COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/lp_base.h @ 1766:6c59b1386fe8

Last change on this file since 1766:6c59b1386fe8 was 1766:6c59b1386fe8, checked in by Alpar Juttner, 18 years ago

Tons of todos have been removed.

File size: 34.2 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      ///Sets all coefficients and the constant component to 0.
297      void clear() {
298        Base::clear();
299        const_comp=0;
300      }
301
302      ///\e
303      Expr &operator+=(const Expr &e) {
304        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
305          (*this)[j->first]+=j->second;
306        const_comp+=e.const_comp;
307        return *this;
308      }
309      ///\e
310      Expr &operator-=(const Expr &e) {
311        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
312          (*this)[j->first]-=j->second;
313        const_comp-=e.const_comp;
314        return *this;
315      }
316      ///\e
317      Expr &operator*=(const Value &c) {
318        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
319          j->second*=c;
320        const_comp*=c;
321        return *this;
322      }
323      ///\e
324      Expr &operator/=(const Value &c) {
325        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
326          j->second/=c;
327        const_comp/=c;
328        return *this;
329      }
330    };
331   
332    ///Linear constraint
333
334    ///This data stucture represents a linear constraint in the LP.
335    ///Basically it is a linear expression with a lower or an upper bound
336    ///(or both). These parts of the constraint can be obtained by the member
337    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
338    ///respectively.
339    ///There are two ways to construct a constraint.
340    ///- You can set the linear expression and the bounds directly
341    ///  by the functions above.
342    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
343    ///  are defined between expressions, or even between constraints whenever
344    ///  it makes sense. Therefore if \c e and \c f are linear expressions and
345    ///  \c s and \c t are numbers, then the followings are valid expressions
346    ///  and thus they can be used directly e.g. in \ref addRow() whenever
347    ///  it makes sense.
348    ///  \code
349    ///  e<=s
350    ///  e<=f
351    ///  s<=e<=t
352    ///  e>=t
353    ///  \endcode
354    ///\warning The validity of a constraint is checked only at run time, so
355    ///e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw a
356    ///\ref LogicError exception.
357    class Constr
358    {
359    public:
360      typedef LpSolverBase::Expr Expr;
361      typedef Expr::Key Key;
362      typedef Expr::Value Value;
363     
364//       static const Value INF;
365//       static const Value NaN;
366
367    protected:
368      Expr _expr;
369      Value _lb,_ub;
370    public:
371      ///\e
372      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
373      ///\e
374      Constr(Value lb,const Expr &e,Value ub) :
375        _expr(e), _lb(lb), _ub(ub) {}
376      ///\e
377      Constr(const Expr &e,Value ub) :
378        _expr(e), _lb(NaN), _ub(ub) {}
379      ///\e
380      Constr(Value lb,const Expr &e) :
381        _expr(e), _lb(lb), _ub(NaN) {}
382      ///\e
383      Constr(const Expr &e) :
384        _expr(e), _lb(NaN), _ub(NaN) {}
385      ///\e
386      void clear()
387      {
388        _expr.clear();
389        _lb=_ub=NaN;
390      }
391
392      ///Reference to the linear expression
393      Expr &expr() { return _expr; }
394      ///Cont reference to the linear expression
395      const Expr &expr() const { return _expr; }
396      ///Reference to the lower bound.
397
398      ///\return
399      ///- \ref INF "INF": the constraint is lower unbounded.
400      ///- \ref NaN "NaN": lower bound has not been set.
401      ///- finite number: the lower bound
402      Value &lowerBound() { return _lb; }
403      ///The const version of \ref lowerBound()
404      const Value &lowerBound() const { return _lb; }
405      ///Reference to the upper bound.
406
407      ///\return
408      ///- \ref INF "INF": the constraint is upper unbounded.
409      ///- \ref NaN "NaN": upper bound has not been set.
410      ///- finite number: the upper bound
411      Value &upperBound() { return _ub; }
412      ///The const version of \ref upperBound()
413      const Value &upperBound() const { return _ub; }
414      ///Is the constraint lower bounded?
415      bool lowerBounded() const {
416        using namespace std;
417        return finite(_lb);
418      }
419      ///Is the constraint upper bounded?
420      bool upperBounded() const {
421        using namespace std;
422        return finite(_ub);
423      }
424    };
425   
426    ///Linear expression of rows
427   
428    ///This data structure represents a column of the matrix,
429    ///thas is it strores a linear expression of the dual variables
430    ///(\ref Row "Row"s).
431    ///
432    ///There are several ways to access and modify the contents of this
433    ///container.
434    ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
435    ///if \c e is an DualExpr and \c v
436    ///and \c w are of type \ref Row, then you can
437    ///read and modify the coefficients like
438    ///these.
439    ///\code
440    ///e[v]=5;
441    ///e[v]+=12;
442    ///e.erase(v);
443    ///\endcode
444    ///or you can also iterate through its elements.
445    ///\code
446    ///double s=0;
447    ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
448    ///  s+=i->second;
449    ///\endcode
450    ///(This code computes the sum of all coefficients).
451    ///- Numbers (<tt>double</tt>'s)
452    ///and variables (\ref Row "Row"s) directly convert to an
453    ///\ref DualExpr and the usual linear operations are defined so 
454    ///\code
455    ///v+w
456    ///2*v-3.12*(v-w/2)
457    ///v*2.1+(3*v+(v*12+w)*3)/2
458    ///\endcode
459    ///are valid \ref DualExpr "DualExpr"essions.
460    ///The usual assignment operations are also defined.
461    ///\code
462    ///e=v+w;
463    ///e+=2*v-3.12*(v-w/2);
464    ///e*=3.4;
465    ///e/=5;
466    ///\endcode
467    ///
468    ///\sa Expr
469    ///
470    class DualExpr : public std::map<Row,Value>
471    {
472    public:
473      typedef LpSolverBase::Row Key;
474      typedef LpSolverBase::Value Value;
475     
476    protected:
477      typedef std::map<Row,Value> Base;
478     
479    public:
480      typedef True IsLinExpression;
481      ///\e
482      DualExpr() : Base() { }
483      ///\e
484      DualExpr(const Key &v) {
485        Base::insert(std::make_pair(v, 1));
486      }
487      ///\e
488      void set(const Key &v,const Value &c) {
489        Base::insert(std::make_pair(v, c));
490      }
491     
492      ///Removes the components with zero coefficient.
493      void simplify() {
494        for (Base::iterator i=Base::begin(); i!=Base::end();) {
495          Base::iterator j=i;
496          ++j;
497          if ((*i).second==0) Base::erase(i);
498          j=i;
499        }
500      }
501
502      ///Sets all coefficients to 0.
503      void clear() {
504        Base::clear();
505      }
506
507      ///\e
508      DualExpr &operator+=(const DualExpr &e) {
509        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
510          (*this)[j->first]+=j->second;
511        return *this;
512      }
513      ///\e
514      DualExpr &operator-=(const DualExpr &e) {
515        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
516          (*this)[j->first]-=j->second;
517        return *this;
518      }
519      ///\e
520      DualExpr &operator*=(const Value &c) {
521        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
522          j->second*=c;
523        return *this;
524      }
525      ///\e
526      DualExpr &operator/=(const Value &c) {
527        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
528          j->second/=c;
529        return *this;
530      }
531    };
532   
533
534  protected:
535    _FixId rows;
536    _FixId cols;
537
538    //Abstract virtual functions
539    virtual LpSolverBase &_newLp() = 0;
540    virtual LpSolverBase &_copyLp(){
541      ///\todo This should be implemented here, too,  when we have problem retrieving routines. It can be overriden.
542
543      //Starting:
544      LpSolverBase & newlp(_newLp());
545      return newlp;
546      //return *(LpSolverBase*)0;
547    };
548
549    virtual int _addCol() = 0;
550    virtual int _addRow() = 0;
551    virtual void _eraseCol(int col) = 0;
552    virtual void _eraseRow(int row) = 0;
553    virtual void _setRowCoeffs(int i,
554                               int length,
555                               int  const * indices,
556                               Value  const * values ) = 0;
557    virtual void _setColCoeffs(int i,
558                               int length,
559                               int  const * indices,
560                               Value  const * values ) = 0;
561    virtual void _setCoeff(int row, int col, Value value) = 0;
562    virtual void _setColLowerBound(int i, Value value) = 0;
563    virtual void _setColUpperBound(int i, Value value) = 0;
564//     virtual void _setRowLowerBound(int i, Value value) = 0;
565//     virtual void _setRowUpperBound(int i, Value value) = 0;
566    virtual void _setRowBounds(int i, Value lower, Value upper) = 0;
567    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
568    virtual void _clearObj()=0;
569//     virtual void _setObj(int length,
570//                          int  const * indices,
571//                          Value  const * values ) = 0;
572    virtual SolveExitStatus _solve() = 0;
573    virtual Value _getPrimal(int i) = 0;
574    virtual Value _getPrimalValue() = 0;
575    virtual SolutionStatus _getPrimalStatus() = 0;
576    virtual SolutionStatus _getDualStatus() = 0;
577    ///\todo This could be implemented here, too, using _getPrimalStatus() and
578    ///_getDualStatus()
579    virtual ProblemTypes _getProblemType() = 0;
580
581    virtual void _setMax() = 0;
582    virtual void _setMin() = 0;
583   
584    //Own protected stuff
585   
586    //Constant component of the objective function
587    Value obj_const_comp;
588   
589
590
591   
592  public:
593
594    ///\e
595    LpSolverBase() : obj_const_comp(0) {}
596
597    ///\e
598    virtual ~LpSolverBase() {}
599
600    ///Creates a new LP problem
601    LpSolverBase &newLp() {return _newLp();}
602    ///Makes a copy of the LP problem
603    LpSolverBase &copyLp() {return _copyLp();}
604   
605    ///\name Build up and modify the LP
606
607    ///@{
608
609    ///Add a new empty column (i.e a new variable) to the LP
610    Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
611
612    ///\brief Adds several new columns
613    ///(i.e a variables) at once
614    ///
615    ///This magic function takes a container as its argument
616    ///and fills its elements
617    ///with new columns (i.e. variables)
618    ///\param t can be
619    ///- a standard STL compatible iterable container with
620    ///\ref Col as its \c values_type
621    ///like
622    ///\code
623    ///std::vector<LpSolverBase::Col>
624    ///std::list<LpSolverBase::Col>
625    ///\endcode
626    ///- a standard STL compatible iterable container with
627    ///\ref Col as its \c mapped_type
628    ///like
629    ///\code
630    ///std::map<AnyType,LpSolverBase::Col>
631    ///\endcode
632    ///- an iterable lemon \ref concept::WriteMap "write map" like
633    ///\code
634    ///ListGraph::NodeMap<LpSolverBase::Col>
635    ///ListGraph::EdgeMap<LpSolverBase::Col>
636    ///\endcode
637    ///\return The number of the created column.
638#ifdef DOXYGEN
639    template<class T>
640    int addColSet(T &t) { return 0;}
641#else
642    template<class T>
643    typename enable_if<typename T::value_type::LpSolverCol,int>::type
644    addColSet(T &t,dummy<0> = 0) {
645      int s=0;
646      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
647      return s;
648    }
649    template<class T>
650    typename enable_if<typename T::value_type::second_type::LpSolverCol,
651                       int>::type
652    addColSet(T &t,dummy<1> = 1) {
653      int s=0;
654      for(typename T::iterator i=t.begin();i!=t.end();++i) {
655        i->second=addCol();
656        s++;
657      }
658      return s;
659    }
660    template<class T>
661    typename enable_if<typename T::ValueSet::value_type::LpSolverCol,
662                       int>::type
663    addColSet(T &t,dummy<2> = 2) {
664      ///\bug <tt>return addColSet(t.valueSet());</tt> should also work.
665      int s=0;
666      for(typename T::ValueSet::iterator i=t.valueSet().begin();
667          i!=t.valueSet().end();
668          ++i)
669        {
670          *i=addCol();
671          s++;
672        }
673      return s;
674    }
675#endif
676
677    ///Set a column (i.e a dual constraint) of the LP
678
679    ///\param c is the column to be modified
680    ///\param e is a dual linear expression (see \ref DualExpr)
681    ///\bug This is a temporary function. The interface will change to
682    ///a better one.
683    void setCol(Col c,const DualExpr &e) {
684      std::vector<int> indices;
685      std::vector<Value> values;
686      indices.push_back(0);
687      values.push_back(0);
688      for(DualExpr::const_iterator i=e.begin(); i!=e.end(); ++i)
689        if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
690          indices.push_back(cols.floatingId((*i).first.id));
691          values.push_back((*i).second);
692        }
693      _setColCoeffs(cols.floatingId(c.id),indices.size()-1,
694                    &indices[0],&values[0]);
695    }
696
697    ///Add a new column to the LP
698
699    ///\param e is a dual linear expression (see \ref DualExpr)
700    ///\param obj is the corresponding component of the objective
701    ///function. It is 0 by default.
702    ///\return The created column.
703    ///\bug This is a temportary function. The interface will change to
704    ///a better one.
705    Col addCol(const DualExpr &e, Value obj=0) {
706      Col c=addCol();
707      setCol(c,e);
708      objCoeff(c,obj);
709      return c;
710    }
711
712    ///Add a new empty row (i.e a new constraint) to the LP
713
714    ///This function adds a new empty row (i.e a new constraint) to the LP.
715    ///\return The created row
716    Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
717
718    ///\brief Add several new rows
719    ///(i.e a constraints) at once
720    ///
721    ///This magic function takes a container as its argument
722    ///and fills its elements
723    ///with new row (i.e. variables)
724    ///\param t can be
725    ///- a standard STL compatible iterable container with
726    ///\ref Row as its \c values_type
727    ///like
728    ///\code
729    ///std::vector<LpSolverBase::Row>
730    ///std::list<LpSolverBase::Row>
731    ///\endcode
732    ///- a standard STL compatible iterable container with
733    ///\ref Row as its \c mapped_type
734    ///like
735    ///\code
736    ///std::map<AnyType,LpSolverBase::Row>
737    ///\endcode
738    ///- an iterable lemon \ref concept::WriteMap "write map" like
739    ///\code
740    ///ListGraph::NodeMap<LpSolverBase::Row>
741    ///ListGraph::EdgeMap<LpSolverBase::Row>
742    ///\endcode
743    ///\return The number of rows created.
744#ifdef DOXYGEN
745    template<class T>
746    int addRowSet(T &t) { return 0;}
747#else
748    template<class T>
749    typename enable_if<typename T::value_type::LpSolverRow,int>::type
750    addRowSet(T &t,dummy<0> = 0) {
751      int s=0;
752      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
753      return s;
754    }
755    template<class T>
756    typename enable_if<typename T::value_type::second_type::LpSolverRow,
757                       int>::type
758    addRowSet(T &t,dummy<1> = 1) {
759      int s=0;
760      for(typename T::iterator i=t.begin();i!=t.end();++i) {
761        i->second=addRow();
762        s++;
763      }
764      return s;
765    }
766    template<class T>
767    typename enable_if<typename T::ValueSet::value_type::LpSolverRow,
768                       int>::type
769    addRowSet(T &t,dummy<2> = 2) {
770      ///\bug <tt>return addRowSet(t.valueSet());</tt> should also work.
771      int s=0;
772      for(typename T::ValueSet::iterator i=t.valueSet().begin();
773          i!=t.valueSet().end();
774          ++i)
775        {
776          *i=addRow();
777          s++;
778        }
779      return s;
780    }
781#endif
782
783    ///Set a row (i.e a constraint) of the LP
784
785    ///\param r is the row to be modified
786    ///\param l is lower bound (-\ref INF means no bound)
787    ///\param e is a linear expression (see \ref Expr)
788    ///\param u is the upper bound (\ref INF means no bound)
789    ///\bug This is a temportary function. The interface will change to
790    ///a better one.
791    ///\todo Option to control whether a constraint with a single variable is
792    ///added or not.
793    void setRow(Row r, Value l,const Expr &e, Value u) {
794      std::vector<int> indices;
795      std::vector<Value> values;
796      indices.push_back(0);
797      values.push_back(0);
798      for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
799        if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
800          indices.push_back(cols.floatingId((*i).first.id));
801          values.push_back((*i).second);
802        }
803      _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
804                    &indices[0],&values[0]);
805//       _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
806//       _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
807       _setRowBounds(rows.floatingId(r.id),l-e.constComp(),u-e.constComp());
808    }
809
810    ///Set a row (i.e a constraint) of the LP
811
812    ///\param r is the row to be modified
813    ///\param c is a linear expression (see \ref Constr)
814    void setRow(Row r, const Constr &c) {
815      setRow(r,
816             c.lowerBounded()?c.lowerBound():-INF,
817             c.expr(),
818             c.upperBounded()?c.upperBound():INF);
819    }
820
821    ///Add a new row (i.e a new constraint) to the LP
822
823    ///\param l is the lower bound (-\ref INF means no bound)
824    ///\param e is a linear expression (see \ref Expr)
825    ///\param u is the upper bound (\ref INF means no bound)
826    ///\return The created row.
827    ///\bug This is a temportary function. The interface will change to
828    ///a better one.
829    Row addRow(Value l,const Expr &e, Value u) {
830      Row r=addRow();
831      setRow(r,l,e,u);
832      return r;
833    }
834
835    ///Add a new row (i.e a new constraint) to the LP
836
837    ///\param c is a linear expression (see \ref Constr)
838    ///\return The created row.
839    Row addRow(const Constr &c) {
840      Row r=addRow();
841      setRow(r,c);
842      return r;
843    }
844    ///Erase a coloumn (i.e a variable) from the LP
845
846    ///\param c is the coloumn to be deleted
847    ///\todo Please check this
848    void eraseCol(Col c) {
849      _eraseCol(cols.floatingId(c.id));
850      cols.erase(c.id);
851    }
852    ///Erase a  row (i.e a constraint) from the LP
853
854    ///\param r is the row to be deleted
855    ///\todo Please check this
856    void eraseRow(Row r) {
857      _eraseRow(rows.floatingId(r.id));
858      rows.erase(r.id);
859    }
860
861    ///Set an element of the coefficient matrix of the LP
862
863    ///\param r is the row of the element to be modified
864    ///\param c is the coloumn of the element to be modified
865    ///\param val is the new value of the coefficient
866    void setCoeff(Row r, Col c, Value val){
867      _setCoeff(rows.floatingId(r.id),cols.floatingId(c.id), val);
868    }
869
870    /// Set the lower bound of a column (i.e a variable)
871
872    /// The upper bound of a variable (column) has to be given by an
873    /// extended number of type Value, i.e. a finite number of type
874    /// Value or -\ref INF.
875    void colLowerBound(Col c, Value value) {
876      _setColLowerBound(cols.floatingId(c.id),value);
877    }
878    /// Set the upper bound of a column (i.e a variable)
879
880    /// The upper bound of a variable (column) has to be given by an
881    /// extended number of type Value, i.e. a finite number of type
882    /// Value or \ref INF.
883    void colUpperBound(Col c, Value value) {
884      _setColUpperBound(cols.floatingId(c.id),value);
885    };
886    /// Set the lower and the upper bounds of a column (i.e a variable)
887
888    /// The lower and the upper bounds of
889    /// a variable (column) have to be given by an
890    /// extended number of type Value, i.e. a finite number of type
891    /// Value, -\ref INF or \ref INF.
892    void colBounds(Col c, Value lower, Value upper) {
893      _setColLowerBound(cols.floatingId(c.id),lower);
894      _setColUpperBound(cols.floatingId(c.id),upper);
895    }
896   
897//     /// Set the lower bound of a row (i.e a constraint)
898
899//     /// The lower bound of a linear expression (row) has to be given by an
900//     /// extended number of type Value, i.e. a finite number of type
901//     /// Value or -\ref INF.
902//     void rowLowerBound(Row r, Value value) {
903//       _setRowLowerBound(rows.floatingId(r.id),value);
904//     };
905//     /// Set the upper bound of a row (i.e a constraint)
906
907//     /// The upper bound of a linear expression (row) has to be given by an
908//     /// extended number of type Value, i.e. a finite number of type
909//     /// Value or \ref INF.
910//     void rowUpperBound(Row r, Value value) {
911//       _setRowUpperBound(rows.floatingId(r.id),value);
912//     };
913
914    /// Set the lower and the upper bounds of a row (i.e a constraint)
915
916    /// The lower and the upper bounds of
917    /// a constraint (row) have to be given by an
918    /// extended number of type Value, i.e. a finite number of type
919    /// Value, -\ref INF or \ref INF.
920    void rowBounds(Row c, Value lower, Value upper) {
921      _setRowBounds(rows.floatingId(c.id),lower, upper);
922      // _setRowUpperBound(rows.floatingId(c.id),upper);
923    }
924   
925    ///Set an element of the objective function
926    void objCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
927    ///Set the objective function
928   
929    ///\param e is a linear expression of type \ref Expr.
930    ///\bug The previous objective function is not cleared!
931    void setObj(Expr e) {
932      _clearObj();
933      for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
934        objCoeff((*i).first,(*i).second);
935      obj_const_comp=e.constComp();
936    }
937
938    ///Maximize
939    void max() { _setMax(); }
940    ///Minimize
941    void min() { _setMin(); }
942
943   
944    ///@}
945
946
947    ///\name Solve the LP
948
949    ///@{
950
951    ///\e Solve the LP problem at hand
952    ///
953    ///\return The result of the optimization procedure. Possible values and their meanings can be found in the documentation of \ref SolveExitStatus.
954    ///
955    ///\todo Which method is used to solve the problem
956    SolveExitStatus solve() { return _solve(); }
957   
958    ///@}
959   
960    ///\name Obtain the solution
961
962    ///@{
963
964    /// The status of the primal problem (the original LP problem)
965    SolutionStatus primalStatus() {
966      return _getPrimalStatus();
967    }
968
969    /// The status of the dual (of the original LP) problem
970    SolutionStatus dualStatus() {
971      return _getDualStatus();
972    }
973
974    ///The type of the original LP problem
975    ProblemTypes problemType() {
976      return _getProblemType();
977    }
978
979    ///\e
980    Value primal(Col c) { return _getPrimal(cols.floatingId(c.id)); }
981
982    ///\e
983
984    ///\return
985    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
986    /// of the primal problem, depending on whether we minimize or maximize.
987    ///- \ref NaN if no primal solution is found.
988    ///- The (finite) objective value if an optimal solution is found.
989    Value primalValue() { return _getPrimalValue()+obj_const_comp;}
990    ///@}
991   
992  }; 
993
994  ///\e
995 
996  ///\relates LpSolverBase::Expr
997  ///
998  inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
999                                      const LpSolverBase::Expr &b)
1000  {
1001    LpSolverBase::Expr tmp(a);
1002    tmp+=b;
1003    return tmp;
1004  }
1005  ///\e
1006 
1007  ///\relates LpSolverBase::Expr
1008  ///
1009  inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
1010                                      const LpSolverBase::Expr &b)
1011  {
1012    LpSolverBase::Expr tmp(a);
1013    tmp-=b;
1014    return tmp;
1015  }
1016  ///\e
1017 
1018  ///\relates LpSolverBase::Expr
1019  ///
1020  inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
1021                                      const LpSolverBase::Value &b)
1022  {
1023    LpSolverBase::Expr tmp(a);
1024    tmp*=b;
1025    return tmp;
1026  }
1027 
1028  ///\e
1029 
1030  ///\relates LpSolverBase::Expr
1031  ///
1032  inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
1033                                      const LpSolverBase::Expr &b)
1034  {
1035    LpSolverBase::Expr tmp(b);
1036    tmp*=a;
1037    return tmp;
1038  }
1039  ///\e
1040 
1041  ///\relates LpSolverBase::Expr
1042  ///
1043  inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
1044                                      const LpSolverBase::Value &b)
1045  {
1046    LpSolverBase::Expr tmp(a);
1047    tmp/=b;
1048    return tmp;
1049  }
1050 
1051  ///\e
1052 
1053  ///\relates LpSolverBase::Constr
1054  ///
1055  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1056                                         const LpSolverBase::Expr &f)
1057  {
1058    return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
1059  }
1060
1061  ///\e
1062 
1063  ///\relates LpSolverBase::Constr
1064  ///
1065  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
1066                                         const LpSolverBase::Expr &f)
1067  {
1068    return LpSolverBase::Constr(e,f);
1069  }
1070
1071  ///\e
1072 
1073  ///\relates LpSolverBase::Constr
1074  ///
1075  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
1076                                         const LpSolverBase::Value &f)
1077  {
1078    return LpSolverBase::Constr(e,f);
1079  }
1080
1081  ///\e
1082 
1083  ///\relates LpSolverBase::Constr
1084  ///
1085  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1086                                         const LpSolverBase::Expr &f)
1087  {
1088    return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
1089  }
1090
1091
1092  ///\e
1093 
1094  ///\relates LpSolverBase::Constr
1095  ///
1096  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
1097                                         const LpSolverBase::Expr &f)
1098  {
1099    return LpSolverBase::Constr(f,e);
1100  }
1101
1102
1103  ///\e
1104 
1105  ///\relates LpSolverBase::Constr
1106  ///
1107  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
1108                                         const LpSolverBase::Value &f)
1109  {
1110    return LpSolverBase::Constr(f,e);
1111  }
1112
1113  ///\e
1114 
1115  ///\relates LpSolverBase::Constr
1116  ///
1117  inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
1118                                         const LpSolverBase::Expr &f)
1119  {
1120    return LpSolverBase::Constr(0,e-f,0);
1121  }
1122
1123  ///\e
1124 
1125  ///\relates LpSolverBase::Constr
1126  ///
1127  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
1128                                         const LpSolverBase::Constr&c)
1129  {
1130    LpSolverBase::Constr tmp(c);
1131    ///\todo Create an own exception type.
1132    if(!isnan(tmp.lowerBound())) throw LogicError();
1133    else tmp.lowerBound()=n;
1134    return tmp;
1135  }
1136  ///\e
1137 
1138  ///\relates LpSolverBase::Constr
1139  ///
1140  inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
1141                                         const LpSolverBase::Value &n)
1142  {
1143    LpSolverBase::Constr tmp(c);
1144    ///\todo Create an own exception type.
1145    if(!isnan(tmp.upperBound())) throw LogicError();
1146    else tmp.upperBound()=n;
1147    return tmp;
1148  }
1149
1150  ///\e
1151 
1152  ///\relates LpSolverBase::Constr
1153  ///
1154  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
1155                                         const LpSolverBase::Constr&c)
1156  {
1157    LpSolverBase::Constr tmp(c);
1158    ///\todo Create an own exception type.
1159    if(!isnan(tmp.upperBound())) throw LogicError();
1160    else tmp.upperBound()=n;
1161    return tmp;
1162  }
1163  ///\e
1164 
1165  ///\relates LpSolverBase::Constr
1166  ///
1167  inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
1168                                         const LpSolverBase::Value &n)
1169  {
1170    LpSolverBase::Constr tmp(c);
1171    ///\todo Create an own exception type.
1172    if(!isnan(tmp.lowerBound())) throw LogicError();
1173    else tmp.lowerBound()=n;
1174    return tmp;
1175  }
1176
1177  ///\e
1178 
1179  ///\relates LpSolverBase::DualExpr
1180  ///
1181  inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
1182                                      const LpSolverBase::DualExpr &b)
1183  {
1184    LpSolverBase::DualExpr tmp(a);
1185    tmp+=b;
1186    return tmp;
1187  }
1188  ///\e
1189 
1190  ///\relates LpSolverBase::DualExpr
1191  ///
1192  inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
1193                                      const LpSolverBase::DualExpr &b)
1194  {
1195    LpSolverBase::DualExpr tmp(a);
1196    tmp-=b;
1197    return tmp;
1198  }
1199  ///\e
1200 
1201  ///\relates LpSolverBase::DualExpr
1202  ///
1203  inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
1204                                      const LpSolverBase::Value &b)
1205  {
1206    LpSolverBase::DualExpr tmp(a);
1207    tmp*=b;
1208    return tmp;
1209  }
1210 
1211  ///\e
1212 
1213  ///\relates LpSolverBase::DualExpr
1214  ///
1215  inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
1216                                      const LpSolverBase::DualExpr &b)
1217  {
1218    LpSolverBase::DualExpr tmp(b);
1219    tmp*=a;
1220    return tmp;
1221  }
1222  ///\e
1223 
1224  ///\relates LpSolverBase::DualExpr
1225  ///
1226  inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
1227                                      const LpSolverBase::Value &b)
1228  {
1229    LpSolverBase::DualExpr tmp(a);
1230    tmp/=b;
1231    return tmp;
1232  }
1233 
1234
1235} //namespace lemon
1236
1237#endif //LEMON_LP_BASE_H
Note: See TracBrowser for help on using the repository browser.