COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/lp_base.h @ 2151:38ec4a930c05

Last change on this file since 2151:38ec4a930c05 was 2148:ab368e0ab662, checked in by athos, 13 years ago

Modifications to the interface: colType() functions, though I left the old integer() functions, too.

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