Obsolte test removed.
2 * src/lemon/lp_base.h - Part of LEMON, a generic C++ optimization library
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5 * (Egervary Combinatorial Optimization Research Group, EGRES).
7 * Permission to use, modify and distribute this software is granted
8 * provided that this copyright notice appears in all copies. For
9 * precise terms see the accompanying LICENSE file.
11 * This software is provided "AS IS" with no warranty of any kind,
12 * express or implied, and with no claim as to its suitability for any
17 #ifndef LEMON_LP_BASE_H
18 #define LEMON_LP_BASE_H
25 #include<lemon/utility.h>
26 #include<lemon/error.h>
27 #include<lemon/invalid.h>
29 //#include"lin_expr.h"
32 ///\brief The interface of the LP solver interface.
33 ///\ingroup gen_opt_group
36 ///Internal data structure to convert floating id's to fix one's
38 ///\todo This might be implemented to be also usable in other places.
41 std::vector<int> index;
42 std::vector<int> cross;
45 _FixId() : first_free(-1) {};
46 ///Convert a floating id to a fix one
48 ///\param n is a floating id
49 ///\return the corresponding fix id
50 int fixId(int n) {return cross[n];}
51 ///Convert a fix id to a floating one
53 ///\param n is a fix id
54 ///\return the corresponding floating id
55 int floatingId(int n) { return index[n];}
56 ///Add a new floating id.
58 ///\param n is a floating id
59 ///\return the fix id of the new value
60 ///\todo Multiple additions should also be handled.
63 if(n>=int(cross.size())) {
66 cross[n]=index.size();
71 int next=index[first_free];
77 ///\todo Create an own exception type.
78 else throw LogicError(); //floatingId-s must form a continuous range;
82 ///\param n is a fix id
89 for(int i=fl+1;i<int(cross.size());++i) {
95 ///An upper bound on the largest fix id.
97 ///\todo Do we need this?
99 std::size_t maxFixId() { return cross.size()-1; }
103 ///Common base class for LP solvers
105 ///\todo Much more docs
106 ///\ingroup gen_opt_group
112 enum SolveExitStatus {
120 enum SolutionStatus {
121 ///Feasible solution has'n been found (but may exist).
123 ///\todo NOTFOUND might be a better name.
126 ///The problem has no feasible solution
128 ///Feasible solution found
130 ///Optimal solution exists and found
132 ///The cost function is unbounded
134 ///\todo Give a feasible solution and an infinite ray (and the
135 ///corresponding bases)
139 ///The floating point type used by the solver
140 typedef double Value;
141 ///The infinity constant
142 static const Value INF;
143 ///The not a number constant
144 static const Value NaN;
146 ///Refer to a column of the LP.
148 ///This type is used to refer to a column of the LP.
150 ///Its value remains valid and correct even after the addition or erase of
153 ///\todo Document what can one do with a Col (INVALID, comparing,
154 ///it is similar to Node/Edge)
158 friend class LpSolverBase;
160 typedef Value ExprValue;
161 typedef True LpSolverCol;
163 Col(const Invalid&) : id(-1) {}
164 bool operator<(Col c) const {return id<c.id;}
165 bool operator==(Col c) const {return id==c.id;}
166 bool operator!=(Col c) const {return id==c.id;}
169 ///Refer to a row of the LP.
171 ///This type is used to refer to a row of the LP.
173 ///Its value remains valid and correct even after the addition or erase of
176 ///\todo Document what can one do with a Row (INVALID, comparing,
177 ///it is similar to Node/Edge)
181 friend class LpSolverBase;
183 typedef Value ExprValue;
184 typedef True LpSolverRow;
186 Row(const Invalid&) : id(-1) {}
187 typedef True LpSolverRow;
188 bool operator<(Row c) const {return id<c.id;}
189 bool operator==(Row c) const {return id==c.id;}
190 bool operator!=(Row c) const {return id==c.id;}
193 ///Linear expression of variables and a constant component
195 ///This data structure strores a linear expression of the variables
196 ///(\ref Col "Col"s) and also has a constant component.
198 ///There are several ways to access and modify the contents of this
200 ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
201 ///if \c e is an Expr and \c v and \c w are of type \ref Col then you can
202 ///read and modify the coefficients like
209 ///or you can also iterate through its elements.
212 ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
215 ///(This code computes the sum of all coefficients).
216 ///- Numbers (<tt>double</tt>'s)
217 ///and variables (\ref Col "Col"s) directly convert to an
218 ///\ref Expr and the usual linear operations are defined so
221 ///2*v-3.12*(v-w/2)+2
222 ///v*2.1+(3*v+(v*12+w+6)*3)/2
224 ///are valid \ref Expr "Expr"essions.
225 ///The usual assignment operations are also defined.
228 ///e+=2*v-3.12*(v-w/2)+2;
232 ///- The constant member can be set and read by \ref constComp()
235 ///double c=e.constComp();
238 ///\note \ref clear() not only sets all coefficients to 0 but also
239 ///clears the constant components.
243 class Expr : public std::map<Col,Value>
246 typedef LpSolverBase::Col Key;
247 typedef LpSolverBase::Value Value;
250 typedef std::map<Col,Value> Base;
254 typedef True IsLinExpression;
256 Expr() : Base(), const_comp(0) { }
258 Expr(const Key &v) : const_comp(0) {
259 Base::insert(std::make_pair(v, 1));
262 Expr(const Value &v) : const_comp(v) {}
264 void set(const Key &v,const Value &c) {
265 Base::insert(std::make_pair(v, c));
268 Value &constComp() { return const_comp; }
270 const Value &constComp() const { return const_comp; }
272 ///Removes the components with zero coefficient.
274 for (Base::iterator i=Base::begin(); i!=Base::end();) {
277 if ((*i).second==0) Base::erase(i);
282 ///Sets all coefficients and the constant component to 0.
289 Expr &operator+=(const Expr &e) {
290 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
291 (*this)[j->first]+=j->second;
292 ///\todo it might be speeded up using "hints"
293 const_comp+=e.const_comp;
297 Expr &operator-=(const Expr &e) {
298 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
299 (*this)[j->first]-=j->second;
300 const_comp-=e.const_comp;
304 Expr &operator*=(const Value &c) {
305 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
311 Expr &operator/=(const Value &c) {
312 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
321 ///\todo document please
326 typedef LpSolverBase::Expr Expr;
327 typedef Expr::Key Key;
328 typedef Expr::Value Value;
330 static const Value INF;
331 static const Value NaN;
332 // static const Value INF=0;
333 // static const Value NaN=1;
340 Constr() : _expr(), _lb(NaN), _ub(NaN) {}
342 Constr(Value lb,const Expr &e,Value ub) :
343 _expr(e), _lb(lb), _ub(ub) {}
345 Constr(const Expr &e,Value ub) :
346 _expr(e), _lb(NaN), _ub(ub) {}
348 Constr(Value lb,const Expr &e) :
349 _expr(e), _lb(lb), _ub(NaN) {}
351 Constr(const Expr &e) :
352 _expr(e), _lb(NaN), _ub(NaN) {}
360 Expr &expr() { return _expr; }
362 const Expr &expr() const { return _expr; }
364 Value &lowerBound() { return _lb; }
366 const Value &lowerBound() const { return _lb; }
368 Value &upperBound() { return _ub; }
370 const Value &upperBound() const { return _ub; }
372 bool lowerBounded() const {
374 return isfinite(_lb);
377 bool upperBounded() const {
379 return isfinite(_ub);
388 //Abstract virtual functions
389 virtual int _addCol() = 0;
390 virtual int _addRow() = 0;
391 virtual void _setRowCoeffs(int i,
394 Value const * values ) = 0;
395 virtual void _setColCoeffs(int i,
398 Value const * values ) = 0;
399 virtual void _setColLowerBound(int i, Value value) = 0;
400 virtual void _setColUpperBound(int i, Value value) = 0;
401 virtual void _setRowLowerBound(int i, Value value) = 0;
402 virtual void _setRowUpperBound(int i, Value value) = 0;
403 virtual void _setObjCoeff(int i, Value obj_coef) = 0;
404 virtual SolveExitStatus _solve() = 0;
405 virtual Value _getPrimal(int i) = 0;
406 virtual Value _getPrimalValue() = 0;
407 virtual SolutionStatus _getPrimalStatus() = 0;
408 virtual void _setMax() = 0;
409 virtual void _setMin() = 0;
411 //Own protected stuff
413 //Constant component of the objective function
414 Value obj_const_comp;
418 ///\bug Unimplemented
424 LpSolverBase() : obj_const_comp(0) {}
427 virtual ~LpSolverBase() {}
429 ///\name Build up and modify of the LP
433 ///Add a new empty column (i.e a new variable) to the LP
434 Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
436 ///\brief Adds several new columns
437 ///(i.e a variables) at once
439 ///This magic function takes a container as its argument
440 ///and fills its elements
441 ///with new columns (i.e. variables)
443 ///- a standard STL compatible iterable container with
444 ///\ref Col as its \c values_type
447 ///std::vector<LpSolverBase::Col>
448 ///std::list<LpSolverBase::Col>
450 ///- a standard STL compatible iterable container with
451 ///\ref Col as its \c mapped_type
454 ///std::map<AnyStatus,LpSolverBase::Col>
456 ///- an iterable lemon \ref concept::WriteMap "write map" like
458 ///ListGraph::NodeMap<LpSolverBase::Col>
459 ///ListGraph::EdgeMap<LpSolverBase::Col>
461 ///\return The number of the created column.
462 ///\bug Iterable nodemap hasn't been implemented yet.
465 int addColSet(T &t) { return 0;}
468 typename enable_if<typename T::value_type::LpSolverCol,int>::type
469 addColSet(T &t,dummy<0> = 0) {
471 for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
475 typename enable_if<typename T::value_type::second_type::LpSolverCol,
477 addColSet(T &t,dummy<1> = 1) {
479 for(typename T::iterator i=t.begin();i!=t.end();++i) {
486 typename enable_if<typename T::ValueSet::value_type::LpSolverCol,
488 addColSet(T &t,dummy<2> = 2) {
489 ///\bug <tt>return addColSet(t.valueSet());</tt> should also work.
491 for(typename T::ValueSet::iterator i=t.valueSet().begin();
492 i!=t.valueSet().end();
502 ///Add a new empty row (i.e a new constaint) to the LP
504 ///This function adds a new empty row (i.e a new constaint) to the LP.
505 ///\return The created row
506 Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
508 ///Set a row (i.e a constaint) of the LP
510 ///\param r is the row to be modified
511 ///\param l is lower bound (-\ref INF means no bound)
512 ///\param e is a linear expression (see \ref Expr)
513 ///\param u is the upper bound (\ref INF means no bound)
514 ///\bug This is a temportary function. The interface will change to
516 ///\todo Option to control whether a constraint with a single variable is
518 void setRow(Row r, Value l,const Expr &e, Value u) {
519 std::vector<int> indices;
520 std::vector<Value> values;
521 indices.push_back(0);
523 for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
524 if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
525 indices.push_back(cols.floatingId((*i).first.id));
526 values.push_back((*i).second);
528 _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
529 &indices[0],&values[0]);
530 _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
531 _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
534 ///Set a row (i.e a constaint) of the LP
536 ///\param r is the row to be modified
537 ///\param c is a linear expression (see \ref Constr)
538 void setRow(Row r, const Constr &c) {
540 c.lowerBounded()?c.lowerBound():-INF,
542 c.upperBounded()?c.upperBound():INF);
545 ///Add a new row (i.e a new constaint) to the LP
547 ///\param l is the lower bound (-\ref INF means no bound)
548 ///\param e is a linear expression (see \ref Expr)
549 ///\param u is the upper bound (\ref INF means no bound)
550 ///\return The created row.
551 ///\bug This is a temportary function. The interface will change to
553 Row addRow(Value l,const Expr &e, Value u) {
559 ///Add a new row (i.e a new constaint) to the LP
561 ///\param c is a linear expression (see \ref Constr)
562 ///\return The created row.
563 Row addRow(const Constr &c) {
569 /// Set the lower bound of a column (i.e a variable)
571 /// The upper bound of a variable (column) has to be given by an
572 /// extended number of type Value, i.e. a finite number of type
573 /// Value or -\ref INF.
574 void colLowerBound(Col c, Value value) {
575 _setColLowerBound(cols.floatingId(c.id),value);
577 /// Set the upper bound of a column (i.e a variable)
579 /// The upper bound of a variable (column) has to be given by an
580 /// extended number of type Value, i.e. a finite number of type
581 /// Value or \ref INF.
582 void colUpperBound(Col c, Value value) {
583 _setColUpperBound(cols.floatingId(c.id),value);
585 /// Set the lower and the upper bounds of a column (i.e a variable)
587 /// The lower and the upper bounds of
588 /// a variable (column) have to be given by an
589 /// extended number of type Value, i.e. a finite number of type
590 /// Value, -\ref INF or \ref INF.
591 void colBounds(Col c, Value lower, Value upper) {
592 _setColLowerBound(cols.floatingId(c.id),lower);
593 _setColUpperBound(cols.floatingId(c.id),upper);
596 /// Set the lower bound of a row (i.e a constraint)
598 /// The lower bound of a linear expression (row) has to be given by an
599 /// extended number of type Value, i.e. a finite number of type
600 /// Value or -\ref INF.
601 void rowLowerBound(Row r, Value value) {
602 _setRowLowerBound(rows.floatingId(r.id),value);
604 /// Set the upper bound of a row (i.e a constraint)
606 /// The upper bound of a linear expression (row) has to be given by an
607 /// extended number of type Value, i.e. a finite number of type
608 /// Value or \ref INF.
609 void rowUpperBound(Row r, Value value) {
610 _setRowUpperBound(rows.floatingId(r.id),value);
612 /// Set the lower and the upper bounds of a row (i.e a variable)
614 /// The lower and the upper bounds of
615 /// a constraint (row) have to be given by an
616 /// extended number of type Value, i.e. a finite number of type
617 /// Value, -\ref INF or \ref INF.
618 void rowBounds(Row c, Value lower, Value upper) {
619 _setRowLowerBound(rows.floatingId(c.id),lower);
620 _setRowUpperBound(rows.floatingId(c.id),upper);
623 ///Set an element of the objective function
624 void objCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
625 ///Set the objective function
627 ///\param e is a linear expression of type \ref Expr.
628 ///\bug The previous objective function is not cleared!
629 void setObj(Expr e) {
631 for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
632 objCoeff((*i).first,(*i).second);
633 obj_const_comp=e.constComp();
637 void max() { _setMax(); }
639 void min() { _setMin(); }
645 ///\name Solve the LP
650 SolveExitStatus solve() { return _solve(); }
654 ///\name Obtain the solution
659 SolutionStatus primalStatus() {
660 return _getPrimalStatus();
664 Value primal(Col c) { return _getPrimal(cols.floatingId(c.id)); }
669 ///- \ref INF or -\ref INF means either infeasibility or unboundedness
670 /// of the primal problem, depending on whether we minimize or maximize.
671 ///- \ref NAN if no primal solution is found.
672 ///- The (finite) objective value if an optimal solution is found.
673 Value primalValue() { return _getPrimalValue()+obj_const_comp;}
680 ///\relates LpSolverBase::Expr
682 inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
683 const LpSolverBase::Expr &b)
685 LpSolverBase::Expr tmp(a);
686 tmp+=b; ///\todo Don't STL have some special 'merge' algorithm?
691 ///\relates LpSolverBase::Expr
693 inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
694 const LpSolverBase::Expr &b)
696 LpSolverBase::Expr tmp(a);
697 tmp-=b; ///\todo Don't STL have some special 'merge' algorithm?
702 ///\relates LpSolverBase::Expr
704 inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
705 const LpSolverBase::Value &b)
707 LpSolverBase::Expr tmp(a);
708 tmp*=b; ///\todo Don't STL have some special 'merge' algorithm?
714 ///\relates LpSolverBase::Expr
716 inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
717 const LpSolverBase::Expr &b)
719 LpSolverBase::Expr tmp(b);
720 tmp*=a; ///\todo Don't STL have some special 'merge' algorithm?
725 ///\relates LpSolverBase::Expr
727 inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
728 const LpSolverBase::Value &b)
730 LpSolverBase::Expr tmp(a);
731 tmp/=b; ///\todo Don't STL have some special 'merge' algorithm?
737 ///\relates LpSolverBase::Constr
739 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
740 const LpSolverBase::Expr &f)
742 return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
747 ///\relates LpSolverBase::Constr
749 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
750 const LpSolverBase::Expr &f)
752 return LpSolverBase::Constr(e,f);
757 ///\relates LpSolverBase::Constr
759 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
760 const LpSolverBase::Value &f)
762 return LpSolverBase::Constr(e,f);
767 ///\relates LpSolverBase::Constr
769 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
770 const LpSolverBase::Expr &f)
772 return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
778 ///\relates LpSolverBase::Constr
780 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
781 const LpSolverBase::Expr &f)
783 return LpSolverBase::Constr(f,e);
789 ///\relates LpSolverBase::Constr
791 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
792 const LpSolverBase::Value &f)
794 return LpSolverBase::Constr(f,e);
799 ///\relates LpSolverBase::Constr
801 inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
802 const LpSolverBase::Expr &f)
804 return LpSolverBase::Constr(0,e-f,0);
809 ///\relates LpSolverBase::Constr
811 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
812 const LpSolverBase::Constr&c)
814 LpSolverBase::Constr tmp(c);
815 ///\todo Create an own exception type.
816 if(!isnan(tmp.lowerBound())) throw LogicError();
817 else tmp.lowerBound()=n;
822 ///\relates LpSolverBase::Constr
824 inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
825 const LpSolverBase::Value &n)
827 LpSolverBase::Constr tmp(c);
828 ///\todo Create an own exception type.
829 if(!isnan(tmp.upperBound())) throw LogicError();
830 else tmp.upperBound()=n;
836 ///\relates LpSolverBase::Constr
838 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
839 const LpSolverBase::Constr&c)
841 LpSolverBase::Constr tmp(c);
842 ///\todo Create an own exception type.
843 if(!isnan(tmp.upperBound())) throw LogicError();
844 else tmp.upperBound()=n;
849 ///\relates LpSolverBase::Constr
851 inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
852 const LpSolverBase::Value &n)
854 LpSolverBase::Constr tmp(c);
855 ///\todo Create an own exception type.
856 if(!isnan(tmp.lowerBound())) throw LogicError();
857 else tmp.lowerBound()=n;
864 #endif //LEMON_LP_BASE_H