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