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