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.
35 ///Internal data structure to convert floating id's to fix one's
37 ///\todo This might by implemented to be usable in other places.
40 std::vector<int> index;
41 std::vector<int> cross;
44 _FixId() : first_free(-1) {};
45 ///Convert a floating id to a fix one
47 ///\param n is a floating id
48 ///\return the corresponding fix id
49 int fixId(int n) {return cross[n];}
50 ///Convert a fix id to a floating one
52 ///\param n is a fix id
53 ///\return the corresponding floating id
54 int floatingId(int n) { return index[n];}
55 ///Add a new floating id.
57 ///\param n is a floating id
58 ///\return the fix id of the new value
59 ///\todo Multiple additions should also be handled.
62 if(n>=int(cross.size())) {
65 cross[n]=index.size();
70 int next=index[first_free];
76 ///\todo Create an own exception type.
77 else throw LogicError(); //floatingId-s must form a continuous range;
81 ///\param n is a fix id
88 for(int i=fl+1;i<int(cross.size());++i) {
94 ///An upper bound on the largest fix id.
96 ///\todo Do we need this?
98 std::size_t maxFixId() { return cross.size()-1; }
102 ///Common base class for LP solvers
119 ///The floating point type used by the solver
120 typedef double Value;
121 ///The infinity constant
122 static const Value INF;
123 ///The not a number constant
124 static const Value NaN;
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;
140 typedef Value ExprValue;
141 typedef True LpSolverCol;
143 Col(const Invalid&) : id(-1) {}
144 bool operator<(Col c) const {return id<c.id;}
145 bool operator==(Col c) const {return id==c.id;}
146 bool operator!=(Col c) const {return id==c.id;}
149 ///Refer to a row of the LP.
151 ///This type is used to refer to a row of the LP.
153 ///Its value remains valid and correct even after the addition or erase of
156 ///\todo Document what can one do with a Row (INVALID, comparing,
157 ///it is similar to Node/Edge)
161 friend class LpSolverBase;
163 typedef Value ExprValue;
164 typedef True LpSolverRow;
166 Row(const Invalid&) : id(-1) {}
167 typedef True LpSolverRow;
168 bool operator<(Row c) const {return id<c.id;}
169 bool operator==(Row c) const {return id==c.id;}
170 bool operator!=(Row c) const {return id==c.id;}
174 // typedef SparseLinExpr<Col> Expr;
175 class Expr : public std::map<Col,Value>
178 typedef LpSolverBase::Col Key;
179 typedef LpSolverBase::Value Value;
182 typedef std::map<Col,Value> Base;
186 typedef True IsLinExpression;
188 Expr() : Base(), const_comp(0) { }
190 Expr(const Key &v) : const_comp(0) {
191 Base::insert(std::make_pair(v, 1));
194 Expr(const Value &v) : const_comp(v) {}
196 void set(const Key &v,const Value &c) {
197 Base::insert(std::make_pair(v, c));
200 Value &constComp() { return const_comp; }
202 const Value &constComp() const { return const_comp; }
204 ///Removes the components with zero coefficient.
206 for (Base::iterator i=Base::begin(); i!=Base::end();) {
209 if ((*i).second==0) Base::erase(i);
214 ///Sets all coefficients and the constant component to 0.
221 Expr &operator+=(const Expr &e) {
222 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
223 (*this)[j->first]+=j->second;
224 ///\todo it might be speeded up using "hints"
225 const_comp+=e.const_comp;
229 Expr &operator-=(const Expr &e) {
230 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
231 (*this)[j->first]-=j->second;
232 const_comp-=e.const_comp;
236 Expr &operator*=(const Value &c) {
237 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
243 Expr &operator/=(const Value &c) {
244 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
252 //typedef LinConstr<Expr> Constr;
256 typedef LpSolverBase::Expr Expr;
257 typedef Expr::Key Key;
258 typedef Expr::Value Value;
260 static const Value INF;
261 static const Value NaN;
262 // static const Value INF=0;
263 // static const Value NaN=1;
270 Constr() : _expr(), _lb(NaN), _ub(NaN) {}
272 Constr(Value lb,const Expr &e,Value ub) :
273 _expr(e), _lb(lb), _ub(ub) {}
275 Constr(const Expr &e,Value ub) :
276 _expr(e), _lb(NaN), _ub(ub) {}
278 Constr(Value lb,const Expr &e) :
279 _expr(e), _lb(lb), _ub(NaN) {}
281 Constr(const Expr &e) :
282 _expr(e), _lb(NaN), _ub(NaN) {}
290 Expr &expr() { return _expr; }
292 const Expr &expr() const { return _expr; }
294 Value &lowerBound() { return _lb; }
296 const Value &lowerBound() const { return _lb; }
298 Value &upperBound() { return _ub; }
300 const Value &upperBound() const { return _ub; }
302 bool lowerBounded() const { return std::isfinite(_lb); }
304 bool upperBounded() const { return std::isfinite(_ub); }
313 virtual int _addCol() = 0;
315 virtual int _addRow() = 0;
318 /// \warning Arrays are indexed from 1 (datum at index 0 is ignored)
320 virtual void _setRowCoeffs(int i,
323 Value const * values ) = 0;
326 /// \warning Arrays are indexed from 1 (datum at index 0 is ignored)
328 virtual void _setColCoeffs(int i,
331 Value const * values ) = 0;
335 /// The lower bound of a variable (column) have to be given by an
336 /// extended number of type Value, i.e. a finite number of type
337 /// Value or -\ref INF.
338 virtual void _setColLowerBound(int i, Value value) = 0;
341 /// The upper bound of a variable (column) have to be given by an
342 /// extended number of type Value, i.e. a finite number of type
343 /// Value or \ref INF.
344 virtual void _setColUpperBound(int i, Value value) = 0;
347 /// The lower bound of a linear expression (row) have to be given by an
348 /// extended number of type Value, i.e. a finite number of type
349 /// Value or -\ref INF.
350 virtual void _setRowLowerBound(int i, Value value) = 0;
353 /// The upper bound of a linear expression (row) have to be given by an
354 /// extended number of type Value, i.e. a finite number of type
355 /// Value or \ref INF.
356 virtual void _setRowUpperBound(int i, Value value) = 0;
359 virtual void _setObjCoeff(int i, Value obj_coef) = 0;
363 ///\bug Wrong interface
365 virtual SolutionType _solve() = 0;
369 ///\bug Wrong interface
371 virtual Value _getSolution(int i) = 0;
374 ///\bug unimplemented!!!!
380 virtual ~LpSolverBase() {}
382 ///\name Building up and modification of the LP
386 ///Add a new empty column (i.e a new variable) to the LP
387 Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
389 ///\brief Fill the elements of a container with newly created columns
390 ///(i.e a new variables)
392 ///This magic function takes a container as its argument
393 ///and fills its elements
394 ///with new columns (i.e. variables)
396 ///- a standard STL compatible iterable container with
397 ///\ref Col as its \c values_type
400 ///std::vector<LpSolverBase::Col>
401 ///std::list<LpSolverBase::Col>
403 ///- a standard STL compatible iterable container with
404 ///\ref Col as its \c mapped_type
407 ///std::map<AnyType,LpSolverBase::Col>
409 ///- an iterable lemon \ref concept::WriteMap "write map" like
411 ///ListGraph::NodeMap<LpSolverBase::Col>
412 ///ListGraph::EdgeMap<LpSolverBase::Col>
414 ///\return The number of the created column.
415 ///\bug Iterable nodemap hasn't been implemented yet.
418 int addColSet(T &t) { return 0;}
421 typename enable_if<typename T::value_type::LpSolverCol,int>::type
422 addColSet(T &t,dummy<0> = 0) {
424 for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
428 typename enable_if<typename T::value_type::second_type::LpSolverCol,
430 addColSet(T &t,dummy<1> = 1) {
432 for(typename T::iterator i=t.begin();i!=t.end();++i) {
439 typename enable_if<typename T::ValueSet::value_type::LpSolverCol,
441 addColSet(T &t,dummy<2> = 2) {
442 ///\bug <tt>return addColSet(t.valueSet());</tt> should also work.
444 for(typename T::ValueSet::iterator i=t.valueSet().begin();
445 i!=t.valueSet().end();
455 ///Add a new empty row (i.e a new constaint) to the LP
457 ///This function adds a new empty row (i.e a new constaint) to the LP.
458 ///\return The created row
459 Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
461 ///Set a row (i.e a constaint) of the LP
463 ///\param r is the row to be modified
464 ///\param l is lower bound (-\ref INF means no bound)
465 ///\param e is a linear expression (see \ref Expr)
466 ///\param u is the upper bound (\ref INF means no bound)
467 ///\bug This is a temportary function. The interface will change to
469 void setRow(Row r, Value l,const Expr &e, Value u) {
470 std::vector<int> indices;
471 std::vector<Value> values;
472 indices.push_back(0);
474 for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
475 if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
476 indices.push_back(cols.floatingId((*i).first.id));
477 values.push_back((*i).second);
479 _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
480 &indices[0],&values[0]);
481 _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
482 _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
485 ///Set a row (i.e a constaint) of the LP
487 ///\param r is the row to be modified
488 ///\param c is a linear expression (see \ref Constr)
489 void setRow(Row r, const Constr &c) {
491 c.lowerBounded()?c.lowerBound():-INF,
493 c.upperBounded()?c.upperBound():INF);
496 ///Add a new row (i.e a new constaint) to the LP
498 ///\param l is the lower bound (-\ref INF means no bound)
499 ///\param e is a linear expression (see \ref Expr)
500 ///\param u is the upper bound (\ref INF means no bound)
501 ///\return The created row.
502 ///\bug This is a temportary function. The interface will change to
504 Row addRow(Value l,const Expr &e, Value u) {
510 ///Add a new row (i.e a new constaint) to the LP
512 ///\param c is a linear expression (see \ref Constr)
513 ///\return The created row.
514 Row addRow(const Constr &c) {
520 /// Set the lower bound of a column (i.e a variable)
522 /// The upper bound of a variable (column) have to be given by an
523 /// extended number of type Value, i.e. a finite number of type
524 /// Value or -\ref INF.
525 virtual void setColLowerBound(Col c, Value value) {
526 _setColLowerBound(cols.floatingId(c.id),value);
528 /// Set the upper bound of a column (i.e a variable)
530 /// The upper bound of a variable (column) have to be given by an
531 /// extended number of type Value, i.e. a finite number of type
532 /// Value or \ref INF.
533 virtual void setColUpperBound(Col c, Value value) {
534 _setColUpperBound(cols.floatingId(c.id),value);
536 /// Set the lower bound of a row (i.e a constraint)
538 /// The lower bound of a linear expression (row) have to be given by an
539 /// extended number of type Value, i.e. a finite number of type
540 /// Value or -\ref INF.
541 virtual void setRowLowerBound(Row r, Value value) {
542 _setRowLowerBound(rows.floatingId(r.id),value);
544 /// Set the upper bound of a row (i.e a constraint)
546 /// The upper bound of a linear expression (row) have to be given by an
547 /// extended number of type Value, i.e. a finite number of type
548 /// Value or \ref INF.
549 virtual void setRowUpperBound(Row r, Value value) {
550 _setRowUpperBound(rows.floatingId(r.id),value);
552 ///Set an element of the objective function
553 void setObjCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
554 ///Set the objective function
556 ///\param e is a linear expression of type \ref Expr.
557 ///\todo What to do with the constant component?
558 void setObj(Expr e) {
560 for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
561 setObjCoeff((*i).first,(*i).second);
567 ///\name Solving the LP
572 SolutionType solve() { return _solve(); }
576 ///\name Obtaining the solution LP
581 Value solution(Col c) { return _getSolution(cols.floatingId(c.id)); }
589 ///\relates LpSolverBase::Expr
591 inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
592 const LpSolverBase::Expr &b)
594 LpSolverBase::Expr tmp(a);
595 tmp+=b; ///\todo Don't STL have some special 'merge' algorithm?
600 ///\relates LpSolverBase::Expr
602 inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
603 const LpSolverBase::Expr &b)
605 LpSolverBase::Expr tmp(a);
606 tmp-=b; ///\todo Don't STL have some special 'merge' algorithm?
611 ///\relates LpSolverBase::Expr
613 inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
614 const LpSolverBase::Value &b)
616 LpSolverBase::Expr tmp(a);
617 tmp*=b; ///\todo Don't STL have some special 'merge' algorithm?
623 ///\relates LpSolverBase::Expr
625 inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
626 const LpSolverBase::Expr &b)
628 LpSolverBase::Expr tmp(b);
629 tmp*=a; ///\todo Don't STL have some special 'merge' algorithm?
634 ///\relates LpSolverBase::Expr
636 inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
637 const LpSolverBase::Value &b)
639 LpSolverBase::Expr tmp(a);
640 tmp/=b; ///\todo Don't STL have some special 'merge' algorithm?
646 ///\relates LpSolverBase::Constr
648 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
649 const LpSolverBase::Expr &f)
651 return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
656 ///\relates LpSolverBase::Constr
658 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
659 const LpSolverBase::Expr &f)
661 return LpSolverBase::Constr(e,f);
666 ///\relates LpSolverBase::Constr
668 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
669 const LpSolverBase::Value &f)
671 return LpSolverBase::Constr(e,f);
676 ///\relates LpSolverBase::Constr
678 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
679 const LpSolverBase::Expr &f)
681 return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
687 ///\relates LpSolverBase::Constr
689 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
690 const LpSolverBase::Expr &f)
692 return LpSolverBase::Constr(f,e);
698 ///\relates LpSolverBase::Constr
700 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
701 const LpSolverBase::Value &f)
703 return LpSolverBase::Constr(f,e);
708 ///\relates LpSolverBase::Constr
710 inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
711 const LpSolverBase::Expr &f)
713 return LpSolverBase::Constr(0,e-f,0);
718 ///\relates LpSolverBase::Constr
720 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
721 const LpSolverBase::Constr&c)
723 LpSolverBase::Constr tmp(c);
724 ///\todo Create an own exception type.
725 if(!isnan(tmp.lowerBound())) throw LogicError();
726 else tmp.lowerBound()=n;
731 ///\relates LpSolverBase::Constr
733 inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
734 const LpSolverBase::Value &n)
736 LpSolverBase::Constr tmp(c);
737 ///\todo Create an own exception type.
738 if(!isnan(tmp.upperBound())) throw LogicError();
739 else tmp.upperBound()=n;
745 ///\relates LpSolverBase::Constr
747 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
748 const LpSolverBase::Constr&c)
750 LpSolverBase::Constr tmp(c);
751 ///\todo Create an own exception type.
752 if(!isnan(tmp.upperBound())) throw LogicError();
753 else tmp.upperBound()=n;
758 ///\relates LpSolverBase::Constr
760 inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
761 const LpSolverBase::Value &n)
763 LpSolverBase::Constr tmp(c);
764 ///\todo Create an own exception type.
765 if(!isnan(tmp.lowerBound())) throw LogicError();
766 else tmp.lowerBound()=n;
773 #endif //LEMON_LP_BASE_H