COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/lp_base.h @ 2028:d0e8a86a1ff2

Last change on this file since 2028:d0e8a86a1ff2 was 2026:8d49961ec50f, checked in by Balazs Dezso, 14 years ago

NaN checking to be conform to MinGW32

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