1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/lp_base.h Tue Dec 02 21:40:33 2008 +0100
1.3 @@ -0,0 +1,1705 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2008
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_LP_BASE_H
1.23 +#define LEMON_LP_BASE_H
1.24 +
1.25 +#include<iostream>
1.26 +#include<vector>
1.27 +#include<map>
1.28 +#include<limits>
1.29 +#include<lemon/math.h>
1.30 +
1.31 +#include<lemon/core.h>
1.32 +#include<lemon/bits/lp_id.h>
1.33 +
1.34 +///\file
1.35 +///\brief The interface of the LP solver interface.
1.36 +///\ingroup lp_group
1.37 +namespace lemon {
1.38 +
1.39 + /// Function to decide whether a floating point value is finite or not.
1.40 +
1.41 + /// Retruns true if the argument is not infinity, minus infinity or NaN.
1.42 + /// It does the same as the isfinite() function defined by C99.
1.43 + template <typename T>
1.44 + bool isFinite(T value)
1.45 + {
1.46 + typedef std::numeric_limits<T> Lim;
1.47 + if ((Lim::has_infinity && (value == Lim::infinity() || value ==
1.48 + -Lim::infinity())) ||
1.49 + ((Lim::has_quiet_NaN || Lim::has_signaling_NaN) && value != value))
1.50 + {
1.51 + return false;
1.52 + }
1.53 + return true;
1.54 + }
1.55 +
1.56 + ///Common base class for LP solvers
1.57 +
1.58 + ///\todo Much more docs
1.59 + ///\ingroup lp_group
1.60 + class LpSolverBase {
1.61 +
1.62 + protected:
1.63 +
1.64 + _lp_bits::LpId rows;
1.65 + _lp_bits::LpId cols;
1.66 +
1.67 + public:
1.68 +
1.69 + ///Possible outcomes of an LP solving procedure
1.70 + enum SolveExitStatus {
1.71 + ///This means that the problem has been successfully solved: either
1.72 + ///an optimal solution has been found or infeasibility/unboundedness
1.73 + ///has been proved.
1.74 + SOLVED = 0,
1.75 + ///Any other case (including the case when some user specified
1.76 + ///limit has been exceeded)
1.77 + UNSOLVED = 1
1.78 + };
1.79 +
1.80 + ///\e
1.81 + enum SolutionStatus {
1.82 + ///Feasible solution hasn't been found (but may exist).
1.83 +
1.84 + ///\todo NOTFOUND might be a better name.
1.85 + ///
1.86 + UNDEFINED = 0,
1.87 + ///The problem has no feasible solution
1.88 + INFEASIBLE = 1,
1.89 + ///Feasible solution found
1.90 + FEASIBLE = 2,
1.91 + ///Optimal solution exists and found
1.92 + OPTIMAL = 3,
1.93 + ///The cost function is unbounded
1.94 +
1.95 + ///\todo Give a feasible solution and an infinite ray (and the
1.96 + ///corresponding bases)
1.97 + INFINITE = 4
1.98 + };
1.99 +
1.100 + ///\e The type of the investigated LP problem
1.101 + enum ProblemTypes {
1.102 + ///Primal-dual feasible
1.103 + PRIMAL_DUAL_FEASIBLE = 0,
1.104 + ///Primal feasible dual infeasible
1.105 + PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
1.106 + ///Primal infeasible dual feasible
1.107 + PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
1.108 + ///Primal-dual infeasible
1.109 + PRIMAL_DUAL_INFEASIBLE = 3,
1.110 + ///Could not determine so far
1.111 + UNKNOWN = 4
1.112 + };
1.113 +
1.114 + ///The floating point type used by the solver
1.115 + typedef double Value;
1.116 + ///The infinity constant
1.117 + static const Value INF;
1.118 + ///The not a number constant
1.119 + static const Value NaN;
1.120 +
1.121 + static inline bool isNaN(const Value& v) { return v!=v; }
1.122 +
1.123 + friend class Col;
1.124 + friend class ColIt;
1.125 + friend class Row;
1.126 +
1.127 + ///Refer to a column of the LP.
1.128 +
1.129 + ///This type is used to refer to a column of the LP.
1.130 + ///
1.131 + ///Its value remains valid and correct even after the addition or erase of
1.132 + ///other columns.
1.133 + ///
1.134 + ///\todo Document what can one do with a Col (INVALID, comparing,
1.135 + ///it is similar to Node/Edge)
1.136 + class Col {
1.137 + protected:
1.138 + int id;
1.139 + friend class LpSolverBase;
1.140 + friend class MipSolverBase;
1.141 + explicit Col(int _id) : id(_id) {}
1.142 + public:
1.143 + typedef Value ExprValue;
1.144 + typedef True LpSolverCol;
1.145 + Col() {}
1.146 + Col(const Invalid&) : id(-1) {}
1.147 + bool operator< (Col c) const {return id< c.id;}
1.148 + bool operator> (Col c) const {return id> c.id;}
1.149 + bool operator==(Col c) const {return id==c.id;}
1.150 + bool operator!=(Col c) const {return id!=c.id;}
1.151 + };
1.152 +
1.153 + class ColIt : public Col {
1.154 + const LpSolverBase *_lp;
1.155 + public:
1.156 + ColIt() {}
1.157 + ColIt(const LpSolverBase &lp) : _lp(&lp)
1.158 + {
1.159 + _lp->cols.firstFix(id);
1.160 + }
1.161 + ColIt(const Invalid&) : Col(INVALID) {}
1.162 + ColIt &operator++()
1.163 + {
1.164 + _lp->cols.nextFix(id);
1.165 + return *this;
1.166 + }
1.167 + };
1.168 +
1.169 + static int id(const Col& col) { return col.id; }
1.170 +
1.171 +
1.172 + ///Refer to a row of the LP.
1.173 +
1.174 + ///This type is used to refer to a row of the LP.
1.175 + ///
1.176 + ///Its value remains valid and correct even after the addition or erase of
1.177 + ///other rows.
1.178 + ///
1.179 + ///\todo Document what can one do with a Row (INVALID, comparing,
1.180 + ///it is similar to Node/Edge)
1.181 + class Row {
1.182 + protected:
1.183 + int id;
1.184 + friend class LpSolverBase;
1.185 + explicit Row(int _id) : id(_id) {}
1.186 + public:
1.187 + typedef Value ExprValue;
1.188 + typedef True LpSolverRow;
1.189 + Row() {}
1.190 + Row(const Invalid&) : id(-1) {}
1.191 +
1.192 + bool operator< (Row c) const {return id< c.id;}
1.193 + bool operator> (Row c) const {return id> c.id;}
1.194 + bool operator==(Row c) const {return id==c.id;}
1.195 + bool operator!=(Row c) const {return id!=c.id;}
1.196 + };
1.197 +
1.198 + class RowIt : public Row {
1.199 + const LpSolverBase *_lp;
1.200 + public:
1.201 + RowIt() {}
1.202 + RowIt(const LpSolverBase &lp) : _lp(&lp)
1.203 + {
1.204 + _lp->rows.firstFix(id);
1.205 + }
1.206 + RowIt(const Invalid&) : Row(INVALID) {}
1.207 + RowIt &operator++()
1.208 + {
1.209 + _lp->rows.nextFix(id);
1.210 + return *this;
1.211 + }
1.212 + };
1.213 +
1.214 + static int id(const Row& row) { return row.id; }
1.215 +
1.216 + protected:
1.217 +
1.218 + int _lpId(const Col& c) const {
1.219 + return cols.floatingId(id(c));
1.220 + }
1.221 +
1.222 + int _lpId(const Row& r) const {
1.223 + return rows.floatingId(id(r));
1.224 + }
1.225 +
1.226 + Col _item(int i, Col) const {
1.227 + return Col(cols.fixId(i));
1.228 + }
1.229 +
1.230 + Row _item(int i, Row) const {
1.231 + return Row(rows.fixId(i));
1.232 + }
1.233 +
1.234 +
1.235 + public:
1.236 +
1.237 + ///Linear expression of variables and a constant component
1.238 +
1.239 + ///This data structure stores a linear expression of the variables
1.240 + ///(\ref Col "Col"s) and also has a constant component.
1.241 + ///
1.242 + ///There are several ways to access and modify the contents of this
1.243 + ///container.
1.244 + ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
1.245 + ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can
1.246 + ///read and modify the coefficients like
1.247 + ///these.
1.248 + ///\code
1.249 + ///e[v]=5;
1.250 + ///e[v]+=12;
1.251 + ///e.erase(v);
1.252 + ///\endcode
1.253 + ///or you can also iterate through its elements.
1.254 + ///\code
1.255 + ///double s=0;
1.256 + ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
1.257 + /// s+=i->second;
1.258 + ///\endcode
1.259 + ///(This code computes the sum of all coefficients).
1.260 + ///- Numbers (<tt>double</tt>'s)
1.261 + ///and variables (\ref Col "Col"s) directly convert to an
1.262 + ///\ref Expr and the usual linear operations are defined, so
1.263 + ///\code
1.264 + ///v+w
1.265 + ///2*v-3.12*(v-w/2)+2
1.266 + ///v*2.1+(3*v+(v*12+w+6)*3)/2
1.267 + ///\endcode
1.268 + ///are valid \ref Expr "Expr"essions.
1.269 + ///The usual assignment operations are also defined.
1.270 + ///\code
1.271 + ///e=v+w;
1.272 + ///e+=2*v-3.12*(v-w/2)+2;
1.273 + ///e*=3.4;
1.274 + ///e/=5;
1.275 + ///\endcode
1.276 + ///- The constant member can be set and read by \ref constComp()
1.277 + ///\code
1.278 + ///e.constComp()=12;
1.279 + ///double c=e.constComp();
1.280 + ///\endcode
1.281 + ///
1.282 + ///\note \ref clear() not only sets all coefficients to 0 but also
1.283 + ///clears the constant components.
1.284 + ///
1.285 + ///\sa Constr
1.286 + ///
1.287 + class Expr : public std::map<Col,Value>
1.288 + {
1.289 + public:
1.290 + typedef LpSolverBase::Col Key;
1.291 + typedef LpSolverBase::Value Value;
1.292 +
1.293 + protected:
1.294 + typedef std::map<Col,Value> Base;
1.295 +
1.296 + Value const_comp;
1.297 + public:
1.298 + typedef True IsLinExpression;
1.299 + ///\e
1.300 + Expr() : Base(), const_comp(0) { }
1.301 + ///\e
1.302 + Expr(const Key &v) : const_comp(0) {
1.303 + Base::insert(std::make_pair(v, 1));
1.304 + }
1.305 + ///\e
1.306 + Expr(const Value &v) : const_comp(v) {}
1.307 + ///\e
1.308 + void set(const Key &v,const Value &c) {
1.309 + Base::insert(std::make_pair(v, c));
1.310 + }
1.311 + ///\e
1.312 + Value &constComp() { return const_comp; }
1.313 + ///\e
1.314 + const Value &constComp() const { return const_comp; }
1.315 +
1.316 + ///Removes the components with zero coefficient.
1.317 + void simplify() {
1.318 + for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.319 + Base::iterator j=i;
1.320 + ++j;
1.321 + if ((*i).second==0) Base::erase(i);
1.322 + i=j;
1.323 + }
1.324 + }
1.325 +
1.326 + void simplify() const {
1.327 + const_cast<Expr*>(this)->simplify();
1.328 + }
1.329 +
1.330 + ///Removes the coefficients closer to zero than \c tolerance.
1.331 + void simplify(double &tolerance) {
1.332 + for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.333 + Base::iterator j=i;
1.334 + ++j;
1.335 + if (std::fabs((*i).second)<tolerance) Base::erase(i);
1.336 + i=j;
1.337 + }
1.338 + }
1.339 +
1.340 + ///Sets all coefficients and the constant component to 0.
1.341 + void clear() {
1.342 + Base::clear();
1.343 + const_comp=0;
1.344 + }
1.345 +
1.346 + ///\e
1.347 + Expr &operator+=(const Expr &e) {
1.348 + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.349 + (*this)[j->first]+=j->second;
1.350 + const_comp+=e.const_comp;
1.351 + return *this;
1.352 + }
1.353 + ///\e
1.354 + Expr &operator-=(const Expr &e) {
1.355 + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.356 + (*this)[j->first]-=j->second;
1.357 + const_comp-=e.const_comp;
1.358 + return *this;
1.359 + }
1.360 + ///\e
1.361 + Expr &operator*=(const Value &c) {
1.362 + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.363 + j->second*=c;
1.364 + const_comp*=c;
1.365 + return *this;
1.366 + }
1.367 + ///\e
1.368 + Expr &operator/=(const Value &c) {
1.369 + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.370 + j->second/=c;
1.371 + const_comp/=c;
1.372 + return *this;
1.373 + }
1.374 +
1.375 + };
1.376 +
1.377 + ///Linear constraint
1.378 +
1.379 + ///This data stucture represents a linear constraint in the LP.
1.380 + ///Basically it is a linear expression with a lower or an upper bound
1.381 + ///(or both). These parts of the constraint can be obtained by the member
1.382 + ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
1.383 + ///respectively.
1.384 + ///There are two ways to construct a constraint.
1.385 + ///- You can set the linear expression and the bounds directly
1.386 + /// by the functions above.
1.387 + ///- The operators <tt>\<=</tt>, <tt>==</tt> and <tt>\>=</tt>
1.388 + /// are defined between expressions, or even between constraints whenever
1.389 + /// it makes sense. Therefore if \c e and \c f are linear expressions and
1.390 + /// \c s and \c t are numbers, then the followings are valid expressions
1.391 + /// and thus they can be used directly e.g. in \ref addRow() whenever
1.392 + /// it makes sense.
1.393 + ///\code
1.394 + /// e<=s
1.395 + /// e<=f
1.396 + /// e==f
1.397 + /// s<=e<=t
1.398 + /// e>=t
1.399 + ///\endcode
1.400 + ///\warning The validity of a constraint is checked only at run time, so
1.401 + ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw
1.402 + ///an assertion.
1.403 + class Constr
1.404 + {
1.405 + public:
1.406 + typedef LpSolverBase::Expr Expr;
1.407 + typedef Expr::Key Key;
1.408 + typedef Expr::Value Value;
1.409 +
1.410 + protected:
1.411 + Expr _expr;
1.412 + Value _lb,_ub;
1.413 + public:
1.414 + ///\e
1.415 + Constr() : _expr(), _lb(NaN), _ub(NaN) {}
1.416 + ///\e
1.417 + Constr(Value lb,const Expr &e,Value ub) :
1.418 + _expr(e), _lb(lb), _ub(ub) {}
1.419 + ///\e
1.420 + Constr(const Expr &e,Value ub) :
1.421 + _expr(e), _lb(NaN), _ub(ub) {}
1.422 + ///\e
1.423 + Constr(Value lb,const Expr &e) :
1.424 + _expr(e), _lb(lb), _ub(NaN) {}
1.425 + ///\e
1.426 + Constr(const Expr &e) :
1.427 + _expr(e), _lb(NaN), _ub(NaN) {}
1.428 + ///\e
1.429 + void clear()
1.430 + {
1.431 + _expr.clear();
1.432 + _lb=_ub=NaN;
1.433 + }
1.434 +
1.435 + ///Reference to the linear expression
1.436 + Expr &expr() { return _expr; }
1.437 + ///Cont reference to the linear expression
1.438 + const Expr &expr() const { return _expr; }
1.439 + ///Reference to the lower bound.
1.440 +
1.441 + ///\return
1.442 + ///- \ref INF "INF": the constraint is lower unbounded.
1.443 + ///- \ref NaN "NaN": lower bound has not been set.
1.444 + ///- finite number: the lower bound
1.445 + Value &lowerBound() { return _lb; }
1.446 + ///The const version of \ref lowerBound()
1.447 + const Value &lowerBound() const { return _lb; }
1.448 + ///Reference to the upper bound.
1.449 +
1.450 + ///\return
1.451 + ///- \ref INF "INF": the constraint is upper unbounded.
1.452 + ///- \ref NaN "NaN": upper bound has not been set.
1.453 + ///- finite number: the upper bound
1.454 + Value &upperBound() { return _ub; }
1.455 + ///The const version of \ref upperBound()
1.456 + const Value &upperBound() const { return _ub; }
1.457 + ///Is the constraint lower bounded?
1.458 + bool lowerBounded() const {
1.459 + return isFinite(_lb);
1.460 + }
1.461 + ///Is the constraint upper bounded?
1.462 + bool upperBounded() const {
1.463 + return isFinite(_ub);
1.464 + }
1.465 +
1.466 + };
1.467 +
1.468 + ///Linear expression of rows
1.469 +
1.470 + ///This data structure represents a column of the matrix,
1.471 + ///thas is it strores a linear expression of the dual variables
1.472 + ///(\ref Row "Row"s).
1.473 + ///
1.474 + ///There are several ways to access and modify the contents of this
1.475 + ///container.
1.476 + ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
1.477 + ///if \c e is an DualExpr and \c v
1.478 + ///and \c w are of type \ref Row, then you can
1.479 + ///read and modify the coefficients like
1.480 + ///these.
1.481 + ///\code
1.482 + ///e[v]=5;
1.483 + ///e[v]+=12;
1.484 + ///e.erase(v);
1.485 + ///\endcode
1.486 + ///or you can also iterate through its elements.
1.487 + ///\code
1.488 + ///double s=0;
1.489 + ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
1.490 + /// s+=i->second;
1.491 + ///\endcode
1.492 + ///(This code computes the sum of all coefficients).
1.493 + ///- Numbers (<tt>double</tt>'s)
1.494 + ///and variables (\ref Row "Row"s) directly convert to an
1.495 + ///\ref DualExpr and the usual linear operations are defined, so
1.496 + ///\code
1.497 + ///v+w
1.498 + ///2*v-3.12*(v-w/2)
1.499 + ///v*2.1+(3*v+(v*12+w)*3)/2
1.500 + ///\endcode
1.501 + ///are valid \ref DualExpr "DualExpr"essions.
1.502 + ///The usual assignment operations are also defined.
1.503 + ///\code
1.504 + ///e=v+w;
1.505 + ///e+=2*v-3.12*(v-w/2);
1.506 + ///e*=3.4;
1.507 + ///e/=5;
1.508 + ///\endcode
1.509 + ///
1.510 + ///\sa Expr
1.511 + ///
1.512 + class DualExpr : public std::map<Row,Value>
1.513 + {
1.514 + public:
1.515 + typedef LpSolverBase::Row Key;
1.516 + typedef LpSolverBase::Value Value;
1.517 +
1.518 + protected:
1.519 + typedef std::map<Row,Value> Base;
1.520 +
1.521 + public:
1.522 + typedef True IsLinExpression;
1.523 + ///\e
1.524 + DualExpr() : Base() { }
1.525 + ///\e
1.526 + DualExpr(const Key &v) {
1.527 + Base::insert(std::make_pair(v, 1));
1.528 + }
1.529 + ///\e
1.530 + void set(const Key &v,const Value &c) {
1.531 + Base::insert(std::make_pair(v, c));
1.532 + }
1.533 +
1.534 + ///Removes the components with zero coefficient.
1.535 + void simplify() {
1.536 + for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.537 + Base::iterator j=i;
1.538 + ++j;
1.539 + if ((*i).second==0) Base::erase(i);
1.540 + i=j;
1.541 + }
1.542 + }
1.543 +
1.544 + void simplify() const {
1.545 + const_cast<DualExpr*>(this)->simplify();
1.546 + }
1.547 +
1.548 + ///Removes the coefficients closer to zero than \c tolerance.
1.549 + void simplify(double &tolerance) {
1.550 + for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.551 + Base::iterator j=i;
1.552 + ++j;
1.553 + if (std::fabs((*i).second)<tolerance) Base::erase(i);
1.554 + i=j;
1.555 + }
1.556 + }
1.557 +
1.558 + ///Sets all coefficients to 0.
1.559 + void clear() {
1.560 + Base::clear();
1.561 + }
1.562 +
1.563 + ///\e
1.564 + DualExpr &operator+=(const DualExpr &e) {
1.565 + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.566 + (*this)[j->first]+=j->second;
1.567 + return *this;
1.568 + }
1.569 + ///\e
1.570 + DualExpr &operator-=(const DualExpr &e) {
1.571 + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.572 + (*this)[j->first]-=j->second;
1.573 + return *this;
1.574 + }
1.575 + ///\e
1.576 + DualExpr &operator*=(const Value &c) {
1.577 + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.578 + j->second*=c;
1.579 + return *this;
1.580 + }
1.581 + ///\e
1.582 + DualExpr &operator/=(const Value &c) {
1.583 + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.584 + j->second/=c;
1.585 + return *this;
1.586 + }
1.587 + };
1.588 +
1.589 +
1.590 + private:
1.591 +
1.592 + template <typename _Expr>
1.593 + class MappedOutputIterator {
1.594 + public:
1.595 +
1.596 + typedef std::insert_iterator<_Expr> Base;
1.597 +
1.598 + typedef std::output_iterator_tag iterator_category;
1.599 + typedef void difference_type;
1.600 + typedef void value_type;
1.601 + typedef void reference;
1.602 + typedef void pointer;
1.603 +
1.604 + MappedOutputIterator(const Base& _base, const LpSolverBase& _lp)
1.605 + : base(_base), lp(_lp) {}
1.606 +
1.607 + MappedOutputIterator& operator*() {
1.608 + return *this;
1.609 + }
1.610 +
1.611 + MappedOutputIterator& operator=(const std::pair<int, Value>& value) {
1.612 + *base = std::make_pair(lp._item(value.first, typename _Expr::Key()),
1.613 + value.second);
1.614 + return *this;
1.615 + }
1.616 +
1.617 + MappedOutputIterator& operator++() {
1.618 + ++base;
1.619 + return *this;
1.620 + }
1.621 +
1.622 + MappedOutputIterator operator++(int) {
1.623 + MappedOutputIterator tmp(*this);
1.624 + ++base;
1.625 + return tmp;
1.626 + }
1.627 +
1.628 + bool operator==(const MappedOutputIterator& it) const {
1.629 + return base == it.base;
1.630 + }
1.631 +
1.632 + bool operator!=(const MappedOutputIterator& it) const {
1.633 + return base != it.base;
1.634 + }
1.635 +
1.636 + private:
1.637 + Base base;
1.638 + const LpSolverBase& lp;
1.639 + };
1.640 +
1.641 + template <typename Expr>
1.642 + class MappedInputIterator {
1.643 + public:
1.644 +
1.645 + typedef typename Expr::const_iterator Base;
1.646 +
1.647 + typedef typename Base::iterator_category iterator_category;
1.648 + typedef typename Base::difference_type difference_type;
1.649 + typedef const std::pair<int, Value> value_type;
1.650 + typedef value_type reference;
1.651 + class pointer {
1.652 + public:
1.653 + pointer(value_type& _value) : value(_value) {}
1.654 + value_type* operator->() { return &value; }
1.655 + private:
1.656 + value_type value;
1.657 + };
1.658 +
1.659 + MappedInputIterator(const Base& _base, const LpSolverBase& _lp)
1.660 + : base(_base), lp(_lp) {}
1.661 +
1.662 + reference operator*() {
1.663 + return std::make_pair(lp._lpId(base->first), base->second);
1.664 + }
1.665 +
1.666 + pointer operator->() {
1.667 + return pointer(operator*());
1.668 + }
1.669 +
1.670 + MappedInputIterator& operator++() {
1.671 + ++base;
1.672 + return *this;
1.673 + }
1.674 +
1.675 + MappedInputIterator operator++(int) {
1.676 + MappedInputIterator tmp(*this);
1.677 + ++base;
1.678 + return tmp;
1.679 + }
1.680 +
1.681 + bool operator==(const MappedInputIterator& it) const {
1.682 + return base == it.base;
1.683 + }
1.684 +
1.685 + bool operator!=(const MappedInputIterator& it) const {
1.686 + return base != it.base;
1.687 + }
1.688 +
1.689 + private:
1.690 + Base base;
1.691 + const LpSolverBase& lp;
1.692 + };
1.693 +
1.694 + protected:
1.695 +
1.696 + /// STL compatible iterator for lp col
1.697 + typedef MappedInputIterator<Expr> ConstRowIterator;
1.698 + /// STL compatible iterator for lp row
1.699 + typedef MappedInputIterator<DualExpr> ConstColIterator;
1.700 +
1.701 + /// STL compatible iterator for lp col
1.702 + typedef MappedOutputIterator<Expr> RowIterator;
1.703 + /// STL compatible iterator for lp row
1.704 + typedef MappedOutputIterator<DualExpr> ColIterator;
1.705 +
1.706 + //Abstract virtual functions
1.707 + virtual LpSolverBase* _newLp() = 0;
1.708 + virtual LpSolverBase* _copyLp(){
1.709 + LpSolverBase* newlp = _newLp();
1.710 +
1.711 + std::map<Col, Col> ref;
1.712 + for (LpSolverBase::ColIt it(*this); it != INVALID; ++it) {
1.713 + Col ccol = newlp->addCol();
1.714 + ref[it] = ccol;
1.715 + newlp->colName(ccol, colName(it));
1.716 + newlp->colLowerBound(ccol, colLowerBound(it));
1.717 + newlp->colUpperBound(ccol, colUpperBound(it));
1.718 + }
1.719 +
1.720 + for (LpSolverBase::RowIt it(*this); it != INVALID; ++it) {
1.721 + Expr e = row(it), ce;
1.722 + for (Expr::iterator jt = e.begin(); jt != e.end(); ++jt) {
1.723 + ce[ref[jt->first]] = jt->second;
1.724 + }
1.725 + ce += e.constComp();
1.726 + Row r = newlp->addRow(ce);
1.727 +
1.728 + double lower, upper;
1.729 + getRowBounds(it, lower, upper);
1.730 + newlp->rowBounds(r, lower, upper);
1.731 + }
1.732 +
1.733 + return newlp;
1.734 + };
1.735 +
1.736 + virtual int _addCol() = 0;
1.737 + virtual int _addRow() = 0;
1.738 +
1.739 + virtual void _eraseCol(int col) = 0;
1.740 + virtual void _eraseRow(int row) = 0;
1.741 +
1.742 + virtual void _getColName(int col, std::string & name) const = 0;
1.743 + virtual void _setColName(int col, const std::string & name) = 0;
1.744 + virtual int _colByName(const std::string& name) const = 0;
1.745 +
1.746 + virtual void _setRowCoeffs(int i, ConstRowIterator b,
1.747 + ConstRowIterator e) = 0;
1.748 + virtual void _getRowCoeffs(int i, RowIterator b) const = 0;
1.749 + virtual void _setColCoeffs(int i, ConstColIterator b,
1.750 + ConstColIterator e) = 0;
1.751 + virtual void _getColCoeffs(int i, ColIterator b) const = 0;
1.752 + virtual void _setCoeff(int row, int col, Value value) = 0;
1.753 + virtual Value _getCoeff(int row, int col) const = 0;
1.754 + virtual void _setColLowerBound(int i, Value value) = 0;
1.755 + virtual Value _getColLowerBound(int i) const = 0;
1.756 + virtual void _setColUpperBound(int i, Value value) = 0;
1.757 + virtual Value _getColUpperBound(int i) const = 0;
1.758 + virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
1.759 + virtual void _getRowBounds(int i, Value &lower, Value &upper) const = 0;
1.760 +
1.761 + virtual void _setObjCoeff(int i, Value obj_coef) = 0;
1.762 + virtual Value _getObjCoeff(int i) const = 0;
1.763 + virtual void _clearObj()=0;
1.764 +
1.765 + virtual SolveExitStatus _solve() = 0;
1.766 + virtual Value _getPrimal(int i) const = 0;
1.767 + virtual Value _getDual(int i) const = 0;
1.768 + virtual Value _getPrimalValue() const = 0;
1.769 + virtual bool _isBasicCol(int i) const = 0;
1.770 + virtual SolutionStatus _getPrimalStatus() const = 0;
1.771 + virtual SolutionStatus _getDualStatus() const = 0;
1.772 + virtual ProblemTypes _getProblemType() const = 0;
1.773 +
1.774 + virtual void _setMax() = 0;
1.775 + virtual void _setMin() = 0;
1.776 +
1.777 +
1.778 + virtual bool _isMax() const = 0;
1.779 +
1.780 + //Own protected stuff
1.781 +
1.782 + //Constant component of the objective function
1.783 + Value obj_const_comp;
1.784 +
1.785 + public:
1.786 +
1.787 + ///\e
1.788 + LpSolverBase() : obj_const_comp(0) {}
1.789 +
1.790 + ///\e
1.791 + virtual ~LpSolverBase() {}
1.792 +
1.793 + ///Creates a new LP problem
1.794 + LpSolverBase* newLp() {return _newLp();}
1.795 + ///Makes a copy of the LP problem
1.796 + LpSolverBase* copyLp() {return _copyLp();}
1.797 +
1.798 + ///\name Build up and modify the LP
1.799 +
1.800 + ///@{
1.801 +
1.802 + ///Add a new empty column (i.e a new variable) to the LP
1.803 + Col addCol() { Col c; _addCol(); c.id = cols.addId(); return c;}
1.804 +
1.805 + ///\brief Adds several new columns
1.806 + ///(i.e a variables) at once
1.807 + ///
1.808 + ///This magic function takes a container as its argument
1.809 + ///and fills its elements
1.810 + ///with new columns (i.e. variables)
1.811 + ///\param t can be
1.812 + ///- a standard STL compatible iterable container with
1.813 + ///\ref Col as its \c values_type
1.814 + ///like
1.815 + ///\code
1.816 + ///std::vector<LpSolverBase::Col>
1.817 + ///std::list<LpSolverBase::Col>
1.818 + ///\endcode
1.819 + ///- a standard STL compatible iterable container with
1.820 + ///\ref Col as its \c mapped_type
1.821 + ///like
1.822 + ///\code
1.823 + ///std::map<AnyType,LpSolverBase::Col>
1.824 + ///\endcode
1.825 + ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.826 + ///\code
1.827 + ///ListGraph::NodeMap<LpSolverBase::Col>
1.828 + ///ListGraph::EdgeMap<LpSolverBase::Col>
1.829 + ///\endcode
1.830 + ///\return The number of the created column.
1.831 +#ifdef DOXYGEN
1.832 + template<class T>
1.833 + int addColSet(T &t) { return 0;}
1.834 +#else
1.835 + template<class T>
1.836 + typename enable_if<typename T::value_type::LpSolverCol,int>::type
1.837 + addColSet(T &t,dummy<0> = 0) {
1.838 + int s=0;
1.839 + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
1.840 + return s;
1.841 + }
1.842 + template<class T>
1.843 + typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.844 + int>::type
1.845 + addColSet(T &t,dummy<1> = 1) {
1.846 + int s=0;
1.847 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.848 + i->second=addCol();
1.849 + s++;
1.850 + }
1.851 + return s;
1.852 + }
1.853 + template<class T>
1.854 + typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.855 + int>::type
1.856 + addColSet(T &t,dummy<2> = 2) {
1.857 + int s=0;
1.858 + for(typename T::MapIt i(t); i!=INVALID; ++i)
1.859 + {
1.860 + i.set(addCol());
1.861 + s++;
1.862 + }
1.863 + return s;
1.864 + }
1.865 +#endif
1.866 +
1.867 + ///Set a column (i.e a dual constraint) of the LP
1.868 +
1.869 + ///\param c is the column to be modified
1.870 + ///\param e is a dual linear expression (see \ref DualExpr)
1.871 + ///a better one.
1.872 + void col(Col c,const DualExpr &e) {
1.873 + e.simplify();
1.874 + _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this),
1.875 + ConstColIterator(e.end(), *this));
1.876 + }
1.877 +
1.878 + ///Get a column (i.e a dual constraint) of the LP
1.879 +
1.880 + ///\param r is the column to get
1.881 + ///\return the dual expression associated to the column
1.882 + DualExpr col(Col c) const {
1.883 + DualExpr e;
1.884 + _getColCoeffs(_lpId(c), ColIterator(std::inserter(e, e.end()), *this));
1.885 + return e;
1.886 + }
1.887 +
1.888 + ///Add a new column to the LP
1.889 +
1.890 + ///\param e is a dual linear expression (see \ref DualExpr)
1.891 + ///\param obj is the corresponding component of the objective
1.892 + ///function. It is 0 by default.
1.893 + ///\return The created column.
1.894 + Col addCol(const DualExpr &e, Value o = 0) {
1.895 + Col c=addCol();
1.896 + col(c,e);
1.897 + objCoeff(c,o);
1.898 + return c;
1.899 + }
1.900 +
1.901 + ///Add a new empty row (i.e a new constraint) to the LP
1.902 +
1.903 + ///This function adds a new empty row (i.e a new constraint) to the LP.
1.904 + ///\return The created row
1.905 + Row addRow() { Row r; _addRow(); r.id = rows.addId(); return r;}
1.906 +
1.907 + ///\brief Add several new rows
1.908 + ///(i.e a constraints) at once
1.909 + ///
1.910 + ///This magic function takes a container as its argument
1.911 + ///and fills its elements
1.912 + ///with new row (i.e. variables)
1.913 + ///\param t can be
1.914 + ///- a standard STL compatible iterable container with
1.915 + ///\ref Row as its \c values_type
1.916 + ///like
1.917 + ///\code
1.918 + ///std::vector<LpSolverBase::Row>
1.919 + ///std::list<LpSolverBase::Row>
1.920 + ///\endcode
1.921 + ///- a standard STL compatible iterable container with
1.922 + ///\ref Row as its \c mapped_type
1.923 + ///like
1.924 + ///\code
1.925 + ///std::map<AnyType,LpSolverBase::Row>
1.926 + ///\endcode
1.927 + ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.928 + ///\code
1.929 + ///ListGraph::NodeMap<LpSolverBase::Row>
1.930 + ///ListGraph::EdgeMap<LpSolverBase::Row>
1.931 + ///\endcode
1.932 + ///\return The number of rows created.
1.933 +#ifdef DOXYGEN
1.934 + template<class T>
1.935 + int addRowSet(T &t) { return 0;}
1.936 +#else
1.937 + template<class T>
1.938 + typename enable_if<typename T::value_type::LpSolverRow,int>::type
1.939 + addRowSet(T &t,dummy<0> = 0) {
1.940 + int s=0;
1.941 + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1.942 + return s;
1.943 + }
1.944 + template<class T>
1.945 + typename enable_if<typename T::value_type::second_type::LpSolverRow,
1.946 + int>::type
1.947 + addRowSet(T &t,dummy<1> = 1) {
1.948 + int s=0;
1.949 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.950 + i->second=addRow();
1.951 + s++;
1.952 + }
1.953 + return s;
1.954 + }
1.955 + template<class T>
1.956 + typename enable_if<typename T::MapIt::Value::LpSolverRow,
1.957 + int>::type
1.958 + addRowSet(T &t,dummy<2> = 2) {
1.959 + int s=0;
1.960 + for(typename T::MapIt i(t); i!=INVALID; ++i)
1.961 + {
1.962 + i.set(addRow());
1.963 + s++;
1.964 + }
1.965 + return s;
1.966 + }
1.967 +#endif
1.968 +
1.969 + ///Set a row (i.e a constraint) of the LP
1.970 +
1.971 + ///\param r is the row to be modified
1.972 + ///\param l is lower bound (-\ref INF means no bound)
1.973 + ///\param e is a linear expression (see \ref Expr)
1.974 + ///\param u is the upper bound (\ref INF means no bound)
1.975 + ///\bug This is a temporary function. The interface will change to
1.976 + ///a better one.
1.977 + ///\todo Option to control whether a constraint with a single variable is
1.978 + ///added or not.
1.979 + void row(Row r, Value l, const Expr &e, Value u) {
1.980 + e.simplify();
1.981 + _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this),
1.982 + ConstRowIterator(e.end(), *this));
1.983 + _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp());
1.984 + }
1.985 +
1.986 + ///Set a row (i.e a constraint) of the LP
1.987 +
1.988 + ///\param r is the row to be modified
1.989 + ///\param c is a linear expression (see \ref Constr)
1.990 + void row(Row r, const Constr &c) {
1.991 + row(r, c.lowerBounded()?c.lowerBound():-INF,
1.992 + c.expr(), c.upperBounded()?c.upperBound():INF);
1.993 + }
1.994 +
1.995 +
1.996 + ///Get a row (i.e a constraint) of the LP
1.997 +
1.998 + ///\param r is the row to get
1.999 + ///\return the expression associated to the row
1.1000 + Expr row(Row r) const {
1.1001 + Expr e;
1.1002 + _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this));
1.1003 + return e;
1.1004 + }
1.1005 +
1.1006 + ///Add a new row (i.e a new constraint) to the LP
1.1007 +
1.1008 + ///\param l is the lower bound (-\ref INF means no bound)
1.1009 + ///\param e is a linear expression (see \ref Expr)
1.1010 + ///\param u is the upper bound (\ref INF means no bound)
1.1011 + ///\return The created row.
1.1012 + ///\bug This is a temporary function. The interface will change to
1.1013 + ///a better one.
1.1014 + Row addRow(Value l,const Expr &e, Value u) {
1.1015 + Row r=addRow();
1.1016 + row(r,l,e,u);
1.1017 + return r;
1.1018 + }
1.1019 +
1.1020 + ///Add a new row (i.e a new constraint) to the LP
1.1021 +
1.1022 + ///\param c is a linear expression (see \ref Constr)
1.1023 + ///\return The created row.
1.1024 + Row addRow(const Constr &c) {
1.1025 + Row r=addRow();
1.1026 + row(r,c);
1.1027 + return r;
1.1028 + }
1.1029 + ///Erase a coloumn (i.e a variable) from the LP
1.1030 +
1.1031 + ///\param c is the coloumn to be deleted
1.1032 + ///\todo Please check this
1.1033 + void eraseCol(Col c) {
1.1034 + _eraseCol(_lpId(c));
1.1035 + cols.eraseId(c.id);
1.1036 + }
1.1037 + ///Erase a row (i.e a constraint) from the LP
1.1038 +
1.1039 + ///\param r is the row to be deleted
1.1040 + ///\todo Please check this
1.1041 + void eraseRow(Row r) {
1.1042 + _eraseRow(_lpId(r));
1.1043 + rows.eraseId(r.id);
1.1044 + }
1.1045 +
1.1046 + /// Get the name of a column
1.1047 +
1.1048 + ///\param c is the coresponding coloumn
1.1049 + ///\return The name of the colunm
1.1050 + std::string colName(Col c) const {
1.1051 + std::string name;
1.1052 + _getColName(_lpId(c), name);
1.1053 + return name;
1.1054 + }
1.1055 +
1.1056 + /// Set the name of a column
1.1057 +
1.1058 + ///\param c is the coresponding coloumn
1.1059 + ///\param name The name to be given
1.1060 + void colName(Col c, const std::string& name) {
1.1061 + _setColName(_lpId(c), name);
1.1062 + }
1.1063 +
1.1064 + /// Get the column by its name
1.1065 +
1.1066 + ///\param name The name of the column
1.1067 + ///\return the proper column or \c INVALID
1.1068 + Col colByName(const std::string& name) const {
1.1069 + int k = _colByName(name);
1.1070 + return k != -1 ? Col(cols.fixId(k)) : Col(INVALID);
1.1071 + }
1.1072 +
1.1073 + /// Set an element of the coefficient matrix of the LP
1.1074 +
1.1075 + ///\param r is the row of the element to be modified
1.1076 + ///\param c is the coloumn of the element to be modified
1.1077 + ///\param val is the new value of the coefficient
1.1078 +
1.1079 + void coeff(Row r, Col c, Value val) {
1.1080 + _setCoeff(_lpId(r),_lpId(c), val);
1.1081 + }
1.1082 +
1.1083 + /// Get an element of the coefficient matrix of the LP
1.1084 +
1.1085 + ///\param r is the row of the element in question
1.1086 + ///\param c is the coloumn of the element in question
1.1087 + ///\return the corresponding coefficient
1.1088 +
1.1089 + Value coeff(Row r, Col c) const {
1.1090 + return _getCoeff(_lpId(r),_lpId(c));
1.1091 + }
1.1092 +
1.1093 + /// Set the lower bound of a column (i.e a variable)
1.1094 +
1.1095 + /// The lower bound of a variable (column) has to be given by an
1.1096 + /// extended number of type Value, i.e. a finite number of type
1.1097 + /// Value or -\ref INF.
1.1098 + void colLowerBound(Col c, Value value) {
1.1099 + _setColLowerBound(_lpId(c),value);
1.1100 + }
1.1101 +
1.1102 + /// Get the lower bound of a column (i.e a variable)
1.1103 +
1.1104 + /// This function returns the lower bound for column (variable) \t c
1.1105 + /// (this might be -\ref INF as well).
1.1106 + ///\return The lower bound for coloumn \t c
1.1107 + Value colLowerBound(Col c) const {
1.1108 + return _getColLowerBound(_lpId(c));
1.1109 + }
1.1110 +
1.1111 + ///\brief Set the lower bound of several columns
1.1112 + ///(i.e a variables) at once
1.1113 + ///
1.1114 + ///This magic function takes a container as its argument
1.1115 + ///and applies the function on all of its elements.
1.1116 + /// The lower bound of a variable (column) has to be given by an
1.1117 + /// extended number of type Value, i.e. a finite number of type
1.1118 + /// Value or -\ref INF.
1.1119 +#ifdef DOXYGEN
1.1120 + template<class T>
1.1121 + void colLowerBound(T &t, Value value) { return 0;}
1.1122 +#else
1.1123 + template<class T>
1.1124 + typename enable_if<typename T::value_type::LpSolverCol,void>::type
1.1125 + colLowerBound(T &t, Value value,dummy<0> = 0) {
1.1126 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1127 + colLowerBound(*i, value);
1.1128 + }
1.1129 + }
1.1130 + template<class T>
1.1131 + typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.1132 + void>::type
1.1133 + colLowerBound(T &t, Value value,dummy<1> = 1) {
1.1134 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1135 + colLowerBound(i->second, value);
1.1136 + }
1.1137 + }
1.1138 + template<class T>
1.1139 + typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.1140 + void>::type
1.1141 + colLowerBound(T &t, Value value,dummy<2> = 2) {
1.1142 + for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1143 + colLowerBound(*i, value);
1.1144 + }
1.1145 + }
1.1146 +#endif
1.1147 +
1.1148 + /// Set the upper bound of a column (i.e a variable)
1.1149 +
1.1150 + /// The upper bound of a variable (column) has to be given by an
1.1151 + /// extended number of type Value, i.e. a finite number of type
1.1152 + /// Value or \ref INF.
1.1153 + void colUpperBound(Col c, Value value) {
1.1154 + _setColUpperBound(_lpId(c),value);
1.1155 + };
1.1156 +
1.1157 + /// Get the upper bound of a column (i.e a variable)
1.1158 +
1.1159 + /// This function returns the upper bound for column (variable) \t c
1.1160 + /// (this might be \ref INF as well).
1.1161 + ///\return The upper bound for coloumn \t c
1.1162 + Value colUpperBound(Col c) const {
1.1163 + return _getColUpperBound(_lpId(c));
1.1164 + }
1.1165 +
1.1166 + ///\brief Set the upper bound of several columns
1.1167 + ///(i.e a variables) at once
1.1168 + ///
1.1169 + ///This magic function takes a container as its argument
1.1170 + ///and applies the function on all of its elements.
1.1171 + /// The upper bound of a variable (column) has to be given by an
1.1172 + /// extended number of type Value, i.e. a finite number of type
1.1173 + /// Value or \ref INF.
1.1174 +#ifdef DOXYGEN
1.1175 + template<class T>
1.1176 + void colUpperBound(T &t, Value value) { return 0;}
1.1177 +#else
1.1178 + template<class T>
1.1179 + typename enable_if<typename T::value_type::LpSolverCol,void>::type
1.1180 + colUpperBound(T &t, Value value,dummy<0> = 0) {
1.1181 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1182 + colUpperBound(*i, value);
1.1183 + }
1.1184 + }
1.1185 + template<class T>
1.1186 + typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.1187 + void>::type
1.1188 + colUpperBound(T &t, Value value,dummy<1> = 1) {
1.1189 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1190 + colUpperBound(i->second, value);
1.1191 + }
1.1192 + }
1.1193 + template<class T>
1.1194 + typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.1195 + void>::type
1.1196 + colUpperBound(T &t, Value value,dummy<2> = 2) {
1.1197 + for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1198 + colUpperBound(*i, value);
1.1199 + }
1.1200 + }
1.1201 +#endif
1.1202 +
1.1203 + /// Set the lower and the upper bounds of a column (i.e a variable)
1.1204 +
1.1205 + /// The lower and the upper bounds of
1.1206 + /// a variable (column) have to be given by an
1.1207 + /// extended number of type Value, i.e. a finite number of type
1.1208 + /// Value, -\ref INF or \ref INF.
1.1209 + void colBounds(Col c, Value lower, Value upper) {
1.1210 + _setColLowerBound(_lpId(c),lower);
1.1211 + _setColUpperBound(_lpId(c),upper);
1.1212 + }
1.1213 +
1.1214 + ///\brief Set the lower and the upper bound of several columns
1.1215 + ///(i.e a variables) at once
1.1216 + ///
1.1217 + ///This magic function takes a container as its argument
1.1218 + ///and applies the function on all of its elements.
1.1219 + /// The lower and the upper bounds of
1.1220 + /// a variable (column) have to be given by an
1.1221 + /// extended number of type Value, i.e. a finite number of type
1.1222 + /// Value, -\ref INF or \ref INF.
1.1223 +#ifdef DOXYGEN
1.1224 + template<class T>
1.1225 + void colBounds(T &t, Value lower, Value upper) { return 0;}
1.1226 +#else
1.1227 + template<class T>
1.1228 + typename enable_if<typename T::value_type::LpSolverCol,void>::type
1.1229 + colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
1.1230 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1231 + colBounds(*i, lower, upper);
1.1232 + }
1.1233 + }
1.1234 + template<class T>
1.1235 + typename enable_if<typename T::value_type::second_type::LpSolverCol,
1.1236 + void>::type
1.1237 + colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
1.1238 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1239 + colBounds(i->second, lower, upper);
1.1240 + }
1.1241 + }
1.1242 + template<class T>
1.1243 + typename enable_if<typename T::MapIt::Value::LpSolverCol,
1.1244 + void>::type
1.1245 + colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
1.1246 + for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1247 + colBounds(*i, lower, upper);
1.1248 + }
1.1249 + }
1.1250 +#endif
1.1251 +
1.1252 +
1.1253 + /// Set the lower and the upper bounds of a row (i.e a constraint)
1.1254 +
1.1255 + /// The lower and the upper bound of a constraint (row) have to be
1.1256 + /// given by an extended number of type Value, i.e. a finite
1.1257 + /// number of type Value, -\ref INF or \ref INF. There is no
1.1258 + /// separate function for the lower and the upper bound because
1.1259 + /// that would have been hard to implement for CPLEX.
1.1260 + void rowBounds(Row c, Value lower, Value upper) {
1.1261 + _setRowBounds(_lpId(c),lower, upper);
1.1262 + }
1.1263 +
1.1264 + /// Get the lower and the upper bounds of a row (i.e a constraint)
1.1265 +
1.1266 + /// The lower and the upper bound of
1.1267 + /// a constraint (row) are
1.1268 + /// extended numbers of type Value, i.e. finite numbers of type
1.1269 + /// Value, -\ref INF or \ref INF.
1.1270 + /// \todo There is no separate function for the
1.1271 + /// lower and the upper bound because we had problems with the
1.1272 + /// implementation of the setting functions for CPLEX:
1.1273 + /// check out whether this can be done for these functions.
1.1274 + void getRowBounds(Row c, Value &lower, Value &upper) const {
1.1275 + _getRowBounds(_lpId(c),lower, upper);
1.1276 + }
1.1277 +
1.1278 + ///Set an element of the objective function
1.1279 + void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); };
1.1280 +
1.1281 + ///Get an element of the objective function
1.1282 + Value objCoeff(Col c) const { return _getObjCoeff(_lpId(c)); };
1.1283 +
1.1284 + ///Set the objective function
1.1285 +
1.1286 + ///\param e is a linear expression of type \ref Expr.
1.1287 + void obj(Expr e) {
1.1288 + _clearObj();
1.1289 + for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
1.1290 + objCoeff((*i).first,(*i).second);
1.1291 + obj_const_comp=e.constComp();
1.1292 + }
1.1293 +
1.1294 + ///Get the objective function
1.1295 +
1.1296 + ///\return the objective function as a linear expression of type \ref Expr.
1.1297 + Expr obj() const {
1.1298 + Expr e;
1.1299 + for (ColIt it(*this); it != INVALID; ++it) {
1.1300 + double c = objCoeff(it);
1.1301 + if (c != 0.0) {
1.1302 + e.insert(std::make_pair(it, c));
1.1303 + }
1.1304 + }
1.1305 + return e;
1.1306 + }
1.1307 +
1.1308 +
1.1309 + ///Maximize
1.1310 + void max() { _setMax(); }
1.1311 + ///Minimize
1.1312 + void min() { _setMin(); }
1.1313 +
1.1314 + ///Query function: is this a maximization problem?
1.1315 + bool isMax() const {return _isMax(); }
1.1316 +
1.1317 + ///Query function: is this a minimization problem?
1.1318 + bool isMin() const {return !isMax(); }
1.1319 +
1.1320 + ///@}
1.1321 +
1.1322 +
1.1323 + ///\name Solve the LP
1.1324 +
1.1325 + ///@{
1.1326 +
1.1327 + ///\e Solve the LP problem at hand
1.1328 + ///
1.1329 + ///\return The result of the optimization procedure. Possible
1.1330 + ///values and their meanings can be found in the documentation of
1.1331 + ///\ref SolveExitStatus.
1.1332 + ///
1.1333 + ///\todo Which method is used to solve the problem
1.1334 + SolveExitStatus solve() { return _solve(); }
1.1335 +
1.1336 + ///@}
1.1337 +
1.1338 + ///\name Obtain the solution
1.1339 +
1.1340 + ///@{
1.1341 +
1.1342 + /// The status of the primal problem (the original LP problem)
1.1343 + SolutionStatus primalStatus() const {
1.1344 + return _getPrimalStatus();
1.1345 + }
1.1346 +
1.1347 + /// The status of the dual (of the original LP) problem
1.1348 + SolutionStatus dualStatus() const {
1.1349 + return _getDualStatus();
1.1350 + }
1.1351 +
1.1352 + ///The type of the original LP problem
1.1353 + ProblemTypes problemType() const {
1.1354 + return _getProblemType();
1.1355 + }
1.1356 +
1.1357 + ///\e
1.1358 + Value primal(Col c) const { return _getPrimal(_lpId(c)); }
1.1359 + ///\e
1.1360 + Value primal(const Expr& e) const {
1.1361 + double res = e.constComp();
1.1362 + for (std::map<Col, double>::const_iterator it = e.begin();
1.1363 + it != e.end(); ++it) {
1.1364 + res += _getPrimal(_lpId(it->first)) * it->second;
1.1365 + }
1.1366 + return res;
1.1367 + }
1.1368 +
1.1369 + ///\e
1.1370 + Value dual(Row r) const { return _getDual(_lpId(r)); }
1.1371 + ///\e
1.1372 + Value dual(const DualExpr& e) const {
1.1373 + double res = 0.0;
1.1374 + for (std::map<Row, double>::const_iterator it = e.begin();
1.1375 + it != e.end(); ++it) {
1.1376 + res += _getPrimal(_lpId(it->first)) * it->second;
1.1377 + }
1.1378 + return res;
1.1379 + }
1.1380 +
1.1381 + ///\e
1.1382 + bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); }
1.1383 +
1.1384 + ///\e
1.1385 +
1.1386 + ///\return
1.1387 + ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1.1388 + /// of the primal problem, depending on whether we minimize or maximize.
1.1389 + ///- \ref NaN if no primal solution is found.
1.1390 + ///- The (finite) objective value if an optimal solution is found.
1.1391 + Value primalValue() const { return _getPrimalValue()+obj_const_comp;}
1.1392 + ///@}
1.1393 +
1.1394 + };
1.1395 +
1.1396 +
1.1397 + /// \ingroup lp_group
1.1398 + ///
1.1399 + /// \brief Common base class for MIP solvers
1.1400 + /// \todo Much more docs
1.1401 + class MipSolverBase : virtual public LpSolverBase{
1.1402 + public:
1.1403 +
1.1404 + ///Possible variable (coloumn) types (e.g. real, integer, binary etc.)
1.1405 + enum ColTypes {
1.1406 + ///Continuous variable
1.1407 + REAL = 0,
1.1408 + ///Integer variable
1.1409 +
1.1410 + ///Unfortunately, cplex 7.5 somewhere writes something like
1.1411 + ///#define INTEGER 'I'
1.1412 + INT = 1
1.1413 + ///\todo No support for other types yet.
1.1414 + };
1.1415 +
1.1416 + ///Sets the type of the given coloumn to the given type
1.1417 + ///
1.1418 + ///Sets the type of the given coloumn to the given type.
1.1419 + void colType(Col c, ColTypes col_type) {
1.1420 + _colType(_lpId(c),col_type);
1.1421 + }
1.1422 +
1.1423 + ///Gives back the type of the column.
1.1424 + ///
1.1425 + ///Gives back the type of the column.
1.1426 + ColTypes colType(Col c) const {
1.1427 + return _colType(_lpId(c));
1.1428 + }
1.1429 +
1.1430 + ///Sets the type of the given Col to integer or remove that property.
1.1431 + ///
1.1432 + ///Sets the type of the given Col to integer or remove that property.
1.1433 + void integer(Col c, bool enable) {
1.1434 + if (enable)
1.1435 + colType(c,INT);
1.1436 + else
1.1437 + colType(c,REAL);
1.1438 + }
1.1439 +
1.1440 + ///Gives back whether the type of the column is integer or not.
1.1441 + ///
1.1442 + ///Gives back the type of the column.
1.1443 + ///\return true if the column has integer type and false if not.
1.1444 + bool integer(Col c) const {
1.1445 + return (colType(c)==INT);
1.1446 + }
1.1447 +
1.1448 + /// The status of the MIP problem
1.1449 + SolutionStatus mipStatus() const {
1.1450 + return _getMipStatus();
1.1451 + }
1.1452 +
1.1453 + protected:
1.1454 +
1.1455 + virtual ColTypes _colType(int col) const = 0;
1.1456 + virtual void _colType(int col, ColTypes col_type) = 0;
1.1457 + virtual SolutionStatus _getMipStatus() const = 0;
1.1458 +
1.1459 + };
1.1460 +
1.1461 + ///\relates LpSolverBase::Expr
1.1462 + ///
1.1463 + inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
1.1464 + const LpSolverBase::Expr &b)
1.1465 + {
1.1466 + LpSolverBase::Expr tmp(a);
1.1467 + tmp+=b;
1.1468 + return tmp;
1.1469 + }
1.1470 + ///\e
1.1471 +
1.1472 + ///\relates LpSolverBase::Expr
1.1473 + ///
1.1474 + inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
1.1475 + const LpSolverBase::Expr &b)
1.1476 + {
1.1477 + LpSolverBase::Expr tmp(a);
1.1478 + tmp-=b;
1.1479 + return tmp;
1.1480 + }
1.1481 + ///\e
1.1482 +
1.1483 + ///\relates LpSolverBase::Expr
1.1484 + ///
1.1485 + inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
1.1486 + const LpSolverBase::Value &b)
1.1487 + {
1.1488 + LpSolverBase::Expr tmp(a);
1.1489 + tmp*=b;
1.1490 + return tmp;
1.1491 + }
1.1492 +
1.1493 + ///\e
1.1494 +
1.1495 + ///\relates LpSolverBase::Expr
1.1496 + ///
1.1497 + inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
1.1498 + const LpSolverBase::Expr &b)
1.1499 + {
1.1500 + LpSolverBase::Expr tmp(b);
1.1501 + tmp*=a;
1.1502 + return tmp;
1.1503 + }
1.1504 + ///\e
1.1505 +
1.1506 + ///\relates LpSolverBase::Expr
1.1507 + ///
1.1508 + inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
1.1509 + const LpSolverBase::Value &b)
1.1510 + {
1.1511 + LpSolverBase::Expr tmp(a);
1.1512 + tmp/=b;
1.1513 + return tmp;
1.1514 + }
1.1515 +
1.1516 + ///\e
1.1517 +
1.1518 + ///\relates LpSolverBase::Constr
1.1519 + ///
1.1520 + inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1.1521 + const LpSolverBase::Expr &f)
1.1522 + {
1.1523 + return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
1.1524 + }
1.1525 +
1.1526 + ///\e
1.1527 +
1.1528 + ///\relates LpSolverBase::Constr
1.1529 + ///
1.1530 + inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
1.1531 + const LpSolverBase::Expr &f)
1.1532 + {
1.1533 + return LpSolverBase::Constr(e,f);
1.1534 + }
1.1535 +
1.1536 + ///\e
1.1537 +
1.1538 + ///\relates LpSolverBase::Constr
1.1539 + ///
1.1540 + inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1.1541 + const LpSolverBase::Value &f)
1.1542 + {
1.1543 + return LpSolverBase::Constr(-LpSolverBase::INF,e,f);
1.1544 + }
1.1545 +
1.1546 + ///\e
1.1547 +
1.1548 + ///\relates LpSolverBase::Constr
1.1549 + ///
1.1550 + inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1.1551 + const LpSolverBase::Expr &f)
1.1552 + {
1.1553 + return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
1.1554 + }
1.1555 +
1.1556 +
1.1557 + ///\e
1.1558 +
1.1559 + ///\relates LpSolverBase::Constr
1.1560 + ///
1.1561 + inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
1.1562 + const LpSolverBase::Expr &f)
1.1563 + {
1.1564 + return LpSolverBase::Constr(f,e);
1.1565 + }
1.1566 +
1.1567 +
1.1568 + ///\e
1.1569 +
1.1570 + ///\relates LpSolverBase::Constr
1.1571 + ///
1.1572 + inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1.1573 + const LpSolverBase::Value &f)
1.1574 + {
1.1575 + return LpSolverBase::Constr(f,e,LpSolverBase::INF);
1.1576 + }
1.1577 +
1.1578 + ///\e
1.1579 +
1.1580 + ///\relates LpSolverBase::Constr
1.1581 + ///
1.1582 + inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1.1583 + const LpSolverBase::Value &f)
1.1584 + {
1.1585 + return LpSolverBase::Constr(f,e,f);
1.1586 + }
1.1587 +
1.1588 + ///\e
1.1589 +
1.1590 + ///\relates LpSolverBase::Constr
1.1591 + ///
1.1592 + inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1.1593 + const LpSolverBase::Expr &f)
1.1594 + {
1.1595 + return LpSolverBase::Constr(0,e-f,0);
1.1596 + }
1.1597 +
1.1598 + ///\e
1.1599 +
1.1600 + ///\relates LpSolverBase::Constr
1.1601 + ///
1.1602 + inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
1.1603 + const LpSolverBase::Constr&c)
1.1604 + {
1.1605 + LpSolverBase::Constr tmp(c);
1.1606 + LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.1607 + tmp.lowerBound()=n;
1.1608 + return tmp;
1.1609 + }
1.1610 + ///\e
1.1611 +
1.1612 + ///\relates LpSolverBase::Constr
1.1613 + ///
1.1614 + inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
1.1615 + const LpSolverBase::Value &n)
1.1616 + {
1.1617 + LpSolverBase::Constr tmp(c);
1.1618 + LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
1.1619 + tmp.upperBound()=n;
1.1620 + return tmp;
1.1621 + }
1.1622 +
1.1623 + ///\e
1.1624 +
1.1625 + ///\relates LpSolverBase::Constr
1.1626 + ///
1.1627 + inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
1.1628 + const LpSolverBase::Constr&c)
1.1629 + {
1.1630 + LpSolverBase::Constr tmp(c);
1.1631 + LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");
1.1632 + tmp.upperBound()=n;
1.1633 + return tmp;
1.1634 + }
1.1635 + ///\e
1.1636 +
1.1637 + ///\relates LpSolverBase::Constr
1.1638 + ///
1.1639 + inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
1.1640 + const LpSolverBase::Value &n)
1.1641 + {
1.1642 + LpSolverBase::Constr tmp(c);
1.1643 + LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.1644 + tmp.lowerBound()=n;
1.1645 + return tmp;
1.1646 + }
1.1647 +
1.1648 + ///\e
1.1649 +
1.1650 + ///\relates LpSolverBase::DualExpr
1.1651 + ///
1.1652 + inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
1.1653 + const LpSolverBase::DualExpr &b)
1.1654 + {
1.1655 + LpSolverBase::DualExpr tmp(a);
1.1656 + tmp+=b;
1.1657 + return tmp;
1.1658 + }
1.1659 + ///\e
1.1660 +
1.1661 + ///\relates LpSolverBase::DualExpr
1.1662 + ///
1.1663 + inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
1.1664 + const LpSolverBase::DualExpr &b)
1.1665 + {
1.1666 + LpSolverBase::DualExpr tmp(a);
1.1667 + tmp-=b;
1.1668 + return tmp;
1.1669 + }
1.1670 + ///\e
1.1671 +
1.1672 + ///\relates LpSolverBase::DualExpr
1.1673 + ///
1.1674 + inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
1.1675 + const LpSolverBase::Value &b)
1.1676 + {
1.1677 + LpSolverBase::DualExpr tmp(a);
1.1678 + tmp*=b;
1.1679 + return tmp;
1.1680 + }
1.1681 +
1.1682 + ///\e
1.1683 +
1.1684 + ///\relates LpSolverBase::DualExpr
1.1685 + ///
1.1686 + inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
1.1687 + const LpSolverBase::DualExpr &b)
1.1688 + {
1.1689 + LpSolverBase::DualExpr tmp(b);
1.1690 + tmp*=a;
1.1691 + return tmp;
1.1692 + }
1.1693 + ///\e
1.1694 +
1.1695 + ///\relates LpSolverBase::DualExpr
1.1696 + ///
1.1697 + inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
1.1698 + const LpSolverBase::Value &b)
1.1699 + {
1.1700 + LpSolverBase::DualExpr tmp(a);
1.1701 + tmp/=b;
1.1702 + return tmp;
1.1703 + }
1.1704 +
1.1705 +
1.1706 +} //namespace lemon
1.1707 +
1.1708 +#endif //LEMON_LP_BASE_H