1.1 --- a/lemon/lp_base.h Tue Dec 02 21:40:33 2008 +0100
1.2 +++ b/lemon/lp_base.h Tue Dec 02 22:48:28 2008 +0100
1.3 @@ -25,41 +25,28 @@
1.4 #include<limits>
1.5 #include<lemon/math.h>
1.6
1.7 +#include<lemon/error.h>
1.8 +#include<lemon/assert.h>
1.9 +
1.10 #include<lemon/core.h>
1.11 -#include<lemon/bits/lp_id.h>
1.12 +#include<lemon/bits/solver_bits.h>
1.13
1.14 ///\file
1.15 ///\brief The interface of the LP solver interface.
1.16 ///\ingroup lp_group
1.17 namespace lemon {
1.18
1.19 - /// Function to decide whether a floating point value is finite or not.
1.20 + ///Common base class for LP and MIP solvers
1.21
1.22 - /// Retruns true if the argument is not infinity, minus infinity or NaN.
1.23 - /// It does the same as the isfinite() function defined by C99.
1.24 - template <typename T>
1.25 - bool isFinite(T value)
1.26 - {
1.27 - typedef std::numeric_limits<T> Lim;
1.28 - if ((Lim::has_infinity && (value == Lim::infinity() || value ==
1.29 - -Lim::infinity())) ||
1.30 - ((Lim::has_quiet_NaN || Lim::has_signaling_NaN) && value != value))
1.31 - {
1.32 - return false;
1.33 - }
1.34 - return true;
1.35 - }
1.36 -
1.37 - ///Common base class for LP solvers
1.38 -
1.39 - ///\todo Much more docs
1.40 + ///Usually this class is not used directly, please use one of the concrete
1.41 + ///implementations of the solver interface.
1.42 ///\ingroup lp_group
1.43 - class LpSolverBase {
1.44 + class LpBase {
1.45
1.46 protected:
1.47
1.48 - _lp_bits::LpId rows;
1.49 - _lp_bits::LpId cols;
1.50 + _solver_bits::VarIndex rows;
1.51 + _solver_bits::VarIndex cols;
1.52
1.53 public:
1.54
1.55 @@ -74,38 +61,12 @@
1.56 UNSOLVED = 1
1.57 };
1.58
1.59 - ///\e
1.60 - enum SolutionStatus {
1.61 - ///Feasible solution hasn't been found (but may exist).
1.62 -
1.63 - ///\todo NOTFOUND might be a better name.
1.64 - ///
1.65 - UNDEFINED = 0,
1.66 - ///The problem has no feasible solution
1.67 - INFEASIBLE = 1,
1.68 - ///Feasible solution found
1.69 - FEASIBLE = 2,
1.70 - ///Optimal solution exists and found
1.71 - OPTIMAL = 3,
1.72 - ///The cost function is unbounded
1.73 -
1.74 - ///\todo Give a feasible solution and an infinite ray (and the
1.75 - ///corresponding bases)
1.76 - INFINITE = 4
1.77 - };
1.78 -
1.79 - ///\e The type of the investigated LP problem
1.80 - enum ProblemTypes {
1.81 - ///Primal-dual feasible
1.82 - PRIMAL_DUAL_FEASIBLE = 0,
1.83 - ///Primal feasible dual infeasible
1.84 - PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
1.85 - ///Primal infeasible dual feasible
1.86 - PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
1.87 - ///Primal-dual infeasible
1.88 - PRIMAL_DUAL_INFEASIBLE = 3,
1.89 - ///Could not determine so far
1.90 - UNKNOWN = 4
1.91 + ///Direction of the optimization
1.92 + enum Sense {
1.93 + /// Minimization
1.94 + MIN,
1.95 + /// Maximization
1.96 + MAX
1.97 };
1.98
1.99 ///The floating point type used by the solver
1.100 @@ -115,11 +76,10 @@
1.101 ///The not a number constant
1.102 static const Value NaN;
1.103
1.104 - static inline bool isNaN(const Value& v) { return v!=v; }
1.105 -
1.106 friend class Col;
1.107 friend class ColIt;
1.108 friend class Row;
1.109 + friend class RowIt;
1.110
1.111 ///Refer to a column of the LP.
1.112
1.113 @@ -128,43 +88,93 @@
1.114 ///Its value remains valid and correct even after the addition or erase of
1.115 ///other columns.
1.116 ///
1.117 - ///\todo Document what can one do with a Col (INVALID, comparing,
1.118 - ///it is similar to Node/Edge)
1.119 + ///\note This class is similar to other Item types in LEMON, like
1.120 + ///Node and Arc types in digraph.
1.121 class Col {
1.122 + friend class LpBase;
1.123 protected:
1.124 - int id;
1.125 - friend class LpSolverBase;
1.126 - friend class MipSolverBase;
1.127 - explicit Col(int _id) : id(_id) {}
1.128 + int _id;
1.129 + explicit Col(int id) : _id(id) {}
1.130 public:
1.131 typedef Value ExprValue;
1.132 - typedef True LpSolverCol;
1.133 + typedef True LpCol;
1.134 + /// Default constructor
1.135 +
1.136 + /// \warning The default constructor sets the Col to an
1.137 + /// undefined value.
1.138 Col() {}
1.139 - Col(const Invalid&) : id(-1) {}
1.140 - bool operator< (Col c) const {return id< c.id;}
1.141 - bool operator> (Col c) const {return id> c.id;}
1.142 - bool operator==(Col c) const {return id==c.id;}
1.143 - bool operator!=(Col c) const {return id!=c.id;}
1.144 + /// Invalid constructor \& conversion.
1.145 +
1.146 + /// This constructor initializes the Col to be invalid.
1.147 + /// \sa Invalid for more details.
1.148 + Col(const Invalid&) : _id(-1) {}
1.149 + /// Equality operator
1.150 +
1.151 + /// Two \ref Col "Col"s are equal if and only if they point to
1.152 + /// the same LP column or both are invalid.
1.153 + bool operator==(Col c) const {return _id == c._id;}
1.154 + /// Inequality operator
1.155 +
1.156 + /// \sa operator==(Col c)
1.157 + ///
1.158 + bool operator!=(Col c) const {return _id != c._id;}
1.159 + /// Artificial ordering operator.
1.160 +
1.161 + /// To allow the use of this object in std::map or similar
1.162 + /// associative container we require this.
1.163 + ///
1.164 + /// \note This operator only have to define some strict ordering of
1.165 + /// the items; this order has nothing to do with the iteration
1.166 + /// ordering of the items.
1.167 + bool operator<(Col c) const {return _id < c._id;}
1.168 };
1.169
1.170 + ///Iterator for iterate over the columns of an LP problem
1.171 +
1.172 + /// Its usage is quite simple, for example you can count the number
1.173 + /// of columns in an LP \c lp:
1.174 + ///\code
1.175 + /// int count=0;
1.176 + /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
1.177 + ///\endcode
1.178 class ColIt : public Col {
1.179 - const LpSolverBase *_lp;
1.180 + const LpBase *_solver;
1.181 public:
1.182 + /// Default constructor
1.183 +
1.184 + /// \warning The default constructor sets the iterator
1.185 + /// to an undefined value.
1.186 ColIt() {}
1.187 - ColIt(const LpSolverBase &lp) : _lp(&lp)
1.188 + /// Sets the iterator to the first Col
1.189 +
1.190 + /// Sets the iterator to the first Col.
1.191 + ///
1.192 + ColIt(const LpBase &solver) : _solver(&solver)
1.193 {
1.194 - _lp->cols.firstFix(id);
1.195 + _solver->cols.firstItem(_id);
1.196 }
1.197 + /// Invalid constructor \& conversion
1.198 +
1.199 + /// Initialize the iterator to be invalid.
1.200 + /// \sa Invalid for more details.
1.201 ColIt(const Invalid&) : Col(INVALID) {}
1.202 + /// Next column
1.203 +
1.204 + /// Assign the iterator to the next column.
1.205 + ///
1.206 ColIt &operator++()
1.207 {
1.208 - _lp->cols.nextFix(id);
1.209 + _solver->cols.nextItem(_id);
1.210 return *this;
1.211 }
1.212 };
1.213
1.214 - static int id(const Col& col) { return col.id; }
1.215 -
1.216 + /// \brief Returns the ID of the column.
1.217 + static int id(const Col& col) { return col._id; }
1.218 + /// \brief Returns the column with the given ID.
1.219 + ///
1.220 + /// \pre The argument should be a valid column ID in the LP problem.
1.221 + static Col colFromId(int id) { return Col(id); }
1.222
1.223 ///Refer to a row of the LP.
1.224
1.225 @@ -173,61 +183,93 @@
1.226 ///Its value remains valid and correct even after the addition or erase of
1.227 ///other rows.
1.228 ///
1.229 - ///\todo Document what can one do with a Row (INVALID, comparing,
1.230 - ///it is similar to Node/Edge)
1.231 + ///\note This class is similar to other Item types in LEMON, like
1.232 + ///Node and Arc types in digraph.
1.233 class Row {
1.234 + friend class LpBase;
1.235 protected:
1.236 - int id;
1.237 - friend class LpSolverBase;
1.238 - explicit Row(int _id) : id(_id) {}
1.239 + int _id;
1.240 + explicit Row(int id) : _id(id) {}
1.241 public:
1.242 typedef Value ExprValue;
1.243 - typedef True LpSolverRow;
1.244 + typedef True LpRow;
1.245 + /// Default constructor
1.246 +
1.247 + /// \warning The default constructor sets the Row to an
1.248 + /// undefined value.
1.249 Row() {}
1.250 - Row(const Invalid&) : id(-1) {}
1.251 + /// Invalid constructor \& conversion.
1.252 +
1.253 + /// This constructor initializes the Row to be invalid.
1.254 + /// \sa Invalid for more details.
1.255 + Row(const Invalid&) : _id(-1) {}
1.256 + /// Equality operator
1.257
1.258 - bool operator< (Row c) const {return id< c.id;}
1.259 - bool operator> (Row c) const {return id> c.id;}
1.260 - bool operator==(Row c) const {return id==c.id;}
1.261 - bool operator!=(Row c) const {return id!=c.id;}
1.262 + /// Two \ref Row "Row"s are equal if and only if they point to
1.263 + /// the same LP row or both are invalid.
1.264 + bool operator==(Row r) const {return _id == r._id;}
1.265 + /// Inequality operator
1.266 +
1.267 + /// \sa operator==(Row r)
1.268 + ///
1.269 + bool operator!=(Row r) const {return _id != r._id;}
1.270 + /// Artificial ordering operator.
1.271 +
1.272 + /// To allow the use of this object in std::map or similar
1.273 + /// associative container we require this.
1.274 + ///
1.275 + /// \note This operator only have to define some strict ordering of
1.276 + /// the items; this order has nothing to do with the iteration
1.277 + /// ordering of the items.
1.278 + bool operator<(Row r) const {return _id < r._id;}
1.279 };
1.280
1.281 + ///Iterator for iterate over the rows of an LP problem
1.282 +
1.283 + /// Its usage is quite simple, for example you can count the number
1.284 + /// of rows in an LP \c lp:
1.285 + ///\code
1.286 + /// int count=0;
1.287 + /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
1.288 + ///\endcode
1.289 class RowIt : public Row {
1.290 - const LpSolverBase *_lp;
1.291 + const LpBase *_solver;
1.292 public:
1.293 + /// Default constructor
1.294 +
1.295 + /// \warning The default constructor sets the iterator
1.296 + /// to an undefined value.
1.297 RowIt() {}
1.298 - RowIt(const LpSolverBase &lp) : _lp(&lp)
1.299 + /// Sets the iterator to the first Row
1.300 +
1.301 + /// Sets the iterator to the first Row.
1.302 + ///
1.303 + RowIt(const LpBase &solver) : _solver(&solver)
1.304 {
1.305 - _lp->rows.firstFix(id);
1.306 + _solver->rows.firstItem(_id);
1.307 }
1.308 + /// Invalid constructor \& conversion
1.309 +
1.310 + /// Initialize the iterator to be invalid.
1.311 + /// \sa Invalid for more details.
1.312 RowIt(const Invalid&) : Row(INVALID) {}
1.313 + /// Next row
1.314 +
1.315 + /// Assign the iterator to the next row.
1.316 + ///
1.317 RowIt &operator++()
1.318 {
1.319 - _lp->rows.nextFix(id);
1.320 + _solver->rows.nextItem(_id);
1.321 return *this;
1.322 }
1.323 };
1.324
1.325 - static int id(const Row& row) { return row.id; }
1.326 -
1.327 - protected:
1.328 -
1.329 - int _lpId(const Col& c) const {
1.330 - return cols.floatingId(id(c));
1.331 - }
1.332 -
1.333 - int _lpId(const Row& r) const {
1.334 - return rows.floatingId(id(r));
1.335 - }
1.336 -
1.337 - Col _item(int i, Col) const {
1.338 - return Col(cols.fixId(i));
1.339 - }
1.340 -
1.341 - Row _item(int i, Row) const {
1.342 - return Row(rows.fixId(i));
1.343 - }
1.344 -
1.345 + /// \brief Returns the ID of the row.
1.346 + static int id(const Row& row) { return row._id; }
1.347 + /// \brief Returns the row with the given ID.
1.348 + ///
1.349 + /// \pre The argument should be a valid row ID in the LP problem.
1.350 + static Row rowFromId(int id) { return Row(id); }
1.351
1.352 public:
1.353
1.354 @@ -238,10 +280,6 @@
1.355 ///
1.356 ///There are several ways to access and modify the contents of this
1.357 ///container.
1.358 - ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
1.359 - ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can
1.360 - ///read and modify the coefficients like
1.361 - ///these.
1.362 ///\code
1.363 ///e[v]=5;
1.364 ///e[v]+=12;
1.365 @@ -250,10 +288,10 @@
1.366 ///or you can also iterate through its elements.
1.367 ///\code
1.368 ///double s=0;
1.369 - ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
1.370 - /// s+=i->second;
1.371 + ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
1.372 + /// s+=*i * primal(i);
1.373 ///\endcode
1.374 - ///(This code computes the sum of all coefficients).
1.375 + ///(This code computes the primal value of the expression).
1.376 ///- Numbers (<tt>double</tt>'s)
1.377 ///and variables (\ref Col "Col"s) directly convert to an
1.378 ///\ref Expr and the usual linear operations are defined, so
1.379 @@ -262,7 +300,7 @@
1.380 ///2*v-3.12*(v-w/2)+2
1.381 ///v*2.1+(3*v+(v*12+w+6)*3)/2
1.382 ///\endcode
1.383 - ///are valid \ref Expr "Expr"essions.
1.384 + ///are valid expressions.
1.385 ///The usual assignment operations are also defined.
1.386 ///\code
1.387 ///e=v+w;
1.388 @@ -270,105 +308,218 @@
1.389 ///e*=3.4;
1.390 ///e/=5;
1.391 ///\endcode
1.392 - ///- The constant member can be set and read by \ref constComp()
1.393 + ///- The constant member can be set and read by dereference
1.394 + /// operator (unary *)
1.395 + ///
1.396 ///\code
1.397 - ///e.constComp()=12;
1.398 - ///double c=e.constComp();
1.399 + ///*e=12;
1.400 + ///double c=*e;
1.401 ///\endcode
1.402 ///
1.403 - ///\note \ref clear() not only sets all coefficients to 0 but also
1.404 - ///clears the constant components.
1.405 - ///
1.406 ///\sa Constr
1.407 - ///
1.408 - class Expr : public std::map<Col,Value>
1.409 - {
1.410 + class Expr {
1.411 + friend class LpBase;
1.412 public:
1.413 - typedef LpSolverBase::Col Key;
1.414 - typedef LpSolverBase::Value Value;
1.415 + /// The key type of the expression
1.416 + typedef LpBase::Col Key;
1.417 + /// The value type of the expression
1.418 + typedef LpBase::Value Value;
1.419
1.420 protected:
1.421 - typedef std::map<Col,Value> Base;
1.422 + Value const_comp;
1.423 + std::map<int, Value> comps;
1.424
1.425 - Value const_comp;
1.426 public:
1.427 - typedef True IsLinExpression;
1.428 - ///\e
1.429 - Expr() : Base(), const_comp(0) { }
1.430 - ///\e
1.431 - Expr(const Key &v) : const_comp(0) {
1.432 - Base::insert(std::make_pair(v, 1));
1.433 + typedef True SolverExpr;
1.434 + /// Default constructor
1.435 +
1.436 + /// Construct an empty expression, the coefficients and
1.437 + /// the constant component are initialized to zero.
1.438 + Expr() : const_comp(0) {}
1.439 + /// Construct an expression from a column
1.440 +
1.441 + /// Construct an expression, which has a term with \c c variable
1.442 + /// and 1.0 coefficient.
1.443 + Expr(const Col &c) : const_comp(0) {
1.444 + typedef std::map<int, Value>::value_type pair_type;
1.445 + comps.insert(pair_type(id(c), 1));
1.446 }
1.447 - ///\e
1.448 + /// Construct an expression from a constant
1.449 +
1.450 + /// Construct an expression, which's constant component is \c v.
1.451 + ///
1.452 Expr(const Value &v) : const_comp(v) {}
1.453 - ///\e
1.454 - void set(const Key &v,const Value &c) {
1.455 - Base::insert(std::make_pair(v, c));
1.456 - }
1.457 - ///\e
1.458 - Value &constComp() { return const_comp; }
1.459 - ///\e
1.460 - const Value &constComp() const { return const_comp; }
1.461 -
1.462 - ///Removes the components with zero coefficient.
1.463 - void simplify() {
1.464 - for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.465 - Base::iterator j=i;
1.466 - ++j;
1.467 - if ((*i).second==0) Base::erase(i);
1.468 - i=j;
1.469 + /// Returns the coefficient of the column
1.470 + Value operator[](const Col& c) const {
1.471 + std::map<int, Value>::const_iterator it=comps.find(id(c));
1.472 + if (it != comps.end()) {
1.473 + return it->second;
1.474 + } else {
1.475 + return 0;
1.476 }
1.477 }
1.478 -
1.479 - void simplify() const {
1.480 - const_cast<Expr*>(this)->simplify();
1.481 + /// Returns the coefficient of the column
1.482 + Value& operator[](const Col& c) {
1.483 + return comps[id(c)];
1.484 + }
1.485 + /// Sets the coefficient of the column
1.486 + void set(const Col &c, const Value &v) {
1.487 + if (v != 0.0) {
1.488 + typedef std::map<int, Value>::value_type pair_type;
1.489 + comps.insert(pair_type(id(c), v));
1.490 + } else {
1.491 + comps.erase(id(c));
1.492 + }
1.493 + }
1.494 + /// Returns the constant component of the expression
1.495 + Value& operator*() { return const_comp; }
1.496 + /// Returns the constant component of the expression
1.497 + const Value& operator*() const { return const_comp; }
1.498 + /// \brief Removes the coefficients which's absolute value does
1.499 + /// not exceed \c epsilon. It also sets to zero the constant
1.500 + /// component, if it does not exceed epsilon in absolute value.
1.501 + void simplify(Value epsilon = 0.0) {
1.502 + std::map<int, Value>::iterator it=comps.begin();
1.503 + while (it != comps.end()) {
1.504 + std::map<int, Value>::iterator jt=it;
1.505 + ++jt;
1.506 + if (std::fabs((*it).second) <= epsilon) comps.erase(it);
1.507 + it=jt;
1.508 + }
1.509 + if (std::fabs(const_comp) <= epsilon) const_comp = 0;
1.510 }
1.511
1.512 - ///Removes the coefficients closer to zero than \c tolerance.
1.513 - void simplify(double &tolerance) {
1.514 - for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.515 - Base::iterator j=i;
1.516 - ++j;
1.517 - if (std::fabs((*i).second)<tolerance) Base::erase(i);
1.518 - i=j;
1.519 - }
1.520 + void simplify(Value epsilon = 0.0) const {
1.521 + const_cast<Expr*>(this)->simplify(epsilon);
1.522 }
1.523
1.524 ///Sets all coefficients and the constant component to 0.
1.525 void clear() {
1.526 - Base::clear();
1.527 + comps.clear();
1.528 const_comp=0;
1.529 }
1.530
1.531 - ///\e
1.532 + ///Compound assignment
1.533 Expr &operator+=(const Expr &e) {
1.534 - for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.535 - (*this)[j->first]+=j->second;
1.536 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.537 + it!=e.comps.end(); ++it)
1.538 + comps[it->first]+=it->second;
1.539 const_comp+=e.const_comp;
1.540 return *this;
1.541 }
1.542 - ///\e
1.543 + ///Compound assignment
1.544 Expr &operator-=(const Expr &e) {
1.545 - for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.546 - (*this)[j->first]-=j->second;
1.547 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.548 + it!=e.comps.end(); ++it)
1.549 + comps[it->first]-=it->second;
1.550 const_comp-=e.const_comp;
1.551 return *this;
1.552 }
1.553 - ///\e
1.554 - Expr &operator*=(const Value &c) {
1.555 - for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.556 - j->second*=c;
1.557 - const_comp*=c;
1.558 + ///Multiply with a constant
1.559 + Expr &operator*=(const Value &v) {
1.560 + for (std::map<int, Value>::iterator it=comps.begin();
1.561 + it!=comps.end(); ++it)
1.562 + it->second*=v;
1.563 + const_comp*=v;
1.564 return *this;
1.565 }
1.566 - ///\e
1.567 + ///Division with a constant
1.568 Expr &operator/=(const Value &c) {
1.569 - for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.570 - j->second/=c;
1.571 + for (std::map<int, Value>::iterator it=comps.begin();
1.572 + it!=comps.end(); ++it)
1.573 + it->second/=c;
1.574 const_comp/=c;
1.575 return *this;
1.576 }
1.577
1.578 + ///Iterator over the expression
1.579 +
1.580 + ///The iterator iterates over the terms of the expression.
1.581 + ///
1.582 + ///\code
1.583 + ///double s=0;
1.584 + ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
1.585 + /// s+= *i * primal(i);
1.586 + ///\endcode
1.587 + class CoeffIt {
1.588 + private:
1.589 +
1.590 + std::map<int, Value>::iterator _it, _end;
1.591 +
1.592 + public:
1.593 +
1.594 + /// Sets the iterator to the first term
1.595 +
1.596 + /// Sets the iterator to the first term of the expression.
1.597 + ///
1.598 + CoeffIt(Expr& e)
1.599 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.600 +
1.601 + /// Convert the iterator to the column of the term
1.602 + operator Col() const {
1.603 + return colFromId(_it->first);
1.604 + }
1.605 +
1.606 + /// Returns the coefficient of the term
1.607 + Value& operator*() { return _it->second; }
1.608 +
1.609 + /// Returns the coefficient of the term
1.610 + const Value& operator*() const { return _it->second; }
1.611 + /// Next term
1.612 +
1.613 + /// Assign the iterator to the next term.
1.614 + ///
1.615 + CoeffIt& operator++() { ++_it; return *this; }
1.616 +
1.617 + /// Equality operator
1.618 + bool operator==(Invalid) const { return _it == _end; }
1.619 + /// Inequality operator
1.620 + bool operator!=(Invalid) const { return _it != _end; }
1.621 + };
1.622 +
1.623 + /// Const iterator over the expression
1.624 +
1.625 + ///The iterator iterates over the terms of the expression.
1.626 + ///
1.627 + ///\code
1.628 + ///double s=0;
1.629 + ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
1.630 + /// s+=*i * primal(i);
1.631 + ///\endcode
1.632 + class ConstCoeffIt {
1.633 + private:
1.634 +
1.635 + std::map<int, Value>::const_iterator _it, _end;
1.636 +
1.637 + public:
1.638 +
1.639 + /// Sets the iterator to the first term
1.640 +
1.641 + /// Sets the iterator to the first term of the expression.
1.642 + ///
1.643 + ConstCoeffIt(const Expr& e)
1.644 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.645 +
1.646 + /// Convert the iterator to the column of the term
1.647 + operator Col() const {
1.648 + return colFromId(_it->first);
1.649 + }
1.650 +
1.651 + /// Returns the coefficient of the term
1.652 + const Value& operator*() const { return _it->second; }
1.653 +
1.654 + /// Next term
1.655 +
1.656 + /// Assign the iterator to the next term.
1.657 + ///
1.658 + ConstCoeffIt& operator++() { ++_it; return *this; }
1.659 +
1.660 + /// Equality operator
1.661 + bool operator==(Invalid) const { return _it == _end; }
1.662 + /// Inequality operator
1.663 + bool operator!=(Invalid) const { return _it != _end; }
1.664 + };
1.665 +
1.666 };
1.667
1.668 ///Linear constraint
1.669 @@ -394,13 +545,13 @@
1.670 /// s<=e<=t
1.671 /// e>=t
1.672 ///\endcode
1.673 - ///\warning The validity of a constraint is checked only at run time, so
1.674 - ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw
1.675 - ///an assertion.
1.676 + ///\warning The validity of a constraint is checked only at run
1.677 + ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
1.678 + ///compile, but will fail an assertion.
1.679 class Constr
1.680 {
1.681 public:
1.682 - typedef LpSolverBase::Expr Expr;
1.683 + typedef LpBase::Expr Expr;
1.684 typedef Expr::Key Key;
1.685 typedef Expr::Value Value;
1.686
1.687 @@ -411,15 +562,8 @@
1.688 ///\e
1.689 Constr() : _expr(), _lb(NaN), _ub(NaN) {}
1.690 ///\e
1.691 - Constr(Value lb,const Expr &e,Value ub) :
1.692 + Constr(Value lb, const Expr &e, Value ub) :
1.693 _expr(e), _lb(lb), _ub(ub) {}
1.694 - ///\e
1.695 - Constr(const Expr &e,Value ub) :
1.696 - _expr(e), _lb(NaN), _ub(ub) {}
1.697 - ///\e
1.698 - Constr(Value lb,const Expr &e) :
1.699 - _expr(e), _lb(lb), _ub(NaN) {}
1.700 - ///\e
1.701 Constr(const Expr &e) :
1.702 _expr(e), _lb(NaN), _ub(NaN) {}
1.703 ///\e
1.704 @@ -453,11 +597,11 @@
1.705 const Value &upperBound() const { return _ub; }
1.706 ///Is the constraint lower bounded?
1.707 bool lowerBounded() const {
1.708 - return isFinite(_lb);
1.709 + return _lb != -INF && !std::isnan(_lb);
1.710 }
1.711 ///Is the constraint upper bounded?
1.712 bool upperBounded() const {
1.713 - return isFinite(_ub);
1.714 + return _ub != INF && !std::isnan(_ub);
1.715 }
1.716
1.717 };
1.718 @@ -470,11 +614,6 @@
1.719 ///
1.720 ///There are several ways to access and modify the contents of this
1.721 ///container.
1.722 - ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
1.723 - ///if \c e is an DualExpr and \c v
1.724 - ///and \c w are of type \ref Row, then you can
1.725 - ///read and modify the coefficients like
1.726 - ///these.
1.727 ///\code
1.728 ///e[v]=5;
1.729 ///e[v]+=12;
1.730 @@ -483,8 +622,8 @@
1.731 ///or you can also iterate through its elements.
1.732 ///\code
1.733 ///double s=0;
1.734 - ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
1.735 - /// s+=i->second;
1.736 + ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
1.737 + /// s+=*i;
1.738 ///\endcode
1.739 ///(This code computes the sum of all coefficients).
1.740 ///- Numbers (<tt>double</tt>'s)
1.741 @@ -495,7 +634,7 @@
1.742 ///2*v-3.12*(v-w/2)
1.743 ///v*2.1+(3*v+(v*12+w)*3)/2
1.744 ///\endcode
1.745 - ///are valid \ref DualExpr "DualExpr"essions.
1.746 + ///are valid \ref DualExpr dual expressions.
1.747 ///The usual assignment operations are also defined.
1.748 ///\code
1.749 ///e=v+w;
1.750 @@ -505,146 +644,237 @@
1.751 ///\endcode
1.752 ///
1.753 ///\sa Expr
1.754 - ///
1.755 - class DualExpr : public std::map<Row,Value>
1.756 - {
1.757 + class DualExpr {
1.758 + friend class LpBase;
1.759 public:
1.760 - typedef LpSolverBase::Row Key;
1.761 - typedef LpSolverBase::Value Value;
1.762 + /// The key type of the expression
1.763 + typedef LpBase::Row Key;
1.764 + /// The value type of the expression
1.765 + typedef LpBase::Value Value;
1.766
1.767 protected:
1.768 - typedef std::map<Row,Value> Base;
1.769 + std::map<int, Value> comps;
1.770
1.771 public:
1.772 - typedef True IsLinExpression;
1.773 - ///\e
1.774 - DualExpr() : Base() { }
1.775 - ///\e
1.776 - DualExpr(const Key &v) {
1.777 - Base::insert(std::make_pair(v, 1));
1.778 + typedef True SolverExpr;
1.779 + /// Default constructor
1.780 +
1.781 + /// Construct an empty expression, the coefficients are
1.782 + /// initialized to zero.
1.783 + DualExpr() {}
1.784 + /// Construct an expression from a row
1.785 +
1.786 + /// Construct an expression, which has a term with \c r dual
1.787 + /// variable and 1.0 coefficient.
1.788 + DualExpr(const Row &r) {
1.789 + typedef std::map<int, Value>::value_type pair_type;
1.790 + comps.insert(pair_type(id(r), 1));
1.791 }
1.792 - ///\e
1.793 - void set(const Key &v,const Value &c) {
1.794 - Base::insert(std::make_pair(v, c));
1.795 + /// Returns the coefficient of the row
1.796 + Value operator[](const Row& r) const {
1.797 + std::map<int, Value>::const_iterator it = comps.find(id(r));
1.798 + if (it != comps.end()) {
1.799 + return it->second;
1.800 + } else {
1.801 + return 0;
1.802 + }
1.803 }
1.804 -
1.805 - ///Removes the components with zero coefficient.
1.806 - void simplify() {
1.807 - for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.808 - Base::iterator j=i;
1.809 - ++j;
1.810 - if ((*i).second==0) Base::erase(i);
1.811 - i=j;
1.812 + /// Returns the coefficient of the row
1.813 + Value& operator[](const Row& r) {
1.814 + return comps[id(r)];
1.815 + }
1.816 + /// Sets the coefficient of the row
1.817 + void set(const Row &r, const Value &v) {
1.818 + if (v != 0.0) {
1.819 + typedef std::map<int, Value>::value_type pair_type;
1.820 + comps.insert(pair_type(id(r), v));
1.821 + } else {
1.822 + comps.erase(id(r));
1.823 + }
1.824 + }
1.825 + /// \brief Removes the coefficients which's absolute value does
1.826 + /// not exceed \c epsilon.
1.827 + void simplify(Value epsilon = 0.0) {
1.828 + std::map<int, Value>::iterator it=comps.begin();
1.829 + while (it != comps.end()) {
1.830 + std::map<int, Value>::iterator jt=it;
1.831 + ++jt;
1.832 + if (std::fabs((*it).second) <= epsilon) comps.erase(it);
1.833 + it=jt;
1.834 }
1.835 }
1.836
1.837 - void simplify() const {
1.838 - const_cast<DualExpr*>(this)->simplify();
1.839 - }
1.840 -
1.841 - ///Removes the coefficients closer to zero than \c tolerance.
1.842 - void simplify(double &tolerance) {
1.843 - for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.844 - Base::iterator j=i;
1.845 - ++j;
1.846 - if (std::fabs((*i).second)<tolerance) Base::erase(i);
1.847 - i=j;
1.848 - }
1.849 + void simplify(Value epsilon = 0.0) const {
1.850 + const_cast<DualExpr*>(this)->simplify(epsilon);
1.851 }
1.852
1.853 ///Sets all coefficients to 0.
1.854 void clear() {
1.855 - Base::clear();
1.856 + comps.clear();
1.857 + }
1.858 + ///Compound assignment
1.859 + DualExpr &operator+=(const DualExpr &e) {
1.860 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.861 + it!=e.comps.end(); ++it)
1.862 + comps[it->first]+=it->second;
1.863 + return *this;
1.864 + }
1.865 + ///Compound assignment
1.866 + DualExpr &operator-=(const DualExpr &e) {
1.867 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.868 + it!=e.comps.end(); ++it)
1.869 + comps[it->first]-=it->second;
1.870 + return *this;
1.871 + }
1.872 + ///Multiply with a constant
1.873 + DualExpr &operator*=(const Value &v) {
1.874 + for (std::map<int, Value>::iterator it=comps.begin();
1.875 + it!=comps.end(); ++it)
1.876 + it->second*=v;
1.877 + return *this;
1.878 + }
1.879 + ///Division with a constant
1.880 + DualExpr &operator/=(const Value &v) {
1.881 + for (std::map<int, Value>::iterator it=comps.begin();
1.882 + it!=comps.end(); ++it)
1.883 + it->second/=v;
1.884 + return *this;
1.885 }
1.886
1.887 - ///\e
1.888 - DualExpr &operator+=(const DualExpr &e) {
1.889 - for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.890 - (*this)[j->first]+=j->second;
1.891 - return *this;
1.892 - }
1.893 - ///\e
1.894 - DualExpr &operator-=(const DualExpr &e) {
1.895 - for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.896 - (*this)[j->first]-=j->second;
1.897 - return *this;
1.898 - }
1.899 - ///\e
1.900 - DualExpr &operator*=(const Value &c) {
1.901 - for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.902 - j->second*=c;
1.903 - return *this;
1.904 - }
1.905 - ///\e
1.906 - DualExpr &operator/=(const Value &c) {
1.907 - for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.908 - j->second/=c;
1.909 - return *this;
1.910 - }
1.911 + ///Iterator over the expression
1.912 +
1.913 + ///The iterator iterates over the terms of the expression.
1.914 + ///
1.915 + ///\code
1.916 + ///double s=0;
1.917 + ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
1.918 + /// s+= *i * dual(i);
1.919 + ///\endcode
1.920 + class CoeffIt {
1.921 + private:
1.922 +
1.923 + std::map<int, Value>::iterator _it, _end;
1.924 +
1.925 + public:
1.926 +
1.927 + /// Sets the iterator to the first term
1.928 +
1.929 + /// Sets the iterator to the first term of the expression.
1.930 + ///
1.931 + CoeffIt(DualExpr& e)
1.932 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.933 +
1.934 + /// Convert the iterator to the row of the term
1.935 + operator Row() const {
1.936 + return rowFromId(_it->first);
1.937 + }
1.938 +
1.939 + /// Returns the coefficient of the term
1.940 + Value& operator*() { return _it->second; }
1.941 +
1.942 + /// Returns the coefficient of the term
1.943 + const Value& operator*() const { return _it->second; }
1.944 +
1.945 + /// Next term
1.946 +
1.947 + /// Assign the iterator to the next term.
1.948 + ///
1.949 + CoeffIt& operator++() { ++_it; return *this; }
1.950 +
1.951 + /// Equality operator
1.952 + bool operator==(Invalid) const { return _it == _end; }
1.953 + /// Inequality operator
1.954 + bool operator!=(Invalid) const { return _it != _end; }
1.955 + };
1.956 +
1.957 + ///Iterator over the expression
1.958 +
1.959 + ///The iterator iterates over the terms of the expression.
1.960 + ///
1.961 + ///\code
1.962 + ///double s=0;
1.963 + ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
1.964 + /// s+= *i * dual(i);
1.965 + ///\endcode
1.966 + class ConstCoeffIt {
1.967 + private:
1.968 +
1.969 + std::map<int, Value>::const_iterator _it, _end;
1.970 +
1.971 + public:
1.972 +
1.973 + /// Sets the iterator to the first term
1.974 +
1.975 + /// Sets the iterator to the first term of the expression.
1.976 + ///
1.977 + ConstCoeffIt(const DualExpr& e)
1.978 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.979 +
1.980 + /// Convert the iterator to the row of the term
1.981 + operator Row() const {
1.982 + return rowFromId(_it->first);
1.983 + }
1.984 +
1.985 + /// Returns the coefficient of the term
1.986 + const Value& operator*() const { return _it->second; }
1.987 +
1.988 + /// Next term
1.989 +
1.990 + /// Assign the iterator to the next term.
1.991 + ///
1.992 + ConstCoeffIt& operator++() { ++_it; return *this; }
1.993 +
1.994 + /// Equality operator
1.995 + bool operator==(Invalid) const { return _it == _end; }
1.996 + /// Inequality operator
1.997 + bool operator!=(Invalid) const { return _it != _end; }
1.998 + };
1.999 };
1.1000
1.1001
1.1002 - private:
1.1003 + protected:
1.1004
1.1005 - template <typename _Expr>
1.1006 - class MappedOutputIterator {
1.1007 + class InsertIterator {
1.1008 + private:
1.1009 +
1.1010 + std::map<int, Value>& _host;
1.1011 + const _solver_bits::VarIndex& _index;
1.1012 +
1.1013 public:
1.1014
1.1015 - typedef std::insert_iterator<_Expr> Base;
1.1016 -
1.1017 typedef std::output_iterator_tag iterator_category;
1.1018 typedef void difference_type;
1.1019 typedef void value_type;
1.1020 typedef void reference;
1.1021 typedef void pointer;
1.1022
1.1023 - MappedOutputIterator(const Base& _base, const LpSolverBase& _lp)
1.1024 - : base(_base), lp(_lp) {}
1.1025 + InsertIterator(std::map<int, Value>& host,
1.1026 + const _solver_bits::VarIndex& index)
1.1027 + : _host(host), _index(index) {}
1.1028
1.1029 - MappedOutputIterator& operator*() {
1.1030 + InsertIterator& operator=(const std::pair<int, Value>& value) {
1.1031 + typedef std::map<int, Value>::value_type pair_type;
1.1032 + _host.insert(pair_type(_index[value.first], value.second));
1.1033 return *this;
1.1034 }
1.1035
1.1036 - MappedOutputIterator& operator=(const std::pair<int, Value>& value) {
1.1037 - *base = std::make_pair(lp._item(value.first, typename _Expr::Key()),
1.1038 - value.second);
1.1039 - return *this;
1.1040 - }
1.1041 + InsertIterator& operator*() { return *this; }
1.1042 + InsertIterator& operator++() { return *this; }
1.1043 + InsertIterator operator++(int) { return *this; }
1.1044
1.1045 - MappedOutputIterator& operator++() {
1.1046 - ++base;
1.1047 - return *this;
1.1048 - }
1.1049 -
1.1050 - MappedOutputIterator operator++(int) {
1.1051 - MappedOutputIterator tmp(*this);
1.1052 - ++base;
1.1053 - return tmp;
1.1054 - }
1.1055 -
1.1056 - bool operator==(const MappedOutputIterator& it) const {
1.1057 - return base == it.base;
1.1058 - }
1.1059 -
1.1060 - bool operator!=(const MappedOutputIterator& it) const {
1.1061 - return base != it.base;
1.1062 - }
1.1063 -
1.1064 - private:
1.1065 - Base base;
1.1066 - const LpSolverBase& lp;
1.1067 };
1.1068
1.1069 - template <typename Expr>
1.1070 - class MappedInputIterator {
1.1071 + class ExprIterator {
1.1072 + private:
1.1073 + std::map<int, Value>::const_iterator _host_it;
1.1074 + const _solver_bits::VarIndex& _index;
1.1075 public:
1.1076
1.1077 - typedef typename Expr::const_iterator Base;
1.1078 -
1.1079 - typedef typename Base::iterator_category iterator_category;
1.1080 - typedef typename Base::difference_type difference_type;
1.1081 + typedef std::bidirectional_iterator_tag iterator_category;
1.1082 + typedef std::ptrdiff_t difference_type;
1.1083 typedef const std::pair<int, Value> value_type;
1.1084 typedef value_type reference;
1.1085 +
1.1086 class pointer {
1.1087 public:
1.1088 pointer(value_type& _value) : value(_value) {}
1.1089 @@ -653,82 +883,49 @@
1.1090 value_type value;
1.1091 };
1.1092
1.1093 - MappedInputIterator(const Base& _base, const LpSolverBase& _lp)
1.1094 - : base(_base), lp(_lp) {}
1.1095 + ExprIterator(const std::map<int, Value>::const_iterator& host_it,
1.1096 + const _solver_bits::VarIndex& index)
1.1097 + : _host_it(host_it), _index(index) {}
1.1098
1.1099 reference operator*() {
1.1100 - return std::make_pair(lp._lpId(base->first), base->second);
1.1101 + return std::make_pair(_index(_host_it->first), _host_it->second);
1.1102 }
1.1103
1.1104 pointer operator->() {
1.1105 return pointer(operator*());
1.1106 }
1.1107
1.1108 - MappedInputIterator& operator++() {
1.1109 - ++base;
1.1110 - return *this;
1.1111 + ExprIterator& operator++() { ++_host_it; return *this; }
1.1112 + ExprIterator operator++(int) {
1.1113 + ExprIterator tmp(*this); ++_host_it; return tmp;
1.1114 }
1.1115
1.1116 - MappedInputIterator operator++(int) {
1.1117 - MappedInputIterator tmp(*this);
1.1118 - ++base;
1.1119 - return tmp;
1.1120 + ExprIterator& operator--() { --_host_it; return *this; }
1.1121 + ExprIterator operator--(int) {
1.1122 + ExprIterator tmp(*this); --_host_it; return tmp;
1.1123 }
1.1124
1.1125 - bool operator==(const MappedInputIterator& it) const {
1.1126 - return base == it.base;
1.1127 + bool operator==(const ExprIterator& it) const {
1.1128 + return _host_it == it._host_it;
1.1129 }
1.1130
1.1131 - bool operator!=(const MappedInputIterator& it) const {
1.1132 - return base != it.base;
1.1133 + bool operator!=(const ExprIterator& it) const {
1.1134 + return _host_it != it._host_it;
1.1135 }
1.1136
1.1137 - private:
1.1138 - Base base;
1.1139 - const LpSolverBase& lp;
1.1140 };
1.1141
1.1142 protected:
1.1143
1.1144 - /// STL compatible iterator for lp col
1.1145 - typedef MappedInputIterator<Expr> ConstRowIterator;
1.1146 - /// STL compatible iterator for lp row
1.1147 - typedef MappedInputIterator<DualExpr> ConstColIterator;
1.1148 + //Abstract virtual functions
1.1149 + virtual LpBase* _newSolver() const = 0;
1.1150 + virtual LpBase* _cloneSolver() const = 0;
1.1151
1.1152 - /// STL compatible iterator for lp col
1.1153 - typedef MappedOutputIterator<Expr> RowIterator;
1.1154 - /// STL compatible iterator for lp row
1.1155 - typedef MappedOutputIterator<DualExpr> ColIterator;
1.1156 + virtual int _addColId(int col) { return cols.addIndex(col); }
1.1157 + virtual int _addRowId(int row) { return rows.addIndex(row); }
1.1158
1.1159 - //Abstract virtual functions
1.1160 - virtual LpSolverBase* _newLp() = 0;
1.1161 - virtual LpSolverBase* _copyLp(){
1.1162 - LpSolverBase* newlp = _newLp();
1.1163 -
1.1164 - std::map<Col, Col> ref;
1.1165 - for (LpSolverBase::ColIt it(*this); it != INVALID; ++it) {
1.1166 - Col ccol = newlp->addCol();
1.1167 - ref[it] = ccol;
1.1168 - newlp->colName(ccol, colName(it));
1.1169 - newlp->colLowerBound(ccol, colLowerBound(it));
1.1170 - newlp->colUpperBound(ccol, colUpperBound(it));
1.1171 - }
1.1172 -
1.1173 - for (LpSolverBase::RowIt it(*this); it != INVALID; ++it) {
1.1174 - Expr e = row(it), ce;
1.1175 - for (Expr::iterator jt = e.begin(); jt != e.end(); ++jt) {
1.1176 - ce[ref[jt->first]] = jt->second;
1.1177 - }
1.1178 - ce += e.constComp();
1.1179 - Row r = newlp->addRow(ce);
1.1180 -
1.1181 - double lower, upper;
1.1182 - getRowBounds(it, lower, upper);
1.1183 - newlp->rowBounds(r, lower, upper);
1.1184 - }
1.1185 -
1.1186 - return newlp;
1.1187 - };
1.1188 + virtual void _eraseColId(int col) { cols.eraseIndex(col); }
1.1189 + virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
1.1190
1.1191 virtual int _addCol() = 0;
1.1192 virtual int _addRow() = 0;
1.1193 @@ -736,93 +933,95 @@
1.1194 virtual void _eraseCol(int col) = 0;
1.1195 virtual void _eraseRow(int row) = 0;
1.1196
1.1197 - virtual void _getColName(int col, std::string & name) const = 0;
1.1198 - virtual void _setColName(int col, const std::string & name) = 0;
1.1199 + virtual void _getColName(int col, std::string& name) const = 0;
1.1200 + virtual void _setColName(int col, const std::string& name) = 0;
1.1201 virtual int _colByName(const std::string& name) const = 0;
1.1202
1.1203 - virtual void _setRowCoeffs(int i, ConstRowIterator b,
1.1204 - ConstRowIterator e) = 0;
1.1205 - virtual void _getRowCoeffs(int i, RowIterator b) const = 0;
1.1206 - virtual void _setColCoeffs(int i, ConstColIterator b,
1.1207 - ConstColIterator e) = 0;
1.1208 - virtual void _getColCoeffs(int i, ColIterator b) const = 0;
1.1209 + virtual void _getRowName(int row, std::string& name) const = 0;
1.1210 + virtual void _setRowName(int row, const std::string& name) = 0;
1.1211 + virtual int _rowByName(const std::string& name) const = 0;
1.1212 +
1.1213 + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1.1214 + virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
1.1215 +
1.1216 + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1.1217 + virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
1.1218 +
1.1219 virtual void _setCoeff(int row, int col, Value value) = 0;
1.1220 virtual Value _getCoeff(int row, int col) const = 0;
1.1221 +
1.1222 virtual void _setColLowerBound(int i, Value value) = 0;
1.1223 virtual Value _getColLowerBound(int i) const = 0;
1.1224 +
1.1225 virtual void _setColUpperBound(int i, Value value) = 0;
1.1226 virtual Value _getColUpperBound(int i) const = 0;
1.1227 - virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
1.1228 - virtual void _getRowBounds(int i, Value &lower, Value &upper) const = 0;
1.1229 +
1.1230 + virtual void _setRowLowerBound(int i, Value value) = 0;
1.1231 + virtual Value _getRowLowerBound(int i) const = 0;
1.1232 +
1.1233 + virtual void _setRowUpperBound(int i, Value value) = 0;
1.1234 + virtual Value _getRowUpperBound(int i) const = 0;
1.1235 +
1.1236 + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
1.1237 + virtual void _getObjCoeffs(InsertIterator b) const = 0;
1.1238
1.1239 virtual void _setObjCoeff(int i, Value obj_coef) = 0;
1.1240 virtual Value _getObjCoeff(int i) const = 0;
1.1241 - virtual void _clearObj()=0;
1.1242
1.1243 - virtual SolveExitStatus _solve() = 0;
1.1244 - virtual Value _getPrimal(int i) const = 0;
1.1245 - virtual Value _getDual(int i) const = 0;
1.1246 - virtual Value _getPrimalValue() const = 0;
1.1247 - virtual bool _isBasicCol(int i) const = 0;
1.1248 - virtual SolutionStatus _getPrimalStatus() const = 0;
1.1249 - virtual SolutionStatus _getDualStatus() const = 0;
1.1250 - virtual ProblemTypes _getProblemType() const = 0;
1.1251 + virtual void _setSense(Sense) = 0;
1.1252 + virtual Sense _getSense() const = 0;
1.1253
1.1254 - virtual void _setMax() = 0;
1.1255 - virtual void _setMin() = 0;
1.1256 + virtual void _clear() = 0;
1.1257
1.1258 -
1.1259 - virtual bool _isMax() const = 0;
1.1260 + virtual const char* _solverName() const = 0;
1.1261
1.1262 //Own protected stuff
1.1263
1.1264 //Constant component of the objective function
1.1265 Value obj_const_comp;
1.1266
1.1267 + LpBase() : rows(), cols(), obj_const_comp(0) {}
1.1268 +
1.1269 public:
1.1270
1.1271 - ///\e
1.1272 - LpSolverBase() : obj_const_comp(0) {}
1.1273 -
1.1274 - ///\e
1.1275 - virtual ~LpSolverBase() {}
1.1276 + /// Virtual destructor
1.1277 + virtual ~LpBase() {}
1.1278
1.1279 ///Creates a new LP problem
1.1280 - LpSolverBase* newLp() {return _newLp();}
1.1281 + LpBase* newSolver() {return _newSolver();}
1.1282 ///Makes a copy of the LP problem
1.1283 - LpSolverBase* copyLp() {return _copyLp();}
1.1284 + LpBase* cloneSolver() {return _cloneSolver();}
1.1285 +
1.1286 + ///Gives back the name of the solver.
1.1287 + const char* solverName() const {return _solverName();}
1.1288
1.1289 ///\name Build up and modify the LP
1.1290
1.1291 ///@{
1.1292
1.1293 ///Add a new empty column (i.e a new variable) to the LP
1.1294 - Col addCol() { Col c; _addCol(); c.id = cols.addId(); return c;}
1.1295 + Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
1.1296
1.1297 - ///\brief Adds several new columns
1.1298 - ///(i.e a variables) at once
1.1299 + ///\brief Adds several new columns (i.e variables) at once
1.1300 ///
1.1301 - ///This magic function takes a container as its argument
1.1302 - ///and fills its elements
1.1303 - ///with new columns (i.e. variables)
1.1304 + ///This magic function takes a container as its argument and fills
1.1305 + ///its elements with new columns (i.e. variables)
1.1306 ///\param t can be
1.1307 ///- a standard STL compatible iterable container with
1.1308 - ///\ref Col as its \c values_type
1.1309 - ///like
1.1310 + ///\ref Col as its \c values_type like
1.1311 ///\code
1.1312 - ///std::vector<LpSolverBase::Col>
1.1313 - ///std::list<LpSolverBase::Col>
1.1314 + ///std::vector<LpBase::Col>
1.1315 + ///std::list<LpBase::Col>
1.1316 ///\endcode
1.1317 ///- a standard STL compatible iterable container with
1.1318 - ///\ref Col as its \c mapped_type
1.1319 - ///like
1.1320 + ///\ref Col as its \c mapped_type like
1.1321 ///\code
1.1322 - ///std::map<AnyType,LpSolverBase::Col>
1.1323 + ///std::map<AnyType,LpBase::Col>
1.1324 ///\endcode
1.1325 ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.1326 ///\code
1.1327 - ///ListGraph::NodeMap<LpSolverBase::Col>
1.1328 - ///ListGraph::EdgeMap<LpSolverBase::Col>
1.1329 + ///ListGraph::NodeMap<LpBase::Col>
1.1330 + ///ListGraph::ArcMap<LpBase::Col>
1.1331 ///\endcode
1.1332 ///\return The number of the created column.
1.1333 #ifdef DOXYGEN
1.1334 @@ -830,14 +1029,14 @@
1.1335 int addColSet(T &t) { return 0;}
1.1336 #else
1.1337 template<class T>
1.1338 - typename enable_if<typename T::value_type::LpSolverCol,int>::type
1.1339 + typename enable_if<typename T::value_type::LpCol,int>::type
1.1340 addColSet(T &t,dummy<0> = 0) {
1.1341 int s=0;
1.1342 for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
1.1343 return s;
1.1344 }
1.1345 template<class T>
1.1346 - typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.1347 + typename enable_if<typename T::value_type::second_type::LpCol,
1.1348 int>::type
1.1349 addColSet(T &t,dummy<1> = 1) {
1.1350 int s=0;
1.1351 @@ -848,7 +1047,7 @@
1.1352 return s;
1.1353 }
1.1354 template<class T>
1.1355 - typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.1356 + typename enable_if<typename T::MapIt::Value::LpCol,
1.1357 int>::type
1.1358 addColSet(T &t,dummy<2> = 2) {
1.1359 int s=0;
1.1360 @@ -866,26 +1065,26 @@
1.1361 ///\param c is the column to be modified
1.1362 ///\param e is a dual linear expression (see \ref DualExpr)
1.1363 ///a better one.
1.1364 - void col(Col c,const DualExpr &e) {
1.1365 + void col(Col c, const DualExpr &e) {
1.1366 e.simplify();
1.1367 - _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this),
1.1368 - ConstColIterator(e.end(), *this));
1.1369 + _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), cols),
1.1370 + ExprIterator(e.comps.end(), cols));
1.1371 }
1.1372
1.1373 ///Get a column (i.e a dual constraint) of the LP
1.1374
1.1375 - ///\param r is the column to get
1.1376 + ///\param c is the column to get
1.1377 ///\return the dual expression associated to the column
1.1378 DualExpr col(Col c) const {
1.1379 DualExpr e;
1.1380 - _getColCoeffs(_lpId(c), ColIterator(std::inserter(e, e.end()), *this));
1.1381 + _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
1.1382 return e;
1.1383 }
1.1384
1.1385 ///Add a new column to the LP
1.1386
1.1387 ///\param e is a dual linear expression (see \ref DualExpr)
1.1388 - ///\param obj is the corresponding component of the objective
1.1389 + ///\param o is the corresponding component of the objective
1.1390 ///function. It is 0 by default.
1.1391 ///\return The created column.
1.1392 Col addCol(const DualExpr &e, Value o = 0) {
1.1393 @@ -899,32 +1098,28 @@
1.1394
1.1395 ///This function adds a new empty row (i.e a new constraint) to the LP.
1.1396 ///\return The created row
1.1397 - Row addRow() { Row r; _addRow(); r.id = rows.addId(); return r;}
1.1398 + Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
1.1399
1.1400 - ///\brief Add several new rows
1.1401 - ///(i.e a constraints) at once
1.1402 + ///\brief Add several new rows (i.e constraints) at once
1.1403 ///
1.1404 - ///This magic function takes a container as its argument
1.1405 - ///and fills its elements
1.1406 - ///with new row (i.e. variables)
1.1407 + ///This magic function takes a container as its argument and fills
1.1408 + ///its elements with new row (i.e. variables)
1.1409 ///\param t can be
1.1410 ///- a standard STL compatible iterable container with
1.1411 - ///\ref Row as its \c values_type
1.1412 - ///like
1.1413 + ///\ref Row as its \c values_type like
1.1414 ///\code
1.1415 - ///std::vector<LpSolverBase::Row>
1.1416 - ///std::list<LpSolverBase::Row>
1.1417 + ///std::vector<LpBase::Row>
1.1418 + ///std::list<LpBase::Row>
1.1419 ///\endcode
1.1420 ///- a standard STL compatible iterable container with
1.1421 - ///\ref Row as its \c mapped_type
1.1422 - ///like
1.1423 + ///\ref Row as its \c mapped_type like
1.1424 ///\code
1.1425 - ///std::map<AnyType,LpSolverBase::Row>
1.1426 + ///std::map<AnyType,LpBase::Row>
1.1427 ///\endcode
1.1428 ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.1429 ///\code
1.1430 - ///ListGraph::NodeMap<LpSolverBase::Row>
1.1431 - ///ListGraph::EdgeMap<LpSolverBase::Row>
1.1432 + ///ListGraph::NodeMap<LpBase::Row>
1.1433 + ///ListGraph::ArcMap<LpBase::Row>
1.1434 ///\endcode
1.1435 ///\return The number of rows created.
1.1436 #ifdef DOXYGEN
1.1437 @@ -932,16 +1127,15 @@
1.1438 int addRowSet(T &t) { return 0;}
1.1439 #else
1.1440 template<class T>
1.1441 - typename enable_if<typename T::value_type::LpSolverRow,int>::type
1.1442 - addRowSet(T &t,dummy<0> = 0) {
1.1443 + typename enable_if<typename T::value_type::LpRow,int>::type
1.1444 + addRowSet(T &t, dummy<0> = 0) {
1.1445 int s=0;
1.1446 for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1.1447 return s;
1.1448 }
1.1449 template<class T>
1.1450 - typename enable_if<typename T::value_type::second_type::LpSolverRow,
1.1451 - int>::type
1.1452 - addRowSet(T &t,dummy<1> = 1) {
1.1453 + typename enable_if<typename T::value_type::second_type::LpRow, int>::type
1.1454 + addRowSet(T &t, dummy<1> = 1) {
1.1455 int s=0;
1.1456 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1457 i->second=addRow();
1.1458 @@ -950,9 +1144,8 @@
1.1459 return s;
1.1460 }
1.1461 template<class T>
1.1462 - typename enable_if<typename T::MapIt::Value::LpSolverRow,
1.1463 - int>::type
1.1464 - addRowSet(T &t,dummy<2> = 2) {
1.1465 + typename enable_if<typename T::MapIt::Value::LpRow, int>::type
1.1466 + addRowSet(T &t, dummy<2> = 2) {
1.1467 int s=0;
1.1468 for(typename T::MapIt i(t); i!=INVALID; ++i)
1.1469 {
1.1470 @@ -969,15 +1162,12 @@
1.1471 ///\param l is lower bound (-\ref INF means no bound)
1.1472 ///\param e is a linear expression (see \ref Expr)
1.1473 ///\param u is the upper bound (\ref INF means no bound)
1.1474 - ///\bug This is a temporary function. The interface will change to
1.1475 - ///a better one.
1.1476 - ///\todo Option to control whether a constraint with a single variable is
1.1477 - ///added or not.
1.1478 void row(Row r, Value l, const Expr &e, Value u) {
1.1479 e.simplify();
1.1480 - _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this),
1.1481 - ConstRowIterator(e.end(), *this));
1.1482 - _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp());
1.1483 + _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
1.1484 + ExprIterator(e.comps.end(), cols));
1.1485 + _setRowLowerBound(rows(id(r)),l - *e);
1.1486 + _setRowUpperBound(rows(id(r)),u - *e);
1.1487 }
1.1488
1.1489 ///Set a row (i.e a constraint) of the LP
1.1490 @@ -996,7 +1186,7 @@
1.1491 ///\return the expression associated to the row
1.1492 Expr row(Row r) const {
1.1493 Expr e;
1.1494 - _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this));
1.1495 + _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
1.1496 return e;
1.1497 }
1.1498
1.1499 @@ -1006,8 +1196,6 @@
1.1500 ///\param e is a linear expression (see \ref Expr)
1.1501 ///\param u is the upper bound (\ref INF means no bound)
1.1502 ///\return The created row.
1.1503 - ///\bug This is a temporary function. The interface will change to
1.1504 - ///a better one.
1.1505 Row addRow(Value l,const Expr &e, Value u) {
1.1506 Row r=addRow();
1.1507 row(r,l,e,u);
1.1508 @@ -1023,39 +1211,37 @@
1.1509 row(r,c);
1.1510 return r;
1.1511 }
1.1512 - ///Erase a coloumn (i.e a variable) from the LP
1.1513 + ///Erase a column (i.e a variable) from the LP
1.1514
1.1515 - ///\param c is the coloumn to be deleted
1.1516 - ///\todo Please check this
1.1517 - void eraseCol(Col c) {
1.1518 - _eraseCol(_lpId(c));
1.1519 - cols.eraseId(c.id);
1.1520 + ///\param c is the column to be deleted
1.1521 + void erase(Col c) {
1.1522 + _eraseCol(cols(id(c)));
1.1523 + _eraseColId(cols(id(c)));
1.1524 }
1.1525 - ///Erase a row (i.e a constraint) from the LP
1.1526 + ///Erase a row (i.e a constraint) from the LP
1.1527
1.1528 ///\param r is the row to be deleted
1.1529 - ///\todo Please check this
1.1530 - void eraseRow(Row r) {
1.1531 - _eraseRow(_lpId(r));
1.1532 - rows.eraseId(r.id);
1.1533 + void erase(Row r) {
1.1534 + _eraseRow(rows(id(r)));
1.1535 + _eraseRowId(rows(id(r)));
1.1536 }
1.1537
1.1538 /// Get the name of a column
1.1539
1.1540 - ///\param c is the coresponding coloumn
1.1541 + ///\param c is the coresponding column
1.1542 ///\return The name of the colunm
1.1543 std::string colName(Col c) const {
1.1544 std::string name;
1.1545 - _getColName(_lpId(c), name);
1.1546 + _getColName(cols(id(c)), name);
1.1547 return name;
1.1548 }
1.1549
1.1550 /// Set the name of a column
1.1551
1.1552 - ///\param c is the coresponding coloumn
1.1553 + ///\param c is the coresponding column
1.1554 ///\param name The name to be given
1.1555 void colName(Col c, const std::string& name) {
1.1556 - _setColName(_lpId(c), name);
1.1557 + _setColName(cols(id(c)), name);
1.1558 }
1.1559
1.1560 /// Get the column by its name
1.1561 @@ -1064,27 +1250,52 @@
1.1562 ///\return the proper column or \c INVALID
1.1563 Col colByName(const std::string& name) const {
1.1564 int k = _colByName(name);
1.1565 - return k != -1 ? Col(cols.fixId(k)) : Col(INVALID);
1.1566 + return k != -1 ? Col(cols[k]) : Col(INVALID);
1.1567 + }
1.1568 +
1.1569 + /// Get the name of a row
1.1570 +
1.1571 + ///\param r is the coresponding row
1.1572 + ///\return The name of the row
1.1573 + std::string rowName(Row r) const {
1.1574 + std::string name;
1.1575 + _getRowName(rows(id(r)), name);
1.1576 + return name;
1.1577 + }
1.1578 +
1.1579 + /// Set the name of a row
1.1580 +
1.1581 + ///\param r is the coresponding row
1.1582 + ///\param name The name to be given
1.1583 + void rowName(Row r, const std::string& name) {
1.1584 + _setRowName(rows(id(r)), name);
1.1585 + }
1.1586 +
1.1587 + /// Get the row by its name
1.1588 +
1.1589 + ///\param name The name of the row
1.1590 + ///\return the proper row or \c INVALID
1.1591 + Row rowByName(const std::string& name) const {
1.1592 + int k = _rowByName(name);
1.1593 + return k != -1 ? Row(rows[k]) : Row(INVALID);
1.1594 }
1.1595
1.1596 /// Set an element of the coefficient matrix of the LP
1.1597
1.1598 ///\param r is the row of the element to be modified
1.1599 - ///\param c is the coloumn of the element to be modified
1.1600 + ///\param c is the column of the element to be modified
1.1601 ///\param val is the new value of the coefficient
1.1602 -
1.1603 void coeff(Row r, Col c, Value val) {
1.1604 - _setCoeff(_lpId(r),_lpId(c), val);
1.1605 + _setCoeff(rows(id(r)),cols(id(c)), val);
1.1606 }
1.1607
1.1608 /// Get an element of the coefficient matrix of the LP
1.1609
1.1610 - ///\param r is the row of the element in question
1.1611 - ///\param c is the coloumn of the element in question
1.1612 + ///\param r is the row of the element
1.1613 + ///\param c is the column of the element
1.1614 ///\return the corresponding coefficient
1.1615 -
1.1616 Value coeff(Row r, Col c) const {
1.1617 - return _getCoeff(_lpId(r),_lpId(c));
1.1618 + return _getCoeff(rows(id(r)),cols(id(c)));
1.1619 }
1.1620
1.1621 /// Set the lower bound of a column (i.e a variable)
1.1622 @@ -1093,39 +1304,39 @@
1.1623 /// extended number of type Value, i.e. a finite number of type
1.1624 /// Value or -\ref INF.
1.1625 void colLowerBound(Col c, Value value) {
1.1626 - _setColLowerBound(_lpId(c),value);
1.1627 + _setColLowerBound(cols(id(c)),value);
1.1628 }
1.1629
1.1630 /// Get the lower bound of a column (i.e a variable)
1.1631
1.1632 - /// This function returns the lower bound for column (variable) \t c
1.1633 + /// This function returns the lower bound for column (variable) \c c
1.1634 /// (this might be -\ref INF as well).
1.1635 - ///\return The lower bound for coloumn \t c
1.1636 + ///\return The lower bound for column \c c
1.1637 Value colLowerBound(Col c) const {
1.1638 - return _getColLowerBound(_lpId(c));
1.1639 + return _getColLowerBound(cols(id(c)));
1.1640 }
1.1641
1.1642 ///\brief Set the lower bound of several columns
1.1643 - ///(i.e a variables) at once
1.1644 + ///(i.e variables) at once
1.1645 ///
1.1646 ///This magic function takes a container as its argument
1.1647 ///and applies the function on all of its elements.
1.1648 - /// The lower bound of a variable (column) has to be given by an
1.1649 - /// extended number of type Value, i.e. a finite number of type
1.1650 - /// Value or -\ref INF.
1.1651 + ///The lower bound of a variable (column) has to be given by an
1.1652 + ///extended number of type Value, i.e. a finite number of type
1.1653 + ///Value or -\ref INF.
1.1654 #ifdef DOXYGEN
1.1655 template<class T>
1.1656 void colLowerBound(T &t, Value value) { return 0;}
1.1657 #else
1.1658 template<class T>
1.1659 - typename enable_if<typename T::value_type::LpSolverCol,void>::type
1.1660 + typename enable_if<typename T::value_type::LpCol,void>::type
1.1661 colLowerBound(T &t, Value value,dummy<0> = 0) {
1.1662 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1663 colLowerBound(*i, value);
1.1664 }
1.1665 }
1.1666 template<class T>
1.1667 - typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.1668 + typename enable_if<typename T::value_type::second_type::LpCol,
1.1669 void>::type
1.1670 colLowerBound(T &t, Value value,dummy<1> = 1) {
1.1671 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1672 @@ -1133,7 +1344,7 @@
1.1673 }
1.1674 }
1.1675 template<class T>
1.1676 - typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.1677 + typename enable_if<typename T::MapIt::Value::LpCol,
1.1678 void>::type
1.1679 colLowerBound(T &t, Value value,dummy<2> = 2) {
1.1680 for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1681 @@ -1148,39 +1359,39 @@
1.1682 /// extended number of type Value, i.e. a finite number of type
1.1683 /// Value or \ref INF.
1.1684 void colUpperBound(Col c, Value value) {
1.1685 - _setColUpperBound(_lpId(c),value);
1.1686 + _setColUpperBound(cols(id(c)),value);
1.1687 };
1.1688
1.1689 /// Get the upper bound of a column (i.e a variable)
1.1690
1.1691 - /// This function returns the upper bound for column (variable) \t c
1.1692 + /// This function returns the upper bound for column (variable) \c c
1.1693 /// (this might be \ref INF as well).
1.1694 - ///\return The upper bound for coloumn \t c
1.1695 + /// \return The upper bound for column \c c
1.1696 Value colUpperBound(Col c) const {
1.1697 - return _getColUpperBound(_lpId(c));
1.1698 + return _getColUpperBound(cols(id(c)));
1.1699 }
1.1700
1.1701 ///\brief Set the upper bound of several columns
1.1702 - ///(i.e a variables) at once
1.1703 + ///(i.e variables) at once
1.1704 ///
1.1705 ///This magic function takes a container as its argument
1.1706 ///and applies the function on all of its elements.
1.1707 - /// The upper bound of a variable (column) has to be given by an
1.1708 - /// extended number of type Value, i.e. a finite number of type
1.1709 - /// Value or \ref INF.
1.1710 + ///The upper bound of a variable (column) has to be given by an
1.1711 + ///extended number of type Value, i.e. a finite number of type
1.1712 + ///Value or \ref INF.
1.1713 #ifdef DOXYGEN
1.1714 template<class T>
1.1715 void colUpperBound(T &t, Value value) { return 0;}
1.1716 #else
1.1717 template<class T>
1.1718 - typename enable_if<typename T::value_type::LpSolverCol,void>::type
1.1719 + typename enable_if<typename T::value_type::LpCol,void>::type
1.1720 colUpperBound(T &t, Value value,dummy<0> = 0) {
1.1721 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1722 colUpperBound(*i, value);
1.1723 }
1.1724 }
1.1725 template<class T>
1.1726 - typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.1727 + typename enable_if<typename T::value_type::second_type::LpCol,
1.1728 void>::type
1.1729 colUpperBound(T &t, Value value,dummy<1> = 1) {
1.1730 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1731 @@ -1188,7 +1399,7 @@
1.1732 }
1.1733 }
1.1734 template<class T>
1.1735 - typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.1736 + typename enable_if<typename T::MapIt::Value::LpCol,
1.1737 void>::type
1.1738 colUpperBound(T &t, Value value,dummy<2> = 2) {
1.1739 for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1740 @@ -1204,12 +1415,12 @@
1.1741 /// extended number of type Value, i.e. a finite number of type
1.1742 /// Value, -\ref INF or \ref INF.
1.1743 void colBounds(Col c, Value lower, Value upper) {
1.1744 - _setColLowerBound(_lpId(c),lower);
1.1745 - _setColUpperBound(_lpId(c),upper);
1.1746 + _setColLowerBound(cols(id(c)),lower);
1.1747 + _setColUpperBound(cols(id(c)),upper);
1.1748 }
1.1749
1.1750 ///\brief Set the lower and the upper bound of several columns
1.1751 - ///(i.e a variables) at once
1.1752 + ///(i.e variables) at once
1.1753 ///
1.1754 ///This magic function takes a container as its argument
1.1755 ///and applies the function on all of its elements.
1.1756 @@ -1222,23 +1433,21 @@
1.1757 void colBounds(T &t, Value lower, Value upper) { return 0;}
1.1758 #else
1.1759 template<class T>
1.1760 - typename enable_if<typename T::value_type::LpSolverCol,void>::type
1.1761 + typename enable_if<typename T::value_type::LpCol,void>::type
1.1762 colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
1.1763 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1764 colBounds(*i, lower, upper);
1.1765 }
1.1766 }
1.1767 template<class T>
1.1768 - typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.1769 - void>::type
1.1770 + typename enable_if<typename T::value_type::second_type::LpCol, void>::type
1.1771 colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
1.1772 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1773 colBounds(i->second, lower, upper);
1.1774 }
1.1775 }
1.1776 template<class T>
1.1777 - typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.1778 - void>::type
1.1779 + typename enable_if<typename T::MapIt::Value::LpCol, void>::type
1.1780 colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
1.1781 for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1782 colBounds(*i, lower, upper);
1.1783 @@ -1246,76 +1455,371 @@
1.1784 }
1.1785 #endif
1.1786
1.1787 + /// Set the lower bound of a row (i.e a constraint)
1.1788
1.1789 - /// Set the lower and the upper bounds of a row (i.e a constraint)
1.1790 -
1.1791 - /// The lower and the upper bound of a constraint (row) have to be
1.1792 - /// given by an extended number of type Value, i.e. a finite
1.1793 - /// number of type Value, -\ref INF or \ref INF. There is no
1.1794 - /// separate function for the lower and the upper bound because
1.1795 - /// that would have been hard to implement for CPLEX.
1.1796 - void rowBounds(Row c, Value lower, Value upper) {
1.1797 - _setRowBounds(_lpId(c),lower, upper);
1.1798 + /// The lower bound of a constraint (row) has to be given by an
1.1799 + /// extended number of type Value, i.e. a finite number of type
1.1800 + /// Value or -\ref INF.
1.1801 + void rowLowerBound(Row r, Value value) {
1.1802 + _setRowLowerBound(rows(id(r)),value);
1.1803 }
1.1804
1.1805 - /// Get the lower and the upper bounds of a row (i.e a constraint)
1.1806 + /// Get the lower bound of a row (i.e a constraint)
1.1807
1.1808 - /// The lower and the upper bound of
1.1809 - /// a constraint (row) are
1.1810 - /// extended numbers of type Value, i.e. finite numbers of type
1.1811 - /// Value, -\ref INF or \ref INF.
1.1812 - /// \todo There is no separate function for the
1.1813 - /// lower and the upper bound because we had problems with the
1.1814 - /// implementation of the setting functions for CPLEX:
1.1815 - /// check out whether this can be done for these functions.
1.1816 - void getRowBounds(Row c, Value &lower, Value &upper) const {
1.1817 - _getRowBounds(_lpId(c),lower, upper);
1.1818 + /// This function returns the lower bound for row (constraint) \c c
1.1819 + /// (this might be -\ref INF as well).
1.1820 + ///\return The lower bound for row \c r
1.1821 + Value rowLowerBound(Row r) const {
1.1822 + return _getRowLowerBound(rows(id(r)));
1.1823 + }
1.1824 +
1.1825 + /// Set the upper bound of a row (i.e a constraint)
1.1826 +
1.1827 + /// The upper bound of a constraint (row) has to be given by an
1.1828 + /// extended number of type Value, i.e. a finite number of type
1.1829 + /// Value or -\ref INF.
1.1830 + void rowUpperBound(Row r, Value value) {
1.1831 + _setRowUpperBound(rows(id(r)),value);
1.1832 + }
1.1833 +
1.1834 + /// Get the upper bound of a row (i.e a constraint)
1.1835 +
1.1836 + /// This function returns the upper bound for row (constraint) \c c
1.1837 + /// (this might be -\ref INF as well).
1.1838 + ///\return The upper bound for row \c r
1.1839 + Value rowUpperBound(Row r) const {
1.1840 + return _getRowUpperBound(rows(id(r)));
1.1841 }
1.1842
1.1843 ///Set an element of the objective function
1.1844 - void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); };
1.1845 + void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
1.1846
1.1847 ///Get an element of the objective function
1.1848 - Value objCoeff(Col c) const { return _getObjCoeff(_lpId(c)); };
1.1849 + Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
1.1850
1.1851 ///Set the objective function
1.1852
1.1853 ///\param e is a linear expression of type \ref Expr.
1.1854 - void obj(Expr e) {
1.1855 - _clearObj();
1.1856 - for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
1.1857 - objCoeff((*i).first,(*i).second);
1.1858 - obj_const_comp=e.constComp();
1.1859 + ///
1.1860 + void obj(const Expr& e) {
1.1861 + _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
1.1862 + ExprIterator(e.comps.end(), cols));
1.1863 + obj_const_comp = *e;
1.1864 }
1.1865
1.1866 ///Get the objective function
1.1867
1.1868 - ///\return the objective function as a linear expression of type \ref Expr.
1.1869 + ///\return the objective function as a linear expression of type
1.1870 + ///Expr.
1.1871 Expr obj() const {
1.1872 Expr e;
1.1873 - for (ColIt it(*this); it != INVALID; ++it) {
1.1874 - double c = objCoeff(it);
1.1875 - if (c != 0.0) {
1.1876 - e.insert(std::make_pair(it, c));
1.1877 - }
1.1878 - }
1.1879 + _getObjCoeffs(InsertIterator(e.comps, cols));
1.1880 + *e = obj_const_comp;
1.1881 return e;
1.1882 }
1.1883
1.1884
1.1885 - ///Maximize
1.1886 - void max() { _setMax(); }
1.1887 - ///Minimize
1.1888 - void min() { _setMin(); }
1.1889 + ///Set the direction of optimization
1.1890 + void sense(Sense sense) { _setSense(sense); }
1.1891
1.1892 - ///Query function: is this a maximization problem?
1.1893 - bool isMax() const {return _isMax(); }
1.1894 + ///Query the direction of the optimization
1.1895 + Sense sense() const {return _getSense(); }
1.1896
1.1897 - ///Query function: is this a minimization problem?
1.1898 - bool isMin() const {return !isMax(); }
1.1899 + ///Set the sense to maximization
1.1900 + void max() { _setSense(MAX); }
1.1901 +
1.1902 + ///Set the sense to maximization
1.1903 + void min() { _setSense(MIN); }
1.1904 +
1.1905 + ///Clears the problem
1.1906 + void clear() { _clear(); }
1.1907
1.1908 ///@}
1.1909
1.1910 + };
1.1911 +
1.1912 + /// Addition
1.1913 +
1.1914 + ///\relates LpBase::Expr
1.1915 + ///
1.1916 + inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
1.1917 + LpBase::Expr tmp(a);
1.1918 + tmp+=b;
1.1919 + return tmp;
1.1920 + }
1.1921 + ///Substraction
1.1922 +
1.1923 + ///\relates LpBase::Expr
1.1924 + ///
1.1925 + inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
1.1926 + LpBase::Expr tmp(a);
1.1927 + tmp-=b;
1.1928 + return tmp;
1.1929 + }
1.1930 + ///Multiply with constant
1.1931 +
1.1932 + ///\relates LpBase::Expr
1.1933 + ///
1.1934 + inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
1.1935 + LpBase::Expr tmp(a);
1.1936 + tmp*=b;
1.1937 + return tmp;
1.1938 + }
1.1939 +
1.1940 + ///Multiply with constant
1.1941 +
1.1942 + ///\relates LpBase::Expr
1.1943 + ///
1.1944 + inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
1.1945 + LpBase::Expr tmp(b);
1.1946 + tmp*=a;
1.1947 + return tmp;
1.1948 + }
1.1949 + ///Divide with constant
1.1950 +
1.1951 + ///\relates LpBase::Expr
1.1952 + ///
1.1953 + inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
1.1954 + LpBase::Expr tmp(a);
1.1955 + tmp/=b;
1.1956 + return tmp;
1.1957 + }
1.1958 +
1.1959 + ///Create constraint
1.1960 +
1.1961 + ///\relates LpBase::Constr
1.1962 + ///
1.1963 + inline LpBase::Constr operator<=(const LpBase::Expr &e,
1.1964 + const LpBase::Expr &f) {
1.1965 + return LpBase::Constr(0, f - e, LpBase::INF);
1.1966 + }
1.1967 +
1.1968 + ///Create constraint
1.1969 +
1.1970 + ///\relates LpBase::Constr
1.1971 + ///
1.1972 + inline LpBase::Constr operator<=(const LpBase::Value &e,
1.1973 + const LpBase::Expr &f) {
1.1974 + return LpBase::Constr(e, f, LpBase::NaN);
1.1975 + }
1.1976 +
1.1977 + ///Create constraint
1.1978 +
1.1979 + ///\relates LpBase::Constr
1.1980 + ///
1.1981 + inline LpBase::Constr operator<=(const LpBase::Expr &e,
1.1982 + const LpBase::Value &f) {
1.1983 + return LpBase::Constr(- LpBase::INF, e, f);
1.1984 + }
1.1985 +
1.1986 + ///Create constraint
1.1987 +
1.1988 + ///\relates LpBase::Constr
1.1989 + ///
1.1990 + inline LpBase::Constr operator>=(const LpBase::Expr &e,
1.1991 + const LpBase::Expr &f) {
1.1992 + return LpBase::Constr(0, e - f, LpBase::INF);
1.1993 + }
1.1994 +
1.1995 +
1.1996 + ///Create constraint
1.1997 +
1.1998 + ///\relates LpBase::Constr
1.1999 + ///
1.2000 + inline LpBase::Constr operator>=(const LpBase::Value &e,
1.2001 + const LpBase::Expr &f) {
1.2002 + return LpBase::Constr(LpBase::NaN, f, e);
1.2003 + }
1.2004 +
1.2005 +
1.2006 + ///Create constraint
1.2007 +
1.2008 + ///\relates LpBase::Constr
1.2009 + ///
1.2010 + inline LpBase::Constr operator>=(const LpBase::Expr &e,
1.2011 + const LpBase::Value &f) {
1.2012 + return LpBase::Constr(f, e, LpBase::INF);
1.2013 + }
1.2014 +
1.2015 + ///Create constraint
1.2016 +
1.2017 + ///\relates LpBase::Constr
1.2018 + ///
1.2019 + inline LpBase::Constr operator==(const LpBase::Expr &e,
1.2020 + const LpBase::Value &f) {
1.2021 + return LpBase::Constr(f, e, f);
1.2022 + }
1.2023 +
1.2024 + ///Create constraint
1.2025 +
1.2026 + ///\relates LpBase::Constr
1.2027 + ///
1.2028 + inline LpBase::Constr operator==(const LpBase::Expr &e,
1.2029 + const LpBase::Expr &f) {
1.2030 + return LpBase::Constr(0, f - e, 0);
1.2031 + }
1.2032 +
1.2033 + ///Create constraint
1.2034 +
1.2035 + ///\relates LpBase::Constr
1.2036 + ///
1.2037 + inline LpBase::Constr operator<=(const LpBase::Value &n,
1.2038 + const LpBase::Constr &c) {
1.2039 + LpBase::Constr tmp(c);
1.2040 + LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
1.2041 + tmp.lowerBound()=n;
1.2042 + return tmp;
1.2043 + }
1.2044 + ///Create constraint
1.2045 +
1.2046 + ///\relates LpBase::Constr
1.2047 + ///
1.2048 + inline LpBase::Constr operator<=(const LpBase::Constr &c,
1.2049 + const LpBase::Value &n)
1.2050 + {
1.2051 + LpBase::Constr tmp(c);
1.2052 + LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
1.2053 + tmp.upperBound()=n;
1.2054 + return tmp;
1.2055 + }
1.2056 +
1.2057 + ///Create constraint
1.2058 +
1.2059 + ///\relates LpBase::Constr
1.2060 + ///
1.2061 + inline LpBase::Constr operator>=(const LpBase::Value &n,
1.2062 + const LpBase::Constr &c) {
1.2063 + LpBase::Constr tmp(c);
1.2064 + LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint");
1.2065 + tmp.upperBound()=n;
1.2066 + return tmp;
1.2067 + }
1.2068 + ///Create constraint
1.2069 +
1.2070 + ///\relates LpBase::Constr
1.2071 + ///
1.2072 + inline LpBase::Constr operator>=(const LpBase::Constr &c,
1.2073 + const LpBase::Value &n)
1.2074 + {
1.2075 + LpBase::Constr tmp(c);
1.2076 + LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint");
1.2077 + tmp.lowerBound()=n;
1.2078 + return tmp;
1.2079 + }
1.2080 +
1.2081 + ///Addition
1.2082 +
1.2083 + ///\relates LpBase::DualExpr
1.2084 + ///
1.2085 + inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
1.2086 + const LpBase::DualExpr &b) {
1.2087 + LpBase::DualExpr tmp(a);
1.2088 + tmp+=b;
1.2089 + return tmp;
1.2090 + }
1.2091 + ///Substraction
1.2092 +
1.2093 + ///\relates LpBase::DualExpr
1.2094 + ///
1.2095 + inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
1.2096 + const LpBase::DualExpr &b) {
1.2097 + LpBase::DualExpr tmp(a);
1.2098 + tmp-=b;
1.2099 + return tmp;
1.2100 + }
1.2101 + ///Multiply with constant
1.2102 +
1.2103 + ///\relates LpBase::DualExpr
1.2104 + ///
1.2105 + inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
1.2106 + const LpBase::Value &b) {
1.2107 + LpBase::DualExpr tmp(a);
1.2108 + tmp*=b;
1.2109 + return tmp;
1.2110 + }
1.2111 +
1.2112 + ///Multiply with constant
1.2113 +
1.2114 + ///\relates LpBase::DualExpr
1.2115 + ///
1.2116 + inline LpBase::DualExpr operator*(const LpBase::Value &a,
1.2117 + const LpBase::DualExpr &b) {
1.2118 + LpBase::DualExpr tmp(b);
1.2119 + tmp*=a;
1.2120 + return tmp;
1.2121 + }
1.2122 + ///Divide with constant
1.2123 +
1.2124 + ///\relates LpBase::DualExpr
1.2125 + ///
1.2126 + inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
1.2127 + const LpBase::Value &b) {
1.2128 + LpBase::DualExpr tmp(a);
1.2129 + tmp/=b;
1.2130 + return tmp;
1.2131 + }
1.2132 +
1.2133 + /// \ingroup lp_group
1.2134 + ///
1.2135 + /// \brief Common base class for LP solvers
1.2136 + ///
1.2137 + /// This class is an abstract base class for LP solvers. This class
1.2138 + /// provides a full interface for set and modify an LP problem,
1.2139 + /// solve it and retrieve the solution. You can use one of the
1.2140 + /// descendants as a concrete implementation, or the \c Lp
1.2141 + /// default LP solver. However, if you would like to handle LP
1.2142 + /// solvers as reference or pointer in a generic way, you can use
1.2143 + /// this class directly.
1.2144 + class LpSolver : virtual public LpBase {
1.2145 + public:
1.2146 +
1.2147 + /// The problem types for primal and dual problems
1.2148 + enum ProblemType {
1.2149 + ///Feasible solution hasn't been found (but may exist).
1.2150 + UNDEFINED = 0,
1.2151 + ///The problem has no feasible solution
1.2152 + INFEASIBLE = 1,
1.2153 + ///Feasible solution found
1.2154 + FEASIBLE = 2,
1.2155 + ///Optimal solution exists and found
1.2156 + OPTIMAL = 3,
1.2157 + ///The cost function is unbounded
1.2158 + UNBOUNDED = 4
1.2159 + };
1.2160 +
1.2161 + ///The basis status of variables
1.2162 + enum VarStatus {
1.2163 + /// The variable is in the basis
1.2164 + BASIC,
1.2165 + /// The variable is free, but not basic
1.2166 + FREE,
1.2167 + /// The variable has active lower bound
1.2168 + LOWER,
1.2169 + /// The variable has active upper bound
1.2170 + UPPER,
1.2171 + /// The variable is non-basic and fixed
1.2172 + FIXED
1.2173 + };
1.2174 +
1.2175 + protected:
1.2176 +
1.2177 + virtual SolveExitStatus _solve() = 0;
1.2178 +
1.2179 + virtual Value _getPrimal(int i) const = 0;
1.2180 + virtual Value _getDual(int i) const = 0;
1.2181 +
1.2182 + virtual Value _getPrimalRay(int i) const = 0;
1.2183 + virtual Value _getDualRay(int i) const = 0;
1.2184 +
1.2185 + virtual Value _getPrimalValue() const = 0;
1.2186 +
1.2187 + virtual VarStatus _getColStatus(int i) const = 0;
1.2188 + virtual VarStatus _getRowStatus(int i) const = 0;
1.2189 +
1.2190 + virtual ProblemType _getPrimalType() const = 0;
1.2191 + virtual ProblemType _getDualType() const = 0;
1.2192 +
1.2193 + public:
1.2194
1.2195 ///\name Solve the LP
1.2196
1.2197 @@ -1326,8 +1830,6 @@
1.2198 ///\return The result of the optimization procedure. Possible
1.2199 ///values and their meanings can be found in the documentation of
1.2200 ///\ref SolveExitStatus.
1.2201 - ///
1.2202 - ///\todo Which method is used to solve the problem
1.2203 SolveExitStatus solve() { return _solve(); }
1.2204
1.2205 ///@}
1.2206 @@ -1336,368 +1838,241 @@
1.2207
1.2208 ///@{
1.2209
1.2210 - /// The status of the primal problem (the original LP problem)
1.2211 - SolutionStatus primalStatus() const {
1.2212 - return _getPrimalStatus();
1.2213 + /// The type of the primal problem
1.2214 + ProblemType primalType() const {
1.2215 + return _getPrimalType();
1.2216 }
1.2217
1.2218 - /// The status of the dual (of the original LP) problem
1.2219 - SolutionStatus dualStatus() const {
1.2220 - return _getDualStatus();
1.2221 + /// The type of the dual problem
1.2222 + ProblemType dualType() const {
1.2223 + return _getDualType();
1.2224 }
1.2225
1.2226 - ///The type of the original LP problem
1.2227 - ProblemTypes problemType() const {
1.2228 - return _getProblemType();
1.2229 + /// Return the primal value of the column
1.2230 +
1.2231 + /// Return the primal value of the column.
1.2232 + /// \pre The problem is solved.
1.2233 + Value primal(Col c) const { return _getPrimal(cols(id(c))); }
1.2234 +
1.2235 + /// Return the primal value of the expression
1.2236 +
1.2237 + /// Return the primal value of the expression, i.e. the dot
1.2238 + /// product of the primal solution and the expression.
1.2239 + /// \pre The problem is solved.
1.2240 + Value primal(const Expr& e) const {
1.2241 + double res = *e;
1.2242 + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1.2243 + res += *c * primal(c);
1.2244 + }
1.2245 + return res;
1.2246 }
1.2247 + /// Returns a component of the primal ray
1.2248 +
1.2249 + /// The primal ray is solution of the modified primal problem,
1.2250 + /// where we change each finite bound to 0, and we looking for a
1.2251 + /// negative objective value in case of minimization, and positive
1.2252 + /// objective value for maximization. If there is such solution,
1.2253 + /// that proofs the unsolvability of the dual problem, and if a
1.2254 + /// feasible primal solution exists, then the unboundness of
1.2255 + /// primal problem.
1.2256 + ///
1.2257 + /// \pre The problem is solved and the dual problem is infeasible.
1.2258 + /// \note Some solvers does not provide primal ray calculation
1.2259 + /// functions.
1.2260 + Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
1.2261
1.2262 - ///\e
1.2263 - Value primal(Col c) const { return _getPrimal(_lpId(c)); }
1.2264 - ///\e
1.2265 - Value primal(const Expr& e) const {
1.2266 - double res = e.constComp();
1.2267 - for (std::map<Col, double>::const_iterator it = e.begin();
1.2268 - it != e.end(); ++it) {
1.2269 - res += _getPrimal(_lpId(it->first)) * it->second;
1.2270 + /// Return the dual value of the row
1.2271 +
1.2272 + /// Return the dual value of the row.
1.2273 + /// \pre The problem is solved.
1.2274 + Value dual(Row r) const { return _getDual(rows(id(r))); }
1.2275 +
1.2276 + /// Return the dual value of the dual expression
1.2277 +
1.2278 + /// Return the dual value of the dual expression, i.e. the dot
1.2279 + /// product of the dual solution and the dual expression.
1.2280 + /// \pre The problem is solved.
1.2281 + Value dual(const DualExpr& e) const {
1.2282 + double res = 0.0;
1.2283 + for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
1.2284 + res += *r * dual(r);
1.2285 }
1.2286 return res;
1.2287 }
1.2288
1.2289 - ///\e
1.2290 - Value dual(Row r) const { return _getDual(_lpId(r)); }
1.2291 - ///\e
1.2292 - Value dual(const DualExpr& e) const {
1.2293 - double res = 0.0;
1.2294 - for (std::map<Row, double>::const_iterator it = e.begin();
1.2295 - it != e.end(); ++it) {
1.2296 - res += _getPrimal(_lpId(it->first)) * it->second;
1.2297 - }
1.2298 - return res;
1.2299 - }
1.2300 + /// Returns a component of the dual ray
1.2301 +
1.2302 + /// The dual ray is solution of the modified primal problem, where
1.2303 + /// we change each finite bound to 0 (i.e. the objective function
1.2304 + /// coefficients in the primal problem), and we looking for a
1.2305 + /// ositive objective value. If there is such solution, that
1.2306 + /// proofs the unsolvability of the primal problem, and if a
1.2307 + /// feasible dual solution exists, then the unboundness of
1.2308 + /// dual problem.
1.2309 + ///
1.2310 + /// \pre The problem is solved and the primal problem is infeasible.
1.2311 + /// \note Some solvers does not provide dual ray calculation
1.2312 + /// functions.
1.2313 + Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
1.2314
1.2315 - ///\e
1.2316 - bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); }
1.2317 + /// Return the basis status of the column
1.2318
1.2319 - ///\e
1.2320 + /// \see VarStatus
1.2321 + VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
1.2322 +
1.2323 + /// Return the basis status of the row
1.2324 +
1.2325 + /// \see VarStatus
1.2326 + VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
1.2327 +
1.2328 + ///The value of the objective function
1.2329
1.2330 ///\return
1.2331 ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1.2332 /// of the primal problem, depending on whether we minimize or maximize.
1.2333 ///- \ref NaN if no primal solution is found.
1.2334 ///- The (finite) objective value if an optimal solution is found.
1.2335 - Value primalValue() const { return _getPrimalValue()+obj_const_comp;}
1.2336 + Value primal() const { return _getPrimalValue()+obj_const_comp;}
1.2337 ///@}
1.2338
1.2339 + LpSolver* newSolver() {return _newSolver();}
1.2340 + LpSolver* cloneSolver() {return _cloneSolver();}
1.2341 +
1.2342 + protected:
1.2343 +
1.2344 + virtual LpSolver* _newSolver() const = 0;
1.2345 + virtual LpSolver* _cloneSolver() const = 0;
1.2346 };
1.2347
1.2348
1.2349 /// \ingroup lp_group
1.2350 ///
1.2351 /// \brief Common base class for MIP solvers
1.2352 - /// \todo Much more docs
1.2353 - class MipSolverBase : virtual public LpSolverBase{
1.2354 + ///
1.2355 + /// This class is an abstract base class for MIP solvers. This class
1.2356 + /// provides a full interface for set and modify an MIP problem,
1.2357 + /// solve it and retrieve the solution. You can use one of the
1.2358 + /// descendants as a concrete implementation, or the \c Lp
1.2359 + /// default MIP solver. However, if you would like to handle MIP
1.2360 + /// solvers as reference or pointer in a generic way, you can use
1.2361 + /// this class directly.
1.2362 + class MipSolver : virtual public LpBase {
1.2363 public:
1.2364
1.2365 - ///Possible variable (coloumn) types (e.g. real, integer, binary etc.)
1.2366 + /// The problem types for MIP problems
1.2367 + enum ProblemType {
1.2368 + ///Feasible solution hasn't been found (but may exist).
1.2369 + UNDEFINED = 0,
1.2370 + ///The problem has no feasible solution
1.2371 + INFEASIBLE = 1,
1.2372 + ///Feasible solution found
1.2373 + FEASIBLE = 2,
1.2374 + ///Optimal solution exists and found
1.2375 + OPTIMAL = 3,
1.2376 + ///The cost function is unbounded
1.2377 + ///
1.2378 + ///The Mip or at least the relaxed problem is unbounded
1.2379 + UNBOUNDED = 4
1.2380 + };
1.2381 +
1.2382 + ///\name Solve the MIP
1.2383 +
1.2384 + ///@{
1.2385 +
1.2386 + /// Solve the MIP problem at hand
1.2387 + ///
1.2388 + ///\return The result of the optimization procedure. Possible
1.2389 + ///values and their meanings can be found in the documentation of
1.2390 + ///\ref SolveExitStatus.
1.2391 + SolveExitStatus solve() { return _solve(); }
1.2392 +
1.2393 + ///@}
1.2394 +
1.2395 + ///\name Setting column type
1.2396 + ///@{
1.2397 +
1.2398 + ///Possible variable (column) types (e.g. real, integer, binary etc.)
1.2399 enum ColTypes {
1.2400 - ///Continuous variable
1.2401 + ///Continuous variable (default)
1.2402 REAL = 0,
1.2403 ///Integer variable
1.2404 -
1.2405 - ///Unfortunately, cplex 7.5 somewhere writes something like
1.2406 - ///#define INTEGER 'I'
1.2407 - INT = 1
1.2408 - ///\todo No support for other types yet.
1.2409 + INTEGER = 1
1.2410 };
1.2411
1.2412 - ///Sets the type of the given coloumn to the given type
1.2413 + ///Sets the type of the given column to the given type
1.2414 +
1.2415 + ///Sets the type of the given column to the given type.
1.2416 ///
1.2417 - ///Sets the type of the given coloumn to the given type.
1.2418 void colType(Col c, ColTypes col_type) {
1.2419 - _colType(_lpId(c),col_type);
1.2420 + _setColType(cols(id(c)),col_type);
1.2421 }
1.2422
1.2423 ///Gives back the type of the column.
1.2424 +
1.2425 + ///Gives back the type of the column.
1.2426 ///
1.2427 - ///Gives back the type of the column.
1.2428 ColTypes colType(Col c) const {
1.2429 - return _colType(_lpId(c));
1.2430 + return _getColType(cols(id(c)));
1.2431 + }
1.2432 + ///@}
1.2433 +
1.2434 + ///\name Obtain the solution
1.2435 +
1.2436 + ///@{
1.2437 +
1.2438 + /// The type of the MIP problem
1.2439 + ProblemType type() const {
1.2440 + return _getType();
1.2441 }
1.2442
1.2443 - ///Sets the type of the given Col to integer or remove that property.
1.2444 - ///
1.2445 - ///Sets the type of the given Col to integer or remove that property.
1.2446 - void integer(Col c, bool enable) {
1.2447 - if (enable)
1.2448 - colType(c,INT);
1.2449 - else
1.2450 - colType(c,REAL);
1.2451 + /// Return the value of the row in the solution
1.2452 +
1.2453 + /// Return the value of the row in the solution.
1.2454 + /// \pre The problem is solved.
1.2455 + Value sol(Col c) const { return _getSol(cols(id(c))); }
1.2456 +
1.2457 + /// Return the value of the expression in the solution
1.2458 +
1.2459 + /// Return the value of the expression in the solution, i.e. the
1.2460 + /// dot product of the solution and the expression.
1.2461 + /// \pre The problem is solved.
1.2462 + Value sol(const Expr& e) const {
1.2463 + double res = *e;
1.2464 + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1.2465 + res += *c * sol(c);
1.2466 + }
1.2467 + return res;
1.2468 }
1.2469 -
1.2470 - ///Gives back whether the type of the column is integer or not.
1.2471 - ///
1.2472 - ///Gives back the type of the column.
1.2473 - ///\return true if the column has integer type and false if not.
1.2474 - bool integer(Col c) const {
1.2475 - return (colType(c)==INT);
1.2476 - }
1.2477 -
1.2478 - /// The status of the MIP problem
1.2479 - SolutionStatus mipStatus() const {
1.2480 - return _getMipStatus();
1.2481 - }
1.2482 + ///The value of the objective function
1.2483 +
1.2484 + ///\return
1.2485 + ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1.2486 + /// of the problem, depending on whether we minimize or maximize.
1.2487 + ///- \ref NaN if no primal solution is found.
1.2488 + ///- The (finite) objective value if an optimal solution is found.
1.2489 + Value solValue() const { return _getSolValue()+obj_const_comp;}
1.2490 + ///@}
1.2491
1.2492 protected:
1.2493
1.2494 - virtual ColTypes _colType(int col) const = 0;
1.2495 - virtual void _colType(int col, ColTypes col_type) = 0;
1.2496 - virtual SolutionStatus _getMipStatus() const = 0;
1.2497 + virtual SolveExitStatus _solve() = 0;
1.2498 + virtual ColTypes _getColType(int col) const = 0;
1.2499 + virtual void _setColType(int col, ColTypes col_type) = 0;
1.2500 + virtual ProblemType _getType() const = 0;
1.2501 + virtual Value _getSol(int i) const = 0;
1.2502 + virtual Value _getSolValue() const = 0;
1.2503
1.2504 + public:
1.2505 +
1.2506 + MipSolver* newSolver() {return _newSolver();}
1.2507 + MipSolver* cloneSolver() {return _cloneSolver();}
1.2508 +
1.2509 + protected:
1.2510 +
1.2511 + virtual MipSolver* _newSolver() const = 0;
1.2512 + virtual MipSolver* _cloneSolver() const = 0;
1.2513 };
1.2514
1.2515 - ///\relates LpSolverBase::Expr
1.2516 - ///
1.2517 - inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
1.2518 - const LpSolverBase::Expr &b)
1.2519 - {
1.2520 - LpSolverBase::Expr tmp(a);
1.2521 - tmp+=b;
1.2522 - return tmp;
1.2523 - }
1.2524 - ///\e
1.2525 -
1.2526 - ///\relates LpSolverBase::Expr
1.2527 - ///
1.2528 - inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
1.2529 - const LpSolverBase::Expr &b)
1.2530 - {
1.2531 - LpSolverBase::Expr tmp(a);
1.2532 - tmp-=b;
1.2533 - return tmp;
1.2534 - }
1.2535 - ///\e
1.2536 -
1.2537 - ///\relates LpSolverBase::Expr
1.2538 - ///
1.2539 - inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
1.2540 - const LpSolverBase::Value &b)
1.2541 - {
1.2542 - LpSolverBase::Expr tmp(a);
1.2543 - tmp*=b;
1.2544 - return tmp;
1.2545 - }
1.2546 -
1.2547 - ///\e
1.2548 -
1.2549 - ///\relates LpSolverBase::Expr
1.2550 - ///
1.2551 - inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
1.2552 - const LpSolverBase::Expr &b)
1.2553 - {
1.2554 - LpSolverBase::Expr tmp(b);
1.2555 - tmp*=a;
1.2556 - return tmp;
1.2557 - }
1.2558 - ///\e
1.2559 -
1.2560 - ///\relates LpSolverBase::Expr
1.2561 - ///
1.2562 - inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
1.2563 - const LpSolverBase::Value &b)
1.2564 - {
1.2565 - LpSolverBase::Expr tmp(a);
1.2566 - tmp/=b;
1.2567 - return tmp;
1.2568 - }
1.2569 -
1.2570 - ///\e
1.2571 -
1.2572 - ///\relates LpSolverBase::Constr
1.2573 - ///
1.2574 - inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1.2575 - const LpSolverBase::Expr &f)
1.2576 - {
1.2577 - return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
1.2578 - }
1.2579 -
1.2580 - ///\e
1.2581 -
1.2582 - ///\relates LpSolverBase::Constr
1.2583 - ///
1.2584 - inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
1.2585 - const LpSolverBase::Expr &f)
1.2586 - {
1.2587 - return LpSolverBase::Constr(e,f);
1.2588 - }
1.2589 -
1.2590 - ///\e
1.2591 -
1.2592 - ///\relates LpSolverBase::Constr
1.2593 - ///
1.2594 - inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1.2595 - const LpSolverBase::Value &f)
1.2596 - {
1.2597 - return LpSolverBase::Constr(-LpSolverBase::INF,e,f);
1.2598 - }
1.2599 -
1.2600 - ///\e
1.2601 -
1.2602 - ///\relates LpSolverBase::Constr
1.2603 - ///
1.2604 - inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1.2605 - const LpSolverBase::Expr &f)
1.2606 - {
1.2607 - return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
1.2608 - }
1.2609 -
1.2610 -
1.2611 - ///\e
1.2612 -
1.2613 - ///\relates LpSolverBase::Constr
1.2614 - ///
1.2615 - inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
1.2616 - const LpSolverBase::Expr &f)
1.2617 - {
1.2618 - return LpSolverBase::Constr(f,e);
1.2619 - }
1.2620 -
1.2621 -
1.2622 - ///\e
1.2623 -
1.2624 - ///\relates LpSolverBase::Constr
1.2625 - ///
1.2626 - inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1.2627 - const LpSolverBase::Value &f)
1.2628 - {
1.2629 - return LpSolverBase::Constr(f,e,LpSolverBase::INF);
1.2630 - }
1.2631 -
1.2632 - ///\e
1.2633 -
1.2634 - ///\relates LpSolverBase::Constr
1.2635 - ///
1.2636 - inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1.2637 - const LpSolverBase::Value &f)
1.2638 - {
1.2639 - return LpSolverBase::Constr(f,e,f);
1.2640 - }
1.2641 -
1.2642 - ///\e
1.2643 -
1.2644 - ///\relates LpSolverBase::Constr
1.2645 - ///
1.2646 - inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1.2647 - const LpSolverBase::Expr &f)
1.2648 - {
1.2649 - return LpSolverBase::Constr(0,e-f,0);
1.2650 - }
1.2651 -
1.2652 - ///\e
1.2653 -
1.2654 - ///\relates LpSolverBase::Constr
1.2655 - ///
1.2656 - inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
1.2657 - const LpSolverBase::Constr&c)
1.2658 - {
1.2659 - LpSolverBase::Constr tmp(c);
1.2660 - LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.2661 - tmp.lowerBound()=n;
1.2662 - return tmp;
1.2663 - }
1.2664 - ///\e
1.2665 -
1.2666 - ///\relates LpSolverBase::Constr
1.2667 - ///
1.2668 - inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
1.2669 - const LpSolverBase::Value &n)
1.2670 - {
1.2671 - LpSolverBase::Constr tmp(c);
1.2672 - LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
1.2673 - tmp.upperBound()=n;
1.2674 - return tmp;
1.2675 - }
1.2676 -
1.2677 - ///\e
1.2678 -
1.2679 - ///\relates LpSolverBase::Constr
1.2680 - ///
1.2681 - inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
1.2682 - const LpSolverBase::Constr&c)
1.2683 - {
1.2684 - LpSolverBase::Constr tmp(c);
1.2685 - LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
1.2686 - tmp.upperBound()=n;
1.2687 - return tmp;
1.2688 - }
1.2689 - ///\e
1.2690 -
1.2691 - ///\relates LpSolverBase::Constr
1.2692 - ///
1.2693 - inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
1.2694 - const LpSolverBase::Value &n)
1.2695 - {
1.2696 - LpSolverBase::Constr tmp(c);
1.2697 - LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.2698 - tmp.lowerBound()=n;
1.2699 - return tmp;
1.2700 - }
1.2701 -
1.2702 - ///\e
1.2703 -
1.2704 - ///\relates LpSolverBase::DualExpr
1.2705 - ///
1.2706 - inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
1.2707 - const LpSolverBase::DualExpr &b)
1.2708 - {
1.2709 - LpSolverBase::DualExpr tmp(a);
1.2710 - tmp+=b;
1.2711 - return tmp;
1.2712 - }
1.2713 - ///\e
1.2714 -
1.2715 - ///\relates LpSolverBase::DualExpr
1.2716 - ///
1.2717 - inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
1.2718 - const LpSolverBase::DualExpr &b)
1.2719 - {
1.2720 - LpSolverBase::DualExpr tmp(a);
1.2721 - tmp-=b;
1.2722 - return tmp;
1.2723 - }
1.2724 - ///\e
1.2725 -
1.2726 - ///\relates LpSolverBase::DualExpr
1.2727 - ///
1.2728 - inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
1.2729 - const LpSolverBase::Value &b)
1.2730 - {
1.2731 - LpSolverBase::DualExpr tmp(a);
1.2732 - tmp*=b;
1.2733 - return tmp;
1.2734 - }
1.2735 -
1.2736 - ///\e
1.2737 -
1.2738 - ///\relates LpSolverBase::DualExpr
1.2739 - ///
1.2740 - inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
1.2741 - const LpSolverBase::DualExpr &b)
1.2742 - {
1.2743 - LpSolverBase::DualExpr tmp(b);
1.2744 - tmp*=a;
1.2745 - return tmp;
1.2746 - }
1.2747 - ///\e
1.2748 -
1.2749 - ///\relates LpSolverBase::DualExpr
1.2750 - ///
1.2751 - inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
1.2752 - const LpSolverBase::Value &b)
1.2753 - {
1.2754 - LpSolverBase::DualExpr tmp(a);
1.2755 - tmp/=b;
1.2756 - return tmp;
1.2757 - }
1.2758
1.2759
1.2760 } //namespace lemon