marci@1031: // -*- c++ -*-
marci@1031: #ifndef LEMON_LP_SOLVER_WRAPPER_H
marci@1031: #define LEMON_LP_SOLVER_WRAPPER_H
marci@1031: 
marci@1031: ///\ingroup misc
marci@1031: ///\file
marci@1031: ///\brief Dijkstra algorithm.
marci@1031: 
marci@1031: // #include <stdio.h>
marci@1031: #include <stdlib.h>
marci@1097: #include <iostream>
marci@1097: #include <map>
marci@1104: #include <limits>
marci@1031: // #include <stdio>
marci@1031: //#include <stdlib>
marci@1031: extern "C" {
marci@1031: #include "glpk.h"
marci@1031: }
marci@1031: 
marci@1031: #include <iostream>
marci@1031: #include <vector>
marci@1031: #include <string>
marci@1031: #include <list>
marci@1031: #include <memory>
marci@1031: #include <utility>
marci@1031: 
marci@1031: //#include <lemon/list_graph.h>
marci@1031: #include <lemon/invalid.h>
marci@1099: #include <expression.h>
marci@1031: //#include <bfs_dfs.h>
marci@1031: //#include <stp.h>
marci@1031: //#include <lemon/max_flow.h>
marci@1031: //#include <augmenting_flow.h>
marci@1031: //#include <iter_map.h>
marci@1031: 
marci@1031: using std::cout;
marci@1031: using std::cin;
marci@1031: using std::endl;
marci@1031: 
marci@1031: namespace lemon {
marci@1031:   
marci@1031:   /// \addtogroup misc
marci@1031:   /// @{
marci@1031: 
marci@1031:   /// \brief A partitioned vector with iterable classes.
marci@1031:   ///
marci@1031:   /// This class implements a container in which the data is stored in an 
marci@1031:   /// stl vector, the range is partitioned into sets and each set is 
marci@1031:   /// doubly linked in a list. 
marci@1031:   /// That is, each class is iterable by lemon iterators, and any member of 
marci@1031:   /// the vector can bo moved to an other class.
marci@1031:   template <typename T>
marci@1031:   class IterablePartition {
marci@1031:   protected:
marci@1031:     struct Node {
marci@1031:       T data;
marci@1031:       int prev; //invalid az -1
marci@1031:       int next; 
marci@1031:     };
marci@1031:     std::vector<Node> nodes;
marci@1031:     struct Tip {
marci@1031:       int first;
marci@1031:       int last;
marci@1031:     };
marci@1031:     std::vector<Tip> tips;
marci@1031:   public:
marci@1031:     /// The classes are indexed by integers from \c 0 to \c classNum()-1.
marci@1031:     int classNum() const { return tips.size(); }
marci@1031:     /// This lemon style iterator iterates through a class. 
marci@1031:     class ClassIt;
marci@1031:     /// Constructor. The number of classes is to be given which is fixed 
marci@1031:     /// over the life of the container. 
marci@1031:     /// The partition classes are indexed from 0 to class_num-1. 
marci@1031:     IterablePartition(int class_num) { 
marci@1031:       for (int i=0; i<class_num; ++i) {
marci@1031: 	Tip t;
marci@1031: 	t.first=t.last=-1;
marci@1031: 	tips.push_back(t);
marci@1031:       }
marci@1031:     }
marci@1031:   protected:
marci@1031:     void befuz(ClassIt it, int class_id) {
marci@1031:       if (tips[class_id].first==-1) {
marci@1031: 	if (tips[class_id].last==-1) {
marci@1031: 	  nodes[it.i].prev=nodes[it.i].next=-1;
marci@1031: 	  tips[class_id].first=tips[class_id].last=it.i;
marci@1031: 	}
marci@1031:       } else {
marci@1031: 	nodes[it.i].prev=tips[class_id].last;
marci@1031: 	nodes[it.i].next=-1;
marci@1031: 	nodes[tips[class_id].last].next=it.i;
marci@1031: 	tips[class_id].last=it.i;
marci@1031:       }
marci@1031:     }
marci@1031:     void kifuz(ClassIt it, int class_id) {
marci@1031:       if (tips[class_id].first==it.i) {
marci@1031: 	if (tips[class_id].last==it.i) {
marci@1031: 	  tips[class_id].first=tips[class_id].last=-1;
marci@1031: 	} else {
marci@1031: 	  tips[class_id].first=nodes[it.i].next;
marci@1031: 	  nodes[nodes[it.i].next].prev=-1;
marci@1031: 	}
marci@1031:       } else {
marci@1031: 	if (tips[class_id].last==it.i) {
marci@1031: 	  tips[class_id].last=nodes[it.i].prev;
marci@1031: 	  nodes[nodes[it.i].prev].next=-1;
marci@1031: 	} else {
marci@1031: 	  nodes[nodes[it.i].next].prev=nodes[it.i].prev;
marci@1031: 	  nodes[nodes[it.i].prev].next=nodes[it.i].next;
marci@1031: 	}
marci@1031:       }
marci@1031:     }
marci@1031:   public:
marci@1031:     /// A new element with data \c t is pushed into the vector and into class 
marci@1031:     /// \c class_id.
marci@1031:     ClassIt push_back(const T& t, int class_id) { 
marci@1031:       Node n;
marci@1031:       n.data=t;
marci@1031:       nodes.push_back(n);
marci@1031:       int i=nodes.size()-1;
marci@1031:       befuz(i, class_id);
marci@1031:       return i;
marci@1031:     }
marci@1031:     /// A member is moved to an other class.
marci@1031:     void set(ClassIt it, int old_class_id, int new_class_id) {
marci@1031:       kifuz(it.i, old_class_id);
marci@1031:       befuz(it.i, new_class_id);
marci@1031:     }
marci@1031:     /// Returns the data pointed by \c it.
marci@1031:     T& operator[](ClassIt it) { return nodes[it.i].data; }
marci@1031:     /// Returns the data pointed by \c it.
marci@1031:     const T& operator[](ClassIt it) const { return nodes[it.i].data; }
marci@1031:     ///.
marci@1031:     class ClassIt {
marci@1031:       friend class IterablePartition;
marci@1031:     protected:
marci@1031:       int i;
marci@1031:     public:
marci@1031:       /// Default constructor.
marci@1031:       ClassIt() { }
marci@1031:       /// This constructor constructs an iterator which points
marci@1031:       /// to the member of th container indexed by the integer _i.
marci@1031:       ClassIt(const int& _i) : i(_i) { }
marci@1031:       /// Invalid constructor.
marci@1031:       ClassIt(const Invalid&) : i(-1) { }
marci@1099:       friend bool operator<(const ClassIt& x, const ClassIt& y);
marci@1099:       friend std::ostream& operator<<(std::ostream& os, 
marci@1099: 				      const ClassIt& it);
marci@1031:     };
marci@1099:     friend bool operator<(const ClassIt& x, const ClassIt& y) {
marci@1099:       return (x.i < y.i);
marci@1099:     }
marci@1099:     friend std::ostream& operator<<(std::ostream& os, 
marci@1099: 				    const ClassIt& it) {
marci@1099:       os << it.i;
marci@1099:       return os;
marci@1099:     }
marci@1031:     /// First member of class \c class_id.
marci@1031:     ClassIt& first(ClassIt& it, int class_id) const {
marci@1031:       it.i=tips[class_id].first;
marci@1031:       return it;
marci@1031:     }
marci@1031:     /// Next member.
marci@1031:     ClassIt& next(ClassIt& it) const {
marci@1031:       it.i=nodes[it.i].next;
marci@1031:       return it;
marci@1031:     }
marci@1031:     /// True iff the iterator is valid.
marci@1031:     bool valid(const ClassIt& it) const { return it.i!=-1; }
marci@1031:   };
marci@1031: 
marci@1097: 
marci@1031:   /*! \e
marci@1143:     \todo kellenene uj iterable structure bele, mert ez nem az igazi
marci@1111:     \todo A[x,y]-t cserel. Jobboldal, baloldal csere.
marci@1111:     \todo LEKERDEZESEK!!!
marci@1111:     \todo DOKSI!!!! Doxygen group!!!
marci@1111:     The aim of this class is to give a general surface to different 
marci@1111:     solvers, i.e. it makes possible to write algorithms using LP's, 
marci@1111:     in which the solver can be changed to an other one easily.
marci@1112:     \nosubgrouping
marci@1111:   */
marci@1048:   template <typename _Value>
marci@1031:   class LPSolverBase {
marci@1112:     
marci@1113:     /*! @name Uncategorized functions and types (public members)
marci@1112:     */
marci@1112:     //@{
marci@1031:   public:
marci@1112: 
marci@1112:     //UNCATEGORIZED
marci@1112: 
marci@1031:     /// \e
marci@1048:     typedef _Value Value;
marci@1048:     /// \e
marci@1144:     typedef IterablePartition<int>::ClassIt Row;
marci@1031:     /// \e
marci@1144:     typedef IterablePartition<int>::ClassIt Col;
marci@1074:   public:
marci@1031:     /// \e
marci@1031:     IterablePartition<int> row_iter_map;
marci@1031:     /// \e
marci@1031:     IterablePartition<int> col_iter_map;
marci@1031:     /// \e
marci@1144:     std::vector<Row> int_row_map;
marci@1143:     /// \e
marci@1144:     std::vector<Col> int_col_map;
marci@1143:     /// \e
marci@1074:     const int VALID_CLASS;
marci@1031:     /// \e
marci@1074:     const int INVALID_CLASS;
marci@1104:     /// \e 
marci@1104:     static const _Value INF;
marci@1031:   public:
marci@1031:     /// \e
marci@1031:     LPSolverBase() : row_iter_map(2), 
marci@1031: 		     col_iter_map(2), 
marci@1074: 		     VALID_CLASS(0), INVALID_CLASS(1) { }
marci@1031:     /// \e
marci@1031:     virtual ~LPSolverBase() { }
marci@1112:     //@}
marci@1081: 
marci@1112:     /*! @name Medium level interface (public members)
marci@1112:       These functions appear in the low level and also in the high level 
marci@1112:       interfaces thus these each of these functions have to be implemented 
marci@1112:       only once in the different interfaces.
marci@1112:       This means that these functions have to be reimplemented for all of the 
marci@1112:       different lp solvers. These are basic functions, and have the same 
marci@1112:       parameter lists in the low and high level interfaces. 
marci@1112:     */
marci@1112:     //@{
marci@1112:   public:
marci@1081: 
marci@1112:     //UNCATEGORIZED FUNCTIONS
marci@1112: 
marci@1031:     /// \e
marci@1031:     virtual void setMinimize() = 0;
marci@1031:     /// \e
marci@1031:     virtual void setMaximize() = 0;
marci@1081: 
marci@1112:     //SOLVER FUNCTIONS
marci@1081: 
marci@1112:     /// \e
marci@1112:     virtual void solveSimplex() = 0;
marci@1112:     /// \e
marci@1112:     virtual void solvePrimalSimplex() = 0;
marci@1112:     /// \e
marci@1112:     virtual void solveDualSimplex() = 0;
marci@1112:     /// \e
marci@1112: 
marci@1112:     //SOLUTION RETRIEVING
marci@1112: 
marci@1112:     /// \e
marci@1112:     virtual _Value getObjVal() = 0;
marci@1112: 
marci@1112:     //OTHER FUNCTIONS
marci@1112: 
marci@1112:     /// \e
marci@1112:     virtual int rowNum() const = 0;
marci@1112:     /// \e
marci@1112:     virtual int colNum() const = 0;
marci@1112:     /// \e
marci@1112:     virtual int warmUp() = 0;
marci@1112:     /// \e
marci@1112:     virtual void printWarmUpStatus(int i) = 0;
marci@1112:     /// \e
marci@1112:     virtual int getPrimalStatus() = 0;
marci@1112:     /// \e
marci@1112:     virtual void printPrimalStatus(int i) = 0;
marci@1112:     /// \e
marci@1112:     virtual int getDualStatus() = 0;
marci@1112:     /// \e
marci@1112:     virtual void printDualStatus(int i) = 0;
marci@1144:     /// Returns the status of the slack variable assigned to row \c row.
marci@1144:     virtual int getRowStat(const Row& row) = 0;
marci@1112:     /// \e
marci@1112:     virtual void printRowStatus(int i) = 0;
marci@1144:     /// Returns the status of the variable assigned to column \c col.
marci@1144:     virtual int getColStat(const Col& col) = 0;
marci@1112:     /// \e
marci@1112:     virtual void printColStatus(int i) = 0;
marci@1112: 
marci@1112:     //@}
marci@1112: 
marci@1112:     /*! @name Low level interface (protected members)
marci@1112:       Problem manipulating functions in the low level interface
marci@1112:     */
marci@1112:     //@{
marci@1074:   protected:
marci@1112: 
marci@1112:     //MATRIX MANIPULATING FUNCTIONS
marci@1112: 
marci@1031:     /// \e
marci@1111:     virtual int _addCol() = 0;
marci@1111:     /// \e
marci@1074:     virtual int _addRow() = 0;
marci@1031:     /// \e
marci@1111:     virtual void _eraseCol(int i) = 0;
marci@1111:     /// \e
marci@1111:     virtual void _eraseRow(int i) = 0;
marci@1081:     /// \e
marci@1081:     virtual void _setRowCoeffs(int i, 
marci@1104: 			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
marci@1081:     /// \e
marci@1143:     /// This routine modifies \c coeffs only by the \c push_back method.
marci@1143:     virtual void _getRowCoeffs(int i, 
marci@1143: 			       std::vector<std::pair<int, _Value> >& coeffs) = 0;
marci@1081:     virtual void _setColCoeffs(int i, 
marci@1104: 			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;
marci@1143:     /// \e
marci@1143:     /// This routine modifies \c coeffs only by the \c push_back method.
marci@1143:     virtual void _getColCoeffs(int i, 
marci@1143: 			       std::vector<std::pair<int, _Value> >& coeffs) = 0;
marci@1081:   public:
marci@1081:     /// \e
marci@1081:     enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };
marci@1081:   protected:
marci@1081:     /// \e
marci@1110:     /// The lower bound of a variable (column) have to be given by an 
marci@1110:     /// extended number of type _Value, i.e. a finite number of type 
marci@1110:     /// _Value or -INF.
marci@1110:     virtual void _setColLowerBound(int i, _Value value) = 0;
marci@1110:     /// \e
marci@1111:     /// The lower bound of a variable (column) is an 
marci@1111:     /// extended number of type _Value, i.e. a finite number of type 
marci@1111:     /// _Value or -INF.
marci@1111:     virtual _Value _getColLowerBound(int i) = 0;
marci@1111:     /// \e
marci@1110:     /// The upper bound of a variable (column) have to be given by an 
marci@1110:     /// extended number of type _Value, i.e. a finite number of type 
marci@1110:     /// _Value or INF.
marci@1110:     virtual void _setColUpperBound(int i, _Value value) = 0;
marci@1110:     /// \e
marci@1110:     /// The upper bound of a variable (column) is an 
marci@1110:     /// extended number of type _Value, i.e. a finite number of type 
marci@1110:     /// _Value or INF.
marci@1110:     virtual _Value _getColUpperBound(int i) = 0;
marci@1110:     /// \e
marci@1111:     /// The lower bound of a linear expression (row) have to be given by an 
marci@1111:     /// extended number of type _Value, i.e. a finite number of type 
marci@1111:     /// _Value or -INF.
marci@1111:     virtual void _setRowLowerBound(int i, _Value value) = 0;
marci@1081:     /// \e
marci@1111:     /// The lower bound of a linear expression (row) is an 
marci@1111:     /// extended number of type _Value, i.e. a finite number of type 
marci@1111:     /// _Value or -INF.
marci@1111:     virtual _Value _getRowLowerBound(int i) = 0;
marci@1111:     /// \e
marci@1111:     /// The upper bound of a linear expression (row) have to be given by an 
marci@1111:     /// extended number of type _Value, i.e. a finite number of type 
marci@1111:     /// _Value or INF.
marci@1111:     virtual void _setRowUpperBound(int i, _Value value) = 0;
marci@1111:     /// \e
marci@1111:     /// The upper bound of a linear expression (row) is an 
marci@1111:     /// extended number of type _Value, i.e. a finite number of type 
marci@1111:     /// _Value or INF.
marci@1111:     virtual _Value _getRowUpperBound(int i) = 0;
marci@1081:     /// \e
marci@1081:     virtual void _setObjCoef(int i, _Value obj_coef) = 0;
marci@1081:     /// \e
marci@1081:     virtual _Value _getObjCoef(int i) = 0;
marci@1112:     
marci@1112:     //SOLUTION RETRIEVING
marci@1081: 
marci@1111:     /// \e
marci@1081:     virtual _Value _getPrimal(int i) = 0;
marci@1112:     //@}
marci@1112:     
marci@1112:     /*! @name High level interface (public members)
marci@1112:       Problem manipulating functions in the high level interface
marci@1112:     */
marci@1112:     //@{
marci@1112:   public:
marci@1081: 
marci@1112:     //MATRIX MANIPULATING FUNCTIONS
marci@1081: 
marci@1074:     /// \e
marci@1144:     Col addCol() {
marci@1074:       int i=_addCol();  
marci@1144:       Col col;
marci@1144:       col_iter_map.first(col, INVALID_CLASS);
marci@1144:       if (col_iter_map.valid(col)) { //van hasznalhato hely
marci@1144: 	col_iter_map.set(col, INVALID_CLASS, VALID_CLASS);
marci@1144: 	col_iter_map[col]=i;
marci@1074:       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
marci@1144: 	col=col_iter_map.push_back(i, VALID_CLASS);
marci@1074:       }
marci@1144:       int_col_map.push_back(col);
marci@1144:       return col;
marci@1074:     }
marci@1074:     /// \e
marci@1144:     Row addRow() {
marci@1111:       int i=_addRow();
marci@1144:       Row row;
marci@1144:       row_iter_map.first(row, INVALID_CLASS);
marci@1144:       if (row_iter_map.valid(row)) { //van hasznalhato hely
marci@1144: 	row_iter_map.set(row, INVALID_CLASS, VALID_CLASS);
marci@1144: 	row_iter_map[row]=i;
marci@1111:       } else { //a cucc vegere kell inzertalni mert nincs szabad hely
marci@1144: 	row=row_iter_map.push_back(i, VALID_CLASS);
marci@1031:       }
marci@1144:       int_row_map.push_back(row);
marci@1144:       return row;
marci@1074:     }
marci@1074:     /// \e
marci@1144:     void eraseCol(const Col& col) {
marci@1144:       col_iter_map.set(col, VALID_CLASS, INVALID_CLASS);
marci@1074:       int cols[2];
marci@1144:       cols[1]=col_iter_map[col];
marci@1074:       _eraseCol(cols[1]);
marci@1144:       col_iter_map[col]=0; //glpk specifikus, de kell ez??
marci@1144:       Col it;
marci@1074:       for (col_iter_map.first(it, VALID_CLASS); 
marci@1074: 	   col_iter_map.valid(it); col_iter_map.next(it)) {
marci@1074: 	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
marci@1074:       }
marci@1143:       int_col_map.erase(int_col_map.begin()+cols[1]);
marci@1031:     }
marci@1031:     /// \e
marci@1144:     void eraseRow(const Row& row) {
marci@1144:       row_iter_map.set(row, VALID_CLASS, INVALID_CLASS);
marci@1074:       int rows[2];
marci@1144:       rows[1]=row_iter_map[row];
marci@1074:       _eraseRow(rows[1]);
marci@1144:       row_iter_map[row]=0; //glpk specifikus, de kell ez??
marci@1144:       Row it;
marci@1074:       for (row_iter_map.first(it, VALID_CLASS); 
marci@1074: 	   row_iter_map.valid(it); row_iter_map.next(it)) {
marci@1074: 	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
marci@1074:       }
marci@1143:       int_row_map.erase(int_row_map.begin()+rows[1]);
marci@1074:     }
marci@1031:     /// \e
marci@1144:     void setColLowerBound(Col col, _Value lo) {
marci@1144:       _setColLowerBound(col_iter_map[col], lo);
marci@1111:     }
marci@1111:     /// \e
marci@1144:     _Value getColLowerBound(Col col) {
marci@1144:       return _getColLowerBound(col_iter_map[col]);
marci@1111:     }
marci@1111:     /// \e
marci@1144:     void setColUpperBound(Col col, _Value up) {
marci@1144:       _setColUpperBound(col_iter_map[col], up);
marci@1110:     }
marci@1110:     /// \e
marci@1144:     _Value getColUpperBound(Col col) {      
marci@1144:       return _getColUpperBound(col_iter_map[col]);
marci@1111:     }
marci@1111:     /// \e
marci@1144:     void setRowLowerBound(Row row, _Value lo) {
marci@1144:       _setRowLowerBound(row_iter_map[row], lo);
marci@1110:     }
marci@1110:     /// \e
marci@1144:     _Value getRowLowerBound(Row row) {
marci@1144:       return _getRowLowerBound(row_iter_map[row]);
marci@1110:     }
marci@1110:     /// \e
marci@1144:     void setRowUpperBound(Row row, _Value up) {
marci@1144:       _setRowUpperBound(row_iter_map[row], up);
marci@1081:     }
marci@1031:     /// \e
marci@1144:     _Value getRowUpperBound(Row row) {      
marci@1144:       return _getRowUpperBound(row_iter_map[row]);
marci@1111:     }
marci@1111:     /// \e
marci@1144:     void setObjCoef(const Col& col, _Value obj_coef) {
marci@1144:       _setObjCoef(col_iter_map[col], obj_coef);
marci@1111:     }
marci@1111:     /// \e
marci@1144:     _Value getObjCoef(const Col& col) {
marci@1144:       return _getObjCoef(col_iter_map[col]);
marci@1081:     }
marci@1081: 
marci@1112:     //SOLUTION RETRIEVING FUNCTIONS
marci@1112: 
marci@1112:     /// \e
marci@1144:     _Value getPrimal(const Col& col) {
marci@1144:       return _getPrimal(col_iter_map[col]);
marci@1112:     }    
marci@1112: 
marci@1112:     //@}
marci@1112: 
marci@1112:     /*! @name User friend interface
marci@1112:       Problem manipulating functions in the user friend interface
marci@1112:     */
marci@1112:     //@{
marci@1112: 
marci@1112:     //EXPRESSION TYPES
marci@1099: 
marci@1099:     /// \e
marci@1144:     typedef Expr<Col, _Value> Expression;
marci@1099:     /// \e
marci@1144:     typedef Expr<Row, _Value> DualExpression;
marci@1144:     /// \e
marci@1144:     typedef Constr<Col, _Value> Constraint;
marci@1112: 
marci@1112:     //MATRIX MANIPULATING FUNCTIONS
marci@1112: 
marci@1099:     /// \e
marci@1144:     void setRowCoeffs(Row row, const Expression& expr) {
marci@1099:       std::vector<std::pair<int, _Value> > row_coeffs;
marci@1099:       for(typename Expression::Data::const_iterator i=expr.data.begin(); 
marci@1099: 	  i!=expr.data.end(); ++i) {
marci@1099: 	row_coeffs.push_back(std::make_pair
marci@1099: 			     (col_iter_map[(*i).first], (*i).second));
marci@1099:       }
marci@1144:       _setRowCoeffs(row_iter_map[row], row_coeffs);
marci@1144:     }
marci@1144:     /// \e 
marci@1144:     void setRow(Row row, const Constraint& constr) {
marci@1144:       setRowCoeffs(row, constr.expr);
marci@1144:       setRowLowerBound(row, constr.lo);
marci@1144:       setRowUpperBound(row, constr.up);
marci@1144:     }
marci@1144:     /// \e 
marci@1144:     Row addRow(const Constraint& constr) {
marci@1144:       Row row=addRow();
marci@1144:       setRowCoeffs(row, constr.expr);
marci@1144:       setRowLowerBound(row, constr.lo);
marci@1144:       setRowUpperBound(row, constr.up);
marci@1144:       return row;
marci@1099:     }
marci@1099:     /// \e
marci@1143:     /// This routine modifies \c expr by only adding to it.
marci@1144:     void getRowCoeffs(Row row, Expression& expr) {
marci@1143:       std::vector<std::pair<int, _Value> > row_coeffs;
marci@1144:       _getRowCoeffs(row_iter_map[row], row_coeffs);
marci@1143:       for(typename std::vector<std::pair<int, _Value> >::const_iterator 
marci@1143:  	    i=row_coeffs.begin(); i!=row_coeffs.end(); ++i) {
marci@1143:  	expr+= (*i).second*int_col_map[(*i).first];
marci@1143:       }
marci@1143:     }
marci@1143:     /// \e
marci@1144:     void setColCoeffs(Col col, const DualExpression& expr) {
marci@1099:       std::vector<std::pair<int, _Value> > col_coeffs;
marci@1099:       for(typename DualExpression::Data::const_iterator i=expr.data.begin(); 
marci@1099: 	  i!=expr.data.end(); ++i) {
marci@1099: 	col_coeffs.push_back(std::make_pair
marci@1099: 			     (row_iter_map[(*i).first], (*i).second));
marci@1099:       }
marci@1144:       _setColCoeffs(col_iter_map[col], col_coeffs);
marci@1099:     }
marci@1099:     /// \e
marci@1143:     /// This routine modifies \c expr by only adding to it.
marci@1144:     void getColCoeffs(Col col, DualExpression& expr) {
marci@1143:       std::vector<std::pair<int, _Value> > col_coeffs;
marci@1144:       _getColCoeffs(col_iter_map[col], col_coeffs);
marci@1143:       for(typename std::vector<std::pair<int, _Value> >::const_iterator 
marci@1143:  	    i=col_coeffs.begin(); i!=col_coeffs.end(); ++i) {
marci@1143:  	expr+= (*i).second*int_row_map[(*i).first];
marci@1143:       }
marci@1143:     }
marci@1143:     /// \e
marci@1143:     /// \bug ez igy nem jo
marci@1099:     void setObjCoeffs(const Expression& expr) {
marci@1099:       for(typename Expression::Data::const_iterator i=expr.data.begin(); 
marci@1099: 	  i!=expr.data.end(); ++i) {
marci@1099: 	setObjCoef((*i).first, (*i).second);
marci@1099:       }
marci@1099:     }
marci@1143:     /// \e
marci@1143:     /// This routine modifies \c expr by only adding to it.
marci@1143:     void getObjCoeffs(Expression& expr) {
marci@1143:       /// FIXME not yet implemented
marci@1143:     }
marci@1112:     //@}
marci@1031:   };
marci@1031:   
marci@1104:   template <typename _Value>
marci@1104:   const _Value LPSolverBase<_Value>::INF=std::numeric_limits<_Value>::infinity();
marci@1104: 
marci@1048: 
marci@1111:   /// \brief Wrapper for GLPK solver
marci@1031:   /// 
marci@1111:   /// This class implements a lemon wrapper for GLPK.
marci@1081:   class LPGLPK : public LPSolverBase<double> {
marci@1031:   public:
marci@1048:     typedef LPSolverBase<double> Parent;
marci@1031: 
marci@1031:   public:
marci@1031:     /// \e
marci@1031:     LPX* lp;
marci@1031: 
marci@1031:   public:
marci@1031:     /// \e
marci@1081:     LPGLPK() : Parent(), 
marci@1031: 			lp(lpx_create_prob()) {
marci@1144:       int_row_map.push_back(Row());
marci@1144:       int_col_map.push_back(Col());
marci@1031:       lpx_set_int_parm(lp, LPX_K_DUAL, 1);
marci@1031:     }
marci@1031:     /// \e
marci@1081:     ~LPGLPK() {
marci@1031:       lpx_delete_prob(lp);
marci@1031:     }
marci@1081: 
marci@1081:     //MATRIX INDEPEDENT MANIPULATING FUNCTIONS
marci@1081: 
marci@1031:     /// \e
marci@1031:     void setMinimize() { 
marci@1031:       lpx_set_obj_dir(lp, LPX_MIN);
marci@1031:     }
marci@1031:     /// \e
marci@1031:     void setMaximize() { 
marci@1031:       lpx_set_obj_dir(lp, LPX_MAX);
marci@1031:     }
marci@1081: 
marci@1081:     //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS
marci@1081: 
marci@1074:   protected:
marci@1031:     /// \e
marci@1074:     int _addCol() { 
marci@1110:       int i=lpx_add_cols(lp, 1);
marci@1110:       _setColLowerBound(i, -INF);
marci@1110:       _setColUpperBound(i, INF);
marci@1110:       return i;
marci@1031:     }
marci@1031:     /// \e
marci@1074:     int _addRow() { 
marci@1110:       int i=lpx_add_rows(lp, 1);
marci@1110:       return i;
marci@1074:     }
marci@1074:     /// \e
marci@1081:     virtual void _setRowCoeffs(int i, 
marci@1104: 			       const std::vector<std::pair<int, double> >& coeffs) {
marci@1074:       int mem_length=1+colNum();
marci@1074:       int* indices = new int[mem_length];
marci@1074:       double* doubles = new double[mem_length];
marci@1074:       int length=0;
marci@1074:       for (std::vector<std::pair<int, double> >::
marci@1074: 	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
marci@1074: 	++length;
marci@1074: 	indices[length]=it->first;
marci@1074: 	doubles[length]=it->second;
marci@1031:       }
marci@1074:       lpx_set_mat_row(lp, i, length, indices, doubles);
marci@1074:       delete [] indices;
marci@1074:       delete [] doubles;
marci@1031:     }
marci@1074:     /// \e
marci@1143:     virtual void _getRowCoeffs(int i, 
marci@1143: 			       std::vector<std::pair<int, double> >& coeffs) {
marci@1143:       int mem_length=1+colNum();
marci@1143:       int* indices = new int[mem_length];
marci@1143:       double* doubles = new double[mem_length];
marci@1143:       int length=lpx_get_mat_row(lp, i, indices, doubles);
marci@1143:       for (int i=1; i<=length; ++i) {
marci@1143: 	coeffs.push_back(std::make_pair(indices[i], doubles[i]));
marci@1143:       }
marci@1143:       delete [] indices;
marci@1143:       delete [] doubles;
marci@1143:     }
marci@1143:     /// \e
marci@1081:     virtual void _setColCoeffs(int i, 
marci@1104: 			       const std::vector<std::pair<int, double> >& coeffs) {
marci@1074:       int mem_length=1+rowNum();
marci@1074:       int* indices = new int[mem_length];
marci@1074:       double* doubles = new double[mem_length];
marci@1074:       int length=0;
marci@1074:       for (std::vector<std::pair<int, double> >::
marci@1074: 	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
marci@1074: 	++length;
marci@1074: 	indices[length]=it->first;
marci@1074: 	doubles[length]=it->second;
marci@1074:       }
marci@1074:       lpx_set_mat_col(lp, i, length, indices, doubles);
marci@1074:       delete [] indices;
marci@1074:       delete [] doubles;
marci@1031:     }
marci@1031:     /// \e
marci@1143:     virtual void _getColCoeffs(int i, 
marci@1143: 			       std::vector<std::pair<int, double> >& coeffs) {
marci@1143:       int mem_length=1+rowNum();
marci@1143:       int* indices = new int[mem_length];
marci@1143:       double* doubles = new double[mem_length];
marci@1143:       int length=lpx_get_mat_col(lp, i, indices, doubles);
marci@1143:       for (int i=1; i<=length; ++i) {
marci@1143: 	coeffs.push_back(std::make_pair(indices[i], doubles[i]));
marci@1143:       }
marci@1143:       delete [] indices;
marci@1143:       delete [] doubles;
marci@1143:     }
marci@1143:     /// \e
marci@1074:     virtual void _eraseCol(int i) {
marci@1031:       int cols[2];
marci@1074:       cols[1]=i;
marci@1031:       lpx_del_cols(lp, 1, cols);
marci@1031:     }
marci@1074:     virtual void _eraseRow(int i) {
marci@1031:       int rows[2];
marci@1074:       rows[1]=i;
marci@1031:       lpx_del_rows(lp, 1, rows);
marci@1031:     }
marci@1110:     virtual void _setColLowerBound(int i, double lo) {
marci@1110:       if (lo==INF) {
marci@1110: 	//FIXME error
marci@1110:       }
marci@1110:       int b=lpx_get_col_type(lp, i);
marci@1110:       double up=lpx_get_col_ub(lp, i);	
marci@1110:       if (lo==-INF) {
marci@1110: 	switch (b) {
marci@1110: 	case LPX_FR:
marci@1110: 	case LPX_LO:
marci@1110: 	  lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
marci@1110: 	  break;
marci@1110: 	case LPX_UP:
marci@1110: 	  break;
marci@1110: 	case LPX_DB:
marci@1110: 	case LPX_FX:
marci@1110: 	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
marci@1110: 	  break;
marci@1110: 	default: ;
marci@1110: 	  //FIXME error
marci@1110: 	}
marci@1110:       } else {
marci@1110: 	switch (b) {
marci@1110: 	case LPX_FR:
marci@1110: 	case LPX_LO:
marci@1110: 	  lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
marci@1110: 	  break;
marci@1110: 	case LPX_UP:	  
marci@1110: 	case LPX_DB:
marci@1110: 	case LPX_FX:
marci@1110: 	  if (lo==up) 
marci@1110: 	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
marci@1110: 	  else 
marci@1110: 	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
marci@1110: 	  break;
marci@1110: 	default: ;
marci@1110: 	  //FIXME error
marci@1110: 	}
marci@1110:       }
marci@1110:     }
marci@1111:     virtual double _getColLowerBound(int i) {
marci@1111:       int b=lpx_get_col_type(lp, i);
marci@1111:       switch (b) {
marci@1111:       case LPX_FR:
marci@1111: 	return -INF;
marci@1111:       case LPX_LO:
marci@1111: 	return lpx_get_col_lb(lp, i);
marci@1111:       case LPX_UP:
marci@1111: 	return -INF;
marci@1111:       case LPX_DB:
marci@1111:       case LPX_FX:
marci@1111: 	return lpx_get_col_lb(lp, i);
marci@1111:       default: ;
marci@1111: 	//FIXME error
marci@1111: 	return 0.0;
marci@1111:       }
marci@1111:     }
marci@1110:     virtual void _setColUpperBound(int i, double up) {
marci@1110:       if (up==-INF) {
marci@1110: 	//FIXME error
marci@1110:       }
marci@1110:       int b=lpx_get_col_type(lp, i);
marci@1110:       double lo=lpx_get_col_lb(lp, i);
marci@1110:       if (up==INF) {
marci@1110: 	switch (b) {
marci@1110: 	case LPX_FR:
marci@1110: 	case LPX_LO:
marci@1110: 	  break;
marci@1110: 	case LPX_UP:
marci@1110: 	  lpx_set_col_bnds(lp, i, LPX_FR, lo, up);
marci@1110: 	  break;
marci@1110: 	case LPX_DB:
marci@1110: 	case LPX_FX:
marci@1110: 	  lpx_set_col_bnds(lp, i, LPX_LO, lo, up);
marci@1110: 	  break;
marci@1110: 	default: ;
marci@1110: 	  //FIXME error
marci@1110: 	}
marci@1110:       } else {
marci@1110: 	switch (b) {
marci@1110: 	case LPX_FR:
marci@1110: 	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
marci@1110: 	case LPX_LO:
marci@1110: 	  if (lo==up) 
marci@1110: 	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
marci@1110: 	  else
marci@1110: 	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
marci@1110: 	  break;
marci@1110: 	case LPX_UP:
marci@1110: 	  lpx_set_col_bnds(lp, i, LPX_UP, lo, up);
marci@1110: 	  break;
marci@1110: 	case LPX_DB:
marci@1110: 	case LPX_FX:
marci@1110: 	  if (lo==up) 
marci@1110: 	    lpx_set_col_bnds(lp, i, LPX_FX, lo, up);
marci@1110: 	  else 
marci@1110: 	    lpx_set_col_bnds(lp, i, LPX_DB, lo, up);
marci@1110: 	  break;
marci@1110: 	default: ;
marci@1110: 	  //FIXME error
marci@1110: 	}
marci@1110:       }
marci@1110:     }
marci@1110:     virtual double _getColUpperBound(int i) {
marci@1110:       int b=lpx_get_col_type(lp, i);
marci@1110:       switch (b) {
marci@1110:       case LPX_FR:
marci@1110:       case LPX_LO:
marci@1110: 	return INF;
marci@1110:       case LPX_UP:
marci@1110:       case LPX_DB:
marci@1110:       case LPX_FX:
marci@1110: 	return lpx_get_col_ub(lp, i);
marci@1110:       default: ;
marci@1110: 	//FIXME error
marci@1110: 	return 0.0;
marci@1110:       }
marci@1110:     }
marci@1111:     virtual void _setRowLowerBound(int i, double lo) {
marci@1111:       if (lo==INF) {
marci@1111: 	//FIXME error
marci@1081:       }
marci@1111:       int b=lpx_get_row_type(lp, i);
marci@1111:       double up=lpx_get_row_ub(lp, i);	
marci@1111:       if (lo==-INF) {
marci@1111: 	switch (b) {
marci@1111: 	case LPX_FR:
marci@1111: 	case LPX_LO:
marci@1111: 	  lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
marci@1111: 	  break;
marci@1111: 	case LPX_UP:
marci@1111: 	  break;
marci@1111: 	case LPX_DB:
marci@1111: 	case LPX_FX:
marci@1111: 	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
marci@1111: 	  break;
marci@1111: 	default: ;
marci@1111: 	  //FIXME error
marci@1111: 	}
marci@1111:       } else {
marci@1111: 	switch (b) {
marci@1111: 	case LPX_FR:
marci@1111: 	case LPX_LO:
marci@1111: 	  lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
marci@1111: 	  break;
marci@1111: 	case LPX_UP:	  
marci@1111: 	case LPX_DB:
marci@1111: 	case LPX_FX:
marci@1111: 	  if (lo==up) 
marci@1111: 	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
marci@1111: 	  else 
marci@1111: 	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
marci@1111: 	  break;
marci@1111: 	default: ;
marci@1111: 	  //FIXME error
marci@1111: 	}
marci@1081:       }
marci@1111:     }
marci@1111:     virtual double _getRowLowerBound(int i) {
marci@1111:       int b=lpx_get_row_type(lp, i);
marci@1111:       switch (b) {
marci@1111:       case LPX_FR:
marci@1111: 	return -INF;
marci@1111:       case LPX_LO:
marci@1111: 	return lpx_get_row_lb(lp, i);
marci@1111:       case LPX_UP:
marci@1111: 	return -INF;
marci@1111:       case LPX_DB:
marci@1111:       case LPX_FX:
marci@1111: 	return lpx_get_row_lb(lp, i);
marci@1111:       default: ;
marci@1111: 	//FIXME error
marci@1111: 	return 0.0;
marci@1111:       }
marci@1111:     }
marci@1111:     virtual void _setRowUpperBound(int i, double up) {
marci@1111:       if (up==-INF) {
marci@1111: 	//FIXME error
marci@1111:       }
marci@1111:       int b=lpx_get_row_type(lp, i);
marci@1111:       double lo=lpx_get_row_lb(lp, i);
marci@1111:       if (up==INF) {
marci@1111: 	switch (b) {
marci@1111: 	case LPX_FR:
marci@1111: 	case LPX_LO:
marci@1111: 	  break;
marci@1111: 	case LPX_UP:
marci@1111: 	  lpx_set_row_bnds(lp, i, LPX_FR, lo, up);
marci@1111: 	  break;
marci@1111: 	case LPX_DB:
marci@1111: 	case LPX_FX:
marci@1111: 	  lpx_set_row_bnds(lp, i, LPX_LO, lo, up);
marci@1111: 	  break;
marci@1111: 	default: ;
marci@1111: 	  //FIXME error
marci@1111: 	}
marci@1111:       } else {
marci@1111: 	switch (b) {
marci@1111: 	case LPX_FR:
marci@1111: 	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
marci@1111: 	case LPX_LO:
marci@1111: 	  if (lo==up) 
marci@1111: 	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
marci@1111: 	  else
marci@1111: 	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
marci@1111: 	  break;
marci@1111: 	case LPX_UP:
marci@1111: 	  lpx_set_row_bnds(lp, i, LPX_UP, lo, up);
marci@1111: 	  break;
marci@1111: 	case LPX_DB:
marci@1111: 	case LPX_FX:
marci@1111: 	  if (lo==up) 
marci@1111: 	    lpx_set_row_bnds(lp, i, LPX_FX, lo, up);
marci@1111: 	  else 
marci@1111: 	    lpx_set_row_bnds(lp, i, LPX_DB, lo, up);
marci@1111: 	  break;
marci@1111: 	default: ;
marci@1111: 	  //FIXME error
marci@1111: 	}
marci@1111:       }
marci@1111:     }
marci@1111:     virtual double _getRowUpperBound(int i) {
marci@1111:       int b=lpx_get_row_type(lp, i);
marci@1111:       switch (b) {
marci@1111:       case LPX_FR:
marci@1111:       case LPX_LO:
marci@1111: 	return INF;
marci@1111:       case LPX_UP:
marci@1111:       case LPX_DB:
marci@1111:       case LPX_FX:
marci@1111: 	return lpx_get_row_ub(lp, i);
marci@1111:       default: ;
marci@1111: 	//FIXME error
marci@1111: 	return 0.0;
marci@1111:       }
marci@1111:     }
marci@1031:     /// \e
marci@1081:     virtual double _getObjCoef(int i) { 
marci@1081:       return lpx_get_obj_coef(lp, i);
marci@1031:     }
marci@1031:     /// \e
marci@1081:     virtual void _setObjCoef(int i, double obj_coef) { 
marci@1081:       lpx_set_obj_coef(lp, i, obj_coef);
marci@1031:     }
marci@1081:   public:
marci@1031:     /// \e
marci@1031:     void solveSimplex() { lpx_simplex(lp); }
marci@1031:     /// \e
marci@1031:     void solvePrimalSimplex() { lpx_simplex(lp); }
marci@1031:     /// \e
marci@1031:     void solveDualSimplex() { lpx_simplex(lp); }
marci@1031:     /// \e
marci@1081:   protected:
marci@1081:     virtual double _getPrimal(int i) {
marci@1081:       return lpx_get_col_prim(lp, i);
marci@1031:     }
marci@1081:   public:
marci@1031:     /// \e
marci@1031:     double getObjVal() { return lpx_get_obj_val(lp); }
marci@1031:     /// \e
marci@1031:     int rowNum() const { return lpx_get_num_rows(lp); }
marci@1031:     /// \e
marci@1031:     int colNum() const { return lpx_get_num_cols(lp); }
marci@1031:     /// \e
marci@1031:     int warmUp() { return lpx_warm_up(lp); }
marci@1031:     /// \e
marci@1031:     void printWarmUpStatus(int i) {
marci@1031:       switch (i) {
marci@1031:       case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
marci@1031:       case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;	
marci@1031:       case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
marci@1031:       case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
marci@1031:       }
marci@1031:     }
marci@1031:     /// \e
marci@1031:     int getPrimalStatus() { return lpx_get_prim_stat(lp); }
marci@1031:     /// \e
marci@1031:     void printPrimalStatus(int i) {
marci@1031:       switch (i) {
marci@1031:       case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
marci@1031:       case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;	
marci@1031:       case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
marci@1031:       case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
marci@1031:       }
marci@1031:     }
marci@1031:     /// \e
marci@1031:     int getDualStatus() { return lpx_get_dual_stat(lp); }
marci@1031:     /// \e
marci@1031:     void printDualStatus(int i) {
marci@1031:       switch (i) {
marci@1031:       case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
marci@1031:       case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;	
marci@1031:       case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
marci@1031:       case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
marci@1031:       }
marci@1031:     }
marci@1144:     /// Returns the status of the slack variable assigned to row \c row.
marci@1144:     int getRowStat(const Row& row) { 
marci@1144:       return lpx_get_row_stat(lp, row_iter_map[row]); 
marci@1031:     }
marci@1031:     /// \e
marci@1031:     void printRowStatus(int i) {
marci@1031:       switch (i) {
marci@1031:       case LPX_BS: cout << "LPX_BS" << endl; break;
marci@1031:       case LPX_NL: cout << "LPX_NL" << endl; break;	
marci@1031:       case LPX_NU: cout << "LPX_NU" << endl; break;
marci@1031:       case LPX_NF: cout << "LPX_NF" << endl; break;
marci@1031:       case LPX_NS: cout << "LPX_NS" << endl; break;
marci@1031:       }
marci@1031:     }
marci@1144:     /// Returns the status of the variable assigned to column \c col.
marci@1144:     int getColStat(const Col& col) { 
marci@1144:       return lpx_get_col_stat(lp, col_iter_map[col]); 
marci@1031:     }
marci@1031:     /// \e
marci@1031:     void printColStatus(int i) {
marci@1031:       switch (i) {
marci@1031:       case LPX_BS: cout << "LPX_BS" << endl; break;
marci@1031:       case LPX_NL: cout << "LPX_NL" << endl; break;	
marci@1031:       case LPX_NU: cout << "LPX_NU" << endl; break;
marci@1031:       case LPX_NF: cout << "LPX_NF" << endl; break;
marci@1031:       case LPX_NS: cout << "LPX_NS" << endl; break;
marci@1031:       }
marci@1031:     }
marci@1031:   };
marci@1031:   
marci@1031:   /// @}
marci@1031: 
marci@1031: } //namespace lemon
marci@1031: 
marci@1031: #endif //LEMON_LP_SOLVER_WRAPPER_H