src/work/marci/lp/lp_solver_wrapper_2.h
author marci
Tue, 11 Jan 2005 17:16:29 +0000
changeset 1075 789bad021e2d
parent 1031 0b7169db694f
permissions -rw-r--r--
:-}
     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 <stdio>
    12 //#include <stdlib>
    13 extern "C" {
    14 #include "glpk.h"
    15 }
    16 
    17 #include <iostream>
    18 #include <vector>
    19 #include <string>
    20 #include <list>
    21 #include <memory>
    22 #include <utility>
    23 
    24 //#include <sage_graph.h>
    25 //#include <lemon/list_graph.h>
    26 //#include <lemon/graph_wrapper.h>
    27 #include <lemon/invalid.h>
    28 //#include <bfs_dfs.h>
    29 //#include <stp.h>
    30 //#include <lemon/max_flow.h>
    31 //#include <augmenting_flow.h>
    32 //#include <iter_map.h>
    33 
    34 using std::cout;
    35 using std::cin;
    36 using std::endl;
    37 
    38 namespace lemon {
    39 
    40   
    41   /// \addtogroup misc
    42   /// @{
    43 
    44   /// \brief A partitioned vector with iterable classes.
    45   ///
    46   /// This class implements a container in which the data is stored in an 
    47   /// stl vector, the range is partitioned into sets and each set is 
    48   /// doubly linked in a list. 
    49   /// That is, each class is iterable by lemon iterators, and any member of 
    50   /// the vector can bo moved to an other class.
    51   template <typename T>
    52   class IterablePartition {
    53   protected:
    54     struct Node {
    55       T data;
    56       int prev; //invalid az -1
    57       int next; 
    58     };
    59     std::vector<Node> nodes;
    60     struct Tip {
    61       int first;
    62       int last;
    63     };
    64     std::vector<Tip> tips;
    65   public:
    66     /// The classes are indexed by integers from \c 0 to \c classNum()-1.
    67     int classNum() const { return tips.size(); }
    68     /// This lemon style iterator iterates through a class. 
    69     class ClassIt;
    70     /// Constructor. The number of classes is to be given which is fixed 
    71     /// over the life of the container. 
    72     /// The partition classes are indexed from 0 to class_num-1. 
    73     IterablePartition(int class_num) { 
    74       for (int i=0; i<class_num; ++i) {
    75 	Tip t;
    76 	t.first=t.last=-1;
    77 	tips.push_back(t);
    78       }
    79     }
    80   protected:
    81     void befuz(ClassIt it, int class_id) {
    82       if (tips[class_id].first==-1) {
    83 	if (tips[class_id].last==-1) {
    84 	  nodes[it.i].prev=nodes[it.i].next=-1;
    85 	  tips[class_id].first=tips[class_id].last=it.i;
    86 	}
    87       } else {
    88 	nodes[it.i].prev=tips[class_id].last;
    89 	nodes[it.i].next=-1;
    90 	nodes[tips[class_id].last].next=it.i;
    91 	tips[class_id].last=it.i;
    92       }
    93     }
    94     void kifuz(ClassIt it, int class_id) {
    95       if (tips[class_id].first==it.i) {
    96 	if (tips[class_id].last==it.i) {
    97 	  tips[class_id].first=tips[class_id].last=-1;
    98 	} else {
    99 	  tips[class_id].first=nodes[it.i].next;
   100 	  nodes[nodes[it.i].next].prev=-1;
   101 	}
   102       } else {
   103 	if (tips[class_id].last==it.i) {
   104 	  tips[class_id].last=nodes[it.i].prev;
   105 	  nodes[nodes[it.i].prev].next=-1;
   106 	} else {
   107 	  nodes[nodes[it.i].next].prev=nodes[it.i].prev;
   108 	  nodes[nodes[it.i].prev].next=nodes[it.i].next;
   109 	}
   110       }
   111     }
   112   public:
   113     /// A new element with data \c t is pushed into the vector and into class 
   114     /// \c class_id.
   115     ClassIt push_back(const T& t, int class_id) { 
   116       Node n;
   117       n.data=t;
   118       nodes.push_back(n);
   119       int i=nodes.size()-1;
   120       befuz(i, class_id);
   121       return i;
   122     }
   123     /// A member is moved to an other class.
   124     void set(ClassIt it, int old_class_id, int new_class_id) {
   125       kifuz(it.i, old_class_id);
   126       befuz(it.i, new_class_id);
   127     }
   128     /// Returns the data pointed by \c it.
   129     T& operator[](ClassIt it) { return nodes[it.i].data; }
   130     /// Returns the data pointed by \c it.
   131     const T& operator[](ClassIt it) const { return nodes[it.i].data; }
   132     ///.
   133     class ClassIt {
   134       friend class IterablePartition;
   135     protected:
   136       int i;
   137     public:
   138       /// Default constructor.
   139       ClassIt() { }
   140       /// This constructor constructs an iterator which points
   141       /// to the member of th container indexed by the integer _i.
   142       ClassIt(const int& _i) : i(_i) { }
   143       /// Invalid constructor.
   144       ClassIt(const Invalid&) : i(-1) { }
   145     };
   146     /// First member of class \c class_id.
   147     ClassIt& first(ClassIt& it, int class_id) const {
   148       it.i=tips[class_id].first;
   149       return it;
   150     }
   151     /// Next member.
   152     ClassIt& next(ClassIt& it) const {
   153       it.i=nodes[it.i].next;
   154       return it;
   155     }
   156     /// True iff the iterator is valid.
   157     bool valid(const ClassIt& it) const { return it.i!=-1; }
   158   };
   159 
   160   /*! \e
   161    */
   162   template <typename _Value>
   163   class LPSolverBase {
   164   public:
   165     /// \e
   166     typedef _Value Value;
   167     /// \e
   168     typedef IterablePartition<int>::ClassIt RowIt;
   169     /// \e
   170     typedef IterablePartition<int>::ClassIt ColIt;
   171   protected:
   172     /// \e
   173     IterablePartition<int> row_iter_map;
   174     /// \e
   175     IterablePartition<int> col_iter_map;
   176     /// \e
   177     const int VALID_ID;
   178     /// \e
   179     const int INVALID_ID;
   180   public:
   181     /// \e
   182     LPSolverBase() : row_iter_map(2), 
   183 		     col_iter_map(2), 
   184 		     VALID_ID(0), INVALID_ID(1) { }
   185     /// \e
   186     virtual ~LPSolverBase() { }
   187     /// \e
   188     virtual void setMinimize() = 0;
   189     /// \e
   190     virtual void setMaximize() = 0;
   191     /// \e
   192     virtual RowIt addRow() = 0;
   193     /// \e
   194     virtual ColIt addCol() = 0;
   195     /// temporally, glpk style indexing
   196     virtual void setRowCoeffs(RowIt row_it, int num, 
   197 			      int* indices, _Value* doubles) = 0;
   198     //pair<RowIt, _Value>-bol kell megadni egy std range-et
   199     /// \e
   200     template <typename Begin, typename End>
   201     void setRowCoeffs(RowIt row_it, Begin begin, End end) {
   202       int mem_length=1+colNum();
   203       int* indices = new int[mem_length];
   204       _Value* doubles = new _Value[mem_length];
   205       int length=0;
   206       for ( ; begin!=end; ++begin) {
   207 	++length;
   208 	indices[length]=col_iter_map[begin->first];
   209 	doubles[length]=begin->second;
   210       }
   211       setRowCoeffs(row_it, length, indices, doubles);
   212       delete [] indices;
   213       delete [] doubles;
   214     }
   215     /// temporally, glpk style indexing
   216     virtual void setColCoeffs(ColIt col_it, int num, 
   217 			      int* indices, _Value* doubles) = 0;
   218     //pair<ColIt, _Value>-bol kell megadni egy std range-et
   219     /// \e
   220     template <typename Begin, typename End>
   221     void setColCoeffs(ColIt col_it, Begin begin, End end) {
   222       int mem_length=1+rowNum();
   223       int* indices = new int[mem_length];
   224       _Value* doubles = new _Value[mem_length];
   225       int length=0;
   226       for ( ; begin!=end; ++begin) {
   227 	++length;
   228 	indices[length]=row_iter_map[begin->first];
   229 	doubles[length]=begin->second;
   230       }
   231       setColCoeffs(col_it, length, indices, doubles);
   232       delete [] indices;
   233       delete [] doubles;
   234     }
   235     /// \e
   236     virtual void eraseCol(const ColIt& col_it) = 0;
   237     /// \e
   238     virtual void eraseRow(const RowIt& row_it) = 0;
   239     /// \e
   240     virtual void setColBounds(const ColIt& col_it, int bound_type, 
   241 			      _Value lo, _Value up) =0; 
   242     /// \e
   243     virtual _Value getObjCoef(const ColIt& col_it) = 0;
   244     /// \e
   245     virtual void setRowBounds(const RowIt& row_it, int bound_type, 
   246 			      _Value lo, _Value up) = 0;
   247     /// \e
   248     virtual void setObjCoef(const ColIt& col_it, _Value obj_coef) = 0;
   249     /// \e
   250     virtual void solveSimplex() = 0;
   251     /// \e
   252     virtual void solvePrimalSimplex() = 0;
   253     /// \e
   254     virtual void solveDualSimplex() = 0;
   255     /// \e
   256     virtual _Value getPrimal(const ColIt& col_it) = 0;
   257     /// \e
   258     virtual _Value getObjVal() = 0;
   259     /// \e
   260     virtual int rowNum() const = 0;
   261     /// \e
   262     virtual int colNum() const = 0;
   263     /// \e
   264     virtual int warmUp() = 0;
   265     /// \e
   266     virtual void printWarmUpStatus(int i) = 0;
   267     /// \e
   268     virtual int getPrimalStatus() = 0;
   269     /// \e
   270     virtual void printPrimalStatus(int i) = 0;
   271     /// \e
   272     virtual int getDualStatus() = 0;
   273     /// \e
   274     virtual void printDualStatus(int i) = 0;
   275     /// Returns the status of the slack variable assigned to row \c row_it.
   276     virtual int getRowStat(const RowIt& row_it) = 0;
   277     /// \e
   278     virtual void printRowStatus(int i) = 0;
   279     /// Returns the status of the variable assigned to column \c col_it.
   280     virtual int getColStat(const ColIt& col_it) = 0;
   281     /// \e
   282     virtual void printColStatus(int i) = 0;
   283   };
   284   
   285 
   286   /// \brief Wrappers for LP solvers
   287   /// 
   288   /// This class implements a lemon wrapper for glpk.
   289   /// Later other LP-solvers will be wrapped into lemon.
   290   /// The aim of this class is to give a general surface to different 
   291   /// solvers, i.e. it makes possible to write algorithms using LP's, 
   292   /// in which the solver can be changed to an other one easily.
   293   class LPSolverWrapper : public LPSolverBase<double> {
   294   public:
   295     typedef LPSolverBase<double> Parent;
   296 
   297     //   class Row {
   298     //   protected:
   299     //     int i;
   300     //   public:
   301     //     Row() { }
   302     //     Row(const Invalid&) : i(0) { }
   303     //     Row(const int& _i) : i(_i) { }
   304     //     operator int() const { return i; }
   305     //   };
   306     //   class RowIt : public Row {
   307     //   public:
   308     //     RowIt(const Row& row) : Row(row) { }
   309     //   };
   310 
   311     //   class Col {
   312     //   protected:
   313     //     int i;
   314     //   public:
   315     //     Col() { }
   316     //     Col(const Invalid&) : i(0) { }
   317     //     Col(const int& _i) : i(_i) { }
   318     //     operator int() const { return i; }
   319     //   };
   320     //   class ColIt : public Col {
   321     //     ColIt(const Col& col) : Col(col) { }
   322     //   };
   323 
   324   public:
   325     /// \e
   326     LPX* lp;
   327 
   328   public:
   329     /// \e
   330     LPSolverWrapper() : Parent(), 
   331 			lp(lpx_create_prob()) {
   332       lpx_set_int_parm(lp, LPX_K_DUAL, 1);
   333     }
   334     /// \e
   335     ~LPSolverWrapper() {
   336       lpx_delete_prob(lp);
   337     }
   338     /// \e
   339     void setMinimize() { 
   340       lpx_set_obj_dir(lp, LPX_MIN);
   341     }
   342     /// \e
   343     void setMaximize() { 
   344       lpx_set_obj_dir(lp, LPX_MAX);
   345     }
   346     /// \e
   347     ColIt addCol() {
   348       int i=lpx_add_cols(lp, 1);  
   349       ColIt col_it;
   350       col_iter_map.first(col_it, INVALID_ID);
   351       if (col_iter_map.valid(col_it)) { //van hasznalhato hely
   352 	col_iter_map.set(col_it, INVALID_ID, VALID_ID);
   353 	col_iter_map[col_it]=i;
   354 	//col_id_to_lp_col_id[col_iter_map[col_it]]=i;
   355       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
   356 	//col_id_to_lp_col_id.push_back(i);
   357 	//int j=col_id_to_lp_col_id.size()-1;
   358 	col_it=col_iter_map.push_back(i, VALID_ID);
   359       }
   360       //    edge_index_map.set(e, i);
   361       //    lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
   362       //    lpx_set_obj_coef(lp, i, cost[e]);    
   363       return col_it;
   364     }
   365     /// \e
   366     RowIt addRow() {
   367       int i=lpx_add_rows(lp, 1);  
   368       RowIt row_it;
   369       row_iter_map.first(row_it, INVALID_ID);
   370       if (row_iter_map.valid(row_it)) { //van hasznalhato hely
   371 	row_iter_map.set(row_it, INVALID_ID, VALID_ID);
   372 	row_iter_map[row_it]=i;
   373       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
   374 	row_it=row_iter_map.push_back(i, VALID_ID);
   375       }
   376       return row_it;
   377     }
   378     using Parent::setRowCoeffs;
   379     void setRowCoeffs(RowIt row_it, int length, 
   380 		      int* indices, double* doubles) {
   381       lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
   382     }
   383     using Parent::setColCoeffs;
   384     void setColCoeffs(ColIt col_it, int length, 
   385 		      int* indices, double* doubles) {
   386       lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
   387     }
   388     //     //pair<RowIt, double>-bol kell megadni egy std range-et
   389     //     /// \e
   390     //     template <typename Begin, typename End>
   391     //     void setColCoeffs(const ColIt& col_it, 
   392     // 		      Begin begin, End end) {
   393     //       int mem_length=1+lpx_get_num_rows(lp);
   394     //       int* indices = new int[mem_length];
   395     //       double* doubles = new double[mem_length];
   396     //       int length=0;
   397     //       for ( ; begin!=end; ++begin) {
   398     // 	++length;
   399     // 	indices[length]=row_iter_map[begin->first];
   400     // 	doubles[length]=begin->second;
   401     //       }
   402     //       lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
   403     //       delete [] indices;
   404     //       delete [] doubles;
   405     //     }
   406     //     //pair<ColIt, double>-bol kell megadni egy std range-et
   407     //     /// \e
   408     //     template <typename Begin, typename End>
   409     //     void setRowCoeffs(const RowIt& row_it, 
   410     // 		      Begin begin, End end) {
   411     //       int mem_length=1+lpx_get_num_cols(lp);
   412     //       int* indices = new int[mem_length];
   413     //       double* doubles = new double[mem_length];
   414     //       int length=0;
   415     //       for ( ; begin!=end; ++begin) {
   416     // 	++length;
   417     // 	indices[length]=col_iter_map[begin->first];
   418     // 	doubles[length]=begin->second;
   419     //       }
   420     //       lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
   421     //       delete [] indices;
   422     //       delete [] doubles;
   423     //     }
   424     /// \e
   425     void eraseCol(const ColIt& col_it) {
   426       col_iter_map.set(col_it, VALID_ID, INVALID_ID);
   427       int cols[2];
   428       cols[1]=col_iter_map[col_it];
   429       lpx_del_cols(lp, 1, cols);
   430       col_iter_map[col_it]=0; //glpk specifikus
   431       ColIt it;
   432       for (col_iter_map.first(it, VALID_ID); 
   433 	   col_iter_map.valid(it); col_iter_map.next(it)) {
   434 	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
   435       }
   436     }
   437     /// \e
   438     void eraseRow(const RowIt& row_it) {
   439       row_iter_map.set(row_it, VALID_ID, INVALID_ID);
   440       int rows[2];
   441       rows[1]=row_iter_map[row_it];
   442       lpx_del_rows(lp, 1, rows);
   443       row_iter_map[row_it]=0; //glpk specifikus
   444       RowIt it;
   445       for (row_iter_map.first(it, VALID_ID); 
   446 	   row_iter_map.valid(it); row_iter_map.next(it)) {
   447 	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
   448       }
   449     }
   450     /// \e
   451     void setColBounds(const ColIt& col_it, int bound_type, 
   452 		      double lo, double up) {
   453       lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
   454     }
   455     /// \e
   456     double getObjCoef(const ColIt& col_it) { 
   457       return lpx_get_obj_coef(lp, col_iter_map[col_it]);
   458     }
   459     /// \e
   460     void setRowBounds(const RowIt& row_it, int bound_type, 
   461 		      double lo, double up) {
   462       lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
   463     }
   464     /// \e
   465     void setObjCoef(const ColIt& col_it, double obj_coef) { 
   466       lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
   467     }
   468     /// \e
   469     void solveSimplex() { lpx_simplex(lp); }
   470     /// \e
   471     void solvePrimalSimplex() { lpx_simplex(lp); }
   472     /// \e
   473     void solveDualSimplex() { lpx_simplex(lp); }
   474     /// \e
   475     double getPrimal(const ColIt& col_it) {
   476       return lpx_get_col_prim(lp, col_iter_map[col_it]);
   477     }
   478     /// \e
   479     double getObjVal() { return lpx_get_obj_val(lp); }
   480     /// \e
   481     int rowNum() const { return lpx_get_num_rows(lp); }
   482     /// \e
   483     int colNum() const { return lpx_get_num_cols(lp); }
   484     /// \e
   485     int warmUp() { return lpx_warm_up(lp); }
   486     /// \e
   487     void printWarmUpStatus(int i) {
   488       switch (i) {
   489       case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
   490       case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;	
   491       case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
   492       case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
   493       }
   494     }
   495     /// \e
   496     int getPrimalStatus() { return lpx_get_prim_stat(lp); }
   497     /// \e
   498     void printPrimalStatus(int i) {
   499       switch (i) {
   500       case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
   501       case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;	
   502       case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
   503       case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
   504       }
   505     }
   506     /// \e
   507     int getDualStatus() { return lpx_get_dual_stat(lp); }
   508     /// \e
   509     void printDualStatus(int i) {
   510       switch (i) {
   511       case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
   512       case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;	
   513       case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
   514       case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
   515       }
   516     }
   517     /// Returns the status of the slack variable assigned to row \c row_it.
   518     int getRowStat(const RowIt& row_it) { 
   519       return lpx_get_row_stat(lp, row_iter_map[row_it]); 
   520     }
   521     /// \e
   522     void printRowStatus(int i) {
   523       switch (i) {
   524       case LPX_BS: cout << "LPX_BS" << endl; break;
   525       case LPX_NL: cout << "LPX_NL" << endl; break;	
   526       case LPX_NU: cout << "LPX_NU" << endl; break;
   527       case LPX_NF: cout << "LPX_NF" << endl; break;
   528       case LPX_NS: cout << "LPX_NS" << endl; break;
   529       }
   530     }
   531     /// Returns the status of the variable assigned to column \c col_it.
   532     int getColStat(const ColIt& col_it) { 
   533       return lpx_get_col_stat(lp, col_iter_map[col_it]); 
   534     }
   535     /// \e
   536     void printColStatus(int i) {
   537       switch (i) {
   538       case LPX_BS: cout << "LPX_BS" << endl; break;
   539       case LPX_NL: cout << "LPX_NL" << endl; break;	
   540       case LPX_NU: cout << "LPX_NU" << endl; break;
   541       case LPX_NF: cout << "LPX_NF" << endl; break;
   542       case LPX_NS: cout << "LPX_NS" << endl; break;
   543       }
   544     }
   545   };
   546   
   547   /// @}
   548 
   549 } //namespace lemon
   550 
   551 #endif //LEMON_LP_SOLVER_WRAPPER_H