1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/lp_base.h Thu Dec 10 17:05:35 2009 +0100
1.3 @@ -0,0 +1,2088 @@
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 void _eraseCol(int col) = 0;
1.950 + virtual void _eraseRow(int row) = 0;
1.951 +
1.952 + virtual void _getColName(int col, std::string& name) const = 0;
1.953 + virtual void _setColName(int col, const std::string& name) = 0;
1.954 + virtual int _colByName(const std::string& name) const = 0;
1.955 +
1.956 + virtual void _getRowName(int row, std::string& name) const = 0;
1.957 + virtual void _setRowName(int row, const std::string& name) = 0;
1.958 + virtual int _rowByName(const std::string& name) const = 0;
1.959 +
1.960 + virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1.961 + virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
1.962 +
1.963 + virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1.964 + virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
1.965 +
1.966 + virtual void _setCoeff(int row, int col, Value value) = 0;
1.967 + virtual Value _getCoeff(int row, int col) const = 0;
1.968 +
1.969 + virtual void _setColLowerBound(int i, Value value) = 0;
1.970 + virtual Value _getColLowerBound(int i) const = 0;
1.971 +
1.972 + virtual void _setColUpperBound(int i, Value value) = 0;
1.973 + virtual Value _getColUpperBound(int i) const = 0;
1.974 +
1.975 + virtual void _setRowLowerBound(int i, Value value) = 0;
1.976 + virtual Value _getRowLowerBound(int i) const = 0;
1.977 +
1.978 + virtual void _setRowUpperBound(int i, Value value) = 0;
1.979 + virtual Value _getRowUpperBound(int i) const = 0;
1.980 +
1.981 + virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
1.982 + virtual void _getObjCoeffs(InsertIterator b) const = 0;
1.983 +
1.984 + virtual void _setObjCoeff(int i, Value obj_coef) = 0;
1.985 + virtual Value _getObjCoeff(int i) const = 0;
1.986 +
1.987 + virtual void _setSense(Sense) = 0;
1.988 + virtual Sense _getSense() const = 0;
1.989 +
1.990 + virtual void _clear() = 0;
1.991 +
1.992 + virtual const char* _solverName() const = 0;
1.993 +
1.994 + virtual void _messageLevel(MessageLevel level) = 0;
1.995 +
1.996 + //Own protected stuff
1.997 +
1.998 + //Constant component of the objective function
1.999 + Value obj_const_comp;
1.1000 +
1.1001 + LpBase() : rows(), cols(), obj_const_comp(0) {}
1.1002 +
1.1003 + public:
1.1004 +
1.1005 + /// Virtual destructor
1.1006 + virtual ~LpBase() {}
1.1007 +
1.1008 + ///Gives back the name of the solver.
1.1009 + const char* solverName() const {return _solverName();}
1.1010 +
1.1011 + ///\name Build Up and Modify the LP
1.1012 +
1.1013 + ///@{
1.1014 +
1.1015 + ///Add a new empty column (i.e a new variable) to the LP
1.1016 + Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
1.1017 +
1.1018 + ///\brief Adds several new columns (i.e variables) at once
1.1019 + ///
1.1020 + ///This magic function takes a container as its argument and fills
1.1021 + ///its elements with new columns (i.e. variables)
1.1022 + ///\param t can be
1.1023 + ///- a standard STL compatible iterable container with
1.1024 + ///\ref Col as its \c values_type like
1.1025 + ///\code
1.1026 + ///std::vector<LpBase::Col>
1.1027 + ///std::list<LpBase::Col>
1.1028 + ///\endcode
1.1029 + ///- a standard STL compatible iterable container with
1.1030 + ///\ref Col as its \c mapped_type like
1.1031 + ///\code
1.1032 + ///std::map<AnyType,LpBase::Col>
1.1033 + ///\endcode
1.1034 + ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.1035 + ///\code
1.1036 + ///ListGraph::NodeMap<LpBase::Col>
1.1037 + ///ListGraph::ArcMap<LpBase::Col>
1.1038 + ///\endcode
1.1039 + ///\return The number of the created column.
1.1040 +#ifdef DOXYGEN
1.1041 + template<class T>
1.1042 + int addColSet(T &t) { return 0;}
1.1043 +#else
1.1044 + template<class T>
1.1045 + typename enable_if<typename T::value_type::LpCol,int>::type
1.1046 + addColSet(T &t,dummy<0> = 0) {
1.1047 + int s=0;
1.1048 + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
1.1049 + return s;
1.1050 + }
1.1051 + template<class T>
1.1052 + typename enable_if<typename T::value_type::second_type::LpCol,
1.1053 + int>::type
1.1054 + addColSet(T &t,dummy<1> = 1) {
1.1055 + int s=0;
1.1056 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1057 + i->second=addCol();
1.1058 + s++;
1.1059 + }
1.1060 + return s;
1.1061 + }
1.1062 + template<class T>
1.1063 + typename enable_if<typename T::MapIt::Value::LpCol,
1.1064 + int>::type
1.1065 + addColSet(T &t,dummy<2> = 2) {
1.1066 + int s=0;
1.1067 + for(typename T::MapIt i(t); i!=INVALID; ++i)
1.1068 + {
1.1069 + i.set(addCol());
1.1070 + s++;
1.1071 + }
1.1072 + return s;
1.1073 + }
1.1074 +#endif
1.1075 +
1.1076 + ///Set a column (i.e a dual constraint) of the LP
1.1077 +
1.1078 + ///\param c is the column to be modified
1.1079 + ///\param e is a dual linear expression (see \ref DualExpr)
1.1080 + ///a better one.
1.1081 + void col(Col c, const DualExpr &e) {
1.1082 + e.simplify();
1.1083 + _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), rows),
1.1084 + ExprIterator(e.comps.end(), rows));
1.1085 + }
1.1086 +
1.1087 + ///Get a column (i.e a dual constraint) of the LP
1.1088 +
1.1089 + ///\param c is the column to get
1.1090 + ///\return the dual expression associated to the column
1.1091 + DualExpr col(Col c) const {
1.1092 + DualExpr e;
1.1093 + _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
1.1094 + return e;
1.1095 + }
1.1096 +
1.1097 + ///Add a new column to the LP
1.1098 +
1.1099 + ///\param e is a dual linear expression (see \ref DualExpr)
1.1100 + ///\param o is the corresponding component of the objective
1.1101 + ///function. It is 0 by default.
1.1102 + ///\return The created column.
1.1103 + Col addCol(const DualExpr &e, Value o = 0) {
1.1104 + Col c=addCol();
1.1105 + col(c,e);
1.1106 + objCoeff(c,o);
1.1107 + return c;
1.1108 + }
1.1109 +
1.1110 + ///Add a new empty row (i.e a new constraint) to the LP
1.1111 +
1.1112 + ///This function adds a new empty row (i.e a new constraint) to the LP.
1.1113 + ///\return The created row
1.1114 + Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
1.1115 +
1.1116 + ///\brief Add several new rows (i.e constraints) at once
1.1117 + ///
1.1118 + ///This magic function takes a container as its argument and fills
1.1119 + ///its elements with new row (i.e. variables)
1.1120 + ///\param t can be
1.1121 + ///- a standard STL compatible iterable container with
1.1122 + ///\ref Row as its \c values_type like
1.1123 + ///\code
1.1124 + ///std::vector<LpBase::Row>
1.1125 + ///std::list<LpBase::Row>
1.1126 + ///\endcode
1.1127 + ///- a standard STL compatible iterable container with
1.1128 + ///\ref Row as its \c mapped_type like
1.1129 + ///\code
1.1130 + ///std::map<AnyType,LpBase::Row>
1.1131 + ///\endcode
1.1132 + ///- an iterable lemon \ref concepts::WriteMap "write map" like
1.1133 + ///\code
1.1134 + ///ListGraph::NodeMap<LpBase::Row>
1.1135 + ///ListGraph::ArcMap<LpBase::Row>
1.1136 + ///\endcode
1.1137 + ///\return The number of rows created.
1.1138 +#ifdef DOXYGEN
1.1139 + template<class T>
1.1140 + int addRowSet(T &t) { return 0;}
1.1141 +#else
1.1142 + template<class T>
1.1143 + typename enable_if<typename T::value_type::LpRow,int>::type
1.1144 + addRowSet(T &t, dummy<0> = 0) {
1.1145 + int s=0;
1.1146 + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1.1147 + return s;
1.1148 + }
1.1149 + template<class T>
1.1150 + typename enable_if<typename T::value_type::second_type::LpRow, int>::type
1.1151 + addRowSet(T &t, dummy<1> = 1) {
1.1152 + int s=0;
1.1153 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1154 + i->second=addRow();
1.1155 + s++;
1.1156 + }
1.1157 + return s;
1.1158 + }
1.1159 + template<class T>
1.1160 + typename enable_if<typename T::MapIt::Value::LpRow, int>::type
1.1161 + addRowSet(T &t, dummy<2> = 2) {
1.1162 + int s=0;
1.1163 + for(typename T::MapIt i(t); i!=INVALID; ++i)
1.1164 + {
1.1165 + i.set(addRow());
1.1166 + s++;
1.1167 + }
1.1168 + return s;
1.1169 + }
1.1170 +#endif
1.1171 +
1.1172 + ///Set a row (i.e a constraint) of the LP
1.1173 +
1.1174 + ///\param r is the row to be modified
1.1175 + ///\param l is lower bound (-\ref INF means no bound)
1.1176 + ///\param e is a linear expression (see \ref Expr)
1.1177 + ///\param u is the upper bound (\ref INF means no bound)
1.1178 + void row(Row r, Value l, const Expr &e, Value u) {
1.1179 + e.simplify();
1.1180 + _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
1.1181 + ExprIterator(e.comps.end(), cols));
1.1182 + _setRowLowerBound(rows(id(r)),l - *e);
1.1183 + _setRowUpperBound(rows(id(r)),u - *e);
1.1184 + }
1.1185 +
1.1186 + ///Set a row (i.e a constraint) of the LP
1.1187 +
1.1188 + ///\param r is the row to be modified
1.1189 + ///\param c is a linear expression (see \ref Constr)
1.1190 + void row(Row r, const Constr &c) {
1.1191 + row(r, c.lowerBounded()?c.lowerBound():-INF,
1.1192 + c.expr(), c.upperBounded()?c.upperBound():INF);
1.1193 + }
1.1194 +
1.1195 +
1.1196 + ///Get a row (i.e a constraint) of the LP
1.1197 +
1.1198 + ///\param r is the row to get
1.1199 + ///\return the expression associated to the row
1.1200 + Expr row(Row r) const {
1.1201 + Expr e;
1.1202 + _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
1.1203 + return e;
1.1204 + }
1.1205 +
1.1206 + ///Add a new row (i.e a new constraint) to the LP
1.1207 +
1.1208 + ///\param l is the lower bound (-\ref INF means no bound)
1.1209 + ///\param e is a linear expression (see \ref Expr)
1.1210 + ///\param u is the upper bound (\ref INF means no bound)
1.1211 + ///\return The created row.
1.1212 + Row addRow(Value l,const Expr &e, Value u) {
1.1213 + Row r=addRow();
1.1214 + row(r,l,e,u);
1.1215 + return r;
1.1216 + }
1.1217 +
1.1218 + ///Add a new row (i.e a new constraint) to the LP
1.1219 +
1.1220 + ///\param c is a linear expression (see \ref Constr)
1.1221 + ///\return The created row.
1.1222 + Row addRow(const Constr &c) {
1.1223 + Row r=addRow();
1.1224 + row(r,c);
1.1225 + return r;
1.1226 + }
1.1227 + ///Erase a column (i.e a variable) from the LP
1.1228 +
1.1229 + ///\param c is the column to be deleted
1.1230 + void erase(Col c) {
1.1231 + _eraseCol(cols(id(c)));
1.1232 + _eraseColId(cols(id(c)));
1.1233 + }
1.1234 + ///Erase a row (i.e a constraint) from the LP
1.1235 +
1.1236 + ///\param r is the row to be deleted
1.1237 + void erase(Row r) {
1.1238 + _eraseRow(rows(id(r)));
1.1239 + _eraseRowId(rows(id(r)));
1.1240 + }
1.1241 +
1.1242 + /// Get the name of a column
1.1243 +
1.1244 + ///\param c is the coresponding column
1.1245 + ///\return The name of the colunm
1.1246 + std::string colName(Col c) const {
1.1247 + std::string name;
1.1248 + _getColName(cols(id(c)), name);
1.1249 + return name;
1.1250 + }
1.1251 +
1.1252 + /// Set the name of a column
1.1253 +
1.1254 + ///\param c is the coresponding column
1.1255 + ///\param name The name to be given
1.1256 + void colName(Col c, const std::string& name) {
1.1257 + _setColName(cols(id(c)), name);
1.1258 + }
1.1259 +
1.1260 + /// Get the column by its name
1.1261 +
1.1262 + ///\param name The name of the column
1.1263 + ///\return the proper column or \c INVALID
1.1264 + Col colByName(const std::string& name) const {
1.1265 + int k = _colByName(name);
1.1266 + return k != -1 ? Col(cols[k]) : Col(INVALID);
1.1267 + }
1.1268 +
1.1269 + /// Get the name of a row
1.1270 +
1.1271 + ///\param r is the coresponding row
1.1272 + ///\return The name of the row
1.1273 + std::string rowName(Row r) const {
1.1274 + std::string name;
1.1275 + _getRowName(rows(id(r)), name);
1.1276 + return name;
1.1277 + }
1.1278 +
1.1279 + /// Set the name of a row
1.1280 +
1.1281 + ///\param r is the coresponding row
1.1282 + ///\param name The name to be given
1.1283 + void rowName(Row r, const std::string& name) {
1.1284 + _setRowName(rows(id(r)), name);
1.1285 + }
1.1286 +
1.1287 + /// Get the row by its name
1.1288 +
1.1289 + ///\param name The name of the row
1.1290 + ///\return the proper row or \c INVALID
1.1291 + Row rowByName(const std::string& name) const {
1.1292 + int k = _rowByName(name);
1.1293 + return k != -1 ? Row(rows[k]) : Row(INVALID);
1.1294 + }
1.1295 +
1.1296 + /// Set an element of the coefficient matrix of the LP
1.1297 +
1.1298 + ///\param r is the row of the element to be modified
1.1299 + ///\param c is the column of the element to be modified
1.1300 + ///\param val is the new value of the coefficient
1.1301 + void coeff(Row r, Col c, Value val) {
1.1302 + _setCoeff(rows(id(r)),cols(id(c)), val);
1.1303 + }
1.1304 +
1.1305 + /// Get an element of the coefficient matrix of the LP
1.1306 +
1.1307 + ///\param r is the row of the element
1.1308 + ///\param c is the column of the element
1.1309 + ///\return the corresponding coefficient
1.1310 + Value coeff(Row r, Col c) const {
1.1311 + return _getCoeff(rows(id(r)),cols(id(c)));
1.1312 + }
1.1313 +
1.1314 + /// Set the lower bound of a column (i.e a variable)
1.1315 +
1.1316 + /// The lower bound of a variable (column) has to be given by an
1.1317 + /// extended number of type Value, i.e. a finite number of type
1.1318 + /// Value or -\ref INF.
1.1319 + void colLowerBound(Col c, Value value) {
1.1320 + _setColLowerBound(cols(id(c)),value);
1.1321 + }
1.1322 +
1.1323 + /// Get the lower bound of a column (i.e a variable)
1.1324 +
1.1325 + /// This function returns the lower bound for column (variable) \c c
1.1326 + /// (this might be -\ref INF as well).
1.1327 + ///\return The lower bound for column \c c
1.1328 + Value colLowerBound(Col c) const {
1.1329 + return _getColLowerBound(cols(id(c)));
1.1330 + }
1.1331 +
1.1332 + ///\brief Set the lower bound of several columns
1.1333 + ///(i.e variables) at once
1.1334 + ///
1.1335 + ///This magic function takes a container as its argument
1.1336 + ///and applies the function on all of its elements.
1.1337 + ///The lower bound of a variable (column) has to be given by an
1.1338 + ///extended number of type Value, i.e. a finite number of type
1.1339 + ///Value or -\ref INF.
1.1340 +#ifdef DOXYGEN
1.1341 + template<class T>
1.1342 + void colLowerBound(T &t, Value value) { return 0;}
1.1343 +#else
1.1344 + template<class T>
1.1345 + typename enable_if<typename T::value_type::LpCol,void>::type
1.1346 + colLowerBound(T &t, Value value,dummy<0> = 0) {
1.1347 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1348 + colLowerBound(*i, value);
1.1349 + }
1.1350 + }
1.1351 + template<class T>
1.1352 + typename enable_if<typename T::value_type::second_type::LpCol,
1.1353 + void>::type
1.1354 + colLowerBound(T &t, Value value,dummy<1> = 1) {
1.1355 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.1356 + colLowerBound(i->second, value);
1.1357 + }
1.1358 + }
1.1359 + template<class T>
1.1360 + typename enable_if<typename T::MapIt::Value::LpCol,
1.1361 + void>::type
1.1362 + colLowerBound(T &t, Value value,dummy<2> = 2) {
1.1363 + for(typename T::MapIt i(t); i!=INVALID; ++i){
1.1364 + colLowerBound(*i, value);
1.1365 + }
1.1366 + }
1.1367 +#endif
1.1368 +
1.1369 + /// Set the upper bound of a column (i.e a variable)
1.1370 +
1.1371 + /// The upper bound of a variable (column) has to be given by an
1.1372 + /// extended number of type Value, i.e. a finite number of type
1.1373 + /// Value or \ref INF.
1.1374 + void colUpperBound(Col c, Value value) {
1.1375 + _setColUpperBound(cols(id(c)),value);
1.1376 + };
1.1377 +
1.1378 + /// Get the upper bound of a column (i.e a variable)
1.1379 +
1.1380 + /// This function returns the upper bound for column (variable) \c c
1.1381 + /// (this might be \ref INF as well).
1.1382 + /// \return The upper bound for column \c c
1.1383 + Value colUpperBound(Col c) const {
1.1384 + return _getColUpperBound(cols(id(c)));
1.1385 + }
1.1386 +
1.1387 + ///\brief Set the upper bound of several columns
1.1388 + ///(i.e variables) at once
1.1389 + ///
1.1390 + ///This magic function takes a container as its argument
1.1391 + ///and applies the function on all of its elements.
1.1392 + ///The upper bound of a variable (column) has to be given by an
1.1393 + ///extended number of type Value, i.e. a finite number of type
1.1394 + ///Value or \ref INF.
1.1395 +#ifdef DOXYGEN
1.1396 + template<class T>
1.1397 + void colUpperBound(T &t, Value value) { return 0;}
1.1398 +#else
1.1399 + template<class T1>
1.1400 + typename enable_if<typename T1::value_type::LpCol,void>::type
1.1401 + colUpperBound(T1 &t, Value value,dummy<0> = 0) {
1.1402 + for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1.1403 + colUpperBound(*i, value);
1.1404 + }
1.1405 + }
1.1406 + template<class T1>
1.1407 + typename enable_if<typename T1::value_type::second_type::LpCol,
1.1408 + void>::type
1.1409 + colUpperBound(T1 &t, Value value,dummy<1> = 1) {
1.1410 + for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1.1411 + colUpperBound(i->second, value);
1.1412 + }
1.1413 + }
1.1414 + template<class T1>
1.1415 + typename enable_if<typename T1::MapIt::Value::LpCol,
1.1416 + void>::type
1.1417 + colUpperBound(T1 &t, Value value,dummy<2> = 2) {
1.1418 + for(typename T1::MapIt i(t); i!=INVALID; ++i){
1.1419 + colUpperBound(*i, value);
1.1420 + }
1.1421 + }
1.1422 +#endif
1.1423 +
1.1424 + /// Set the lower and the upper bounds of a column (i.e a variable)
1.1425 +
1.1426 + /// The lower and the upper bounds of
1.1427 + /// a variable (column) have to be given by an
1.1428 + /// extended number of type Value, i.e. a finite number of type
1.1429 + /// Value, -\ref INF or \ref INF.
1.1430 + void colBounds(Col c, Value lower, Value upper) {
1.1431 + _setColLowerBound(cols(id(c)),lower);
1.1432 + _setColUpperBound(cols(id(c)),upper);
1.1433 + }
1.1434 +
1.1435 + ///\brief Set the lower and the upper bound of several columns
1.1436 + ///(i.e variables) at once
1.1437 + ///
1.1438 + ///This magic function takes a container as its argument
1.1439 + ///and applies the function on all of its elements.
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 +#ifdef DOXYGEN
1.1445 + template<class T>
1.1446 + void colBounds(T &t, Value lower, Value upper) { return 0;}
1.1447 +#else
1.1448 + template<class T2>
1.1449 + typename enable_if<typename T2::value_type::LpCol,void>::type
1.1450 + colBounds(T2 &t, Value lower, Value upper,dummy<0> = 0) {
1.1451 + for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1.1452 + colBounds(*i, lower, upper);
1.1453 + }
1.1454 + }
1.1455 + template<class T2>
1.1456 + typename enable_if<typename T2::value_type::second_type::LpCol, void>::type
1.1457 + colBounds(T2 &t, Value lower, Value upper,dummy<1> = 1) {
1.1458 + for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1.1459 + colBounds(i->second, lower, upper);
1.1460 + }
1.1461 + }
1.1462 + template<class T2>
1.1463 + typename enable_if<typename T2::MapIt::Value::LpCol, void>::type
1.1464 + colBounds(T2 &t, Value lower, Value upper,dummy<2> = 2) {
1.1465 + for(typename T2::MapIt i(t); i!=INVALID; ++i){
1.1466 + colBounds(*i, lower, upper);
1.1467 + }
1.1468 + }
1.1469 +#endif
1.1470 +
1.1471 + /// Set the lower bound of a row (i.e a constraint)
1.1472 +
1.1473 + /// The lower bound of a constraint (row) has to be given by an
1.1474 + /// extended number of type Value, i.e. a finite number of type
1.1475 + /// Value or -\ref INF.
1.1476 + void rowLowerBound(Row r, Value value) {
1.1477 + _setRowLowerBound(rows(id(r)),value);
1.1478 + }
1.1479 +
1.1480 + /// Get the lower bound of a row (i.e a constraint)
1.1481 +
1.1482 + /// This function returns the lower bound for row (constraint) \c c
1.1483 + /// (this might be -\ref INF as well).
1.1484 + ///\return The lower bound for row \c r
1.1485 + Value rowLowerBound(Row r) const {
1.1486 + return _getRowLowerBound(rows(id(r)));
1.1487 + }
1.1488 +
1.1489 + /// Set the upper bound of a row (i.e a constraint)
1.1490 +
1.1491 + /// The upper bound of a constraint (row) has to be given by an
1.1492 + /// extended number of type Value, i.e. a finite number of type
1.1493 + /// Value or -\ref INF.
1.1494 + void rowUpperBound(Row r, Value value) {
1.1495 + _setRowUpperBound(rows(id(r)),value);
1.1496 + }
1.1497 +
1.1498 + /// Get the upper bound of a row (i.e a constraint)
1.1499 +
1.1500 + /// This function returns the upper bound for row (constraint) \c c
1.1501 + /// (this might be -\ref INF as well).
1.1502 + ///\return The upper bound for row \c r
1.1503 + Value rowUpperBound(Row r) const {
1.1504 + return _getRowUpperBound(rows(id(r)));
1.1505 + }
1.1506 +
1.1507 + ///Set an element of the objective function
1.1508 + void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
1.1509 +
1.1510 + ///Get an element of the objective function
1.1511 + Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
1.1512 +
1.1513 + ///Set the objective function
1.1514 +
1.1515 + ///\param e is a linear expression of type \ref Expr.
1.1516 + ///
1.1517 + void obj(const Expr& e) {
1.1518 + _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
1.1519 + ExprIterator(e.comps.end(), cols));
1.1520 + obj_const_comp = *e;
1.1521 + }
1.1522 +
1.1523 + ///Get the objective function
1.1524 +
1.1525 + ///\return the objective function as a linear expression of type
1.1526 + ///Expr.
1.1527 + Expr obj() const {
1.1528 + Expr e;
1.1529 + _getObjCoeffs(InsertIterator(e.comps, cols));
1.1530 + *e = obj_const_comp;
1.1531 + return e;
1.1532 + }
1.1533 +
1.1534 +
1.1535 + ///Set the direction of optimization
1.1536 + void sense(Sense sense) { _setSense(sense); }
1.1537 +
1.1538 + ///Query the direction of the optimization
1.1539 + Sense sense() const {return _getSense(); }
1.1540 +
1.1541 + ///Set the sense to maximization
1.1542 + void max() { _setSense(MAX); }
1.1543 +
1.1544 + ///Set the sense to maximization
1.1545 + void min() { _setSense(MIN); }
1.1546 +
1.1547 + ///Clears the problem
1.1548 + void clear() { _clear(); }
1.1549 +
1.1550 + /// Sets the message level of the solver
1.1551 + void messageLevel(MessageLevel level) { _messageLevel(level); }
1.1552 +
1.1553 + ///@}
1.1554 +
1.1555 + };
1.1556 +
1.1557 + /// Addition
1.1558 +
1.1559 + ///\relates LpBase::Expr
1.1560 + ///
1.1561 + inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
1.1562 + LpBase::Expr tmp(a);
1.1563 + tmp+=b;
1.1564 + return tmp;
1.1565 + }
1.1566 + ///Substraction
1.1567 +
1.1568 + ///\relates LpBase::Expr
1.1569 + ///
1.1570 + inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
1.1571 + LpBase::Expr tmp(a);
1.1572 + tmp-=b;
1.1573 + return tmp;
1.1574 + }
1.1575 + ///Multiply with constant
1.1576 +
1.1577 + ///\relates LpBase::Expr
1.1578 + ///
1.1579 + inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
1.1580 + LpBase::Expr tmp(a);
1.1581 + tmp*=b;
1.1582 + return tmp;
1.1583 + }
1.1584 +
1.1585 + ///Multiply with constant
1.1586 +
1.1587 + ///\relates LpBase::Expr
1.1588 + ///
1.1589 + inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
1.1590 + LpBase::Expr tmp(b);
1.1591 + tmp*=a;
1.1592 + return tmp;
1.1593 + }
1.1594 + ///Divide with constant
1.1595 +
1.1596 + ///\relates LpBase::Expr
1.1597 + ///
1.1598 + inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
1.1599 + LpBase::Expr tmp(a);
1.1600 + tmp/=b;
1.1601 + return tmp;
1.1602 + }
1.1603 +
1.1604 + ///Create constraint
1.1605 +
1.1606 + ///\relates LpBase::Constr
1.1607 + ///
1.1608 + inline LpBase::Constr operator<=(const LpBase::Expr &e,
1.1609 + const LpBase::Expr &f) {
1.1610 + return LpBase::Constr(0, f - e, LpBase::INF);
1.1611 + }
1.1612 +
1.1613 + ///Create constraint
1.1614 +
1.1615 + ///\relates LpBase::Constr
1.1616 + ///
1.1617 + inline LpBase::Constr operator<=(const LpBase::Value &e,
1.1618 + const LpBase::Expr &f) {
1.1619 + return LpBase::Constr(e, f, LpBase::NaN);
1.1620 + }
1.1621 +
1.1622 + ///Create constraint
1.1623 +
1.1624 + ///\relates LpBase::Constr
1.1625 + ///
1.1626 + inline LpBase::Constr operator<=(const LpBase::Expr &e,
1.1627 + const LpBase::Value &f) {
1.1628 + return LpBase::Constr(- LpBase::INF, e, f);
1.1629 + }
1.1630 +
1.1631 + ///Create constraint
1.1632 +
1.1633 + ///\relates LpBase::Constr
1.1634 + ///
1.1635 + inline LpBase::Constr operator>=(const LpBase::Expr &e,
1.1636 + const LpBase::Expr &f) {
1.1637 + return LpBase::Constr(0, e - f, LpBase::INF);
1.1638 + }
1.1639 +
1.1640 +
1.1641 + ///Create constraint
1.1642 +
1.1643 + ///\relates LpBase::Constr
1.1644 + ///
1.1645 + inline LpBase::Constr operator>=(const LpBase::Value &e,
1.1646 + const LpBase::Expr &f) {
1.1647 + return LpBase::Constr(LpBase::NaN, f, e);
1.1648 + }
1.1649 +
1.1650 +
1.1651 + ///Create constraint
1.1652 +
1.1653 + ///\relates LpBase::Constr
1.1654 + ///
1.1655 + inline LpBase::Constr operator>=(const LpBase::Expr &e,
1.1656 + const LpBase::Value &f) {
1.1657 + return LpBase::Constr(f, e, LpBase::INF);
1.1658 + }
1.1659 +
1.1660 + ///Create constraint
1.1661 +
1.1662 + ///\relates LpBase::Constr
1.1663 + ///
1.1664 + inline LpBase::Constr operator==(const LpBase::Expr &e,
1.1665 + const LpBase::Value &f) {
1.1666 + return LpBase::Constr(f, e, f);
1.1667 + }
1.1668 +
1.1669 + ///Create constraint
1.1670 +
1.1671 + ///\relates LpBase::Constr
1.1672 + ///
1.1673 + inline LpBase::Constr operator==(const LpBase::Expr &e,
1.1674 + const LpBase::Expr &f) {
1.1675 + return LpBase::Constr(0, f - e, 0);
1.1676 + }
1.1677 +
1.1678 + ///Create constraint
1.1679 +
1.1680 + ///\relates LpBase::Constr
1.1681 + ///
1.1682 + inline LpBase::Constr operator<=(const LpBase::Value &n,
1.1683 + const LpBase::Constr &c) {
1.1684 + LpBase::Constr tmp(c);
1.1685 + LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.1686 + tmp.lowerBound()=n;
1.1687 + return tmp;
1.1688 + }
1.1689 + ///Create constraint
1.1690 +
1.1691 + ///\relates LpBase::Constr
1.1692 + ///
1.1693 + inline LpBase::Constr operator<=(const LpBase::Constr &c,
1.1694 + const LpBase::Value &n)
1.1695 + {
1.1696 + LpBase::Constr tmp(c);
1.1697 + LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1.1698 + tmp.upperBound()=n;
1.1699 + return tmp;
1.1700 + }
1.1701 +
1.1702 + ///Create constraint
1.1703 +
1.1704 + ///\relates LpBase::Constr
1.1705 + ///
1.1706 + inline LpBase::Constr operator>=(const LpBase::Value &n,
1.1707 + const LpBase::Constr &c) {
1.1708 + LpBase::Constr tmp(c);
1.1709 + LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1.1710 + tmp.upperBound()=n;
1.1711 + return tmp;
1.1712 + }
1.1713 + ///Create constraint
1.1714 +
1.1715 + ///\relates LpBase::Constr
1.1716 + ///
1.1717 + inline LpBase::Constr operator>=(const LpBase::Constr &c,
1.1718 + const LpBase::Value &n)
1.1719 + {
1.1720 + LpBase::Constr tmp(c);
1.1721 + LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1.1722 + tmp.lowerBound()=n;
1.1723 + return tmp;
1.1724 + }
1.1725 +
1.1726 + ///Addition
1.1727 +
1.1728 + ///\relates LpBase::DualExpr
1.1729 + ///
1.1730 + inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
1.1731 + const LpBase::DualExpr &b) {
1.1732 + LpBase::DualExpr tmp(a);
1.1733 + tmp+=b;
1.1734 + return tmp;
1.1735 + }
1.1736 + ///Substraction
1.1737 +
1.1738 + ///\relates LpBase::DualExpr
1.1739 + ///
1.1740 + inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
1.1741 + const LpBase::DualExpr &b) {
1.1742 + LpBase::DualExpr tmp(a);
1.1743 + tmp-=b;
1.1744 + return tmp;
1.1745 + }
1.1746 + ///Multiply with constant
1.1747 +
1.1748 + ///\relates LpBase::DualExpr
1.1749 + ///
1.1750 + inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
1.1751 + const LpBase::Value &b) {
1.1752 + LpBase::DualExpr tmp(a);
1.1753 + tmp*=b;
1.1754 + return tmp;
1.1755 + }
1.1756 +
1.1757 + ///Multiply with constant
1.1758 +
1.1759 + ///\relates LpBase::DualExpr
1.1760 + ///
1.1761 + inline LpBase::DualExpr operator*(const LpBase::Value &a,
1.1762 + const LpBase::DualExpr &b) {
1.1763 + LpBase::DualExpr tmp(b);
1.1764 + tmp*=a;
1.1765 + return tmp;
1.1766 + }
1.1767 + ///Divide with constant
1.1768 +
1.1769 + ///\relates LpBase::DualExpr
1.1770 + ///
1.1771 + inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
1.1772 + const LpBase::Value &b) {
1.1773 + LpBase::DualExpr tmp(a);
1.1774 + tmp/=b;
1.1775 + return tmp;
1.1776 + }
1.1777 +
1.1778 + /// \ingroup lp_group
1.1779 + ///
1.1780 + /// \brief Common base class for LP solvers
1.1781 + ///
1.1782 + /// This class is an abstract base class for LP solvers. This class
1.1783 + /// provides a full interface for set and modify an LP problem,
1.1784 + /// solve it and retrieve the solution. You can use one of the
1.1785 + /// descendants as a concrete implementation, or the \c Lp
1.1786 + /// default LP solver. However, if you would like to handle LP
1.1787 + /// solvers as reference or pointer in a generic way, you can use
1.1788 + /// this class directly.
1.1789 + class LpSolver : virtual public LpBase {
1.1790 + public:
1.1791 +
1.1792 + /// The problem types for primal and dual problems
1.1793 + enum ProblemType {
1.1794 + /// = 0. Feasible solution hasn't been found (but may exist).
1.1795 + UNDEFINED = 0,
1.1796 + /// = 1. The problem has no feasible solution.
1.1797 + INFEASIBLE = 1,
1.1798 + /// = 2. Feasible solution found.
1.1799 + FEASIBLE = 2,
1.1800 + /// = 3. Optimal solution exists and found.
1.1801 + OPTIMAL = 3,
1.1802 + /// = 4. The cost function is unbounded.
1.1803 + UNBOUNDED = 4
1.1804 + };
1.1805 +
1.1806 + ///The basis status of variables
1.1807 + enum VarStatus {
1.1808 + /// The variable is in the basis
1.1809 + BASIC,
1.1810 + /// The variable is free, but not basic
1.1811 + FREE,
1.1812 + /// The variable has active lower bound
1.1813 + LOWER,
1.1814 + /// The variable has active upper bound
1.1815 + UPPER,
1.1816 + /// The variable is non-basic and fixed
1.1817 + FIXED
1.1818 + };
1.1819 +
1.1820 + protected:
1.1821 +
1.1822 + virtual SolveExitStatus _solve() = 0;
1.1823 +
1.1824 + virtual Value _getPrimal(int i) const = 0;
1.1825 + virtual Value _getDual(int i) const = 0;
1.1826 +
1.1827 + virtual Value _getPrimalRay(int i) const = 0;
1.1828 + virtual Value _getDualRay(int i) const = 0;
1.1829 +
1.1830 + virtual Value _getPrimalValue() const = 0;
1.1831 +
1.1832 + virtual VarStatus _getColStatus(int i) const = 0;
1.1833 + virtual VarStatus _getRowStatus(int i) const = 0;
1.1834 +
1.1835 + virtual ProblemType _getPrimalType() const = 0;
1.1836 + virtual ProblemType _getDualType() const = 0;
1.1837 +
1.1838 + public:
1.1839 +
1.1840 + ///Allocate a new LP problem instance
1.1841 + virtual LpSolver* newSolver() const = 0;
1.1842 + ///Make a copy of the LP problem
1.1843 + virtual LpSolver* cloneSolver() const = 0;
1.1844 +
1.1845 + ///\name Solve the LP
1.1846 +
1.1847 + ///@{
1.1848 +
1.1849 + ///\e Solve the LP problem at hand
1.1850 + ///
1.1851 + ///\return The result of the optimization procedure. Possible
1.1852 + ///values and their meanings can be found in the documentation of
1.1853 + ///\ref SolveExitStatus.
1.1854 + SolveExitStatus solve() { return _solve(); }
1.1855 +
1.1856 + ///@}
1.1857 +
1.1858 + ///\name Obtain the Solution
1.1859 +
1.1860 + ///@{
1.1861 +
1.1862 + /// The type of the primal problem
1.1863 + ProblemType primalType() const {
1.1864 + return _getPrimalType();
1.1865 + }
1.1866 +
1.1867 + /// The type of the dual problem
1.1868 + ProblemType dualType() const {
1.1869 + return _getDualType();
1.1870 + }
1.1871 +
1.1872 + /// Return the primal value of the column
1.1873 +
1.1874 + /// Return the primal value of the column.
1.1875 + /// \pre The problem is solved.
1.1876 + Value primal(Col c) const { return _getPrimal(cols(id(c))); }
1.1877 +
1.1878 + /// Return the primal value of the expression
1.1879 +
1.1880 + /// Return the primal value of the expression, i.e. the dot
1.1881 + /// product of the primal solution and the expression.
1.1882 + /// \pre The problem is solved.
1.1883 + Value primal(const Expr& e) const {
1.1884 + double res = *e;
1.1885 + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1.1886 + res += *c * primal(c);
1.1887 + }
1.1888 + return res;
1.1889 + }
1.1890 + /// Returns a component of the primal ray
1.1891 +
1.1892 + /// The primal ray is solution of the modified primal problem,
1.1893 + /// where we change each finite bound to 0, and we looking for a
1.1894 + /// negative objective value in case of minimization, and positive
1.1895 + /// objective value for maximization. If there is such solution,
1.1896 + /// that proofs the unsolvability of the dual problem, and if a
1.1897 + /// feasible primal solution exists, then the unboundness of
1.1898 + /// primal problem.
1.1899 + ///
1.1900 + /// \pre The problem is solved and the dual problem is infeasible.
1.1901 + /// \note Some solvers does not provide primal ray calculation
1.1902 + /// functions.
1.1903 + Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
1.1904 +
1.1905 + /// Return the dual value of the row
1.1906 +
1.1907 + /// Return the dual value of the row.
1.1908 + /// \pre The problem is solved.
1.1909 + Value dual(Row r) const { return _getDual(rows(id(r))); }
1.1910 +
1.1911 + /// Return the dual value of the dual expression
1.1912 +
1.1913 + /// Return the dual value of the dual expression, i.e. the dot
1.1914 + /// product of the dual solution and the dual expression.
1.1915 + /// \pre The problem is solved.
1.1916 + Value dual(const DualExpr& e) const {
1.1917 + double res = 0.0;
1.1918 + for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
1.1919 + res += *r * dual(r);
1.1920 + }
1.1921 + return res;
1.1922 + }
1.1923 +
1.1924 + /// Returns a component of the dual ray
1.1925 +
1.1926 + /// The dual ray is solution of the modified primal problem, where
1.1927 + /// we change each finite bound to 0 (i.e. the objective function
1.1928 + /// coefficients in the primal problem), and we looking for a
1.1929 + /// ositive objective value. If there is such solution, that
1.1930 + /// proofs the unsolvability of the primal problem, and if a
1.1931 + /// feasible dual solution exists, then the unboundness of
1.1932 + /// dual problem.
1.1933 + ///
1.1934 + /// \pre The problem is solved and the primal problem is infeasible.
1.1935 + /// \note Some solvers does not provide dual ray calculation
1.1936 + /// functions.
1.1937 + Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
1.1938 +
1.1939 + /// Return the basis status of the column
1.1940 +
1.1941 + /// \see VarStatus
1.1942 + VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
1.1943 +
1.1944 + /// Return the basis status of the row
1.1945 +
1.1946 + /// \see VarStatus
1.1947 + VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
1.1948 +
1.1949 + ///The value of the objective function
1.1950 +
1.1951 + ///\return
1.1952 + ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1.1953 + /// of the primal problem, depending on whether we minimize or maximize.
1.1954 + ///- \ref NaN if no primal solution is found.
1.1955 + ///- The (finite) objective value if an optimal solution is found.
1.1956 + Value primal() const { return _getPrimalValue()+obj_const_comp;}
1.1957 + ///@}
1.1958 +
1.1959 + protected:
1.1960 +
1.1961 + };
1.1962 +
1.1963 +
1.1964 + /// \ingroup lp_group
1.1965 + ///
1.1966 + /// \brief Common base class for MIP solvers
1.1967 + ///
1.1968 + /// This class is an abstract base class for MIP solvers. This class
1.1969 + /// provides a full interface for set and modify an MIP problem,
1.1970 + /// solve it and retrieve the solution. You can use one of the
1.1971 + /// descendants as a concrete implementation, or the \c Lp
1.1972 + /// default MIP solver. However, if you would like to handle MIP
1.1973 + /// solvers as reference or pointer in a generic way, you can use
1.1974 + /// this class directly.
1.1975 + class MipSolver : virtual public LpBase {
1.1976 + public:
1.1977 +
1.1978 + /// The problem types for MIP problems
1.1979 + enum ProblemType {
1.1980 + /// = 0. Feasible solution hasn't been found (but may exist).
1.1981 + UNDEFINED = 0,
1.1982 + /// = 1. The problem has no feasible solution.
1.1983 + INFEASIBLE = 1,
1.1984 + /// = 2. Feasible solution found.
1.1985 + FEASIBLE = 2,
1.1986 + /// = 3. Optimal solution exists and found.
1.1987 + OPTIMAL = 3,
1.1988 + /// = 4. The cost function is unbounded.
1.1989 + ///The Mip or at least the relaxed problem is unbounded.
1.1990 + UNBOUNDED = 4
1.1991 + };
1.1992 +
1.1993 + ///Allocate a new MIP problem instance
1.1994 + virtual MipSolver* newSolver() const = 0;
1.1995 + ///Make a copy of the MIP problem
1.1996 + virtual MipSolver* cloneSolver() const = 0;
1.1997 +
1.1998 + ///\name Solve the MIP
1.1999 +
1.2000 + ///@{
1.2001 +
1.2002 + /// Solve the MIP problem at hand
1.2003 + ///
1.2004 + ///\return The result of the optimization procedure. Possible
1.2005 + ///values and their meanings can be found in the documentation of
1.2006 + ///\ref SolveExitStatus.
1.2007 + SolveExitStatus solve() { return _solve(); }
1.2008 +
1.2009 + ///@}
1.2010 +
1.2011 + ///\name Set Column Type
1.2012 + ///@{
1.2013 +
1.2014 + ///Possible variable (column) types (e.g. real, integer, binary etc.)
1.2015 + enum ColTypes {
1.2016 + /// = 0. Continuous variable (default).
1.2017 + REAL = 0,
1.2018 + /// = 1. Integer variable.
1.2019 + INTEGER = 1
1.2020 + };
1.2021 +
1.2022 + ///Sets the type of the given column to the given type
1.2023 +
1.2024 + ///Sets the type of the given column to the given type.
1.2025 + ///
1.2026 + void colType(Col c, ColTypes col_type) {
1.2027 + _setColType(cols(id(c)),col_type);
1.2028 + }
1.2029 +
1.2030 + ///Gives back the type of the column.
1.2031 +
1.2032 + ///Gives back the type of the column.
1.2033 + ///
1.2034 + ColTypes colType(Col c) const {
1.2035 + return _getColType(cols(id(c)));
1.2036 + }
1.2037 + ///@}
1.2038 +
1.2039 + ///\name Obtain the Solution
1.2040 +
1.2041 + ///@{
1.2042 +
1.2043 + /// The type of the MIP problem
1.2044 + ProblemType type() const {
1.2045 + return _getType();
1.2046 + }
1.2047 +
1.2048 + /// Return the value of the row in the solution
1.2049 +
1.2050 + /// Return the value of the row in the solution.
1.2051 + /// \pre The problem is solved.
1.2052 + Value sol(Col c) const { return _getSol(cols(id(c))); }
1.2053 +
1.2054 + /// Return the value of the expression in the solution
1.2055 +
1.2056 + /// Return the value of the expression in the solution, i.e. the
1.2057 + /// dot product of the solution and the expression.
1.2058 + /// \pre The problem is solved.
1.2059 + Value sol(const Expr& e) const {
1.2060 + double res = *e;
1.2061 + for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1.2062 + res += *c * sol(c);
1.2063 + }
1.2064 + return res;
1.2065 + }
1.2066 + ///The value of the objective function
1.2067 +
1.2068 + ///\return
1.2069 + ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1.2070 + /// of the problem, depending on whether we minimize or maximize.
1.2071 + ///- \ref NaN if no primal solution is found.
1.2072 + ///- The (finite) objective value if an optimal solution is found.
1.2073 + Value solValue() const { return _getSolValue()+obj_const_comp;}
1.2074 + ///@}
1.2075 +
1.2076 + protected:
1.2077 +
1.2078 + virtual SolveExitStatus _solve() = 0;
1.2079 + virtual ColTypes _getColType(int col) const = 0;
1.2080 + virtual void _setColType(int col, ColTypes col_type) = 0;
1.2081 + virtual ProblemType _getType() const = 0;
1.2082 + virtual Value _getSol(int i) const = 0;
1.2083 + virtual Value _getSolValue() const = 0;
1.2084 +
1.2085 + };
1.2086 +
1.2087 +
1.2088 +
1.2089 +} //namespace lemon
1.2090 +
1.2091 +#endif //LEMON_LP_BASE_H