deba@458: /* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@458:  *
deba@458:  * This file is a part of LEMON, a generic C++ optimization library.
deba@458:  *
deba@458:  * Copyright (C) 2003-2008
deba@458:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@458:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@458:  *
deba@458:  * Permission to use, modify and distribute this software is granted
deba@458:  * provided that this copyright notice appears in all copies. For
deba@458:  * precise terms see the accompanying LICENSE file.
deba@458:  *
deba@458:  * This software is provided "AS IS" with no warranty of any kind,
deba@458:  * express or implied, and with no claim as to its suitability for any
deba@458:  * purpose.
deba@458:  *
deba@458:  */
deba@458: 
deba@458: #ifndef LEMON_LP_BASE_H
deba@458: #define LEMON_LP_BASE_H
deba@458: 
deba@458: #include<iostream>
deba@458: #include<vector>
deba@458: #include<map>
deba@458: #include<limits>
deba@458: #include<lemon/math.h>
deba@458: 
deba@459: #include<lemon/error.h>
deba@459: #include<lemon/assert.h>
deba@459: 
deba@458: #include<lemon/core.h>
deba@459: #include<lemon/bits/solver_bits.h>
deba@458: 
deba@458: ///\file
deba@458: ///\brief The interface of the LP solver interface.
deba@458: ///\ingroup lp_group
deba@458: namespace lemon {
deba@458: 
deba@459:   ///Common base class for LP and MIP solvers
deba@458: 
deba@459:   ///Usually this class is not used directly, please use one of the concrete
deba@459:   ///implementations of the solver interface.
deba@458:   ///\ingroup lp_group
deba@459:   class LpBase {
deba@458: 
deba@458:   protected:
deba@458: 
deba@459:     _solver_bits::VarIndex rows;
deba@459:     _solver_bits::VarIndex cols;
deba@458: 
deba@458:   public:
deba@458: 
deba@458:     ///Possible outcomes of an LP solving procedure
deba@458:     enum SolveExitStatus {
kpeter@576:       /// = 0. It means that the problem has been successfully solved: either
deba@458:       ///an optimal solution has been found or infeasibility/unboundedness
deba@458:       ///has been proved.
deba@458:       SOLVED = 0,
kpeter@576:       /// = 1. Any other case (including the case when some user specified
kpeter@576:       ///limit has been exceeded).
deba@458:       UNSOLVED = 1
deba@458:     };
deba@458: 
deba@459:     ///Direction of the optimization
deba@459:     enum Sense {
deba@459:       /// Minimization
deba@459:       MIN,
deba@459:       /// Maximization
deba@459:       MAX
deba@458:     };
deba@458: 
deba@568:     ///Enum for \c messageLevel() parameter
deba@568:     enum MessageLevel {
kpeter@576:       /// No output (default value).
deba@568:       MESSAGE_NOTHING,
kpeter@576:       /// Error messages only.
deba@568:       MESSAGE_ERROR,
kpeter@576:       /// Warnings.
deba@568:       MESSAGE_WARNING,
kpeter@576:       /// Normal output.
deba@568:       MESSAGE_NORMAL,
kpeter@576:       /// Verbose output.
deba@568:       MESSAGE_VERBOSE
deba@568:     };
deba@568:     
deba@568: 
deba@458:     ///The floating point type used by the solver
deba@458:     typedef double Value;
deba@458:     ///The infinity constant
deba@458:     static const Value INF;
deba@458:     ///The not a number constant
deba@458:     static const Value NaN;
deba@458: 
deba@458:     friend class Col;
deba@458:     friend class ColIt;
deba@458:     friend class Row;
deba@459:     friend class RowIt;
deba@458: 
deba@458:     ///Refer to a column of the LP.
deba@458: 
deba@458:     ///This type is used to refer to a column of the LP.
deba@458:     ///
deba@458:     ///Its value remains valid and correct even after the addition or erase of
deba@458:     ///other columns.
deba@458:     ///
deba@459:     ///\note This class is similar to other Item types in LEMON, like
deba@459:     ///Node and Arc types in digraph.
deba@458:     class Col {
deba@459:       friend class LpBase;
deba@458:     protected:
deba@459:       int _id;
deba@459:       explicit Col(int id) : _id(id) {}
deba@458:     public:
deba@458:       typedef Value ExprValue;
deba@459:       typedef True LpCol;
deba@459:       /// Default constructor
deba@459:       
deba@459:       /// \warning The default constructor sets the Col to an
deba@459:       /// undefined value.
deba@458:       Col() {}
deba@459:       /// Invalid constructor \& conversion.
deba@459:       
deba@459:       /// This constructor initializes the Col to be invalid.
deba@459:       /// \sa Invalid for more details.      
deba@459:       Col(const Invalid&) : _id(-1) {}
deba@459:       /// Equality operator
deba@459: 
deba@459:       /// Two \ref Col "Col"s are equal if and only if they point to
deba@459:       /// the same LP column or both are invalid.
deba@459:       bool operator==(Col c) const  {return _id == c._id;}
deba@459:       /// Inequality operator
deba@459: 
deba@459:       /// \sa operator==(Col c)
deba@459:       ///
deba@459:       bool operator!=(Col c) const  {return _id != c._id;}
deba@459:       /// Artificial ordering operator.
deba@459: 
deba@459:       /// To allow the use of this object in std::map or similar
deba@459:       /// associative container we require this.
deba@459:       ///
deba@459:       /// \note This operator only have to define some strict ordering of
deba@459:       /// the items; this order has nothing to do with the iteration
deba@459:       /// ordering of the items.
deba@459:       bool operator<(Col c) const  {return _id < c._id;}
deba@458:     };
deba@458: 
deba@459:     ///Iterator for iterate over the columns of an LP problem
deba@459: 
deba@459:     /// Its usage is quite simple, for example you can count the number
deba@459:     /// of columns in an LP \c lp:
deba@459:     ///\code
deba@459:     /// int count=0;
deba@459:     /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
deba@459:     ///\endcode
deba@458:     class ColIt : public Col {
deba@459:       const LpBase *_solver;
deba@458:     public:
deba@459:       /// Default constructor
deba@459:       
deba@459:       /// \warning The default constructor sets the iterator
deba@459:       /// to an undefined value.
deba@458:       ColIt() {}
deba@459:       /// Sets the iterator to the first Col
deba@459:       
deba@459:       /// Sets the iterator to the first Col.
deba@459:       ///
deba@459:       ColIt(const LpBase &solver) : _solver(&solver)
deba@458:       {
deba@459:         _solver->cols.firstItem(_id);
deba@458:       }
deba@459:       /// Invalid constructor \& conversion
deba@459:       
deba@459:       /// Initialize the iterator to be invalid.
deba@459:       /// \sa Invalid for more details.
deba@458:       ColIt(const Invalid&) : Col(INVALID) {}
deba@459:       /// Next column
deba@459:       
deba@459:       /// Assign the iterator to the next column.
deba@459:       ///
deba@458:       ColIt &operator++()
deba@458:       {
deba@459:         _solver->cols.nextItem(_id);
deba@458:         return *this;
deba@458:       }
deba@458:     };
deba@458: 
deba@459:     /// \brief Returns the ID of the column.
deba@459:     static int id(const Col& col) { return col._id; }
deba@459:     /// \brief Returns the column with the given ID.
deba@459:     ///
deba@459:     /// \pre The argument should be a valid column ID in the LP problem.
deba@459:     static Col colFromId(int id) { return Col(id); }
deba@458: 
deba@458:     ///Refer to a row of the LP.
deba@458: 
deba@458:     ///This type is used to refer to a row of the LP.
deba@458:     ///
deba@458:     ///Its value remains valid and correct even after the addition or erase of
deba@458:     ///other rows.
deba@458:     ///
deba@459:     ///\note This class is similar to other Item types in LEMON, like
deba@459:     ///Node and Arc types in digraph.
deba@458:     class Row {
deba@459:       friend class LpBase;
deba@458:     protected:
deba@459:       int _id;
deba@459:       explicit Row(int id) : _id(id) {}
deba@458:     public:
deba@458:       typedef Value ExprValue;
deba@459:       typedef True LpRow;
deba@459:       /// Default constructor
deba@459:       
deba@459:       /// \warning The default constructor sets the Row to an
deba@459:       /// undefined value.
deba@458:       Row() {}
deba@459:       /// Invalid constructor \& conversion.
deba@459:       
deba@459:       /// This constructor initializes the Row to be invalid.
deba@459:       /// \sa Invalid for more details.      
deba@459:       Row(const Invalid&) : _id(-1) {}
deba@459:       /// Equality operator
deba@458: 
deba@459:       /// Two \ref Row "Row"s are equal if and only if they point to
deba@459:       /// the same LP row or both are invalid.
deba@459:       bool operator==(Row r) const  {return _id == r._id;}
deba@459:       /// Inequality operator
deba@459:       
deba@459:       /// \sa operator==(Row r)
deba@459:       ///
deba@459:       bool operator!=(Row r) const  {return _id != r._id;}
deba@459:       /// Artificial ordering operator.
deba@459: 
deba@459:       /// To allow the use of this object in std::map or similar
deba@459:       /// associative container we require this.
deba@459:       ///
deba@459:       /// \note This operator only have to define some strict ordering of
deba@459:       /// the items; this order has nothing to do with the iteration
deba@459:       /// ordering of the items.
deba@459:       bool operator<(Row r) const  {return _id < r._id;}
deba@458:     };
deba@458: 
deba@459:     ///Iterator for iterate over the rows of an LP problem
deba@459: 
deba@459:     /// Its usage is quite simple, for example you can count the number
deba@459:     /// of rows in an LP \c lp:
deba@459:     ///\code
deba@459:     /// int count=0;
deba@459:     /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
deba@459:     ///\endcode
deba@458:     class RowIt : public Row {
deba@459:       const LpBase *_solver;
deba@458:     public:
deba@459:       /// Default constructor
deba@459:       
deba@459:       /// \warning The default constructor sets the iterator
deba@459:       /// to an undefined value.
deba@458:       RowIt() {}
deba@459:       /// Sets the iterator to the first Row
deba@459:       
deba@459:       /// Sets the iterator to the first Row.
deba@459:       ///
deba@459:       RowIt(const LpBase &solver) : _solver(&solver)
deba@458:       {
deba@459:         _solver->rows.firstItem(_id);
deba@458:       }
deba@459:       /// Invalid constructor \& conversion
deba@459:       
deba@459:       /// Initialize the iterator to be invalid.
deba@459:       /// \sa Invalid for more details.
deba@458:       RowIt(const Invalid&) : Row(INVALID) {}
deba@459:       /// Next row
deba@459:       
deba@459:       /// Assign the iterator to the next row.
deba@459:       ///
deba@458:       RowIt &operator++()
deba@458:       {
deba@459:         _solver->rows.nextItem(_id);
deba@458:         return *this;
deba@458:       }
deba@458:     };
deba@458: 
deba@459:     /// \brief Returns the ID of the row.
deba@459:     static int id(const Row& row) { return row._id; }
deba@459:     /// \brief Returns the row with the given ID.
deba@459:     ///
deba@459:     /// \pre The argument should be a valid row ID in the LP problem.
deba@459:     static Row rowFromId(int id) { return Row(id); }
deba@458: 
deba@458:   public:
deba@458: 
deba@458:     ///Linear expression of variables and a constant component
deba@458: 
deba@458:     ///This data structure stores a linear expression of the variables
deba@458:     ///(\ref Col "Col"s) and also has a constant component.
deba@458:     ///
deba@458:     ///There are several ways to access and modify the contents of this
deba@458:     ///container.
deba@458:     ///\code
deba@458:     ///e[v]=5;
deba@458:     ///e[v]+=12;
deba@458:     ///e.erase(v);
deba@458:     ///\endcode
deba@458:     ///or you can also iterate through its elements.
deba@458:     ///\code
deba@458:     ///double s=0;
deba@459:     ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
deba@459:     ///  s+=*i * primal(i);
deba@458:     ///\endcode
deba@459:     ///(This code computes the primal value of the expression).
deba@458:     ///- Numbers (<tt>double</tt>'s)
deba@458:     ///and variables (\ref Col "Col"s) directly convert to an
deba@458:     ///\ref Expr and the usual linear operations are defined, so
deba@458:     ///\code
deba@458:     ///v+w
deba@458:     ///2*v-3.12*(v-w/2)+2
deba@458:     ///v*2.1+(3*v+(v*12+w+6)*3)/2
deba@458:     ///\endcode
deba@459:     ///are valid expressions.
deba@458:     ///The usual assignment operations are also defined.
deba@458:     ///\code
deba@458:     ///e=v+w;
deba@458:     ///e+=2*v-3.12*(v-w/2)+2;
deba@458:     ///e*=3.4;
deba@458:     ///e/=5;
deba@458:     ///\endcode
deba@459:     ///- The constant member can be set and read by dereference
deba@459:     ///  operator (unary *)
deba@459:     ///
deba@458:     ///\code
deba@459:     ///*e=12;
deba@459:     ///double c=*e;
deba@458:     ///\endcode
deba@458:     ///
deba@458:     ///\sa Constr
deba@459:     class Expr {
deba@459:       friend class LpBase;
deba@458:     public:
deba@459:       /// The key type of the expression
deba@459:       typedef LpBase::Col Key;
deba@459:       /// The value type of the expression
deba@459:       typedef LpBase::Value Value;
deba@458: 
deba@458:     protected:
deba@459:       Value const_comp;
deba@459:       std::map<int, Value> comps;
deba@458: 
deba@458:     public:
deba@459:       typedef True SolverExpr;
deba@459:       /// Default constructor
deba@459:       
deba@459:       /// Construct an empty expression, the coefficients and
deba@459:       /// the constant component are initialized to zero.
deba@459:       Expr() : const_comp(0) {}
deba@459:       /// Construct an expression from a column
deba@459: 
deba@459:       /// Construct an expression, which has a term with \c c variable
deba@459:       /// and 1.0 coefficient.
deba@459:       Expr(const Col &c) : const_comp(0) {
deba@459:         typedef std::map<int, Value>::value_type pair_type;
deba@459:         comps.insert(pair_type(id(c), 1));
deba@458:       }
deba@459:       /// Construct an expression from a constant
deba@459: 
deba@459:       /// Construct an expression, which's constant component is \c v.
deba@459:       ///
deba@458:       Expr(const Value &v) : const_comp(v) {}
deba@459:       /// Returns the coefficient of the column
deba@459:       Value operator[](const Col& c) const {
deba@459:         std::map<int, Value>::const_iterator it=comps.find(id(c));
deba@459:         if (it != comps.end()) {
deba@459:           return it->second;
deba@459:         } else {
deba@459:           return 0;
deba@458:         }
deba@458:       }
deba@459:       /// Returns the coefficient of the column
deba@459:       Value& operator[](const Col& c) {
deba@459:         return comps[id(c)];
deba@459:       }
deba@459:       /// Sets the coefficient of the column
deba@459:       void set(const Col &c, const Value &v) {
deba@459:         if (v != 0.0) {
deba@459:           typedef std::map<int, Value>::value_type pair_type;
deba@459:           comps.insert(pair_type(id(c), v));
deba@459:         } else {
deba@459:           comps.erase(id(c));
deba@459:         }
deba@459:       }
deba@459:       /// Returns the constant component of the expression
deba@459:       Value& operator*() { return const_comp; }
deba@459:       /// Returns the constant component of the expression
deba@459:       const Value& operator*() const { return const_comp; }
deba@459:       /// \brief Removes the coefficients which's absolute value does
deba@459:       /// not exceed \c epsilon. It also sets to zero the constant
deba@459:       /// component, if it does not exceed epsilon in absolute value.
deba@459:       void simplify(Value epsilon = 0.0) {
deba@459:         std::map<int, Value>::iterator it=comps.begin();
deba@459:         while (it != comps.end()) {
deba@459:           std::map<int, Value>::iterator jt=it;
deba@459:           ++jt;
deba@459:           if (std::fabs((*it).second) <= epsilon) comps.erase(it);
deba@459:           it=jt;
deba@459:         }
deba@459:         if (std::fabs(const_comp) <= epsilon) const_comp = 0;
deba@458:       }
deba@458: 
deba@459:       void simplify(Value epsilon = 0.0) const {
deba@459:         const_cast<Expr*>(this)->simplify(epsilon);
deba@458:       }
deba@458: 
deba@458:       ///Sets all coefficients and the constant component to 0.
deba@458:       void clear() {
deba@459:         comps.clear();
deba@458:         const_comp=0;
deba@458:       }
deba@458: 
deba@459:       ///Compound assignment
deba@458:       Expr &operator+=(const Expr &e) {
deba@459:         for (std::map<int, Value>::const_iterator it=e.comps.begin();
deba@459:              it!=e.comps.end(); ++it)
deba@459:           comps[it->first]+=it->second;
deba@458:         const_comp+=e.const_comp;
deba@458:         return *this;
deba@458:       }
deba@459:       ///Compound assignment
deba@458:       Expr &operator-=(const Expr &e) {
deba@459:         for (std::map<int, Value>::const_iterator it=e.comps.begin();
deba@459:              it!=e.comps.end(); ++it)
deba@459:           comps[it->first]-=it->second;
deba@458:         const_comp-=e.const_comp;
deba@458:         return *this;
deba@458:       }
deba@459:       ///Multiply with a constant
deba@459:       Expr &operator*=(const Value &v) {
deba@459:         for (std::map<int, Value>::iterator it=comps.begin();
deba@459:              it!=comps.end(); ++it)
deba@459:           it->second*=v;
deba@459:         const_comp*=v;
deba@458:         return *this;
deba@458:       }
deba@459:       ///Division with a constant
deba@458:       Expr &operator/=(const Value &c) {
deba@459:         for (std::map<int, Value>::iterator it=comps.begin();
deba@459:              it!=comps.end(); ++it)
deba@459:           it->second/=c;
deba@458:         const_comp/=c;
deba@458:         return *this;
deba@458:       }
deba@458: 
deba@459:       ///Iterator over the expression
deba@459:       
deba@459:       ///The iterator iterates over the terms of the expression. 
deba@459:       /// 
deba@459:       ///\code
deba@459:       ///double s=0;
deba@459:       ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
deba@459:       ///  s+= *i * primal(i);
deba@459:       ///\endcode
deba@459:       class CoeffIt {
deba@459:       private:
deba@459: 
deba@459:         std::map<int, Value>::iterator _it, _end;
deba@459: 
deba@459:       public:
deba@459: 
deba@459:         /// Sets the iterator to the first term
deba@459:         
deba@459:         /// Sets the iterator to the first term of the expression.
deba@459:         ///
deba@459:         CoeffIt(Expr& e)
deba@459:           : _it(e.comps.begin()), _end(e.comps.end()){}
deba@459: 
deba@459:         /// Convert the iterator to the column of the term
deba@459:         operator Col() const {
deba@459:           return colFromId(_it->first);
deba@459:         }
deba@459: 
deba@459:         /// Returns the coefficient of the term
deba@459:         Value& operator*() { return _it->second; }
deba@459: 
deba@459:         /// Returns the coefficient of the term
deba@459:         const Value& operator*() const { return _it->second; }
deba@459:         /// Next term
deba@459:         
deba@459:         /// Assign the iterator to the next term.
deba@459:         ///
deba@459:         CoeffIt& operator++() { ++_it; return *this; }
deba@459: 
deba@459:         /// Equality operator
deba@459:         bool operator==(Invalid) const { return _it == _end; }
deba@459:         /// Inequality operator
deba@459:         bool operator!=(Invalid) const { return _it != _end; }
deba@459:       };
deba@459: 
deba@459:       /// Const iterator over the expression
deba@459:       
deba@459:       ///The iterator iterates over the terms of the expression. 
deba@459:       /// 
deba@459:       ///\code
deba@459:       ///double s=0;
deba@459:       ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
deba@459:       ///  s+=*i * primal(i);
deba@459:       ///\endcode
deba@459:       class ConstCoeffIt {
deba@459:       private:
deba@459: 
deba@459:         std::map<int, Value>::const_iterator _it, _end;
deba@459: 
deba@459:       public:
deba@459: 
deba@459:         /// Sets the iterator to the first term
deba@459:         
deba@459:         /// Sets the iterator to the first term of the expression.
deba@459:         ///
deba@459:         ConstCoeffIt(const Expr& e)
deba@459:           : _it(e.comps.begin()), _end(e.comps.end()){}
deba@459: 
deba@459:         /// Convert the iterator to the column of the term
deba@459:         operator Col() const {
deba@459:           return colFromId(_it->first);
deba@459:         }
deba@459: 
deba@459:         /// Returns the coefficient of the term
deba@459:         const Value& operator*() const { return _it->second; }
deba@459: 
deba@459:         /// Next term
deba@459:         
deba@459:         /// Assign the iterator to the next term.
deba@459:         ///
deba@459:         ConstCoeffIt& operator++() { ++_it; return *this; }
deba@459: 
deba@459:         /// Equality operator
deba@459:         bool operator==(Invalid) const { return _it == _end; }
deba@459:         /// Inequality operator
deba@459:         bool operator!=(Invalid) const { return _it != _end; }
deba@459:       };
deba@459: 
deba@458:     };
deba@458: 
deba@458:     ///Linear constraint
deba@458: 
deba@458:     ///This data stucture represents a linear constraint in the LP.
deba@458:     ///Basically it is a linear expression with a lower or an upper bound
deba@458:     ///(or both). These parts of the constraint can be obtained by the member
deba@458:     ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
deba@458:     ///respectively.
deba@458:     ///There are two ways to construct a constraint.
deba@458:     ///- You can set the linear expression and the bounds directly
deba@458:     ///  by the functions above.
deba@458:     ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
deba@458:     ///  are defined between expressions, or even between constraints whenever
deba@458:     ///  it makes sense. Therefore if \c e and \c f are linear expressions and
deba@458:     ///  \c s and \c t are numbers, then the followings are valid expressions
deba@458:     ///  and thus they can be used directly e.g. in \ref addRow() whenever
deba@458:     ///  it makes sense.
deba@458:     ///\code
deba@458:     ///  e<=s
deba@458:     ///  e<=f
deba@458:     ///  e==f
deba@458:     ///  s<=e<=t
deba@458:     ///  e>=t
deba@458:     ///\endcode
deba@459:     ///\warning The validity of a constraint is checked only at run
deba@459:     ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
deba@459:     ///compile, but will fail an assertion.
deba@458:     class Constr
deba@458:     {
deba@458:     public:
deba@459:       typedef LpBase::Expr Expr;
deba@458:       typedef Expr::Key Key;
deba@458:       typedef Expr::Value Value;
deba@458: 
deba@458:     protected:
deba@458:       Expr _expr;
deba@458:       Value _lb,_ub;
deba@458:     public:
deba@458:       ///\e
deba@458:       Constr() : _expr(), _lb(NaN), _ub(NaN) {}
deba@458:       ///\e
deba@459:       Constr(Value lb, const Expr &e, Value ub) :
deba@458:         _expr(e), _lb(lb), _ub(ub) {}
deba@458:       Constr(const Expr &e) :
deba@458:         _expr(e), _lb(NaN), _ub(NaN) {}
deba@458:       ///\e
deba@458:       void clear()
deba@458:       {
deba@458:         _expr.clear();
deba@458:         _lb=_ub=NaN;
deba@458:       }
deba@458: 
deba@458:       ///Reference to the linear expression
deba@458:       Expr &expr() { return _expr; }
deba@458:       ///Cont reference to the linear expression
deba@458:       const Expr &expr() const { return _expr; }
deba@458:       ///Reference to the lower bound.
deba@458: 
deba@458:       ///\return
deba@458:       ///- \ref INF "INF": the constraint is lower unbounded.
deba@458:       ///- \ref NaN "NaN": lower bound has not been set.
deba@458:       ///- finite number: the lower bound
deba@458:       Value &lowerBound() { return _lb; }
deba@458:       ///The const version of \ref lowerBound()
deba@458:       const Value &lowerBound() const { return _lb; }
deba@458:       ///Reference to the upper bound.
deba@458: 
deba@458:       ///\return
deba@458:       ///- \ref INF "INF": the constraint is upper unbounded.
deba@458:       ///- \ref NaN "NaN": upper bound has not been set.
deba@458:       ///- finite number: the upper bound
deba@458:       Value &upperBound() { return _ub; }
deba@458:       ///The const version of \ref upperBound()
deba@458:       const Value &upperBound() const { return _ub; }
deba@458:       ///Is the constraint lower bounded?
deba@458:       bool lowerBounded() const {
alpar@487:         return _lb != -INF && !isNaN(_lb);
deba@458:       }
deba@458:       ///Is the constraint upper bounded?
deba@458:       bool upperBounded() const {
alpar@487:         return _ub != INF && !isNaN(_ub);
deba@458:       }
deba@458: 
deba@458:     };
deba@458: 
deba@458:     ///Linear expression of rows
deba@458: 
deba@458:     ///This data structure represents a column of the matrix,
deba@458:     ///thas is it strores a linear expression of the dual variables
deba@458:     ///(\ref Row "Row"s).
deba@458:     ///
deba@458:     ///There are several ways to access and modify the contents of this
deba@458:     ///container.
deba@458:     ///\code
deba@458:     ///e[v]=5;
deba@458:     ///e[v]+=12;
deba@458:     ///e.erase(v);
deba@458:     ///\endcode
deba@458:     ///or you can also iterate through its elements.
deba@458:     ///\code
deba@458:     ///double s=0;
deba@459:     ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
deba@459:     ///  s+=*i;
deba@458:     ///\endcode
deba@458:     ///(This code computes the sum of all coefficients).
deba@458:     ///- Numbers (<tt>double</tt>'s)
deba@458:     ///and variables (\ref Row "Row"s) directly convert to an
deba@458:     ///\ref DualExpr and the usual linear operations are defined, so
deba@458:     ///\code
deba@458:     ///v+w
deba@458:     ///2*v-3.12*(v-w/2)
deba@458:     ///v*2.1+(3*v+(v*12+w)*3)/2
deba@458:     ///\endcode
deba@459:     ///are valid \ref DualExpr dual expressions.
deba@458:     ///The usual assignment operations are also defined.
deba@458:     ///\code
deba@458:     ///e=v+w;
deba@458:     ///e+=2*v-3.12*(v-w/2);
deba@458:     ///e*=3.4;
deba@458:     ///e/=5;
deba@458:     ///\endcode
deba@458:     ///
deba@458:     ///\sa Expr
deba@459:     class DualExpr {
deba@459:       friend class LpBase;
deba@458:     public:
deba@459:       /// The key type of the expression
deba@459:       typedef LpBase::Row Key;
deba@459:       /// The value type of the expression
deba@459:       typedef LpBase::Value Value;
deba@458: 
deba@458:     protected:
deba@459:       std::map<int, Value> comps;
deba@458: 
deba@458:     public:
deba@459:       typedef True SolverExpr;
deba@459:       /// Default constructor
deba@459:       
deba@459:       /// Construct an empty expression, the coefficients are
deba@459:       /// initialized to zero.
deba@459:       DualExpr() {}
deba@459:       /// Construct an expression from a row
deba@459: 
deba@459:       /// Construct an expression, which has a term with \c r dual
deba@459:       /// variable and 1.0 coefficient.
deba@459:       DualExpr(const Row &r) {
deba@459:         typedef std::map<int, Value>::value_type pair_type;
deba@459:         comps.insert(pair_type(id(r), 1));
deba@458:       }
deba@459:       /// Returns the coefficient of the row
deba@459:       Value operator[](const Row& r) const {
deba@459:         std::map<int, Value>::const_iterator it = comps.find(id(r));
deba@459:         if (it != comps.end()) {
deba@459:           return it->second;
deba@459:         } else {
deba@459:           return 0;
deba@459:         }
deba@458:       }
deba@459:       /// Returns the coefficient of the row
deba@459:       Value& operator[](const Row& r) {
deba@459:         return comps[id(r)];
deba@459:       }
deba@459:       /// Sets the coefficient of the row
deba@459:       void set(const Row &r, const Value &v) {
deba@459:         if (v != 0.0) {
deba@459:           typedef std::map<int, Value>::value_type pair_type;
deba@459:           comps.insert(pair_type(id(r), v));
deba@459:         } else {
deba@459:           comps.erase(id(r));
deba@459:         }
deba@459:       }
deba@459:       /// \brief Removes the coefficients which's absolute value does
deba@459:       /// not exceed \c epsilon. 
deba@459:       void simplify(Value epsilon = 0.0) {
deba@459:         std::map<int, Value>::iterator it=comps.begin();
deba@459:         while (it != comps.end()) {
deba@459:           std::map<int, Value>::iterator jt=it;
deba@459:           ++jt;
deba@459:           if (std::fabs((*it).second) <= epsilon) comps.erase(it);
deba@459:           it=jt;
deba@458:         }
deba@458:       }
deba@458: 
deba@459:       void simplify(Value epsilon = 0.0) const {
deba@459:         const_cast<DualExpr*>(this)->simplify(epsilon);
deba@458:       }
deba@458: 
deba@458:       ///Sets all coefficients to 0.
deba@458:       void clear() {
deba@459:         comps.clear();
deba@459:       }
deba@459:       ///Compound assignment
deba@459:       DualExpr &operator+=(const DualExpr &e) {
deba@459:         for (std::map<int, Value>::const_iterator it=e.comps.begin();
deba@459:              it!=e.comps.end(); ++it)
deba@459:           comps[it->first]+=it->second;
deba@459:         return *this;
deba@459:       }
deba@459:       ///Compound assignment
deba@459:       DualExpr &operator-=(const DualExpr &e) {
deba@459:         for (std::map<int, Value>::const_iterator it=e.comps.begin();
deba@459:              it!=e.comps.end(); ++it)
deba@459:           comps[it->first]-=it->second;
deba@459:         return *this;
deba@459:       }
deba@459:       ///Multiply with a constant
deba@459:       DualExpr &operator*=(const Value &v) {
deba@459:         for (std::map<int, Value>::iterator it=comps.begin();
deba@459:              it!=comps.end(); ++it)
deba@459:           it->second*=v;
deba@459:         return *this;
deba@459:       }
deba@459:       ///Division with a constant
deba@459:       DualExpr &operator/=(const Value &v) {
deba@459:         for (std::map<int, Value>::iterator it=comps.begin();
deba@459:              it!=comps.end(); ++it)
deba@459:           it->second/=v;
deba@459:         return *this;
deba@458:       }
deba@458: 
deba@459:       ///Iterator over the expression
deba@459:       
deba@459:       ///The iterator iterates over the terms of the expression. 
deba@459:       /// 
deba@459:       ///\code
deba@459:       ///double s=0;
deba@459:       ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
deba@459:       ///  s+= *i * dual(i);
deba@459:       ///\endcode
deba@459:       class CoeffIt {
deba@459:       private:
deba@459: 
deba@459:         std::map<int, Value>::iterator _it, _end;
deba@459: 
deba@459:       public:
deba@459: 
deba@459:         /// Sets the iterator to the first term
deba@459:         
deba@459:         /// Sets the iterator to the first term of the expression.
deba@459:         ///
deba@459:         CoeffIt(DualExpr& e)
deba@459:           : _it(e.comps.begin()), _end(e.comps.end()){}
deba@459: 
deba@459:         /// Convert the iterator to the row of the term
deba@459:         operator Row() const {
deba@459:           return rowFromId(_it->first);
deba@459:         }
deba@459: 
deba@459:         /// Returns the coefficient of the term
deba@459:         Value& operator*() { return _it->second; }
deba@459: 
deba@459:         /// Returns the coefficient of the term
deba@459:         const Value& operator*() const { return _it->second; }
deba@459: 
deba@459:         /// Next term
deba@459:         
deba@459:         /// Assign the iterator to the next term.
deba@459:         ///
deba@459:         CoeffIt& operator++() { ++_it; return *this; }
deba@459: 
deba@459:         /// Equality operator
deba@459:         bool operator==(Invalid) const { return _it == _end; }
deba@459:         /// Inequality operator
deba@459:         bool operator!=(Invalid) const { return _it != _end; }
deba@459:       };
deba@459: 
deba@459:       ///Iterator over the expression
deba@459:       
deba@459:       ///The iterator iterates over the terms of the expression. 
deba@459:       /// 
deba@459:       ///\code
deba@459:       ///double s=0;
deba@459:       ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
deba@459:       ///  s+= *i * dual(i);
deba@459:       ///\endcode
deba@459:       class ConstCoeffIt {
deba@459:       private:
deba@459: 
deba@459:         std::map<int, Value>::const_iterator _it, _end;
deba@459: 
deba@459:       public:
deba@459: 
deba@459:         /// Sets the iterator to the first term
deba@459:         
deba@459:         /// Sets the iterator to the first term of the expression.
deba@459:         ///
deba@459:         ConstCoeffIt(const DualExpr& e)
deba@459:           : _it(e.comps.begin()), _end(e.comps.end()){}
deba@459: 
deba@459:         /// Convert the iterator to the row of the term
deba@459:         operator Row() const {
deba@459:           return rowFromId(_it->first);
deba@459:         }
deba@459: 
deba@459:         /// Returns the coefficient of the term
deba@459:         const Value& operator*() const { return _it->second; }
deba@459: 
deba@459:         /// Next term
deba@459:         
deba@459:         /// Assign the iterator to the next term.
deba@459:         ///
deba@459:         ConstCoeffIt& operator++() { ++_it; return *this; }
deba@459: 
deba@459:         /// Equality operator
deba@459:         bool operator==(Invalid) const { return _it == _end; }
deba@459:         /// Inequality operator
deba@459:         bool operator!=(Invalid) const { return _it != _end; }
deba@459:       };
deba@458:     };
deba@458: 
deba@458: 
deba@459:   protected:
deba@458: 
deba@459:     class InsertIterator {
deba@459:     private:
deba@459: 
deba@459:       std::map<int, Value>& _host;
deba@459:       const _solver_bits::VarIndex& _index;
deba@459: 
deba@458:     public:
deba@458: 
deba@458:       typedef std::output_iterator_tag iterator_category;
deba@458:       typedef void difference_type;
deba@458:       typedef void value_type;
deba@458:       typedef void reference;
deba@458:       typedef void pointer;
deba@458: 
deba@459:       InsertIterator(std::map<int, Value>& host,
deba@459:                    const _solver_bits::VarIndex& index)
deba@459:         : _host(host), _index(index) {}
deba@458: 
deba@459:       InsertIterator& operator=(const std::pair<int, Value>& value) {
deba@459:         typedef std::map<int, Value>::value_type pair_type;
deba@459:         _host.insert(pair_type(_index[value.first], value.second));
deba@458:         return *this;
deba@458:       }
deba@458: 
deba@459:       InsertIterator& operator*() { return *this; }
deba@459:       InsertIterator& operator++() { return *this; }
deba@459:       InsertIterator operator++(int) { return *this; }
deba@458: 
deba@458:     };
deba@458: 
deba@459:     class ExprIterator {
deba@459:     private:
deba@459:       std::map<int, Value>::const_iterator _host_it;
deba@459:       const _solver_bits::VarIndex& _index;
deba@458:     public:
deba@458: 
deba@459:       typedef std::bidirectional_iterator_tag iterator_category;
deba@459:       typedef std::ptrdiff_t difference_type;
deba@458:       typedef const std::pair<int, Value> value_type;
deba@458:       typedef value_type reference;
deba@459: 
deba@458:       class pointer {
deba@458:       public:
deba@458:         pointer(value_type& _value) : value(_value) {}
deba@458:         value_type* operator->() { return &value; }
deba@458:       private:
deba@458:         value_type value;
deba@458:       };
deba@458: 
deba@459:       ExprIterator(const std::map<int, Value>::const_iterator& host_it,
deba@459:                    const _solver_bits::VarIndex& index)
deba@459:         : _host_it(host_it), _index(index) {}
deba@458: 
deba@458:       reference operator*() {
deba@459:         return std::make_pair(_index(_host_it->first), _host_it->second);
deba@458:       }
deba@458: 
deba@458:       pointer operator->() {
deba@458:         return pointer(operator*());
deba@458:       }
deba@458: 
deba@459:       ExprIterator& operator++() { ++_host_it; return *this; }
deba@459:       ExprIterator operator++(int) {
deba@459:         ExprIterator tmp(*this); ++_host_it; return tmp;
deba@458:       }
deba@458: 
deba@459:       ExprIterator& operator--() { --_host_it; return *this; }
deba@459:       ExprIterator operator--(int) {
deba@459:         ExprIterator tmp(*this); --_host_it; return tmp;
deba@458:       }
deba@458: 
deba@459:       bool operator==(const ExprIterator& it) const {
deba@459:         return _host_it == it._host_it;
deba@458:       }
deba@458: 
deba@459:       bool operator!=(const ExprIterator& it) const {
deba@459:         return _host_it != it._host_it;
deba@458:       }
deba@458: 
deba@458:     };
deba@458: 
deba@458:   protected:
deba@458: 
deba@459:     //Abstract virtual functions
deba@458: 
deba@459:     virtual int _addColId(int col) { return cols.addIndex(col); }
deba@459:     virtual int _addRowId(int row) { return rows.addIndex(row); }
deba@458: 
deba@459:     virtual void _eraseColId(int col) { cols.eraseIndex(col); }
deba@459:     virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
deba@458: 
deba@458:     virtual int _addCol() = 0;
deba@458:     virtual int _addRow() = 0;
deba@458: 
deba@458:     virtual void _eraseCol(int col) = 0;
deba@458:     virtual void _eraseRow(int row) = 0;
deba@458: 
deba@459:     virtual void _getColName(int col, std::string& name) const = 0;
deba@459:     virtual void _setColName(int col, const std::string& name) = 0;
deba@458:     virtual int _colByName(const std::string& name) const = 0;
deba@458: 
deba@459:     virtual void _getRowName(int row, std::string& name) const = 0;
deba@459:     virtual void _setRowName(int row, const std::string& name) = 0;
deba@459:     virtual int _rowByName(const std::string& name) const = 0;
deba@459: 
deba@459:     virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
deba@459:     virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
deba@459: 
deba@459:     virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
deba@459:     virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
deba@459: 
deba@458:     virtual void _setCoeff(int row, int col, Value value) = 0;
deba@458:     virtual Value _getCoeff(int row, int col) const = 0;
deba@459: 
deba@458:     virtual void _setColLowerBound(int i, Value value) = 0;
deba@458:     virtual Value _getColLowerBound(int i) const = 0;
deba@459: 
deba@458:     virtual void _setColUpperBound(int i, Value value) = 0;
deba@458:     virtual Value _getColUpperBound(int i) const = 0;
deba@459: 
deba@459:     virtual void _setRowLowerBound(int i, Value value) = 0;
deba@459:     virtual Value _getRowLowerBound(int i) const = 0;
deba@459: 
deba@459:     virtual void _setRowUpperBound(int i, Value value) = 0;
deba@459:     virtual Value _getRowUpperBound(int i) const = 0;
deba@459: 
deba@459:     virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
deba@459:     virtual void _getObjCoeffs(InsertIterator b) const = 0;
deba@458: 
deba@458:     virtual void _setObjCoeff(int i, Value obj_coef) = 0;
deba@458:     virtual Value _getObjCoeff(int i) const = 0;
deba@458: 
deba@459:     virtual void _setSense(Sense) = 0;
deba@459:     virtual Sense _getSense() const = 0;
deba@458: 
deba@459:     virtual void _clear() = 0;
deba@458: 
deba@459:     virtual const char* _solverName() const = 0;
deba@458: 
deba@568:     virtual void _messageLevel(MessageLevel level) = 0;
deba@568: 
deba@458:     //Own protected stuff
deba@458: 
deba@458:     //Constant component of the objective function
deba@458:     Value obj_const_comp;
deba@458: 
deba@459:     LpBase() : rows(), cols(), obj_const_comp(0) {}
deba@459: 
deba@458:   public:
deba@458: 
deba@459:     /// Virtual destructor
deba@459:     virtual ~LpBase() {}
deba@458: 
deba@459:     ///Gives back the name of the solver.
deba@459:     const char* solverName() const {return _solverName();}
deba@458: 
kpeter@576:     ///\name Build Up and Modify the LP
deba@458: 
deba@458:     ///@{
deba@458: 
deba@458:     ///Add a new empty column (i.e a new variable) to the LP
deba@459:     Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
deba@458: 
deba@459:     ///\brief Adds several new columns (i.e variables) at once
deba@458:     ///
deba@459:     ///This magic function takes a container as its argument and fills
deba@459:     ///its elements with new columns (i.e. variables)
deba@458:     ///\param t can be
deba@458:     ///- a standard STL compatible iterable container with
deba@459:     ///\ref Col as its \c values_type like
deba@458:     ///\code
deba@459:     ///std::vector<LpBase::Col>
deba@459:     ///std::list<LpBase::Col>
deba@458:     ///\endcode
deba@458:     ///- a standard STL compatible iterable container with
deba@459:     ///\ref Col as its \c mapped_type like
deba@458:     ///\code
deba@459:     ///std::map<AnyType,LpBase::Col>
deba@458:     ///\endcode
deba@458:     ///- an iterable lemon \ref concepts::WriteMap "write map" like
deba@458:     ///\code
deba@459:     ///ListGraph::NodeMap<LpBase::Col>
deba@459:     ///ListGraph::ArcMap<LpBase::Col>
deba@458:     ///\endcode
deba@458:     ///\return The number of the created column.
deba@458: #ifdef DOXYGEN
deba@458:     template<class T>
deba@458:     int addColSet(T &t) { return 0;}
deba@458: #else
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::value_type::LpCol,int>::type
deba@458:     addColSet(T &t,dummy<0> = 0) {
deba@458:       int s=0;
deba@458:       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
deba@458:       return s;
deba@458:     }
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::value_type::second_type::LpCol,
deba@458:                        int>::type
deba@458:     addColSet(T &t,dummy<1> = 1) {
deba@458:       int s=0;
deba@458:       for(typename T::iterator i=t.begin();i!=t.end();++i) {
deba@458:         i->second=addCol();
deba@458:         s++;
deba@458:       }
deba@458:       return s;
deba@458:     }
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::MapIt::Value::LpCol,
deba@458:                        int>::type
deba@458:     addColSet(T &t,dummy<2> = 2) {
deba@458:       int s=0;
deba@458:       for(typename T::MapIt i(t); i!=INVALID; ++i)
deba@458:         {
deba@458:           i.set(addCol());
deba@458:           s++;
deba@458:         }
deba@458:       return s;
deba@458:     }
deba@458: #endif
deba@458: 
deba@458:     ///Set a column (i.e a dual constraint) of the LP
deba@458: 
deba@458:     ///\param c is the column to be modified
deba@458:     ///\param e is a dual linear expression (see \ref DualExpr)
deba@458:     ///a better one.
deba@459:     void col(Col c, const DualExpr &e) {
deba@458:       e.simplify();
deba@471:       _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), rows),
deba@471:                     ExprIterator(e.comps.end(), rows));
deba@458:     }
deba@458: 
deba@458:     ///Get a column (i.e a dual constraint) of the LP
deba@458: 
deba@459:     ///\param c is the column to get
deba@458:     ///\return the dual expression associated to the column
deba@458:     DualExpr col(Col c) const {
deba@458:       DualExpr e;
deba@459:       _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
deba@458:       return e;
deba@458:     }
deba@458: 
deba@458:     ///Add a new column to the LP
deba@458: 
deba@458:     ///\param e is a dual linear expression (see \ref DualExpr)
deba@459:     ///\param o is the corresponding component of the objective
deba@458:     ///function. It is 0 by default.
deba@458:     ///\return The created column.
deba@458:     Col addCol(const DualExpr &e, Value o = 0) {
deba@458:       Col c=addCol();
deba@458:       col(c,e);
deba@458:       objCoeff(c,o);
deba@458:       return c;
deba@458:     }
deba@458: 
deba@458:     ///Add a new empty row (i.e a new constraint) to the LP
deba@458: 
deba@458:     ///This function adds a new empty row (i.e a new constraint) to the LP.
deba@458:     ///\return The created row
deba@459:     Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
deba@458: 
deba@459:     ///\brief Add several new rows (i.e constraints) at once
deba@458:     ///
deba@459:     ///This magic function takes a container as its argument and fills
deba@459:     ///its elements with new row (i.e. variables)
deba@458:     ///\param t can be
deba@458:     ///- a standard STL compatible iterable container with
deba@459:     ///\ref Row as its \c values_type like
deba@458:     ///\code
deba@459:     ///std::vector<LpBase::Row>
deba@459:     ///std::list<LpBase::Row>
deba@458:     ///\endcode
deba@458:     ///- a standard STL compatible iterable container with
deba@459:     ///\ref Row as its \c mapped_type like
deba@458:     ///\code
deba@459:     ///std::map<AnyType,LpBase::Row>
deba@458:     ///\endcode
deba@458:     ///- an iterable lemon \ref concepts::WriteMap "write map" like
deba@458:     ///\code
deba@459:     ///ListGraph::NodeMap<LpBase::Row>
deba@459:     ///ListGraph::ArcMap<LpBase::Row>
deba@458:     ///\endcode
deba@458:     ///\return The number of rows created.
deba@458: #ifdef DOXYGEN
deba@458:     template<class T>
deba@458:     int addRowSet(T &t) { return 0;}
deba@458: #else
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::value_type::LpRow,int>::type
deba@459:     addRowSet(T &t, dummy<0> = 0) {
deba@458:       int s=0;
deba@458:       for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
deba@458:       return s;
deba@458:     }
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::value_type::second_type::LpRow, int>::type
deba@459:     addRowSet(T &t, dummy<1> = 1) {
deba@458:       int s=0;
deba@458:       for(typename T::iterator i=t.begin();i!=t.end();++i) {
deba@458:         i->second=addRow();
deba@458:         s++;
deba@458:       }
deba@458:       return s;
deba@458:     }
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::MapIt::Value::LpRow, int>::type
deba@459:     addRowSet(T &t, dummy<2> = 2) {
deba@458:       int s=0;
deba@458:       for(typename T::MapIt i(t); i!=INVALID; ++i)
deba@458:         {
deba@458:           i.set(addRow());
deba@458:           s++;
deba@458:         }
deba@458:       return s;
deba@458:     }
deba@458: #endif
deba@458: 
deba@458:     ///Set a row (i.e a constraint) of the LP
deba@458: 
deba@458:     ///\param r is the row to be modified
deba@458:     ///\param l is lower bound (-\ref INF means no bound)
deba@458:     ///\param e is a linear expression (see \ref Expr)
deba@458:     ///\param u is the upper bound (\ref INF means no bound)
deba@458:     void row(Row r, Value l, const Expr &e, Value u) {
deba@458:       e.simplify();
deba@459:       _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
deba@459:                     ExprIterator(e.comps.end(), cols));
deba@459:       _setRowLowerBound(rows(id(r)),l - *e);
deba@459:       _setRowUpperBound(rows(id(r)),u - *e);
deba@458:     }
deba@458: 
deba@458:     ///Set a row (i.e a constraint) of the LP
deba@458: 
deba@458:     ///\param r is the row to be modified
deba@458:     ///\param c is a linear expression (see \ref Constr)
deba@458:     void row(Row r, const Constr &c) {
deba@458:       row(r, c.lowerBounded()?c.lowerBound():-INF,
deba@458:           c.expr(), c.upperBounded()?c.upperBound():INF);
deba@458:     }
deba@458: 
deba@458: 
deba@458:     ///Get a row (i.e a constraint) of the LP
deba@458: 
deba@458:     ///\param r is the row to get
deba@458:     ///\return the expression associated to the row
deba@458:     Expr row(Row r) const {
deba@458:       Expr e;
deba@459:       _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
deba@458:       return e;
deba@458:     }
deba@458: 
deba@458:     ///Add a new row (i.e a new constraint) to the LP
deba@458: 
deba@458:     ///\param l is the lower bound (-\ref INF means no bound)
deba@458:     ///\param e is a linear expression (see \ref Expr)
deba@458:     ///\param u is the upper bound (\ref INF means no bound)
deba@458:     ///\return The created row.
deba@458:     Row addRow(Value l,const Expr &e, Value u) {
deba@458:       Row r=addRow();
deba@458:       row(r,l,e,u);
deba@458:       return r;
deba@458:     }
deba@458: 
deba@458:     ///Add a new row (i.e a new constraint) to the LP
deba@458: 
deba@458:     ///\param c is a linear expression (see \ref Constr)
deba@458:     ///\return The created row.
deba@458:     Row addRow(const Constr &c) {
deba@458:       Row r=addRow();
deba@458:       row(r,c);
deba@458:       return r;
deba@458:     }
deba@459:     ///Erase a column (i.e a variable) from the LP
deba@458: 
deba@459:     ///\param c is the column to be deleted
deba@459:     void erase(Col c) {
deba@459:       _eraseCol(cols(id(c)));
deba@459:       _eraseColId(cols(id(c)));
deba@458:     }
deba@459:     ///Erase a row (i.e a constraint) from the LP
deba@458: 
deba@458:     ///\param r is the row to be deleted
deba@459:     void erase(Row r) {
deba@459:       _eraseRow(rows(id(r)));
deba@459:       _eraseRowId(rows(id(r)));
deba@458:     }
deba@458: 
deba@458:     /// Get the name of a column
deba@458: 
deba@459:     ///\param c is the coresponding column
deba@458:     ///\return The name of the colunm
deba@458:     std::string colName(Col c) const {
deba@458:       std::string name;
deba@459:       _getColName(cols(id(c)), name);
deba@458:       return name;
deba@458:     }
deba@458: 
deba@458:     /// Set the name of a column
deba@458: 
deba@459:     ///\param c is the coresponding column
deba@458:     ///\param name The name to be given
deba@458:     void colName(Col c, const std::string& name) {
deba@459:       _setColName(cols(id(c)), name);
deba@458:     }
deba@458: 
deba@458:     /// Get the column by its name
deba@458: 
deba@458:     ///\param name The name of the column
deba@458:     ///\return the proper column or \c INVALID
deba@458:     Col colByName(const std::string& name) const {
deba@458:       int k = _colByName(name);
deba@459:       return k != -1 ? Col(cols[k]) : Col(INVALID);
deba@459:     }
deba@459: 
deba@459:     /// Get the name of a row
deba@459: 
deba@459:     ///\param r is the coresponding row
deba@459:     ///\return The name of the row
deba@459:     std::string rowName(Row r) const {
deba@459:       std::string name;
deba@459:       _getRowName(rows(id(r)), name);
deba@459:       return name;
deba@459:     }
deba@459: 
deba@459:     /// Set the name of a row
deba@459: 
deba@459:     ///\param r is the coresponding row
deba@459:     ///\param name The name to be given
deba@459:     void rowName(Row r, const std::string& name) {
deba@459:       _setRowName(rows(id(r)), name);
deba@459:     }
deba@459: 
deba@459:     /// Get the row by its name
deba@459: 
deba@459:     ///\param name The name of the row
deba@459:     ///\return the proper row or \c INVALID
deba@459:     Row rowByName(const std::string& name) const {
deba@459:       int k = _rowByName(name);
deba@459:       return k != -1 ? Row(rows[k]) : Row(INVALID);
deba@458:     }
deba@458: 
deba@458:     /// Set an element of the coefficient matrix of the LP
deba@458: 
deba@458:     ///\param r is the row of the element to be modified
deba@459:     ///\param c is the column of the element to be modified
deba@458:     ///\param val is the new value of the coefficient
deba@458:     void coeff(Row r, Col c, Value val) {
deba@459:       _setCoeff(rows(id(r)),cols(id(c)), val);
deba@458:     }
deba@458: 
deba@458:     /// Get an element of the coefficient matrix of the LP
deba@458: 
deba@459:     ///\param r is the row of the element
deba@459:     ///\param c is the column of the element
deba@458:     ///\return the corresponding coefficient
deba@458:     Value coeff(Row r, Col c) const {
deba@459:       return _getCoeff(rows(id(r)),cols(id(c)));
deba@458:     }
deba@458: 
deba@458:     /// Set the lower bound of a column (i.e a variable)
deba@458: 
deba@458:     /// The lower bound of a variable (column) has to be given by an
deba@458:     /// extended number of type Value, i.e. a finite number of type
deba@458:     /// Value or -\ref INF.
deba@458:     void colLowerBound(Col c, Value value) {
deba@459:       _setColLowerBound(cols(id(c)),value);
deba@458:     }
deba@458: 
deba@458:     /// Get the lower bound of a column (i.e a variable)
deba@458: 
deba@459:     /// This function returns the lower bound for column (variable) \c c
deba@458:     /// (this might be -\ref INF as well).
deba@459:     ///\return The lower bound for column \c c
deba@458:     Value colLowerBound(Col c) const {
deba@459:       return _getColLowerBound(cols(id(c)));
deba@458:     }
deba@458: 
deba@458:     ///\brief Set the lower bound of  several columns
deba@459:     ///(i.e variables) at once
deba@458:     ///
deba@458:     ///This magic function takes a container as its argument
deba@458:     ///and applies the function on all of its elements.
deba@459:     ///The lower bound of a variable (column) has to be given by an
deba@459:     ///extended number of type Value, i.e. a finite number of type
deba@459:     ///Value or -\ref INF.
deba@458: #ifdef DOXYGEN
deba@458:     template<class T>
deba@458:     void colLowerBound(T &t, Value value) { return 0;}
deba@458: #else
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::value_type::LpCol,void>::type
deba@458:     colLowerBound(T &t, Value value,dummy<0> = 0) {
deba@458:       for(typename T::iterator i=t.begin();i!=t.end();++i) {
deba@458:         colLowerBound(*i, value);
deba@458:       }
deba@458:     }
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::value_type::second_type::LpCol,
deba@458:                        void>::type
deba@458:     colLowerBound(T &t, Value value,dummy<1> = 1) {
deba@458:       for(typename T::iterator i=t.begin();i!=t.end();++i) {
deba@458:         colLowerBound(i->second, value);
deba@458:       }
deba@458:     }
deba@458:     template<class T>
deba@459:     typename enable_if<typename T::MapIt::Value::LpCol,
deba@458:                        void>::type
deba@458:     colLowerBound(T &t, Value value,dummy<2> = 2) {
deba@458:       for(typename T::MapIt i(t); i!=INVALID; ++i){
deba@458:         colLowerBound(*i, value);
deba@458:       }
deba@458:     }
deba@458: #endif
deba@458: 
deba@458:     /// Set the upper bound of a column (i.e a variable)
deba@458: 
deba@458:     /// The upper bound of a variable (column) has to be given by an
deba@458:     /// extended number of type Value, i.e. a finite number of type
deba@458:     /// Value or \ref INF.
deba@458:     void colUpperBound(Col c, Value value) {
deba@459:       _setColUpperBound(cols(id(c)),value);
deba@458:     };
deba@458: 
deba@458:     /// Get the upper bound of a column (i.e a variable)
deba@458: 
deba@459:     /// This function returns the upper bound for column (variable) \c c
deba@458:     /// (this might be \ref INF as well).
deba@459:     /// \return The upper bound for column \c c
deba@458:     Value colUpperBound(Col c) const {
deba@459:       return _getColUpperBound(cols(id(c)));
deba@458:     }
deba@458: 
deba@458:     ///\brief Set the upper bound of  several columns
deba@459:     ///(i.e variables) at once
deba@458:     ///
deba@458:     ///This magic function takes a container as its argument
deba@458:     ///and applies the function on all of its elements.
deba@459:     ///The upper bound of a variable (column) has to be given by an
deba@459:     ///extended number of type Value, i.e. a finite number of type
deba@459:     ///Value or \ref INF.
deba@458: #ifdef DOXYGEN
deba@458:     template<class T>
deba@458:     void colUpperBound(T &t, Value value) { return 0;}
deba@458: #else
tapolcai@490:     template<class T1>
tapolcai@490:     typename enable_if<typename T1::value_type::LpCol,void>::type
tapolcai@490:     colUpperBound(T1 &t, Value value,dummy<0> = 0) {
tapolcai@490:       for(typename T1::iterator i=t.begin();i!=t.end();++i) {
deba@458:         colUpperBound(*i, value);
deba@458:       }
deba@458:     }
tapolcai@490:     template<class T1>
tapolcai@490:     typename enable_if<typename T1::value_type::second_type::LpCol,
deba@458:                        void>::type
tapolcai@490:     colUpperBound(T1 &t, Value value,dummy<1> = 1) {
tapolcai@490:       for(typename T1::iterator i=t.begin();i!=t.end();++i) {
deba@458:         colUpperBound(i->second, value);
deba@458:       }
deba@458:     }
tapolcai@490:     template<class T1>
tapolcai@490:     typename enable_if<typename T1::MapIt::Value::LpCol,
deba@458:                        void>::type
tapolcai@490:     colUpperBound(T1 &t, Value value,dummy<2> = 2) {
tapolcai@490:       for(typename T1::MapIt i(t); i!=INVALID; ++i){
deba@458:         colUpperBound(*i, value);
deba@458:       }
deba@458:     }
deba@458: #endif
deba@458: 
deba@458:     /// Set the lower and the upper bounds of a column (i.e a variable)
deba@458: 
deba@458:     /// The lower and the upper bounds of
deba@458:     /// a variable (column) have to be given by an
deba@458:     /// extended number of type Value, i.e. a finite number of type
deba@458:     /// Value, -\ref INF or \ref INF.
deba@458:     void colBounds(Col c, Value lower, Value upper) {
deba@459:       _setColLowerBound(cols(id(c)),lower);
deba@459:       _setColUpperBound(cols(id(c)),upper);
deba@458:     }
deba@458: 
deba@458:     ///\brief Set the lower and the upper bound of several columns
deba@459:     ///(i.e variables) at once
deba@458:     ///
deba@458:     ///This magic function takes a container as its argument
deba@458:     ///and applies the function on all of its elements.
deba@458:     /// The lower and the upper bounds of
deba@458:     /// a variable (column) have to be given by an
deba@458:     /// extended number of type Value, i.e. a finite number of type
deba@458:     /// Value, -\ref INF or \ref INF.
deba@458: #ifdef DOXYGEN
deba@458:     template<class T>
deba@458:     void colBounds(T &t, Value lower, Value upper) { return 0;}
deba@458: #else
tapolcai@490:     template<class T2>
tapolcai@490:     typename enable_if<typename T2::value_type::LpCol,void>::type
tapolcai@490:     colBounds(T2 &t, Value lower, Value upper,dummy<0> = 0) {
tapolcai@490:       for(typename T2::iterator i=t.begin();i!=t.end();++i) {
deba@458:         colBounds(*i, lower, upper);
deba@458:       }
deba@458:     }
tapolcai@490:     template<class T2>
tapolcai@490:     typename enable_if<typename T2::value_type::second_type::LpCol, void>::type
tapolcai@490:     colBounds(T2 &t, Value lower, Value upper,dummy<1> = 1) {
tapolcai@490:       for(typename T2::iterator i=t.begin();i!=t.end();++i) {
deba@458:         colBounds(i->second, lower, upper);
deba@458:       }
deba@458:     }
tapolcai@490:     template<class T2>
tapolcai@490:     typename enable_if<typename T2::MapIt::Value::LpCol, void>::type
tapolcai@490:     colBounds(T2 &t, Value lower, Value upper,dummy<2> = 2) {
tapolcai@490:       for(typename T2::MapIt i(t); i!=INVALID; ++i){
deba@458:         colBounds(*i, lower, upper);
deba@458:       }
deba@458:     }
deba@458: #endif
deba@458: 
deba@459:     /// Set the lower bound of a row (i.e a constraint)
deba@458: 
deba@459:     /// The lower bound of a constraint (row) has to be given by an
deba@459:     /// extended number of type Value, i.e. a finite number of type
deba@459:     /// Value or -\ref INF.
deba@459:     void rowLowerBound(Row r, Value value) {
deba@459:       _setRowLowerBound(rows(id(r)),value);
deba@458:     }
deba@458: 
deba@459:     /// Get the lower bound of a row (i.e a constraint)
deba@458: 
deba@459:     /// This function returns the lower bound for row (constraint) \c c
deba@459:     /// (this might be -\ref INF as well).
deba@459:     ///\return The lower bound for row \c r
deba@459:     Value rowLowerBound(Row r) const {
deba@459:       return _getRowLowerBound(rows(id(r)));
deba@459:     }
deba@459: 
deba@459:     /// Set the upper bound of a row (i.e a constraint)
deba@459: 
deba@459:     /// The upper bound of a constraint (row) has to be given by an
deba@459:     /// extended number of type Value, i.e. a finite number of type
deba@459:     /// Value or -\ref INF.
deba@459:     void rowUpperBound(Row r, Value value) {
deba@459:       _setRowUpperBound(rows(id(r)),value);
deba@459:     }
deba@459: 
deba@459:     /// Get the upper bound of a row (i.e a constraint)
deba@459: 
deba@459:     /// This function returns the upper bound for row (constraint) \c c
deba@459:     /// (this might be -\ref INF as well).
deba@459:     ///\return The upper bound for row \c r
deba@459:     Value rowUpperBound(Row r) const {
deba@459:       return _getRowUpperBound(rows(id(r)));
deba@458:     }
deba@458: 
deba@458:     ///Set an element of the objective function
deba@459:     void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
deba@458: 
deba@458:     ///Get an element of the objective function
deba@459:     Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
deba@458: 
deba@458:     ///Set the objective function
deba@458: 
deba@458:     ///\param e is a linear expression of type \ref Expr.
deba@459:     ///
deba@459:     void obj(const Expr& e) {
deba@459:       _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
deba@459:                     ExprIterator(e.comps.end(), cols));
deba@459:       obj_const_comp = *e;
deba@458:     }
deba@458: 
deba@458:     ///Get the objective function
deba@458: 
deba@459:     ///\return the objective function as a linear expression of type
deba@459:     ///Expr.
deba@458:     Expr obj() const {
deba@458:       Expr e;
deba@459:       _getObjCoeffs(InsertIterator(e.comps, cols));
deba@459:       *e = obj_const_comp;
deba@458:       return e;
deba@458:     }
deba@458: 
deba@458: 
deba@459:     ///Set the direction of optimization
deba@459:     void sense(Sense sense) { _setSense(sense); }
deba@458: 
deba@459:     ///Query the direction of the optimization
deba@459:     Sense sense() const {return _getSense(); }
deba@458: 
deba@459:     ///Set the sense to maximization
deba@459:     void max() { _setSense(MAX); }
deba@459: 
deba@459:     ///Set the sense to maximization
deba@459:     void min() { _setSense(MIN); }
deba@459: 
deba@459:     ///Clears the problem
deba@459:     void clear() { _clear(); }
deba@458: 
deba@568:     /// Sets the message level of the solver
deba@568:     void messageLevel(MessageLevel level) { _messageLevel(level); }
deba@568: 
deba@458:     ///@}
deba@458: 
deba@459:   };
deba@459: 
deba@459:   /// Addition
deba@459: 
deba@459:   ///\relates LpBase::Expr
deba@459:   ///
deba@459:   inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
deba@459:     LpBase::Expr tmp(a);
deba@459:     tmp+=b;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Substraction
deba@459: 
deba@459:   ///\relates LpBase::Expr
deba@459:   ///
deba@459:   inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
deba@459:     LpBase::Expr tmp(a);
deba@459:     tmp-=b;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Multiply with constant
deba@459: 
deba@459:   ///\relates LpBase::Expr
deba@459:   ///
deba@459:   inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
deba@459:     LpBase::Expr tmp(a);
deba@459:     tmp*=b;
deba@459:     return tmp;
deba@459:   }
deba@459: 
deba@459:   ///Multiply with constant
deba@459: 
deba@459:   ///\relates LpBase::Expr
deba@459:   ///
deba@459:   inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
deba@459:     LpBase::Expr tmp(b);
deba@459:     tmp*=a;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Divide with constant
deba@459: 
deba@459:   ///\relates LpBase::Expr
deba@459:   ///
deba@459:   inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
deba@459:     LpBase::Expr tmp(a);
deba@459:     tmp/=b;
deba@459:     return tmp;
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator<=(const LpBase::Expr &e,
deba@459:                                    const LpBase::Expr &f) {
retvari@766:     return LpBase::Constr(0, f - e, LpBase::NaN);
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator<=(const LpBase::Value &e,
deba@459:                                    const LpBase::Expr &f) {
deba@459:     return LpBase::Constr(e, f, LpBase::NaN);
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator<=(const LpBase::Expr &e,
deba@459:                                    const LpBase::Value &f) {
retvari@766:     return LpBase::Constr(LpBase::NaN, e, f);
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator>=(const LpBase::Expr &e,
deba@459:                                    const LpBase::Expr &f) {
retvari@766:     return LpBase::Constr(0, e - f, LpBase::NaN);
deba@459:   }
deba@459: 
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator>=(const LpBase::Value &e,
deba@459:                                    const LpBase::Expr &f) {
deba@459:     return LpBase::Constr(LpBase::NaN, f, e);
deba@459:   }
deba@459: 
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator>=(const LpBase::Expr &e,
deba@459:                                    const LpBase::Value &f) {
retvari@766:     return LpBase::Constr(f, e, LpBase::NaN);
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator==(const LpBase::Expr &e,
deba@459:                                    const LpBase::Value &f) {
deba@459:     return LpBase::Constr(f, e, f);
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator==(const LpBase::Expr &e,
deba@459:                                    const LpBase::Expr &f) {
deba@459:     return LpBase::Constr(0, f - e, 0);
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator<=(const LpBase::Value &n,
deba@459:                                    const LpBase::Constr &c) {
deba@459:     LpBase::Constr tmp(c);
alpar@487:     LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
deba@459:     tmp.lowerBound()=n;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator<=(const LpBase::Constr &c,
deba@459:                                    const LpBase::Value &n)
deba@459:   {
deba@459:     LpBase::Constr tmp(c);
alpar@487:     LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
deba@459:     tmp.upperBound()=n;
deba@459:     return tmp;
deba@459:   }
deba@459: 
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator>=(const LpBase::Value &n,
deba@459:                                    const LpBase::Constr &c) {
deba@459:     LpBase::Constr tmp(c);
alpar@487:     LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
deba@459:     tmp.upperBound()=n;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Create constraint
deba@459: 
deba@459:   ///\relates LpBase::Constr
deba@459:   ///
deba@459:   inline LpBase::Constr operator>=(const LpBase::Constr &c,
deba@459:                                    const LpBase::Value &n)
deba@459:   {
deba@459:     LpBase::Constr tmp(c);
alpar@487:     LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
deba@459:     tmp.lowerBound()=n;
deba@459:     return tmp;
deba@459:   }
deba@459: 
deba@459:   ///Addition
deba@459: 
deba@459:   ///\relates LpBase::DualExpr
deba@459:   ///
deba@459:   inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
deba@459:                                     const LpBase::DualExpr &b) {
deba@459:     LpBase::DualExpr tmp(a);
deba@459:     tmp+=b;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Substraction
deba@459: 
deba@459:   ///\relates LpBase::DualExpr
deba@459:   ///
deba@459:   inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
deba@459:                                     const LpBase::DualExpr &b) {
deba@459:     LpBase::DualExpr tmp(a);
deba@459:     tmp-=b;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Multiply with constant
deba@459: 
deba@459:   ///\relates LpBase::DualExpr
deba@459:   ///
deba@459:   inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
deba@459:                                     const LpBase::Value &b) {
deba@459:     LpBase::DualExpr tmp(a);
deba@459:     tmp*=b;
deba@459:     return tmp;
deba@459:   }
deba@459: 
deba@459:   ///Multiply with constant
deba@459: 
deba@459:   ///\relates LpBase::DualExpr
deba@459:   ///
deba@459:   inline LpBase::DualExpr operator*(const LpBase::Value &a,
deba@459:                                     const LpBase::DualExpr &b) {
deba@459:     LpBase::DualExpr tmp(b);
deba@459:     tmp*=a;
deba@459:     return tmp;
deba@459:   }
deba@459:   ///Divide with constant
deba@459: 
deba@459:   ///\relates LpBase::DualExpr
deba@459:   ///
deba@459:   inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
deba@459:                                     const LpBase::Value &b) {
deba@459:     LpBase::DualExpr tmp(a);
deba@459:     tmp/=b;
deba@459:     return tmp;
deba@459:   }
deba@459: 
deba@459:   /// \ingroup lp_group
deba@459:   ///
deba@459:   /// \brief Common base class for LP solvers
deba@459:   ///
deba@459:   /// This class is an abstract base class for LP solvers. This class
deba@459:   /// provides a full interface for set and modify an LP problem,
deba@459:   /// solve it and retrieve the solution. You can use one of the
deba@459:   /// descendants as a concrete implementation, or the \c Lp
deba@459:   /// default LP solver. However, if you would like to handle LP
deba@459:   /// solvers as reference or pointer in a generic way, you can use
deba@459:   /// this class directly.
deba@459:   class LpSolver : virtual public LpBase {
deba@459:   public:
deba@459: 
deba@459:     /// The problem types for primal and dual problems
deba@459:     enum ProblemType {
kpeter@576:       /// = 0. Feasible solution hasn't been found (but may exist).
deba@459:       UNDEFINED = 0,
kpeter@576:       /// = 1. The problem has no feasible solution.
deba@459:       INFEASIBLE = 1,
kpeter@576:       /// = 2. Feasible solution found.
deba@459:       FEASIBLE = 2,
kpeter@576:       /// = 3. Optimal solution exists and found.
deba@459:       OPTIMAL = 3,
kpeter@576:       /// = 4. The cost function is unbounded.
deba@459:       UNBOUNDED = 4
deba@459:     };
deba@459: 
deba@459:     ///The basis status of variables
deba@459:     enum VarStatus {
deba@459:       /// The variable is in the basis
deba@459:       BASIC, 
deba@459:       /// The variable is free, but not basic
deba@459:       FREE,
deba@459:       /// The variable has active lower bound 
deba@459:       LOWER,
deba@459:       /// The variable has active upper bound
deba@459:       UPPER,
deba@459:       /// The variable is non-basic and fixed
deba@459:       FIXED
deba@459:     };
deba@459: 
deba@459:   protected:
deba@459: 
deba@459:     virtual SolveExitStatus _solve() = 0;
deba@459: 
deba@459:     virtual Value _getPrimal(int i) const = 0;
deba@459:     virtual Value _getDual(int i) const = 0;
deba@459: 
deba@459:     virtual Value _getPrimalRay(int i) const = 0;
deba@459:     virtual Value _getDualRay(int i) const = 0;
deba@459: 
deba@459:     virtual Value _getPrimalValue() const = 0;
deba@459: 
deba@459:     virtual VarStatus _getColStatus(int i) const = 0;
deba@459:     virtual VarStatus _getRowStatus(int i) const = 0;
deba@459: 
deba@459:     virtual ProblemType _getPrimalType() const = 0;
deba@459:     virtual ProblemType _getDualType() const = 0;
deba@459: 
deba@459:   public:
deba@458: 
alpar@528:     ///Allocate a new LP problem instance
alpar@528:     virtual LpSolver* newSolver() const = 0;
alpar@528:     ///Make a copy of the LP problem
alpar@528:     virtual LpSolver* cloneSolver() const = 0;
alpar@528: 
deba@458:     ///\name Solve the LP
deba@458: 
deba@458:     ///@{
deba@458: 
deba@458:     ///\e Solve the LP problem at hand
deba@458:     ///
deba@458:     ///\return The result of the optimization procedure. Possible
deba@458:     ///values and their meanings can be found in the documentation of
deba@458:     ///\ref SolveExitStatus.
deba@458:     SolveExitStatus solve() { return _solve(); }
deba@458: 
deba@458:     ///@}
deba@458: 
kpeter@576:     ///\name Obtain the Solution
deba@458: 
deba@458:     ///@{
deba@458: 
deba@459:     /// The type of the primal problem
deba@459:     ProblemType primalType() const {
deba@459:       return _getPrimalType();
deba@458:     }
deba@458: 
deba@459:     /// The type of the dual problem
deba@459:     ProblemType dualType() const {
deba@459:       return _getDualType();
deba@458:     }
deba@458: 
deba@459:     /// Return the primal value of the column
deba@459: 
deba@459:     /// Return the primal value of the column.
deba@459:     /// \pre The problem is solved.
deba@459:     Value primal(Col c) const { return _getPrimal(cols(id(c))); }
deba@459: 
deba@459:     /// Return the primal value of the expression
deba@459: 
deba@459:     /// Return the primal value of the expression, i.e. the dot
deba@459:     /// product of the primal solution and the expression.
deba@459:     /// \pre The problem is solved.
deba@459:     Value primal(const Expr& e) const {
deba@459:       double res = *e;
deba@459:       for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
deba@459:         res += *c * primal(c);
deba@459:       }
deba@459:       return res;
deba@458:     }
deba@459:     /// Returns a component of the primal ray
deba@459:     
deba@459:     /// The primal ray is solution of the modified primal problem,
deba@459:     /// where we change each finite bound to 0, and we looking for a
deba@459:     /// negative objective value in case of minimization, and positive
deba@459:     /// objective value for maximization. If there is such solution,
deba@459:     /// that proofs the unsolvability of the dual problem, and if a
deba@459:     /// feasible primal solution exists, then the unboundness of
deba@459:     /// primal problem.
deba@459:     ///
deba@459:     /// \pre The problem is solved and the dual problem is infeasible.
deba@459:     /// \note Some solvers does not provide primal ray calculation
deba@459:     /// functions.
deba@459:     Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
deba@458: 
deba@459:     /// Return the dual value of the row
deba@459: 
deba@459:     /// Return the dual value of the row.
deba@459:     /// \pre The problem is solved.
deba@459:     Value dual(Row r) const { return _getDual(rows(id(r))); }
deba@459: 
deba@459:     /// Return the dual value of the dual expression
deba@459: 
deba@459:     /// Return the dual value of the dual expression, i.e. the dot
deba@459:     /// product of the dual solution and the dual expression.
deba@459:     /// \pre The problem is solved.
deba@459:     Value dual(const DualExpr& e) const {
deba@459:       double res = 0.0;
deba@459:       for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
deba@459:         res += *r * dual(r);
deba@458:       }
deba@458:       return res;
deba@458:     }
deba@458: 
deba@459:     /// Returns a component of the dual ray
deba@459:     
deba@459:     /// The dual ray is solution of the modified primal problem, where
deba@459:     /// we change each finite bound to 0 (i.e. the objective function
deba@459:     /// coefficients in the primal problem), and we looking for a
deba@459:     /// ositive objective value. If there is such solution, that
deba@459:     /// proofs the unsolvability of the primal problem, and if a
deba@459:     /// feasible dual solution exists, then the unboundness of
deba@459:     /// dual problem.
deba@459:     ///
deba@459:     /// \pre The problem is solved and the primal problem is infeasible.
deba@459:     /// \note Some solvers does not provide dual ray calculation
deba@459:     /// functions.
deba@459:     Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
deba@458: 
deba@459:     /// Return the basis status of the column
deba@458: 
deba@459:     /// \see VarStatus
deba@459:     VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
deba@459: 
deba@459:     /// Return the basis status of the row
deba@459: 
deba@459:     /// \see VarStatus
deba@459:     VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
deba@459: 
deba@459:     ///The value of the objective function
deba@458: 
deba@458:     ///\return
deba@458:     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
deba@458:     /// of the primal problem, depending on whether we minimize or maximize.
deba@458:     ///- \ref NaN if no primal solution is found.
deba@458:     ///- The (finite) objective value if an optimal solution is found.
deba@459:     Value primal() const { return _getPrimalValue()+obj_const_comp;}
deba@458:     ///@}
deba@458: 
deba@459:   protected:
deba@459: 
deba@458:   };
deba@458: 
deba@458: 
deba@458:   /// \ingroup lp_group
deba@458:   ///
deba@458:   /// \brief Common base class for MIP solvers
deba@459:   ///
deba@459:   /// This class is an abstract base class for MIP solvers. This class
deba@459:   /// provides a full interface for set and modify an MIP problem,
deba@459:   /// solve it and retrieve the solution. You can use one of the
deba@459:   /// descendants as a concrete implementation, or the \c Lp
deba@459:   /// default MIP solver. However, if you would like to handle MIP
deba@459:   /// solvers as reference or pointer in a generic way, you can use
deba@459:   /// this class directly.
deba@459:   class MipSolver : virtual public LpBase {
deba@458:   public:
deba@458: 
deba@459:     /// The problem types for MIP problems
deba@459:     enum ProblemType {
kpeter@576:       /// = 0. Feasible solution hasn't been found (but may exist).
deba@459:       UNDEFINED = 0,
kpeter@576:       /// = 1. The problem has no feasible solution.
deba@459:       INFEASIBLE = 1,
kpeter@576:       /// = 2. Feasible solution found.
deba@459:       FEASIBLE = 2,
kpeter@576:       /// = 3. Optimal solution exists and found.
deba@459:       OPTIMAL = 3,
kpeter@576:       /// = 4. The cost function is unbounded.
kpeter@576:       ///The Mip or at least the relaxed problem is unbounded.
deba@459:       UNBOUNDED = 4
deba@459:     };
deba@459: 
alpar@528:     ///Allocate a new MIP problem instance
alpar@528:     virtual MipSolver* newSolver() const = 0;
alpar@528:     ///Make a copy of the MIP problem
alpar@528:     virtual MipSolver* cloneSolver() const = 0;
alpar@528: 
deba@459:     ///\name Solve the MIP
deba@459: 
deba@459:     ///@{
deba@459: 
deba@459:     /// Solve the MIP problem at hand
deba@459:     ///
deba@459:     ///\return The result of the optimization procedure. Possible
deba@459:     ///values and their meanings can be found in the documentation of
deba@459:     ///\ref SolveExitStatus.
deba@459:     SolveExitStatus solve() { return _solve(); }
deba@459: 
deba@459:     ///@}
deba@459: 
kpeter@576:     ///\name Set Column Type
deba@459:     ///@{
deba@459: 
deba@459:     ///Possible variable (column) types (e.g. real, integer, binary etc.)
deba@458:     enum ColTypes {
kpeter@576:       /// = 0. Continuous variable (default).
deba@458:       REAL = 0,
kpeter@576:       /// = 1. Integer variable.
deba@459:       INTEGER = 1
deba@458:     };
deba@458: 
deba@459:     ///Sets the type of the given column to the given type
deba@459: 
deba@459:     ///Sets the type of the given column to the given type.
deba@458:     ///
deba@458:     void colType(Col c, ColTypes col_type) {
deba@459:       _setColType(cols(id(c)),col_type);
deba@458:     }
deba@458: 
deba@458:     ///Gives back the type of the column.
deba@459: 
deba@459:     ///Gives back the type of the column.
deba@458:     ///
deba@458:     ColTypes colType(Col c) const {
deba@459:       return _getColType(cols(id(c)));
deba@459:     }
deba@459:     ///@}
deba@459: 
kpeter@576:     ///\name Obtain the Solution
deba@459: 
deba@459:     ///@{
deba@459: 
deba@459:     /// The type of the MIP problem
deba@459:     ProblemType type() const {
deba@459:       return _getType();
deba@458:     }
deba@458: 
deba@459:     /// Return the value of the row in the solution
deba@459: 
deba@459:     ///  Return the value of the row in the solution.
deba@459:     /// \pre The problem is solved.
deba@459:     Value sol(Col c) const { return _getSol(cols(id(c))); }
deba@459: 
deba@459:     /// Return the value of the expression in the solution
deba@459: 
deba@459:     /// Return the value of the expression in the solution, i.e. the
deba@459:     /// dot product of the solution and the expression.
deba@459:     /// \pre The problem is solved.
deba@459:     Value sol(const Expr& e) const {
deba@459:       double res = *e;
deba@459:       for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
deba@459:         res += *c * sol(c);
deba@459:       }
deba@459:       return res;
deba@458:     }
deba@459:     ///The value of the objective function
deba@459:     
deba@459:     ///\return
deba@459:     ///- \ref INF or -\ref INF means either infeasibility or unboundedness
deba@459:     /// of the problem, depending on whether we minimize or maximize.
deba@459:     ///- \ref NaN if no primal solution is found.
deba@459:     ///- The (finite) objective value if an optimal solution is found.
deba@459:     Value solValue() const { return _getSolValue()+obj_const_comp;}
deba@459:     ///@}
deba@458: 
deba@458:   protected:
deba@458: 
deba@459:     virtual SolveExitStatus _solve() = 0;
deba@459:     virtual ColTypes _getColType(int col) const = 0;
deba@459:     virtual void _setColType(int col, ColTypes col_type) = 0;
deba@459:     virtual ProblemType _getType() const = 0;
deba@459:     virtual Value _getSol(int i) const = 0;
deba@459:     virtual Value _getSolValue() const = 0;
deba@458: 
deba@458:   };
deba@458: 
deba@458: 
deba@458: 
deba@458: } //namespace lemon
deba@458: 
deba@458: #endif //LEMON_LP_BASE_H