lemon/lp_base.h
author hegyi
Tue, 29 Nov 2005 19:31:58 +0000
changeset 1837 8dd6160ff699
parent 1787 932b8490caf0
child 1840 173b53b28d7c
permissions -rw-r--r--
Structure of GUI is now more clear-cut than before.
     1 /* -*- C++ -*-
     2  * lemon/lp_base.h - Part of LEMON, a generic C++ optimization library
     3  *
     4  * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     5  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     6  *
     7  * Permission to use, modify and distribute this software is granted
     8  * provided that this copyright notice appears in all copies. For
     9  * precise terms see the accompanying LICENSE file.
    10  *
    11  * This software is provided "AS IS" with no warranty of any kind,
    12  * express or implied, and with no claim as to its suitability for any
    13  * purpose.
    14  *
    15  */
    16 
    17 #ifndef LEMON_LP_BASE_H
    18 #define LEMON_LP_BASE_H
    19 
    20 #include<vector>
    21 #include<map>
    22 #include<limits>
    23 #include<cmath>
    24 
    25 #include<lemon/utility.h>
    26 #include<lemon/error.h>
    27 #include<lemon/invalid.h>
    28 
    29 ///\file
    30 ///\brief The interface of the LP solver interface.
    31 ///\ingroup gen_opt_group
    32 namespace lemon {
    33   
    34   ///Internal data structure to convert floating id's to fix one's
    35     
    36   ///\todo This might be implemented to be also usable in other places.
    37   class _FixId 
    38   {
    39   protected:
    40     std::vector<int> index;
    41     std::vector<int> cross;
    42     int first_free;
    43   public:
    44     _FixId() : first_free(-1) {};
    45     ///Convert a floating id to a fix one
    46 
    47     ///\param n is a floating id
    48     ///\return the corresponding fix id
    49     int fixId(int n) const {return cross[n];}
    50     ///Convert a fix id to a floating one
    51 
    52     ///\param n is a fix id
    53     ///\return the corresponding floating id
    54     int floatingId(int n) const { return index[n];}
    55     ///Add a new floating id.
    56 
    57     ///\param n is a floating id
    58     ///\return the fix id of the new value
    59     ///\todo Multiple additions should also be handled.
    60     int insert(int n)
    61     {
    62       if(n>=int(cross.size())) {
    63 	cross.resize(n+1);
    64 	if(first_free==-1) {
    65 	  cross[n]=index.size();
    66 	  index.push_back(n);
    67 	}
    68 	else {
    69 	  cross[n]=first_free;
    70 	  int next=index[first_free];
    71 	  index[first_free]=n;
    72 	  first_free=next;
    73 	}
    74 	return cross[n];
    75       }
    76       ///\todo Create an own exception type.
    77       else throw LogicError(); //floatingId-s must form a continuous range;
    78     }
    79     ///Remove a fix id.
    80 
    81     ///\param n is a fix id
    82     ///
    83     void erase(int n) 
    84     {
    85       int fl=index[n];
    86       index[n]=first_free;
    87       first_free=n;
    88       for(int i=fl+1;i<int(cross.size());++i) {
    89 	cross[i-1]=cross[i];
    90 	index[cross[i]]--;
    91       }
    92       cross.pop_back();
    93     }
    94     ///An upper bound on the largest fix id.
    95 
    96     ///\todo Do we need this?
    97     ///
    98     std::size_t maxFixId() { return cross.size()-1; }
    99   
   100   };
   101     
   102   ///Common base class for LP solvers
   103   
   104   ///\todo Much more docs
   105   ///\ingroup gen_opt_group
   106   class LpSolverBase {
   107 
   108   public:
   109 
   110     ///Possible outcomes of an LP solving procedure
   111     enum SolveExitStatus {
   112       ///This means that the problem has been successfully solved: either
   113       ///an optimal solution has been found or infeasibility/unboundedness
   114       ///has been proved.
   115       SOLVED = 0,
   116       ///Any other case (including the case when some user specified limit has been exceeded)
   117       UNSOLVED = 1
   118     };
   119       
   120       ///\e
   121     enum SolutionStatus {
   122       ///Feasible solution has'n been found (but may exist).
   123 
   124       ///\todo NOTFOUND might be a better name.
   125       ///
   126       UNDEFINED = 0,
   127       ///The problem has no feasible solution
   128       INFEASIBLE = 1,
   129       ///Feasible solution found
   130       FEASIBLE = 2,
   131       ///Optimal solution exists and found
   132       OPTIMAL = 3,
   133       ///The cost function is unbounded
   134 
   135       ///\todo Give a feasible solution and an infinite ray (and the
   136       ///corresponding bases)
   137       INFINITE = 4
   138     };
   139 
   140     ///\e The type of the investigated LP problem
   141     enum ProblemTypes {
   142       ///Primal-dual feasible
   143       PRIMAL_DUAL_FEASIBLE = 0,
   144       ///Primal feasible dual infeasible
   145       PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
   146       ///Primal infeasible dual feasible
   147       PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
   148       ///Primal-dual infeasible
   149       PRIMAL_DUAL_INFEASIBLE = 3,
   150       ///Could not determine so far
   151       UNKNOWN = 4
   152     };
   153 
   154     ///The floating point type used by the solver
   155     typedef double Value;
   156     ///The infinity constant
   157     static const Value INF;
   158     ///The not a number constant
   159     static const Value NaN;
   160     
   161     ///Refer to a column of the LP.
   162 
   163     ///This type is used to refer to a column of the LP.
   164     ///
   165     ///Its value remains valid and correct even after the addition or erase of
   166     ///other columns.
   167     ///
   168     ///\todo Document what can one do with a Col (INVALID, comparing,
   169     ///it is similar to Node/Edge)
   170     class Col {
   171     protected:
   172       int id;
   173       friend class LpSolverBase;
   174     public:
   175       typedef Value ExprValue;
   176       typedef True LpSolverCol;
   177       Col() {}
   178       Col(const Invalid&) : id(-1) {}
   179       bool operator<(Col c) const  {return id<c.id;}
   180       bool operator==(Col c) const  {return id==c.id;}
   181       bool operator!=(Col c) const  {return id==c.id;}
   182     };
   183 
   184     ///Refer to a row of the LP.
   185 
   186     ///This type is used to refer to a row of the LP.
   187     ///
   188     ///Its value remains valid and correct even after the addition or erase of
   189     ///other rows.
   190     ///
   191     ///\todo Document what can one do with a Row (INVALID, comparing,
   192     ///it is similar to Node/Edge)
   193     class Row {
   194     protected:
   195       int id;
   196       friend class LpSolverBase;
   197     public:
   198       typedef Value ExprValue;
   199       typedef True LpSolverRow;
   200       Row() {}
   201       Row(const Invalid&) : id(-1) {}
   202 
   203       bool operator<(Row c) const  {return id<c.id;}
   204       bool operator==(Row c) const  {return id==c.id;}
   205       bool operator!=(Row c) const  {return id==c.id;} 
   206    };
   207     
   208     ///Linear expression of variables and a constant component
   209     
   210     ///This data structure strores a linear expression of the variables
   211     ///(\ref Col "Col"s) and also has a constant component.
   212     ///
   213     ///There are several ways to access and modify the contents of this
   214     ///container.
   215     ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
   216     ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can
   217     ///read and modify the coefficients like
   218     ///these.
   219     ///\code
   220     ///e[v]=5;
   221     ///e[v]+=12;
   222     ///e.erase(v);
   223     ///\endcode
   224     ///or you can also iterate through its elements.
   225     ///\code
   226     ///double s=0;
   227     ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
   228     ///  s+=i->second;
   229     ///\endcode
   230     ///(This code computes the sum of all coefficients).
   231     ///- Numbers (<tt>double</tt>'s)
   232     ///and variables (\ref Col "Col"s) directly convert to an
   233     ///\ref Expr and the usual linear operations are defined so  
   234     ///\code
   235     ///v+w
   236     ///2*v-3.12*(v-w/2)+2
   237     ///v*2.1+(3*v+(v*12+w+6)*3)/2
   238     ///\endcode
   239     ///are valid \ref Expr "Expr"essions.
   240     ///The usual assignment operations are also defined.
   241     ///\code
   242     ///e=v+w;
   243     ///e+=2*v-3.12*(v-w/2)+2;
   244     ///e*=3.4;
   245     ///e/=5;
   246     ///\endcode
   247     ///- The constant member can be set and read by \ref constComp()
   248     ///\code
   249     ///e.constComp()=12;
   250     ///double c=e.constComp();
   251     ///\endcode
   252     ///
   253     ///\note \ref clear() not only sets all coefficients to 0 but also
   254     ///clears the constant components.
   255     ///
   256     ///\sa Constr
   257     ///
   258     class Expr : public std::map<Col,Value>
   259     {
   260     public:
   261       typedef LpSolverBase::Col Key; 
   262       typedef LpSolverBase::Value Value;
   263       
   264     protected:
   265       typedef std::map<Col,Value> Base;
   266       
   267       Value const_comp;
   268   public:
   269       typedef True IsLinExpression;
   270       ///\e
   271       Expr() : Base(), const_comp(0) { }
   272       ///\e
   273       Expr(const Key &v) : const_comp(0) {
   274 	Base::insert(std::make_pair(v, 1));
   275       }
   276       ///\e
   277       Expr(const Value &v) : const_comp(v) {}
   278       ///\e
   279       void set(const Key &v,const Value &c) {
   280 	Base::insert(std::make_pair(v, c));
   281       }
   282       ///\e
   283       Value &constComp() { return const_comp; }
   284       ///\e
   285       const Value &constComp() const { return const_comp; }
   286       
   287       ///Removes the components with zero coefficient.
   288       void simplify() {
   289 	for (Base::iterator i=Base::begin(); i!=Base::end();) {
   290 	  Base::iterator j=i;
   291 	  ++j;
   292 	  if ((*i).second==0) Base::erase(i);
   293 	  j=i;
   294 	}
   295       }
   296 
   297       ///Removes the coefficients closer to zero than \c tolerance.
   298       void simplify(double &tolerance) {
   299 	for (Base::iterator i=Base::begin(); i!=Base::end();) {
   300 	  Base::iterator j=i;
   301 	  ++j;
   302 	  if (std::fabs((*i).second)<tolerance) Base::erase(i);
   303 	  j=i;
   304 	}
   305       }
   306 
   307       ///Sets all coefficients and the constant component to 0.
   308       void clear() {
   309 	Base::clear();
   310 	const_comp=0;
   311       }
   312 
   313       ///\e
   314       Expr &operator+=(const Expr &e) {
   315 	for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
   316 	  (*this)[j->first]+=j->second;
   317 	const_comp+=e.const_comp;
   318 	return *this;
   319       }
   320       ///\e
   321       Expr &operator-=(const Expr &e) {
   322 	for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
   323 	  (*this)[j->first]-=j->second;
   324 	const_comp-=e.const_comp;
   325 	return *this;
   326       }
   327       ///\e
   328       Expr &operator*=(const Value &c) {
   329 	for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
   330 	  j->second*=c;
   331 	const_comp*=c;
   332 	return *this;
   333       }
   334       ///\e
   335       Expr &operator/=(const Value &c) {
   336 	for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
   337 	  j->second/=c;
   338 	const_comp/=c;
   339 	return *this;
   340       }
   341     };
   342     
   343     ///Linear constraint
   344 
   345     ///This data stucture represents a linear constraint in the LP.
   346     ///Basically it is a linear expression with a lower or an upper bound
   347     ///(or both). These parts of the constraint can be obtained by the member
   348     ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
   349     ///respectively.
   350     ///There are two ways to construct a constraint.
   351     ///- You can set the linear expression and the bounds directly
   352     ///  by the functions above.
   353     ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
   354     ///  are defined between expressions, or even between constraints whenever
   355     ///  it makes sense. Therefore if \c e and \c f are linear expressions and
   356     ///  \c s and \c t are numbers, then the followings are valid expressions
   357     ///  and thus they can be used directly e.g. in \ref addRow() whenever
   358     ///  it makes sense.
   359     ///  \code
   360     ///  e<=s
   361     ///  e<=f
   362     ///  s<=e<=t
   363     ///  e>=t
   364     ///  \endcode
   365     ///\warning The validity of a constraint is checked only at run time, so
   366     ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw a
   367     ///\ref LogicError exception.
   368     class Constr
   369     {
   370     public:
   371       typedef LpSolverBase::Expr Expr;
   372       typedef Expr::Key Key;
   373       typedef Expr::Value Value;
   374       
   375 //       static const Value INF;
   376 //       static const Value NaN;
   377 
   378     protected:
   379       Expr _expr;
   380       Value _lb,_ub;
   381     public:
   382       ///\e
   383       Constr() : _expr(), _lb(NaN), _ub(NaN) {}
   384       ///\e
   385       Constr(Value lb,const Expr &e,Value ub) :
   386 	_expr(e), _lb(lb), _ub(ub) {}
   387       ///\e
   388       Constr(const Expr &e,Value ub) : 
   389 	_expr(e), _lb(NaN), _ub(ub) {}
   390       ///\e
   391       Constr(Value lb,const Expr &e) :
   392 	_expr(e), _lb(lb), _ub(NaN) {}
   393       ///\e
   394       Constr(const Expr &e) : 
   395 	_expr(e), _lb(NaN), _ub(NaN) {}
   396       ///\e
   397       void clear() 
   398       {
   399 	_expr.clear();
   400 	_lb=_ub=NaN;
   401       }
   402 
   403       ///Reference to the linear expression 
   404       Expr &expr() { return _expr; }
   405       ///Cont reference to the linear expression 
   406       const Expr &expr() const { return _expr; }
   407       ///Reference to the lower bound.
   408 
   409       ///\return
   410       ///- \ref INF "INF": the constraint is lower unbounded.
   411       ///- \ref NaN "NaN": lower bound has not been set.
   412       ///- finite number: the lower bound
   413       Value &lowerBound() { return _lb; }
   414       ///The const version of \ref lowerBound()
   415       const Value &lowerBound() const { return _lb; }
   416       ///Reference to the upper bound.
   417 
   418       ///\return
   419       ///- \ref INF "INF": the constraint is upper unbounded.
   420       ///- \ref NaN "NaN": upper bound has not been set.
   421       ///- finite number: the upper bound
   422       Value &upperBound() { return _ub; }
   423       ///The const version of \ref upperBound()
   424       const Value &upperBound() const { return _ub; }
   425       ///Is the constraint lower bounded?
   426       bool lowerBounded() const { 
   427 	using namespace std;
   428 	return finite(_lb);
   429       }
   430       ///Is the constraint upper bounded?
   431       bool upperBounded() const {
   432 	using namespace std;
   433 	return finite(_ub);
   434       }
   435     };
   436     
   437     ///Linear expression of rows
   438     
   439     ///This data structure represents a column of the matrix,
   440     ///thas is it strores a linear expression of the dual variables
   441     ///(\ref Row "Row"s).
   442     ///
   443     ///There are several ways to access and modify the contents of this
   444     ///container.
   445     ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
   446     ///if \c e is an DualExpr and \c v
   447     ///and \c w are of type \ref Row, then you can
   448     ///read and modify the coefficients like
   449     ///these.
   450     ///\code
   451     ///e[v]=5;
   452     ///e[v]+=12;
   453     ///e.erase(v);
   454     ///\endcode
   455     ///or you can also iterate through its elements.
   456     ///\code
   457     ///double s=0;
   458     ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
   459     ///  s+=i->second;
   460     ///\endcode
   461     ///(This code computes the sum of all coefficients).
   462     ///- Numbers (<tt>double</tt>'s)
   463     ///and variables (\ref Row "Row"s) directly convert to an
   464     ///\ref DualExpr and the usual linear operations are defined so  
   465     ///\code
   466     ///v+w
   467     ///2*v-3.12*(v-w/2)
   468     ///v*2.1+(3*v+(v*12+w)*3)/2
   469     ///\endcode
   470     ///are valid \ref DualExpr "DualExpr"essions.
   471     ///The usual assignment operations are also defined.
   472     ///\code
   473     ///e=v+w;
   474     ///e+=2*v-3.12*(v-w/2);
   475     ///e*=3.4;
   476     ///e/=5;
   477     ///\endcode
   478     ///
   479     ///\sa Expr
   480     ///
   481     class DualExpr : public std::map<Row,Value>
   482     {
   483     public:
   484       typedef LpSolverBase::Row Key; 
   485       typedef LpSolverBase::Value Value;
   486       
   487     protected:
   488       typedef std::map<Row,Value> Base;
   489       
   490     public:
   491       typedef True IsLinExpression;
   492       ///\e
   493       DualExpr() : Base() { }
   494       ///\e
   495       DualExpr(const Key &v) {
   496 	Base::insert(std::make_pair(v, 1));
   497       }
   498       ///\e
   499       void set(const Key &v,const Value &c) {
   500 	Base::insert(std::make_pair(v, c));
   501       }
   502       
   503       ///Removes the components with zero coefficient.
   504       void simplify() {
   505 	for (Base::iterator i=Base::begin(); i!=Base::end();) {
   506 	  Base::iterator j=i;
   507 	  ++j;
   508 	  if ((*i).second==0) Base::erase(i);
   509 	  j=i;
   510 	}
   511       }
   512 
   513       ///Removes the coefficients closer to zero than \c tolerance.
   514       void simplify(double &tolerance) {
   515 	for (Base::iterator i=Base::begin(); i!=Base::end();) {
   516 	  Base::iterator j=i;
   517 	  ++j;
   518 	  if (std::fabs((*i).second)<tolerance) Base::erase(i);
   519 	  j=i;
   520 	}
   521       }
   522 
   523 
   524       ///Sets all coefficients to 0.
   525       void clear() {
   526 	Base::clear();
   527       }
   528 
   529       ///\e
   530       DualExpr &operator+=(const DualExpr &e) {
   531 	for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
   532 	  (*this)[j->first]+=j->second;
   533 	return *this;
   534       }
   535       ///\e
   536       DualExpr &operator-=(const DualExpr &e) {
   537 	for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
   538 	  (*this)[j->first]-=j->second;
   539 	return *this;
   540       }
   541       ///\e
   542       DualExpr &operator*=(const Value &c) {
   543 	for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
   544 	  j->second*=c;
   545 	return *this;
   546       }
   547       ///\e
   548       DualExpr &operator/=(const Value &c) {
   549 	for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
   550 	  j->second/=c;
   551 	return *this;
   552       }
   553     };
   554     
   555 
   556   protected:
   557     _FixId rows;
   558     _FixId cols;
   559 
   560     //Abstract virtual functions
   561     virtual LpSolverBase &_newLp() = 0;
   562     virtual LpSolverBase &_copyLp(){
   563       ///\todo This should be implemented here, too,  when we have problem retrieving routines. It can be overriden.
   564 
   565       //Starting:
   566       LpSolverBase & newlp(_newLp());
   567       return newlp;
   568       //return *(LpSolverBase*)0;
   569     };
   570 
   571     virtual int _addCol() = 0;
   572     virtual int _addRow() = 0;
   573     virtual void _eraseCol(int col) = 0;
   574     virtual void _eraseRow(int row) = 0;
   575     virtual void _setRowCoeffs(int i, 
   576 			       int length,
   577                                int  const * indices, 
   578                                Value  const * values ) = 0;
   579     virtual void _setColCoeffs(int i, 
   580 			       int length,
   581                                int  const * indices, 
   582                                Value  const * values ) = 0;
   583     virtual void _setCoeff(int row, int col, Value value) = 0;
   584     virtual void _setColLowerBound(int i, Value value) = 0;
   585     virtual void _setColUpperBound(int i, Value value) = 0;
   586 //     virtual void _setRowLowerBound(int i, Value value) = 0;
   587 //     virtual void _setRowUpperBound(int i, Value value) = 0;
   588     virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
   589     virtual void _setObjCoeff(int i, Value obj_coef) = 0;
   590     virtual void _clearObj()=0;
   591 //     virtual void _setObj(int length,
   592 //                          int  const * indices, 
   593 //                          Value  const * values ) = 0;
   594     virtual SolveExitStatus _solve() = 0;
   595     virtual Value _getPrimal(int i) = 0;
   596     virtual Value _getDual(int i) = 0;
   597     virtual Value _getPrimalValue() = 0;
   598     virtual SolutionStatus _getPrimalStatus() = 0;
   599     virtual SolutionStatus _getDualStatus() = 0;
   600     ///\todo This could be implemented here, too, using _getPrimalStatus() and
   601     ///_getDualStatus()
   602     virtual ProblemTypes _getProblemType() = 0;
   603 
   604     virtual void _setMax() = 0;
   605     virtual void _setMin() = 0;
   606     
   607     //Own protected stuff
   608     
   609     //Constant component of the objective function
   610     Value obj_const_comp;
   611     
   612 
   613 
   614     
   615   public:
   616 
   617     ///\e
   618     LpSolverBase() : obj_const_comp(0) {}
   619 
   620     ///\e
   621     virtual ~LpSolverBase() {}
   622 
   623     ///Creates a new LP problem
   624     LpSolverBase &newLp() {return _newLp();}
   625     ///Makes a copy of the LP problem
   626     LpSolverBase &copyLp() {return _copyLp();}
   627     
   628     ///\name Build up and modify the LP
   629 
   630     ///@{
   631 
   632     ///Add a new empty column (i.e a new variable) to the LP
   633     Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
   634 
   635     ///\brief Adds several new columns
   636     ///(i.e a variables) at once
   637     ///
   638     ///This magic function takes a container as its argument
   639     ///and fills its elements
   640     ///with new columns (i.e. variables)
   641     ///\param t can be
   642     ///- a standard STL compatible iterable container with
   643     ///\ref Col as its \c values_type
   644     ///like
   645     ///\code
   646     ///std::vector<LpSolverBase::Col>
   647     ///std::list<LpSolverBase::Col>
   648     ///\endcode
   649     ///- a standard STL compatible iterable container with
   650     ///\ref Col as its \c mapped_type
   651     ///like
   652     ///\code
   653     ///std::map<AnyType,LpSolverBase::Col>
   654     ///\endcode
   655     ///- an iterable lemon \ref concept::WriteMap "write map" like 
   656     ///\code
   657     ///ListGraph::NodeMap<LpSolverBase::Col>
   658     ///ListGraph::EdgeMap<LpSolverBase::Col>
   659     ///\endcode
   660     ///\return The number of the created column.
   661 #ifdef DOXYGEN
   662     template<class T>
   663     int addColSet(T &t) { return 0;} 
   664 #else
   665     template<class T>
   666     typename enable_if<typename T::value_type::LpSolverCol,int>::type
   667     addColSet(T &t,dummy<0> = 0) {
   668       int s=0;
   669       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
   670       return s;
   671     }
   672     template<class T>
   673     typename enable_if<typename T::value_type::second_type::LpSolverCol,
   674 		       int>::type
   675     addColSet(T &t,dummy<1> = 1) { 
   676       int s=0;
   677       for(typename T::iterator i=t.begin();i!=t.end();++i) {
   678 	i->second=addCol();
   679 	s++;
   680       }
   681       return s;
   682     }
   683     template<class T>
   684     typename enable_if<typename T::MapIt::Value::LpSolverCol,
   685 		       int>::type
   686     addColSet(T &t,dummy<2> = 2) { 
   687       int s=0;
   688       for(typename T::MapIt i(t); i!=INVALID; ++i)
   689 	{
   690 	  i.set(addCol());
   691 	  s++;
   692 	}
   693       return s;
   694     }
   695 #endif
   696 
   697     ///Set a column (i.e a dual constraint) of the LP
   698 
   699     ///\param c is the column to be modified
   700     ///\param e is a dual linear expression (see \ref DualExpr)
   701     ///\bug This is a temporary function. The interface will change to
   702     ///a better one.
   703     void setCol(Col c,const DualExpr &e) {
   704       std::vector<int> indices;
   705       std::vector<Value> values;
   706       indices.push_back(0);
   707       values.push_back(0);
   708       for(DualExpr::const_iterator i=e.begin(); i!=e.end(); ++i)
   709 	if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
   710 	  indices.push_back(rows.floatingId((*i).first.id));
   711 	  values.push_back((*i).second);
   712 	}
   713       _setColCoeffs(cols.floatingId(c.id),indices.size()-1,
   714 		    &indices[0],&values[0]);
   715     }
   716 
   717     ///Add a new column to the LP
   718 
   719     ///\param e is a dual linear expression (see \ref DualExpr)
   720     ///\param obj is the corresponding component of the objective
   721     ///function. It is 0 by default.
   722     ///\return The created column.
   723     ///\bug This is a temportary function. The interface will change to
   724     ///a better one.
   725     Col addCol(const DualExpr &e, Value obj=0) {
   726       Col c=addCol();
   727       setCol(c,e);
   728       objCoeff(c,obj);
   729       return c;
   730     }
   731 
   732     ///Add a new empty row (i.e a new constraint) to the LP
   733 
   734     ///This function adds a new empty row (i.e a new constraint) to the LP.
   735     ///\return The created row
   736     Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
   737 
   738     ///\brief Add several new rows
   739     ///(i.e a constraints) at once
   740     ///
   741     ///This magic function takes a container as its argument
   742     ///and fills its elements
   743     ///with new row (i.e. variables)
   744     ///\param t can be
   745     ///- a standard STL compatible iterable container with
   746     ///\ref Row as its \c values_type
   747     ///like
   748     ///\code
   749     ///std::vector<LpSolverBase::Row>
   750     ///std::list<LpSolverBase::Row>
   751     ///\endcode
   752     ///- a standard STL compatible iterable container with
   753     ///\ref Row as its \c mapped_type
   754     ///like
   755     ///\code
   756     ///std::map<AnyType,LpSolverBase::Row>
   757     ///\endcode
   758     ///- an iterable lemon \ref concept::WriteMap "write map" like 
   759     ///\code
   760     ///ListGraph::NodeMap<LpSolverBase::Row>
   761     ///ListGraph::EdgeMap<LpSolverBase::Row>
   762     ///\endcode
   763     ///\return The number of rows created.
   764 #ifdef DOXYGEN
   765     template<class T>
   766     int addRowSet(T &t) { return 0;} 
   767 #else
   768     template<class T>
   769     typename enable_if<typename T::value_type::LpSolverRow,int>::type
   770     addRowSet(T &t,dummy<0> = 0) {
   771       int s=0;
   772       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
   773       return s;
   774     }
   775     template<class T>
   776     typename enable_if<typename T::value_type::second_type::LpSolverRow,
   777 		       int>::type
   778     addRowSet(T &t,dummy<1> = 1) { 
   779       int s=0;
   780       for(typename T::iterator i=t.begin();i!=t.end();++i) {
   781 	i->second=addRow();
   782 	s++;
   783       }
   784       return s;
   785     }
   786     template<class T>
   787     typename enable_if<typename T::MapIt::Value::LpSolverRow,
   788 		       int>::type
   789     addRowSet(T &t,dummy<2> = 2) { 
   790       int s=0;
   791       for(typename T::MapIt i(t); i!=INVALID; ++i)
   792 	{
   793 	  i.set(addRow());
   794 	  s++;
   795 	}
   796       return s;
   797     }
   798 #endif
   799 
   800     ///Set a row (i.e a constraint) of the LP
   801 
   802     ///\param r is the row to be modified
   803     ///\param l is lower bound (-\ref INF means no bound)
   804     ///\param e is a linear expression (see \ref Expr)
   805     ///\param u is the upper bound (\ref INF means no bound)
   806     ///\bug This is a temportary function. The interface will change to
   807     ///a better one.
   808     ///\todo Option to control whether a constraint with a single variable is
   809     ///added or not.
   810     void setRow(Row r, Value l,const Expr &e, Value u) {
   811       std::vector<int> indices;
   812       std::vector<Value> values;
   813       indices.push_back(0);
   814       values.push_back(0);
   815       for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
   816 	if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
   817 	  indices.push_back(cols.floatingId((*i).first.id));
   818 	  values.push_back((*i).second);
   819 	}
   820       _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
   821 		    &indices[0],&values[0]);
   822 //       _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
   823 //       _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
   824        _setRowBounds(rows.floatingId(r.id),l-e.constComp(),u-e.constComp());
   825     }
   826 
   827     ///Set a row (i.e a constraint) of the LP
   828 
   829     ///\param r is the row to be modified
   830     ///\param c is a linear expression (see \ref Constr)
   831     void setRow(Row r, const Constr &c) {
   832       setRow(r,
   833 	     c.lowerBounded()?c.lowerBound():-INF,
   834 	     c.expr(),
   835 	     c.upperBounded()?c.upperBound():INF);
   836     }
   837 
   838     ///Add a new row (i.e a new constraint) to the LP
   839 
   840     ///\param l is the lower bound (-\ref INF means no bound)
   841     ///\param e is a linear expression (see \ref Expr)
   842     ///\param u is the upper bound (\ref INF means no bound)
   843     ///\return The created row.
   844     ///\bug This is a temportary function. The interface will change to
   845     ///a better one.
   846     Row addRow(Value l,const Expr &e, Value u) {
   847       Row r=addRow();
   848       setRow(r,l,e,u);
   849       return r;
   850     }
   851 
   852     ///Add a new row (i.e a new constraint) to the LP
   853 
   854     ///\param c is a linear expression (see \ref Constr)
   855     ///\return The created row.
   856     Row addRow(const Constr &c) {
   857       Row r=addRow();
   858       setRow(r,c);
   859       return r;
   860     }
   861     ///Erase a coloumn (i.e a variable) from the LP
   862 
   863     ///\param c is the coloumn to be deleted
   864     ///\todo Please check this
   865     void eraseCol(Col c) {
   866       _eraseCol(cols.floatingId(c.id));
   867       cols.erase(c.id);
   868     }
   869     ///Erase a  row (i.e a constraint) from the LP
   870 
   871     ///\param r is the row to be deleted
   872     ///\todo Please check this
   873     void eraseRow(Row r) {
   874       _eraseRow(rows.floatingId(r.id));
   875       rows.erase(r.id);
   876     }
   877 
   878     ///Set an element of the coefficient matrix of the LP
   879 
   880     ///\param r is the row of the element to be modified
   881     ///\param c is the coloumn of the element to be modified
   882     ///\param val is the new value of the coefficient
   883     void setCoeff(Row r, Col c, Value val){
   884       _setCoeff(rows.floatingId(r.id),cols.floatingId(c.id), val);
   885     }
   886 
   887     /// Set the lower bound of a column (i.e a variable)
   888 
   889     /// The upper bound of a variable (column) has to be given by an 
   890     /// extended number of type Value, i.e. a finite number of type 
   891     /// Value or -\ref INF.
   892     void colLowerBound(Col c, Value value) {
   893       _setColLowerBound(cols.floatingId(c.id),value);
   894     }
   895     /// Set the upper bound of a column (i.e a variable)
   896 
   897     /// The upper bound of a variable (column) has to be given by an 
   898     /// extended number of type Value, i.e. a finite number of type 
   899     /// Value or \ref INF.
   900     void colUpperBound(Col c, Value value) {
   901       _setColUpperBound(cols.floatingId(c.id),value);
   902     };
   903     /// Set the lower and the upper bounds of a column (i.e a variable)
   904 
   905     /// The lower and the upper bounds of
   906     /// a variable (column) have to be given by an 
   907     /// extended number of type Value, i.e. a finite number of type 
   908     /// Value, -\ref INF or \ref INF.
   909     void colBounds(Col c, Value lower, Value upper) {
   910       _setColLowerBound(cols.floatingId(c.id),lower);
   911       _setColUpperBound(cols.floatingId(c.id),upper);
   912     }
   913     
   914 //     /// Set the lower bound of a row (i.e a constraint)
   915 
   916 //     /// The lower bound of a linear expression (row) has to be given by an 
   917 //     /// extended number of type Value, i.e. a finite number of type 
   918 //     /// Value or -\ref INF.
   919 //     void rowLowerBound(Row r, Value value) {
   920 //       _setRowLowerBound(rows.floatingId(r.id),value);
   921 //     };
   922 //     /// Set the upper bound of a row (i.e a constraint)
   923 
   924 //     /// The upper bound of a linear expression (row) has to be given by an 
   925 //     /// extended number of type Value, i.e. a finite number of type 
   926 //     /// Value or \ref INF.
   927 //     void rowUpperBound(Row r, Value value) {
   928 //       _setRowUpperBound(rows.floatingId(r.id),value);
   929 //     };
   930 
   931     /// Set the lower and the upper bounds of a row (i.e a constraint)
   932 
   933     /// The lower and the upper bounds of
   934     /// a constraint (row) have to be given by an 
   935     /// extended number of type Value, i.e. a finite number of type 
   936     /// Value, -\ref INF or \ref INF.
   937     void rowBounds(Row c, Value lower, Value upper) {
   938       _setRowBounds(rows.floatingId(c.id),lower, upper);
   939       // _setRowUpperBound(rows.floatingId(c.id),upper);
   940     }
   941     
   942     ///Set an element of the objective function
   943     void objCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
   944     ///Set the objective function
   945     
   946     ///\param e is a linear expression of type \ref Expr.
   947     ///\bug The previous objective function is not cleared!
   948     void setObj(Expr e) {
   949       _clearObj();
   950       for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
   951 	objCoeff((*i).first,(*i).second);
   952       obj_const_comp=e.constComp();
   953     }
   954 
   955     ///Maximize
   956     void max() { _setMax(); }
   957     ///Minimize
   958     void min() { _setMin(); }
   959 
   960     
   961     ///@}
   962 
   963 
   964     ///\name Solve the LP
   965 
   966     ///@{
   967 
   968     ///\e Solve the LP problem at hand
   969     ///
   970     ///\return The result of the optimization procedure. Possible values and their meanings can be found in the documentation of \ref SolveExitStatus.
   971     ///
   972     ///\todo Which method is used to solve the problem
   973     SolveExitStatus solve() { return _solve(); }
   974     
   975     ///@}
   976     
   977     ///\name Obtain the solution
   978 
   979     ///@{
   980 
   981     /// The status of the primal problem (the original LP problem)
   982     SolutionStatus primalStatus() {
   983       return _getPrimalStatus();
   984     }
   985 
   986     /// The status of the dual (of the original LP) problem 
   987     SolutionStatus dualStatus() {
   988       return _getDualStatus();
   989     }
   990 
   991     ///The type of the original LP problem
   992     ProblemTypes problemType() {
   993       return _getProblemType();
   994     }
   995 
   996     ///\e
   997     Value primal(Col c) { return _getPrimal(cols.floatingId(c.id)); }
   998 
   999     ///\e
  1000     Value dual(Row r) { return _getDual(rows.floatingId(r.id)); }
  1001 
  1002     ///\e
  1003 
  1004     ///\return
  1005     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
  1006     /// of the primal problem, depending on whether we minimize or maximize.
  1007     ///- \ref NaN if no primal solution is found.
  1008     ///- The (finite) objective value if an optimal solution is found.
  1009     Value primalValue() { return _getPrimalValue()+obj_const_comp;}
  1010     ///@}
  1011     
  1012   };  
  1013 
  1014   ///\e
  1015   
  1016   ///\relates LpSolverBase::Expr
  1017   ///
  1018   inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
  1019 				      const LpSolverBase::Expr &b) 
  1020   {
  1021     LpSolverBase::Expr tmp(a);
  1022     tmp+=b;
  1023     return tmp;
  1024   }
  1025   ///\e
  1026   
  1027   ///\relates LpSolverBase::Expr
  1028   ///
  1029   inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
  1030 				      const LpSolverBase::Expr &b) 
  1031   {
  1032     LpSolverBase::Expr tmp(a);
  1033     tmp-=b;
  1034     return tmp;
  1035   }
  1036   ///\e
  1037   
  1038   ///\relates LpSolverBase::Expr
  1039   ///
  1040   inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
  1041 				      const LpSolverBase::Value &b) 
  1042   {
  1043     LpSolverBase::Expr tmp(a);
  1044     tmp*=b;
  1045     return tmp;
  1046   }
  1047   
  1048   ///\e
  1049   
  1050   ///\relates LpSolverBase::Expr
  1051   ///
  1052   inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
  1053 				      const LpSolverBase::Expr &b) 
  1054   {
  1055     LpSolverBase::Expr tmp(b);
  1056     tmp*=a;
  1057     return tmp;
  1058   }
  1059   ///\e
  1060   
  1061   ///\relates LpSolverBase::Expr
  1062   ///
  1063   inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
  1064 				      const LpSolverBase::Value &b) 
  1065   {
  1066     LpSolverBase::Expr tmp(a);
  1067     tmp/=b;
  1068     return tmp;
  1069   }
  1070   
  1071   ///\e
  1072   
  1073   ///\relates LpSolverBase::Constr
  1074   ///
  1075   inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
  1076 					 const LpSolverBase::Expr &f) 
  1077   {
  1078     return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
  1079   }
  1080 
  1081   ///\e
  1082   
  1083   ///\relates LpSolverBase::Constr
  1084   ///
  1085   inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
  1086 					 const LpSolverBase::Expr &f) 
  1087   {
  1088     return LpSolverBase::Constr(e,f);
  1089   }
  1090 
  1091   ///\e
  1092   
  1093   ///\relates LpSolverBase::Constr
  1094   ///
  1095   inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
  1096 					 const LpSolverBase::Value &f) 
  1097   {
  1098     return LpSolverBase::Constr(e,f);
  1099   }
  1100 
  1101   ///\e
  1102   
  1103   ///\relates LpSolverBase::Constr
  1104   ///
  1105   inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
  1106 					 const LpSolverBase::Expr &f) 
  1107   {
  1108     return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
  1109   }
  1110 
  1111 
  1112   ///\e
  1113   
  1114   ///\relates LpSolverBase::Constr
  1115   ///
  1116   inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
  1117 					 const LpSolverBase::Expr &f) 
  1118   {
  1119     return LpSolverBase::Constr(f,e);
  1120   }
  1121 
  1122 
  1123   ///\e
  1124   
  1125   ///\relates LpSolverBase::Constr
  1126   ///
  1127   inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
  1128 					 const LpSolverBase::Value &f) 
  1129   {
  1130     return LpSolverBase::Constr(f,e);
  1131   }
  1132 
  1133   ///\e
  1134   
  1135   ///\relates LpSolverBase::Constr
  1136   ///
  1137   inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
  1138 					 const LpSolverBase::Expr &f) 
  1139   {
  1140     return LpSolverBase::Constr(0,e-f,0);
  1141   }
  1142 
  1143   ///\e
  1144   
  1145   ///\relates LpSolverBase::Constr
  1146   ///
  1147   inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
  1148 					 const LpSolverBase::Constr&c) 
  1149   {
  1150     LpSolverBase::Constr tmp(c);
  1151     ///\todo Create an own exception type.
  1152     if(!isnan(tmp.lowerBound())) throw LogicError();
  1153     else tmp.lowerBound()=n;
  1154     return tmp;
  1155   }
  1156   ///\e
  1157   
  1158   ///\relates LpSolverBase::Constr
  1159   ///
  1160   inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
  1161 					 const LpSolverBase::Value &n)
  1162   {
  1163     LpSolverBase::Constr tmp(c);
  1164     ///\todo Create an own exception type.
  1165     if(!isnan(tmp.upperBound())) throw LogicError();
  1166     else tmp.upperBound()=n;
  1167     return tmp;
  1168   }
  1169 
  1170   ///\e
  1171   
  1172   ///\relates LpSolverBase::Constr
  1173   ///
  1174   inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
  1175 					 const LpSolverBase::Constr&c) 
  1176   {
  1177     LpSolverBase::Constr tmp(c);
  1178     ///\todo Create an own exception type.
  1179     if(!isnan(tmp.upperBound())) throw LogicError();
  1180     else tmp.upperBound()=n;
  1181     return tmp;
  1182   }
  1183   ///\e
  1184   
  1185   ///\relates LpSolverBase::Constr
  1186   ///
  1187   inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
  1188 					 const LpSolverBase::Value &n)
  1189   {
  1190     LpSolverBase::Constr tmp(c);
  1191     ///\todo Create an own exception type.
  1192     if(!isnan(tmp.lowerBound())) throw LogicError();
  1193     else tmp.lowerBound()=n;
  1194     return tmp;
  1195   }
  1196 
  1197   ///\e
  1198   
  1199   ///\relates LpSolverBase::DualExpr
  1200   ///
  1201   inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
  1202 				      const LpSolverBase::DualExpr &b) 
  1203   {
  1204     LpSolverBase::DualExpr tmp(a);
  1205     tmp+=b;
  1206     return tmp;
  1207   }
  1208   ///\e
  1209   
  1210   ///\relates LpSolverBase::DualExpr
  1211   ///
  1212   inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
  1213 				      const LpSolverBase::DualExpr &b) 
  1214   {
  1215     LpSolverBase::DualExpr tmp(a);
  1216     tmp-=b;
  1217     return tmp;
  1218   }
  1219   ///\e
  1220   
  1221   ///\relates LpSolverBase::DualExpr
  1222   ///
  1223   inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
  1224 				      const LpSolverBase::Value &b) 
  1225   {
  1226     LpSolverBase::DualExpr tmp(a);
  1227     tmp*=b;
  1228     return tmp;
  1229   }
  1230   
  1231   ///\e
  1232   
  1233   ///\relates LpSolverBase::DualExpr
  1234   ///
  1235   inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
  1236 				      const LpSolverBase::DualExpr &b) 
  1237   {
  1238     LpSolverBase::DualExpr tmp(b);
  1239     tmp*=a;
  1240     return tmp;
  1241   }
  1242   ///\e
  1243   
  1244   ///\relates LpSolverBase::DualExpr
  1245   ///
  1246   inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
  1247 				      const LpSolverBase::Value &b) 
  1248   {
  1249     LpSolverBase::DualExpr tmp(a);
  1250     tmp/=b;
  1251     return tmp;
  1252   }
  1253   
  1254 
  1255 } //namespace lemon
  1256 
  1257 #endif //LEMON_LP_BASE_H