lemon/lp_base.h
author Balazs Dezso <deba@inf.elte.hu>
Thu, 04 Mar 2010 15:20:59 +0100
changeset 951 41d7ac528c3a
parent 623 745e182d0139
child 793 e4554cd6b2bf
child 1081 f1398882a928
child 1092 2eebc8f7dca5
permissions -rw-r--r--
Uniforming primal scale to 2 (#314)
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2008
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_LP_BASE_H
    20 #define LEMON_LP_BASE_H
    21 
    22 #include<iostream>
    23 #include<vector>
    24 #include<map>
    25 #include<limits>
    26 #include<lemon/math.h>
    27 
    28 #include<lemon/error.h>
    29 #include<lemon/assert.h>
    30 
    31 #include<lemon/core.h>
    32 #include<lemon/bits/solver_bits.h>
    33 
    34 ///\file
    35 ///\brief The interface of the LP solver interface.
    36 ///\ingroup lp_group
    37 namespace lemon {
    38 
    39   ///Common base class for LP and MIP solvers
    40 
    41   ///Usually this class is not used directly, please use one of the concrete
    42   ///implementations of the solver interface.
    43   ///\ingroup lp_group
    44   class LpBase {
    45 
    46   protected:
    47 
    48     _solver_bits::VarIndex rows;
    49     _solver_bits::VarIndex cols;
    50 
    51   public:
    52 
    53     ///Possible outcomes of an LP solving procedure
    54     enum SolveExitStatus {
    55       /// = 0. It means that the problem has been successfully solved: either
    56       ///an optimal solution has been found or infeasibility/unboundedness
    57       ///has been proved.
    58       SOLVED = 0,
    59       /// = 1. Any other case (including the case when some user specified
    60       ///limit has been exceeded).
    61       UNSOLVED = 1
    62     };
    63 
    64     ///Direction of the optimization
    65     enum Sense {
    66       /// Minimization
    67       MIN,
    68       /// Maximization
    69       MAX
    70     };
    71 
    72     ///Enum for \c messageLevel() parameter
    73     enum MessageLevel {
    74       /// No output (default value).
    75       MESSAGE_NOTHING,
    76       /// Error messages only.
    77       MESSAGE_ERROR,
    78       /// Warnings.
    79       MESSAGE_WARNING,
    80       /// Normal output.
    81       MESSAGE_NORMAL,
    82       /// Verbose output.
    83       MESSAGE_VERBOSE
    84     };
    85     
    86 
    87     ///The floating point type used by the solver
    88     typedef double Value;
    89     ///The infinity constant
    90     static const Value INF;
    91     ///The not a number constant
    92     static const Value NaN;
    93 
    94     friend class Col;
    95     friend class ColIt;
    96     friend class Row;
    97     friend class RowIt;
    98 
    99     ///Refer to a column of the LP.
   100 
   101     ///This type is used to refer to a column of the LP.
   102     ///
   103     ///Its value remains valid and correct even after the addition or erase of
   104     ///other columns.
   105     ///
   106     ///\note This class is similar to other Item types in LEMON, like
   107     ///Node and Arc types in digraph.
   108     class Col {
   109       friend class LpBase;
   110     protected:
   111       int _id;
   112       explicit Col(int id) : _id(id) {}
   113     public:
   114       typedef Value ExprValue;
   115       typedef True LpCol;
   116       /// Default constructor
   117       
   118       /// \warning The default constructor sets the Col to an
   119       /// undefined value.
   120       Col() {}
   121       /// Invalid constructor \& conversion.
   122       
   123       /// This constructor initializes the Col to be invalid.
   124       /// \sa Invalid for more details.      
   125       Col(const Invalid&) : _id(-1) {}
   126       /// Equality operator
   127 
   128       /// Two \ref Col "Col"s are equal if and only if they point to
   129       /// the same LP column or both are invalid.
   130       bool operator==(Col c) const  {return _id == c._id;}
   131       /// Inequality operator
   132 
   133       /// \sa operator==(Col c)
   134       ///
   135       bool operator!=(Col c) const  {return _id != c._id;}
   136       /// Artificial ordering operator.
   137 
   138       /// To allow the use of this object in std::map or similar
   139       /// associative container we require this.
   140       ///
   141       /// \note This operator only have to define some strict ordering of
   142       /// the items; this order has nothing to do with the iteration
   143       /// ordering of the items.
   144       bool operator<(Col c) const  {return _id < c._id;}
   145     };
   146 
   147     ///Iterator for iterate over the columns of an LP problem
   148 
   149     /// Its usage is quite simple, for example you can count the number
   150     /// of columns in an LP \c lp:
   151     ///\code
   152     /// int count=0;
   153     /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
   154     ///\endcode
   155     class ColIt : public Col {
   156       const LpBase *_solver;
   157     public:
   158       /// Default constructor
   159       
   160       /// \warning The default constructor sets the iterator
   161       /// to an undefined value.
   162       ColIt() {}
   163       /// Sets the iterator to the first Col
   164       
   165       /// Sets the iterator to the first Col.
   166       ///
   167       ColIt(const LpBase &solver) : _solver(&solver)
   168       {
   169         _solver->cols.firstItem(_id);
   170       }
   171       /// Invalid constructor \& conversion
   172       
   173       /// Initialize the iterator to be invalid.
   174       /// \sa Invalid for more details.
   175       ColIt(const Invalid&) : Col(INVALID) {}
   176       /// Next column
   177       
   178       /// Assign the iterator to the next column.
   179       ///
   180       ColIt &operator++()
   181       {
   182         _solver->cols.nextItem(_id);
   183         return *this;
   184       }
   185     };
   186 
   187     /// \brief Returns the ID of the column.
   188     static int id(const Col& col) { return col._id; }
   189     /// \brief Returns the column with the given ID.
   190     ///
   191     /// \pre The argument should be a valid column ID in the LP problem.
   192     static Col colFromId(int id) { return Col(id); }
   193 
   194     ///Refer to a row of the LP.
   195 
   196     ///This type is used to refer to a row of the LP.
   197     ///
   198     ///Its value remains valid and correct even after the addition or erase of
   199     ///other rows.
   200     ///
   201     ///\note This class is similar to other Item types in LEMON, like
   202     ///Node and Arc types in digraph.
   203     class Row {
   204       friend class LpBase;
   205     protected:
   206       int _id;
   207       explicit Row(int id) : _id(id) {}
   208     public:
   209       typedef Value ExprValue;
   210       typedef True LpRow;
   211       /// Default constructor
   212       
   213       /// \warning The default constructor sets the Row to an
   214       /// undefined value.
   215       Row() {}
   216       /// Invalid constructor \& conversion.
   217       
   218       /// This constructor initializes the Row to be invalid.
   219       /// \sa Invalid for more details.      
   220       Row(const Invalid&) : _id(-1) {}
   221       /// Equality operator
   222 
   223       /// Two \ref Row "Row"s are equal if and only if they point to
   224       /// the same LP row or both are invalid.
   225       bool operator==(Row r) const  {return _id == r._id;}
   226       /// Inequality operator
   227       
   228       /// \sa operator==(Row r)
   229       ///
   230       bool operator!=(Row r) const  {return _id != r._id;}
   231       /// Artificial ordering operator.
   232 
   233       /// To allow the use of this object in std::map or similar
   234       /// associative container we require this.
   235       ///
   236       /// \note This operator only have to define some strict ordering of
   237       /// the items; this order has nothing to do with the iteration
   238       /// ordering of the items.
   239       bool operator<(Row r) const  {return _id < r._id;}
   240     };
   241 
   242     ///Iterator for iterate over the rows of an LP problem
   243 
   244     /// Its usage is quite simple, for example you can count the number
   245     /// of rows in an LP \c lp:
   246     ///\code
   247     /// int count=0;
   248     /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
   249     ///\endcode
   250     class RowIt : public Row {
   251       const LpBase *_solver;
   252     public:
   253       /// Default constructor
   254       
   255       /// \warning The default constructor sets the iterator
   256       /// to an undefined value.
   257       RowIt() {}
   258       /// Sets the iterator to the first Row
   259       
   260       /// Sets the iterator to the first Row.
   261       ///
   262       RowIt(const LpBase &solver) : _solver(&solver)
   263       {
   264         _solver->rows.firstItem(_id);
   265       }
   266       /// Invalid constructor \& conversion
   267       
   268       /// Initialize the iterator to be invalid.
   269       /// \sa Invalid for more details.
   270       RowIt(const Invalid&) : Row(INVALID) {}
   271       /// Next row
   272       
   273       /// Assign the iterator to the next row.
   274       ///
   275       RowIt &operator++()
   276       {
   277         _solver->rows.nextItem(_id);
   278         return *this;
   279       }
   280     };
   281 
   282     /// \brief Returns the ID of the row.
   283     static int id(const Row& row) { return row._id; }
   284     /// \brief Returns the row with the given ID.
   285     ///
   286     /// \pre The argument should be a valid row ID in the LP problem.
   287     static Row rowFromId(int id) { return Row(id); }
   288 
   289   public:
   290 
   291     ///Linear expression of variables and a constant component
   292 
   293     ///This data structure stores a linear expression of the variables
   294     ///(\ref Col "Col"s) and also has a constant component.
   295     ///
   296     ///There are several ways to access and modify the contents of this
   297     ///container.
   298     ///\code
   299     ///e[v]=5;
   300     ///e[v]+=12;
   301     ///e.erase(v);
   302     ///\endcode
   303     ///or you can also iterate through its elements.
   304     ///\code
   305     ///double s=0;
   306     ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
   307     ///  s+=*i * primal(i);
   308     ///\endcode
   309     ///(This code computes the primal value of the expression).
   310     ///- Numbers (<tt>double</tt>'s)
   311     ///and variables (\ref Col "Col"s) directly convert to an
   312     ///\ref Expr and the usual linear operations are defined, so
   313     ///\code
   314     ///v+w
   315     ///2*v-3.12*(v-w/2)+2
   316     ///v*2.1+(3*v+(v*12+w+6)*3)/2
   317     ///\endcode
   318     ///are valid expressions.
   319     ///The usual assignment operations are also defined.
   320     ///\code
   321     ///e=v+w;
   322     ///e+=2*v-3.12*(v-w/2)+2;
   323     ///e*=3.4;
   324     ///e/=5;
   325     ///\endcode
   326     ///- The constant member can be set and read by dereference
   327     ///  operator (unary *)
   328     ///
   329     ///\code
   330     ///*e=12;
   331     ///double c=*e;
   332     ///\endcode
   333     ///
   334     ///\sa Constr
   335     class Expr {
   336       friend class LpBase;
   337     public:
   338       /// The key type of the expression
   339       typedef LpBase::Col Key;
   340       /// The value type of the expression
   341       typedef LpBase::Value Value;
   342 
   343     protected:
   344       Value const_comp;
   345       std::map<int, Value> comps;
   346 
   347     public:
   348       typedef True SolverExpr;
   349       /// Default constructor
   350       
   351       /// Construct an empty expression, the coefficients and
   352       /// the constant component are initialized to zero.
   353       Expr() : const_comp(0) {}
   354       /// Construct an expression from a column
   355 
   356       /// Construct an expression, which has a term with \c c variable
   357       /// and 1.0 coefficient.
   358       Expr(const Col &c) : const_comp(0) {
   359         typedef std::map<int, Value>::value_type pair_type;
   360         comps.insert(pair_type(id(c), 1));
   361       }
   362       /// Construct an expression from a constant
   363 
   364       /// Construct an expression, which's constant component is \c v.
   365       ///
   366       Expr(const Value &v) : const_comp(v) {}
   367       /// Returns the coefficient of the column
   368       Value operator[](const Col& c) const {
   369         std::map<int, Value>::const_iterator it=comps.find(id(c));
   370         if (it != comps.end()) {
   371           return it->second;
   372         } else {
   373           return 0;
   374         }
   375       }
   376       /// Returns the coefficient of the column
   377       Value& operator[](const Col& c) {
   378         return comps[id(c)];
   379       }
   380       /// Sets the coefficient of the column
   381       void set(const Col &c, const Value &v) {
   382         if (v != 0.0) {
   383           typedef std::map<int, Value>::value_type pair_type;
   384           comps.insert(pair_type(id(c), v));
   385         } else {
   386           comps.erase(id(c));
   387         }
   388       }
   389       /// Returns the constant component of the expression
   390       Value& operator*() { return const_comp; }
   391       /// Returns the constant component of the expression
   392       const Value& operator*() const { return const_comp; }
   393       /// \brief Removes the coefficients which's absolute value does
   394       /// not exceed \c epsilon. It also sets to zero the constant
   395       /// component, if it does not exceed epsilon in absolute value.
   396       void simplify(Value epsilon = 0.0) {
   397         std::map<int, Value>::iterator it=comps.begin();
   398         while (it != comps.end()) {
   399           std::map<int, Value>::iterator jt=it;
   400           ++jt;
   401           if (std::fabs((*it).second) <= epsilon) comps.erase(it);
   402           it=jt;
   403         }
   404         if (std::fabs(const_comp) <= epsilon) const_comp = 0;
   405       }
   406 
   407       void simplify(Value epsilon = 0.0) const {
   408         const_cast<Expr*>(this)->simplify(epsilon);
   409       }
   410 
   411       ///Sets all coefficients and the constant component to 0.
   412       void clear() {
   413         comps.clear();
   414         const_comp=0;
   415       }
   416 
   417       ///Compound assignment
   418       Expr &operator+=(const Expr &e) {
   419         for (std::map<int, Value>::const_iterator it=e.comps.begin();
   420              it!=e.comps.end(); ++it)
   421           comps[it->first]+=it->second;
   422         const_comp+=e.const_comp;
   423         return *this;
   424       }
   425       ///Compound assignment
   426       Expr &operator-=(const Expr &e) {
   427         for (std::map<int, Value>::const_iterator it=e.comps.begin();
   428              it!=e.comps.end(); ++it)
   429           comps[it->first]-=it->second;
   430         const_comp-=e.const_comp;
   431         return *this;
   432       }
   433       ///Multiply with a constant
   434       Expr &operator*=(const Value &v) {
   435         for (std::map<int, Value>::iterator it=comps.begin();
   436              it!=comps.end(); ++it)
   437           it->second*=v;
   438         const_comp*=v;
   439         return *this;
   440       }
   441       ///Division with a constant
   442       Expr &operator/=(const Value &c) {
   443         for (std::map<int, Value>::iterator it=comps.begin();
   444              it!=comps.end(); ++it)
   445           it->second/=c;
   446         const_comp/=c;
   447         return *this;
   448       }
   449 
   450       ///Iterator over the expression
   451       
   452       ///The iterator iterates over the terms of the expression. 
   453       /// 
   454       ///\code
   455       ///double s=0;
   456       ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
   457       ///  s+= *i * primal(i);
   458       ///\endcode
   459       class CoeffIt {
   460       private:
   461 
   462         std::map<int, Value>::iterator _it, _end;
   463 
   464       public:
   465 
   466         /// Sets the iterator to the first term
   467         
   468         /// Sets the iterator to the first term of the expression.
   469         ///
   470         CoeffIt(Expr& e)
   471           : _it(e.comps.begin()), _end(e.comps.end()){}
   472 
   473         /// Convert the iterator to the column of the term
   474         operator Col() const {
   475           return colFromId(_it->first);
   476         }
   477 
   478         /// Returns the coefficient of the term
   479         Value& operator*() { return _it->second; }
   480 
   481         /// Returns the coefficient of the term
   482         const Value& operator*() const { return _it->second; }
   483         /// Next term
   484         
   485         /// Assign the iterator to the next term.
   486         ///
   487         CoeffIt& operator++() { ++_it; return *this; }
   488 
   489         /// Equality operator
   490         bool operator==(Invalid) const { return _it == _end; }
   491         /// Inequality operator
   492         bool operator!=(Invalid) const { return _it != _end; }
   493       };
   494 
   495       /// Const iterator over the expression
   496       
   497       ///The iterator iterates over the terms of the expression. 
   498       /// 
   499       ///\code
   500       ///double s=0;
   501       ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
   502       ///  s+=*i * primal(i);
   503       ///\endcode
   504       class ConstCoeffIt {
   505       private:
   506 
   507         std::map<int, Value>::const_iterator _it, _end;
   508 
   509       public:
   510 
   511         /// Sets the iterator to the first term
   512         
   513         /// Sets the iterator to the first term of the expression.
   514         ///
   515         ConstCoeffIt(const Expr& e)
   516           : _it(e.comps.begin()), _end(e.comps.end()){}
   517 
   518         /// Convert the iterator to the column of the term
   519         operator Col() const {
   520           return colFromId(_it->first);
   521         }
   522 
   523         /// Returns the coefficient of the term
   524         const Value& operator*() const { return _it->second; }
   525 
   526         /// Next term
   527         
   528         /// Assign the iterator to the next term.
   529         ///
   530         ConstCoeffIt& operator++() { ++_it; return *this; }
   531 
   532         /// Equality operator
   533         bool operator==(Invalid) const { return _it == _end; }
   534         /// Inequality operator
   535         bool operator!=(Invalid) const { return _it != _end; }
   536       };
   537 
   538     };
   539 
   540     ///Linear constraint
   541 
   542     ///This data stucture represents a linear constraint in the LP.
   543     ///Basically it is a linear expression with a lower or an upper bound
   544     ///(or both). These parts of the constraint can be obtained by the member
   545     ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
   546     ///respectively.
   547     ///There are two ways to construct a constraint.
   548     ///- You can set the linear expression and the bounds directly
   549     ///  by the functions above.
   550     ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
   551     ///  are defined between expressions, or even between constraints whenever
   552     ///  it makes sense. Therefore if \c e and \c f are linear expressions and
   553     ///  \c s and \c t are numbers, then the followings are valid expressions
   554     ///  and thus they can be used directly e.g. in \ref addRow() whenever
   555     ///  it makes sense.
   556     ///\code
   557     ///  e<=s
   558     ///  e<=f
   559     ///  e==f
   560     ///  s<=e<=t
   561     ///  e>=t
   562     ///\endcode
   563     ///\warning The validity of a constraint is checked only at run
   564     ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
   565     ///compile, but will fail an assertion.
   566     class Constr
   567     {
   568     public:
   569       typedef LpBase::Expr Expr;
   570       typedef Expr::Key Key;
   571       typedef Expr::Value Value;
   572 
   573     protected:
   574       Expr _expr;
   575       Value _lb,_ub;
   576     public:
   577       ///\e
   578       Constr() : _expr(), _lb(NaN), _ub(NaN) {}
   579       ///\e
   580       Constr(Value lb, const Expr &e, Value ub) :
   581         _expr(e), _lb(lb), _ub(ub) {}
   582       Constr(const Expr &e) :
   583         _expr(e), _lb(NaN), _ub(NaN) {}
   584       ///\e
   585       void clear()
   586       {
   587         _expr.clear();
   588         _lb=_ub=NaN;
   589       }
   590 
   591       ///Reference to the linear expression
   592       Expr &expr() { return _expr; }
   593       ///Cont reference to the linear expression
   594       const Expr &expr() const { return _expr; }
   595       ///Reference to the lower bound.
   596 
   597       ///\return
   598       ///- \ref INF "INF": the constraint is lower unbounded.
   599       ///- \ref NaN "NaN": lower bound has not been set.
   600       ///- finite number: the lower bound
   601       Value &lowerBound() { return _lb; }
   602       ///The const version of \ref lowerBound()
   603       const Value &lowerBound() const { return _lb; }
   604       ///Reference to the upper bound.
   605 
   606       ///\return
   607       ///- \ref INF "INF": the constraint is upper unbounded.
   608       ///- \ref NaN "NaN": upper bound has not been set.
   609       ///- finite number: the upper bound
   610       Value &upperBound() { return _ub; }
   611       ///The const version of \ref upperBound()
   612       const Value &upperBound() const { return _ub; }
   613       ///Is the constraint lower bounded?
   614       bool lowerBounded() const {
   615         return _lb != -INF && !isNaN(_lb);
   616       }
   617       ///Is the constraint upper bounded?
   618       bool upperBounded() const {
   619         return _ub != INF && !isNaN(_ub);
   620       }
   621 
   622     };
   623 
   624     ///Linear expression of rows
   625 
   626     ///This data structure represents a column of the matrix,
   627     ///thas is it strores a linear expression of the dual variables
   628     ///(\ref Row "Row"s).
   629     ///
   630     ///There are several ways to access and modify the contents of this
   631     ///container.
   632     ///\code
   633     ///e[v]=5;
   634     ///e[v]+=12;
   635     ///e.erase(v);
   636     ///\endcode
   637     ///or you can also iterate through its elements.
   638     ///\code
   639     ///double s=0;
   640     ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
   641     ///  s+=*i;
   642     ///\endcode
   643     ///(This code computes the sum of all coefficients).
   644     ///- Numbers (<tt>double</tt>'s)
   645     ///and variables (\ref Row "Row"s) directly convert to an
   646     ///\ref DualExpr and the usual linear operations are defined, so
   647     ///\code
   648     ///v+w
   649     ///2*v-3.12*(v-w/2)
   650     ///v*2.1+(3*v+(v*12+w)*3)/2
   651     ///\endcode
   652     ///are valid \ref DualExpr dual expressions.
   653     ///The usual assignment operations are also defined.
   654     ///\code
   655     ///e=v+w;
   656     ///e+=2*v-3.12*(v-w/2);
   657     ///e*=3.4;
   658     ///e/=5;
   659     ///\endcode
   660     ///
   661     ///\sa Expr
   662     class DualExpr {
   663       friend class LpBase;
   664     public:
   665       /// The key type of the expression
   666       typedef LpBase::Row Key;
   667       /// The value type of the expression
   668       typedef LpBase::Value Value;
   669 
   670     protected:
   671       std::map<int, Value> comps;
   672 
   673     public:
   674       typedef True SolverExpr;
   675       /// Default constructor
   676       
   677       /// Construct an empty expression, the coefficients are
   678       /// initialized to zero.
   679       DualExpr() {}
   680       /// Construct an expression from a row
   681 
   682       /// Construct an expression, which has a term with \c r dual
   683       /// variable and 1.0 coefficient.
   684       DualExpr(const Row &r) {
   685         typedef std::map<int, Value>::value_type pair_type;
   686         comps.insert(pair_type(id(r), 1));
   687       }
   688       /// Returns the coefficient of the row
   689       Value operator[](const Row& r) const {
   690         std::map<int, Value>::const_iterator it = comps.find(id(r));
   691         if (it != comps.end()) {
   692           return it->second;
   693         } else {
   694           return 0;
   695         }
   696       }
   697       /// Returns the coefficient of the row
   698       Value& operator[](const Row& r) {
   699         return comps[id(r)];
   700       }
   701       /// Sets the coefficient of the row
   702       void set(const Row &r, const Value &v) {
   703         if (v != 0.0) {
   704           typedef std::map<int, Value>::value_type pair_type;
   705           comps.insert(pair_type(id(r), v));
   706         } else {
   707           comps.erase(id(r));
   708         }
   709       }
   710       /// \brief Removes the coefficients which's absolute value does
   711       /// not exceed \c epsilon. 
   712       void simplify(Value epsilon = 0.0) {
   713         std::map<int, Value>::iterator it=comps.begin();
   714         while (it != comps.end()) {
   715           std::map<int, Value>::iterator jt=it;
   716           ++jt;
   717           if (std::fabs((*it).second) <= epsilon) comps.erase(it);
   718           it=jt;
   719         }
   720       }
   721 
   722       void simplify(Value epsilon = 0.0) const {
   723         const_cast<DualExpr*>(this)->simplify(epsilon);
   724       }
   725 
   726       ///Sets all coefficients to 0.
   727       void clear() {
   728         comps.clear();
   729       }
   730       ///Compound assignment
   731       DualExpr &operator+=(const DualExpr &e) {
   732         for (std::map<int, Value>::const_iterator it=e.comps.begin();
   733              it!=e.comps.end(); ++it)
   734           comps[it->first]+=it->second;
   735         return *this;
   736       }
   737       ///Compound assignment
   738       DualExpr &operator-=(const DualExpr &e) {
   739         for (std::map<int, Value>::const_iterator it=e.comps.begin();
   740              it!=e.comps.end(); ++it)
   741           comps[it->first]-=it->second;
   742         return *this;
   743       }
   744       ///Multiply with a constant
   745       DualExpr &operator*=(const Value &v) {
   746         for (std::map<int, Value>::iterator it=comps.begin();
   747              it!=comps.end(); ++it)
   748           it->second*=v;
   749         return *this;
   750       }
   751       ///Division with a constant
   752       DualExpr &operator/=(const Value &v) {
   753         for (std::map<int, Value>::iterator it=comps.begin();
   754              it!=comps.end(); ++it)
   755           it->second/=v;
   756         return *this;
   757       }
   758 
   759       ///Iterator over the expression
   760       
   761       ///The iterator iterates over the terms of the expression. 
   762       /// 
   763       ///\code
   764       ///double s=0;
   765       ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
   766       ///  s+= *i * dual(i);
   767       ///\endcode
   768       class CoeffIt {
   769       private:
   770 
   771         std::map<int, Value>::iterator _it, _end;
   772 
   773       public:
   774 
   775         /// Sets the iterator to the first term
   776         
   777         /// Sets the iterator to the first term of the expression.
   778         ///
   779         CoeffIt(DualExpr& e)
   780           : _it(e.comps.begin()), _end(e.comps.end()){}
   781 
   782         /// Convert the iterator to the row of the term
   783         operator Row() const {
   784           return rowFromId(_it->first);
   785         }
   786 
   787         /// Returns the coefficient of the term
   788         Value& operator*() { return _it->second; }
   789 
   790         /// Returns the coefficient of the term
   791         const Value& operator*() const { return _it->second; }
   792 
   793         /// Next term
   794         
   795         /// Assign the iterator to the next term.
   796         ///
   797         CoeffIt& operator++() { ++_it; return *this; }
   798 
   799         /// Equality operator
   800         bool operator==(Invalid) const { return _it == _end; }
   801         /// Inequality operator
   802         bool operator!=(Invalid) const { return _it != _end; }
   803       };
   804 
   805       ///Iterator over the expression
   806       
   807       ///The iterator iterates over the terms of the expression. 
   808       /// 
   809       ///\code
   810       ///double s=0;
   811       ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
   812       ///  s+= *i * dual(i);
   813       ///\endcode
   814       class ConstCoeffIt {
   815       private:
   816 
   817         std::map<int, Value>::const_iterator _it, _end;
   818 
   819       public:
   820 
   821         /// Sets the iterator to the first term
   822         
   823         /// Sets the iterator to the first term of the expression.
   824         ///
   825         ConstCoeffIt(const DualExpr& e)
   826           : _it(e.comps.begin()), _end(e.comps.end()){}
   827 
   828         /// Convert the iterator to the row of the term
   829         operator Row() const {
   830           return rowFromId(_it->first);
   831         }
   832 
   833         /// Returns the coefficient of the term
   834         const Value& operator*() const { return _it->second; }
   835 
   836         /// Next term
   837         
   838         /// Assign the iterator to the next term.
   839         ///
   840         ConstCoeffIt& operator++() { ++_it; return *this; }
   841 
   842         /// Equality operator
   843         bool operator==(Invalid) const { return _it == _end; }
   844         /// Inequality operator
   845         bool operator!=(Invalid) const { return _it != _end; }
   846       };
   847     };
   848 
   849 
   850   protected:
   851 
   852     class InsertIterator {
   853     private:
   854 
   855       std::map<int, Value>& _host;
   856       const _solver_bits::VarIndex& _index;
   857 
   858     public:
   859 
   860       typedef std::output_iterator_tag iterator_category;
   861       typedef void difference_type;
   862       typedef void value_type;
   863       typedef void reference;
   864       typedef void pointer;
   865 
   866       InsertIterator(std::map<int, Value>& host,
   867                    const _solver_bits::VarIndex& index)
   868         : _host(host), _index(index) {}
   869 
   870       InsertIterator& operator=(const std::pair<int, Value>& value) {
   871         typedef std::map<int, Value>::value_type pair_type;
   872         _host.insert(pair_type(_index[value.first], value.second));
   873         return *this;
   874       }
   875 
   876       InsertIterator& operator*() { return *this; }
   877       InsertIterator& operator++() { return *this; }
   878       InsertIterator operator++(int) { return *this; }
   879 
   880     };
   881 
   882     class ExprIterator {
   883     private:
   884       std::map<int, Value>::const_iterator _host_it;
   885       const _solver_bits::VarIndex& _index;
   886     public:
   887 
   888       typedef std::bidirectional_iterator_tag iterator_category;
   889       typedef std::ptrdiff_t difference_type;
   890       typedef const std::pair<int, Value> value_type;
   891       typedef value_type reference;
   892 
   893       class pointer {
   894       public:
   895         pointer(value_type& _value) : value(_value) {}
   896         value_type* operator->() { return &value; }
   897       private:
   898         value_type value;
   899       };
   900 
   901       ExprIterator(const std::map<int, Value>::const_iterator& host_it,
   902                    const _solver_bits::VarIndex& index)
   903         : _host_it(host_it), _index(index) {}
   904 
   905       reference operator*() {
   906         return std::make_pair(_index(_host_it->first), _host_it->second);
   907       }
   908 
   909       pointer operator->() {
   910         return pointer(operator*());
   911       }
   912 
   913       ExprIterator& operator++() { ++_host_it; return *this; }
   914       ExprIterator operator++(int) {
   915         ExprIterator tmp(*this); ++_host_it; return tmp;
   916       }
   917 
   918       ExprIterator& operator--() { --_host_it; return *this; }
   919       ExprIterator operator--(int) {
   920         ExprIterator tmp(*this); --_host_it; return tmp;
   921       }
   922 
   923       bool operator==(const ExprIterator& it) const {
   924         return _host_it == it._host_it;
   925       }
   926 
   927       bool operator!=(const ExprIterator& it) const {
   928         return _host_it != it._host_it;
   929       }
   930 
   931     };
   932 
   933   protected:
   934 
   935     //Abstract virtual functions
   936 
   937     virtual int _addColId(int col) { return cols.addIndex(col); }
   938     virtual int _addRowId(int row) { return rows.addIndex(row); }
   939 
   940     virtual void _eraseColId(int col) { cols.eraseIndex(col); }
   941     virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
   942 
   943     virtual int _addCol() = 0;
   944     virtual int _addRow() = 0;
   945 
   946     virtual void _eraseCol(int col) = 0;
   947     virtual void _eraseRow(int row) = 0;
   948 
   949     virtual void _getColName(int col, std::string& name) const = 0;
   950     virtual void _setColName(int col, const std::string& name) = 0;
   951     virtual int _colByName(const std::string& name) const = 0;
   952 
   953     virtual void _getRowName(int row, std::string& name) const = 0;
   954     virtual void _setRowName(int row, const std::string& name) = 0;
   955     virtual int _rowByName(const std::string& name) const = 0;
   956 
   957     virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
   958     virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
   959 
   960     virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
   961     virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
   962 
   963     virtual void _setCoeff(int row, int col, Value value) = 0;
   964     virtual Value _getCoeff(int row, int col) const = 0;
   965 
   966     virtual void _setColLowerBound(int i, Value value) = 0;
   967     virtual Value _getColLowerBound(int i) const = 0;
   968 
   969     virtual void _setColUpperBound(int i, Value value) = 0;
   970     virtual Value _getColUpperBound(int i) const = 0;
   971 
   972     virtual void _setRowLowerBound(int i, Value value) = 0;
   973     virtual Value _getRowLowerBound(int i) const = 0;
   974 
   975     virtual void _setRowUpperBound(int i, Value value) = 0;
   976     virtual Value _getRowUpperBound(int i) const = 0;
   977 
   978     virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
   979     virtual void _getObjCoeffs(InsertIterator b) const = 0;
   980 
   981     virtual void _setObjCoeff(int i, Value obj_coef) = 0;
   982     virtual Value _getObjCoeff(int i) const = 0;
   983 
   984     virtual void _setSense(Sense) = 0;
   985     virtual Sense _getSense() const = 0;
   986 
   987     virtual void _clear() = 0;
   988 
   989     virtual const char* _solverName() const = 0;
   990 
   991     virtual void _messageLevel(MessageLevel level) = 0;
   992 
   993     //Own protected stuff
   994 
   995     //Constant component of the objective function
   996     Value obj_const_comp;
   997 
   998     LpBase() : rows(), cols(), obj_const_comp(0) {}
   999 
  1000   public:
  1001 
  1002     /// Virtual destructor
  1003     virtual ~LpBase() {}
  1004 
  1005     ///Gives back the name of the solver.
  1006     const char* solverName() const {return _solverName();}
  1007 
  1008     ///\name Build Up and Modify the LP
  1009 
  1010     ///@{
  1011 
  1012     ///Add a new empty column (i.e a new variable) to the LP
  1013     Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
  1014 
  1015     ///\brief Adds several new columns (i.e variables) at once
  1016     ///
  1017     ///This magic function takes a container as its argument and fills
  1018     ///its elements with new columns (i.e. variables)
  1019     ///\param t can be
  1020     ///- a standard STL compatible iterable container with
  1021     ///\ref Col as its \c values_type like
  1022     ///\code
  1023     ///std::vector<LpBase::Col>
  1024     ///std::list<LpBase::Col>
  1025     ///\endcode
  1026     ///- a standard STL compatible iterable container with
  1027     ///\ref Col as its \c mapped_type like
  1028     ///\code
  1029     ///std::map<AnyType,LpBase::Col>
  1030     ///\endcode
  1031     ///- an iterable lemon \ref concepts::WriteMap "write map" like
  1032     ///\code
  1033     ///ListGraph::NodeMap<LpBase::Col>
  1034     ///ListGraph::ArcMap<LpBase::Col>
  1035     ///\endcode
  1036     ///\return The number of the created column.
  1037 #ifdef DOXYGEN
  1038     template<class T>
  1039     int addColSet(T &t) { return 0;}
  1040 #else
  1041     template<class T>
  1042     typename enable_if<typename T::value_type::LpCol,int>::type
  1043     addColSet(T &t,dummy<0> = 0) {
  1044       int s=0;
  1045       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
  1046       return s;
  1047     }
  1048     template<class T>
  1049     typename enable_if<typename T::value_type::second_type::LpCol,
  1050                        int>::type
  1051     addColSet(T &t,dummy<1> = 1) {
  1052       int s=0;
  1053       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1054         i->second=addCol();
  1055         s++;
  1056       }
  1057       return s;
  1058     }
  1059     template<class T>
  1060     typename enable_if<typename T::MapIt::Value::LpCol,
  1061                        int>::type
  1062     addColSet(T &t,dummy<2> = 2) {
  1063       int s=0;
  1064       for(typename T::MapIt i(t); i!=INVALID; ++i)
  1065         {
  1066           i.set(addCol());
  1067           s++;
  1068         }
  1069       return s;
  1070     }
  1071 #endif
  1072 
  1073     ///Set a column (i.e a dual constraint) of the LP
  1074 
  1075     ///\param c is the column to be modified
  1076     ///\param e is a dual linear expression (see \ref DualExpr)
  1077     ///a better one.
  1078     void col(Col c, const DualExpr &e) {
  1079       e.simplify();
  1080       _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), rows),
  1081                     ExprIterator(e.comps.end(), rows));
  1082     }
  1083 
  1084     ///Get a column (i.e a dual constraint) of the LP
  1085 
  1086     ///\param c is the column to get
  1087     ///\return the dual expression associated to the column
  1088     DualExpr col(Col c) const {
  1089       DualExpr e;
  1090       _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
  1091       return e;
  1092     }
  1093 
  1094     ///Add a new column to the LP
  1095 
  1096     ///\param e is a dual linear expression (see \ref DualExpr)
  1097     ///\param o is the corresponding component of the objective
  1098     ///function. It is 0 by default.
  1099     ///\return The created column.
  1100     Col addCol(const DualExpr &e, Value o = 0) {
  1101       Col c=addCol();
  1102       col(c,e);
  1103       objCoeff(c,o);
  1104       return c;
  1105     }
  1106 
  1107     ///Add a new empty row (i.e a new constraint) to the LP
  1108 
  1109     ///This function adds a new empty row (i.e a new constraint) to the LP.
  1110     ///\return The created row
  1111     Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
  1112 
  1113     ///\brief Add several new rows (i.e constraints) at once
  1114     ///
  1115     ///This magic function takes a container as its argument and fills
  1116     ///its elements with new row (i.e. variables)
  1117     ///\param t can be
  1118     ///- a standard STL compatible iterable container with
  1119     ///\ref Row as its \c values_type like
  1120     ///\code
  1121     ///std::vector<LpBase::Row>
  1122     ///std::list<LpBase::Row>
  1123     ///\endcode
  1124     ///- a standard STL compatible iterable container with
  1125     ///\ref Row as its \c mapped_type like
  1126     ///\code
  1127     ///std::map<AnyType,LpBase::Row>
  1128     ///\endcode
  1129     ///- an iterable lemon \ref concepts::WriteMap "write map" like
  1130     ///\code
  1131     ///ListGraph::NodeMap<LpBase::Row>
  1132     ///ListGraph::ArcMap<LpBase::Row>
  1133     ///\endcode
  1134     ///\return The number of rows created.
  1135 #ifdef DOXYGEN
  1136     template<class T>
  1137     int addRowSet(T &t) { return 0;}
  1138 #else
  1139     template<class T>
  1140     typename enable_if<typename T::value_type::LpRow,int>::type
  1141     addRowSet(T &t, dummy<0> = 0) {
  1142       int s=0;
  1143       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
  1144       return s;
  1145     }
  1146     template<class T>
  1147     typename enable_if<typename T::value_type::second_type::LpRow, int>::type
  1148     addRowSet(T &t, dummy<1> = 1) {
  1149       int s=0;
  1150       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1151         i->second=addRow();
  1152         s++;
  1153       }
  1154       return s;
  1155     }
  1156     template<class T>
  1157     typename enable_if<typename T::MapIt::Value::LpRow, int>::type
  1158     addRowSet(T &t, dummy<2> = 2) {
  1159       int s=0;
  1160       for(typename T::MapIt i(t); i!=INVALID; ++i)
  1161         {
  1162           i.set(addRow());
  1163           s++;
  1164         }
  1165       return s;
  1166     }
  1167 #endif
  1168 
  1169     ///Set a row (i.e a constraint) of the LP
  1170 
  1171     ///\param r is the row to be modified
  1172     ///\param l is lower bound (-\ref INF means no bound)
  1173     ///\param e is a linear expression (see \ref Expr)
  1174     ///\param u is the upper bound (\ref INF means no bound)
  1175     void row(Row r, Value l, const Expr &e, Value u) {
  1176       e.simplify();
  1177       _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
  1178                     ExprIterator(e.comps.end(), cols));
  1179       _setRowLowerBound(rows(id(r)),l - *e);
  1180       _setRowUpperBound(rows(id(r)),u - *e);
  1181     }
  1182 
  1183     ///Set a row (i.e a constraint) of the LP
  1184 
  1185     ///\param r is the row to be modified
  1186     ///\param c is a linear expression (see \ref Constr)
  1187     void row(Row r, const Constr &c) {
  1188       row(r, c.lowerBounded()?c.lowerBound():-INF,
  1189           c.expr(), c.upperBounded()?c.upperBound():INF);
  1190     }
  1191 
  1192 
  1193     ///Get a row (i.e a constraint) of the LP
  1194 
  1195     ///\param r is the row to get
  1196     ///\return the expression associated to the row
  1197     Expr row(Row r) const {
  1198       Expr e;
  1199       _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
  1200       return e;
  1201     }
  1202 
  1203     ///Add a new row (i.e a new constraint) to the LP
  1204 
  1205     ///\param l is the lower bound (-\ref INF means no bound)
  1206     ///\param e is a linear expression (see \ref Expr)
  1207     ///\param u is the upper bound (\ref INF means no bound)
  1208     ///\return The created row.
  1209     Row addRow(Value l,const Expr &e, Value u) {
  1210       Row r=addRow();
  1211       row(r,l,e,u);
  1212       return r;
  1213     }
  1214 
  1215     ///Add a new row (i.e a new constraint) to the LP
  1216 
  1217     ///\param c is a linear expression (see \ref Constr)
  1218     ///\return The created row.
  1219     Row addRow(const Constr &c) {
  1220       Row r=addRow();
  1221       row(r,c);
  1222       return r;
  1223     }
  1224     ///Erase a column (i.e a variable) from the LP
  1225 
  1226     ///\param c is the column to be deleted
  1227     void erase(Col c) {
  1228       _eraseCol(cols(id(c)));
  1229       _eraseColId(cols(id(c)));
  1230     }
  1231     ///Erase a row (i.e a constraint) from the LP
  1232 
  1233     ///\param r is the row to be deleted
  1234     void erase(Row r) {
  1235       _eraseRow(rows(id(r)));
  1236       _eraseRowId(rows(id(r)));
  1237     }
  1238 
  1239     /// Get the name of a column
  1240 
  1241     ///\param c is the coresponding column
  1242     ///\return The name of the colunm
  1243     std::string colName(Col c) const {
  1244       std::string name;
  1245       _getColName(cols(id(c)), name);
  1246       return name;
  1247     }
  1248 
  1249     /// Set the name of a column
  1250 
  1251     ///\param c is the coresponding column
  1252     ///\param name The name to be given
  1253     void colName(Col c, const std::string& name) {
  1254       _setColName(cols(id(c)), name);
  1255     }
  1256 
  1257     /// Get the column by its name
  1258 
  1259     ///\param name The name of the column
  1260     ///\return the proper column or \c INVALID
  1261     Col colByName(const std::string& name) const {
  1262       int k = _colByName(name);
  1263       return k != -1 ? Col(cols[k]) : Col(INVALID);
  1264     }
  1265 
  1266     /// Get the name of a row
  1267 
  1268     ///\param r is the coresponding row
  1269     ///\return The name of the row
  1270     std::string rowName(Row r) const {
  1271       std::string name;
  1272       _getRowName(rows(id(r)), name);
  1273       return name;
  1274     }
  1275 
  1276     /// Set the name of a row
  1277 
  1278     ///\param r is the coresponding row
  1279     ///\param name The name to be given
  1280     void rowName(Row r, const std::string& name) {
  1281       _setRowName(rows(id(r)), name);
  1282     }
  1283 
  1284     /// Get the row by its name
  1285 
  1286     ///\param name The name of the row
  1287     ///\return the proper row or \c INVALID
  1288     Row rowByName(const std::string& name) const {
  1289       int k = _rowByName(name);
  1290       return k != -1 ? Row(rows[k]) : Row(INVALID);
  1291     }
  1292 
  1293     /// Set an element of the coefficient matrix of the LP
  1294 
  1295     ///\param r is the row of the element to be modified
  1296     ///\param c is the column of the element to be modified
  1297     ///\param val is the new value of the coefficient
  1298     void coeff(Row r, Col c, Value val) {
  1299       _setCoeff(rows(id(r)),cols(id(c)), val);
  1300     }
  1301 
  1302     /// Get an element of the coefficient matrix of the LP
  1303 
  1304     ///\param r is the row of the element
  1305     ///\param c is the column of the element
  1306     ///\return the corresponding coefficient
  1307     Value coeff(Row r, Col c) const {
  1308       return _getCoeff(rows(id(r)),cols(id(c)));
  1309     }
  1310 
  1311     /// Set the lower bound of a column (i.e a variable)
  1312 
  1313     /// The lower bound of a variable (column) has to be given by an
  1314     /// extended number of type Value, i.e. a finite number of type
  1315     /// Value or -\ref INF.
  1316     void colLowerBound(Col c, Value value) {
  1317       _setColLowerBound(cols(id(c)),value);
  1318     }
  1319 
  1320     /// Get the lower bound of a column (i.e a variable)
  1321 
  1322     /// This function returns the lower bound for column (variable) \c c
  1323     /// (this might be -\ref INF as well).
  1324     ///\return The lower bound for column \c c
  1325     Value colLowerBound(Col c) const {
  1326       return _getColLowerBound(cols(id(c)));
  1327     }
  1328 
  1329     ///\brief Set the lower bound of  several columns
  1330     ///(i.e variables) at once
  1331     ///
  1332     ///This magic function takes a container as its argument
  1333     ///and applies the function on all of its elements.
  1334     ///The lower bound of a variable (column) has to be given by an
  1335     ///extended number of type Value, i.e. a finite number of type
  1336     ///Value or -\ref INF.
  1337 #ifdef DOXYGEN
  1338     template<class T>
  1339     void colLowerBound(T &t, Value value) { return 0;}
  1340 #else
  1341     template<class T>
  1342     typename enable_if<typename T::value_type::LpCol,void>::type
  1343     colLowerBound(T &t, Value value,dummy<0> = 0) {
  1344       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1345         colLowerBound(*i, value);
  1346       }
  1347     }
  1348     template<class T>
  1349     typename enable_if<typename T::value_type::second_type::LpCol,
  1350                        void>::type
  1351     colLowerBound(T &t, Value value,dummy<1> = 1) {
  1352       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1353         colLowerBound(i->second, value);
  1354       }
  1355     }
  1356     template<class T>
  1357     typename enable_if<typename T::MapIt::Value::LpCol,
  1358                        void>::type
  1359     colLowerBound(T &t, Value value,dummy<2> = 2) {
  1360       for(typename T::MapIt i(t); i!=INVALID; ++i){
  1361         colLowerBound(*i, value);
  1362       }
  1363     }
  1364 #endif
  1365 
  1366     /// Set the upper bound of a column (i.e a variable)
  1367 
  1368     /// The upper bound of a variable (column) has to be given by an
  1369     /// extended number of type Value, i.e. a finite number of type
  1370     /// Value or \ref INF.
  1371     void colUpperBound(Col c, Value value) {
  1372       _setColUpperBound(cols(id(c)),value);
  1373     };
  1374 
  1375     /// Get the upper bound of a column (i.e a variable)
  1376 
  1377     /// This function returns the upper bound for column (variable) \c c
  1378     /// (this might be \ref INF as well).
  1379     /// \return The upper bound for column \c c
  1380     Value colUpperBound(Col c) const {
  1381       return _getColUpperBound(cols(id(c)));
  1382     }
  1383 
  1384     ///\brief Set the upper bound of  several columns
  1385     ///(i.e variables) at once
  1386     ///
  1387     ///This magic function takes a container as its argument
  1388     ///and applies the function on all of its elements.
  1389     ///The upper bound of a variable (column) has to be given by an
  1390     ///extended number of type Value, i.e. a finite number of type
  1391     ///Value or \ref INF.
  1392 #ifdef DOXYGEN
  1393     template<class T>
  1394     void colUpperBound(T &t, Value value) { return 0;}
  1395 #else
  1396     template<class T1>
  1397     typename enable_if<typename T1::value_type::LpCol,void>::type
  1398     colUpperBound(T1 &t, Value value,dummy<0> = 0) {
  1399       for(typename T1::iterator i=t.begin();i!=t.end();++i) {
  1400         colUpperBound(*i, value);
  1401       }
  1402     }
  1403     template<class T1>
  1404     typename enable_if<typename T1::value_type::second_type::LpCol,
  1405                        void>::type
  1406     colUpperBound(T1 &t, Value value,dummy<1> = 1) {
  1407       for(typename T1::iterator i=t.begin();i!=t.end();++i) {
  1408         colUpperBound(i->second, value);
  1409       }
  1410     }
  1411     template<class T1>
  1412     typename enable_if<typename T1::MapIt::Value::LpCol,
  1413                        void>::type
  1414     colUpperBound(T1 &t, Value value,dummy<2> = 2) {
  1415       for(typename T1::MapIt i(t); i!=INVALID; ++i){
  1416         colUpperBound(*i, value);
  1417       }
  1418     }
  1419 #endif
  1420 
  1421     /// Set the lower and the upper bounds of a column (i.e a variable)
  1422 
  1423     /// The lower and the upper bounds of
  1424     /// a variable (column) have to be given by an
  1425     /// extended number of type Value, i.e. a finite number of type
  1426     /// Value, -\ref INF or \ref INF.
  1427     void colBounds(Col c, Value lower, Value upper) {
  1428       _setColLowerBound(cols(id(c)),lower);
  1429       _setColUpperBound(cols(id(c)),upper);
  1430     }
  1431 
  1432     ///\brief Set the lower and the upper bound of several columns
  1433     ///(i.e variables) at once
  1434     ///
  1435     ///This magic function takes a container as its argument
  1436     ///and applies the function on all of its elements.
  1437     /// The lower and the upper bounds of
  1438     /// a variable (column) have to be given by an
  1439     /// extended number of type Value, i.e. a finite number of type
  1440     /// Value, -\ref INF or \ref INF.
  1441 #ifdef DOXYGEN
  1442     template<class T>
  1443     void colBounds(T &t, Value lower, Value upper) { return 0;}
  1444 #else
  1445     template<class T2>
  1446     typename enable_if<typename T2::value_type::LpCol,void>::type
  1447     colBounds(T2 &t, Value lower, Value upper,dummy<0> = 0) {
  1448       for(typename T2::iterator i=t.begin();i!=t.end();++i) {
  1449         colBounds(*i, lower, upper);
  1450       }
  1451     }
  1452     template<class T2>
  1453     typename enable_if<typename T2::value_type::second_type::LpCol, void>::type
  1454     colBounds(T2 &t, Value lower, Value upper,dummy<1> = 1) {
  1455       for(typename T2::iterator i=t.begin();i!=t.end();++i) {
  1456         colBounds(i->second, lower, upper);
  1457       }
  1458     }
  1459     template<class T2>
  1460     typename enable_if<typename T2::MapIt::Value::LpCol, void>::type
  1461     colBounds(T2 &t, Value lower, Value upper,dummy<2> = 2) {
  1462       for(typename T2::MapIt i(t); i!=INVALID; ++i){
  1463         colBounds(*i, lower, upper);
  1464       }
  1465     }
  1466 #endif
  1467 
  1468     /// Set the lower bound of a row (i.e a constraint)
  1469 
  1470     /// The lower bound of a constraint (row) has to be given by an
  1471     /// extended number of type Value, i.e. a finite number of type
  1472     /// Value or -\ref INF.
  1473     void rowLowerBound(Row r, Value value) {
  1474       _setRowLowerBound(rows(id(r)),value);
  1475     }
  1476 
  1477     /// Get the lower bound of a row (i.e a constraint)
  1478 
  1479     /// This function returns the lower bound for row (constraint) \c c
  1480     /// (this might be -\ref INF as well).
  1481     ///\return The lower bound for row \c r
  1482     Value rowLowerBound(Row r) const {
  1483       return _getRowLowerBound(rows(id(r)));
  1484     }
  1485 
  1486     /// Set the upper bound of a row (i.e a constraint)
  1487 
  1488     /// The upper bound of a constraint (row) has to be given by an
  1489     /// extended number of type Value, i.e. a finite number of type
  1490     /// Value or -\ref INF.
  1491     void rowUpperBound(Row r, Value value) {
  1492       _setRowUpperBound(rows(id(r)),value);
  1493     }
  1494 
  1495     /// Get the upper bound of a row (i.e a constraint)
  1496 
  1497     /// This function returns the upper bound for row (constraint) \c c
  1498     /// (this might be -\ref INF as well).
  1499     ///\return The upper bound for row \c r
  1500     Value rowUpperBound(Row r) const {
  1501       return _getRowUpperBound(rows(id(r)));
  1502     }
  1503 
  1504     ///Set an element of the objective function
  1505     void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
  1506 
  1507     ///Get an element of the objective function
  1508     Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
  1509 
  1510     ///Set the objective function
  1511 
  1512     ///\param e is a linear expression of type \ref Expr.
  1513     ///
  1514     void obj(const Expr& e) {
  1515       _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
  1516                     ExprIterator(e.comps.end(), cols));
  1517       obj_const_comp = *e;
  1518     }
  1519 
  1520     ///Get the objective function
  1521 
  1522     ///\return the objective function as a linear expression of type
  1523     ///Expr.
  1524     Expr obj() const {
  1525       Expr e;
  1526       _getObjCoeffs(InsertIterator(e.comps, cols));
  1527       *e = obj_const_comp;
  1528       return e;
  1529     }
  1530 
  1531 
  1532     ///Set the direction of optimization
  1533     void sense(Sense sense) { _setSense(sense); }
  1534 
  1535     ///Query the direction of the optimization
  1536     Sense sense() const {return _getSense(); }
  1537 
  1538     ///Set the sense to maximization
  1539     void max() { _setSense(MAX); }
  1540 
  1541     ///Set the sense to maximization
  1542     void min() { _setSense(MIN); }
  1543 
  1544     ///Clears the problem
  1545     void clear() { _clear(); }
  1546 
  1547     /// Sets the message level of the solver
  1548     void messageLevel(MessageLevel level) { _messageLevel(level); }
  1549 
  1550     ///@}
  1551 
  1552   };
  1553 
  1554   /// Addition
  1555 
  1556   ///\relates LpBase::Expr
  1557   ///
  1558   inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
  1559     LpBase::Expr tmp(a);
  1560     tmp+=b;
  1561     return tmp;
  1562   }
  1563   ///Substraction
  1564 
  1565   ///\relates LpBase::Expr
  1566   ///
  1567   inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
  1568     LpBase::Expr tmp(a);
  1569     tmp-=b;
  1570     return tmp;
  1571   }
  1572   ///Multiply with constant
  1573 
  1574   ///\relates LpBase::Expr
  1575   ///
  1576   inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
  1577     LpBase::Expr tmp(a);
  1578     tmp*=b;
  1579     return tmp;
  1580   }
  1581 
  1582   ///Multiply with constant
  1583 
  1584   ///\relates LpBase::Expr
  1585   ///
  1586   inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
  1587     LpBase::Expr tmp(b);
  1588     tmp*=a;
  1589     return tmp;
  1590   }
  1591   ///Divide with constant
  1592 
  1593   ///\relates LpBase::Expr
  1594   ///
  1595   inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
  1596     LpBase::Expr tmp(a);
  1597     tmp/=b;
  1598     return tmp;
  1599   }
  1600 
  1601   ///Create constraint
  1602 
  1603   ///\relates LpBase::Constr
  1604   ///
  1605   inline LpBase::Constr operator<=(const LpBase::Expr &e,
  1606                                    const LpBase::Expr &f) {
  1607     return LpBase::Constr(0, f - e, LpBase::INF);
  1608   }
  1609 
  1610   ///Create constraint
  1611 
  1612   ///\relates LpBase::Constr
  1613   ///
  1614   inline LpBase::Constr operator<=(const LpBase::Value &e,
  1615                                    const LpBase::Expr &f) {
  1616     return LpBase::Constr(e, f, LpBase::NaN);
  1617   }
  1618 
  1619   ///Create constraint
  1620 
  1621   ///\relates LpBase::Constr
  1622   ///
  1623   inline LpBase::Constr operator<=(const LpBase::Expr &e,
  1624                                    const LpBase::Value &f) {
  1625     return LpBase::Constr(- LpBase::INF, e, f);
  1626   }
  1627 
  1628   ///Create constraint
  1629 
  1630   ///\relates LpBase::Constr
  1631   ///
  1632   inline LpBase::Constr operator>=(const LpBase::Expr &e,
  1633                                    const LpBase::Expr &f) {
  1634     return LpBase::Constr(0, e - f, LpBase::INF);
  1635   }
  1636 
  1637 
  1638   ///Create constraint
  1639 
  1640   ///\relates LpBase::Constr
  1641   ///
  1642   inline LpBase::Constr operator>=(const LpBase::Value &e,
  1643                                    const LpBase::Expr &f) {
  1644     return LpBase::Constr(LpBase::NaN, f, e);
  1645   }
  1646 
  1647 
  1648   ///Create constraint
  1649 
  1650   ///\relates LpBase::Constr
  1651   ///
  1652   inline LpBase::Constr operator>=(const LpBase::Expr &e,
  1653                                    const LpBase::Value &f) {
  1654     return LpBase::Constr(f, e, LpBase::INF);
  1655   }
  1656 
  1657   ///Create constraint
  1658 
  1659   ///\relates LpBase::Constr
  1660   ///
  1661   inline LpBase::Constr operator==(const LpBase::Expr &e,
  1662                                    const LpBase::Value &f) {
  1663     return LpBase::Constr(f, e, f);
  1664   }
  1665 
  1666   ///Create constraint
  1667 
  1668   ///\relates LpBase::Constr
  1669   ///
  1670   inline LpBase::Constr operator==(const LpBase::Expr &e,
  1671                                    const LpBase::Expr &f) {
  1672     return LpBase::Constr(0, f - e, 0);
  1673   }
  1674 
  1675   ///Create constraint
  1676 
  1677   ///\relates LpBase::Constr
  1678   ///
  1679   inline LpBase::Constr operator<=(const LpBase::Value &n,
  1680                                    const LpBase::Constr &c) {
  1681     LpBase::Constr tmp(c);
  1682     LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
  1683     tmp.lowerBound()=n;
  1684     return tmp;
  1685   }
  1686   ///Create constraint
  1687 
  1688   ///\relates LpBase::Constr
  1689   ///
  1690   inline LpBase::Constr operator<=(const LpBase::Constr &c,
  1691                                    const LpBase::Value &n)
  1692   {
  1693     LpBase::Constr tmp(c);
  1694     LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
  1695     tmp.upperBound()=n;
  1696     return tmp;
  1697   }
  1698 
  1699   ///Create constraint
  1700 
  1701   ///\relates LpBase::Constr
  1702   ///
  1703   inline LpBase::Constr operator>=(const LpBase::Value &n,
  1704                                    const LpBase::Constr &c) {
  1705     LpBase::Constr tmp(c);
  1706     LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
  1707     tmp.upperBound()=n;
  1708     return tmp;
  1709   }
  1710   ///Create constraint
  1711 
  1712   ///\relates LpBase::Constr
  1713   ///
  1714   inline LpBase::Constr operator>=(const LpBase::Constr &c,
  1715                                    const LpBase::Value &n)
  1716   {
  1717     LpBase::Constr tmp(c);
  1718     LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
  1719     tmp.lowerBound()=n;
  1720     return tmp;
  1721   }
  1722 
  1723   ///Addition
  1724 
  1725   ///\relates LpBase::DualExpr
  1726   ///
  1727   inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
  1728                                     const LpBase::DualExpr &b) {
  1729     LpBase::DualExpr tmp(a);
  1730     tmp+=b;
  1731     return tmp;
  1732   }
  1733   ///Substraction
  1734 
  1735   ///\relates LpBase::DualExpr
  1736   ///
  1737   inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
  1738                                     const LpBase::DualExpr &b) {
  1739     LpBase::DualExpr tmp(a);
  1740     tmp-=b;
  1741     return tmp;
  1742   }
  1743   ///Multiply with constant
  1744 
  1745   ///\relates LpBase::DualExpr
  1746   ///
  1747   inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
  1748                                     const LpBase::Value &b) {
  1749     LpBase::DualExpr tmp(a);
  1750     tmp*=b;
  1751     return tmp;
  1752   }
  1753 
  1754   ///Multiply with constant
  1755 
  1756   ///\relates LpBase::DualExpr
  1757   ///
  1758   inline LpBase::DualExpr operator*(const LpBase::Value &a,
  1759                                     const LpBase::DualExpr &b) {
  1760     LpBase::DualExpr tmp(b);
  1761     tmp*=a;
  1762     return tmp;
  1763   }
  1764   ///Divide with constant
  1765 
  1766   ///\relates LpBase::DualExpr
  1767   ///
  1768   inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
  1769                                     const LpBase::Value &b) {
  1770     LpBase::DualExpr tmp(a);
  1771     tmp/=b;
  1772     return tmp;
  1773   }
  1774 
  1775   /// \ingroup lp_group
  1776   ///
  1777   /// \brief Common base class for LP solvers
  1778   ///
  1779   /// This class is an abstract base class for LP solvers. This class
  1780   /// provides a full interface for set and modify an LP problem,
  1781   /// solve it and retrieve the solution. You can use one of the
  1782   /// descendants as a concrete implementation, or the \c Lp
  1783   /// default LP solver. However, if you would like to handle LP
  1784   /// solvers as reference or pointer in a generic way, you can use
  1785   /// this class directly.
  1786   class LpSolver : virtual public LpBase {
  1787   public:
  1788 
  1789     /// The problem types for primal and dual problems
  1790     enum ProblemType {
  1791       /// = 0. Feasible solution hasn't been found (but may exist).
  1792       UNDEFINED = 0,
  1793       /// = 1. The problem has no feasible solution.
  1794       INFEASIBLE = 1,
  1795       /// = 2. Feasible solution found.
  1796       FEASIBLE = 2,
  1797       /// = 3. Optimal solution exists and found.
  1798       OPTIMAL = 3,
  1799       /// = 4. The cost function is unbounded.
  1800       UNBOUNDED = 4
  1801     };
  1802 
  1803     ///The basis status of variables
  1804     enum VarStatus {
  1805       /// The variable is in the basis
  1806       BASIC, 
  1807       /// The variable is free, but not basic
  1808       FREE,
  1809       /// The variable has active lower bound 
  1810       LOWER,
  1811       /// The variable has active upper bound
  1812       UPPER,
  1813       /// The variable is non-basic and fixed
  1814       FIXED
  1815     };
  1816 
  1817   protected:
  1818 
  1819     virtual SolveExitStatus _solve() = 0;
  1820 
  1821     virtual Value _getPrimal(int i) const = 0;
  1822     virtual Value _getDual(int i) const = 0;
  1823 
  1824     virtual Value _getPrimalRay(int i) const = 0;
  1825     virtual Value _getDualRay(int i) const = 0;
  1826 
  1827     virtual Value _getPrimalValue() const = 0;
  1828 
  1829     virtual VarStatus _getColStatus(int i) const = 0;
  1830     virtual VarStatus _getRowStatus(int i) const = 0;
  1831 
  1832     virtual ProblemType _getPrimalType() const = 0;
  1833     virtual ProblemType _getDualType() const = 0;
  1834 
  1835   public:
  1836 
  1837     ///Allocate a new LP problem instance
  1838     virtual LpSolver* newSolver() const = 0;
  1839     ///Make a copy of the LP problem
  1840     virtual LpSolver* cloneSolver() const = 0;
  1841 
  1842     ///\name Solve the LP
  1843 
  1844     ///@{
  1845 
  1846     ///\e Solve the LP problem at hand
  1847     ///
  1848     ///\return The result of the optimization procedure. Possible
  1849     ///values and their meanings can be found in the documentation of
  1850     ///\ref SolveExitStatus.
  1851     SolveExitStatus solve() { return _solve(); }
  1852 
  1853     ///@}
  1854 
  1855     ///\name Obtain the Solution
  1856 
  1857     ///@{
  1858 
  1859     /// The type of the primal problem
  1860     ProblemType primalType() const {
  1861       return _getPrimalType();
  1862     }
  1863 
  1864     /// The type of the dual problem
  1865     ProblemType dualType() const {
  1866       return _getDualType();
  1867     }
  1868 
  1869     /// Return the primal value of the column
  1870 
  1871     /// Return the primal value of the column.
  1872     /// \pre The problem is solved.
  1873     Value primal(Col c) const { return _getPrimal(cols(id(c))); }
  1874 
  1875     /// Return the primal value of the expression
  1876 
  1877     /// Return the primal value of the expression, i.e. the dot
  1878     /// product of the primal solution and the expression.
  1879     /// \pre The problem is solved.
  1880     Value primal(const Expr& e) const {
  1881       double res = *e;
  1882       for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
  1883         res += *c * primal(c);
  1884       }
  1885       return res;
  1886     }
  1887     /// Returns a component of the primal ray
  1888     
  1889     /// The primal ray is solution of the modified primal problem,
  1890     /// where we change each finite bound to 0, and we looking for a
  1891     /// negative objective value in case of minimization, and positive
  1892     /// objective value for maximization. If there is such solution,
  1893     /// that proofs the unsolvability of the dual problem, and if a
  1894     /// feasible primal solution exists, then the unboundness of
  1895     /// primal problem.
  1896     ///
  1897     /// \pre The problem is solved and the dual problem is infeasible.
  1898     /// \note Some solvers does not provide primal ray calculation
  1899     /// functions.
  1900     Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
  1901 
  1902     /// Return the dual value of the row
  1903 
  1904     /// Return the dual value of the row.
  1905     /// \pre The problem is solved.
  1906     Value dual(Row r) const { return _getDual(rows(id(r))); }
  1907 
  1908     /// Return the dual value of the dual expression
  1909 
  1910     /// Return the dual value of the dual expression, i.e. the dot
  1911     /// product of the dual solution and the dual expression.
  1912     /// \pre The problem is solved.
  1913     Value dual(const DualExpr& e) const {
  1914       double res = 0.0;
  1915       for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
  1916         res += *r * dual(r);
  1917       }
  1918       return res;
  1919     }
  1920 
  1921     /// Returns a component of the dual ray
  1922     
  1923     /// The dual ray is solution of the modified primal problem, where
  1924     /// we change each finite bound to 0 (i.e. the objective function
  1925     /// coefficients in the primal problem), and we looking for a
  1926     /// ositive objective value. If there is such solution, that
  1927     /// proofs the unsolvability of the primal problem, and if a
  1928     /// feasible dual solution exists, then the unboundness of
  1929     /// dual problem.
  1930     ///
  1931     /// \pre The problem is solved and the primal problem is infeasible.
  1932     /// \note Some solvers does not provide dual ray calculation
  1933     /// functions.
  1934     Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
  1935 
  1936     /// Return the basis status of the column
  1937 
  1938     /// \see VarStatus
  1939     VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
  1940 
  1941     /// Return the basis status of the row
  1942 
  1943     /// \see VarStatus
  1944     VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
  1945 
  1946     ///The value of the objective function
  1947 
  1948     ///\return
  1949     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
  1950     /// of the primal problem, depending on whether we minimize or maximize.
  1951     ///- \ref NaN if no primal solution is found.
  1952     ///- The (finite) objective value if an optimal solution is found.
  1953     Value primal() const { return _getPrimalValue()+obj_const_comp;}
  1954     ///@}
  1955 
  1956   protected:
  1957 
  1958   };
  1959 
  1960 
  1961   /// \ingroup lp_group
  1962   ///
  1963   /// \brief Common base class for MIP solvers
  1964   ///
  1965   /// This class is an abstract base class for MIP solvers. This class
  1966   /// provides a full interface for set and modify an MIP problem,
  1967   /// solve it and retrieve the solution. You can use one of the
  1968   /// descendants as a concrete implementation, or the \c Lp
  1969   /// default MIP solver. However, if you would like to handle MIP
  1970   /// solvers as reference or pointer in a generic way, you can use
  1971   /// this class directly.
  1972   class MipSolver : virtual public LpBase {
  1973   public:
  1974 
  1975     /// The problem types for MIP problems
  1976     enum ProblemType {
  1977       /// = 0. Feasible solution hasn't been found (but may exist).
  1978       UNDEFINED = 0,
  1979       /// = 1. The problem has no feasible solution.
  1980       INFEASIBLE = 1,
  1981       /// = 2. Feasible solution found.
  1982       FEASIBLE = 2,
  1983       /// = 3. Optimal solution exists and found.
  1984       OPTIMAL = 3,
  1985       /// = 4. The cost function is unbounded.
  1986       ///The Mip or at least the relaxed problem is unbounded.
  1987       UNBOUNDED = 4
  1988     };
  1989 
  1990     ///Allocate a new MIP problem instance
  1991     virtual MipSolver* newSolver() const = 0;
  1992     ///Make a copy of the MIP problem
  1993     virtual MipSolver* cloneSolver() const = 0;
  1994 
  1995     ///\name Solve the MIP
  1996 
  1997     ///@{
  1998 
  1999     /// Solve the MIP problem at hand
  2000     ///
  2001     ///\return The result of the optimization procedure. Possible
  2002     ///values and their meanings can be found in the documentation of
  2003     ///\ref SolveExitStatus.
  2004     SolveExitStatus solve() { return _solve(); }
  2005 
  2006     ///@}
  2007 
  2008     ///\name Set Column Type
  2009     ///@{
  2010 
  2011     ///Possible variable (column) types (e.g. real, integer, binary etc.)
  2012     enum ColTypes {
  2013       /// = 0. Continuous variable (default).
  2014       REAL = 0,
  2015       /// = 1. Integer variable.
  2016       INTEGER = 1
  2017     };
  2018 
  2019     ///Sets the type of the given column to the given type
  2020 
  2021     ///Sets the type of the given column to the given type.
  2022     ///
  2023     void colType(Col c, ColTypes col_type) {
  2024       _setColType(cols(id(c)),col_type);
  2025     }
  2026 
  2027     ///Gives back the type of the column.
  2028 
  2029     ///Gives back the type of the column.
  2030     ///
  2031     ColTypes colType(Col c) const {
  2032       return _getColType(cols(id(c)));
  2033     }
  2034     ///@}
  2035 
  2036     ///\name Obtain the Solution
  2037 
  2038     ///@{
  2039 
  2040     /// The type of the MIP problem
  2041     ProblemType type() const {
  2042       return _getType();
  2043     }
  2044 
  2045     /// Return the value of the row in the solution
  2046 
  2047     ///  Return the value of the row in the solution.
  2048     /// \pre The problem is solved.
  2049     Value sol(Col c) const { return _getSol(cols(id(c))); }
  2050 
  2051     /// Return the value of the expression in the solution
  2052 
  2053     /// Return the value of the expression in the solution, i.e. the
  2054     /// dot product of the solution and the expression.
  2055     /// \pre The problem is solved.
  2056     Value sol(const Expr& e) const {
  2057       double res = *e;
  2058       for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
  2059         res += *c * sol(c);
  2060       }
  2061       return res;
  2062     }
  2063     ///The value of the objective function
  2064     
  2065     ///\return
  2066     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
  2067     /// of the problem, depending on whether we minimize or maximize.
  2068     ///- \ref NaN if no primal solution is found.
  2069     ///- The (finite) objective value if an optimal solution is found.
  2070     Value solValue() const { return _getSolValue()+obj_const_comp;}
  2071     ///@}
  2072 
  2073   protected:
  2074 
  2075     virtual SolveExitStatus _solve() = 0;
  2076     virtual ColTypes _getColType(int col) const = 0;
  2077     virtual void _setColType(int col, ColTypes col_type) = 0;
  2078     virtual ProblemType _getType() const = 0;
  2079     virtual Value _getSol(int i) const = 0;
  2080     virtual Value _getSolValue() const = 0;
  2081 
  2082   };
  2083 
  2084 
  2085 
  2086 } //namespace lemon
  2087 
  2088 #endif //LEMON_LP_BASE_H