COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/lp_base.h @ 1993:2115143eceea

Last change on this file since 1993:2115143eceea was 1993:2115143eceea, checked in by Balazs Dezso, 14 years ago

utility, invalid and traits moved to bits

File size: 39.0 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2006
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
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.
12 *
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
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_LP_BASE_H
20#define LEMON_LP_BASE_H
21
22#include<vector>
23#include<map>
24#include<limits>
25#include<cmath>
26
27#include<lemon/bits/utility.h>
28#include<lemon/error.h>
29#include<lemon/bits/invalid.h>
30
31///\file
32///\brief The interface of the LP solver interface.
33///\ingroup gen_opt_group
34namespace lemon {
35 
36  ///Internal data structure to convert floating id's to fix one's
37   
38  ///\todo This might be implemented to be also usable in other places.
39  class _FixId
40  {
41  protected:
42    std::vector<int> index;
43    std::vector<int> cross;
44    int first_free;
45  public:
46    _FixId() : first_free(-1) {};
47    ///Convert a floating id to a fix one
48
49    ///\param n is a floating id
50    ///\return the corresponding fix id
51    int fixId(int n) const {return cross[n];}
52    ///Convert a fix id to a floating one
53
54    ///\param n is a fix id
55    ///\return the corresponding floating id
56    int floatingId(int n) const { return index[n];}
57    ///Add a new floating id.
58
59    ///\param n is a floating id
60    ///\return the fix id of the new value
61    ///\todo Multiple additions should also be handled.
62    int insert(int n)
63    {
64      if(n>=int(cross.size())) {
65        cross.resize(n+1);
66        if(first_free==-1) {
67          cross[n]=index.size();
68          index.push_back(n);
69        }
70        else {
71          cross[n]=first_free;
72          int next=index[first_free];
73          index[first_free]=n;
74          first_free=next;
75        }
76        return cross[n];
77      }
78      ///\todo Create an own exception type.
79      else throw LogicError(); //floatingId-s must form a continuous range;
80    }
81    ///Remove a fix id.
82
83    ///\param n is a fix id
84    ///
85    void erase(int n)
86    {
87      int fl=index[n];
88      index[n]=first_free;
89      first_free=n;
90      for(int i=fl+1;i<int(cross.size());++i) {
91        cross[i-1]=cross[i];
92        index[cross[i]]--;
93      }
94      cross.pop_back();
95    }
96    ///An upper bound on the largest fix id.
97
98    ///\todo Do we need this?
99    ///
100    std::size_t maxFixId() { return cross.size()-1; }
101 
102  };
103   
104  ///Common base class for LP solvers
105 
106  ///\todo Much more docs
107  ///\ingroup gen_opt_group
108  class LpSolverBase {
109
110  public:
111
112    ///Possible outcomes of an LP solving procedure
113    enum SolveExitStatus {
114      ///This means that the problem has been successfully solved: either
115      ///an optimal solution has been found or infeasibility/unboundedness
116      ///has been proved.
117      SOLVED = 0,
118      ///Any other case (including the case when some user specified limit has been exceeded)
119      UNSOLVED = 1
120    };
121     
122      ///\e
123    enum SolutionStatus {
124      ///Feasible solution has'n been found (but may exist).
125
126      ///\todo NOTFOUND might be a better name.
127      ///
128      UNDEFINED = 0,
129      ///The problem has no feasible solution
130      INFEASIBLE = 1,
131      ///Feasible solution found
132      FEASIBLE = 2,
133      ///Optimal solution exists and found
134      OPTIMAL = 3,
135      ///The cost function is unbounded
136
137      ///\todo Give a feasible solution and an infinite ray (and the
138      ///corresponding bases)
139      INFINITE = 4
140    };
141
142    ///\e The type of the investigated LP problem
143    enum ProblemTypes {
144      ///Primal-dual feasible
145      PRIMAL_DUAL_FEASIBLE = 0,
146      ///Primal feasible dual infeasible
147      PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1,
148      ///Primal infeasible dual feasible
149      PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2,
150      ///Primal-dual infeasible
151      PRIMAL_DUAL_INFEASIBLE = 3,
152      ///Could not determine so far
153      UNKNOWN = 4
154    };
155
156    ///The floating point type used by the solver
157    typedef double Value;
158    ///The infinity constant
159    static const Value INF;
160    ///The not a number constant
161    static const Value NaN;
162   
163    ///Refer to a column of the LP.
164
165    ///This type is used to refer to a column of the LP.
166    ///
167    ///Its value remains valid and correct even after the addition or erase of
168    ///other columns.
169    ///
170    ///\todo Document what can one do with a Col (INVALID, comparing,
171    ///it is similar to Node/Edge)
172    class Col {
173    protected:
174      int id;
175      friend class LpSolverBase;
176    public:
177      typedef Value ExprValue;
178      typedef True LpSolverCol;
179      Col() {}
180      Col(const Invalid&) : id(-1) {}
181      bool operator< (Col c) const  {return id< c.id;}
182      bool operator> (Col c) const  {return id> c.id;}
183      bool operator==(Col c) const  {return id==c.id;}
184      bool operator!=(Col c) const  {return id!=c.id;}
185    };
186
187    ///Refer to a row of the LP.
188
189    ///This type is used to refer to a row of the LP.
190    ///
191    ///Its value remains valid and correct even after the addition or erase of
192    ///other rows.
193    ///
194    ///\todo Document what can one do with a Row (INVALID, comparing,
195    ///it is similar to Node/Edge)
196    class Row {
197    protected:
198      int id;
199      friend class LpSolverBase;
200    public:
201      typedef Value ExprValue;
202      typedef True LpSolverRow;
203      Row() {}
204      Row(const Invalid&) : id(-1) {}
205
206      bool operator< (Row c) const  {return id< c.id;}
207      bool operator> (Row c) const  {return id> c.id;}
208      bool operator==(Row c) const  {return id==c.id;}
209      bool operator!=(Row c) const  {return id!=c.id;}
210   };
211   
212    ///Linear expression of variables and a constant component
213   
214    ///This data structure strores a linear expression of the variables
215    ///(\ref Col "Col"s) and also has a constant component.
216    ///
217    ///There are several ways to access and modify the contents of this
218    ///container.
219    ///- Its it fully compatible with \c std::map<Col,double>, so for expamle
220    ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can
221    ///read and modify the coefficients like
222    ///these.
223    ///\code
224    ///e[v]=5;
225    ///e[v]+=12;
226    ///e.erase(v);
227    ///\endcode
228    ///or you can also iterate through its elements.
229    ///\code
230    ///double s=0;
231    ///for(LpSolverBase::Expr::iterator i=e.begin();i!=e.end();++i)
232    ///  s+=i->second;
233    ///\endcode
234    ///(This code computes the sum of all coefficients).
235    ///- Numbers (<tt>double</tt>'s)
236    ///and variables (\ref Col "Col"s) directly convert to an
237    ///\ref Expr and the usual linear operations are defined, so 
238    ///\code
239    ///v+w
240    ///2*v-3.12*(v-w/2)+2
241    ///v*2.1+(3*v+(v*12+w+6)*3)/2
242    ///\endcode
243    ///are valid \ref Expr "Expr"essions.
244    ///The usual assignment operations are also defined.
245    ///\code
246    ///e=v+w;
247    ///e+=2*v-3.12*(v-w/2)+2;
248    ///e*=3.4;
249    ///e/=5;
250    ///\endcode
251    ///- The constant member can be set and read by \ref constComp()
252    ///\code
253    ///e.constComp()=12;
254    ///double c=e.constComp();
255    ///\endcode
256    ///
257    ///\note \ref clear() not only sets all coefficients to 0 but also
258    ///clears the constant components.
259    ///
260    ///\sa Constr
261    ///
262    class Expr : public std::map<Col,Value>
263    {
264    public:
265      typedef LpSolverBase::Col Key;
266      typedef LpSolverBase::Value Value;
267     
268    protected:
269      typedef std::map<Col,Value> Base;
270     
271      Value const_comp;
272  public:
273      typedef True IsLinExpression;
274      ///\e
275      Expr() : Base(), const_comp(0) { }
276      ///\e
277      Expr(const Key &v) : const_comp(0) {
278        Base::insert(std::make_pair(v, 1));
279      }
280      ///\e
281      Expr(const Value &v) : const_comp(v) {}
282      ///\e
283      void set(const Key &v,const Value &c) {
284        Base::insert(std::make_pair(v, c));
285      }
286      ///\e
287      Value &constComp() { return const_comp; }
288      ///\e
289      const Value &constComp() const { return const_comp; }
290     
291      ///Removes the components with zero coefficient.
292      void simplify() {
293        for (Base::iterator i=Base::begin(); i!=Base::end();) {
294          Base::iterator j=i;
295          ++j;
296          if ((*i).second==0) Base::erase(i);
297          j=i;
298        }
299      }
300
301      ///Removes the coefficients closer to zero than \c tolerance.
302      void simplify(double &tolerance) {
303        for (Base::iterator i=Base::begin(); i!=Base::end();) {
304          Base::iterator j=i;
305          ++j;
306          if (std::fabs((*i).second)<tolerance) Base::erase(i);
307          j=i;
308        }
309      }
310
311      ///Sets all coefficients and the constant component to 0.
312      void clear() {
313        Base::clear();
314        const_comp=0;
315      }
316
317      ///\e
318      Expr &operator+=(const Expr &e) {
319        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
320          (*this)[j->first]+=j->second;
321        const_comp+=e.const_comp;
322        return *this;
323      }
324      ///\e
325      Expr &operator-=(const Expr &e) {
326        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
327          (*this)[j->first]-=j->second;
328        const_comp-=e.const_comp;
329        return *this;
330      }
331      ///\e
332      Expr &operator*=(const Value &c) {
333        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
334          j->second*=c;
335        const_comp*=c;
336        return *this;
337      }
338      ///\e
339      Expr &operator/=(const Value &c) {
340        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
341          j->second/=c;
342        const_comp/=c;
343        return *this;
344      }
345    };
346   
347    ///Linear constraint
348
349    ///This data stucture represents a linear constraint in the LP.
350    ///Basically it is a linear expression with a lower or an upper bound
351    ///(or both). These parts of the constraint can be obtained by the member
352    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
353    ///respectively.
354    ///There are two ways to construct a constraint.
355    ///- You can set the linear expression and the bounds directly
356    ///  by the functions above.
357    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
358    ///  are defined between expressions, or even between constraints whenever
359    ///  it makes sense. Therefore if \c e and \c f are linear expressions and
360    ///  \c s and \c t are numbers, then the followings are valid expressions
361    ///  and thus they can be used directly e.g. in \ref addRow() whenever
362    ///  it makes sense.
363    ///\code
364    ///  e<=s
365    ///  e<=f
366    ///  e==f
367    ///  s<=e<=t
368    ///  e>=t
369    ///\endcode
370    ///\warning The validity of a constraint is checked only at run time, so
371    ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw a
372    ///\ref LogicError exception.
373    class Constr
374    {
375    public:
376      typedef LpSolverBase::Expr Expr;
377      typedef Expr::Key Key;
378      typedef Expr::Value Value;
379     
380//       static const Value INF;
381//       static const Value NaN;
382
383    protected:
384      Expr _expr;
385      Value _lb,_ub;
386    public:
387      ///\e
388      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
389      ///\e
390      Constr(Value lb,const Expr &e,Value ub) :
391        _expr(e), _lb(lb), _ub(ub) {}
392      ///\e
393      Constr(const Expr &e,Value ub) :
394        _expr(e), _lb(NaN), _ub(ub) {}
395      ///\e
396      Constr(Value lb,const Expr &e) :
397        _expr(e), _lb(lb), _ub(NaN) {}
398      ///\e
399      Constr(const Expr &e) :
400        _expr(e), _lb(NaN), _ub(NaN) {}
401      ///\e
402      void clear()
403      {
404        _expr.clear();
405        _lb=_ub=NaN;
406      }
407
408      ///Reference to the linear expression
409      Expr &expr() { return _expr; }
410      ///Cont reference to the linear expression
411      const Expr &expr() const { return _expr; }
412      ///Reference to the lower bound.
413
414      ///\return
415      ///- \ref INF "INF": the constraint is lower unbounded.
416      ///- \ref NaN "NaN": lower bound has not been set.
417      ///- finite number: the lower bound
418      Value &lowerBound() { return _lb; }
419      ///The const version of \ref lowerBound()
420      const Value &lowerBound() const { return _lb; }
421      ///Reference to the upper bound.
422
423      ///\return
424      ///- \ref INF "INF": the constraint is upper unbounded.
425      ///- \ref NaN "NaN": upper bound has not been set.
426      ///- finite number: the upper bound
427      Value &upperBound() { return _ub; }
428      ///The const version of \ref upperBound()
429      const Value &upperBound() const { return _ub; }
430      ///Is the constraint lower bounded?
431      bool lowerBounded() const {
432        using namespace std;
433        return finite(_lb);
434      }
435      ///Is the constraint upper bounded?
436      bool upperBounded() const {
437        using namespace std;
438        return finite(_ub);
439      }
440    };
441   
442    ///Linear expression of rows
443   
444    ///This data structure represents a column of the matrix,
445    ///thas is it strores a linear expression of the dual variables
446    ///(\ref Row "Row"s).
447    ///
448    ///There are several ways to access and modify the contents of this
449    ///container.
450    ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
451    ///if \c e is an DualExpr and \c v
452    ///and \c w are of type \ref Row, then you can
453    ///read and modify the coefficients like
454    ///these.
455    ///\code
456    ///e[v]=5;
457    ///e[v]+=12;
458    ///e.erase(v);
459    ///\endcode
460    ///or you can also iterate through its elements.
461    ///\code
462    ///double s=0;
463    ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
464    ///  s+=i->second;
465    ///\endcode
466    ///(This code computes the sum of all coefficients).
467    ///- Numbers (<tt>double</tt>'s)
468    ///and variables (\ref Row "Row"s) directly convert to an
469    ///\ref DualExpr and the usual linear operations are defined, so
470    ///\code
471    ///v+w
472    ///2*v-3.12*(v-w/2)
473    ///v*2.1+(3*v+(v*12+w)*3)/2
474    ///\endcode
475    ///are valid \ref DualExpr "DualExpr"essions.
476    ///The usual assignment operations are also defined.
477    ///\code
478    ///e=v+w;
479    ///e+=2*v-3.12*(v-w/2);
480    ///e*=3.4;
481    ///e/=5;
482    ///\endcode
483    ///
484    ///\sa Expr
485    ///
486    class DualExpr : public std::map<Row,Value>
487    {
488    public:
489      typedef LpSolverBase::Row Key;
490      typedef LpSolverBase::Value Value;
491     
492    protected:
493      typedef std::map<Row,Value> Base;
494     
495    public:
496      typedef True IsLinExpression;
497      ///\e
498      DualExpr() : Base() { }
499      ///\e
500      DualExpr(const Key &v) {
501        Base::insert(std::make_pair(v, 1));
502      }
503      ///\e
504      void set(const Key &v,const Value &c) {
505        Base::insert(std::make_pair(v, c));
506      }
507     
508      ///Removes the components with zero coefficient.
509      void simplify() {
510        for (Base::iterator i=Base::begin(); i!=Base::end();) {
511          Base::iterator j=i;
512          ++j;
513          if ((*i).second==0) Base::erase(i);
514          j=i;
515        }
516      }
517
518      ///Removes the coefficients closer to zero than \c tolerance.
519      void simplify(double &tolerance) {
520        for (Base::iterator i=Base::begin(); i!=Base::end();) {
521          Base::iterator j=i;
522          ++j;
523          if (std::fabs((*i).second)<tolerance) Base::erase(i);
524          j=i;
525        }
526      }
527
528
529      ///Sets all coefficients to 0.
530      void clear() {
531        Base::clear();
532      }
533
534      ///\e
535      DualExpr &operator+=(const DualExpr &e) {
536        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
537          (*this)[j->first]+=j->second;
538        return *this;
539      }
540      ///\e
541      DualExpr &operator-=(const DualExpr &e) {
542        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
543          (*this)[j->first]-=j->second;
544        return *this;
545      }
546      ///\e
547      DualExpr &operator*=(const Value &c) {
548        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
549          j->second*=c;
550        return *this;
551      }
552      ///\e
553      DualExpr &operator/=(const Value &c) {
554        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
555          j->second/=c;
556        return *this;
557      }
558    };
559   
560
561  protected:
562    _FixId rows;
563    _FixId cols;
564
565    //Abstract virtual functions
566    virtual LpSolverBase &_newLp() = 0;
567    virtual LpSolverBase &_copyLp(){
568      ///\todo This should be implemented here, too,  when we have problem retrieving routines. It can be overriden.
569
570      //Starting:
571      LpSolverBase & newlp(_newLp());
572      return newlp;
573      //return *(LpSolverBase*)0;
574    };
575
576    virtual int _addCol() = 0;
577    virtual int _addRow() = 0;
578    virtual void _eraseCol(int col) = 0;
579    virtual void _eraseRow(int row) = 0;
580    virtual void _getColName(int col,       std::string & name) = 0;
581    virtual void _setColName(int col, const std::string & name) = 0;
582    virtual void _setRowCoeffs(int i,
583                               int length,
584                               int  const * indices,
585                               Value  const * values ) = 0;
586    virtual void _setColCoeffs(int i,
587                               int length,
588                               int  const * indices,
589                               Value  const * values ) = 0;
590    virtual void _setCoeff(int row, int col, Value value) = 0;
591    virtual void _setColLowerBound(int i, Value value) = 0;
592    virtual void _setColUpperBound(int i, Value value) = 0;
593//     virtual void _setRowLowerBound(int i, Value value) = 0;
594//     virtual void _setRowUpperBound(int i, Value value) = 0;
595    virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
596    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
597    virtual void _clearObj()=0;
598//     virtual void _setObj(int length,
599//                          int  const * indices,
600//                          Value  const * values ) = 0;
601    virtual SolveExitStatus _solve() = 0;
602    virtual Value _getPrimal(int i) = 0;
603    virtual Value _getDual(int i) = 0;
604    virtual Value _getPrimalValue() = 0;
605    virtual bool _isBasicCol(int i) = 0;
606    virtual SolutionStatus _getPrimalStatus() = 0;
607    virtual SolutionStatus _getDualStatus() = 0;
608    ///\todo This could be implemented here, too, using _getPrimalStatus() and
609    ///_getDualStatus()
610    virtual ProblemTypes _getProblemType() = 0;
611
612    virtual void _setMax() = 0;
613    virtual void _setMin() = 0;
614   
615    //Own protected stuff
616   
617    //Constant component of the objective function
618    Value obj_const_comp;
619   
620
621
622   
623  public:
624
625    ///\e
626    LpSolverBase() : obj_const_comp(0) {}
627
628    ///\e
629    virtual ~LpSolverBase() {}
630
631    ///Creates a new LP problem
632    LpSolverBase &newLp() {return _newLp();}
633    ///Makes a copy of the LP problem
634    LpSolverBase &copyLp() {return _copyLp();}
635   
636    ///\name Build up and modify the LP
637
638    ///@{
639
640    ///Add a new empty column (i.e a new variable) to the LP
641    Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
642
643    ///\brief Adds several new columns
644    ///(i.e a variables) at once
645    ///
646    ///This magic function takes a container as its argument
647    ///and fills its elements
648    ///with new columns (i.e. variables)
649    ///\param t can be
650    ///- a standard STL compatible iterable container with
651    ///\ref Col as its \c values_type
652    ///like
653    ///\code
654    ///std::vector<LpSolverBase::Col>
655    ///std::list<LpSolverBase::Col>
656    ///\endcode
657    ///- a standard STL compatible iterable container with
658    ///\ref Col as its \c mapped_type
659    ///like
660    ///\code
661    ///std::map<AnyType,LpSolverBase::Col>
662    ///\endcode
663    ///- an iterable lemon \ref concept::WriteMap "write map" like
664    ///\code
665    ///ListGraph::NodeMap<LpSolverBase::Col>
666    ///ListGraph::EdgeMap<LpSolverBase::Col>
667    ///\endcode
668    ///\return The number of the created column.
669#ifdef DOXYGEN
670    template<class T>
671    int addColSet(T &t) { return 0;}
672#else
673    template<class T>
674    typename enable_if<typename T::value_type::LpSolverCol,int>::type
675    addColSet(T &t,dummy<0> = 0) {
676      int s=0;
677      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
678      return s;
679    }
680    template<class T>
681    typename enable_if<typename T::value_type::second_type::LpSolverCol,
682                       int>::type
683    addColSet(T &t,dummy<1> = 1) {
684      int s=0;
685      for(typename T::iterator i=t.begin();i!=t.end();++i) {
686        i->second=addCol();
687        s++;
688      }
689      return s;
690    }
691    template<class T>
692    typename enable_if<typename T::MapIt::Value::LpSolverCol,
693                       int>::type
694    addColSet(T &t,dummy<2> = 2) {
695      int s=0;
696      for(typename T::MapIt i(t); i!=INVALID; ++i)
697        {
698          i.set(addCol());
699          s++;
700        }
701      return s;
702    }
703#endif
704
705    ///Set a column (i.e a dual constraint) of the LP
706
707    ///\param c is the column to be modified
708    ///\param e is a dual linear expression (see \ref DualExpr)
709    ///a better one.
710    void col(Col c,const DualExpr &e) {
711      std::vector<int> indices;
712      std::vector<Value> values;
713      indices.push_back(0);
714      values.push_back(0);
715      for(DualExpr::const_iterator i=e.begin(); i!=e.end(); ++i)
716        if((*i).second!=0) {
717          indices.push_back(rows.floatingId((*i).first.id));
718          values.push_back((*i).second);
719        }
720      _setColCoeffs(cols.floatingId(c.id),indices.size()-1,
721                    &indices[0],&values[0]);
722    }
723
724    ///Add a new column to the LP
725
726    ///\param e is a dual linear expression (see \ref DualExpr)
727    ///\param obj is the corresponding component of the objective
728    ///function. It is 0 by default.
729    ///\return The created column.
730    Col addCol(const DualExpr &e, Value obj=0) {
731      Col c=addCol();
732      col(c,e);
733      objCoeff(c,obj);
734      return c;
735    }
736
737    ///Add a new empty row (i.e a new constraint) to the LP
738
739    ///This function adds a new empty row (i.e a new constraint) to the LP.
740    ///\return The created row
741    Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
742
743    ///\brief Add several new rows
744    ///(i.e a constraints) at once
745    ///
746    ///This magic function takes a container as its argument
747    ///and fills its elements
748    ///with new row (i.e. variables)
749    ///\param t can be
750    ///- a standard STL compatible iterable container with
751    ///\ref Row as its \c values_type
752    ///like
753    ///\code
754    ///std::vector<LpSolverBase::Row>
755    ///std::list<LpSolverBase::Row>
756    ///\endcode
757    ///- a standard STL compatible iterable container with
758    ///\ref Row as its \c mapped_type
759    ///like
760    ///\code
761    ///std::map<AnyType,LpSolverBase::Row>
762    ///\endcode
763    ///- an iterable lemon \ref concept::WriteMap "write map" like
764    ///\code
765    ///ListGraph::NodeMap<LpSolverBase::Row>
766    ///ListGraph::EdgeMap<LpSolverBase::Row>
767    ///\endcode
768    ///\return The number of rows created.
769#ifdef DOXYGEN
770    template<class T>
771    int addRowSet(T &t) { return 0;}
772#else
773    template<class T>
774    typename enable_if<typename T::value_type::LpSolverRow,int>::type
775    addRowSet(T &t,dummy<0> = 0) {
776      int s=0;
777      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
778      return s;
779    }
780    template<class T>
781    typename enable_if<typename T::value_type::second_type::LpSolverRow,
782                       int>::type
783    addRowSet(T &t,dummy<1> = 1) {
784      int s=0;
785      for(typename T::iterator i=t.begin();i!=t.end();++i) {
786        i->second=addRow();
787        s++;
788      }
789      return s;
790    }
791    template<class T>
792    typename enable_if<typename T::MapIt::Value::LpSolverRow,
793                       int>::type
794    addRowSet(T &t,dummy<2> = 2) {
795      int s=0;
796      for(typename T::MapIt i(t); i!=INVALID; ++i)
797        {
798          i.set(addRow());
799          s++;
800        }
801      return s;
802    }
803#endif
804
805    ///Set a row (i.e a constraint) of the LP
806
807    ///\param r is the row to be modified
808    ///\param l is lower bound (-\ref INF means no bound)
809    ///\param e is a linear expression (see \ref Expr)
810    ///\param u is the upper bound (\ref INF means no bound)
811    ///\bug This is a temportary function. The interface will change to
812    ///a better one.
813    ///\todo Option to control whether a constraint with a single variable is
814    ///added or not.
815    void row(Row r, Value l,const Expr &e, Value u) {
816      std::vector<int> indices;
817      std::vector<Value> values;
818      indices.push_back(0);
819      values.push_back(0);
820      for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
821        if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
822          indices.push_back(cols.floatingId((*i).first.id));
823          values.push_back((*i).second);
824        }
825      _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
826                    &indices[0],&values[0]);
827//       _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
828//       _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
829       _setRowBounds(rows.floatingId(r.id),l-e.constComp(),u-e.constComp());
830    }
831
832    ///Set a row (i.e a constraint) of the LP
833
834    ///\param r is the row to be modified
835    ///\param c is a linear expression (see \ref Constr)
836    void row(Row r, const Constr &c) {
837      row(r,
838             c.lowerBounded()?c.lowerBound():-INF,
839             c.expr(),
840             c.upperBounded()?c.upperBound():INF);
841    }
842
843    ///Add a new row (i.e a new constraint) to the LP
844
845    ///\param l is the lower bound (-\ref INF means no bound)
846    ///\param e is a linear expression (see \ref Expr)
847    ///\param u is the upper bound (\ref INF means no bound)
848    ///\return The created row.
849    ///\bug This is a temportary function. The interface will change to
850    ///a better one.
851    Row addRow(Value l,const Expr &e, Value u) {
852      Row r=addRow();
853      row(r,l,e,u);
854      return r;
855    }
856
857    ///Add a new row (i.e a new constraint) to the LP
858
859    ///\param c is a linear expression (see \ref Constr)
860    ///\return The created row.
861    Row addRow(const Constr &c) {
862      Row r=addRow();
863      row(r,c);
864      return r;
865    }
866    ///Erase a coloumn (i.e a variable) from the LP
867
868    ///\param c is the coloumn to be deleted
869    ///\todo Please check this
870    void eraseCol(Col c) {
871      _eraseCol(cols.floatingId(c.id));
872      cols.erase(c.id);
873    }
874    ///Erase a  row (i.e a constraint) from the LP
875
876    ///\param r is the row to be deleted
877    ///\todo Please check this
878    void eraseRow(Row r) {
879      _eraseRow(rows.floatingId(r.id));
880      rows.erase(r.id);
881    }
882
883    /// Get the name of a column
884   
885    ///\param c is the coresponding coloumn
886    ///\return The name of the colunm
887    std::string ColName(Col c){
888      std::string name;
889      _getColName(cols.floatingId(c.id), name);
890      return name;
891    }
892   
893    /// Set the name of a column
894   
895    ///\param c is the coresponding coloumn
896    ///\param name The name to be given
897    void ColName(Col c, const std::string & name){
898      _setColName(cols.floatingId(c.id), name);
899    }
900   
901    /// Set an element of the coefficient matrix of the LP
902
903    ///\param r is the row of the element to be modified
904    ///\param c is the coloumn of the element to be modified
905    ///\param val is the new value of the coefficient
906
907    void Coeff(Row r, Col c, Value val){
908      _setCoeff(rows.floatingId(r.id),cols.floatingId(c.id), val);
909    }
910
911    /// Set the lower bound of a column (i.e a variable)
912
913    /// The lower bound of a variable (column) has to be given by an
914    /// extended number of type Value, i.e. a finite number of type
915    /// Value or -\ref INF.
916    void colLowerBound(Col c, Value value) {
917      _setColLowerBound(cols.floatingId(c.id),value);
918    }
919   
920    ///\brief Set the lower bound of  several columns
921    ///(i.e a variables) at once
922    ///
923    ///This magic function takes a container as its argument
924    ///and applies the function on all of its elements.
925    /// The lower bound of a variable (column) has to be given by an
926    /// extended number of type Value, i.e. a finite number of type
927    /// Value or -\ref INF.
928#ifdef DOXYGEN
929    template<class T>
930    void colLowerBound(T &t, Value value) { return 0;}
931#else
932    template<class T>
933    typename enable_if<typename T::value_type::LpSolverCol,void>::type
934    colLowerBound(T &t, Value value,dummy<0> = 0) {
935      for(typename T::iterator i=t.begin();i!=t.end();++i) {
936        colLowerBound(*i, value);
937      }
938    }
939    template<class T>
940    typename enable_if<typename T::value_type::second_type::LpSolverCol,
941                       void>::type
942    colLowerBound(T &t, Value value,dummy<1> = 1) {
943      for(typename T::iterator i=t.begin();i!=t.end();++i) {
944        colLowerBound(i->second, value);
945      }
946    }
947    template<class T>
948    typename enable_if<typename T::MapIt::Value::LpSolverCol,
949                       void>::type
950    colLowerBound(T &t, Value value,dummy<2> = 2) {
951      for(typename T::MapIt i(t); i!=INVALID; ++i){
952        colLowerBound(*i, value);
953      }
954    }
955#endif
956   
957    /// Set the upper bound of a column (i.e a variable)
958
959    /// The upper bound of a variable (column) has to be given by an
960    /// extended number of type Value, i.e. a finite number of type
961    /// Value or \ref INF.
962    void colUpperBound(Col c, Value value) {
963      _setColUpperBound(cols.floatingId(c.id),value);
964    };
965
966    ///\brief Set the lower bound of  several columns
967    ///(i.e a variables) at once
968    ///
969    ///This magic function takes a container as its argument
970    ///and applies the function on all of its elements.
971    /// The upper bound of a variable (column) has to be given by an
972    /// extended number of type Value, i.e. a finite number of type
973    /// Value or \ref INF.
974#ifdef DOXYGEN
975    template<class T>
976    void colUpperBound(T &t, Value value) { return 0;}
977#else
978    template<class T>
979    typename enable_if<typename T::value_type::LpSolverCol,void>::type
980    colUpperBound(T &t, Value value,dummy<0> = 0) {
981      for(typename T::iterator i=t.begin();i!=t.end();++i) {
982        colUpperBound(*i, value);
983      }
984    }
985    template<class T>
986    typename enable_if<typename T::value_type::second_type::LpSolverCol,
987                       void>::type
988    colUpperBound(T &t, Value value,dummy<1> = 1) {
989      for(typename T::iterator i=t.begin();i!=t.end();++i) {
990        colUpperBound(i->second, value);
991      }
992    }
993    template<class T>
994    typename enable_if<typename T::MapIt::Value::LpSolverCol,
995                       void>::type
996    colUpperBound(T &t, Value value,dummy<2> = 2) {
997      for(typename T::MapIt i(t); i!=INVALID; ++i){
998        colUpperBound(*i, value);
999      }
1000    }
1001#endif
1002
1003    /// Set the lower and the upper bounds of a column (i.e a variable)
1004
1005    /// The lower and the upper bounds of
1006    /// a variable (column) have to be given by an
1007    /// extended number of type Value, i.e. a finite number of type
1008    /// Value, -\ref INF or \ref INF.
1009    void colBounds(Col c, Value lower, Value upper) {
1010      _setColLowerBound(cols.floatingId(c.id),lower);
1011      _setColUpperBound(cols.floatingId(c.id),upper);
1012    }
1013   
1014    ///\brief Set the lower and the upper bound of several columns
1015    ///(i.e a variables) at once
1016    ///
1017    ///This magic function takes a container as its argument
1018    ///and applies the function on all of its elements.
1019    /// The lower and the upper bounds of
1020    /// a variable (column) have to be given by an
1021    /// extended number of type Value, i.e. a finite number of type
1022    /// Value, -\ref INF or \ref INF.
1023#ifdef DOXYGEN
1024    template<class T>
1025    void colBounds(T &t, Value lower, Value upper) { return 0;}
1026#else
1027    template<class T>
1028    typename enable_if<typename T::value_type::LpSolverCol,void>::type
1029    colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
1030      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1031        colBounds(*i, lower, upper);
1032      }
1033    }
1034    template<class T>
1035    typename enable_if<typename T::value_type::second_type::LpSolverCol,
1036                       void>::type
1037    colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
1038      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1039        colBounds(i->second, lower, upper);
1040      }
1041    }
1042    template<class T>
1043    typename enable_if<typename T::MapIt::Value::LpSolverCol,
1044                       void>::type
1045    colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
1046      for(typename T::MapIt i(t); i!=INVALID; ++i){
1047        colBounds(*i, lower, upper);
1048      }
1049    }
1050#endif
1051   
1052//     /// Set the lower bound of a row (i.e a constraint)
1053
1054//     /// The lower bound of a linear expression (row) has to be given by an
1055//     /// extended number of type Value, i.e. a finite number of type
1056//     /// Value or -\ref INF.
1057//     void rowLowerBound(Row r, Value value) {
1058//       _setRowLowerBound(rows.floatingId(r.id),value);
1059//     };
1060//     /// Set the upper bound of a row (i.e a constraint)
1061
1062//     /// The upper bound of a linear expression (row) has to be given by an
1063//     /// extended number of type Value, i.e. a finite number of type
1064//     /// Value or \ref INF.
1065//     void rowUpperBound(Row r, Value value) {
1066//       _setRowUpperBound(rows.floatingId(r.id),value);
1067//     };
1068
1069    /// Set the lower and the upper bounds of a row (i.e a constraint)
1070
1071    /// The lower and the upper bounds of
1072    /// a constraint (row) have to be given by an
1073    /// extended number of type Value, i.e. a finite number of type
1074    /// Value, -\ref INF or \ref INF.
1075    void rowBounds(Row c, Value lower, Value upper) {
1076      _setRowBounds(rows.floatingId(c.id),lower, upper);
1077      // _setRowUpperBound(rows.floatingId(c.id),upper);
1078    }
1079   
1080    ///Set an element of the objective function
1081    void objCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
1082    ///Set the objective function
1083   
1084    ///\param e is a linear expression of type \ref Expr.
1085    ///\bug Is should be called obj()
1086    void setObj(Expr e) {
1087      _clearObj();
1088      for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
1089        objCoeff((*i).first,(*i).second);
1090      obj_const_comp=e.constComp();
1091    }
1092
1093    ///Maximize
1094    void max() { _setMax(); }
1095    ///Minimize
1096    void min() { _setMin(); }
1097
1098   
1099    ///@}
1100
1101
1102    ///\name Solve the LP
1103
1104    ///@{
1105
1106    ///\e Solve the LP problem at hand
1107    ///
1108    ///\return The result of the optimization procedure. Possible values and their meanings can be found in the documentation of \ref SolveExitStatus.
1109    ///
1110    ///\todo Which method is used to solve the problem
1111    SolveExitStatus solve() { return _solve(); }
1112   
1113    ///@}
1114   
1115    ///\name Obtain the solution
1116
1117    ///@{
1118
1119    /// The status of the primal problem (the original LP problem)
1120    SolutionStatus primalStatus() {
1121      return _getPrimalStatus();
1122    }
1123
1124    /// The status of the dual (of the original LP) problem
1125    SolutionStatus dualStatus() {
1126      return _getDualStatus();
1127    }
1128
1129    ///The type of the original LP problem
1130    ProblemTypes problemType() {
1131      return _getProblemType();
1132    }
1133
1134    ///\e
1135    Value primal(Col c) { return _getPrimal(cols.floatingId(c.id)); }
1136
1137    ///\e
1138    Value dual(Row r) { return _getDual(rows.floatingId(r.id)); }
1139
1140    ///\e
1141    bool isBasicCol(Col c) { return _isBasicCol(cols.floatingId(c.id)); }
1142
1143    ///\e
1144
1145    ///\return
1146    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1147    /// of the primal problem, depending on whether we minimize or maximize.
1148    ///- \ref NaN if no primal solution is found.
1149    ///- The (finite) objective value if an optimal solution is found.
1150    Value primalValue() { return _getPrimalValue()+obj_const_comp;}
1151    ///@}
1152   
1153  }; 
1154
1155  ///\e
1156 
1157  ///\relates LpSolverBase::Expr
1158  ///
1159  inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
1160                                      const LpSolverBase::Expr &b)
1161  {
1162    LpSolverBase::Expr tmp(a);
1163    tmp+=b;
1164    return tmp;
1165  }
1166  ///\e
1167 
1168  ///\relates LpSolverBase::Expr
1169  ///
1170  inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
1171                                      const LpSolverBase::Expr &b)
1172  {
1173    LpSolverBase::Expr tmp(a);
1174    tmp-=b;
1175    return tmp;
1176  }
1177  ///\e
1178 
1179  ///\relates LpSolverBase::Expr
1180  ///
1181  inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
1182                                      const LpSolverBase::Value &b)
1183  {
1184    LpSolverBase::Expr tmp(a);
1185    tmp*=b;
1186    return tmp;
1187  }
1188 
1189  ///\e
1190 
1191  ///\relates LpSolverBase::Expr
1192  ///
1193  inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
1194                                      const LpSolverBase::Expr &b)
1195  {
1196    LpSolverBase::Expr tmp(b);
1197    tmp*=a;
1198    return tmp;
1199  }
1200  ///\e
1201 
1202  ///\relates LpSolverBase::Expr
1203  ///
1204  inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
1205                                      const LpSolverBase::Value &b)
1206  {
1207    LpSolverBase::Expr tmp(a);
1208    tmp/=b;
1209    return tmp;
1210  }
1211 
1212  ///\e
1213 
1214  ///\relates LpSolverBase::Constr
1215  ///
1216  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1217                                         const LpSolverBase::Expr &f)
1218  {
1219    return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
1220  }
1221
1222  ///\e
1223 
1224  ///\relates LpSolverBase::Constr
1225  ///
1226  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
1227                                         const LpSolverBase::Expr &f)
1228  {
1229    return LpSolverBase::Constr(e,f);
1230  }
1231
1232  ///\e
1233 
1234  ///\relates LpSolverBase::Constr
1235  ///
1236  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1237                                         const LpSolverBase::Value &f)
1238  {
1239    return LpSolverBase::Constr(e,f);
1240  }
1241
1242  ///\e
1243 
1244  ///\relates LpSolverBase::Constr
1245  ///
1246  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1247                                         const LpSolverBase::Expr &f)
1248  {
1249    return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
1250  }
1251
1252
1253  ///\e
1254 
1255  ///\relates LpSolverBase::Constr
1256  ///
1257  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
1258                                         const LpSolverBase::Expr &f)
1259  {
1260    return LpSolverBase::Constr(f,e);
1261  }
1262
1263
1264  ///\e
1265 
1266  ///\relates LpSolverBase::Constr
1267  ///
1268  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1269                                         const LpSolverBase::Value &f)
1270  {
1271    return LpSolverBase::Constr(f,e);
1272  }
1273
1274  ///\e
1275 
1276  ///\relates LpSolverBase::Constr
1277  ///
1278  inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1279                                         const LpSolverBase::Expr &f)
1280  {
1281    return LpSolverBase::Constr(0,e-f,0);
1282  }
1283
1284  ///\e
1285 
1286  ///\relates LpSolverBase::Constr
1287  ///
1288  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
1289                                         const LpSolverBase::Constr&c)
1290  {
1291    LpSolverBase::Constr tmp(c);
1292    ///\todo Create an own exception type.
1293    if(!isnan(tmp.lowerBound())) throw LogicError();
1294    else tmp.lowerBound()=n;
1295    return tmp;
1296  }
1297  ///\e
1298 
1299  ///\relates LpSolverBase::Constr
1300  ///
1301  inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
1302                                         const LpSolverBase::Value &n)
1303  {
1304    LpSolverBase::Constr tmp(c);
1305    ///\todo Create an own exception type.
1306    if(!isnan(tmp.upperBound())) throw LogicError();
1307    else tmp.upperBound()=n;
1308    return tmp;
1309  }
1310
1311  ///\e
1312 
1313  ///\relates LpSolverBase::Constr
1314  ///
1315  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
1316                                         const LpSolverBase::Constr&c)
1317  {
1318    LpSolverBase::Constr tmp(c);
1319    ///\todo Create an own exception type.
1320    if(!isnan(tmp.upperBound())) throw LogicError();
1321    else tmp.upperBound()=n;
1322    return tmp;
1323  }
1324  ///\e
1325 
1326  ///\relates LpSolverBase::Constr
1327  ///
1328  inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
1329                                         const LpSolverBase::Value &n)
1330  {
1331    LpSolverBase::Constr tmp(c);
1332    ///\todo Create an own exception type.
1333    if(!isnan(tmp.lowerBound())) throw LogicError();
1334    else tmp.lowerBound()=n;
1335    return tmp;
1336  }
1337
1338  ///\e
1339 
1340  ///\relates LpSolverBase::DualExpr
1341  ///
1342  inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
1343                                      const LpSolverBase::DualExpr &b)
1344  {
1345    LpSolverBase::DualExpr tmp(a);
1346    tmp+=b;
1347    return tmp;
1348  }
1349  ///\e
1350 
1351  ///\relates LpSolverBase::DualExpr
1352  ///
1353  inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
1354                                      const LpSolverBase::DualExpr &b)
1355  {
1356    LpSolverBase::DualExpr tmp(a);
1357    tmp-=b;
1358    return tmp;
1359  }
1360  ///\e
1361 
1362  ///\relates LpSolverBase::DualExpr
1363  ///
1364  inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
1365                                      const LpSolverBase::Value &b)
1366  {
1367    LpSolverBase::DualExpr tmp(a);
1368    tmp*=b;
1369    return tmp;
1370  }
1371 
1372  ///\e
1373 
1374  ///\relates LpSolverBase::DualExpr
1375  ///
1376  inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
1377                                      const LpSolverBase::DualExpr &b)
1378  {
1379    LpSolverBase::DualExpr tmp(b);
1380    tmp*=a;
1381    return tmp;
1382  }
1383  ///\e
1384 
1385  ///\relates LpSolverBase::DualExpr
1386  ///
1387  inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
1388                                      const LpSolverBase::Value &b)
1389  {
1390    LpSolverBase::DualExpr tmp(a);
1391    tmp/=b;
1392    return tmp;
1393  }
1394 
1395
1396} //namespace lemon
1397
1398#endif //LEMON_LP_BASE_H
Note: See TracBrowser for help on using the repository browser.