COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/lp/lp_solver_wrapper_3.h @ 1104:23a54f889272

Last change on this file since 1104:23a54f889272 was 1104:23a54f889272, checked in by marci, 19 years ago

small changes, a try for max flow using expression

File size: 18.2 KB
Line 
1// -*- c++ -*-
2#ifndef LEMON_LP_SOLVER_WRAPPER_H
3#define LEMON_LP_SOLVER_WRAPPER_H
4
5///\ingroup misc
6///\file
7///\brief Dijkstra algorithm.
8
9// #include <stdio.h>
10#include <stdlib.h>
11#include <iostream>
12#include <map>
13#include <limits>
14// #include <stdio>
15//#include <stdlib>
16extern "C" {
17#include "glpk.h"
18}
19
20#include <iostream>
21#include <vector>
22#include <string>
23#include <list>
24#include <memory>
25#include <utility>
26
27//#include <sage_graph.h>
28//#include <lemon/list_graph.h>
29//#include <lemon/graph_wrapper.h>
30#include <lemon/invalid.h>
31#include <expression.h>
32//#include <bfs_dfs.h>
33//#include <stp.h>
34//#include <lemon/max_flow.h>
35//#include <augmenting_flow.h>
36//#include <iter_map.h>
37
38using std::cout;
39using std::cin;
40using std::endl;
41
42namespace lemon {
43 
44  /// \addtogroup misc
45  /// @{
46
47  /// \brief A partitioned vector with iterable classes.
48  ///
49  /// This class implements a container in which the data is stored in an
50  /// stl vector, the range is partitioned into sets and each set is
51  /// doubly linked in a list.
52  /// That is, each class is iterable by lemon iterators, and any member of
53  /// the vector can bo moved to an other class.
54  template <typename T>
55  class IterablePartition {
56  protected:
57    struct Node {
58      T data;
59      int prev; //invalid az -1
60      int next;
61    };
62    std::vector<Node> nodes;
63    struct Tip {
64      int first;
65      int last;
66    };
67    std::vector<Tip> tips;
68  public:
69    /// The classes are indexed by integers from \c 0 to \c classNum()-1.
70    int classNum() const { return tips.size(); }
71    /// This lemon style iterator iterates through a class.
72    class ClassIt;
73    /// Constructor. The number of classes is to be given which is fixed
74    /// over the life of the container.
75    /// The partition classes are indexed from 0 to class_num-1.
76    IterablePartition(int class_num) {
77      for (int i=0; i<class_num; ++i) {
78        Tip t;
79        t.first=t.last=-1;
80        tips.push_back(t);
81      }
82    }
83  protected:
84    void befuz(ClassIt it, int class_id) {
85      if (tips[class_id].first==-1) {
86        if (tips[class_id].last==-1) {
87          nodes[it.i].prev=nodes[it.i].next=-1;
88          tips[class_id].first=tips[class_id].last=it.i;
89        }
90      } else {
91        nodes[it.i].prev=tips[class_id].last;
92        nodes[it.i].next=-1;
93        nodes[tips[class_id].last].next=it.i;
94        tips[class_id].last=it.i;
95      }
96    }
97    void kifuz(ClassIt it, int class_id) {
98      if (tips[class_id].first==it.i) {
99        if (tips[class_id].last==it.i) {
100          tips[class_id].first=tips[class_id].last=-1;
101        } else {
102          tips[class_id].first=nodes[it.i].next;
103          nodes[nodes[it.i].next].prev=-1;
104        }
105      } else {
106        if (tips[class_id].last==it.i) {
107          tips[class_id].last=nodes[it.i].prev;
108          nodes[nodes[it.i].prev].next=-1;
109        } else {
110          nodes[nodes[it.i].next].prev=nodes[it.i].prev;
111          nodes[nodes[it.i].prev].next=nodes[it.i].next;
112        }
113      }
114    }
115  public:
116    /// A new element with data \c t is pushed into the vector and into class
117    /// \c class_id.
118    ClassIt push_back(const T& t, int class_id) {
119      Node n;
120      n.data=t;
121      nodes.push_back(n);
122      int i=nodes.size()-1;
123      befuz(i, class_id);
124      return i;
125    }
126    /// A member is moved to an other class.
127    void set(ClassIt it, int old_class_id, int new_class_id) {
128      kifuz(it.i, old_class_id);
129      befuz(it.i, new_class_id);
130    }
131    /// Returns the data pointed by \c it.
132    T& operator[](ClassIt it) { return nodes[it.i].data; }
133    /// Returns the data pointed by \c it.
134    const T& operator[](ClassIt it) const { return nodes[it.i].data; }
135    ///.
136    class ClassIt {
137      friend class IterablePartition;
138    protected:
139      int i;
140    public:
141      /// Default constructor.
142      ClassIt() { }
143      /// This constructor constructs an iterator which points
144      /// to the member of th container indexed by the integer _i.
145      ClassIt(const int& _i) : i(_i) { }
146      /// Invalid constructor.
147      ClassIt(const Invalid&) : i(-1) { }
148      friend bool operator<(const ClassIt& x, const ClassIt& y);
149      friend std::ostream& operator<<(std::ostream& os,
150                                      const ClassIt& it);
151    };
152    friend bool operator<(const ClassIt& x, const ClassIt& y) {
153      return (x.i < y.i);
154    }
155    friend std::ostream& operator<<(std::ostream& os,
156                                    const ClassIt& it) {
157      os << it.i;
158      return os;
159    }
160    /// First member of class \c class_id.
161    ClassIt& first(ClassIt& it, int class_id) const {
162      it.i=tips[class_id].first;
163      return it;
164    }
165    /// Next member.
166    ClassIt& next(ClassIt& it) const {
167      it.i=nodes[it.i].next;
168      return it;
169    }
170    /// True iff the iterator is valid.
171    bool valid(const ClassIt& it) const { return it.i!=-1; }
172  };
173
174
175  /*! \e
176
177  \todo A[x,y]-t cserel. Jobboldal, baloldal csere.
178  \todo LEKERDEZESEK!!!
179  \todo DOKSI!!!! Doxygen group!!!
180
181   */
182  template <typename _Value>
183  class LPSolverBase {
184  public:
185    /// \e
186    typedef _Value Value;
187    /// \e
188    typedef IterablePartition<int>::ClassIt RowIt;
189    /// \e
190    typedef IterablePartition<int>::ClassIt ColIt;
191  public:
192    /// \e
193    IterablePartition<int> row_iter_map;
194    /// \e
195    IterablePartition<int> col_iter_map;
196    /// \e
197    const int VALID_CLASS;
198    /// \e
199    const int INVALID_CLASS;
200    /// \e
201    static const _Value INF;
202  public:
203    /// \e
204    LPSolverBase() : row_iter_map(2),
205                     col_iter_map(2),
206                     VALID_CLASS(0), INVALID_CLASS(1) { }
207    /// \e
208    virtual ~LPSolverBase() { }
209
210    //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
211
212  public:
213    /// \e
214    virtual void setMinimize() = 0;
215    /// \e
216    virtual void setMaximize() = 0;
217
218    //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
219
220  protected:
221    /// \e
222    virtual int _addRow() = 0;
223    /// \e
224    virtual int _addCol() = 0;
225    /// \e
226    virtual void _setRowCoeffs(int i,
227                               const std::vector<std::pair<int, _Value> >& coeffs) = 0;
228    /// \e
229    virtual void _setColCoeffs(int i,
230                               const std::vector<std::pair<int, _Value> >& coeffs) = 0;
231    /// \e
232    virtual void _eraseCol(int i) = 0;
233    /// \e
234    virtual void _eraseRow(int i) = 0;
235  public:
236    /// \e
237    enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };
238  protected:
239    /// \e
240    virtual void _setColBounds(int i, Bound bound,
241                               _Value lo, _Value up) = 0;
242    /// \e
243    virtual void _setRowBounds(int i, Bound bound,
244                               _Value lo, _Value up) = 0;
245    /// \e
246    virtual void _setObjCoef(int i, _Value obj_coef) = 0;
247    /// \e
248    virtual _Value _getObjCoef(int i) = 0;
249
250    //LOW LEVEL, SOLUTION RETRIEVING FUNCTIONS
251
252  protected:
253    virtual _Value _getPrimal(int i) = 0;
254
255    //HIGH LEVEL INTERFACE, MATRIX MANIPULATING FUNTIONS
256
257  public:
258    /// \e
259    RowIt addRow() {
260      int i=_addRow();
261      RowIt row_it;
262      row_iter_map.first(row_it, INVALID_CLASS);
263      if (row_iter_map.valid(row_it)) { //van hasznalhato hely
264        row_iter_map.set(row_it, INVALID_CLASS, VALID_CLASS);
265        row_iter_map[row_it]=i;
266      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
267        row_it=row_iter_map.push_back(i, VALID_CLASS);
268      }
269      return row_it;
270    }
271    /// \e
272    ColIt addCol() {
273      int i=_addCol(); 
274      ColIt col_it;
275      col_iter_map.first(col_it, INVALID_CLASS);
276      if (col_iter_map.valid(col_it)) { //van hasznalhato hely
277        col_iter_map.set(col_it, INVALID_CLASS, VALID_CLASS);
278        col_iter_map[col_it]=i;
279      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
280        col_it=col_iter_map.push_back(i, VALID_CLASS);
281      }
282      return col_it;
283    }
284    /// \e
285    template <typename Begin, typename End>
286    void setRowCoeffs(RowIt row_it, Begin begin, End end) {
287      std::vector<std::pair<int, double> > coeffs;
288      for ( ; begin!=end; ++begin) {
289        coeffs.push_back(std::
290                         make_pair(col_iter_map[begin->first], begin->second));
291      }
292      _setRowCoeffs(row_iter_map[row_it], coeffs);
293    }
294    /// \e
295    template <typename Begin, typename End>
296    void setColCoeffs(ColIt col_it, Begin begin, End end) {
297      std::vector<std::pair<int, double> > coeffs;
298      for ( ; begin!=end; ++begin) {
299        coeffs.push_back(std::
300                         make_pair(row_iter_map[begin->first], begin->second));
301      }
302      _setColCoeffs(col_iter_map[col_it], coeffs);
303    }
304    /// \e
305    void eraseCol(const ColIt& col_it) {
306      col_iter_map.set(col_it, VALID_CLASS, INVALID_CLASS);
307      int cols[2];
308      cols[1]=col_iter_map[col_it];
309      _eraseCol(cols[1]);
310      col_iter_map[col_it]=0; //glpk specifikus, de kell ez??
311      ColIt it;
312      for (col_iter_map.first(it, VALID_CLASS);
313           col_iter_map.valid(it); col_iter_map.next(it)) {
314        if (col_iter_map[it]>cols[1]) --col_iter_map[it];
315      }
316    }
317    /// \e
318    void eraseRow(const RowIt& row_it) {
319      row_iter_map.set(row_it, VALID_CLASS, INVALID_CLASS);
320      int rows[2];
321      rows[1]=row_iter_map[row_it];
322      _eraseRow(rows[1]);
323      row_iter_map[row_it]=0; //glpk specifikus, de kell ez??
324      RowIt it;
325      for (row_iter_map.first(it, VALID_CLASS);
326           row_iter_map.valid(it); row_iter_map.next(it)) {
327        if (row_iter_map[it]>rows[1]) --row_iter_map[it];
328      }
329    }
330    /// \e
331    void setColBounds(const ColIt& col_it, Bound bound,
332                      _Value lo, _Value up) {
333      _setColBounds(col_iter_map[col_it], bound, lo, up);
334    }
335    /// \e
336    void setRowBounds(const RowIt& row_it, Bound bound,
337                      _Value lo, _Value up) {
338      _setRowBounds(row_iter_map[row_it], bound, lo, up);
339    }
340    /// \e
341    void setObjCoef(const ColIt& col_it, _Value obj_coef) {
342      _setObjCoef(col_iter_map[col_it], obj_coef);
343    }
344    /// \e
345    _Value getObjCoef(const ColIt& col_it) {
346      return _getObjCoef(col_iter_map[col_it]);
347    }
348
349    //MOST HIGH LEVEL, USER FRIEND FUNCTIONS
350
351    /// \e
352    typedef Expr<ColIt, _Value> Expression;
353    /// \e
354    typedef Expr<RowIt, _Value> DualExpression;
355    /// \e
356    void setRowCoeffs(RowIt row_it, const Expression& expr) {
357      std::vector<std::pair<int, _Value> > row_coeffs;
358      for(typename Expression::Data::const_iterator i=expr.data.begin();
359          i!=expr.data.end(); ++i) {
360        row_coeffs.push_back(std::make_pair
361                             (col_iter_map[(*i).first], (*i).second));
362      }
363      _setRowCoeffs(row_iter_map[row_it], row_coeffs);
364    }
365    /// \e
366    void setColCoeffs(ColIt col_it, const DualExpression& expr) {
367      std::vector<std::pair<int, _Value> > col_coeffs;
368      for(typename DualExpression::Data::const_iterator i=expr.data.begin();
369          i!=expr.data.end(); ++i) {
370        col_coeffs.push_back(std::make_pair
371                             (row_iter_map[(*i).first], (*i).second));
372      }
373      _setColCoeffs(col_iter_map[col_it], col_coeffs);
374    }
375    /// \e
376    void setObjCoeffs(const Expression& expr) {
377      for(typename Expression::Data::const_iterator i=expr.data.begin();
378          i!=expr.data.end(); ++i) {
379        setObjCoef((*i).first, (*i).second);
380      }
381    }
382    //SOLVER FUNCTIONS
383
384    /// \e
385    virtual void solveSimplex() = 0;
386    /// \e
387    virtual void solvePrimalSimplex() = 0;
388    /// \e
389    virtual void solveDualSimplex() = 0;
390    /// \e
391
392    //HIGH LEVEL, SOLUTION RETRIEVING FUNCTIONS
393
394  public:
395    _Value getPrimal(const ColIt& col_it) {
396      return _getPrimal(col_iter_map[col_it]);
397    }
398    /// \e
399    virtual _Value getObjVal() = 0;
400
401    //OTHER FUNCTIONS
402
403    /// \e
404    virtual int rowNum() const = 0;
405    /// \e
406    virtual int colNum() const = 0;
407    /// \e
408    virtual int warmUp() = 0;
409    /// \e
410    virtual void printWarmUpStatus(int i) = 0;
411    /// \e
412    virtual int getPrimalStatus() = 0;
413    /// \e
414    virtual void printPrimalStatus(int i) = 0;
415    /// \e
416    virtual int getDualStatus() = 0;
417    /// \e
418    virtual void printDualStatus(int i) = 0;
419    /// Returns the status of the slack variable assigned to row \c row_it.
420    virtual int getRowStat(const RowIt& row_it) = 0;
421    /// \e
422    virtual void printRowStatus(int i) = 0;
423    /// Returns the status of the variable assigned to column \c col_it.
424    virtual int getColStat(const ColIt& col_it) = 0;
425    /// \e
426    virtual void printColStatus(int i) = 0;
427  };
428 
429  template <typename _Value>
430  const _Value LPSolverBase<_Value>::INF=std::numeric_limits<_Value>::infinity();
431
432
433  /// \brief Wrappers for LP solvers
434  ///
435  /// This class implements a lemon wrapper for glpk.
436  /// Later other LP-solvers will be wrapped into lemon.
437  /// The aim of this class is to give a general surface to different
438  /// solvers, i.e. it makes possible to write algorithms using LP's,
439  /// in which the solver can be changed to an other one easily.
440  class LPGLPK : public LPSolverBase<double> {
441  public:
442    typedef LPSolverBase<double> Parent;
443
444  public:
445    /// \e
446    LPX* lp;
447
448  public:
449    /// \e
450    LPGLPK() : Parent(),
451                        lp(lpx_create_prob()) {
452      lpx_set_int_parm(lp, LPX_K_DUAL, 1);
453    }
454    /// \e
455    ~LPGLPK() {
456      lpx_delete_prob(lp);
457    }
458
459    //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
460
461    /// \e
462    void setMinimize() {
463      lpx_set_obj_dir(lp, LPX_MIN);
464    }
465    /// \e
466    void setMaximize() {
467      lpx_set_obj_dir(lp, LPX_MAX);
468    }
469
470    //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
471
472  protected:
473    /// \e
474    int _addCol() {
475      return lpx_add_cols(lp, 1);
476    }
477    /// \e
478    int _addRow() {
479      return lpx_add_rows(lp, 1);
480    }
481    /// \e
482    virtual void _setRowCoeffs(int i,
483                               const std::vector<std::pair<int, double> >& coeffs) {
484      int mem_length=1+colNum();
485      int* indices = new int[mem_length];
486      double* doubles = new double[mem_length];
487      int length=0;
488      for (std::vector<std::pair<int, double> >::
489             const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
490        ++length;
491        indices[length]=it->first;
492        doubles[length]=it->second;
493//      std::cout << "  " << indices[length] << " "
494//                << doubles[length] << std::endl;
495      }
496//      std::cout << i << " " << length << std::endl;
497      lpx_set_mat_row(lp, i, length, indices, doubles);
498      delete [] indices;
499      delete [] doubles;
500    }
501    /// \e
502    virtual void _setColCoeffs(int i,
503                               const std::vector<std::pair<int, double> >& coeffs) {
504      int mem_length=1+rowNum();
505      int* indices = new int[mem_length];
506      double* doubles = new double[mem_length];
507      int length=0;
508      for (std::vector<std::pair<int, double> >::
509             const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
510        ++length;
511        indices[length]=it->first;
512        doubles[length]=it->second;
513      }
514      lpx_set_mat_col(lp, i, length, indices, doubles);
515      delete [] indices;
516      delete [] doubles;
517    }
518    /// \e
519    virtual void _eraseCol(int i) {
520      int cols[2];
521      cols[1]=i;
522      lpx_del_cols(lp, 1, cols);
523    }
524    virtual void _eraseRow(int i) {
525      int rows[2];
526      rows[1]=i;
527      lpx_del_rows(lp, 1, rows);
528    }
529    virtual void _setColBounds(int i, Bound bound,
530                               double lo, double up) {
531      switch (bound) {
532      case FREE:
533        lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
534        break;
535      case LOWER:
536        lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
537        break;
538      case UPPER:
539        lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
540        break;
541      case DOUBLE:
542        lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
543        break;
544      case FIXED:
545        lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
546        break;
547      }
548    }
549    virtual void _setRowBounds(int i, Bound bound,
550                               double lo, double up) {
551      switch (bound) {
552      case FREE:
553        lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
554        break;
555      case LOWER:
556        lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
557        break;
558      case UPPER:
559        lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
560        break;
561      case DOUBLE:
562        lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
563        break;
564      case FIXED:
565        lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
566        break;
567      }
568    }
569  protected:
570    /// \e
571    virtual double _getObjCoef(int i) {
572      return lpx_get_obj_coef(lp, i);
573    }
574    /// \e
575    virtual void _setObjCoef(int i, double obj_coef) {
576      lpx_set_obj_coef(lp, i, obj_coef);
577    }
578  public:
579    /// \e
580    void solveSimplex() { lpx_simplex(lp); }
581    /// \e
582    void solvePrimalSimplex() { lpx_simplex(lp); }
583    /// \e
584    void solveDualSimplex() { lpx_simplex(lp); }
585    /// \e
586  protected:
587    virtual double _getPrimal(int i) {
588      return lpx_get_col_prim(lp, i);
589    }
590  public:
591    /// \e
592    double getObjVal() { return lpx_get_obj_val(lp); }
593    /// \e
594    int rowNum() const { return lpx_get_num_rows(lp); }
595    /// \e
596    int colNum() const { return lpx_get_num_cols(lp); }
597    /// \e
598    int warmUp() { return lpx_warm_up(lp); }
599    /// \e
600    void printWarmUpStatus(int i) {
601      switch (i) {
602      case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
603      case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;   
604      case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
605      case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
606      }
607    }
608    /// \e
609    int getPrimalStatus() { return lpx_get_prim_stat(lp); }
610    /// \e
611    void printPrimalStatus(int i) {
612      switch (i) {
613      case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
614      case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;     
615      case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
616      case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
617      }
618    }
619    /// \e
620    int getDualStatus() { return lpx_get_dual_stat(lp); }
621    /// \e
622    void printDualStatus(int i) {
623      switch (i) {
624      case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
625      case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;     
626      case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
627      case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
628      }
629    }
630    /// Returns the status of the slack variable assigned to row \c row_it.
631    int getRowStat(const RowIt& row_it) {
632      return lpx_get_row_stat(lp, row_iter_map[row_it]);
633    }
634    /// \e
635    void printRowStatus(int i) {
636      switch (i) {
637      case LPX_BS: cout << "LPX_BS" << endl; break;
638      case LPX_NL: cout << "LPX_NL" << endl; break;     
639      case LPX_NU: cout << "LPX_NU" << endl; break;
640      case LPX_NF: cout << "LPX_NF" << endl; break;
641      case LPX_NS: cout << "LPX_NS" << endl; break;
642      }
643    }
644    /// Returns the status of the variable assigned to column \c col_it.
645    int getColStat(const ColIt& col_it) {
646      return lpx_get_col_stat(lp, col_iter_map[col_it]);
647    }
648    /// \e
649    void printColStatus(int i) {
650      switch (i) {
651      case LPX_BS: cout << "LPX_BS" << endl; break;
652      case LPX_NL: cout << "LPX_NL" << endl; break;     
653      case LPX_NU: cout << "LPX_NU" << endl; break;
654      case LPX_NF: cout << "LPX_NF" << endl; break;
655      case LPX_NS: cout << "LPX_NS" << endl; break;
656      }
657    }
658  };
659 
660  /// @}
661
662} //namespace lemon
663
664#endif //LEMON_LP_SOLVER_WRAPPER_H
Note: See TracBrowser for help on using the repository browser.