1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/lp_base.h Thu Nov 05 15:50:01 2009 +0100
1.3 @@ -0,0 +1,2102 @@
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/error.h>
1.32 +#include<lemon/assert.h>
1.33 +
1.34 +#include<lemon/core.h>
1.35 +#include<lemon/bits/solver_bits.h>
1.36 +
1.37 +///\file
1.38 +///\brief The interface of the LP solver interface.
1.39 +///\ingroup lp_group
1.40 +namespace lemon {
1.41 +
1.42 + ///Common base class for LP and MIP solvers
1.43 +
1.44 + ///Usually this class is not used directly, please use one of the concrete
1.45 + ///implementations of the solver interface.
1.46 + ///\ingroup lp_group
1.47 + class LpBase {
1.48 +
1.49 + protected:
1.50 +
1.51 + _solver_bits::VarIndex rows;
1.52 + _solver_bits::VarIndex cols;
1.53 +
1.54 + public:
1.55 +
1.56 + ///Possible outcomes of an LP solving procedure
1.57 + enum SolveExitStatus {
1.58 + /// = 0. It means that the problem has been successfully solved: either
1.59 + ///an optimal solution has been found or infeasibility/unboundedness
1.60 + ///has been proved.
1.61 + SOLVED = 0,
1.62 + /// = 1. Any other case (including the case when some user specified
1.63 + ///limit has been exceeded).
1.64 + UNSOLVED = 1
1.65 + };
1.66 +
1.67 + ///Direction of the optimization
1.68 + enum Sense {
1.69 + /// Minimization
1.70 + MIN,
1.71 + /// Maximization
1.72 + MAX
1.73 + };
1.74 +
1.75 + ///Enum for \c messageLevel() parameter
1.76 + enum MessageLevel {
1.77 + /// No output (default value).
1.78 + MESSAGE_NOTHING,
1.79 + /// Error messages only.
1.80 + MESSAGE_ERROR,
1.81 + /// Warnings.
1.82 + MESSAGE_WARNING,
1.83 + /// Normal output.
1.84 + MESSAGE_NORMAL,
1.85 + /// Verbose output.
1.86 + MESSAGE_VERBOSE
1.87 + };
1.88 +
1.89 +
1.90 + ///The floating point type used by the solver
1.91 + typedef double Value;
1.92 + ///The infinity constant
1.93 + static const Value INF;
1.94 + ///The not a number constant
1.95 + static const Value NaN;
1.96 +
1.97 + friend class Col;
1.98 + friend class ColIt;
1.99 + friend class Row;
1.100 + friend class RowIt;
1.101 +
1.102 + ///Refer to a column of the LP.
1.103 +
1.104 + ///This type is used to refer to a column of the LP.
1.105 + ///
1.106 + ///Its value remains valid and correct even after the addition or erase of
1.107 + ///other columns.
1.108 + ///
1.109 + ///\note This class is similar to other Item types in LEMON, like
1.110 + ///Node and Arc types in digraph.
1.111 + class Col {
1.112 + friend class LpBase;
1.113 + protected:
1.114 + int _id;
1.115 + explicit Col(int id) : _id(id) {}
1.116 + public:
1.117 + typedef Value ExprValue;
1.118 + typedef True LpCol;
1.119 + /// Default constructor
1.120 +
1.121 + /// \warning The default constructor sets the Col to an
1.122 + /// undefined value.
1.123 + Col() {}
1.124 + /// Invalid constructor \& conversion.
1.125 +
1.126 + /// This constructor initializes the Col to be invalid.
1.127 + /// \sa Invalid for more details.
1.128 + Col(const Invalid&) : _id(-1) {}
1.129 + /// Equality operator
1.130 +
1.131 + /// Two \ref Col "Col"s are equal if and only if they point to
1.132 + /// the same LP column or both are invalid.
1.133 + bool operator==(Col c) const {return _id == c._id;}
1.134 + /// Inequality operator
1.135 +
1.136 + /// \sa operator==(Col c)
1.137 + ///
1.138 + bool operator!=(Col c) const {return _id != c._id;}
1.139 + /// Artificial ordering operator.
1.140 +
1.141 + /// To allow the use of this object in std::map or similar
1.142 + /// associative container we require this.
1.143 + ///
1.144 + /// \note This operator only have to define some strict ordering of
1.145 + /// the items; this order has nothing to do with the iteration
1.146 + /// ordering of the items.
1.147 + bool operator<(Col c) const {return _id < c._id;}
1.148 + };
1.149 +
1.150 + ///Iterator for iterate over the columns of an LP problem
1.151 +
1.152 + /// Its usage is quite simple, for example you can count the number
1.153 + /// of columns in an LP \c lp:
1.154 + ///\code
1.155 + /// int count=0;
1.156 + /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
1.157 + ///\endcode
1.158 + class ColIt : public Col {
1.159 + const LpBase *_solver;
1.160 + public:
1.161 + /// Default constructor
1.162 +
1.163 + /// \warning The default constructor sets the iterator
1.164 + /// to an undefined value.
1.165 + ColIt() {}
1.166 + /// Sets the iterator to the first Col
1.167 +
1.168 + /// Sets the iterator to the first Col.
1.169 + ///
1.170 + ColIt(const LpBase &solver) : _solver(&solver)
1.171 + {
1.172 + _solver->cols.firstItem(_id);
1.173 + }
1.174 + /// Invalid constructor \& conversion
1.175 +
1.176 + /// Initialize the iterator to be invalid.
1.177 + /// \sa Invalid for more details.
1.178 + ColIt(const Invalid&) : Col(INVALID) {}
1.179 + /// Next column
1.180 +
1.181 + /// Assign the iterator to the next column.
1.182 + ///
1.183 + ColIt &operator++()
1.184 + {
1.185 + _solver->cols.nextItem(_id);
1.186 + return *this;
1.187 + }
1.188 + };
1.189 +
1.190 + /// \brief Returns the ID of the column.
1.191 + static int id(const Col& col) { return col._id; }
1.192 + /// \brief Returns the column with the given ID.
1.193 + ///
1.194 + /// \pre The argument should be a valid column ID in the LP problem.
1.195 + static Col colFromId(int id) { return Col(id); }
1.196 +
1.197 + ///Refer to a row of the LP.
1.198 +
1.199 + ///This type is used to refer to a row of the LP.
1.200 + ///
1.201 + ///Its value remains valid and correct even after the addition or erase of
1.202 + ///other rows.
1.203 + ///
1.204 + ///\note This class is similar to other Item types in LEMON, like
1.205 + ///Node and Arc types in digraph.
1.206 + class Row {
1.207 + friend class LpBase;
1.208 + protected:
1.209 + int _id;
1.210 + explicit Row(int id) : _id(id) {}
1.211 + public:
1.212 + typedef Value ExprValue;
1.213 + typedef True LpRow;
1.214 + /// Default constructor
1.215 +
1.216 + /// \warning The default constructor sets the Row to an
1.217 + /// undefined value.
1.218 + Row() {}
1.219 + /// Invalid constructor \& conversion.
1.220 +
1.221 + /// This constructor initializes the Row to be invalid.
1.222 + /// \sa Invalid for more details.
1.223 + Row(const Invalid&) : _id(-1) {}
1.224 + /// Equality operator
1.225 +
1.226 + /// Two \ref Row "Row"s are equal if and only if they point to
1.227 + /// the same LP row or both are invalid.
1.228 + bool operator==(Row r) const {return _id == r._id;}
1.229 + /// Inequality operator
1.230 +
1.231 + /// \sa operator==(Row r)
1.232 + ///
1.233 + bool operator!=(Row r) const {return _id != r._id;}
1.234 + /// Artificial ordering operator.
1.235 +
1.236 + /// To allow the use of this object in std::map or similar
1.237 + /// associative container we require this.
1.238 + ///
1.239 + /// \note This operator only have to define some strict ordering of
1.240 + /// the items; this order has nothing to do with the iteration
1.241 + /// ordering of the items.
1.242 + bool operator<(Row r) const {return _id < r._id;}
1.243 + };
1.244 +
1.245 + ///Iterator for iterate over the rows of an LP problem
1.246 +
1.247 + /// Its usage is quite simple, for example you can count the number
1.248 + /// of rows in an LP \c lp:
1.249 + ///\code
1.250 + /// int count=0;
1.251 + /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
1.252 + ///\endcode
1.253 + class RowIt : public Row {
1.254 + const LpBase *_solver;
1.255 + public:
1.256 + /// Default constructor
1.257 +
1.258 + /// \warning The default constructor sets the iterator
1.259 + /// to an undefined value.
1.260 + RowIt() {}
1.261 + /// Sets the iterator to the first Row
1.262 +
1.263 + /// Sets the iterator to the first Row.
1.264 + ///
1.265 + RowIt(const LpBase &solver) : _solver(&solver)
1.266 + {
1.267 + _solver->rows.firstItem(_id);
1.268 + }
1.269 + /// Invalid constructor \& conversion
1.270 +
1.271 + /// Initialize the iterator to be invalid.
1.272 + /// \sa Invalid for more details.
1.273 + RowIt(const Invalid&) : Row(INVALID) {}
1.274 + /// Next row
1.275 +
1.276 + /// Assign the iterator to the next row.
1.277 + ///
1.278 + RowIt &operator++()
1.279 + {
1.280 + _solver->rows.nextItem(_id);
1.281 + return *this;
1.282 + }
1.283 + };
1.284 +
1.285 + /// \brief Returns the ID of the row.
1.286 + static int id(const Row& row) { return row._id; }
1.287 + /// \brief Returns the row with the given ID.
1.288 + ///
1.289 + /// \pre The argument should be a valid row ID in the LP problem.
1.290 + static Row rowFromId(int id) { return Row(id); }
1.291 +
1.292 + public:
1.293 +
1.294 + ///Linear expression of variables and a constant component
1.295 +
1.296 + ///This data structure stores a linear expression of the variables
1.297 + ///(\ref Col "Col"s) and also has a constant component.
1.298 + ///
1.299 + ///There are several ways to access and modify the contents of this
1.300 + ///container.
1.301 + ///\code
1.302 + ///e[v]=5;
1.303 + ///e[v]+=12;
1.304 + ///e.erase(v);
1.305 + ///\endcode
1.306 + ///or you can also iterate through its elements.
1.307 + ///\code
1.308 + ///double s=0;
1.309 + ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
1.310 + /// s+=*i * primal(i);
1.311 + ///\endcode
1.312 + ///(This code computes the primal value of the expression).
1.313 + ///- Numbers (<tt>double</tt>'s)
1.314 + ///and variables (\ref Col "Col"s) directly convert to an
1.315 + ///\ref Expr and the usual linear operations are defined, so
1.316 + ///\code
1.317 + ///v+w
1.318 + ///2*v-3.12*(v-w/2)+2
1.319 + ///v*2.1+(3*v+(v*12+w+6)*3)/2
1.320 + ///\endcode
1.321 + ///are valid expressions.
1.322 + ///The usual assignment operations are also defined.
1.323 + ///\code
1.324 + ///e=v+w;
1.325 + ///e+=2*v-3.12*(v-w/2)+2;
1.326 + ///e*=3.4;
1.327 + ///e/=5;
1.328 + ///\endcode
1.329 + ///- The constant member can be set and read by dereference
1.330 + /// operator (unary *)
1.331 + ///
1.332 + ///\code
1.333 + ///*e=12;
1.334 + ///double c=*e;
1.335 + ///\endcode
1.336 + ///
1.337 + ///\sa Constr
1.338 + class Expr {
1.339 + friend class LpBase;
1.340 + public:
1.341 + /// The key type of the expression
1.342 + typedef LpBase::Col Key;
1.343 + /// The value type of the expression
1.344 + typedef LpBase::Value Value;
1.345 +
1.346 + protected:
1.347 + Value const_comp;
1.348 + std::map<int, Value> comps;
1.349 +
1.350 + public:
1.351 + typedef True SolverExpr;
1.352 + /// Default constructor
1.353 +
1.354 + /// Construct an empty expression, the coefficients and
1.355 + /// the constant component are initialized to zero.
1.356 + Expr() : const_comp(0) {}
1.357 + /// Construct an expression from a column
1.358 +
1.359 + /// Construct an expression, which has a term with \c c variable
1.360 + /// and 1.0 coefficient.
1.361 + Expr(const Col &c) : const_comp(0) {
1.362 + typedef std::map<int, Value>::value_type pair_type;
1.363 + comps.insert(pair_type(id(c), 1));
1.364 + }
1.365 + /// Construct an expression from a constant
1.366 +
1.367 + /// Construct an expression, which's constant component is \c v.
1.368 + ///
1.369 + Expr(const Value &v) : const_comp(v) {}
1.370 + /// Returns the coefficient of the column
1.371 + Value operator[](const Col& c) const {
1.372 + std::map<int, Value>::const_iterator it=comps.find(id(c));
1.373 + if (it != comps.end()) {
1.374 + return it->second;
1.375 + } else {
1.376 + return 0;
1.377 + }
1.378 + }
1.379 + /// Returns the coefficient of the column
1.380 + Value& operator[](const Col& c) {
1.381 + return comps[id(c)];
1.382 + }
1.383 + /// Sets the coefficient of the column
1.384 + void set(const Col &c, const Value &v) {
1.385 + if (v != 0.0) {
1.386 + typedef std::map<int, Value>::value_type pair_type;
1.387 + comps.insert(pair_type(id(c), v));
1.388 + } else {
1.389 + comps.erase(id(c));
1.390 + }
1.391 + }
1.392 + /// Returns the constant component of the expression
1.393 + Value& operator*() { return const_comp; }
1.394 + /// Returns the constant component of the expression
1.395 + const Value& operator*() const { return const_comp; }
1.396 + /// \brief Removes the coefficients which's absolute value does
1.397 + /// not exceed \c epsilon. It also sets to zero the constant
1.398 + /// component, if it does not exceed epsilon in absolute value.
1.399 + void simplify(Value epsilon = 0.0) {
1.400 + std::map<int, Value>::iterator it=comps.begin();
1.401 + while (it != comps.end()) {
1.402 + std::map<int, Value>::iterator jt=it;
1.403 + ++jt;
1.404 + if (std::fabs((*it).second) <= epsilon) comps.erase(it);
1.405 + it=jt;
1.406 + }
1.407 + if (std::fabs(const_comp) <= epsilon) const_comp = 0;
1.408 + }
1.409 +
1.410 + void simplify(Value epsilon = 0.0) const {
1.411 + const_cast<Expr*>(this)->simplify(epsilon);
1.412 + }
1.413 +
1.414 + ///Sets all coefficients and the constant component to 0.
1.415 + void clear() {
1.416 + comps.clear();
1.417 + const_comp=0;
1.418 + }
1.419 +
1.420 + ///Compound assignment
1.421 + Expr &operator+=(const Expr &e) {
1.422 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.423 + it!=e.comps.end(); ++it)
1.424 + comps[it->first]+=it->second;
1.425 + const_comp+=e.const_comp;
1.426 + return *this;
1.427 + }
1.428 + ///Compound assignment
1.429 + Expr &operator-=(const Expr &e) {
1.430 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.431 + it!=e.comps.end(); ++it)
1.432 + comps[it->first]-=it->second;
1.433 + const_comp-=e.const_comp;
1.434 + return *this;
1.435 + }
1.436 + ///Multiply with a constant
1.437 + Expr &operator*=(const Value &v) {
1.438 + for (std::map<int, Value>::iterator it=comps.begin();
1.439 + it!=comps.end(); ++it)
1.440 + it->second*=v;
1.441 + const_comp*=v;
1.442 + return *this;
1.443 + }
1.444 + ///Division with a constant
1.445 + Expr &operator/=(const Value &c) {
1.446 + for (std::map<int, Value>::iterator it=comps.begin();
1.447 + it!=comps.end(); ++it)
1.448 + it->second/=c;
1.449 + const_comp/=c;
1.450 + return *this;
1.451 + }
1.452 +
1.453 + ///Iterator over the expression
1.454 +
1.455 + ///The iterator iterates over the terms of the expression.
1.456 + ///
1.457 + ///\code
1.458 + ///double s=0;
1.459 + ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
1.460 + /// s+= *i * primal(i);
1.461 + ///\endcode
1.462 + class CoeffIt {
1.463 + private:
1.464 +
1.465 + std::map<int, Value>::iterator _it, _end;
1.466 +
1.467 + public:
1.468 +
1.469 + /// Sets the iterator to the first term
1.470 +
1.471 + /// Sets the iterator to the first term of the expression.
1.472 + ///
1.473 + CoeffIt(Expr& e)
1.474 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.475 +
1.476 + /// Convert the iterator to the column of the term
1.477 + operator Col() const {
1.478 + return colFromId(_it->first);
1.479 + }
1.480 +
1.481 + /// Returns the coefficient of the term
1.482 + Value& operator*() { return _it->second; }
1.483 +
1.484 + /// Returns the coefficient of the term
1.485 + const Value& operator*() const { return _it->second; }
1.486 + /// Next term
1.487 +
1.488 + /// Assign the iterator to the next term.
1.489 + ///
1.490 + CoeffIt& operator++() { ++_it; return *this; }
1.491 +
1.492 + /// Equality operator
1.493 + bool operator==(Invalid) const { return _it == _end; }
1.494 + /// Inequality operator
1.495 + bool operator!=(Invalid) const { return _it != _end; }
1.496 + };
1.497 +
1.498 + /// Const iterator over the expression
1.499 +
1.500 + ///The iterator iterates over the terms of the expression.
1.501 + ///
1.502 + ///\code
1.503 + ///double s=0;
1.504 + ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
1.505 + /// s+=*i * primal(i);
1.506 + ///\endcode
1.507 + class ConstCoeffIt {
1.508 + private:
1.509 +
1.510 + std::map<int, Value>::const_iterator _it, _end;
1.511 +
1.512 + public:
1.513 +
1.514 + /// Sets the iterator to the first term
1.515 +
1.516 + /// Sets the iterator to the first term of the expression.
1.517 + ///
1.518 + ConstCoeffIt(const Expr& e)
1.519 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.520 +
1.521 + /// Convert the iterator to the column of the term
1.522 + operator Col() const {
1.523 + return colFromId(_it->first);
1.524 + }
1.525 +
1.526 + /// Returns the coefficient of the term
1.527 + const Value& operator*() const { return _it->second; }
1.528 +
1.529 + /// Next term
1.530 +
1.531 + /// Assign the iterator to the next term.
1.532 + ///
1.533 + ConstCoeffIt& operator++() { ++_it; return *this; }
1.534 +
1.535 + /// Equality operator
1.536 + bool operator==(Invalid) const { return _it == _end; }
1.537 + /// Inequality operator
1.538 + bool operator!=(Invalid) const { return _it != _end; }
1.539 + };
1.540 +
1.541 + };
1.542 +
1.543 + ///Linear constraint
1.544 +
1.545 + ///This data stucture represents a linear constraint in the LP.
1.546 + ///Basically it is a linear expression with a lower or an upper bound
1.547 + ///(or both). These parts of the constraint can be obtained by the member
1.548 + ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
1.549 + ///respectively.
1.550 + ///There are two ways to construct a constraint.
1.551 + ///- You can set the linear expression and the bounds directly
1.552 + /// by the functions above.
1.553 + ///- The operators <tt>\<=</tt>, <tt>==</tt> and <tt>\>=</tt>
1.554 + /// are defined between expressions, or even between constraints whenever
1.555 + /// it makes sense. Therefore if \c e and \c f are linear expressions and
1.556 + /// \c s and \c t are numbers, then the followings are valid expressions
1.557 + /// and thus they can be used directly e.g. in \ref addRow() whenever
1.558 + /// it makes sense.
1.559 + ///\code
1.560 + /// e<=s
1.561 + /// e<=f
1.562 + /// e==f
1.563 + /// s<=e<=t
1.564 + /// e>=t
1.565 + ///\endcode
1.566 + ///\warning The validity of a constraint is checked only at run
1.567 + ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
1.568 + ///compile, but will fail an assertion.
1.569 + class Constr
1.570 + {
1.571 + public:
1.572 + typedef LpBase::Expr Expr;
1.573 + typedef Expr::Key Key;
1.574 + typedef Expr::Value Value;
1.575 +
1.576 + protected:
1.577 + Expr _expr;
1.578 + Value _lb,_ub;
1.579 + public:
1.580 + ///\e
1.581 + Constr() : _expr(), _lb(NaN), _ub(NaN) {}
1.582 + ///\e
1.583 + Constr(Value lb, const Expr &e, Value ub) :
1.584 + _expr(e), _lb(lb), _ub(ub) {}
1.585 + Constr(const Expr &e) :
1.586 + _expr(e), _lb(NaN), _ub(NaN) {}
1.587 + ///\e
1.588 + void clear()
1.589 + {
1.590 + _expr.clear();
1.591 + _lb=_ub=NaN;
1.592 + }
1.593 +
1.594 + ///Reference to the linear expression
1.595 + Expr &expr() { return _expr; }
1.596 + ///Cont reference to the linear expression
1.597 + const Expr &expr() const { return _expr; }
1.598 + ///Reference to the lower bound.
1.599 +
1.600 + ///\return
1.601 + ///- \ref INF "INF": the constraint is lower unbounded.
1.602 + ///- \ref NaN "NaN": lower bound has not been set.
1.603 + ///- finite number: the lower bound
1.604 + Value &lowerBound() { return _lb; }
1.605 + ///The const version of \ref lowerBound()
1.606 + const Value &lowerBound() const { return _lb; }
1.607 + ///Reference to the upper bound.
1.608 +
1.609 + ///\return
1.610 + ///- \ref INF "INF": the constraint is upper unbounded.
1.611 + ///- \ref NaN "NaN": upper bound has not been set.
1.612 + ///- finite number: the upper bound
1.613 + Value &upperBound() { return _ub; }
1.614 + ///The const version of \ref upperBound()
1.615 + const Value &upperBound() const { return _ub; }
1.616 + ///Is the constraint lower bounded?
1.617 + bool lowerBounded() const {
1.618 + return _lb != -INF && !isNaN(_lb);
1.619 + }
1.620 + ///Is the constraint upper bounded?
1.621 + bool upperBounded() const {
1.622 + return _ub != INF && !isNaN(_ub);
1.623 + }
1.624 +
1.625 + };
1.626 +
1.627 + ///Linear expression of rows
1.628 +
1.629 + ///This data structure represents a column of the matrix,
1.630 + ///thas is it strores a linear expression of the dual variables
1.631 + ///(\ref Row "Row"s).
1.632 + ///
1.633 + ///There are several ways to access and modify the contents of this
1.634 + ///container.
1.635 + ///\code
1.636 + ///e[v]=5;
1.637 + ///e[v]+=12;
1.638 + ///e.erase(v);
1.639 + ///\endcode
1.640 + ///or you can also iterate through its elements.
1.641 + ///\code
1.642 + ///double s=0;
1.643 + ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
1.644 + /// s+=*i;
1.645 + ///\endcode
1.646 + ///(This code computes the sum of all coefficients).
1.647 + ///- Numbers (<tt>double</tt>'s)
1.648 + ///and variables (\ref Row "Row"s) directly convert to an
1.649 + ///\ref DualExpr and the usual linear operations are defined, so
1.650 + ///\code
1.651 + ///v+w
1.652 + ///2*v-3.12*(v-w/2)
1.653 + ///v*2.1+(3*v+(v*12+w)*3)/2
1.654 + ///\endcode
1.655 + ///are valid \ref DualExpr dual expressions.
1.656 + ///The usual assignment operations are also defined.
1.657 + ///\code
1.658 + ///e=v+w;
1.659 + ///e+=2*v-3.12*(v-w/2);
1.660 + ///e*=3.4;
1.661 + ///e/=5;
1.662 + ///\endcode
1.663 + ///
1.664 + ///\sa Expr
1.665 + class DualExpr {
1.666 + friend class LpBase;
1.667 + public:
1.668 + /// The key type of the expression
1.669 + typedef LpBase::Row Key;
1.670 + /// The value type of the expression
1.671 + typedef LpBase::Value Value;
1.672 +
1.673 + protected:
1.674 + std::map<int, Value> comps;
1.675 +
1.676 + public:
1.677 + typedef True SolverExpr;
1.678 + /// Default constructor
1.679 +
1.680 + /// Construct an empty expression, the coefficients are
1.681 + /// initialized to zero.
1.682 + DualExpr() {}
1.683 + /// Construct an expression from a row
1.684 +
1.685 + /// Construct an expression, which has a term with \c r dual
1.686 + /// variable and 1.0 coefficient.
1.687 + DualExpr(const Row &r) {
1.688 + typedef std::map<int, Value>::value_type pair_type;
1.689 + comps.insert(pair_type(id(r), 1));
1.690 + }
1.691 + /// Returns the coefficient of the row
1.692 + Value operator[](const Row& r) const {
1.693 + std::map<int, Value>::const_iterator it = comps.find(id(r));
1.694 + if (it != comps.end()) {
1.695 + return it->second;
1.696 + } else {
1.697 + return 0;
1.698 + }
1.699 + }
1.700 + /// Returns the coefficient of the row
1.701 + Value& operator[](const Row& r) {
1.702 + return comps[id(r)];
1.703 + }
1.704 + /// Sets the coefficient of the row
1.705 + void set(const Row &r, const Value &v) {
1.706 + if (v != 0.0) {
1.707 + typedef std::map<int, Value>::value_type pair_type;
1.708 + comps.insert(pair_type(id(r), v));
1.709 + } else {
1.710 + comps.erase(id(r));
1.711 + }
1.712 + }
1.713 + /// \brief Removes the coefficients which's absolute value does
1.714 + /// not exceed \c epsilon.
1.715 + void simplify(Value epsilon = 0.0) {
1.716 + std::map<int, Value>::iterator it=comps.begin();
1.717 + while (it != comps.end()) {
1.718 + std::map<int, Value>::iterator jt=it;
1.719 + ++jt;
1.720 + if (std::fabs((*it).second) <= epsilon) comps.erase(it);
1.721 + it=jt;
1.722 + }
1.723 + }
1.724 +
1.725 + void simplify(Value epsilon = 0.0) const {
1.726 + const_cast<DualExpr*>(this)->simplify(epsilon);
1.727 + }
1.728 +
1.729 + ///Sets all coefficients to 0.
1.730 + void clear() {
1.731 + comps.clear();
1.732 + }
1.733 + ///Compound assignment
1.734 + DualExpr &operator+=(const DualExpr &e) {
1.735 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.736 + it!=e.comps.end(); ++it)
1.737 + comps[it->first]+=it->second;
1.738 + return *this;
1.739 + }
1.740 + ///Compound assignment
1.741 + DualExpr &operator-=(const DualExpr &e) {
1.742 + for (std::map<int, Value>::const_iterator it=e.comps.begin();
1.743 + it!=e.comps.end(); ++it)
1.744 + comps[it->first]-=it->second;
1.745 + return *this;
1.746 + }
1.747 + ///Multiply with a constant
1.748 + DualExpr &operator*=(const Value &v) {
1.749 + for (std::map<int, Value>::iterator it=comps.begin();
1.750 + it!=comps.end(); ++it)
1.751 + it->second*=v;
1.752 + return *this;
1.753 + }
1.754 + ///Division with a constant
1.755 + DualExpr &operator/=(const Value &v) {
1.756 + for (std::map<int, Value>::iterator it=comps.begin();
1.757 + it!=comps.end(); ++it)
1.758 + it->second/=v;
1.759 + return *this;
1.760 + }
1.761 +
1.762 + ///Iterator over the expression
1.763 +
1.764 + ///The iterator iterates over the terms of the expression.
1.765 + ///
1.766 + ///\code
1.767 + ///double s=0;
1.768 + ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
1.769 + /// s+= *i * dual(i);
1.770 + ///\endcode
1.771 + class CoeffIt {
1.772 + private:
1.773 +
1.774 + std::map<int, Value>::iterator _it, _end;
1.775 +
1.776 + public:
1.777 +
1.778 + /// Sets the iterator to the first term
1.779 +
1.780 + /// Sets the iterator to the first term of the expression.
1.781 + ///
1.782 + CoeffIt(DualExpr& e)
1.783 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.784 +
1.785 + /// Convert the iterator to the row of the term
1.786 + operator Row() const {
1.787 + return rowFromId(_it->first);
1.788 + }
1.789 +
1.790 + /// Returns the coefficient of the term
1.791 + Value& operator*() { return _it->second; }
1.792 +
1.793 + /// Returns the coefficient of the term
1.794 + const Value& operator*() const { return _it->second; }
1.795 +
1.796 + /// Next term
1.797 +
1.798 + /// Assign the iterator to the next term.
1.799 + ///
1.800 + CoeffIt& operator++() { ++_it; return *this; }
1.801 +
1.802 + /// Equality operator
1.803 + bool operator==(Invalid) const { return _it == _end; }
1.804 + /// Inequality operator
1.805 + bool operator!=(Invalid) const { return _it != _end; }
1.806 + };
1.807 +
1.808 + ///Iterator over the expression
1.809 +
1.810 + ///The iterator iterates over the terms of the expression.
1.811 + ///
1.812 + ///\code
1.813 + ///double s=0;
1.814 + ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
1.815 + /// s+= *i * dual(i);
1.816 + ///\endcode
1.817 + class ConstCoeffIt {
1.818 + private:
1.819 +
1.820 + std::map<int, Value>::const_iterator _it, _end;
1.821 +
1.822 + public:
1.823 +
1.824 + /// Sets the iterator to the first term
1.825 +
1.826 + /// Sets the iterator to the first term of the expression.
1.827 + ///
1.828 + ConstCoeffIt(const DualExpr& e)
1.829 + : _it(e.comps.begin()), _end(e.comps.end()){}
1.830 +
1.831 + /// Convert the iterator to the row of the term
1.832 + operator Row() const {
1.833 + return rowFromId(_it->first);
1.834 + }
1.835 +
1.836 + /// Returns the coefficient of the term
1.837 + const Value& operator*() const { return _it->second; }
1.838 +
1.839 + /// Next term
1.840 +
1.841 + /// Assign the iterator to the next term.
1.842 + ///
1.843 + ConstCoeffIt& operator++() { ++_it; return *this; }
1.844 +
1.845 + /// Equality operator
1.846 + bool operator==(Invalid) const { return _it == _end; }
1.847 + /// Inequality operator
1.848 + bool operator!=(Invalid) const { return _it != _end; }
1.849 + };
1.850 + };
1.851 +
1.852 +
1.853 + protected:
1.854 +
1.855 + class InsertIterator {
1.856 + private:
1.857 +
1.858 + std::map<int, Value>& _host;
1.859 + const _solver_bits::VarIndex& _index;
1.860 +
1.861 + public:
1.862 +
1.863 + typedef std::output_iterator_tag iterator_category;
1.864 + typedef void difference_type;
1.865 + typedef void value_type;
1.866 + typedef void reference;
1.867 + typedef void pointer;
1.868 +
1.869 + InsertIterator(std::map<int, Value>& host,
1.870 + const _solver_bits::VarIndex& index)
1.871 + : _host(host), _index(index) {}
1.872 +
1.873 + InsertIterator& operator=(const std::pair<int, Value>& value) {
1.874 + typedef std::map<int, Value>::value_type pair_type;
1.875 + _host.insert(pair_type(_index[value.first], value.second));
1.876 + return *this;
1.877 + }
1.878 +
1.879 + InsertIterator& operator*() { return *this; }
1.880 + InsertIterator& operator++() { return *this; }
1.881 + InsertIterator operator++(int) { return *this; }
1.882 +
1.883 + };
1.884 +
1.885 + class ExprIterator {
1.886 + private:
1.887 + std::map<int, Value>::const_iterator _host_it;
1.888 + const _solver_bits::VarIndex& _index;
1.889 + public:
1.890 +
1.891 + typedef std::bidirectional_iterator_tag iterator_category;
1.892 + typedef std::ptrdiff_t difference_type;
1.893 + typedef const std::pair<int, Value> value_type;
1.894 + typedef value_type reference;
1.895 +
1.896 + class pointer {
1.897 + public:
1.898 + pointer(value_type& _value) : value(_value) {}
1.899 + value_type* operator->() { return &value; }
1.900 + private:
1.901 + value_type value;
1.902 + };
1.903 +
1.904 + ExprIterator(const std::map<int, Value>::const_iterator& host_it,
1.905 + const _solver_bits::VarIndex& index)
1.906 + : _host_it(host_it), _index(index) {}
1.907 +
1.908 + reference operator*() {
1.909 + return std::make_pair(_index(_host_it->first), _host_it->second);
1.910 + }
1.911 +
1.912 + pointer operator->() {
1.913 + return pointer(operator*());
1.914 + }
1.915 +
1.916 + ExprIterator& operator++() { ++_host_it; return *this; }
1.917 + ExprIterator operator++(int) {
1.918 + ExprIterator tmp(*this); ++_host_it; return tmp;
1.919 + }
1.920 +
1.921 + ExprIterator& operator--() { --_host_it; return *this; }
1.922 + ExprIterator operator--(int) {
1.923 + ExprIterator tmp(*this); --_host_it; return tmp;
1.924 + }
1.925 +
1.926 + bool operator==(const ExprIterator& it) const {
1.927 + return _host_it == it._host_it;
1.928 + }
1.929 +
1.930 + bool operator!=(const ExprIterator& it) const {
1.931 + return _host_it != it._host_it;
1.932 + }
1.933 +
1.934 + };
1.935 +
1.936 + protected:
1.937 +
1.938 + //Abstract virtual functions
1.939 +
1.940 + virtual int _addColId(int col) { return cols.addIndex(col); }
1.941 + virtual int _addRowId(int row) { return rows.addIndex(row); }
1.942 +
1.943 + virtual void _eraseColId(int col) { cols.eraseIndex(col); }
1.944 + virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
1.945 +
1.946 + virtual int _addCol() = 0;
1.947 + virtual int _addRow() = 0;
1.948 +
1.949 + virtual int _addRow(Value l, ExprIterator b, ExprIterator e, Value u) {
1.950 + int row = _addRow();
1.951 + _setRowCoeffs(row, b, e);
1.952 + _setRowLowerBound(row, l);
1.953 + _setRowUpperBound(row, u);
1.954 + return row;
1.955 + }
1.956 +
1.957 + virtual void _eraseCol(int col) = 0;
1.958 + virtual void _eraseRow(int row) = 0;
1.959 +
1.960 + virtual void _getColName(int col, std::string& name) const = 0;
1.961 + virtual void _setColName(int col, const std::string& name) = 0;
1.962 + virtual int _colByName(const std::string& name) const = 0;
1.963 +
1.964 + virtual void _getRowName(int row, std::string& name) const = 0;
1.965 + virtual void _setRowName(int row, const std::string& name) = 0;
1.966 + virtual int _rowByName(const std::string& name) const = 0;
1.967 +
1.968 + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1.969 + virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
1.970 +
1.971 + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1.972 + virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
1.973 +
1.974 + virtual void _setCoeff(int row, int col, Value value) = 0;
1.975 + virtual Value _getCoeff(int row, int col) const = 0;
1.976 +
1.977 + virtual void _setColLowerBound(int i, Value value) = 0;
1.978 + virtual Value _getColLowerBound(int i) const = 0;
1.979 +
1.980 + virtual void _setColUpperBound(int i, Value value) = 0;
1.981 + virtual Value _getColUpperBound(int i) const = 0;
1.982 +
1.983 + virtual void _setRowLowerBound(int i, Value value) = 0;
1.984 + virtual Value _getRowLowerBound(int i) const = 0;
1.985 +
1.986 + virtual void _setRowUpperBound(int i, Value value) = 0;
1.987 + virtual Value _getRowUpperBound(int i) const = 0;
1.988 +
1.989 + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
1.990 + virtual void _getObjCoeffs(InsertIterator b) const = 0;
1.991 +
1.992 + virtual void _setObjCoeff(int i, Value obj_coef) = 0;
1.993 + virtual Value _getObjCoeff(int i) const = 0;
1.994 +
1.995 + virtual void _setSense(Sense) = 0;
1.996 + virtual Sense _getSense() const = 0;
1.997 +
1.998 + virtual void _clear() = 0;
1.999 +
1.1000 + virtual const char* _solverName() const = 0;
1.1001 +
1.1002 + virtual void _messageLevel(MessageLevel level) = 0;
1.1003 +
1.1004 + //Own protected stuff
1.1005 +
1.1006 + //Constant component of the objective function
1.1007 + Value obj_const_comp;
1.1008 +
1.1009 + LpBase() : rows(), cols(), obj_const_comp(0) {}
1.1010 +
1.1011 + public:
1.1012 +
1.1013 + /// Virtual destructor
1.1014 + virtual ~LpBase() {}
1.1015 +
1.1016 + ///Gives back the name of the solver.
1.1017 + const char* solverName() const {return _solverName();}
1.1018 +
1.1019 + ///\name Build Up and Modify the LP
1.1020 +
1.1021 + ///@{
1.1022 +
1.1023 + ///Add a new empty column (i.e a new variable) to the LP
1.1024 + Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
1.1025 +
1.1026 + ///\brief Adds several new columns (i.e variables) at once
1.1027 + ///
1.1028 + ///This magic function takes a container as its argument and fills
1.1029 + ///its elements with new columns (i.e. variables)
1.1030 + ///\param t can be
1.1031 + ///- a standard STL compatible iterable container with
1.1032 + ///\ref Col as its \c values_type like
1.1033 + ///\code
1.1034 + ///std::vector<LpBase::Col>
1.1035 + ///std::list<LpBase::Col>
1.1036 + ///\endcode
1.1037 + ///- a standard STL compatible iterable container with
1.1038 + ///\ref Col as its \c mapped_type like
1.1039 + ///\code
1.1040 + ///std::map<AnyType,LpBase::Col>
1.1041 + ///\endcode
1.1042 + ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.1043 + ///\code
1.1044 + ///ListGraph::NodeMap<LpBase::Col>
1.1045 + ///ListGraph::ArcMap<LpBase::Col>
1.1046 + ///\endcode
1.1047 + ///\return The number of the created column.
1.1048 +#ifdef DOXYGEN
1.1049 + template<class T>
1.1050 + int addColSet(T &t) { return 0;}
1.1051 +#else
1.1052 + template<class T>
1.1053 + typename enable_if<typename T::value_type::LpCol,int>::type
1.1054 + addColSet(T &t,dummy<0> = 0) {
1.1055 + int s=0;
1.1056 + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
1.1057 + return s;
1.1058 + }
1.1059 + template<class T>
1.1060 + typename enable_if<typename T::value_type::second_type::LpCol,
1.1061 + int>::type
1.1062 + addColSet(T &t,dummy<1> = 1) {
1.1063 + int s=0;
1.1064 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1065 + i->second=addCol();
1.1066 + s++;
1.1067 + }
1.1068 + return s;
1.1069 + }
1.1070 + template<class T>
1.1071 + typename enable_if<typename T::MapIt::Value::LpCol,
1.1072 + int>::type
1.1073 + addColSet(T &t,dummy<2> = 2) {
1.1074 + int s=0;
1.1075 + for(typename T::MapIt i(t); i!=INVALID; ++i)
1.1076 + {
1.1077 + i.set(addCol());
1.1078 + s++;
1.1079 + }
1.1080 + return s;
1.1081 + }
1.1082 +#endif
1.1083 +
1.1084 + ///Set a column (i.e a dual constraint) of the LP
1.1085 +
1.1086 + ///\param c is the column to be modified
1.1087 + ///\param e is a dual linear expression (see \ref DualExpr)
1.1088 + ///a better one.
1.1089 + void col(Col c, const DualExpr &e) {
1.1090 + e.simplify();
1.1091 + _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), rows),
1.1092 + ExprIterator(e.comps.end(), rows));
1.1093 + }
1.1094 +
1.1095 + ///Get a column (i.e a dual constraint) of the LP
1.1096 +
1.1097 + ///\param c is the column to get
1.1098 + ///\return the dual expression associated to the column
1.1099 + DualExpr col(Col c) const {
1.1100 + DualExpr e;
1.1101 + _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
1.1102 + return e;
1.1103 + }
1.1104 +
1.1105 + ///Add a new column to the LP
1.1106 +
1.1107 + ///\param e is a dual linear expression (see \ref DualExpr)
1.1108 + ///\param o is the corresponding component of the objective
1.1109 + ///function. It is 0 by default.
1.1110 + ///\return The created column.
1.1111 + Col addCol(const DualExpr &e, Value o = 0) {
1.1112 + Col c=addCol();
1.1113 + col(c,e);
1.1114 + objCoeff(c,o);
1.1115 + return c;
1.1116 + }
1.1117 +
1.1118 + ///Add a new empty row (i.e a new constraint) to the LP
1.1119 +
1.1120 + ///This function adds a new empty row (i.e a new constraint) to the LP.
1.1121 + ///\return The created row
1.1122 + Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
1.1123 +
1.1124 + ///\brief Add several new rows (i.e constraints) at once
1.1125 + ///
1.1126 + ///This magic function takes a container as its argument and fills
1.1127 + ///its elements with new row (i.e. variables)
1.1128 + ///\param t can be
1.1129 + ///- a standard STL compatible iterable container with
1.1130 + ///\ref Row as its \c values_type like
1.1131 + ///\code
1.1132 + ///std::vector<LpBase::Row>
1.1133 + ///std::list<LpBase::Row>
1.1134 + ///\endcode
1.1135 + ///- a standard STL compatible iterable container with
1.1136 + ///\ref Row as its \c mapped_type like
1.1137 + ///\code
1.1138 + ///std::map<AnyType,LpBase::Row>
1.1139 + ///\endcode
1.1140 + ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.1141 + ///\code
1.1142 + ///ListGraph::NodeMap<LpBase::Row>
1.1143 + ///ListGraph::ArcMap<LpBase::Row>
1.1144 + ///\endcode
1.1145 + ///\return The number of rows created.
1.1146 +#ifdef DOXYGEN
1.1147 + template<class T>
1.1148 + int addRowSet(T &t) { return 0;}
1.1149 +#else
1.1150 + template<class T>
1.1151 + typename enable_if<typename T::value_type::LpRow,int>::type
1.1152 + addRowSet(T &t, dummy<0> = 0) {
1.1153 + int s=0;
1.1154 + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1.1155 + return s;
1.1156 + }
1.1157 + template<class T>
1.1158 + typename enable_if<typename T::value_type::second_type::LpRow, int>::type
1.1159 + addRowSet(T &t, dummy<1> = 1) {
1.1160 + int s=0;
1.1161 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1162 + i->second=addRow();
1.1163 + s++;
1.1164 + }
1.1165 + return s;
1.1166 + }
1.1167 + template<class T>
1.1168 + typename enable_if<typename T::MapIt::Value::LpRow, int>::type
1.1169 + addRowSet(T &t, dummy<2> = 2) {
1.1170 + int s=0;
1.1171 + for(typename T::MapIt i(t); i!=INVALID; ++i)
1.1172 + {
1.1173 + i.set(addRow());
1.1174 + s++;
1.1175 + }
1.1176 + return s;
1.1177 + }
1.1178 +#endif
1.1179 +
1.1180 + ///Set a row (i.e a constraint) of the LP
1.1181 +
1.1182 + ///\param r is the row to be modified
1.1183 + ///\param l is lower bound (-\ref INF means no bound)
1.1184 + ///\param e is a linear expression (see \ref Expr)
1.1185 + ///\param u is the upper bound (\ref INF means no bound)
1.1186 + void row(Row r, Value l, const Expr &e, Value u) {
1.1187 + e.simplify();
1.1188 + _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
1.1189 + ExprIterator(e.comps.end(), cols));
1.1190 + _setRowLowerBound(rows(id(r)),l - *e);
1.1191 + _setRowUpperBound(rows(id(r)),u - *e);
1.1192 + }
1.1193 +
1.1194 + ///Set a row (i.e a constraint) of the LP
1.1195 +
1.1196 + ///\param r is the row to be modified
1.1197 + ///\param c is a linear expression (see \ref Constr)
1.1198 + void row(Row r, const Constr &c) {
1.1199 + row(r, c.lowerBounded()?c.lowerBound():-INF,
1.1200 + c.expr(), c.upperBounded()?c.upperBound():INF);
1.1201 + }
1.1202 +
1.1203 +
1.1204 + ///Get a row (i.e a constraint) of the LP
1.1205 +
1.1206 + ///\param r is the row to get
1.1207 + ///\return the expression associated to the row
1.1208 + Expr row(Row r) const {
1.1209 + Expr e;
1.1210 + _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
1.1211 + return e;
1.1212 + }
1.1213 +
1.1214 + ///Add a new row (i.e a new constraint) to the LP
1.1215 +
1.1216 + ///\param l is the lower bound (-\ref INF means no bound)
1.1217 + ///\param e is a linear expression (see \ref Expr)
1.1218 + ///\param u is the upper bound (\ref INF means no bound)
1.1219 + ///\return The created row.
1.1220 + Row addRow(Value l,const Expr &e, Value u) {
1.1221 + Row r;
1.1222 + e.simplify();
1.1223 + r._id = _addRowId(_addRow(l - *e, ExprIterator(e.comps.begin(), cols),
1.1224 + ExprIterator(e.comps.end(), cols), u - *e));
1.1225 + return r;
1.1226 + }
1.1227 +
1.1228 + ///Add a new row (i.e a new constraint) to the LP
1.1229 +
1.1230 + ///\param c is a linear expression (see \ref Constr)
1.1231 + ///\return The created row.
1.1232 + Row addRow(const Constr &c) {
1.1233 + Row r;
1.1234 + c.expr().simplify();
1.1235 + r._id = _addRowId(_addRow(c.lowerBounded()?c.lowerBound():-INF,
1.1236 + ExprIterator(c.expr().comps.begin(), cols),
1.1237 + ExprIterator(c.expr().comps.end(), cols),
1.1238 + c.upperBounded()?c.upperBound():INF));
1.1239 + return r;
1.1240 + }
1.1241 + ///Erase a column (i.e a variable) from the LP
1.1242 +
1.1243 + ///\param c is the column to be deleted
1.1244 + void erase(Col c) {
1.1245 + _eraseCol(cols(id(c)));
1.1246 + _eraseColId(cols(id(c)));
1.1247 + }
1.1248 + ///Erase a row (i.e a constraint) from the LP
1.1249 +
1.1250 + ///\param r is the row to be deleted
1.1251 + void erase(Row r) {
1.1252 + _eraseRow(rows(id(r)));
1.1253 + _eraseRowId(rows(id(r)));
1.1254 + }
1.1255 +
1.1256 + /// Get the name of a column
1.1257 +
1.1258 + ///\param c is the coresponding column
1.1259 + ///\return The name of the colunm
1.1260 + std::string colName(Col c) const {
1.1261 + std::string name;
1.1262 + _getColName(cols(id(c)), name);
1.1263 + return name;
1.1264 + }
1.1265 +
1.1266 + /// Set the name of a column
1.1267 +
1.1268 + ///\param c is the coresponding column
1.1269 + ///\param name The name to be given
1.1270 + void colName(Col c, const std::string& name) {
1.1271 + _setColName(cols(id(c)), name);
1.1272 + }
1.1273 +
1.1274 + /// Get the column by its name
1.1275 +
1.1276 + ///\param name The name of the column
1.1277 + ///\return the proper column or \c INVALID
1.1278 + Col colByName(const std::string& name) const {
1.1279 + int k = _colByName(name);
1.1280 + return k != -1 ? Col(cols[k]) : Col(INVALID);
1.1281 + }
1.1282 +
1.1283 + /// Get the name of a row
1.1284 +
1.1285 + ///\param r is the coresponding row
1.1286 + ///\return The name of the row
1.1287 + std::string rowName(Row r) const {
1.1288 + std::string name;
1.1289 + _getRowName(rows(id(r)), name);
1.1290 + return name;
1.1291 + }
1.1292 +
1.1293 + /// Set the name of a row
1.1294 +
1.1295 + ///\param r is the coresponding row
1.1296 + ///\param name The name to be given
1.1297 + void rowName(Row r, const std::string& name) {
1.1298 + _setRowName(rows(id(r)), name);
1.1299 + }
1.1300 +
1.1301 + /// Get the row by its name
1.1302 +
1.1303 + ///\param name The name of the row
1.1304 + ///\return the proper row or \c INVALID
1.1305 + Row rowByName(const std::string& name) const {
1.1306 + int k = _rowByName(name);
1.1307 + return k != -1 ? Row(rows[k]) : Row(INVALID);
1.1308 + }
1.1309 +
1.1310 + /// Set an element of the coefficient matrix of the LP
1.1311 +
1.1312 + ///\param r is the row of the element to be modified
1.1313 + ///\param c is the column of the element to be modified
1.1314 + ///\param val is the new value of the coefficient
1.1315 + void coeff(Row r, Col c, Value val) {
1.1316 + _setCoeff(rows(id(r)),cols(id(c)), val);
1.1317 + }
1.1318 +
1.1319 + /// Get an element of the coefficient matrix of the LP
1.1320 +
1.1321 + ///\param r is the row of the element
1.1322 + ///\param c is the column of the element
1.1323 + ///\return the corresponding coefficient
1.1324 + Value coeff(Row r, Col c) const {
1.1325 + return _getCoeff(rows(id(r)),cols(id(c)));
1.1326 + }
1.1327 +
1.1328 + /// Set the lower bound of a column (i.e a variable)
1.1329 +
1.1330 + /// The lower bound of a variable (column) has to be given by an
1.1331 + /// extended number of type Value, i.e. a finite number of type
1.1332 + /// Value or -\ref INF.
1.1333 + void colLowerBound(Col c, Value value) {
1.1334 + _setColLowerBound(cols(id(c)),value);
1.1335 + }
1.1336 +
1.1337 + /// Get the lower bound of a column (i.e a variable)
1.1338 +
1.1339 + /// This function returns the lower bound for column (variable) \c c
1.1340 + /// (this might be -\ref INF as well).
1.1341 + ///\return The lower bound for column \c c
1.1342 + Value colLowerBound(Col c) const {
1.1343 + return _getColLowerBound(cols(id(c)));
1.1344 + }
1.1345 +
1.1346 + ///\brief Set the lower bound of several columns
1.1347 + ///(i.e variables) at once
1.1348 + ///
1.1349 + ///This magic function takes a container as its argument
1.1350 + ///and applies the function on all of its elements.
1.1351 + ///The lower bound of a variable (column) has to be given by an
1.1352 + ///extended number of type Value, i.e. a finite number of type
1.1353 + ///Value or -\ref INF.
1.1354 +#ifdef DOXYGEN
1.1355 + template<class T>
1.1356 + void colLowerBound(T &t, Value value) { return 0;}
1.1357 +#else
1.1358 + template<class T>
1.1359 + typename enable_if<typename T::value_type::LpCol,void>::type
1.1360 + colLowerBound(T &t, Value value,dummy<0> = 0) {
1.1361 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1362 + colLowerBound(*i, value);
1.1363 + }
1.1364 + }
1.1365 + template<class T>
1.1366 + typename enable_if<typename T::value_type::second_type::LpCol,
1.1367 + void>::type
1.1368 + colLowerBound(T &t, Value value,dummy<1> = 1) {
1.1369 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1370 + colLowerBound(i->second, value);
1.1371 + }
1.1372 + }
1.1373 + template<class T>
1.1374 + typename enable_if<typename T::MapIt::Value::LpCol,
1.1375 + void>::type
1.1376 + colLowerBound(T &t, Value value,dummy<2> = 2) {
1.1377 + for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1378 + colLowerBound(*i, value);
1.1379 + }
1.1380 + }
1.1381 +#endif
1.1382 +
1.1383 + /// Set the upper bound of a column (i.e a variable)
1.1384 +
1.1385 + /// The upper bound of a variable (column) has to be given by an
1.1386 + /// extended number of type Value, i.e. a finite number of type
1.1387 + /// Value or \ref INF.
1.1388 + void colUpperBound(Col c, Value value) {
1.1389 + _setColUpperBound(cols(id(c)),value);
1.1390 + };
1.1391 +
1.1392 + /// Get the upper bound of a column (i.e a variable)
1.1393 +
1.1394 + /// This function returns the upper bound for column (variable) \c c
1.1395 + /// (this might be \ref INF as well).
1.1396 + /// \return The upper bound for column \c c
1.1397 + Value colUpperBound(Col c) const {
1.1398 + return _getColUpperBound(cols(id(c)));
1.1399 + }
1.1400 +
1.1401 + ///\brief Set the upper bound of several columns
1.1402 + ///(i.e variables) at once
1.1403 + ///
1.1404 + ///This magic function takes a container as its argument
1.1405 + ///and applies the function on all of its elements.
1.1406 + ///The upper bound of a variable (column) has to be given by an
1.1407 + ///extended number of type Value, i.e. a finite number of type
1.1408 + ///Value or \ref INF.
1.1409 +#ifdef DOXYGEN
1.1410 + template<class T>
1.1411 + void colUpperBound(T &t, Value value) { return 0;}
1.1412 +#else
1.1413 + template<class T1>
1.1414 + typename enable_if<typename T1::value_type::LpCol,void>::type
1.1415 + colUpperBound(T1 &t, Value value,dummy<0> = 0) {
1.1416 + for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1.1417 + colUpperBound(*i, value);
1.1418 + }
1.1419 + }
1.1420 + template<class T1>
1.1421 + typename enable_if<typename T1::value_type::second_type::LpCol,
1.1422 + void>::type
1.1423 + colUpperBound(T1 &t, Value value,dummy<1> = 1) {
1.1424 + for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1.1425 + colUpperBound(i->second, value);
1.1426 + }
1.1427 + }
1.1428 + template<class T1>
1.1429 + typename enable_if<typename T1::MapIt::Value::LpCol,
1.1430 + void>::type
1.1431 + colUpperBound(T1 &t, Value value,dummy<2> = 2) {
1.1432 + for(typename T1::MapIt i(t); i!=INVALID; ++i){
1.1433 + colUpperBound(*i, value);
1.1434 + }
1.1435 + }
1.1436 +#endif
1.1437 +
1.1438 + /// Set the lower and the upper bounds of a column (i.e a variable)
1.1439 +
1.1440 + /// The lower and the upper bounds of
1.1441 + /// a variable (column) have to be given by an
1.1442 + /// extended number of type Value, i.e. a finite number of type
1.1443 + /// Value, -\ref INF or \ref INF.
1.1444 + void colBounds(Col c, Value lower, Value upper) {
1.1445 + _setColLowerBound(cols(id(c)),lower);
1.1446 + _setColUpperBound(cols(id(c)),upper);
1.1447 + }
1.1448 +
1.1449 + ///\brief Set the lower and the upper bound of several columns
1.1450 + ///(i.e variables) at once
1.1451 + ///
1.1452 + ///This magic function takes a container as its argument
1.1453 + ///and applies the function on all of its elements.
1.1454 + /// The lower and the upper bounds of
1.1455 + /// a variable (column) have to be given by an
1.1456 + /// extended number of type Value, i.e. a finite number of type
1.1457 + /// Value, -\ref INF or \ref INF.
1.1458 +#ifdef DOXYGEN
1.1459 + template<class T>
1.1460 + void colBounds(T &t, Value lower, Value upper) { return 0;}
1.1461 +#else
1.1462 + template<class T2>
1.1463 + typename enable_if<typename T2::value_type::LpCol,void>::type
1.1464 + colBounds(T2 &t, Value lower, Value upper,dummy<0> = 0) {
1.1465 + for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1.1466 + colBounds(*i, lower, upper);
1.1467 + }
1.1468 + }
1.1469 + template<class T2>
1.1470 + typename enable_if<typename T2::value_type::second_type::LpCol, void>::type
1.1471 + colBounds(T2 &t, Value lower, Value upper,dummy<1> = 1) {
1.1472 + for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1.1473 + colBounds(i->second, lower, upper);
1.1474 + }
1.1475 + }
1.1476 + template<class T2>
1.1477 + typename enable_if<typename T2::MapIt::Value::LpCol, void>::type
1.1478 + colBounds(T2 &t, Value lower, Value upper,dummy<2> = 2) {
1.1479 + for(typename T2::MapIt i(t); i!=INVALID; ++i){
1.1480 + colBounds(*i, lower, upper);
1.1481 + }
1.1482 + }
1.1483 +#endif
1.1484 +
1.1485 + /// Set the lower bound of a row (i.e a constraint)
1.1486 +
1.1487 + /// The lower bound of a constraint (row) has to be given by an
1.1488 + /// extended number of type Value, i.e. a finite number of type
1.1489 + /// Value or -\ref INF.
1.1490 + void rowLowerBound(Row r, Value value) {
1.1491 + _setRowLowerBound(rows(id(r)),value);
1.1492 + }
1.1493 +
1.1494 + /// Get the lower bound of a row (i.e a constraint)
1.1495 +
1.1496 + /// This function returns the lower bound for row (constraint) \c c
1.1497 + /// (this might be -\ref INF as well).
1.1498 + ///\return The lower bound for row \c r
1.1499 + Value rowLowerBound(Row r) const {
1.1500 + return _getRowLowerBound(rows(id(r)));
1.1501 + }
1.1502 +
1.1503 + /// Set the upper bound of a row (i.e a constraint)
1.1504 +
1.1505 + /// The upper bound of a constraint (row) has to be given by an
1.1506 + /// extended number of type Value, i.e. a finite number of type
1.1507 + /// Value or -\ref INF.
1.1508 + void rowUpperBound(Row r, Value value) {
1.1509 + _setRowUpperBound(rows(id(r)),value);
1.1510 + }
1.1511 +
1.1512 + /// Get the upper bound of a row (i.e a constraint)
1.1513 +
1.1514 + /// This function returns the upper bound for row (constraint) \c c
1.1515 + /// (this might be -\ref INF as well).
1.1516 + ///\return The upper bound for row \c r
1.1517 + Value rowUpperBound(Row r) const {
1.1518 + return _getRowUpperBound(rows(id(r)));
1.1519 + }
1.1520 +
1.1521 + ///Set an element of the objective function
1.1522 + void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
1.1523 +
1.1524 + ///Get an element of the objective function
1.1525 + Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
1.1526 +
1.1527 + ///Set the objective function
1.1528 +
1.1529 + ///\param e is a linear expression of type \ref Expr.
1.1530 + ///
1.1531 + void obj(const Expr& e) {
1.1532 + _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
1.1533 + ExprIterator(e.comps.end(), cols));
1.1534 + obj_const_comp = *e;
1.1535 + }
1.1536 +
1.1537 + ///Get the objective function
1.1538 +
1.1539 + ///\return the objective function as a linear expression of type
1.1540 + ///Expr.
1.1541 + Expr obj() const {
1.1542 + Expr e;
1.1543 + _getObjCoeffs(InsertIterator(e.comps, cols));
1.1544 + *e = obj_const_comp;
1.1545 + return e;
1.1546 + }
1.1547 +
1.1548 +
1.1549 + ///Set the direction of optimization
1.1550 + void sense(Sense sense) { _setSense(sense); }
1.1551 +
1.1552 + ///Query the direction of the optimization
1.1553 + Sense sense() const {return _getSense(); }
1.1554 +
1.1555 + ///Set the sense to maximization
1.1556 + void max() { _setSense(MAX); }
1.1557 +
1.1558 + ///Set the sense to maximization
1.1559 + void min() { _setSense(MIN); }
1.1560 +
1.1561 + ///Clears the problem
1.1562 + void clear() { _clear(); }
1.1563 +
1.1564 + /// Sets the message level of the solver
1.1565 + void messageLevel(MessageLevel level) { _messageLevel(level); }
1.1566 +
1.1567 + ///@}
1.1568 +
1.1569 + };
1.1570 +
1.1571 + /// Addition
1.1572 +
1.1573 + ///\relates LpBase::Expr
1.1574 + ///
1.1575 + inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
1.1576 + LpBase::Expr tmp(a);
1.1577 + tmp+=b;
1.1578 + return tmp;
1.1579 + }
1.1580 + ///Substraction
1.1581 +
1.1582 + ///\relates LpBase::Expr
1.1583 + ///
1.1584 + inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
1.1585 + LpBase::Expr tmp(a);
1.1586 + tmp-=b;
1.1587 + return tmp;
1.1588 + }
1.1589 + ///Multiply with constant
1.1590 +
1.1591 + ///\relates LpBase::Expr
1.1592 + ///
1.1593 + inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
1.1594 + LpBase::Expr tmp(a);
1.1595 + tmp*=b;
1.1596 + return tmp;
1.1597 + }
1.1598 +
1.1599 + ///Multiply with constant
1.1600 +
1.1601 + ///\relates LpBase::Expr
1.1602 + ///
1.1603 + inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
1.1604 + LpBase::Expr tmp(b);
1.1605 + tmp*=a;
1.1606 + return tmp;
1.1607 + }
1.1608 + ///Divide with constant
1.1609 +
1.1610 + ///\relates LpBase::Expr
1.1611 + ///
1.1612 + inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
1.1613 + LpBase::Expr tmp(a);
1.1614 + tmp/=b;
1.1615 + return tmp;
1.1616 + }
1.1617 +
1.1618 + ///Create constraint
1.1619 +
1.1620 + ///\relates LpBase::Constr
1.1621 + ///
1.1622 + inline LpBase::Constr operator<=(const LpBase::Expr &e,
1.1623 + const LpBase::Expr &f) {
1.1624 + return LpBase::Constr(0, f - e, LpBase::INF);
1.1625 + }
1.1626 +
1.1627 + ///Create constraint
1.1628 +
1.1629 + ///\relates LpBase::Constr
1.1630 + ///
1.1631 + inline LpBase::Constr operator<=(const LpBase::Value &e,
1.1632 + const LpBase::Expr &f) {
1.1633 + return LpBase::Constr(e, f, LpBase::NaN);
1.1634 + }
1.1635 +
1.1636 + ///Create constraint
1.1637 +
1.1638 + ///\relates LpBase::Constr
1.1639 + ///
1.1640 + inline LpBase::Constr operator<=(const LpBase::Expr &e,
1.1641 + const LpBase::Value &f) {
1.1642 + return LpBase::Constr(- LpBase::INF, e, f);
1.1643 + }
1.1644 +
1.1645 + ///Create constraint
1.1646 +
1.1647 + ///\relates LpBase::Constr
1.1648 + ///
1.1649 + inline LpBase::Constr operator>=(const LpBase::Expr &e,
1.1650 + const LpBase::Expr &f) {
1.1651 + return LpBase::Constr(0, e - f, LpBase::INF);
1.1652 + }
1.1653 +
1.1654 +
1.1655 + ///Create constraint
1.1656 +
1.1657 + ///\relates LpBase::Constr
1.1658 + ///
1.1659 + inline LpBase::Constr operator>=(const LpBase::Value &e,
1.1660 + const LpBase::Expr &f) {
1.1661 + return LpBase::Constr(LpBase::NaN, f, e);
1.1662 + }
1.1663 +
1.1664 +
1.1665 + ///Create constraint
1.1666 +
1.1667 + ///\relates LpBase::Constr
1.1668 + ///
1.1669 + inline LpBase::Constr operator>=(const LpBase::Expr &e,
1.1670 + const LpBase::Value &f) {
1.1671 + return LpBase::Constr(f, e, LpBase::INF);
1.1672 + }
1.1673 +
1.1674 + ///Create constraint
1.1675 +
1.1676 + ///\relates LpBase::Constr
1.1677 + ///
1.1678 + inline LpBase::Constr operator==(const LpBase::Expr &e,
1.1679 + const LpBase::Value &f) {
1.1680 + return LpBase::Constr(f, e, f);
1.1681 + }
1.1682 +
1.1683 + ///Create constraint
1.1684 +
1.1685 + ///\relates LpBase::Constr
1.1686 + ///
1.1687 + inline LpBase::Constr operator==(const LpBase::Expr &e,
1.1688 + const LpBase::Expr &f) {
1.1689 + return LpBase::Constr(0, f - e, 0);
1.1690 + }
1.1691 +
1.1692 + ///Create constraint
1.1693 +
1.1694 + ///\relates LpBase::Constr
1.1695 + ///
1.1696 + inline LpBase::Constr operator<=(const LpBase::Value &n,
1.1697 + const LpBase::Constr &c) {
1.1698 + LpBase::Constr tmp(c);
1.1699 + LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.1700 + tmp.lowerBound()=n;
1.1701 + return tmp;
1.1702 + }
1.1703 + ///Create constraint
1.1704 +
1.1705 + ///\relates LpBase::Constr
1.1706 + ///
1.1707 + inline LpBase::Constr operator<=(const LpBase::Constr &c,
1.1708 + const LpBase::Value &n)
1.1709 + {
1.1710 + LpBase::Constr tmp(c);
1.1711 + LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1.1712 + tmp.upperBound()=n;
1.1713 + return tmp;
1.1714 + }
1.1715 +
1.1716 + ///Create constraint
1.1717 +
1.1718 + ///\relates LpBase::Constr
1.1719 + ///
1.1720 + inline LpBase::Constr operator>=(const LpBase::Value &n,
1.1721 + const LpBase::Constr &c) {
1.1722 + LpBase::Constr tmp(c);
1.1723 + LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1.1724 + tmp.upperBound()=n;
1.1725 + return tmp;
1.1726 + }
1.1727 + ///Create constraint
1.1728 +
1.1729 + ///\relates LpBase::Constr
1.1730 + ///
1.1731 + inline LpBase::Constr operator>=(const LpBase::Constr &c,
1.1732 + const LpBase::Value &n)
1.1733 + {
1.1734 + LpBase::Constr tmp(c);
1.1735 + LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.1736 + tmp.lowerBound()=n;
1.1737 + return tmp;
1.1738 + }
1.1739 +
1.1740 + ///Addition
1.1741 +
1.1742 + ///\relates LpBase::DualExpr
1.1743 + ///
1.1744 + inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
1.1745 + const LpBase::DualExpr &b) {
1.1746 + LpBase::DualExpr tmp(a);
1.1747 + tmp+=b;
1.1748 + return tmp;
1.1749 + }
1.1750 + ///Substraction
1.1751 +
1.1752 + ///\relates LpBase::DualExpr
1.1753 + ///
1.1754 + inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
1.1755 + const LpBase::DualExpr &b) {
1.1756 + LpBase::DualExpr tmp(a);
1.1757 + tmp-=b;
1.1758 + return tmp;
1.1759 + }
1.1760 + ///Multiply with constant
1.1761 +
1.1762 + ///\relates LpBase::DualExpr
1.1763 + ///
1.1764 + inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
1.1765 + const LpBase::Value &b) {
1.1766 + LpBase::DualExpr tmp(a);
1.1767 + tmp*=b;
1.1768 + return tmp;
1.1769 + }
1.1770 +
1.1771 + ///Multiply with constant
1.1772 +
1.1773 + ///\relates LpBase::DualExpr
1.1774 + ///
1.1775 + inline LpBase::DualExpr operator*(const LpBase::Value &a,
1.1776 + const LpBase::DualExpr &b) {
1.1777 + LpBase::DualExpr tmp(b);
1.1778 + tmp*=a;
1.1779 + return tmp;
1.1780 + }
1.1781 + ///Divide with constant
1.1782 +
1.1783 + ///\relates LpBase::DualExpr
1.1784 + ///
1.1785 + inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
1.1786 + const LpBase::Value &b) {
1.1787 + LpBase::DualExpr tmp(a);
1.1788 + tmp/=b;
1.1789 + return tmp;
1.1790 + }
1.1791 +
1.1792 + /// \ingroup lp_group
1.1793 + ///
1.1794 + /// \brief Common base class for LP solvers
1.1795 + ///
1.1796 + /// This class is an abstract base class for LP solvers. This class
1.1797 + /// provides a full interface for set and modify an LP problem,
1.1798 + /// solve it and retrieve the solution. You can use one of the
1.1799 + /// descendants as a concrete implementation, or the \c Lp
1.1800 + /// default LP solver. However, if you would like to handle LP
1.1801 + /// solvers as reference or pointer in a generic way, you can use
1.1802 + /// this class directly.
1.1803 + class LpSolver : virtual public LpBase {
1.1804 + public:
1.1805 +
1.1806 + /// The problem types for primal and dual problems
1.1807 + enum ProblemType {
1.1808 + /// = 0. Feasible solution hasn't been found (but may exist).
1.1809 + UNDEFINED = 0,
1.1810 + /// = 1. The problem has no feasible solution.
1.1811 + INFEASIBLE = 1,
1.1812 + /// = 2. Feasible solution found.
1.1813 + FEASIBLE = 2,
1.1814 + /// = 3. Optimal solution exists and found.
1.1815 + OPTIMAL = 3,
1.1816 + /// = 4. The cost function is unbounded.
1.1817 + UNBOUNDED = 4
1.1818 + };
1.1819 +
1.1820 + ///The basis status of variables
1.1821 + enum VarStatus {
1.1822 + /// The variable is in the basis
1.1823 + BASIC,
1.1824 + /// The variable is free, but not basic
1.1825 + FREE,
1.1826 + /// The variable has active lower bound
1.1827 + LOWER,
1.1828 + /// The variable has active upper bound
1.1829 + UPPER,
1.1830 + /// The variable is non-basic and fixed
1.1831 + FIXED
1.1832 + };
1.1833 +
1.1834 + protected:
1.1835 +
1.1836 + virtual SolveExitStatus _solve() = 0;
1.1837 +
1.1838 + virtual Value _getPrimal(int i) const = 0;
1.1839 + virtual Value _getDual(int i) const = 0;
1.1840 +
1.1841 + virtual Value _getPrimalRay(int i) const = 0;
1.1842 + virtual Value _getDualRay(int i) const = 0;
1.1843 +
1.1844 + virtual Value _getPrimalValue() const = 0;
1.1845 +
1.1846 + virtual VarStatus _getColStatus(int i) const = 0;
1.1847 + virtual VarStatus _getRowStatus(int i) const = 0;
1.1848 +
1.1849 + virtual ProblemType _getPrimalType() const = 0;
1.1850 + virtual ProblemType _getDualType() const = 0;
1.1851 +
1.1852 + public:
1.1853 +
1.1854 + ///Allocate a new LP problem instance
1.1855 + virtual LpSolver* newSolver() const = 0;
1.1856 + ///Make a copy of the LP problem
1.1857 + virtual LpSolver* cloneSolver() const = 0;
1.1858 +
1.1859 + ///\name Solve the LP
1.1860 +
1.1861 + ///@{
1.1862 +
1.1863 + ///\e Solve the LP problem at hand
1.1864 + ///
1.1865 + ///\return The result of the optimization procedure. Possible
1.1866 + ///values and their meanings can be found in the documentation of
1.1867 + ///\ref SolveExitStatus.
1.1868 + SolveExitStatus solve() { return _solve(); }
1.1869 +
1.1870 + ///@}
1.1871 +
1.1872 + ///\name Obtain the Solution
1.1873 +
1.1874 + ///@{
1.1875 +
1.1876 + /// The type of the primal problem
1.1877 + ProblemType primalType() const {
1.1878 + return _getPrimalType();
1.1879 + }
1.1880 +
1.1881 + /// The type of the dual problem
1.1882 + ProblemType dualType() const {
1.1883 + return _getDualType();
1.1884 + }
1.1885 +
1.1886 + /// Return the primal value of the column
1.1887 +
1.1888 + /// Return the primal value of the column.
1.1889 + /// \pre The problem is solved.
1.1890 + Value primal(Col c) const { return _getPrimal(cols(id(c))); }
1.1891 +
1.1892 + /// Return the primal value of the expression
1.1893 +
1.1894 + /// Return the primal value of the expression, i.e. the dot
1.1895 + /// product of the primal solution and the expression.
1.1896 + /// \pre The problem is solved.
1.1897 + Value primal(const Expr& e) const {
1.1898 + double res = *e;
1.1899 + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1.1900 + res += *c * primal(c);
1.1901 + }
1.1902 + return res;
1.1903 + }
1.1904 + /// Returns a component of the primal ray
1.1905 +
1.1906 + /// The primal ray is solution of the modified primal problem,
1.1907 + /// where we change each finite bound to 0, and we looking for a
1.1908 + /// negative objective value in case of minimization, and positive
1.1909 + /// objective value for maximization. If there is such solution,
1.1910 + /// that proofs the unsolvability of the dual problem, and if a
1.1911 + /// feasible primal solution exists, then the unboundness of
1.1912 + /// primal problem.
1.1913 + ///
1.1914 + /// \pre The problem is solved and the dual problem is infeasible.
1.1915 + /// \note Some solvers does not provide primal ray calculation
1.1916 + /// functions.
1.1917 + Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
1.1918 +
1.1919 + /// Return the dual value of the row
1.1920 +
1.1921 + /// Return the dual value of the row.
1.1922 + /// \pre The problem is solved.
1.1923 + Value dual(Row r) const { return _getDual(rows(id(r))); }
1.1924 +
1.1925 + /// Return the dual value of the dual expression
1.1926 +
1.1927 + /// Return the dual value of the dual expression, i.e. the dot
1.1928 + /// product of the dual solution and the dual expression.
1.1929 + /// \pre The problem is solved.
1.1930 + Value dual(const DualExpr& e) const {
1.1931 + double res = 0.0;
1.1932 + for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
1.1933 + res += *r * dual(r);
1.1934 + }
1.1935 + return res;
1.1936 + }
1.1937 +
1.1938 + /// Returns a component of the dual ray
1.1939 +
1.1940 + /// The dual ray is solution of the modified primal problem, where
1.1941 + /// we change each finite bound to 0 (i.e. the objective function
1.1942 + /// coefficients in the primal problem), and we looking for a
1.1943 + /// ositive objective value. If there is such solution, that
1.1944 + /// proofs the unsolvability of the primal problem, and if a
1.1945 + /// feasible dual solution exists, then the unboundness of
1.1946 + /// dual problem.
1.1947 + ///
1.1948 + /// \pre The problem is solved and the primal problem is infeasible.
1.1949 + /// \note Some solvers does not provide dual ray calculation
1.1950 + /// functions.
1.1951 + Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
1.1952 +
1.1953 + /// Return the basis status of the column
1.1954 +
1.1955 + /// \see VarStatus
1.1956 + VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
1.1957 +
1.1958 + /// Return the basis status of the row
1.1959 +
1.1960 + /// \see VarStatus
1.1961 + VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
1.1962 +
1.1963 + ///The value of the objective function
1.1964 +
1.1965 + ///\return
1.1966 + ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1.1967 + /// of the primal problem, depending on whether we minimize or maximize.
1.1968 + ///- \ref NaN if no primal solution is found.
1.1969 + ///- The (finite) objective value if an optimal solution is found.
1.1970 + Value primal() const { return _getPrimalValue()+obj_const_comp;}
1.1971 + ///@}
1.1972 +
1.1973 + protected:
1.1974 +
1.1975 + };
1.1976 +
1.1977 +
1.1978 + /// \ingroup lp_group
1.1979 + ///
1.1980 + /// \brief Common base class for MIP solvers
1.1981 + ///
1.1982 + /// This class is an abstract base class for MIP solvers. This class
1.1983 + /// provides a full interface for set and modify an MIP problem,
1.1984 + /// solve it and retrieve the solution. You can use one of the
1.1985 + /// descendants as a concrete implementation, or the \c Lp
1.1986 + /// default MIP solver. However, if you would like to handle MIP
1.1987 + /// solvers as reference or pointer in a generic way, you can use
1.1988 + /// this class directly.
1.1989 + class MipSolver : virtual public LpBase {
1.1990 + public:
1.1991 +
1.1992 + /// The problem types for MIP problems
1.1993 + enum ProblemType {
1.1994 + /// = 0. Feasible solution hasn't been found (but may exist).
1.1995 + UNDEFINED = 0,
1.1996 + /// = 1. The problem has no feasible solution.
1.1997 + INFEASIBLE = 1,
1.1998 + /// = 2. Feasible solution found.
1.1999 + FEASIBLE = 2,
1.2000 + /// = 3. Optimal solution exists and found.
1.2001 + OPTIMAL = 3,
1.2002 + /// = 4. The cost function is unbounded.
1.2003 + ///The Mip or at least the relaxed problem is unbounded.
1.2004 + UNBOUNDED = 4
1.2005 + };
1.2006 +
1.2007 + ///Allocate a new MIP problem instance
1.2008 + virtual MipSolver* newSolver() const = 0;
1.2009 + ///Make a copy of the MIP problem
1.2010 + virtual MipSolver* cloneSolver() const = 0;
1.2011 +
1.2012 + ///\name Solve the MIP
1.2013 +
1.2014 + ///@{
1.2015 +
1.2016 + /// Solve the MIP problem at hand
1.2017 + ///
1.2018 + ///\return The result of the optimization procedure. Possible
1.2019 + ///values and their meanings can be found in the documentation of
1.2020 + ///\ref SolveExitStatus.
1.2021 + SolveExitStatus solve() { return _solve(); }
1.2022 +
1.2023 + ///@}
1.2024 +
1.2025 + ///\name Set Column Type
1.2026 + ///@{
1.2027 +
1.2028 + ///Possible variable (column) types (e.g. real, integer, binary etc.)
1.2029 + enum ColTypes {
1.2030 + /// = 0. Continuous variable (default).
1.2031 + REAL = 0,
1.2032 + /// = 1. Integer variable.
1.2033 + INTEGER = 1
1.2034 + };
1.2035 +
1.2036 + ///Sets the type of the given column to the given type
1.2037 +
1.2038 + ///Sets the type of the given column to the given type.
1.2039 + ///
1.2040 + void colType(Col c, ColTypes col_type) {
1.2041 + _setColType(cols(id(c)),col_type);
1.2042 + }
1.2043 +
1.2044 + ///Gives back the type of the column.
1.2045 +
1.2046 + ///Gives back the type of the column.
1.2047 + ///
1.2048 + ColTypes colType(Col c) const {
1.2049 + return _getColType(cols(id(c)));
1.2050 + }
1.2051 + ///@}
1.2052 +
1.2053 + ///\name Obtain the Solution
1.2054 +
1.2055 + ///@{
1.2056 +
1.2057 + /// The type of the MIP problem
1.2058 + ProblemType type() const {
1.2059 + return _getType();
1.2060 + }
1.2061 +
1.2062 + /// Return the value of the row in the solution
1.2063 +
1.2064 + /// Return the value of the row in the solution.
1.2065 + /// \pre The problem is solved.
1.2066 + Value sol(Col c) const { return _getSol(cols(id(c))); }
1.2067 +
1.2068 + /// Return the value of the expression in the solution
1.2069 +
1.2070 + /// Return the value of the expression in the solution, i.e. the
1.2071 + /// dot product of the solution and the expression.
1.2072 + /// \pre The problem is solved.
1.2073 + Value sol(const Expr& e) const {
1.2074 + double res = *e;
1.2075 + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1.2076 + res += *c * sol(c);
1.2077 + }
1.2078 + return res;
1.2079 + }
1.2080 + ///The value of the objective function
1.2081 +
1.2082 + ///\return
1.2083 + ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1.2084 + /// of the problem, depending on whether we minimize or maximize.
1.2085 + ///- \ref NaN if no primal solution is found.
1.2086 + ///- The (finite) objective value if an optimal solution is found.
1.2087 + Value solValue() const { return _getSolValue()+obj_const_comp;}
1.2088 + ///@}
1.2089 +
1.2090 + protected:
1.2091 +
1.2092 + virtual SolveExitStatus _solve() = 0;
1.2093 + virtual ColTypes _getColType(int col) const = 0;
1.2094 + virtual void _setColType(int col, ColTypes col_type) = 0;
1.2095 + virtual ProblemType _getType() const = 0;
1.2096 + virtual Value _getSol(int i) const = 0;
1.2097 + virtual Value _getSolValue() const = 0;
1.2098 +
1.2099 + };
1.2100 +
1.2101 +
1.2102 +
1.2103 +} //namespace lemon
1.2104 +
1.2105 +#endif //LEMON_LP_BASE_H