lemon/lp_base.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 19 Feb 2010 14:08:32 +0100
changeset 844 a6eb9698c321
parent 786 e20173729589
child 877 141f9c0db4a3
permissions -rw-r--r--
Support tolerance technique for BellmanFord (#51)

A new operation traits class BellmanFordToleranceOperationTraits
is introduced, which uses the tolerance technique in its less()
function. This class can be used with the SetOperationTraits
named template parameter.
     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 int _addRow(Value l, ExprIterator b, ExprIterator e, Value u) {
   947       int row = _addRow();
   948       _setRowCoeffs(row, b, e);
   949       _setRowLowerBound(row, l);
   950       _setRowUpperBound(row, u);
   951       return row;
   952     }
   953 
   954     virtual void _eraseCol(int col) = 0;
   955     virtual void _eraseRow(int row) = 0;
   956 
   957     virtual void _getColName(int col, std::string& name) const = 0;
   958     virtual void _setColName(int col, const std::string& name) = 0;
   959     virtual int _colByName(const std::string& name) const = 0;
   960 
   961     virtual void _getRowName(int row, std::string& name) const = 0;
   962     virtual void _setRowName(int row, const std::string& name) = 0;
   963     virtual int _rowByName(const std::string& name) const = 0;
   964 
   965     virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
   966     virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
   967 
   968     virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
   969     virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
   970 
   971     virtual void _setCoeff(int row, int col, Value value) = 0;
   972     virtual Value _getCoeff(int row, int col) const = 0;
   973 
   974     virtual void _setColLowerBound(int i, Value value) = 0;
   975     virtual Value _getColLowerBound(int i) const = 0;
   976 
   977     virtual void _setColUpperBound(int i, Value value) = 0;
   978     virtual Value _getColUpperBound(int i) const = 0;
   979 
   980     virtual void _setRowLowerBound(int i, Value value) = 0;
   981     virtual Value _getRowLowerBound(int i) const = 0;
   982 
   983     virtual void _setRowUpperBound(int i, Value value) = 0;
   984     virtual Value _getRowUpperBound(int i) const = 0;
   985 
   986     virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
   987     virtual void _getObjCoeffs(InsertIterator b) const = 0;
   988 
   989     virtual void _setObjCoeff(int i, Value obj_coef) = 0;
   990     virtual Value _getObjCoeff(int i) const = 0;
   991 
   992     virtual void _setSense(Sense) = 0;
   993     virtual Sense _getSense() const = 0;
   994 
   995     virtual void _clear() = 0;
   996 
   997     virtual const char* _solverName() const = 0;
   998 
   999     virtual void _messageLevel(MessageLevel level) = 0;
  1000 
  1001     //Own protected stuff
  1002 
  1003     //Constant component of the objective function
  1004     Value obj_const_comp;
  1005 
  1006     LpBase() : rows(), cols(), obj_const_comp(0) {}
  1007 
  1008   public:
  1009 
  1010     /// Virtual destructor
  1011     virtual ~LpBase() {}
  1012 
  1013     ///Gives back the name of the solver.
  1014     const char* solverName() const {return _solverName();}
  1015 
  1016     ///\name Build Up and Modify the LP
  1017 
  1018     ///@{
  1019 
  1020     ///Add a new empty column (i.e a new variable) to the LP
  1021     Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
  1022 
  1023     ///\brief Adds several new columns (i.e variables) at once
  1024     ///
  1025     ///This magic function takes a container as its argument and fills
  1026     ///its elements with new columns (i.e. variables)
  1027     ///\param t can be
  1028     ///- a standard STL compatible iterable container with
  1029     ///\ref Col as its \c values_type like
  1030     ///\code
  1031     ///std::vector<LpBase::Col>
  1032     ///std::list<LpBase::Col>
  1033     ///\endcode
  1034     ///- a standard STL compatible iterable container with
  1035     ///\ref Col as its \c mapped_type like
  1036     ///\code
  1037     ///std::map<AnyType,LpBase::Col>
  1038     ///\endcode
  1039     ///- an iterable lemon \ref concepts::WriteMap "write map" like
  1040     ///\code
  1041     ///ListGraph::NodeMap<LpBase::Col>
  1042     ///ListGraph::ArcMap<LpBase::Col>
  1043     ///\endcode
  1044     ///\return The number of the created column.
  1045 #ifdef DOXYGEN
  1046     template<class T>
  1047     int addColSet(T &t) { return 0;}
  1048 #else
  1049     template<class T>
  1050     typename enable_if<typename T::value_type::LpCol,int>::type
  1051     addColSet(T &t,dummy<0> = 0) {
  1052       int s=0;
  1053       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
  1054       return s;
  1055     }
  1056     template<class T>
  1057     typename enable_if<typename T::value_type::second_type::LpCol,
  1058                        int>::type
  1059     addColSet(T &t,dummy<1> = 1) {
  1060       int s=0;
  1061       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1062         i->second=addCol();
  1063         s++;
  1064       }
  1065       return s;
  1066     }
  1067     template<class T>
  1068     typename enable_if<typename T::MapIt::Value::LpCol,
  1069                        int>::type
  1070     addColSet(T &t,dummy<2> = 2) {
  1071       int s=0;
  1072       for(typename T::MapIt i(t); i!=INVALID; ++i)
  1073         {
  1074           i.set(addCol());
  1075           s++;
  1076         }
  1077       return s;
  1078     }
  1079 #endif
  1080 
  1081     ///Set a column (i.e a dual constraint) of the LP
  1082 
  1083     ///\param c is the column to be modified
  1084     ///\param e is a dual linear expression (see \ref DualExpr)
  1085     ///a better one.
  1086     void col(Col c, const DualExpr &e) {
  1087       e.simplify();
  1088       _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), rows),
  1089                     ExprIterator(e.comps.end(), rows));
  1090     }
  1091 
  1092     ///Get a column (i.e a dual constraint) of the LP
  1093 
  1094     ///\param c is the column to get
  1095     ///\return the dual expression associated to the column
  1096     DualExpr col(Col c) const {
  1097       DualExpr e;
  1098       _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
  1099       return e;
  1100     }
  1101 
  1102     ///Add a new column to the LP
  1103 
  1104     ///\param e is a dual linear expression (see \ref DualExpr)
  1105     ///\param o is the corresponding component of the objective
  1106     ///function. It is 0 by default.
  1107     ///\return The created column.
  1108     Col addCol(const DualExpr &e, Value o = 0) {
  1109       Col c=addCol();
  1110       col(c,e);
  1111       objCoeff(c,o);
  1112       return c;
  1113     }
  1114 
  1115     ///Add a new empty row (i.e a new constraint) to the LP
  1116 
  1117     ///This function adds a new empty row (i.e a new constraint) to the LP.
  1118     ///\return The created row
  1119     Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
  1120 
  1121     ///\brief Add several new rows (i.e constraints) at once
  1122     ///
  1123     ///This magic function takes a container as its argument and fills
  1124     ///its elements with new row (i.e. variables)
  1125     ///\param t can be
  1126     ///- a standard STL compatible iterable container with
  1127     ///\ref Row as its \c values_type like
  1128     ///\code
  1129     ///std::vector<LpBase::Row>
  1130     ///std::list<LpBase::Row>
  1131     ///\endcode
  1132     ///- a standard STL compatible iterable container with
  1133     ///\ref Row as its \c mapped_type like
  1134     ///\code
  1135     ///std::map<AnyType,LpBase::Row>
  1136     ///\endcode
  1137     ///- an iterable lemon \ref concepts::WriteMap "write map" like
  1138     ///\code
  1139     ///ListGraph::NodeMap<LpBase::Row>
  1140     ///ListGraph::ArcMap<LpBase::Row>
  1141     ///\endcode
  1142     ///\return The number of rows created.
  1143 #ifdef DOXYGEN
  1144     template<class T>
  1145     int addRowSet(T &t) { return 0;}
  1146 #else
  1147     template<class T>
  1148     typename enable_if<typename T::value_type::LpRow,int>::type
  1149     addRowSet(T &t, dummy<0> = 0) {
  1150       int s=0;
  1151       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
  1152       return s;
  1153     }
  1154     template<class T>
  1155     typename enable_if<typename T::value_type::second_type::LpRow, int>::type
  1156     addRowSet(T &t, dummy<1> = 1) {
  1157       int s=0;
  1158       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1159         i->second=addRow();
  1160         s++;
  1161       }
  1162       return s;
  1163     }
  1164     template<class T>
  1165     typename enable_if<typename T::MapIt::Value::LpRow, int>::type
  1166     addRowSet(T &t, dummy<2> = 2) {
  1167       int s=0;
  1168       for(typename T::MapIt i(t); i!=INVALID; ++i)
  1169         {
  1170           i.set(addRow());
  1171           s++;
  1172         }
  1173       return s;
  1174     }
  1175 #endif
  1176 
  1177     ///Set a row (i.e a constraint) of the LP
  1178 
  1179     ///\param r is the row to be modified
  1180     ///\param l is lower bound (-\ref INF means no bound)
  1181     ///\param e is a linear expression (see \ref Expr)
  1182     ///\param u is the upper bound (\ref INF means no bound)
  1183     void row(Row r, Value l, const Expr &e, Value u) {
  1184       e.simplify();
  1185       _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
  1186                     ExprIterator(e.comps.end(), cols));
  1187       _setRowLowerBound(rows(id(r)),l - *e);
  1188       _setRowUpperBound(rows(id(r)),u - *e);
  1189     }
  1190 
  1191     ///Set a row (i.e a constraint) of the LP
  1192 
  1193     ///\param r is the row to be modified
  1194     ///\param c is a linear expression (see \ref Constr)
  1195     void row(Row r, const Constr &c) {
  1196       row(r, c.lowerBounded()?c.lowerBound():-INF,
  1197           c.expr(), c.upperBounded()?c.upperBound():INF);
  1198     }
  1199 
  1200 
  1201     ///Get a row (i.e a constraint) of the LP
  1202 
  1203     ///\param r is the row to get
  1204     ///\return the expression associated to the row
  1205     Expr row(Row r) const {
  1206       Expr e;
  1207       _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
  1208       return e;
  1209     }
  1210 
  1211     ///Add a new row (i.e a new constraint) to the LP
  1212 
  1213     ///\param l is the lower bound (-\ref INF means no bound)
  1214     ///\param e is a linear expression (see \ref Expr)
  1215     ///\param u is the upper bound (\ref INF means no bound)
  1216     ///\return The created row.
  1217     Row addRow(Value l,const Expr &e, Value u) {
  1218       Row r;
  1219       e.simplify();
  1220       r._id = _addRowId(_addRow(l - *e, ExprIterator(e.comps.begin(), cols),
  1221                                 ExprIterator(e.comps.end(), cols), u - *e));
  1222       return r;
  1223     }
  1224 
  1225     ///Add a new row (i.e a new constraint) to the LP
  1226 
  1227     ///\param c is a linear expression (see \ref Constr)
  1228     ///\return The created row.
  1229     Row addRow(const Constr &c) {
  1230       Row r;
  1231       c.expr().simplify();
  1232       r._id = _addRowId(_addRow(c.lowerBounded()?c.lowerBound()-*c.expr():-INF, 
  1233                                 ExprIterator(c.expr().comps.begin(), cols),
  1234                                 ExprIterator(c.expr().comps.end(), cols),
  1235                                 c.upperBounded()?c.upperBound()-*c.expr():INF));
  1236       return r;
  1237     }
  1238     ///Erase a column (i.e a variable) from the LP
  1239 
  1240     ///\param c is the column to be deleted
  1241     void erase(Col c) {
  1242       _eraseCol(cols(id(c)));
  1243       _eraseColId(cols(id(c)));
  1244     }
  1245     ///Erase a row (i.e a constraint) from the LP
  1246 
  1247     ///\param r is the row to be deleted
  1248     void erase(Row r) {
  1249       _eraseRow(rows(id(r)));
  1250       _eraseRowId(rows(id(r)));
  1251     }
  1252 
  1253     /// Get the name of a column
  1254 
  1255     ///\param c is the coresponding column
  1256     ///\return The name of the colunm
  1257     std::string colName(Col c) const {
  1258       std::string name;
  1259       _getColName(cols(id(c)), name);
  1260       return name;
  1261     }
  1262 
  1263     /// Set the name of a column
  1264 
  1265     ///\param c is the coresponding column
  1266     ///\param name The name to be given
  1267     void colName(Col c, const std::string& name) {
  1268       _setColName(cols(id(c)), name);
  1269     }
  1270 
  1271     /// Get the column by its name
  1272 
  1273     ///\param name The name of the column
  1274     ///\return the proper column or \c INVALID
  1275     Col colByName(const std::string& name) const {
  1276       int k = _colByName(name);
  1277       return k != -1 ? Col(cols[k]) : Col(INVALID);
  1278     }
  1279 
  1280     /// Get the name of a row
  1281 
  1282     ///\param r is the coresponding row
  1283     ///\return The name of the row
  1284     std::string rowName(Row r) const {
  1285       std::string name;
  1286       _getRowName(rows(id(r)), name);
  1287       return name;
  1288     }
  1289 
  1290     /// Set the name of a row
  1291 
  1292     ///\param r is the coresponding row
  1293     ///\param name The name to be given
  1294     void rowName(Row r, const std::string& name) {
  1295       _setRowName(rows(id(r)), name);
  1296     }
  1297 
  1298     /// Get the row by its name
  1299 
  1300     ///\param name The name of the row
  1301     ///\return the proper row or \c INVALID
  1302     Row rowByName(const std::string& name) const {
  1303       int k = _rowByName(name);
  1304       return k != -1 ? Row(rows[k]) : Row(INVALID);
  1305     }
  1306 
  1307     /// Set an element of the coefficient matrix of the LP
  1308 
  1309     ///\param r is the row of the element to be modified
  1310     ///\param c is the column of the element to be modified
  1311     ///\param val is the new value of the coefficient
  1312     void coeff(Row r, Col c, Value val) {
  1313       _setCoeff(rows(id(r)),cols(id(c)), val);
  1314     }
  1315 
  1316     /// Get an element of the coefficient matrix of the LP
  1317 
  1318     ///\param r is the row of the element
  1319     ///\param c is the column of the element
  1320     ///\return the corresponding coefficient
  1321     Value coeff(Row r, Col c) const {
  1322       return _getCoeff(rows(id(r)),cols(id(c)));
  1323     }
  1324 
  1325     /// Set the lower bound of a column (i.e a variable)
  1326 
  1327     /// The lower bound of a variable (column) has to be given by an
  1328     /// extended number of type Value, i.e. a finite number of type
  1329     /// Value or -\ref INF.
  1330     void colLowerBound(Col c, Value value) {
  1331       _setColLowerBound(cols(id(c)),value);
  1332     }
  1333 
  1334     /// Get the lower bound of a column (i.e a variable)
  1335 
  1336     /// This function returns the lower bound for column (variable) \c c
  1337     /// (this might be -\ref INF as well).
  1338     ///\return The lower bound for column \c c
  1339     Value colLowerBound(Col c) const {
  1340       return _getColLowerBound(cols(id(c)));
  1341     }
  1342 
  1343     ///\brief Set the lower bound of  several columns
  1344     ///(i.e variables) at once
  1345     ///
  1346     ///This magic function takes a container as its argument
  1347     ///and applies the function on all of its elements.
  1348     ///The lower bound of a variable (column) has to be given by an
  1349     ///extended number of type Value, i.e. a finite number of type
  1350     ///Value or -\ref INF.
  1351 #ifdef DOXYGEN
  1352     template<class T>
  1353     void colLowerBound(T &t, Value value) { return 0;}
  1354 #else
  1355     template<class T>
  1356     typename enable_if<typename T::value_type::LpCol,void>::type
  1357     colLowerBound(T &t, Value value,dummy<0> = 0) {
  1358       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1359         colLowerBound(*i, value);
  1360       }
  1361     }
  1362     template<class T>
  1363     typename enable_if<typename T::value_type::second_type::LpCol,
  1364                        void>::type
  1365     colLowerBound(T &t, Value value,dummy<1> = 1) {
  1366       for(typename T::iterator i=t.begin();i!=t.end();++i) {
  1367         colLowerBound(i->second, value);
  1368       }
  1369     }
  1370     template<class T>
  1371     typename enable_if<typename T::MapIt::Value::LpCol,
  1372                        void>::type
  1373     colLowerBound(T &t, Value value,dummy<2> = 2) {
  1374       for(typename T::MapIt i(t); i!=INVALID; ++i){
  1375         colLowerBound(*i, value);
  1376       }
  1377     }
  1378 #endif
  1379 
  1380     /// Set the upper bound of a column (i.e a variable)
  1381 
  1382     /// The upper bound of a variable (column) has to be given by an
  1383     /// extended number of type Value, i.e. a finite number of type
  1384     /// Value or \ref INF.
  1385     void colUpperBound(Col c, Value value) {
  1386       _setColUpperBound(cols(id(c)),value);
  1387     };
  1388 
  1389     /// Get the upper bound of a column (i.e a variable)
  1390 
  1391     /// This function returns the upper bound for column (variable) \c c
  1392     /// (this might be \ref INF as well).
  1393     /// \return The upper bound for column \c c
  1394     Value colUpperBound(Col c) const {
  1395       return _getColUpperBound(cols(id(c)));
  1396     }
  1397 
  1398     ///\brief Set the upper bound of  several columns
  1399     ///(i.e variables) at once
  1400     ///
  1401     ///This magic function takes a container as its argument
  1402     ///and applies the function on all of its elements.
  1403     ///The upper bound of a variable (column) has to be given by an
  1404     ///extended number of type Value, i.e. a finite number of type
  1405     ///Value or \ref INF.
  1406 #ifdef DOXYGEN
  1407     template<class T>
  1408     void colUpperBound(T &t, Value value) { return 0;}
  1409 #else
  1410     template<class T1>
  1411     typename enable_if<typename T1::value_type::LpCol,void>::type
  1412     colUpperBound(T1 &t, Value value,dummy<0> = 0) {
  1413       for(typename T1::iterator i=t.begin();i!=t.end();++i) {
  1414         colUpperBound(*i, value);
  1415       }
  1416     }
  1417     template<class T1>
  1418     typename enable_if<typename T1::value_type::second_type::LpCol,
  1419                        void>::type
  1420     colUpperBound(T1 &t, Value value,dummy<1> = 1) {
  1421       for(typename T1::iterator i=t.begin();i!=t.end();++i) {
  1422         colUpperBound(i->second, value);
  1423       }
  1424     }
  1425     template<class T1>
  1426     typename enable_if<typename T1::MapIt::Value::LpCol,
  1427                        void>::type
  1428     colUpperBound(T1 &t, Value value,dummy<2> = 2) {
  1429       for(typename T1::MapIt i(t); i!=INVALID; ++i){
  1430         colUpperBound(*i, value);
  1431       }
  1432     }
  1433 #endif
  1434 
  1435     /// Set the lower and the upper bounds of a column (i.e a variable)
  1436 
  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     void colBounds(Col c, Value lower, Value upper) {
  1442       _setColLowerBound(cols(id(c)),lower);
  1443       _setColUpperBound(cols(id(c)),upper);
  1444     }
  1445 
  1446     ///\brief Set the lower and the upper bound of several columns
  1447     ///(i.e variables) at once
  1448     ///
  1449     ///This magic function takes a container as its argument
  1450     ///and applies the function on all of its elements.
  1451     /// The lower and the upper bounds of
  1452     /// a variable (column) have to be given by an
  1453     /// extended number of type Value, i.e. a finite number of type
  1454     /// Value, -\ref INF or \ref INF.
  1455 #ifdef DOXYGEN
  1456     template<class T>
  1457     void colBounds(T &t, Value lower, Value upper) { return 0;}
  1458 #else
  1459     template<class T2>
  1460     typename enable_if<typename T2::value_type::LpCol,void>::type
  1461     colBounds(T2 &t, Value lower, Value upper,dummy<0> = 0) {
  1462       for(typename T2::iterator i=t.begin();i!=t.end();++i) {
  1463         colBounds(*i, lower, upper);
  1464       }
  1465     }
  1466     template<class T2>
  1467     typename enable_if<typename T2::value_type::second_type::LpCol, void>::type
  1468     colBounds(T2 &t, Value lower, Value upper,dummy<1> = 1) {
  1469       for(typename T2::iterator i=t.begin();i!=t.end();++i) {
  1470         colBounds(i->second, lower, upper);
  1471       }
  1472     }
  1473     template<class T2>
  1474     typename enable_if<typename T2::MapIt::Value::LpCol, void>::type
  1475     colBounds(T2 &t, Value lower, Value upper,dummy<2> = 2) {
  1476       for(typename T2::MapIt i(t); i!=INVALID; ++i){
  1477         colBounds(*i, lower, upper);
  1478       }
  1479     }
  1480 #endif
  1481 
  1482     /// Set the lower bound of a row (i.e a constraint)
  1483 
  1484     /// The lower bound of a constraint (row) has to be given by an
  1485     /// extended number of type Value, i.e. a finite number of type
  1486     /// Value or -\ref INF.
  1487     void rowLowerBound(Row r, Value value) {
  1488       _setRowLowerBound(rows(id(r)),value);
  1489     }
  1490 
  1491     /// Get the lower bound of a row (i.e a constraint)
  1492 
  1493     /// This function returns the lower bound for row (constraint) \c c
  1494     /// (this might be -\ref INF as well).
  1495     ///\return The lower bound for row \c r
  1496     Value rowLowerBound(Row r) const {
  1497       return _getRowLowerBound(rows(id(r)));
  1498     }
  1499 
  1500     /// Set the upper bound of a row (i.e a constraint)
  1501 
  1502     /// The upper bound of a constraint (row) has to be given by an
  1503     /// extended number of type Value, i.e. a finite number of type
  1504     /// Value or -\ref INF.
  1505     void rowUpperBound(Row r, Value value) {
  1506       _setRowUpperBound(rows(id(r)),value);
  1507     }
  1508 
  1509     /// Get the upper bound of a row (i.e a constraint)
  1510 
  1511     /// This function returns the upper bound for row (constraint) \c c
  1512     /// (this might be -\ref INF as well).
  1513     ///\return The upper bound for row \c r
  1514     Value rowUpperBound(Row r) const {
  1515       return _getRowUpperBound(rows(id(r)));
  1516     }
  1517 
  1518     ///Set an element of the objective function
  1519     void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
  1520 
  1521     ///Get an element of the objective function
  1522     Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
  1523 
  1524     ///Set the objective function
  1525 
  1526     ///\param e is a linear expression of type \ref Expr.
  1527     ///
  1528     void obj(const Expr& e) {
  1529       _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
  1530                     ExprIterator(e.comps.end(), cols));
  1531       obj_const_comp = *e;
  1532     }
  1533 
  1534     ///Get the objective function
  1535 
  1536     ///\return the objective function as a linear expression of type
  1537     ///Expr.
  1538     Expr obj() const {
  1539       Expr e;
  1540       _getObjCoeffs(InsertIterator(e.comps, cols));
  1541       *e = obj_const_comp;
  1542       return e;
  1543     }
  1544 
  1545 
  1546     ///Set the direction of optimization
  1547     void sense(Sense sense) { _setSense(sense); }
  1548 
  1549     ///Query the direction of the optimization
  1550     Sense sense() const {return _getSense(); }
  1551 
  1552     ///Set the sense to maximization
  1553     void max() { _setSense(MAX); }
  1554 
  1555     ///Set the sense to maximization
  1556     void min() { _setSense(MIN); }
  1557 
  1558     ///Clears the problem
  1559     void clear() { _clear(); }
  1560 
  1561     /// Sets the message level of the solver
  1562     void messageLevel(MessageLevel level) { _messageLevel(level); }
  1563 
  1564     ///@}
  1565 
  1566   };
  1567 
  1568   /// Addition
  1569 
  1570   ///\relates LpBase::Expr
  1571   ///
  1572   inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
  1573     LpBase::Expr tmp(a);
  1574     tmp+=b;
  1575     return tmp;
  1576   }
  1577   ///Substraction
  1578 
  1579   ///\relates LpBase::Expr
  1580   ///
  1581   inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
  1582     LpBase::Expr tmp(a);
  1583     tmp-=b;
  1584     return tmp;
  1585   }
  1586   ///Multiply with constant
  1587 
  1588   ///\relates LpBase::Expr
  1589   ///
  1590   inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
  1591     LpBase::Expr tmp(a);
  1592     tmp*=b;
  1593     return tmp;
  1594   }
  1595 
  1596   ///Multiply with constant
  1597 
  1598   ///\relates LpBase::Expr
  1599   ///
  1600   inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
  1601     LpBase::Expr tmp(b);
  1602     tmp*=a;
  1603     return tmp;
  1604   }
  1605   ///Divide with constant
  1606 
  1607   ///\relates LpBase::Expr
  1608   ///
  1609   inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
  1610     LpBase::Expr tmp(a);
  1611     tmp/=b;
  1612     return tmp;
  1613   }
  1614 
  1615   ///Create constraint
  1616 
  1617   ///\relates LpBase::Constr
  1618   ///
  1619   inline LpBase::Constr operator<=(const LpBase::Expr &e,
  1620                                    const LpBase::Expr &f) {
  1621     return LpBase::Constr(0, f - e, LpBase::INF);
  1622   }
  1623 
  1624   ///Create constraint
  1625 
  1626   ///\relates LpBase::Constr
  1627   ///
  1628   inline LpBase::Constr operator<=(const LpBase::Value &e,
  1629                                    const LpBase::Expr &f) {
  1630     return LpBase::Constr(e, f, LpBase::NaN);
  1631   }
  1632 
  1633   ///Create constraint
  1634 
  1635   ///\relates LpBase::Constr
  1636   ///
  1637   inline LpBase::Constr operator<=(const LpBase::Expr &e,
  1638                                    const LpBase::Value &f) {
  1639     return LpBase::Constr(- LpBase::INF, e, f);
  1640   }
  1641 
  1642   ///Create constraint
  1643 
  1644   ///\relates LpBase::Constr
  1645   ///
  1646   inline LpBase::Constr operator>=(const LpBase::Expr &e,
  1647                                    const LpBase::Expr &f) {
  1648     return LpBase::Constr(0, e - f, LpBase::INF);
  1649   }
  1650 
  1651 
  1652   ///Create constraint
  1653 
  1654   ///\relates LpBase::Constr
  1655   ///
  1656   inline LpBase::Constr operator>=(const LpBase::Value &e,
  1657                                    const LpBase::Expr &f) {
  1658     return LpBase::Constr(LpBase::NaN, f, e);
  1659   }
  1660 
  1661 
  1662   ///Create constraint
  1663 
  1664   ///\relates LpBase::Constr
  1665   ///
  1666   inline LpBase::Constr operator>=(const LpBase::Expr &e,
  1667                                    const LpBase::Value &f) {
  1668     return LpBase::Constr(f, e, LpBase::INF);
  1669   }
  1670 
  1671   ///Create constraint
  1672 
  1673   ///\relates LpBase::Constr
  1674   ///
  1675   inline LpBase::Constr operator==(const LpBase::Expr &e,
  1676                                    const LpBase::Value &f) {
  1677     return LpBase::Constr(f, e, f);
  1678   }
  1679 
  1680   ///Create constraint
  1681 
  1682   ///\relates LpBase::Constr
  1683   ///
  1684   inline LpBase::Constr operator==(const LpBase::Expr &e,
  1685                                    const LpBase::Expr &f) {
  1686     return LpBase::Constr(0, f - e, 0);
  1687   }
  1688 
  1689   ///Create constraint
  1690 
  1691   ///\relates LpBase::Constr
  1692   ///
  1693   inline LpBase::Constr operator<=(const LpBase::Value &n,
  1694                                    const LpBase::Constr &c) {
  1695     LpBase::Constr tmp(c);
  1696     LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
  1697     tmp.lowerBound()=n;
  1698     return tmp;
  1699   }
  1700   ///Create constraint
  1701 
  1702   ///\relates LpBase::Constr
  1703   ///
  1704   inline LpBase::Constr operator<=(const LpBase::Constr &c,
  1705                                    const LpBase::Value &n)
  1706   {
  1707     LpBase::Constr tmp(c);
  1708     LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
  1709     tmp.upperBound()=n;
  1710     return tmp;
  1711   }
  1712 
  1713   ///Create constraint
  1714 
  1715   ///\relates LpBase::Constr
  1716   ///
  1717   inline LpBase::Constr operator>=(const LpBase::Value &n,
  1718                                    const LpBase::Constr &c) {
  1719     LpBase::Constr tmp(c);
  1720     LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
  1721     tmp.upperBound()=n;
  1722     return tmp;
  1723   }
  1724   ///Create constraint
  1725 
  1726   ///\relates LpBase::Constr
  1727   ///
  1728   inline LpBase::Constr operator>=(const LpBase::Constr &c,
  1729                                    const LpBase::Value &n)
  1730   {
  1731     LpBase::Constr tmp(c);
  1732     LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
  1733     tmp.lowerBound()=n;
  1734     return tmp;
  1735   }
  1736 
  1737   ///Addition
  1738 
  1739   ///\relates LpBase::DualExpr
  1740   ///
  1741   inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
  1742                                     const LpBase::DualExpr &b) {
  1743     LpBase::DualExpr tmp(a);
  1744     tmp+=b;
  1745     return tmp;
  1746   }
  1747   ///Substraction
  1748 
  1749   ///\relates LpBase::DualExpr
  1750   ///
  1751   inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
  1752                                     const LpBase::DualExpr &b) {
  1753     LpBase::DualExpr tmp(a);
  1754     tmp-=b;
  1755     return tmp;
  1756   }
  1757   ///Multiply with constant
  1758 
  1759   ///\relates LpBase::DualExpr
  1760   ///
  1761   inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
  1762                                     const LpBase::Value &b) {
  1763     LpBase::DualExpr tmp(a);
  1764     tmp*=b;
  1765     return tmp;
  1766   }
  1767 
  1768   ///Multiply with constant
  1769 
  1770   ///\relates LpBase::DualExpr
  1771   ///
  1772   inline LpBase::DualExpr operator*(const LpBase::Value &a,
  1773                                     const LpBase::DualExpr &b) {
  1774     LpBase::DualExpr tmp(b);
  1775     tmp*=a;
  1776     return tmp;
  1777   }
  1778   ///Divide with constant
  1779 
  1780   ///\relates LpBase::DualExpr
  1781   ///
  1782   inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
  1783                                     const LpBase::Value &b) {
  1784     LpBase::DualExpr tmp(a);
  1785     tmp/=b;
  1786     return tmp;
  1787   }
  1788 
  1789   /// \ingroup lp_group
  1790   ///
  1791   /// \brief Common base class for LP solvers
  1792   ///
  1793   /// This class is an abstract base class for LP solvers. This class
  1794   /// provides a full interface for set and modify an LP problem,
  1795   /// solve it and retrieve the solution. You can use one of the
  1796   /// descendants as a concrete implementation, or the \c Lp
  1797   /// default LP solver. However, if you would like to handle LP
  1798   /// solvers as reference or pointer in a generic way, you can use
  1799   /// this class directly.
  1800   class LpSolver : virtual public LpBase {
  1801   public:
  1802 
  1803     /// The problem types for primal and dual problems
  1804     enum ProblemType {
  1805       /// = 0. Feasible solution hasn't been found (but may exist).
  1806       UNDEFINED = 0,
  1807       /// = 1. The problem has no feasible solution.
  1808       INFEASIBLE = 1,
  1809       /// = 2. Feasible solution found.
  1810       FEASIBLE = 2,
  1811       /// = 3. Optimal solution exists and found.
  1812       OPTIMAL = 3,
  1813       /// = 4. The cost function is unbounded.
  1814       UNBOUNDED = 4
  1815     };
  1816 
  1817     ///The basis status of variables
  1818     enum VarStatus {
  1819       /// The variable is in the basis
  1820       BASIC, 
  1821       /// The variable is free, but not basic
  1822       FREE,
  1823       /// The variable has active lower bound 
  1824       LOWER,
  1825       /// The variable has active upper bound
  1826       UPPER,
  1827       /// The variable is non-basic and fixed
  1828       FIXED
  1829     };
  1830 
  1831   protected:
  1832 
  1833     virtual SolveExitStatus _solve() = 0;
  1834 
  1835     virtual Value _getPrimal(int i) const = 0;
  1836     virtual Value _getDual(int i) const = 0;
  1837 
  1838     virtual Value _getPrimalRay(int i) const = 0;
  1839     virtual Value _getDualRay(int i) const = 0;
  1840 
  1841     virtual Value _getPrimalValue() const = 0;
  1842 
  1843     virtual VarStatus _getColStatus(int i) const = 0;
  1844     virtual VarStatus _getRowStatus(int i) const = 0;
  1845 
  1846     virtual ProblemType _getPrimalType() const = 0;
  1847     virtual ProblemType _getDualType() const = 0;
  1848 
  1849   public:
  1850 
  1851     ///Allocate a new LP problem instance
  1852     virtual LpSolver* newSolver() const = 0;
  1853     ///Make a copy of the LP problem
  1854     virtual LpSolver* cloneSolver() const = 0;
  1855 
  1856     ///\name Solve the LP
  1857 
  1858     ///@{
  1859 
  1860     ///\e Solve the LP problem at hand
  1861     ///
  1862     ///\return The result of the optimization procedure. Possible
  1863     ///values and their meanings can be found in the documentation of
  1864     ///\ref SolveExitStatus.
  1865     SolveExitStatus solve() { return _solve(); }
  1866 
  1867     ///@}
  1868 
  1869     ///\name Obtain the Solution
  1870 
  1871     ///@{
  1872 
  1873     /// The type of the primal problem
  1874     ProblemType primalType() const {
  1875       return _getPrimalType();
  1876     }
  1877 
  1878     /// The type of the dual problem
  1879     ProblemType dualType() const {
  1880       return _getDualType();
  1881     }
  1882 
  1883     /// Return the primal value of the column
  1884 
  1885     /// Return the primal value of the column.
  1886     /// \pre The problem is solved.
  1887     Value primal(Col c) const { return _getPrimal(cols(id(c))); }
  1888 
  1889     /// Return the primal value of the expression
  1890 
  1891     /// Return the primal value of the expression, i.e. the dot
  1892     /// product of the primal solution and the expression.
  1893     /// \pre The problem is solved.
  1894     Value primal(const Expr& e) const {
  1895       double res = *e;
  1896       for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
  1897         res += *c * primal(c);
  1898       }
  1899       return res;
  1900     }
  1901     /// Returns a component of the primal ray
  1902     
  1903     /// The primal ray is solution of the modified primal problem,
  1904     /// where we change each finite bound to 0, and we looking for a
  1905     /// negative objective value in case of minimization, and positive
  1906     /// objective value for maximization. If there is such solution,
  1907     /// that proofs the unsolvability of the dual problem, and if a
  1908     /// feasible primal solution exists, then the unboundness of
  1909     /// primal problem.
  1910     ///
  1911     /// \pre The problem is solved and the dual problem is infeasible.
  1912     /// \note Some solvers does not provide primal ray calculation
  1913     /// functions.
  1914     Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
  1915 
  1916     /// Return the dual value of the row
  1917 
  1918     /// Return the dual value of the row.
  1919     /// \pre The problem is solved.
  1920     Value dual(Row r) const { return _getDual(rows(id(r))); }
  1921 
  1922     /// Return the dual value of the dual expression
  1923 
  1924     /// Return the dual value of the dual expression, i.e. the dot
  1925     /// product of the dual solution and the dual expression.
  1926     /// \pre The problem is solved.
  1927     Value dual(const DualExpr& e) const {
  1928       double res = 0.0;
  1929       for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
  1930         res += *r * dual(r);
  1931       }
  1932       return res;
  1933     }
  1934 
  1935     /// Returns a component of the dual ray
  1936     
  1937     /// The dual ray is solution of the modified primal problem, where
  1938     /// we change each finite bound to 0 (i.e. the objective function
  1939     /// coefficients in the primal problem), and we looking for a
  1940     /// ositive objective value. If there is such solution, that
  1941     /// proofs the unsolvability of the primal problem, and if a
  1942     /// feasible dual solution exists, then the unboundness of
  1943     /// dual problem.
  1944     ///
  1945     /// \pre The problem is solved and the primal problem is infeasible.
  1946     /// \note Some solvers does not provide dual ray calculation
  1947     /// functions.
  1948     Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
  1949 
  1950     /// Return the basis status of the column
  1951 
  1952     /// \see VarStatus
  1953     VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
  1954 
  1955     /// Return the basis status of the row
  1956 
  1957     /// \see VarStatus
  1958     VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
  1959 
  1960     ///The value of the objective function
  1961 
  1962     ///\return
  1963     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
  1964     /// of the primal problem, depending on whether we minimize or maximize.
  1965     ///- \ref NaN if no primal solution is found.
  1966     ///- The (finite) objective value if an optimal solution is found.
  1967     Value primal() const { return _getPrimalValue()+obj_const_comp;}
  1968     ///@}
  1969 
  1970   protected:
  1971 
  1972   };
  1973 
  1974 
  1975   /// \ingroup lp_group
  1976   ///
  1977   /// \brief Common base class for MIP solvers
  1978   ///
  1979   /// This class is an abstract base class for MIP solvers. This class
  1980   /// provides a full interface for set and modify an MIP problem,
  1981   /// solve it and retrieve the solution. You can use one of the
  1982   /// descendants as a concrete implementation, or the \c Lp
  1983   /// default MIP solver. However, if you would like to handle MIP
  1984   /// solvers as reference or pointer in a generic way, you can use
  1985   /// this class directly.
  1986   class MipSolver : virtual public LpBase {
  1987   public:
  1988 
  1989     /// The problem types for MIP problems
  1990     enum ProblemType {
  1991       /// = 0. Feasible solution hasn't been found (but may exist).
  1992       UNDEFINED = 0,
  1993       /// = 1. The problem has no feasible solution.
  1994       INFEASIBLE = 1,
  1995       /// = 2. Feasible solution found.
  1996       FEASIBLE = 2,
  1997       /// = 3. Optimal solution exists and found.
  1998       OPTIMAL = 3,
  1999       /// = 4. The cost function is unbounded.
  2000       ///The Mip or at least the relaxed problem is unbounded.
  2001       UNBOUNDED = 4
  2002     };
  2003 
  2004     ///Allocate a new MIP problem instance
  2005     virtual MipSolver* newSolver() const = 0;
  2006     ///Make a copy of the MIP problem
  2007     virtual MipSolver* cloneSolver() const = 0;
  2008 
  2009     ///\name Solve the MIP
  2010 
  2011     ///@{
  2012 
  2013     /// Solve the MIP problem at hand
  2014     ///
  2015     ///\return The result of the optimization procedure. Possible
  2016     ///values and their meanings can be found in the documentation of
  2017     ///\ref SolveExitStatus.
  2018     SolveExitStatus solve() { return _solve(); }
  2019 
  2020     ///@}
  2021 
  2022     ///\name Set Column Type
  2023     ///@{
  2024 
  2025     ///Possible variable (column) types (e.g. real, integer, binary etc.)
  2026     enum ColTypes {
  2027       /// = 0. Continuous variable (default).
  2028       REAL = 0,
  2029       /// = 1. Integer variable.
  2030       INTEGER = 1
  2031     };
  2032 
  2033     ///Sets the type of the given column to the given type
  2034 
  2035     ///Sets the type of the given column to the given type.
  2036     ///
  2037     void colType(Col c, ColTypes col_type) {
  2038       _setColType(cols(id(c)),col_type);
  2039     }
  2040 
  2041     ///Gives back the type of the column.
  2042 
  2043     ///Gives back the type of the column.
  2044     ///
  2045     ColTypes colType(Col c) const {
  2046       return _getColType(cols(id(c)));
  2047     }
  2048     ///@}
  2049 
  2050     ///\name Obtain the Solution
  2051 
  2052     ///@{
  2053 
  2054     /// The type of the MIP problem
  2055     ProblemType type() const {
  2056       return _getType();
  2057     }
  2058 
  2059     /// Return the value of the row in the solution
  2060 
  2061     ///  Return the value of the row in the solution.
  2062     /// \pre The problem is solved.
  2063     Value sol(Col c) const { return _getSol(cols(id(c))); }
  2064 
  2065     /// Return the value of the expression in the solution
  2066 
  2067     /// Return the value of the expression in the solution, i.e. the
  2068     /// dot product of the solution and the expression.
  2069     /// \pre The problem is solved.
  2070     Value sol(const Expr& e) const {
  2071       double res = *e;
  2072       for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
  2073         res += *c * sol(c);
  2074       }
  2075       return res;
  2076     }
  2077     ///The value of the objective function
  2078     
  2079     ///\return
  2080     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
  2081     /// of the problem, depending on whether we minimize or maximize.
  2082     ///- \ref NaN if no primal solution is found.
  2083     ///- The (finite) objective value if an optimal solution is found.
  2084     Value solValue() const { return _getSolValue()+obj_const_comp;}
  2085     ///@}
  2086 
  2087   protected:
  2088 
  2089     virtual SolveExitStatus _solve() = 0;
  2090     virtual ColTypes _getColType(int col) const = 0;
  2091     virtual void _setColType(int col, ColTypes col_type) = 0;
  2092     virtual ProblemType _getType() const = 0;
  2093     virtual Value _getSol(int i) const = 0;
  2094     virtual Value _getSolValue() const = 0;
  2095 
  2096   };
  2097 
  2098 
  2099 
  2100 } //namespace lemon
  2101 
  2102 #endif //LEMON_LP_BASE_H