src/work/marci/lp/lp_solver_wrapper_2.h
author marci
Mon, 06 Dec 2004 11:56:10 +0000
changeset 1031 0b7169db694f
child 1048 38a49245a701
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   class LPSolverBase {
   163   public:
   164     /// \e
   165     typedef IterablePartition<int>::ClassIt RowIt;
   166     /// \e
   167     typedef IterablePartition<int>::ClassIt ColIt;
   168   protected:
   169     /// \e
   170     IterablePartition<int> row_iter_map;
   171     /// \e
   172     IterablePartition<int> col_iter_map;
   173     /// \e
   174     const int VALID_ID;
   175     /// \e
   176     const int INVALID_ID;
   177   public:
   178     /// \e
   179     LPSolverBase() : row_iter_map(2), 
   180 		     col_iter_map(2), 
   181 		     VALID_ID(0), INVALID_ID(1) { }
   182     /// \e
   183     virtual ~LPSolverBase() { }
   184     /// \e
   185     virtual void setMinimize() = 0;
   186     /// \e
   187     virtual void setMaximize() = 0;
   188     /// \e
   189     virtual RowIt addRow() = 0;
   190     /// \e
   191     virtual ColIt addCol() = 0;
   192     /// temporally, glpk style indexing
   193     virtual void setRowCoeffs(RowIt row_it, int num, 
   194 			      int* indices, double* doubles) = 0;
   195     //pair<RowIt, double>-bol kell megadni egy std range-et
   196     /// \e
   197     template <typename Begin, typename End>
   198     void setRowCoeffs(RowIt row_it, Begin begin, End end) {
   199       int mem_length=1+colNum();
   200       int* indices = new int[mem_length];
   201       double* doubles = new double[mem_length];
   202       int length=0;
   203       for ( ; begin!=end; ++begin) {
   204 	++length;
   205 	indices[length]=col_iter_map[begin->first];
   206 	doubles[length]=begin->second;
   207       }
   208       setRowCoeffs(row_it, length, indices, doubles);
   209       delete [] indices;
   210       delete [] doubles;
   211     }
   212     /// temporally, glpk style indexing
   213     virtual void setColCoeffs(ColIt col_it, int num, 
   214 			      int* indices, double* doubles) = 0;
   215     //pair<ColIt, double>-bol kell megadni egy std range-et
   216     /// \e
   217     template <typename Begin, typename End>
   218     void setColCoeffs(ColIt col_it, Begin begin, End end) {
   219       int mem_length=1+rowNum();
   220       int* indices = new int[mem_length];
   221       double* doubles = new double[mem_length];
   222       int length=0;
   223       for ( ; begin!=end; ++begin) {
   224 	++length;
   225 	indices[length]=row_iter_map[begin->first];
   226 	doubles[length]=begin->second;
   227       }
   228       setColCoeffs(col_it, length, indices, doubles);
   229       delete [] indices;
   230       delete [] doubles;
   231     }
   232     /// \e
   233     virtual void eraseCol(const ColIt& col_it) = 0;
   234     /// \e
   235     virtual void eraseRow(const RowIt& row_it) = 0;
   236     /// \e
   237     virtual void setColBounds(const ColIt& col_it, int bound_type, 
   238 			      double lo, double up) =0; 
   239     /// \e
   240     virtual double getObjCoef(const ColIt& col_it) = 0;
   241     /// \e
   242     virtual void setRowBounds(const RowIt& row_it, int bound_type, 
   243 			      double lo, double up) = 0;
   244     /// \e
   245     virtual void setObjCoef(const ColIt& col_it, double obj_coef) = 0;
   246     /// \e
   247     virtual void solveSimplex() = 0;
   248     /// \e
   249     virtual void solvePrimalSimplex() = 0;
   250     /// \e
   251     virtual void solveDualSimplex() = 0;
   252     /// \e
   253     virtual double getPrimal(const ColIt& col_it) = 0;
   254     /// \e
   255     virtual double getObjVal() = 0;
   256     /// \e
   257     virtual int rowNum() const = 0;
   258     /// \e
   259     virtual int colNum() const = 0;
   260     /// \e
   261     virtual int warmUp() = 0;
   262     /// \e
   263     virtual void printWarmUpStatus(int i) = 0;
   264     /// \e
   265     virtual int getPrimalStatus() = 0;
   266     /// \e
   267     virtual void printPrimalStatus(int i) = 0;
   268     /// \e
   269     virtual int getDualStatus() = 0;
   270     /// \e
   271     virtual void printDualStatus(int i) = 0;
   272     /// Returns the status of the slack variable assigned to row \c row_it.
   273     virtual int getRowStat(const RowIt& row_it) = 0;
   274     /// \e
   275     virtual void printRowStatus(int i) = 0;
   276     /// Returns the status of the variable assigned to column \c col_it.
   277     virtual int getColStat(const ColIt& col_it) = 0;
   278     /// \e
   279     virtual void printColStatus(int i) = 0;
   280   };
   281   
   282   /// \brief Wrappers for LP solvers
   283   /// 
   284   /// This class implements a lemon wrapper for glpk.
   285   /// Later other LP-solvers will be wrapped into lemon.
   286   /// The aim of this class is to give a general surface to different 
   287   /// solvers, i.e. it makes possible to write algorithms using LP's, 
   288   /// in which the solver can be changed to an other one easily.
   289   class LPSolverWrapper : public LPSolverBase {
   290   public:
   291     typedef LPSolverBase Parent;
   292 
   293     //   class Row {
   294     //   protected:
   295     //     int i;
   296     //   public:
   297     //     Row() { }
   298     //     Row(const Invalid&) : i(0) { }
   299     //     Row(const int& _i) : i(_i) { }
   300     //     operator int() const { return i; }
   301     //   };
   302     //   class RowIt : public Row {
   303     //   public:
   304     //     RowIt(const Row& row) : Row(row) { }
   305     //   };
   306 
   307     //   class Col {
   308     //   protected:
   309     //     int i;
   310     //   public:
   311     //     Col() { }
   312     //     Col(const Invalid&) : i(0) { }
   313     //     Col(const int& _i) : i(_i) { }
   314     //     operator int() const { return i; }
   315     //   };
   316     //   class ColIt : public Col {
   317     //     ColIt(const Col& col) : Col(col) { }
   318     //   };
   319 
   320   public:
   321     /// \e
   322     LPX* lp;
   323 
   324   public:
   325     /// \e
   326     LPSolverWrapper() : LPSolverBase(), 
   327 			lp(lpx_create_prob()) {
   328       lpx_set_int_parm(lp, LPX_K_DUAL, 1);
   329     }
   330     /// \e
   331     ~LPSolverWrapper() {
   332       lpx_delete_prob(lp);
   333     }
   334     /// \e
   335     void setMinimize() { 
   336       lpx_set_obj_dir(lp, LPX_MIN);
   337     }
   338     /// \e
   339     void setMaximize() { 
   340       lpx_set_obj_dir(lp, LPX_MAX);
   341     }
   342     /// \e
   343     ColIt addCol() {
   344       int i=lpx_add_cols(lp, 1);  
   345       ColIt col_it;
   346       col_iter_map.first(col_it, INVALID_ID);
   347       if (col_iter_map.valid(col_it)) { //van hasznalhato hely
   348 	col_iter_map.set(col_it, INVALID_ID, VALID_ID);
   349 	col_iter_map[col_it]=i;
   350 	//col_id_to_lp_col_id[col_iter_map[col_it]]=i;
   351       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
   352 	//col_id_to_lp_col_id.push_back(i);
   353 	//int j=col_id_to_lp_col_id.size()-1;
   354 	col_it=col_iter_map.push_back(i, VALID_ID);
   355       }
   356       //    edge_index_map.set(e, i);
   357       //    lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
   358       //    lpx_set_obj_coef(lp, i, cost[e]);    
   359       return col_it;
   360     }
   361     /// \e
   362     RowIt addRow() {
   363       int i=lpx_add_rows(lp, 1);  
   364       RowIt row_it;
   365       row_iter_map.first(row_it, INVALID_ID);
   366       if (row_iter_map.valid(row_it)) { //van hasznalhato hely
   367 	row_iter_map.set(row_it, INVALID_ID, VALID_ID);
   368 	row_iter_map[row_it]=i;
   369       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
   370 	row_it=row_iter_map.push_back(i, VALID_ID);
   371       }
   372       return row_it;
   373     }
   374     using Parent::setRowCoeffs;
   375     void setRowCoeffs(RowIt row_it, int length, 
   376 		      int* indices, double* doubles) {
   377       lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
   378     }
   379     using Parent::setColCoeffs;
   380     void setColCoeffs(ColIt col_it, int length, 
   381 		      int* indices, double* doubles) {
   382       lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
   383     }
   384     //     //pair<RowIt, double>-bol kell megadni egy std range-et
   385     //     /// \e
   386     //     template <typename Begin, typename End>
   387     //     void setColCoeffs(const ColIt& col_it, 
   388     // 		      Begin begin, End end) {
   389     //       int mem_length=1+lpx_get_num_rows(lp);
   390     //       int* indices = new int[mem_length];
   391     //       double* doubles = new double[mem_length];
   392     //       int length=0;
   393     //       for ( ; begin!=end; ++begin) {
   394     // 	++length;
   395     // 	indices[length]=row_iter_map[begin->first];
   396     // 	doubles[length]=begin->second;
   397     //       }
   398     //       lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
   399     //       delete [] indices;
   400     //       delete [] doubles;
   401     //     }
   402     //     //pair<ColIt, double>-bol kell megadni egy std range-et
   403     //     /// \e
   404     //     template <typename Begin, typename End>
   405     //     void setRowCoeffs(const RowIt& row_it, 
   406     // 		      Begin begin, End end) {
   407     //       int mem_length=1+lpx_get_num_cols(lp);
   408     //       int* indices = new int[mem_length];
   409     //       double* doubles = new double[mem_length];
   410     //       int length=0;
   411     //       for ( ; begin!=end; ++begin) {
   412     // 	++length;
   413     // 	indices[length]=col_iter_map[begin->first];
   414     // 	doubles[length]=begin->second;
   415     //       }
   416     //       lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
   417     //       delete [] indices;
   418     //       delete [] doubles;
   419     //     }
   420     /// \e
   421     void eraseCol(const ColIt& col_it) {
   422       col_iter_map.set(col_it, VALID_ID, INVALID_ID);
   423       int cols[2];
   424       cols[1]=col_iter_map[col_it];
   425       lpx_del_cols(lp, 1, cols);
   426       col_iter_map[col_it]=0; //glpk specifikus
   427       ColIt it;
   428       for (col_iter_map.first(it, VALID_ID); 
   429 	   col_iter_map.valid(it); col_iter_map.next(it)) {
   430 	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
   431       }
   432     }
   433     /// \e
   434     void eraseRow(const RowIt& row_it) {
   435       row_iter_map.set(row_it, VALID_ID, INVALID_ID);
   436       int rows[2];
   437       rows[1]=row_iter_map[row_it];
   438       lpx_del_rows(lp, 1, rows);
   439       row_iter_map[row_it]=0; //glpk specifikus
   440       RowIt it;
   441       for (row_iter_map.first(it, VALID_ID); 
   442 	   row_iter_map.valid(it); row_iter_map.next(it)) {
   443 	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
   444       }
   445     }
   446     /// \e
   447     void setColBounds(const ColIt& col_it, int bound_type, 
   448 		      double lo, double up) {
   449       lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
   450     }
   451     /// \e
   452     double getObjCoef(const ColIt& col_it) { 
   453       return lpx_get_obj_coef(lp, col_iter_map[col_it]);
   454     }
   455     /// \e
   456     void setRowBounds(const RowIt& row_it, int bound_type, 
   457 		      double lo, double up) {
   458       lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
   459     }
   460     /// \e
   461     void setObjCoef(const ColIt& col_it, double obj_coef) { 
   462       lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
   463     }
   464     /// \e
   465     void solveSimplex() { lpx_simplex(lp); }
   466     /// \e
   467     void solvePrimalSimplex() { lpx_simplex(lp); }
   468     /// \e
   469     void solveDualSimplex() { lpx_simplex(lp); }
   470     /// \e
   471     double getPrimal(const ColIt& col_it) {
   472       return lpx_get_col_prim(lp, col_iter_map[col_it]);
   473     }
   474     /// \e
   475     double getObjVal() { return lpx_get_obj_val(lp); }
   476     /// \e
   477     int rowNum() const { return lpx_get_num_rows(lp); }
   478     /// \e
   479     int colNum() const { return lpx_get_num_cols(lp); }
   480     /// \e
   481     int warmUp() { return lpx_warm_up(lp); }
   482     /// \e
   483     void printWarmUpStatus(int i) {
   484       switch (i) {
   485       case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
   486       case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;	
   487       case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
   488       case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
   489       }
   490     }
   491     /// \e
   492     int getPrimalStatus() { return lpx_get_prim_stat(lp); }
   493     /// \e
   494     void printPrimalStatus(int i) {
   495       switch (i) {
   496       case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
   497       case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;	
   498       case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
   499       case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
   500       }
   501     }
   502     /// \e
   503     int getDualStatus() { return lpx_get_dual_stat(lp); }
   504     /// \e
   505     void printDualStatus(int i) {
   506       switch (i) {
   507       case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
   508       case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;	
   509       case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
   510       case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
   511       }
   512     }
   513     /// Returns the status of the slack variable assigned to row \c row_it.
   514     int getRowStat(const RowIt& row_it) { 
   515       return lpx_get_row_stat(lp, row_iter_map[row_it]); 
   516     }
   517     /// \e
   518     void printRowStatus(int i) {
   519       switch (i) {
   520       case LPX_BS: cout << "LPX_BS" << endl; break;
   521       case LPX_NL: cout << "LPX_NL" << endl; break;	
   522       case LPX_NU: cout << "LPX_NU" << endl; break;
   523       case LPX_NF: cout << "LPX_NF" << endl; break;
   524       case LPX_NS: cout << "LPX_NS" << endl; break;
   525       }
   526     }
   527     /// Returns the status of the variable assigned to column \c col_it.
   528     int getColStat(const ColIt& col_it) { 
   529       return lpx_get_col_stat(lp, col_iter_map[col_it]); 
   530     }
   531     /// \e
   532     void printColStatus(int i) {
   533       switch (i) {
   534       case LPX_BS: cout << "LPX_BS" << endl; break;
   535       case LPX_NL: cout << "LPX_NL" << endl; break;	
   536       case LPX_NU: cout << "LPX_NU" << endl; break;
   537       case LPX_NF: cout << "LPX_NF" << endl; break;
   538       case LPX_NS: cout << "LPX_NS" << endl; break;
   539       }
   540     }
   541   };
   542   
   543   /// @}
   544 
   545 } //namespace lemon
   546 
   547 #endif //LEMON_LP_SOLVER_WRAPPER_H