src/lemon/lp_base.h
author athos
Thu, 07 Apr 2005 15:22:03 +0000
changeset 1319 6e277ba3fc76
parent 1305 c3dc75d4af24
child 1323 3aaadfb7de3d
permissions -rw-r--r--
Cplex interface has improved a lot.
     1 /* -*- C++ -*-
     2  * src/lemon/lp_base.h - Part of LEMON, a generic C++ optimization library
     3  *
     4  * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     5  * (Egervary Combinatorial Optimization Research Group, 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<math.h>
    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 namespace lemon {
    34   
    35   ///Internal data structure to convert floating id's to fix one's
    36     
    37   ///\todo This might be implemented to be also usable in other places.
    38   class _FixId 
    39   {
    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) {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) { 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   class LpSolverBase {
   104     
   105   public:
   106 
   107     ///\e
   108     enum SolveExitStatus {
   109       ///\e
   110       SOLVED = 0,
   111       ///\e
   112       UNSOLVED = 1
   113     };
   114       
   115     ///\e
   116     enum SolutionStatus {
   117       ///Feasible solution has'n been found (but may exist).
   118 
   119       ///\todo NOTFOUND might be a better name.
   120       ///
   121       UNDEFINED = 0,
   122       ///The problem has no feasible solution
   123       INFEASIBLE = 1,
   124       ///Feasible solution found
   125       FEASIBLE = 2,
   126       ///Optimal solution exists and found
   127       OPTIMAL = 3,
   128       ///The cost function is unbounded
   129 
   130       ///\todo Give a feasible solution and an infinite ray (and the
   131       ///corresponding bases)
   132       INFINITE = 4
   133     };
   134       
   135     ///The floating point type used by the solver
   136     typedef double Value;
   137     ///The infinity constant
   138     static const Value INF;
   139     ///The not a number constant
   140     static const Value NaN;
   141     
   142     ///Refer to a column of the LP.
   143 
   144     ///This type is used to refer to a column of the LP.
   145     ///
   146     ///Its value remains valid and correct even after the addition or erase of
   147     ///other columns.
   148     ///
   149     ///\todo Document what can one do with a Col (INVALID, comparing,
   150     ///it is similar to Node/Edge)
   151     class Col {
   152     protected:
   153       int id;
   154       friend class LpSolverBase;
   155     public:
   156       typedef Value ExprValue;
   157       typedef True LpSolverCol;
   158       Col() {}
   159       Col(const Invalid&) : id(-1) {}
   160       bool operator<(Col c) const  {return id<c.id;}
   161       bool operator==(Col c) const  {return id==c.id;}
   162       bool operator!=(Col c) const  {return id==c.id;}
   163     };
   164 
   165     ///Refer to a row of the LP.
   166 
   167     ///This type is used to refer to a row of the LP.
   168     ///
   169     ///Its value remains valid and correct even after the addition or erase of
   170     ///other rows.
   171     ///
   172     ///\todo Document what can one do with a Row (INVALID, comparing,
   173     ///it is similar to Node/Edge)
   174     class Row {
   175     protected:
   176       int id;
   177       friend class LpSolverBase;
   178     public:
   179       typedef Value ExprValue;
   180       typedef True LpSolverRow;
   181       Row() {}
   182       Row(const Invalid&) : id(-1) {}
   183       typedef True LpSolverRow;
   184       bool operator<(Row c) const  {return id<c.id;}
   185       bool operator==(Row c) const  {return id==c.id;}
   186       bool operator!=(Row c) const  {return id==c.id;} 
   187    };
   188     
   189     ///Linear expression of variables and a constant component
   190     
   191     ///This data structure strores a linear expression of the variables
   192     ///(\ref Col "Col"s) and also has a constant component.
   193     ///
   194     ///There are several ways to access and modify the contents of this
   195     ///container.
   196     ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
   197     ///if \c e is an Expr and \c v and \c w are of type \ref Col then you can
   198     ///read and modify the coefficients like
   199     ///these.
   200     ///\code
   201     ///e[v]=5;
   202     ///e[v]+=12;
   203     ///e.erase(v);
   204     ///\endcode
   205     ///or you can also iterate through its elements.
   206     ///\code
   207     ///double s=0;
   208     ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
   209     ///  s+=i->second;
   210     ///\endcode
   211     ///(This code computes the sum of all coefficients).
   212     ///- Numbers (<tt>double</tt>'s)
   213     ///and variables (\ref Col "Col"s) directly convert to an
   214     ///\ref Expr and the usual linear operations are defined so  
   215     ///\code
   216     ///v+w
   217     ///2*v-3.12*(v-w/2)+2
   218     ///v*2.1+(3*v+(v*12+w+6)*3)/2
   219     ///\endcode
   220     ///are valid expressions. The usual assignment operations are also defined.
   221     ///\code
   222     ///e=v+w;
   223     ///e+=2*v-3.12*(v-w/2)+2;
   224     ///e*=3.4;
   225     ///e/=5;
   226     ///\endcode
   227     ///- The constant member can be set and read by \ref constComp()
   228     ///\code
   229     ///e.constComp()=12;
   230     ///double c=e.constComp();
   231     ///\endcode
   232     ///
   233     ///\note that \ref clear() not only sets all coefficients to 0 but also
   234     ///clears the constant components.
   235     class Expr : public std::map<Col,Value>
   236     {
   237     public:
   238       typedef LpSolverBase::Col Key; 
   239       typedef LpSolverBase::Value Value;
   240       
   241     protected:
   242       typedef std::map<Col,Value> Base;
   243       
   244       Value const_comp;
   245   public:
   246       typedef True IsLinExpression;
   247       ///\e
   248       Expr() : Base(), const_comp(0) { }
   249       ///\e
   250       Expr(const Key &v) : const_comp(0) {
   251 	Base::insert(std::make_pair(v, 1));
   252       }
   253       ///\e
   254       Expr(const Value &v) : const_comp(v) {}
   255       ///\e
   256       void set(const Key &v,const Value &c) {
   257 	Base::insert(std::make_pair(v, c));
   258       }
   259       ///\e
   260       Value &constComp() { return const_comp; }
   261       ///\e
   262       const Value &constComp() const { return const_comp; }
   263       
   264       ///Removes the components with zero coefficient.
   265       void simplify() {
   266 	for (Base::iterator i=Base::begin(); i!=Base::end();) {
   267 	  Base::iterator j=i;
   268 	  ++j;
   269 	  if ((*i).second==0) Base::erase(i);
   270 	  j=i;
   271 	}
   272       }
   273 
   274       ///Sets all coefficients and the constant component to 0.
   275       void clear() {
   276 	Base::clear();
   277 	const_comp=0;
   278       }
   279 
   280       ///\e
   281       Expr &operator+=(const Expr &e) {
   282 	for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
   283 	  (*this)[j->first]+=j->second;
   284 	///\todo it might be speeded up using "hints"
   285 	const_comp+=e.const_comp;
   286 	return *this;
   287       }
   288       ///\e
   289       Expr &operator-=(const Expr &e) {
   290 	for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
   291 	  (*this)[j->first]-=j->second;
   292 	const_comp-=e.const_comp;
   293 	return *this;
   294       }
   295       ///\e
   296       Expr &operator*=(const Value &c) {
   297 	for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
   298 	  j->second*=c;
   299 	const_comp*=c;
   300 	return *this;
   301       }
   302       ///\e
   303       Expr &operator/=(const Value &c) {
   304 	for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
   305 	  j->second/=c;
   306 	const_comp/=c;
   307 	return *this;
   308       }
   309     };
   310     
   311     ///Linear constraint
   312     //typedef LinConstr<Expr> Constr;
   313     class Constr
   314     {
   315     public:
   316       typedef LpSolverBase::Expr Expr;
   317       typedef Expr::Key Key;
   318       typedef Expr::Value Value;
   319       
   320       static const Value INF;
   321       static const Value NaN;
   322       //     static const Value INF=0;
   323       //     static const Value NaN=1;
   324       
   325     protected:
   326       Expr _expr;
   327       Value _lb,_ub;
   328     public:
   329       ///\e
   330       Constr() : _expr(), _lb(NaN), _ub(NaN) {}
   331       ///\e
   332       Constr(Value lb,const Expr &e,Value ub) :
   333 	_expr(e), _lb(lb), _ub(ub) {}
   334       ///\e
   335       Constr(const Expr &e,Value ub) : 
   336 	_expr(e), _lb(NaN), _ub(ub) {}
   337       ///\e
   338       Constr(Value lb,const Expr &e) :
   339 	_expr(e), _lb(lb), _ub(NaN) {}
   340       ///\e
   341       Constr(const Expr &e) : 
   342 	_expr(e), _lb(NaN), _ub(NaN) {}
   343       ///\e
   344       void clear() 
   345       {
   346 	_expr.clear();
   347 	_lb=_ub=NaN;
   348       }
   349       ///\e
   350       Expr &expr() { return _expr; }
   351       ///\e
   352       const Expr &expr() const { return _expr; }
   353       ///\e
   354       Value &lowerBound() { return _lb; }
   355       ///\e
   356       const Value &lowerBound() const { return _lb; }
   357       ///\e
   358       Value &upperBound() { return _ub; }
   359       ///\e
   360       const Value &upperBound() const { return _ub; }
   361       ///\e
   362       bool lowerBounded() const { 
   363 	using namespace std;
   364 	return isfinite(_lb);
   365       }
   366       ///\e
   367       bool upperBounded() const {
   368 	using namespace std;
   369 	return isfinite(_ub);
   370       }
   371     };
   372     
   373 
   374   protected:
   375     _FixId rows;
   376     _FixId cols;
   377 
   378     virtual int _addCol() = 0;
   379     virtual int _addRow() = 0;
   380     virtual void _setRowCoeffs(int i, 
   381 			       int length,
   382                                int  const * indices, 
   383                                Value  const * values ) = 0;
   384     virtual void _setColCoeffs(int i, 
   385 			       int length,
   386                                int  const * indices, 
   387                                Value  const * values ) = 0;
   388     virtual void _setColLowerBound(int i, Value value) = 0;
   389     virtual void _setColUpperBound(int i, Value value) = 0;
   390     virtual void _setRowLowerBound(int i, Value value) = 0;
   391     virtual void _setRowUpperBound(int i, Value value) = 0;
   392     virtual void _setObjCoeff(int i, Value obj_coef) = 0;
   393     virtual SolveExitStatus _solve() = 0;
   394     virtual Value _getPrimal(int i) = 0;
   395     virtual Value _getPrimalValue() = 0;
   396     virtual SolutionStatus _getPrimalStatus() = 0;
   397     virtual void _setMax() = 0;
   398     virtual void _setMin() = 0;
   399     
   400 
   401     void clearObj() {}
   402   public:
   403 
   404 
   405     ///\e
   406     virtual ~LpSolverBase() {}
   407 
   408     ///\name Build up and modify of the LP
   409 
   410     ///@{
   411 
   412     ///Add a new empty column (i.e a new variable) to the LP
   413     Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
   414 
   415     ///\brief Adds several new columns
   416     ///(i.e a variables) at once
   417     ///
   418     ///This magic function takes a container as its argument
   419     ///and fills its elements
   420     ///with new columns (i.e. variables)
   421     ///\param t can be
   422     ///- a standard STL compatible iterable container with
   423     ///\ref Col as its \c values_type
   424     ///like
   425     ///\code
   426     ///std::vector<LpSolverBase::Col>
   427     ///std::list<LpSolverBase::Col>
   428     ///\endcode
   429     ///- a standard STL compatible iterable container with
   430     ///\ref Col as its \c mapped_type
   431     ///like
   432     ///\code
   433     ///std::map<AnyStatus,LpSolverBase::Col>
   434     ///\endcode
   435     ///- an iterable lemon \ref concept::WriteMap "write map" like 
   436     ///\code
   437     ///ListGraph::NodeMap<LpSolverBase::Col>
   438     ///ListGraph::EdgeMap<LpSolverBase::Col>
   439     ///\endcode
   440     ///\return The number of the created column.
   441     ///\bug Iterable nodemap hasn't been implemented yet.
   442 #ifdef DOXYGEN
   443     template<class T>
   444     int addColSet(T &t) { return 0;} 
   445 #else
   446     template<class T>
   447     typename enable_if<typename T::value_type::LpSolverCol,int>::type
   448     addColSet(T &t,dummy<0> = 0) {
   449       int s=0;
   450       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
   451       return s;
   452     }
   453     template<class T>
   454     typename enable_if<typename T::value_type::second_type::LpSolverCol,
   455 		       int>::type
   456     addColSet(T &t,dummy<1> = 1) { 
   457       int s=0;
   458       for(typename T::iterator i=t.begin();i!=t.end();++i) {
   459 	i->second=addCol();
   460 	s++;
   461       }
   462       return s;
   463     }
   464     template<class T>
   465     typename enable_if<typename T::ValueSet::value_type::LpSolverCol,
   466 		       int>::type
   467     addColSet(T &t,dummy<2> = 2) { 
   468       ///\bug <tt>return addColSet(t.valueSet());</tt> should also work.
   469       int s=0;
   470       for(typename T::ValueSet::iterator i=t.valueSet().begin();
   471 	  i!=t.valueSet().end();
   472 	  ++i)
   473 	{
   474 	  *i=addCol();
   475 	  s++;
   476 	}
   477       return s;
   478     }
   479 #endif
   480 
   481     ///Add a new empty row (i.e a new constaint) to the LP
   482 
   483     ///This function adds a new empty row (i.e a new constaint) to the LP.
   484     ///\return The created row
   485     Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
   486 
   487     ///Set a row (i.e a constaint) of the LP
   488 
   489     ///\param r is the row to be modified
   490     ///\param l is lower bound (-\ref INF means no bound)
   491     ///\param e is a linear expression (see \ref Expr)
   492     ///\param u is the upper bound (\ref INF means no bound)
   493     ///\bug This is a temportary function. The interface will change to
   494     ///a better one.
   495     void setRow(Row r, Value l,const Expr &e, Value u) {
   496       std::vector<int> indices;
   497       std::vector<Value> values;
   498       indices.push_back(0);
   499       values.push_back(0);
   500       for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
   501 	if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
   502 	  indices.push_back(cols.floatingId((*i).first.id));
   503 	  values.push_back((*i).second);
   504 	}
   505       _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
   506 		    &indices[0],&values[0]);
   507       _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
   508       _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
   509     }
   510 
   511     ///Set a row (i.e a constaint) of the LP
   512 
   513     ///\param r is the row to be modified
   514     ///\param c is a linear expression (see \ref Constr)
   515     void setRow(Row r, const Constr &c) {
   516       setRow(r,
   517 	     c.lowerBounded()?c.lowerBound():-INF,
   518 	     c.expr(),
   519 	     c.upperBounded()?c.upperBound():INF);
   520     }
   521 
   522     ///Add a new row (i.e a new constaint) to the LP
   523 
   524     ///\param l is the lower bound (-\ref INF means no bound)
   525     ///\param e is a linear expression (see \ref Expr)
   526     ///\param u is the upper bound (\ref INF means no bound)
   527     ///\return The created row.
   528     ///\bug This is a temportary function. The interface will change to
   529     ///a better one.
   530     Row addRow(Value l,const Expr &e, Value u) {
   531       Row r=addRow();
   532       setRow(r,l,e,u);
   533       return r;
   534     }
   535 
   536     ///Add a new row (i.e a new constaint) to the LP
   537 
   538     ///\param c is a linear expression (see \ref Constr)
   539     ///\return The created row.
   540     Row addRow(const Constr &c) {
   541       Row r=addRow();
   542       setRow(r,c);
   543       return r;
   544     }
   545 
   546     /// Set the lower bound of a column (i.e a variable)
   547 
   548     /// The upper bound of a variable (column) has to be given by an 
   549     /// extended number of type Value, i.e. a finite number of type 
   550     /// Value or -\ref INF.
   551     void colLowerBound(Col c, Value value) {
   552       _setColLowerBound(cols.floatingId(c.id),value);
   553     }
   554     /// Set the upper bound of a column (i.e a variable)
   555 
   556     /// The upper bound of a variable (column) has to be given by an 
   557     /// extended number of type Value, i.e. a finite number of type 
   558     /// Value or \ref INF.
   559     void colUpperBound(Col c, Value value) {
   560       _setColUpperBound(cols.floatingId(c.id),value);
   561     };
   562     /// Set the lower and the upper bounds of a column (i.e a variable)
   563 
   564     /// The lower and the upper bounds of
   565     /// a variable (column) have to be given by an 
   566     /// extended number of type Value, i.e. a finite number of type 
   567     /// Value, -\ref INF or \ref INF.
   568     void colBounds(Col c, Value lower, Value upper) {
   569       _setColLowerBound(cols.floatingId(c.id),lower);
   570       _setColUpperBound(cols.floatingId(c.id),upper);
   571     }
   572     
   573     /// Set the lower bound of a row (i.e a constraint)
   574 
   575     /// The lower bound of a linear expression (row) has to be given by an 
   576     /// extended number of type Value, i.e. a finite number of type 
   577     /// Value or -\ref INF.
   578     void rowLowerBound(Row r, Value value) {
   579       _setRowLowerBound(rows.floatingId(r.id),value);
   580     };
   581     /// Set the upper bound of a row (i.e a constraint)
   582 
   583     /// The upper bound of a linear expression (row) has to be given by an 
   584     /// extended number of type Value, i.e. a finite number of type 
   585     /// Value or \ref INF.
   586     void rowUpperBound(Row r, Value value) {
   587       _setRowUpperBound(rows.floatingId(r.id),value);
   588     };
   589     /// Set the lower and the upper bounds of a row (i.e a variable)
   590 
   591     /// The lower and the upper bounds of
   592     /// a constraint (row) have to be given by an 
   593     /// extended number of type Value, i.e. a finite number of type 
   594     /// Value, -\ref INF or \ref INF.
   595     void rowBounds(Row c, Value lower, Value upper) {
   596       _setRowLowerBound(rows.floatingId(c.id),lower);
   597       _setRowUpperBound(rows.floatingId(c.id),upper);
   598     }
   599     
   600     ///Set an element of the objective function
   601     void objCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
   602     ///Set the objective function
   603     
   604     ///\param e is a linear expression of type \ref Expr.
   605     ///\todo What to do with the constant component?
   606     void setObj(Expr e) {
   607       clearObj();
   608       for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
   609 	objCoeff((*i).first,(*i).second);
   610     }
   611 
   612     ///Maximize
   613     void max() { _setMax(); }
   614     ///Minimize
   615     void min() { _setMin(); }
   616 
   617     
   618     ///@}
   619 
   620 
   621     ///\name Solve the LP
   622 
   623     ///@{
   624 
   625     ///\e
   626     SolveExitStatus solve() { return _solve(); }
   627     
   628     ///@}
   629     
   630     ///\name Obtain the solution
   631 
   632     ///@{
   633 
   634     ///\e
   635     SolutionStatus primalStatus() {
   636       return _getPrimalStatus();
   637     }
   638 
   639     ///\e
   640     Value primal(Col c) { return _getPrimal(cols.floatingId(c.id)); }
   641 
   642     ///\e
   643 
   644     ///\return
   645     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
   646     /// of the primal problem, depending on whether we minimize or maximize.
   647     ///- \ref NAN if no primal solution is found.
   648     ///- The (finite) objective value if an optimal solution is found.
   649     Value primalValue() { return _getPrimalValue();}
   650     ///@}
   651     
   652   };  
   653 
   654   ///\e
   655   
   656   ///\relates LpSolverBase::Expr
   657   ///
   658   inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
   659 				      const LpSolverBase::Expr &b) 
   660   {
   661     LpSolverBase::Expr tmp(a);
   662     tmp+=b; ///\todo Don't STL have some special 'merge' algorithm?
   663     return tmp;
   664   }
   665   ///\e
   666   
   667   ///\relates LpSolverBase::Expr
   668   ///
   669   inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
   670 				      const LpSolverBase::Expr &b) 
   671   {
   672     LpSolverBase::Expr tmp(a);
   673     tmp-=b; ///\todo Don't STL have some special 'merge' algorithm?
   674     return tmp;
   675   }
   676   ///\e
   677   
   678   ///\relates LpSolverBase::Expr
   679   ///
   680   inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
   681 				      const LpSolverBase::Value &b) 
   682   {
   683     LpSolverBase::Expr tmp(a);
   684     tmp*=b; ///\todo Don't STL have some special 'merge' algorithm?
   685     return tmp;
   686   }
   687   
   688   ///\e
   689   
   690   ///\relates LpSolverBase::Expr
   691   ///
   692   inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
   693 				      const LpSolverBase::Expr &b) 
   694   {
   695     LpSolverBase::Expr tmp(b);
   696     tmp*=a; ///\todo Don't STL have some special 'merge' algorithm?
   697     return tmp;
   698   }
   699   ///\e
   700   
   701   ///\relates LpSolverBase::Expr
   702   ///
   703   inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
   704 				      const LpSolverBase::Value &b) 
   705   {
   706     LpSolverBase::Expr tmp(a);
   707     tmp/=b; ///\todo Don't STL have some special 'merge' algorithm?
   708     return tmp;
   709   }
   710   
   711   ///\e
   712   
   713   ///\relates LpSolverBase::Constr
   714   ///
   715   inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
   716 					 const LpSolverBase::Expr &f) 
   717   {
   718     return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
   719   }
   720 
   721   ///\e
   722   
   723   ///\relates LpSolverBase::Constr
   724   ///
   725   inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
   726 					 const LpSolverBase::Expr &f) 
   727   {
   728     return LpSolverBase::Constr(e,f);
   729   }
   730 
   731   ///\e
   732   
   733   ///\relates LpSolverBase::Constr
   734   ///
   735   inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
   736 					 const LpSolverBase::Value &f) 
   737   {
   738     return LpSolverBase::Constr(e,f);
   739   }
   740 
   741   ///\e
   742   
   743   ///\relates LpSolverBase::Constr
   744   ///
   745   inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
   746 					 const LpSolverBase::Expr &f) 
   747   {
   748     return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
   749   }
   750 
   751 
   752   ///\e
   753   
   754   ///\relates LpSolverBase::Constr
   755   ///
   756   inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
   757 					 const LpSolverBase::Expr &f) 
   758   {
   759     return LpSolverBase::Constr(f,e);
   760   }
   761 
   762 
   763   ///\e
   764   
   765   ///\relates LpSolverBase::Constr
   766   ///
   767   inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
   768 					 const LpSolverBase::Value &f) 
   769   {
   770     return LpSolverBase::Constr(f,e);
   771   }
   772 
   773   ///\e
   774   
   775   ///\relates LpSolverBase::Constr
   776   ///
   777   inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
   778 					 const LpSolverBase::Expr &f) 
   779   {
   780     return LpSolverBase::Constr(0,e-f,0);
   781   }
   782 
   783   ///\e
   784   
   785   ///\relates LpSolverBase::Constr
   786   ///
   787   inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
   788 					 const LpSolverBase::Constr&c) 
   789   {
   790     LpSolverBase::Constr tmp(c);
   791     ///\todo Create an own exception type.
   792     if(!isnan(tmp.lowerBound())) throw LogicError();
   793     else tmp.lowerBound()=n;
   794     return tmp;
   795   }
   796   ///\e
   797   
   798   ///\relates LpSolverBase::Constr
   799   ///
   800   inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
   801 					 const LpSolverBase::Value &n)
   802   {
   803     LpSolverBase::Constr tmp(c);
   804     ///\todo Create an own exception type.
   805     if(!isnan(tmp.upperBound())) throw LogicError();
   806     else tmp.upperBound()=n;
   807     return tmp;
   808   }
   809 
   810   ///\e
   811   
   812   ///\relates LpSolverBase::Constr
   813   ///
   814   inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
   815 					 const LpSolverBase::Constr&c) 
   816   {
   817     LpSolverBase::Constr tmp(c);
   818     ///\todo Create an own exception type.
   819     if(!isnan(tmp.upperBound())) throw LogicError();
   820     else tmp.upperBound()=n;
   821     return tmp;
   822   }
   823   ///\e
   824   
   825   ///\relates LpSolverBase::Constr
   826   ///
   827   inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
   828 					 const LpSolverBase::Value &n)
   829   {
   830     LpSolverBase::Constr tmp(c);
   831     ///\todo Create an own exception type.
   832     if(!isnan(tmp.lowerBound())) throw LogicError();
   833     else tmp.lowerBound()=n;
   834     return tmp;
   835   }
   836 
   837 
   838 } //namespace lemon
   839 
   840 #endif //LEMON_LP_BASE_H