COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/lp_base.h @ 1771:5faaa9880d4d

Last change on this file since 1771:5faaa9880d4d was 1771:5faaa9880d4d, checked in by Alpar Juttner, 14 years ago

(Dual)Expr::simplify(double tolerance) added

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