Redesigned CapacityScaling algorithm with almost the same interface.
The new version does not use the ResidualGraphAdaptor for performance reasons.
Scaling can be enabled and disabled with a parameter of the run() function.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2007
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_LP_BASE_H
20 #define LEMON_LP_BASE_H
28 #include<lemon/error.h>
29 #include<lemon/bits/invalid.h>
30 #include<lemon/bits/utility.h>
31 #include<lemon/bits/lp_id.h>
34 ///\brief The interface of the LP solver interface.
38 /// Function to decide whether a floating point value is finite or not.
40 /// Retruns true if the argument is not infinity, minus infinity or NaN.
41 /// It does the same as the isfinite() function defined by C99.
43 bool isFinite(T value)
45 typedef std::numeric_limits<T> Lim;
46 if (Lim::has_infinity && (value == Lim::infinity() || value ==
48 (Lim::has_quiet_NaN || Lim::has_signaling_NaN) && value != value)
55 ///Common base class for LP solvers
57 ///\todo Much more docs
68 ///Possible outcomes of an LP solving procedure
69 enum SolveExitStatus {
70 ///This means that the problem has been successfully solved: either
71 ///an optimal solution has been found or infeasibility/unboundedness
74 ///Any other case (including the case when some user specified
75 ///limit has been exceeded)
81 ///Feasible solution hasn't been found (but may exist).
83 ///\todo NOTFOUND might be a better name.
86 ///The problem has no feasible solution
88 ///Feasible solution found
90 ///Optimal solution exists and found
92 ///The cost function is unbounded
94 ///\todo Give a feasible solution and an infinite ray (and the
95 ///corresponding bases)
99 ///\e The type of the investigated LP problem
101 ///Primal-dual feasible
102 PRIMAL_DUAL_FEASIBLE = 0,
103 ///Primal feasible dual infeasible
104 PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
105 ///Primal infeasible dual feasible
106 PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
107 ///Primal-dual infeasible
108 PRIMAL_DUAL_INFEASIBLE = 3,
109 ///Could not determine so far
113 ///The floating point type used by the solver
114 typedef double Value;
115 ///The infinity constant
116 static const Value INF;
117 ///The not a number constant
118 static const Value NaN;
120 static inline bool isNaN(const Value& v) { return v!=v; }
126 ///Refer to a column of the LP.
128 ///This type is used to refer to a column of the LP.
130 ///Its value remains valid and correct even after the addition or erase of
133 ///\todo Document what can one do with a Col (INVALID, comparing,
134 ///it is similar to Node/Edge)
138 friend class LpSolverBase;
139 friend class MipSolverBase;
140 explicit Col(int _id) : id(_id) {}
142 typedef Value ExprValue;
143 typedef True LpSolverCol;
145 Col(const Invalid&) : id(-1) {}
146 bool operator< (Col c) const {return id< c.id;}
147 bool operator> (Col c) const {return id> c.id;}
148 bool operator==(Col c) const {return id==c.id;}
149 bool operator!=(Col c) const {return id!=c.id;}
152 class ColIt : public Col {
153 const LpSolverBase *_lp;
156 ColIt(const LpSolverBase &lp) : _lp(&lp)
158 _lp->cols.firstFix(id);
160 ColIt(const Invalid&) : Col(INVALID) {}
163 _lp->cols.nextFix(id);
168 static int id(const Col& col) { return col.id; }
171 ///Refer to a row of the LP.
173 ///This type is used to refer to a row of the LP.
175 ///Its value remains valid and correct even after the addition or erase of
178 ///\todo Document what can one do with a Row (INVALID, comparing,
179 ///it is similar to Node/Edge)
183 friend class LpSolverBase;
184 explicit Row(int _id) : id(_id) {}
186 typedef Value ExprValue;
187 typedef True LpSolverRow;
189 Row(const Invalid&) : id(-1) {}
191 bool operator< (Row c) const {return id< c.id;}
192 bool operator> (Row c) const {return id> c.id;}
193 bool operator==(Row c) const {return id==c.id;}
194 bool operator!=(Row c) const {return id!=c.id;}
197 class RowIt : public Row {
198 const LpSolverBase *_lp;
201 RowIt(const LpSolverBase &lp) : _lp(&lp)
203 _lp->rows.firstFix(id);
205 RowIt(const Invalid&) : Row(INVALID) {}
208 _lp->rows.nextFix(id);
213 static int id(const Row& row) { return row.id; }
217 int _lpId(const Col& c) const {
218 return cols.floatingId(id(c));
221 int _lpId(const Row& r) const {
222 return rows.floatingId(id(r));
225 Col _item(int i, Col) const {
226 return Col(cols.fixId(i));
229 Row _item(int i, Row) const {
230 return Row(rows.fixId(i));
236 ///Linear expression of variables and a constant component
238 ///This data structure stores a linear expression of the variables
239 ///(\ref Col "Col"s) and also has a constant component.
241 ///There are several ways to access and modify the contents of this
243 ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
244 ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can
245 ///read and modify the coefficients like
252 ///or you can also iterate through its elements.
255 ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
258 ///(This code computes the sum of all coefficients).
259 ///- Numbers (<tt>double</tt>'s)
260 ///and variables (\ref Col "Col"s) directly convert to an
261 ///\ref Expr and the usual linear operations are defined, so
264 ///2*v-3.12*(v-w/2)+2
265 ///v*2.1+(3*v+(v*12+w+6)*3)/2
267 ///are valid \ref Expr "Expr"essions.
268 ///The usual assignment operations are also defined.
271 ///e+=2*v-3.12*(v-w/2)+2;
275 ///- The constant member can be set and read by \ref constComp()
278 ///double c=e.constComp();
281 ///\note \ref clear() not only sets all coefficients to 0 but also
282 ///clears the constant components.
286 class Expr : public std::map<Col,Value>
289 typedef LpSolverBase::Col Key;
290 typedef LpSolverBase::Value Value;
293 typedef std::map<Col,Value> Base;
297 typedef True IsLinExpression;
299 Expr() : Base(), const_comp(0) { }
301 Expr(const Key &v) : const_comp(0) {
302 Base::insert(std::make_pair(v, 1));
305 Expr(const Value &v) : const_comp(v) {}
307 void set(const Key &v,const Value &c) {
308 Base::insert(std::make_pair(v, c));
311 Value &constComp() { return const_comp; }
313 const Value &constComp() const { return const_comp; }
315 ///Removes the components with zero coefficient.
317 for (Base::iterator i=Base::begin(); i!=Base::end();) {
320 if ((*i).second==0) Base::erase(i);
325 void simplify() const {
326 const_cast<Expr*>(this)->simplify();
329 ///Removes the coefficients closer to zero than \c tolerance.
330 void simplify(double &tolerance) {
331 for (Base::iterator i=Base::begin(); i!=Base::end();) {
334 if (std::fabs((*i).second)<tolerance) Base::erase(i);
339 ///Sets all coefficients and the constant component to 0.
346 Expr &operator+=(const Expr &e) {
347 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
348 (*this)[j->first]+=j->second;
349 const_comp+=e.const_comp;
353 Expr &operator-=(const Expr &e) {
354 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
355 (*this)[j->first]-=j->second;
356 const_comp-=e.const_comp;
360 Expr &operator*=(const Value &c) {
361 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
367 Expr &operator/=(const Value &c) {
368 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
375 void prettyPrint(std::ostream &os) {
376 //std::fmtflags os.flags();
377 //os.setf(std::ios::showpos);
378 Base::iterator j=Base::begin();
380 os<<j->second<<"*x["<<id(j->first)<<"]";
382 for (; j!=Base::end(); ++j){
385 os<<j->second<<"*x["<<id(j->first)<<"]";
387 //Nem valami korrekt, de nem talaltam meg, hogy kell
388 //os.unsetf(std::ios::showpos);
397 ///This data stucture represents a linear constraint in the LP.
398 ///Basically it is a linear expression with a lower or an upper bound
399 ///(or both). These parts of the constraint can be obtained by the member
400 ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
402 ///There are two ways to construct a constraint.
403 ///- You can set the linear expression and the bounds directly
404 /// by the functions above.
405 ///- The operators <tt>\<=</tt>, <tt>==</tt> and <tt>\>=</tt>
406 /// are defined between expressions, or even between constraints whenever
407 /// it makes sense. Therefore if \c e and \c f are linear expressions and
408 /// \c s and \c t are numbers, then the followings are valid expressions
409 /// and thus they can be used directly e.g. in \ref addRow() whenever
418 ///\warning The validity of a constraint is checked only at run time, so
419 ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw a
420 ///\ref LogicError exception.
424 typedef LpSolverBase::Expr Expr;
425 typedef Expr::Key Key;
426 typedef Expr::Value Value;
433 Constr() : _expr(), _lb(NaN), _ub(NaN) {}
435 Constr(Value lb,const Expr &e,Value ub) :
436 _expr(e), _lb(lb), _ub(ub) {}
438 Constr(const Expr &e,Value ub) :
439 _expr(e), _lb(NaN), _ub(ub) {}
441 Constr(Value lb,const Expr &e) :
442 _expr(e), _lb(lb), _ub(NaN) {}
444 Constr(const Expr &e) :
445 _expr(e), _lb(NaN), _ub(NaN) {}
453 ///Reference to the linear expression
454 Expr &expr() { return _expr; }
455 ///Cont reference to the linear expression
456 const Expr &expr() const { return _expr; }
457 ///Reference to the lower bound.
460 ///- \ref INF "INF": the constraint is lower unbounded.
461 ///- \ref NaN "NaN": lower bound has not been set.
462 ///- finite number: the lower bound
463 Value &lowerBound() { return _lb; }
464 ///The const version of \ref lowerBound()
465 const Value &lowerBound() const { return _lb; }
466 ///Reference to the upper bound.
469 ///- \ref INF "INF": the constraint is upper unbounded.
470 ///- \ref NaN "NaN": upper bound has not been set.
471 ///- finite number: the upper bound
472 Value &upperBound() { return _ub; }
473 ///The const version of \ref upperBound()
474 const Value &upperBound() const { return _ub; }
475 ///Is the constraint lower bounded?
476 bool lowerBounded() const {
477 return isFinite(_lb);
479 ///Is the constraint upper bounded?
480 bool upperBounded() const {
481 return isFinite(_ub);
484 void prettyPrint(std::ostream &os) {
485 if (_lb==-LpSolverBase::INF||isNaN(_lb))
489 _expr.prettyPrint(os);
490 if (_ub==LpSolverBase::INF)
499 ///Linear expression of rows
501 ///This data structure represents a column of the matrix,
502 ///thas is it strores a linear expression of the dual variables
503 ///(\ref Row "Row"s).
505 ///There are several ways to access and modify the contents of this
507 ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
508 ///if \c e is an DualExpr and \c v
509 ///and \c w are of type \ref Row, then you can
510 ///read and modify the coefficients like
517 ///or you can also iterate through its elements.
520 ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
523 ///(This code computes the sum of all coefficients).
524 ///- Numbers (<tt>double</tt>'s)
525 ///and variables (\ref Row "Row"s) directly convert to an
526 ///\ref DualExpr and the usual linear operations are defined, so
530 ///v*2.1+(3*v+(v*12+w)*3)/2
532 ///are valid \ref DualExpr "DualExpr"essions.
533 ///The usual assignment operations are also defined.
536 ///e+=2*v-3.12*(v-w/2);
543 class DualExpr : public std::map<Row,Value>
546 typedef LpSolverBase::Row Key;
547 typedef LpSolverBase::Value Value;
550 typedef std::map<Row,Value> Base;
553 typedef True IsLinExpression;
555 DualExpr() : Base() { }
557 DualExpr(const Key &v) {
558 Base::insert(std::make_pair(v, 1));
561 void set(const Key &v,const Value &c) {
562 Base::insert(std::make_pair(v, c));
565 ///Removes the components with zero coefficient.
567 for (Base::iterator i=Base::begin(); i!=Base::end();) {
570 if ((*i).second==0) Base::erase(i);
575 void simplify() const {
576 const_cast<DualExpr*>(this)->simplify();
579 ///Removes the coefficients closer to zero than \c tolerance.
580 void simplify(double &tolerance) {
581 for (Base::iterator i=Base::begin(); i!=Base::end();) {
584 if (std::fabs((*i).second)<tolerance) Base::erase(i);
589 ///Sets all coefficients to 0.
595 DualExpr &operator+=(const DualExpr &e) {
596 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
597 (*this)[j->first]+=j->second;
601 DualExpr &operator-=(const DualExpr &e) {
602 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
603 (*this)[j->first]-=j->second;
607 DualExpr &operator*=(const Value &c) {
608 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
613 DualExpr &operator/=(const Value &c) {
614 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
623 template <typename _Expr>
624 class MappedOutputIterator {
627 typedef std::insert_iterator<_Expr> Base;
629 typedef std::output_iterator_tag iterator_category;
630 typedef void difference_type;
631 typedef void value_type;
632 typedef void reference;
633 typedef void pointer;
635 MappedOutputIterator(const Base& _base, const LpSolverBase& _lp)
636 : base(_base), lp(_lp) {}
638 MappedOutputIterator& operator*() {
642 MappedOutputIterator& operator=(const std::pair<int, Value>& value) {
643 *base = std::make_pair(lp._item(value.first, typename _Expr::Key()),
648 MappedOutputIterator& operator++() {
653 MappedOutputIterator operator++(int) {
654 MappedOutputIterator tmp(*this);
659 bool operator==(const MappedOutputIterator& it) const {
660 return base == it.base;
663 bool operator!=(const MappedOutputIterator& it) const {
664 return base != it.base;
669 const LpSolverBase& lp;
672 template <typename Expr>
673 class MappedInputIterator {
676 typedef typename Expr::const_iterator Base;
678 typedef typename Base::iterator_category iterator_category;
679 typedef typename Base::difference_type difference_type;
680 typedef const std::pair<int, Value> value_type;
681 typedef value_type reference;
684 pointer(value_type& _value) : value(_value) {}
685 value_type* operator->() { return &value; }
690 MappedInputIterator(const Base& _base, const LpSolverBase& _lp)
691 : base(_base), lp(_lp) {}
693 reference operator*() {
694 return std::make_pair(lp._lpId(base->first), base->second);
697 pointer operator->() {
698 return pointer(operator*());
701 MappedInputIterator& operator++() {
706 MappedInputIterator operator++(int) {
707 MappedInputIterator tmp(*this);
712 bool operator==(const MappedInputIterator& it) const {
713 return base == it.base;
716 bool operator!=(const MappedInputIterator& it) const {
717 return base != it.base;
722 const LpSolverBase& lp;
727 /// STL compatible iterator for lp col
728 typedef MappedInputIterator<Expr> ConstRowIterator;
729 /// STL compatible iterator for lp row
730 typedef MappedInputIterator<DualExpr> ConstColIterator;
732 /// STL compatible iterator for lp col
733 typedef MappedOutputIterator<Expr> RowIterator;
734 /// STL compatible iterator for lp row
735 typedef MappedOutputIterator<DualExpr> ColIterator;
737 //Abstract virtual functions
738 virtual LpSolverBase &_newLp() = 0;
739 virtual LpSolverBase &_copyLp(){
740 ///\todo This should be implemented here, too, when we have
741 ///problem retrieving routines. It can be overriden.
744 LpSolverBase & newlp(_newLp());
746 //return *(LpSolverBase*)0;
749 virtual int _addCol() = 0;
750 virtual int _addRow() = 0;
752 virtual void _eraseCol(int col) = 0;
753 virtual void _eraseRow(int row) = 0;
755 virtual void _getColName(int col, std::string & name) const = 0;
756 virtual void _setColName(int col, const std::string & name) = 0;
757 virtual int _colByName(const std::string& name) const = 0;
759 virtual void _setRowCoeffs(int i, ConstRowIterator b,
760 ConstRowIterator e) = 0;
761 virtual void _getRowCoeffs(int i, RowIterator b) const = 0;
762 virtual void _setColCoeffs(int i, ConstColIterator b,
763 ConstColIterator e) = 0;
764 virtual void _getColCoeffs(int i, ColIterator b) const = 0;
765 virtual void _setCoeff(int row, int col, Value value) = 0;
766 virtual Value _getCoeff(int row, int col) const = 0;
767 virtual void _setColLowerBound(int i, Value value) = 0;
768 virtual Value _getColLowerBound(int i) const = 0;
769 virtual void _setColUpperBound(int i, Value value) = 0;
770 virtual Value _getColUpperBound(int i) const = 0;
771 virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
772 virtual void _getRowBounds(int i, Value &lower, Value &upper) const = 0;
774 virtual void _setObjCoeff(int i, Value obj_coef) = 0;
775 virtual Value _getObjCoeff(int i) const = 0;
776 virtual void _clearObj()=0;
778 virtual SolveExitStatus _solve() = 0;
779 virtual Value _getPrimal(int i) const = 0;
780 virtual Value _getDual(int i) const = 0;
781 virtual Value _getPrimalValue() const = 0;
782 virtual bool _isBasicCol(int i) const = 0;
783 virtual SolutionStatus _getPrimalStatus() const = 0;
784 virtual SolutionStatus _getDualStatus() const = 0;
785 virtual ProblemTypes _getProblemType() const = 0;
787 virtual void _setMax() = 0;
788 virtual void _setMin() = 0;
791 virtual bool _isMax() const = 0;
793 //Own protected stuff
795 //Constant component of the objective function
796 Value obj_const_comp;
801 LpSolverBase() : obj_const_comp(0) {}
804 virtual ~LpSolverBase() {}
806 ///Creates a new LP problem
807 LpSolverBase &newLp() {return _newLp();}
808 ///Makes a copy of the LP problem
809 LpSolverBase ©Lp() {return _copyLp();}
811 ///\name Build up and modify the LP
815 ///Add a new empty column (i.e a new variable) to the LP
816 Col addCol() { Col c; _addCol(); c.id = cols.addId(); return c;}
818 ///\brief Adds several new columns
819 ///(i.e a variables) at once
821 ///This magic function takes a container as its argument
822 ///and fills its elements
823 ///with new columns (i.e. variables)
825 ///- a standard STL compatible iterable container with
826 ///\ref Col as its \c values_type
829 ///std::vector<LpSolverBase::Col>
830 ///std::list<LpSolverBase::Col>
832 ///- a standard STL compatible iterable container with
833 ///\ref Col as its \c mapped_type
836 ///std::map<AnyType,LpSolverBase::Col>
838 ///- an iterable lemon \ref concepts::WriteMap "write map" like
840 ///ListGraph::NodeMap<LpSolverBase::Col>
841 ///ListGraph::EdgeMap<LpSolverBase::Col>
843 ///\return The number of the created column.
846 int addColSet(T &t) { return 0;}
849 typename enable_if<typename T::value_type::LpSolverCol,int>::type
850 addColSet(T &t,dummy<0> = 0) {
852 for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
856 typename enable_if<typename T::value_type::second_type::LpSolverCol,
858 addColSet(T &t,dummy<1> = 1) {
860 for(typename T::iterator i=t.begin();i!=t.end();++i) {
867 typename enable_if<typename T::MapIt::Value::LpSolverCol,
869 addColSet(T &t,dummy<2> = 2) {
871 for(typename T::MapIt i(t); i!=INVALID; ++i)
880 ///Set a column (i.e a dual constraint) of the LP
882 ///\param c is the column to be modified
883 ///\param e is a dual linear expression (see \ref DualExpr)
885 void col(Col c,const DualExpr &e) {
887 _setColCoeffs(_lpId(c), ConstColIterator(e.begin(), *this),
888 ConstColIterator(e.end(), *this));
891 ///Get a column (i.e a dual constraint) of the LP
893 ///\param r is the column to get
894 ///\return the dual expression associated to the column
895 DualExpr col(Col c) const {
897 _getColCoeffs(_lpId(c), ColIterator(std::inserter(e, e.end()), *this));
901 ///Add a new column to the LP
903 ///\param e is a dual linear expression (see \ref DualExpr)
904 ///\param obj is the corresponding component of the objective
905 ///function. It is 0 by default.
906 ///\return The created column.
907 Col addCol(const DualExpr &e, Value o = 0) {
914 ///Add a new empty row (i.e a new constraint) to the LP
916 ///This function adds a new empty row (i.e a new constraint) to the LP.
917 ///\return The created row
918 Row addRow() { Row r; _addRow(); r.id = rows.addId(); return r;}
920 ///\brief Add several new rows
921 ///(i.e a constraints) at once
923 ///This magic function takes a container as its argument
924 ///and fills its elements
925 ///with new row (i.e. variables)
927 ///- a standard STL compatible iterable container with
928 ///\ref Row as its \c values_type
931 ///std::vector<LpSolverBase::Row>
932 ///std::list<LpSolverBase::Row>
934 ///- a standard STL compatible iterable container with
935 ///\ref Row as its \c mapped_type
938 ///std::map<AnyType,LpSolverBase::Row>
940 ///- an iterable lemon \ref concepts::WriteMap "write map" like
942 ///ListGraph::NodeMap<LpSolverBase::Row>
943 ///ListGraph::EdgeMap<LpSolverBase::Row>
945 ///\return The number of rows created.
948 int addRowSet(T &t) { return 0;}
951 typename enable_if<typename T::value_type::LpSolverRow,int>::type
952 addRowSet(T &t,dummy<0> = 0) {
954 for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
958 typename enable_if<typename T::value_type::second_type::LpSolverRow,
960 addRowSet(T &t,dummy<1> = 1) {
962 for(typename T::iterator i=t.begin();i!=t.end();++i) {
969 typename enable_if<typename T::MapIt::Value::LpSolverRow,
971 addRowSet(T &t,dummy<2> = 2) {
973 for(typename T::MapIt i(t); i!=INVALID; ++i)
982 ///Set a row (i.e a constraint) of the LP
984 ///\param r is the row to be modified
985 ///\param l is lower bound (-\ref INF means no bound)
986 ///\param e is a linear expression (see \ref Expr)
987 ///\param u is the upper bound (\ref INF means no bound)
988 ///\bug This is a temporary function. The interface will change to
990 ///\todo Option to control whether a constraint with a single variable is
992 void row(Row r, Value l, const Expr &e, Value u) {
994 _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this),
995 ConstRowIterator(e.end(), *this));
996 _setRowBounds(_lpId(r),l-e.constComp(),u-e.constComp());
999 ///Set a row (i.e a constraint) of the LP
1001 ///\param r is the row to be modified
1002 ///\param c is a linear expression (see \ref Constr)
1003 void row(Row r, const Constr &c) {
1004 row(r, c.lowerBounded()?c.lowerBound():-INF,
1005 c.expr(), c.upperBounded()?c.upperBound():INF);
1009 ///Get a row (i.e a constraint) of the LP
1011 ///\param r is the row to get
1012 ///\return the expression associated to the row
1013 Expr row(Row r) const {
1015 _getRowCoeffs(_lpId(r), RowIterator(std::inserter(e, e.end()), *this));
1019 ///Add a new row (i.e a new constraint) to the LP
1021 ///\param l is the lower bound (-\ref INF means no bound)
1022 ///\param e is a linear expression (see \ref Expr)
1023 ///\param u is the upper bound (\ref INF means no bound)
1024 ///\return The created row.
1025 ///\bug This is a temporary function. The interface will change to
1027 Row addRow(Value l,const Expr &e, Value u) {
1033 ///Add a new row (i.e a new constraint) to the LP
1035 ///\param c is a linear expression (see \ref Constr)
1036 ///\return The created row.
1037 Row addRow(const Constr &c) {
1042 ///Erase a coloumn (i.e a variable) from the LP
1044 ///\param c is the coloumn to be deleted
1045 ///\todo Please check this
1046 void eraseCol(Col c) {
1047 _eraseCol(_lpId(c));
1050 ///Erase a row (i.e a constraint) from the LP
1052 ///\param r is the row to be deleted
1053 ///\todo Please check this
1054 void eraseRow(Row r) {
1055 _eraseRow(_lpId(r));
1059 /// Get the name of a column
1061 ///\param c is the coresponding coloumn
1062 ///\return The name of the colunm
1063 std::string colName(Col c) const {
1065 _getColName(_lpId(c), name);
1069 /// Set the name of a column
1071 ///\param c is the coresponding coloumn
1072 ///\param name The name to be given
1073 void colName(Col c, const std::string& name) {
1074 _setColName(_lpId(c), name);
1077 /// Get the column by its name
1079 ///\param name The name of the column
1080 ///\return the proper column or \c INVALID
1081 Col colByName(const std::string& name) const {
1082 int k = _colByName(name);
1083 return k != -1 ? Col(cols.fixId(k)) : Col(INVALID);
1086 /// Set an element of the coefficient matrix of the LP
1088 ///\param r is the row of the element to be modified
1089 ///\param c is the coloumn of the element to be modified
1090 ///\param val is the new value of the coefficient
1092 void coeff(Row r, Col c, Value val) {
1093 _setCoeff(_lpId(r),_lpId(c), val);
1096 /// Get an element of the coefficient matrix of the LP
1098 ///\param r is the row of the element in question
1099 ///\param c is the coloumn of the element in question
1100 ///\return the corresponding coefficient
1102 Value coeff(Row r, Col c) const {
1103 return _getCoeff(_lpId(r),_lpId(c));
1106 /// Set the lower bound of a column (i.e a variable)
1108 /// The lower bound of a variable (column) has to be given by an
1109 /// extended number of type Value, i.e. a finite number of type
1110 /// Value or -\ref INF.
1111 void colLowerBound(Col c, Value value) {
1112 _setColLowerBound(_lpId(c),value);
1115 /// Get the lower bound of a column (i.e a variable)
1117 /// This function returns the lower bound for column (variable) \t c
1118 /// (this might be -\ref INF as well).
1119 ///\return The lower bound for coloumn \t c
1120 Value colLowerBound(Col c) const {
1121 return _getColLowerBound(_lpId(c));
1124 ///\brief Set the lower bound of several columns
1125 ///(i.e a variables) at once
1127 ///This magic function takes a container as its argument
1128 ///and applies the function on all of its elements.
1129 /// The lower bound of a variable (column) has to be given by an
1130 /// extended number of type Value, i.e. a finite number of type
1131 /// Value or -\ref INF.
1134 void colLowerBound(T &t, Value value) { return 0;}
1137 typename enable_if<typename T::value_type::LpSolverCol,void>::type
1138 colLowerBound(T &t, Value value,dummy<0> = 0) {
1139 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1140 colLowerBound(*i, value);
1144 typename enable_if<typename T::value_type::second_type::LpSolverCol,
1146 colLowerBound(T &t, Value value,dummy<1> = 1) {
1147 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1148 colLowerBound(i->second, value);
1152 typename enable_if<typename T::MapIt::Value::LpSolverCol,
1154 colLowerBound(T &t, Value value,dummy<2> = 2) {
1155 for(typename T::MapIt i(t); i!=INVALID; ++i){
1156 colLowerBound(*i, value);
1161 /// Set the upper bound of a column (i.e a variable)
1163 /// The upper bound of a variable (column) has to be given by an
1164 /// extended number of type Value, i.e. a finite number of type
1165 /// Value or \ref INF.
1166 void colUpperBound(Col c, Value value) {
1167 _setColUpperBound(_lpId(c),value);
1170 /// Get the upper bound of a column (i.e a variable)
1172 /// This function returns the upper bound for column (variable) \t c
1173 /// (this might be \ref INF as well).
1174 ///\return The upper bound for coloumn \t c
1175 Value colUpperBound(Col c) const {
1176 return _getColUpperBound(_lpId(c));
1179 ///\brief Set the upper bound of several columns
1180 ///(i.e a variables) at once
1182 ///This magic function takes a container as its argument
1183 ///and applies the function on all of its elements.
1184 /// The upper bound of a variable (column) has to be given by an
1185 /// extended number of type Value, i.e. a finite number of type
1186 /// Value or \ref INF.
1189 void colUpperBound(T &t, Value value) { return 0;}
1192 typename enable_if<typename T::value_type::LpSolverCol,void>::type
1193 colUpperBound(T &t, Value value,dummy<0> = 0) {
1194 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1195 colUpperBound(*i, value);
1199 typename enable_if<typename T::value_type::second_type::LpSolverCol,
1201 colUpperBound(T &t, Value value,dummy<1> = 1) {
1202 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1203 colUpperBound(i->second, value);
1207 typename enable_if<typename T::MapIt::Value::LpSolverCol,
1209 colUpperBound(T &t, Value value,dummy<2> = 2) {
1210 for(typename T::MapIt i(t); i!=INVALID; ++i){
1211 colUpperBound(*i, value);
1216 /// Set the lower and the upper bounds of a column (i.e a variable)
1218 /// The lower and the upper bounds of
1219 /// a variable (column) have to be given by an
1220 /// extended number of type Value, i.e. a finite number of type
1221 /// Value, -\ref INF or \ref INF.
1222 void colBounds(Col c, Value lower, Value upper) {
1223 _setColLowerBound(_lpId(c),lower);
1224 _setColUpperBound(_lpId(c),upper);
1227 ///\brief Set the lower and the upper bound of several columns
1228 ///(i.e a variables) at once
1230 ///This magic function takes a container as its argument
1231 ///and applies the function on all of its elements.
1232 /// The lower and the upper bounds of
1233 /// a variable (column) have to be given by an
1234 /// extended number of type Value, i.e. a finite number of type
1235 /// Value, -\ref INF or \ref INF.
1238 void colBounds(T &t, Value lower, Value upper) { return 0;}
1241 typename enable_if<typename T::value_type::LpSolverCol,void>::type
1242 colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
1243 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1244 colBounds(*i, lower, upper);
1248 typename enable_if<typename T::value_type::second_type::LpSolverCol,
1250 colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
1251 for(typename T::iterator i=t.begin();i!=t.end();++i) {
1252 colBounds(i->second, lower, upper);
1256 typename enable_if<typename T::MapIt::Value::LpSolverCol,
1258 colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
1259 for(typename T::MapIt i(t); i!=INVALID; ++i){
1260 colBounds(*i, lower, upper);
1266 /// Set the lower and the upper bounds of a row (i.e a constraint)
1268 /// The lower and the upper bound of a constraint (row) have to be
1269 /// given by an extended number of type Value, i.e. a finite
1270 /// number of type Value, -\ref INF or \ref INF. There is no
1271 /// separate function for the lower and the upper bound because
1272 /// that would have been hard to implement for CPLEX.
1273 void rowBounds(Row c, Value lower, Value upper) {
1274 _setRowBounds(_lpId(c),lower, upper);
1277 /// Get the lower and the upper bounds of a row (i.e a constraint)
1279 /// The lower and the upper bound of
1280 /// a constraint (row) are
1281 /// extended numbers of type Value, i.e. finite numbers of type
1282 /// Value, -\ref INF or \ref INF.
1283 /// \todo There is no separate function for the
1284 /// lower and the upper bound because we had problems with the
1285 /// implementation of the setting functions for CPLEX:
1286 /// check out whether this can be done for these functions.
1287 void getRowBounds(Row c, Value &lower, Value &upper) const {
1288 _getRowBounds(_lpId(c),lower, upper);
1291 ///Set an element of the objective function
1292 void objCoeff(Col c, Value v) {_setObjCoeff(_lpId(c),v); };
1294 ///Get an element of the objective function
1295 Value objCoeff(Col c) const { return _getObjCoeff(_lpId(c)); };
1297 ///Set the objective function
1299 ///\param e is a linear expression of type \ref Expr.
1302 for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
1303 objCoeff((*i).first,(*i).second);
1304 obj_const_comp=e.constComp();
1307 ///Get the objective function
1309 ///\return the objective function as a linear expression of type \ref Expr.
1312 for (ColIt it(*this); it != INVALID; ++it) {
1313 double c = objCoeff(it);
1315 e.insert(std::make_pair(it, c));
1323 void max() { _setMax(); }
1325 void min() { _setMin(); }
1327 ///Query function: is this a maximization problem?
1328 bool isMax() const {return _isMax(); }
1330 ///Query function: is this a minimization problem?
1331 bool isMin() const {return !isMax(); }
1336 ///\name Solve the LP
1340 ///\e Solve the LP problem at hand
1342 ///\return The result of the optimization procedure. Possible
1343 ///values and their meanings can be found in the documentation of
1344 ///\ref SolveExitStatus.
1346 ///\todo Which method is used to solve the problem
1347 SolveExitStatus solve() { return _solve(); }
1351 ///\name Obtain the solution
1355 /// The status of the primal problem (the original LP problem)
1356 SolutionStatus primalStatus() const {
1357 return _getPrimalStatus();
1360 /// The status of the dual (of the original LP) problem
1361 SolutionStatus dualStatus() const {
1362 return _getDualStatus();
1365 ///The type of the original LP problem
1366 ProblemTypes problemType() const {
1367 return _getProblemType();
1371 Value primal(Col c) const { return _getPrimal(_lpId(c)); }
1373 Value primal(const Expr& e) const {
1374 double res = e.constComp();
1375 for (std::map<Col, double>::const_iterator it = e.begin();
1376 it != e.end(); ++it) {
1377 res += _getPrimal(_lpId(it->first)) * it->second;
1383 Value dual(Row r) const { return _getDual(_lpId(r)); }
1385 Value dual(const DualExpr& e) const {
1387 for (std::map<Row, double>::const_iterator it = e.begin();
1388 it != e.end(); ++it) {
1389 res += _getPrimal(_lpId(it->first)) * it->second;
1395 bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); }
1400 ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1401 /// of the primal problem, depending on whether we minimize or maximize.
1402 ///- \ref NaN if no primal solution is found.
1403 ///- The (finite) objective value if an optimal solution is found.
1404 Value primalValue() const { return _getPrimalValue()+obj_const_comp;}
1410 /// \ingroup lp_group
1412 /// \brief Common base class for MIP solvers
1413 /// \todo Much more docs
1414 class MipSolverBase : virtual public LpSolverBase{
1417 ///Possible variable (coloumn) types (e.g. real, integer, binary etc.)
1419 ///Continuous variable
1423 ///Unfortunately, cplex 7.5 somewhere writes something like
1424 ///#define INTEGER 'I'
1426 ///\todo No support for other types yet.
1429 ///Sets the type of the given coloumn to the given type
1431 ///Sets the type of the given coloumn to the given type.
1432 void colType(Col c, ColTypes col_type) {
1433 _colType(_lpId(c),col_type);
1436 ///Gives back the type of the column.
1438 ///Gives back the type of the column.
1439 ColTypes colType(Col c) const {
1440 return _colType(_lpId(c));
1443 ///Sets the type of the given Col to integer or remove that property.
1445 ///Sets the type of the given Col to integer or remove that property.
1446 void integer(Col c, bool enable) {
1453 ///Gives back whether the type of the column is integer or not.
1455 ///Gives back the type of the column.
1456 ///\return true if the column has integer type and false if not.
1457 bool integer(Col c) const {
1458 return (colType(c)==INT);
1461 /// The status of the MIP problem
1462 SolutionStatus mipStatus() const {
1463 return _getMipStatus();
1468 virtual ColTypes _colType(int col) const = 0;
1469 virtual void _colType(int col, ColTypes col_type) = 0;
1470 virtual SolutionStatus _getMipStatus() const = 0;
1474 ///\relates LpSolverBase::Expr
1476 inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
1477 const LpSolverBase::Expr &b)
1479 LpSolverBase::Expr tmp(a);
1485 ///\relates LpSolverBase::Expr
1487 inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
1488 const LpSolverBase::Expr &b)
1490 LpSolverBase::Expr tmp(a);
1496 ///\relates LpSolverBase::Expr
1498 inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
1499 const LpSolverBase::Value &b)
1501 LpSolverBase::Expr tmp(a);
1508 ///\relates LpSolverBase::Expr
1510 inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
1511 const LpSolverBase::Expr &b)
1513 LpSolverBase::Expr tmp(b);
1519 ///\relates LpSolverBase::Expr
1521 inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
1522 const LpSolverBase::Value &b)
1524 LpSolverBase::Expr tmp(a);
1531 ///\relates LpSolverBase::Constr
1533 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1534 const LpSolverBase::Expr &f)
1536 return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
1541 ///\relates LpSolverBase::Constr
1543 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
1544 const LpSolverBase::Expr &f)
1546 return LpSolverBase::Constr(e,f);
1551 ///\relates LpSolverBase::Constr
1553 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1554 const LpSolverBase::Value &f)
1556 return LpSolverBase::Constr(e,f);
1561 ///\relates LpSolverBase::Constr
1563 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1564 const LpSolverBase::Expr &f)
1566 return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
1572 ///\relates LpSolverBase::Constr
1574 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
1575 const LpSolverBase::Expr &f)
1577 return LpSolverBase::Constr(f,e);
1583 ///\relates LpSolverBase::Constr
1585 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1586 const LpSolverBase::Value &f)
1588 return LpSolverBase::Constr(f,e);
1593 ///\relates LpSolverBase::Constr
1595 inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1596 const LpSolverBase::Value &f)
1598 return LpSolverBase::Constr(f,e,f);
1603 ///\relates LpSolverBase::Constr
1605 inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1606 const LpSolverBase::Expr &f)
1608 return LpSolverBase::Constr(0,e-f,0);
1613 ///\relates LpSolverBase::Constr
1615 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
1616 const LpSolverBase::Constr&c)
1618 LpSolverBase::Constr tmp(c);
1619 ///\todo Create an own exception type.
1620 if(!LpSolverBase::isNaN(tmp.lowerBound())) throw LogicError();
1621 else tmp.lowerBound()=n;
1626 ///\relates LpSolverBase::Constr
1628 inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
1629 const LpSolverBase::Value &n)
1631 LpSolverBase::Constr tmp(c);
1632 ///\todo Create an own exception type.
1633 if(!LpSolverBase::isNaN(tmp.upperBound())) throw LogicError();
1634 else tmp.upperBound()=n;
1640 ///\relates LpSolverBase::Constr
1642 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
1643 const LpSolverBase::Constr&c)
1645 LpSolverBase::Constr tmp(c);
1646 ///\todo Create an own exception type.
1647 if(!LpSolverBase::isNaN(tmp.upperBound())) throw LogicError();
1648 else tmp.upperBound()=n;
1653 ///\relates LpSolverBase::Constr
1655 inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
1656 const LpSolverBase::Value &n)
1658 LpSolverBase::Constr tmp(c);
1659 ///\todo Create an own exception type.
1660 if(!LpSolverBase::isNaN(tmp.lowerBound())) throw LogicError();
1661 else tmp.lowerBound()=n;
1667 ///\relates LpSolverBase::DualExpr
1669 inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
1670 const LpSolverBase::DualExpr &b)
1672 LpSolverBase::DualExpr tmp(a);
1678 ///\relates LpSolverBase::DualExpr
1680 inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
1681 const LpSolverBase::DualExpr &b)
1683 LpSolverBase::DualExpr tmp(a);
1689 ///\relates LpSolverBase::DualExpr
1691 inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
1692 const LpSolverBase::Value &b)
1694 LpSolverBase::DualExpr tmp(a);
1701 ///\relates LpSolverBase::DualExpr
1703 inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
1704 const LpSolverBase::DualExpr &b)
1706 LpSolverBase::DualExpr tmp(b);
1712 ///\relates LpSolverBase::DualExpr
1714 inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
1715 const LpSolverBase::Value &b)
1717 LpSolverBase::DualExpr tmp(a);
1725 #endif //LEMON_LP_BASE_H