COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/lp_base.h @ 1946:17eb3eaad9f8

Last change on this file since 1946:17eb3eaad9f8 was 1908:e225719bde6b, checked in by Alpar Juttner, 18 years ago

Better doc.

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