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