COIN-OR::LEMON - Graph Library

source: lemon/lemon/lp_base.h @ 1369:9fd86ec2cb81

Last change on this file since 1369:9fd86ec2cb81 was 1353:760a5f690163, checked in by Alpar Juttner <alpar@…>, 4 years ago

Minor doc fixes

File size: 64.1 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2013
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<iostream>
23#include<vector>
24#include<map>
25#include<limits>
26#include<lemon/math.h>
27
28#include<lemon/error.h>
29#include<lemon/assert.h>
30
31#include<lemon/core.h>
32#include<lemon/bits/solver_bits.h>
33
34#include<lemon/bits/stl_iterators.h>
35
36///\file
37///\brief The interface of the LP solver interface.
38///\ingroup lp_group
39namespace lemon {
40
41  ///Common base class for LP and MIP solvers
42
43  ///Usually this class is not used directly, please use one of the concrete
44  ///implementations of the solver interface.
45  ///\ingroup lp_group
46  class LpBase {
47
48  protected:
49
50    _solver_bits::VarIndex _rows;
51    _solver_bits::VarIndex _cols;
52
53  public:
54
55    ///Possible outcomes of an LP solving procedure
56    enum SolveExitStatus {
57      /// = 0. It means that the problem has been successfully solved: either
58      ///an optimal solution has been found or infeasibility/unboundedness
59      ///has been proved.
60      SOLVED = 0,
61      /// = 1. Any other case (including the case when some user specified
62      ///limit has been exceeded).
63      UNSOLVED = 1
64    };
65
66    ///Direction of the optimization
67    enum Sense {
68      /// Minimization
69      MIN,
70      /// Maximization
71      MAX
72    };
73
74    ///Enum for \c messageLevel() parameter
75    enum MessageLevel {
76      /// No output (default value).
77      MESSAGE_NOTHING,
78      /// Error messages only.
79      MESSAGE_ERROR,
80      /// Warnings.
81      MESSAGE_WARNING,
82      /// Normal output.
83      MESSAGE_NORMAL,
84      /// Verbose output.
85      MESSAGE_VERBOSE
86    };
87
88
89    ///The floating point type used by the solver
90    typedef double Value;
91    ///The infinity constant
92    static const Value INF;
93    ///The not a number constant
94    static const Value NaN;
95
96    friend class Col;
97    friend class ColIt;
98    friend class Row;
99    friend class RowIt;
100
101    ///Refer to a column of the LP.
102
103    ///This type is used to refer to a column of the LP.
104    ///
105    ///Its value remains valid and correct even after the addition or erase of
106    ///other columns.
107    ///
108    ///\note This class is similar to other Item types in LEMON, like
109    ///Node and Arc types in digraph.
110    class Col {
111      friend class LpBase;
112    protected:
113      int _id;
114      explicit Col(int id) : _id(id) {}
115    public:
116      typedef Value ExprValue;
117      typedef True LpCol;
118      /// Default constructor
119
120      /// \warning The default constructor sets the Col to an
121      /// undefined value.
122      Col() {}
123      /// Invalid constructor \& conversion.
124
125      /// This constructor initializes the Col to be invalid.
126      /// \sa Invalid for more details.
127      Col(const Invalid&) : _id(-1) {}
128      /// Equality operator
129
130      /// Two \ref Col "Col"s are equal if and only if they point to
131      /// the same LP column or both are invalid.
132      bool operator==(Col c) const  {return _id == c._id;}
133      /// Inequality operator
134
135      /// \sa operator==(Col c)
136      ///
137      bool operator!=(Col c) const  {return _id != c._id;}
138      /// Artificial ordering operator.
139
140      /// To allow the use of this object in std::map or similar
141      /// associative container we require this.
142      ///
143      /// \note This operator only have to define some strict ordering of
144      /// the items; this order has nothing to do with the iteration
145      /// ordering of the items.
146      bool operator<(Col c) const  {return _id < c._id;}
147    };
148
149    ///Iterator for iterate over the columns of an LP problem
150
151    /// Its usage is quite simple, for example, you can count the number
152    /// of columns in an LP \c lp:
153    ///\code
154    /// int count=0;
155    /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
156    ///\endcode
157    class ColIt : public Col {
158      const LpBase *_solver;
159    public:
160      /// Default constructor
161
162      /// \warning The default constructor sets the iterator
163      /// to an undefined value.
164      ColIt() {}
165      /// Sets the iterator to the first Col
166
167      /// Sets the iterator to the first Col.
168      ///
169      ColIt(const LpBase &solver) : _solver(&solver)
170      {
171        _solver->_cols.firstItem(_id);
172      }
173      /// Invalid constructor \& conversion
174
175      /// Initialize the iterator to be invalid.
176      /// \sa Invalid for more details.
177      ColIt(const Invalid&) : Col(INVALID) {}
178      /// Next column
179
180      /// Assign the iterator to the next column.
181      ///
182      ColIt &operator++()
183      {
184        _solver->_cols.nextItem(_id);
185        return *this;
186      }
187    };
188
189    /// \brief Gets the collection of the columns of the LP problem.
190    ///
191    /// This function can be used for iterating on
192    /// the columns of the LP problem. It returns a wrapped ColIt, which looks
193    /// like an STL container (by having begin() and end())
194    /// which you can use in range-based for loops, STL algorithms, etc.
195    /// For example you can write:
196    ///\code
197    /// for(auto c: lp.cols())
198    ///   doSomething(c);
199    ///\endcode
200    LemonRangeWrapper1<ColIt, LpBase> cols() {
201      return LemonRangeWrapper1<ColIt, LpBase>(*this);
202    }
203
204
205    /// \brief Returns the ID of the column.
206    static int id(const Col& col) { return col._id; }
207    /// \brief Returns the column with the given ID.
208    ///
209    /// \pre The argument should be a valid column ID in the LP problem.
210    static Col colFromId(int id) { return Col(id); }
211
212    ///Refer to a row of the LP.
213
214    ///This type is used to refer to a row of the LP.
215    ///
216    ///Its value remains valid and correct even after the addition or erase of
217    ///other rows.
218    ///
219    ///\note This class is similar to other Item types in LEMON, like
220    ///Node and Arc types in digraph.
221    class Row {
222      friend class LpBase;
223    protected:
224      int _id;
225      explicit Row(int id) : _id(id) {}
226    public:
227      typedef Value ExprValue;
228      typedef True LpRow;
229      /// Default constructor
230
231      /// \warning The default constructor sets the Row to an
232      /// undefined value.
233      Row() {}
234      /// Invalid constructor \& conversion.
235
236      /// This constructor initializes the Row to be invalid.
237      /// \sa Invalid for more details.
238      Row(const Invalid&) : _id(-1) {}
239      /// Equality operator
240
241      /// Two \ref Row "Row"s are equal if and only if they point to
242      /// the same LP row or both are invalid.
243      bool operator==(Row r) const  {return _id == r._id;}
244      /// Inequality operator
245
246      /// \sa operator==(Row r)
247      ///
248      bool operator!=(Row r) const  {return _id != r._id;}
249      /// Artificial ordering operator.
250
251      /// To allow the use of this object in std::map or similar
252      /// associative container we require this.
253      ///
254      /// \note This operator only have to define some strict ordering of
255      /// the items; this order has nothing to do with the iteration
256      /// ordering of the items.
257      bool operator<(Row r) const  {return _id < r._id;}
258    };
259
260    ///Iterator for iterate over the rows of an LP problem
261
262    /// Its usage is quite simple, for example, you can count the number
263    /// of rows in an LP \c lp:
264    ///\code
265    /// int count=0;
266    /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
267    ///\endcode
268    class RowIt : public Row {
269      const LpBase *_solver;
270    public:
271      /// Default constructor
272
273      /// \warning The default constructor sets the iterator
274      /// to an undefined value.
275      RowIt() {}
276      /// Sets the iterator to the first Row
277
278      /// Sets the iterator to the first Row.
279      ///
280      RowIt(const LpBase &solver) : _solver(&solver)
281      {
282        _solver->_rows.firstItem(_id);
283      }
284      /// Invalid constructor \& conversion
285
286      /// Initialize the iterator to be invalid.
287      /// \sa Invalid for more details.
288      RowIt(const Invalid&) : Row(INVALID) {}
289      /// Next row
290
291      /// Assign the iterator to the next row.
292      ///
293      RowIt &operator++()
294      {
295        _solver->_rows.nextItem(_id);
296        return *this;
297      }
298    };
299   
300    /// \brief Gets the collection of the rows of the LP problem.
301    ///
302    /// This function can be used for iterating on
303    /// the rows of the LP problem. It returns a wrapped RowIt, which looks
304    /// like an STL container (by having begin() and end())
305    /// which you can use in range-based for loops, STL algorithms, etc.
306    /// For example you can write:
307    ///\code
308    /// for(auto c: lp.rows())
309    ///   doSomething(c);
310    ///\endcode
311    LemonRangeWrapper1<RowIt, LpBase> rows() {
312      return LemonRangeWrapper1<RowIt, LpBase>(*this);
313    }
314   
315
316    /// \brief Returns the ID of the row.
317    static int id(const Row& row) { return row._id; }
318    /// \brief Returns the row with the given ID.
319    ///
320    /// \pre The argument should be a valid row ID in the LP problem.
321    static Row rowFromId(int id) { return Row(id); }
322
323  public:
324
325    ///Linear expression of variables and a constant component
326
327    ///This data structure stores a linear expression of the variables
328    ///(\ref Col "Col"s) and also has a constant component.
329    ///
330    ///There are several ways to access and modify the contents of this
331    ///container.
332    ///\code
333    ///e[v]=5;
334    ///e[v]+=12;
335    ///e.erase(v);
336    ///\endcode
337    ///or you can also iterate through its elements.
338    ///\code
339    ///double s=0;
340    ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
341    ///  s+=*i * primal(i);
342    ///\endcode
343    ///(This code computes the primal value of the expression).
344    ///- Numbers (<tt>double</tt>'s)
345    ///and variables (\ref Col "Col"s) directly convert to an
346    ///\ref Expr and the usual linear operations are defined, so
347    ///\code
348    ///v+w
349    ///2*v-3.12*(v-w/2)+2
350    ///v*2.1+(3*v+(v*12+w+6)*3)/2
351    ///\endcode
352    ///are valid expressions.
353    ///The usual assignment operations are also defined.
354    ///\code
355    ///e=v+w;
356    ///e+=2*v-3.12*(v-w/2)+2;
357    ///e*=3.4;
358    ///e/=5;
359    ///\endcode
360    ///- The constant member can be set and read by dereference
361    ///  operator (unary *)
362    ///
363    ///\code
364    ///*e=12;
365    ///double c=*e;
366    ///\endcode
367    ///
368    ///\sa Constr
369    class Expr {
370      friend class LpBase;
371    public:
372      /// The key type of the expression
373      typedef LpBase::Col Key;
374      /// The value type of the expression
375      typedef LpBase::Value Value;
376
377    protected:
378      Value const_comp;
379      std::map<int, Value> comps;
380
381    public:
382      typedef True SolverExpr;
383      /// Default constructor
384
385      /// Construct an empty expression, the coefficients and
386      /// the constant component are initialized to zero.
387      Expr() : const_comp(0) {}
388      /// Construct an expression from a column
389
390      /// Construct an expression, which has a term with \c c variable
391      /// and 1.0 coefficient.
392      Expr(const Col &c) : const_comp(0) {
393        typedef std::map<int, Value>::value_type pair_type;
394        comps.insert(pair_type(id(c), 1));
395      }
396      /// Construct an expression from a constant
397
398      /// Construct an expression, which's constant component is \c v.
399      ///
400      Expr(const Value &v) : const_comp(v) {}
401      /// Returns the coefficient of the column
402      Value operator[](const Col& c) const {
403        std::map<int, Value>::const_iterator it=comps.find(id(c));
404        if (it != comps.end()) {
405          return it->second;
406        } else {
407          return 0;
408        }
409      }
410      /// Returns the coefficient of the column
411      Value& operator[](const Col& c) {
412        return comps[id(c)];
413      }
414      /// Sets the coefficient of the column
415      void set(const Col &c, const Value &v) {
416        if (v != 0.0) {
417          typedef std::map<int, Value>::value_type pair_type;
418          comps.insert(pair_type(id(c), v));
419        } else {
420          comps.erase(id(c));
421        }
422      }
423      /// Returns the constant component of the expression
424      Value& operator*() { return const_comp; }
425      /// Returns the constant component of the expression
426      const Value& operator*() const { return const_comp; }
427      /// \brief Removes the coefficients which's absolute value does
428      /// not exceed \c epsilon. It also sets to zero the constant
429      /// component, if it does not exceed epsilon in absolute value.
430      void simplify(Value epsilon = 0.0) {
431        std::map<int, Value>::iterator it=comps.begin();
432        while (it != comps.end()) {
433          std::map<int, Value>::iterator jt=it;
434          ++jt;
435          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
436          it=jt;
437        }
438        if (std::fabs(const_comp) <= epsilon) const_comp = 0;
439      }
440
441      void simplify(Value epsilon = 0.0) const {
442        const_cast<Expr*>(this)->simplify(epsilon);
443      }
444
445      ///Sets all coefficients and the constant component to 0.
446      void clear() {
447        comps.clear();
448        const_comp=0;
449      }
450
451      ///Compound assignment
452      Expr &operator+=(const Expr &e) {
453        for (std::map<int, Value>::const_iterator it=e.comps.begin();
454             it!=e.comps.end(); ++it)
455          comps[it->first]+=it->second;
456        const_comp+=e.const_comp;
457        return *this;
458      }
459      ///Compound assignment
460      Expr &operator-=(const Expr &e) {
461        for (std::map<int, Value>::const_iterator it=e.comps.begin();
462             it!=e.comps.end(); ++it)
463          comps[it->first]-=it->second;
464        const_comp-=e.const_comp;
465        return *this;
466      }
467      ///Multiply with a constant
468      Expr &operator*=(const Value &v) {
469        for (std::map<int, Value>::iterator it=comps.begin();
470             it!=comps.end(); ++it)
471          it->second*=v;
472        const_comp*=v;
473        return *this;
474      }
475      ///Division with a constant
476      Expr &operator/=(const Value &c) {
477        for (std::map<int, Value>::iterator it=comps.begin();
478             it!=comps.end(); ++it)
479          it->second/=c;
480        const_comp/=c;
481        return *this;
482      }
483
484      ///Iterator over the expression
485
486      ///The iterator iterates over the terms of the expression.
487      ///
488      ///\code
489      ///double s=0;
490      ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
491      ///  s+= *i * primal(i);
492      ///\endcode
493      class CoeffIt {
494      private:
495
496        std::map<int, Value>::iterator _it, _end;
497
498      public:
499
500        /// Sets the iterator to the first term
501
502        /// Sets the iterator to the first term of the expression.
503        ///
504        CoeffIt(Expr& e)
505          : _it(e.comps.begin()), _end(e.comps.end()){}
506
507        /// Convert the iterator to the column of the term
508        operator Col() const {
509          return colFromId(_it->first);
510        }
511
512        /// Returns the coefficient of the term
513        Value& operator*() { return _it->second; }
514
515        /// Returns the coefficient of the term
516        const Value& operator*() const { return _it->second; }
517        /// Next term
518
519        /// Assign the iterator to the next term.
520        ///
521        CoeffIt& operator++() { ++_it; return *this; }
522
523        /// Equality operator
524        bool operator==(Invalid) const { return _it == _end; }
525        /// Inequality operator
526        bool operator!=(Invalid) const { return _it != _end; }
527      };
528
529      /// Const iterator over the expression
530
531      ///The iterator iterates over the terms of the expression.
532      ///
533      ///\code
534      ///double s=0;
535      ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
536      ///  s+=*i * primal(i);
537      ///\endcode
538      class ConstCoeffIt {
539      private:
540
541        std::map<int, Value>::const_iterator _it, _end;
542
543      public:
544
545        /// Sets the iterator to the first term
546
547        /// Sets the iterator to the first term of the expression.
548        ///
549        ConstCoeffIt(const Expr& e)
550          : _it(e.comps.begin()), _end(e.comps.end()){}
551
552        /// Convert the iterator to the column of the term
553        operator Col() const {
554          return colFromId(_it->first);
555        }
556
557        /// Returns the coefficient of the term
558        const Value& operator*() const { return _it->second; }
559
560        /// Next term
561
562        /// Assign the iterator to the next term.
563        ///
564        ConstCoeffIt& operator++() { ++_it; return *this; }
565
566        /// Equality operator
567        bool operator==(Invalid) const { return _it == _end; }
568        /// Inequality operator
569        bool operator!=(Invalid) const { return _it != _end; }
570      };
571
572    };
573
574    ///Linear constraint
575
576    ///This data stucture represents a linear constraint in the LP.
577    ///Basically it is a linear expression with a lower or an upper bound
578    ///(or both). These parts of the constraint can be obtained by the member
579    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
580    ///respectively.
581    ///There are two ways to construct a constraint.
582    ///- You can set the linear expression and the bounds directly
583    ///  by the functions above.
584    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
585    ///  are defined between expressions, or even between constraints whenever
586    ///  it makes sense. Therefore if \c e and \c f are linear expressions and
587    ///  \c s and \c t are numbers, then the followings are valid expressions
588    ///  and thus they can be used directly e.g. in \ref addRow() whenever
589    ///  it makes sense.
590    ///\code
591    ///  e<=s
592    ///  e<=f
593    ///  e==f
594    ///  s<=e<=t
595    ///  e>=t
596    ///\endcode
597    ///\warning The validity of a constraint is checked only at run
598    ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
599    ///compile, but will fail an assertion.
600    class Constr
601    {
602    public:
603      typedef LpBase::Expr Expr;
604      typedef Expr::Key Key;
605      typedef Expr::Value Value;
606
607    protected:
608      Expr _expr;
609      Value _lb,_ub;
610    public:
611      ///\e
612      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
613      ///\e
614      Constr(Value lb, const Expr &e, Value ub) :
615        _expr(e), _lb(lb), _ub(ub) {}
616      Constr(const Expr &e) :
617        _expr(e), _lb(NaN), _ub(NaN) {}
618      ///\e
619      void clear()
620      {
621        _expr.clear();
622        _lb=_ub=NaN;
623      }
624
625      ///Reference to the linear expression
626      Expr &expr() { return _expr; }
627      ///Cont reference to the linear expression
628      const Expr &expr() const { return _expr; }
629      ///Reference to the lower bound.
630
631      ///\return
632      ///- \ref INF "INF": the constraint is lower unbounded.
633      ///- \ref NaN "NaN": lower bound has not been set.
634      ///- finite number: the lower bound
635      Value &lowerBound() { return _lb; }
636      ///The const version of \ref lowerBound()
637      const Value &lowerBound() const { return _lb; }
638      ///Reference to the upper bound.
639
640      ///\return
641      ///- \ref INF "INF": the constraint is upper unbounded.
642      ///- \ref NaN "NaN": upper bound has not been set.
643      ///- finite number: the upper bound
644      Value &upperBound() { return _ub; }
645      ///The const version of \ref upperBound()
646      const Value &upperBound() const { return _ub; }
647      ///Is the constraint lower bounded?
648      bool lowerBounded() const {
649        return _lb != -INF && !isNaN(_lb);
650      }
651      ///Is the constraint upper bounded?
652      bool upperBounded() const {
653        return _ub != INF && !isNaN(_ub);
654      }
655
656    };
657
658    ///Linear expression of rows
659
660    ///This data structure represents a column of the matrix,
661    ///thas is it strores a linear expression of the dual variables
662    ///(\ref LpBase::Row "Row"s).
663    ///
664    ///There are several ways to access and modify the contents of this
665    ///container.
666    ///\code
667    ///e[v]=5;
668    ///e[v]+=12;
669    ///e.erase(v);
670    ///\endcode
671    ///or you can also iterate through its elements.
672    ///\code
673    ///double s=0;
674    ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
675    ///  s+=*i;
676    ///\endcode
677    ///(This code computes the sum of all coefficients).
678    ///- Numbers (<tt>double</tt>'s)
679    ///and variables (\ref LpBase::Row "Row"s) directly convert to an
680    ///\ref DualExpr and the usual linear operations are defined, so
681    ///\code
682    ///v+w
683    ///2*v-3.12*(v-w/2)
684    ///v*2.1+(3*v+(v*12+w)*3)/2
685    ///\endcode
686    ///are valid \ref DualExpr dual expressions.
687    ///The usual assignment operations are also defined.
688    ///\code
689    ///e=v+w;
690    ///e+=2*v-3.12*(v-w/2);
691    ///e*=3.4;
692    ///e/=5;
693    ///\endcode
694    ///
695    ///\sa Expr
696    class DualExpr {
697      friend class LpBase;
698    public:
699      /// The key type of the expression
700      typedef LpBase::Row Key;
701      /// The value type of the expression
702      typedef LpBase::Value Value;
703
704    protected:
705      std::map<int, Value> comps;
706
707    public:
708      typedef True SolverExpr;
709      /// Default constructor
710
711      /// Construct an empty expression, the coefficients are
712      /// initialized to zero.
713      DualExpr() {}
714      /// Construct an expression from a row
715
716      /// Construct an expression, which has a term with \c r dual
717      /// variable and 1.0 coefficient.
718      DualExpr(const Row &r) {
719        typedef std::map<int, Value>::value_type pair_type;
720        comps.insert(pair_type(id(r), 1));
721      }
722      /// Returns the coefficient of the row
723      Value operator[](const Row& r) const {
724        std::map<int, Value>::const_iterator it = comps.find(id(r));
725        if (it != comps.end()) {
726          return it->second;
727        } else {
728          return 0;
729        }
730      }
731      /// Returns the coefficient of the row
732      Value& operator[](const Row& r) {
733        return comps[id(r)];
734      }
735      /// Sets the coefficient of the row
736      void set(const Row &r, const Value &v) {
737        if (v != 0.0) {
738          typedef std::map<int, Value>::value_type pair_type;
739          comps.insert(pair_type(id(r), v));
740        } else {
741          comps.erase(id(r));
742        }
743      }
744      /// \brief Removes the coefficients which's absolute value does
745      /// not exceed \c epsilon.
746      void simplify(Value epsilon = 0.0) {
747        std::map<int, Value>::iterator it=comps.begin();
748        while (it != comps.end()) {
749          std::map<int, Value>::iterator jt=it;
750          ++jt;
751          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
752          it=jt;
753        }
754      }
755
756      void simplify(Value epsilon = 0.0) const {
757        const_cast<DualExpr*>(this)->simplify(epsilon);
758      }
759
760      ///Sets all coefficients to 0.
761      void clear() {
762        comps.clear();
763      }
764      ///Compound assignment
765      DualExpr &operator+=(const DualExpr &e) {
766        for (std::map<int, Value>::const_iterator it=e.comps.begin();
767             it!=e.comps.end(); ++it)
768          comps[it->first]+=it->second;
769        return *this;
770      }
771      ///Compound assignment
772      DualExpr &operator-=(const DualExpr &e) {
773        for (std::map<int, Value>::const_iterator it=e.comps.begin();
774             it!=e.comps.end(); ++it)
775          comps[it->first]-=it->second;
776        return *this;
777      }
778      ///Multiply with a constant
779      DualExpr &operator*=(const Value &v) {
780        for (std::map<int, Value>::iterator it=comps.begin();
781             it!=comps.end(); ++it)
782          it->second*=v;
783        return *this;
784      }
785      ///Division with a constant
786      DualExpr &operator/=(const Value &v) {
787        for (std::map<int, Value>::iterator it=comps.begin();
788             it!=comps.end(); ++it)
789          it->second/=v;
790        return *this;
791      }
792
793      ///Iterator over the expression
794
795      ///The iterator iterates over the terms of the expression.
796      ///
797      ///\code
798      ///double s=0;
799      ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
800      ///  s+= *i * dual(i);
801      ///\endcode
802      class CoeffIt {
803      private:
804
805        std::map<int, Value>::iterator _it, _end;
806
807      public:
808
809        /// Sets the iterator to the first term
810
811        /// Sets the iterator to the first term of the expression.
812        ///
813        CoeffIt(DualExpr& e)
814          : _it(e.comps.begin()), _end(e.comps.end()){}
815
816        /// Convert the iterator to the row of the term
817        operator Row() const {
818          return rowFromId(_it->first);
819        }
820
821        /// Returns the coefficient of the term
822        Value& operator*() { return _it->second; }
823
824        /// Returns the coefficient of the term
825        const Value& operator*() const { return _it->second; }
826
827        /// Next term
828
829        /// Assign the iterator to the next term.
830        ///
831        CoeffIt& operator++() { ++_it; return *this; }
832
833        /// Equality operator
834        bool operator==(Invalid) const { return _it == _end; }
835        /// Inequality operator
836        bool operator!=(Invalid) const { return _it != _end; }
837      };
838
839      ///Iterator over the expression
840
841      ///The iterator iterates over the terms of the expression.
842      ///
843      ///\code
844      ///double s=0;
845      ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
846      ///  s+= *i * dual(i);
847      ///\endcode
848      class ConstCoeffIt {
849      private:
850
851        std::map<int, Value>::const_iterator _it, _end;
852
853      public:
854
855        /// Sets the iterator to the first term
856
857        /// Sets the iterator to the first term of the expression.
858        ///
859        ConstCoeffIt(const DualExpr& e)
860          : _it(e.comps.begin()), _end(e.comps.end()){}
861
862        /// Convert the iterator to the row of the term
863        operator Row() const {
864          return rowFromId(_it->first);
865        }
866
867        /// Returns the coefficient of the term
868        const Value& operator*() const { return _it->second; }
869
870        /// Next term
871
872        /// Assign the iterator to the next term.
873        ///
874        ConstCoeffIt& operator++() { ++_it; return *this; }
875
876        /// Equality operator
877        bool operator==(Invalid) const { return _it == _end; }
878        /// Inequality operator
879        bool operator!=(Invalid) const { return _it != _end; }
880      };
881    };
882
883
884  protected:
885
886    class InsertIterator {
887    private:
888
889      std::map<int, Value>& _host;
890      const _solver_bits::VarIndex& _index;
891
892    public:
893
894      typedef std::output_iterator_tag iterator_category;
895      typedef void difference_type;
896      typedef void value_type;
897      typedef void reference;
898      typedef void pointer;
899
900      InsertIterator(std::map<int, Value>& host,
901                   const _solver_bits::VarIndex& index)
902        : _host(host), _index(index) {}
903
904      InsertIterator& operator=(const std::pair<int, Value>& value) {
905        typedef std::map<int, Value>::value_type pair_type;
906        _host.insert(pair_type(_index[value.first], value.second));
907        return *this;
908      }
909
910      InsertIterator& operator*() { return *this; }
911      InsertIterator& operator++() { return *this; }
912      InsertIterator operator++(int) { return *this; }
913
914    };
915
916    class ExprIterator {
917    private:
918      std::map<int, Value>::const_iterator _host_it;
919      const _solver_bits::VarIndex& _index;
920    public:
921
922      typedef std::bidirectional_iterator_tag iterator_category;
923      typedef std::ptrdiff_t difference_type;
924      typedef const std::pair<int, Value> value_type;
925      typedef value_type reference;
926
927      class pointer {
928      public:
929        pointer(value_type& _value) : value(_value) {}
930        value_type* operator->() { return &value; }
931      private:
932        value_type value;
933      };
934
935      ExprIterator(const std::map<int, Value>::const_iterator& host_it,
936                   const _solver_bits::VarIndex& index)
937        : _host_it(host_it), _index(index) {}
938
939      reference operator*() {
940        return std::make_pair(_index(_host_it->first), _host_it->second);
941      }
942
943      pointer operator->() {
944        return pointer(operator*());
945      }
946
947      ExprIterator& operator++() { ++_host_it; return *this; }
948      ExprIterator operator++(int) {
949        ExprIterator tmp(*this); ++_host_it; return tmp;
950      }
951
952      ExprIterator& operator--() { --_host_it; return *this; }
953      ExprIterator operator--(int) {
954        ExprIterator tmp(*this); --_host_it; return tmp;
955      }
956
957      bool operator==(const ExprIterator& it) const {
958        return _host_it == it._host_it;
959      }
960
961      bool operator!=(const ExprIterator& it) const {
962        return _host_it != it._host_it;
963      }
964
965    };
966
967  protected:
968
969    //Abstract virtual functions
970
971    virtual int _addColId(int col) { return _cols.addIndex(col); }
972    virtual int _addRowId(int row) { return _rows.addIndex(row); }
973
974    virtual void _eraseColId(int col) { _cols.eraseIndex(col); }
975    virtual void _eraseRowId(int row) { _rows.eraseIndex(row); }
976
977    virtual int _addCol() = 0;
978    virtual int _addRow() = 0;
979
980    virtual int _addRow(Value l, ExprIterator b, ExprIterator e, Value u) {
981      int row = _addRow();
982      _setRowCoeffs(row, b, e);
983      _setRowLowerBound(row, l);
984      _setRowUpperBound(row, u);
985      return row;
986    }
987
988    virtual void _eraseCol(int col) = 0;
989    virtual void _eraseRow(int row) = 0;
990
991    virtual void _getColName(int col, std::string& name) const = 0;
992    virtual void _setColName(int col, const std::string& name) = 0;
993    virtual int _colByName(const std::string& name) const = 0;
994
995    virtual void _getRowName(int row, std::string& name) const = 0;
996    virtual void _setRowName(int row, const std::string& name) = 0;
997    virtual int _rowByName(const std::string& name) const = 0;
998
999    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1000    virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
1001
1002    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
1003    virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
1004
1005    virtual void _setCoeff(int row, int col, Value value) = 0;
1006    virtual Value _getCoeff(int row, int col) const = 0;
1007
1008    virtual void _setColLowerBound(int i, Value value) = 0;
1009    virtual Value _getColLowerBound(int i) const = 0;
1010
1011    virtual void _setColUpperBound(int i, Value value) = 0;
1012    virtual Value _getColUpperBound(int i) const = 0;
1013
1014    virtual void _setRowLowerBound(int i, Value value) = 0;
1015    virtual Value _getRowLowerBound(int i) const = 0;
1016
1017    virtual void _setRowUpperBound(int i, Value value) = 0;
1018    virtual Value _getRowUpperBound(int i) const = 0;
1019
1020    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
1021    virtual void _getObjCoeffs(InsertIterator b) const = 0;
1022
1023    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
1024    virtual Value _getObjCoeff(int i) const = 0;
1025
1026    virtual void _setSense(Sense) = 0;
1027    virtual Sense _getSense() const = 0;
1028
1029    virtual void _clear() = 0;
1030
1031    virtual const char* _solverName() const = 0;
1032
1033    virtual void _messageLevel(MessageLevel level) = 0;
1034
1035    //Own protected stuff
1036
1037    //Constant component of the objective function
1038    Value obj_const_comp;
1039
1040    LpBase() : _rows(), _cols(), obj_const_comp(0) {}
1041
1042  public:
1043
1044    ///Unsupported file format exception
1045    class UnsupportedFormatError : public Exception
1046    {
1047      std::string _format;
1048      mutable std::string _what;
1049    public:
1050      explicit UnsupportedFormatError(std::string format) throw()
1051        : _format(format) { }
1052      virtual ~UnsupportedFormatError() throw() {}
1053      virtual const char* what() const throw() {
1054        try {
1055          _what.clear();
1056          std::ostringstream oss;
1057          oss << "lemon::UnsupportedFormatError: " << _format;
1058          _what = oss.str();
1059        }
1060        catch (...) {}
1061        if (!_what.empty()) return _what.c_str();
1062        else return "lemon::UnsupportedFormatError";
1063      }
1064    };
1065
1066  protected:
1067    virtual void _write(std::string, std::string format) const
1068    {
1069      throw UnsupportedFormatError(format);
1070    }
1071
1072  public:
1073
1074    /// Virtual destructor
1075    virtual ~LpBase() {}
1076
1077    ///Gives back the name of the solver.
1078    const char* solverName() const {return _solverName();}
1079
1080    ///\name Build Up and Modify the LP
1081
1082    ///@{
1083
1084    ///Add a new empty column (i.e a new variable) to the LP
1085    Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
1086
1087    ///\brief Adds several new columns (i.e variables) at once
1088    ///
1089    ///This magic function takes a container as its argument and fills
1090    ///its elements with new columns (i.e. variables)
1091    ///\param t can be
1092    ///- a standard STL compatible iterable container with
1093    ///\ref Col as its \c values_type like
1094    ///\code
1095    ///std::vector<LpBase::Col>
1096    ///std::list<LpBase::Col>
1097    ///\endcode
1098    ///- a standard STL compatible iterable container with
1099    ///\ref Col as its \c mapped_type like
1100    ///\code
1101    ///std::map<AnyType,LpBase::Col>
1102    ///\endcode
1103    ///- an iterable lemon \ref concepts::WriteMap "write map" like
1104    ///\code
1105    ///ListGraph::NodeMap<LpBase::Col>
1106    ///ListGraph::ArcMap<LpBase::Col>
1107    ///\endcode
1108    ///\return The number of the created column.
1109#ifdef DOXYGEN
1110    template<class T>
1111    int addColSet(T &t) { return 0;}
1112#else
1113    template<class T>
1114    typename enable_if<typename T::value_type::LpCol,int>::type
1115    addColSet(T &t,dummy<0> = 0) {
1116      int s=0;
1117      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
1118      return s;
1119    }
1120    template<class T>
1121    typename enable_if<typename T::value_type::second_type::LpCol,
1122                       int>::type
1123    addColSet(T &t,dummy<1> = 1) {
1124      int s=0;
1125      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1126        i->second=addCol();
1127        s++;
1128      }
1129      return s;
1130    }
1131    template<class T>
1132    typename enable_if<typename T::MapIt::Value::LpCol,
1133                       int>::type
1134    addColSet(T &t,dummy<2> = 2) {
1135      int s=0;
1136      for(typename T::MapIt i(t); i!=INVALID; ++i)
1137        {
1138          i.set(addCol());
1139          s++;
1140        }
1141      return s;
1142    }
1143#endif
1144
1145    ///Set a column (i.e a dual constraint) of the LP
1146
1147    ///\param c is the column to be modified
1148    ///\param e is a dual linear expression (see \ref DualExpr)
1149    ///a better one.
1150    void col(Col c, const DualExpr &e) {
1151      e.simplify();
1152      _setColCoeffs(_cols(id(c)), ExprIterator(e.comps.begin(), _rows),
1153                    ExprIterator(e.comps.end(), _rows));
1154    }
1155
1156    ///Get a column (i.e a dual constraint) of the LP
1157
1158    ///\param c is the column to get
1159    ///\return the dual expression associated to the column
1160    DualExpr col(Col c) const {
1161      DualExpr e;
1162      _getColCoeffs(_cols(id(c)), InsertIterator(e.comps, _rows));
1163      return e;
1164    }
1165
1166    ///Add a new column to the LP
1167
1168    ///\param e is a dual linear expression (see \ref DualExpr)
1169    ///\param o is the corresponding component of the objective
1170    ///function. It is 0 by default.
1171    ///\return The created column.
1172    Col addCol(const DualExpr &e, Value o = 0) {
1173      Col c=addCol();
1174      col(c,e);
1175      objCoeff(c,o);
1176      return c;
1177    }
1178
1179    ///Add a new empty row (i.e a new constraint) to the LP
1180
1181    ///This function adds a new empty row (i.e a new constraint) to the LP.
1182    ///\return The created row
1183    Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
1184
1185    ///\brief Add several new rows (i.e constraints) at once
1186    ///
1187    ///This magic function takes a container as its argument and fills
1188    ///its elements with new row (i.e. variables)
1189    ///\param t can be
1190    ///- a standard STL compatible iterable container with
1191    ///\ref LpBase::Row "Row" as its \c values_type like
1192    ///\code
1193    ///std::vector<LpBase::Row>
1194    ///std::list<LpBase::Row>
1195    ///\endcode
1196    ///- a standard STL compatible iterable container with
1197    ///\ref LpBase::Row "Row" as its \c mapped_type like
1198    ///\code
1199    ///std::map<AnyType,LpBase::Row>
1200    ///\endcode
1201    ///- an iterable lemon \ref concepts::WriteMap "write map" like
1202    ///\code
1203    ///ListGraph::NodeMap<LpBase::Row>
1204    ///ListGraph::ArcMap<LpBase::Row>
1205    ///\endcode
1206    ///\return The number of rows created.
1207#ifdef DOXYGEN
1208    template<class T>
1209    int addRowSet(T &t) { return 0;}
1210#else
1211    template<class T>
1212    typename enable_if<typename T::value_type::LpRow,int>::type
1213    addRowSet(T &t, dummy<0> = 0) {
1214      int s=0;
1215      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1216      return s;
1217    }
1218    template<class T>
1219    typename enable_if<typename T::value_type::second_type::LpRow, int>::type
1220    addRowSet(T &t, dummy<1> = 1) {
1221      int s=0;
1222      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1223        i->second=addRow();
1224        s++;
1225      }
1226      return s;
1227    }
1228    template<class T>
1229    typename enable_if<typename T::MapIt::Value::LpRow, int>::type
1230    addRowSet(T &t, dummy<2> = 2) {
1231      int s=0;
1232      for(typename T::MapIt i(t); i!=INVALID; ++i)
1233        {
1234          i.set(addRow());
1235          s++;
1236        }
1237      return s;
1238    }
1239#endif
1240
1241    ///Set a row (i.e a constraint) of the LP
1242
1243    ///\param r is the row to be modified
1244    ///\param l is lower bound (-\ref INF means no bound)
1245    ///\param e is a linear expression (see \ref Expr)
1246    ///\param u is the upper bound (\ref INF means no bound)
1247    void row(Row r, Value l, const Expr &e, Value u) {
1248      e.simplify();
1249      _setRowCoeffs(_rows(id(r)), ExprIterator(e.comps.begin(), _cols),
1250                    ExprIterator(e.comps.end(), _cols));
1251      _setRowLowerBound(_rows(id(r)),l - *e);
1252      _setRowUpperBound(_rows(id(r)),u - *e);
1253    }
1254
1255    ///Set a row (i.e a constraint) of the LP
1256
1257    ///\param r is the row to be modified
1258    ///\param c is a linear expression (see \ref Constr)
1259    void row(Row r, const Constr &c) {
1260      row(r, c.lowerBounded()?c.lowerBound():-INF,
1261          c.expr(), c.upperBounded()?c.upperBound():INF);
1262    }
1263
1264
1265    ///Get a row (i.e a constraint) of the LP
1266
1267    ///\param r is the row to get
1268    ///\return the expression associated to the row
1269    Expr row(Row r) const {
1270      Expr e;
1271      _getRowCoeffs(_rows(id(r)), InsertIterator(e.comps, _cols));
1272      return e;
1273    }
1274
1275    ///Add a new row (i.e a new constraint) to the LP
1276
1277    ///\param l is the lower bound (-\ref INF means no bound)
1278    ///\param e is a linear expression (see \ref Expr)
1279    ///\param u is the upper bound (\ref INF means no bound)
1280    ///\return The created row.
1281    Row addRow(Value l,const Expr &e, Value u) {
1282      Row r;
1283      e.simplify();
1284      r._id = _addRowId(_addRow(l - *e, ExprIterator(e.comps.begin(), _cols),
1285                                ExprIterator(e.comps.end(), _cols), u - *e));
1286      return r;
1287    }
1288
1289    ///Add a new row (i.e a new constraint) to the LP
1290
1291    ///\param c is a linear expression (see \ref Constr)
1292    ///\return The created row.
1293    Row addRow(const Constr &c) {
1294      Row r;
1295      c.expr().simplify();
1296      r._id = _addRowId(_addRow(c.lowerBounded()?c.lowerBound()-*c.expr():-INF,
1297                                ExprIterator(c.expr().comps.begin(), _cols),
1298                                ExprIterator(c.expr().comps.end(), _cols),
1299                                c.upperBounded()?c.upperBound()-*c.expr():INF));
1300      return r;
1301    }
1302    ///Erase a column (i.e a variable) from the LP
1303
1304    ///\param c is the column to be deleted
1305    void erase(Col c) {
1306      _eraseCol(_cols(id(c)));
1307      _eraseColId(_cols(id(c)));
1308    }
1309    ///Erase a row (i.e a constraint) from the LP
1310
1311    ///\param r is the row to be deleted
1312    void erase(Row r) {
1313      _eraseRow(_rows(id(r)));
1314      _eraseRowId(_rows(id(r)));
1315    }
1316
1317    /// Get the name of a column
1318
1319    ///\param c is the coresponding column
1320    ///\return The name of the colunm
1321    std::string colName(Col c) const {
1322      std::string name;
1323      _getColName(_cols(id(c)), name);
1324      return name;
1325    }
1326
1327    /// Set the name of a column
1328
1329    ///\param c is the coresponding column
1330    ///\param name The name to be given
1331    void colName(Col c, const std::string& name) {
1332      _setColName(_cols(id(c)), name);
1333    }
1334
1335    /// Get the column by its name
1336
1337    ///\param name The name of the column
1338    ///\return the proper column or \c INVALID
1339    Col colByName(const std::string& name) const {
1340      int k = _colByName(name);
1341      return k != -1 ? Col(_cols[k]) : Col(INVALID);
1342    }
1343
1344    /// Get the name of a row
1345
1346    ///\param r is the coresponding row
1347    ///\return The name of the row
1348    std::string rowName(Row r) const {
1349      std::string name;
1350      _getRowName(_rows(id(r)), name);
1351      return name;
1352    }
1353
1354    /// Set the name of a row
1355
1356    ///\param r is the coresponding row
1357    ///\param name The name to be given
1358    void rowName(Row r, const std::string& name) {
1359      _setRowName(_rows(id(r)), name);
1360    }
1361
1362    /// Get the row by its name
1363
1364    ///\param name The name of the row
1365    ///\return the proper row or \c INVALID
1366    Row rowByName(const std::string& name) const {
1367      int k = _rowByName(name);
1368      return k != -1 ? Row(_rows[k]) : Row(INVALID);
1369    }
1370
1371    /// Set an element of the coefficient matrix of the LP
1372
1373    ///\param r is the row of the element to be modified
1374    ///\param c is the column of the element to be modified
1375    ///\param val is the new value of the coefficient
1376    void coeff(Row r, Col c, Value val) {
1377      _setCoeff(_rows(id(r)),_cols(id(c)), val);
1378    }
1379
1380    /// Get an element of the coefficient matrix of the LP
1381
1382    ///\param r is the row of the element
1383    ///\param c is the column of the element
1384    ///\return the corresponding coefficient
1385    Value coeff(Row r, Col c) const {
1386      return _getCoeff(_rows(id(r)),_cols(id(c)));
1387    }
1388
1389    /// Set the lower bound of a column (i.e a variable)
1390
1391    /// The lower bound of a variable (column) has to be given by an
1392    /// extended number of type Value, i.e. a finite number of type
1393    /// Value or -\ref INF.
1394    void colLowerBound(Col c, Value value) {
1395      _setColLowerBound(_cols(id(c)),value);
1396    }
1397
1398    /// Get the lower bound of a column (i.e a variable)
1399
1400    /// This function returns the lower bound for column (variable) \c c
1401    /// (this might be -\ref INF as well).
1402    ///\return The lower bound for column \c c
1403    Value colLowerBound(Col c) const {
1404      return _getColLowerBound(_cols(id(c)));
1405    }
1406
1407    ///\brief Set the lower bound of  several columns
1408    ///(i.e variables) at once
1409    ///
1410    ///This magic function takes a container as its argument
1411    ///and applies the function on all of its elements.
1412    ///The lower bound of a variable (column) has to be given by an
1413    ///extended number of type Value, i.e. a finite number of type
1414    ///Value or -\ref INF.
1415#ifdef DOXYGEN
1416    template<class T>
1417    void colLowerBound(T &t, Value value) { return 0;}
1418#else
1419    template<class T>
1420    typename enable_if<typename T::value_type::LpCol,void>::type
1421    colLowerBound(T &t, Value value,dummy<0> = 0) {
1422      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1423        colLowerBound(*i, value);
1424      }
1425    }
1426    template<class T>
1427    typename enable_if<typename T::value_type::second_type::LpCol,
1428                       void>::type
1429    colLowerBound(T &t, Value value,dummy<1> = 1) {
1430      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1431        colLowerBound(i->second, value);
1432      }
1433    }
1434    template<class T>
1435    typename enable_if<typename T::MapIt::Value::LpCol,
1436                       void>::type
1437    colLowerBound(T &t, Value value,dummy<2> = 2) {
1438      for(typename T::MapIt i(t); i!=INVALID; ++i){
1439        colLowerBound(*i, value);
1440      }
1441    }
1442#endif
1443
1444    /// Set the upper bound of a column (i.e a variable)
1445
1446    /// The upper bound of a variable (column) has to be given by an
1447    /// extended number of type Value, i.e. a finite number of type
1448    /// Value or \ref INF.
1449    void colUpperBound(Col c, Value value) {
1450      _setColUpperBound(_cols(id(c)),value);
1451    };
1452
1453    /// Get the upper bound of a column (i.e a variable)
1454
1455    /// This function returns the upper bound for column (variable) \c c
1456    /// (this might be \ref INF as well).
1457    /// \return The upper bound for column \c c
1458    Value colUpperBound(Col c) const {
1459      return _getColUpperBound(_cols(id(c)));
1460    }
1461
1462    ///\brief Set the upper bound of  several columns
1463    ///(i.e variables) at once
1464    ///
1465    ///This magic function takes a container as its argument
1466    ///and applies the function on all of its elements.
1467    ///The upper bound of a variable (column) has to be given by an
1468    ///extended number of type Value, i.e. a finite number of type
1469    ///Value or \ref INF.
1470#ifdef DOXYGEN
1471    template<class T>
1472    void colUpperBound(T &t, Value value) { return 0;}
1473#else
1474    template<class T1>
1475    typename enable_if<typename T1::value_type::LpCol,void>::type
1476    colUpperBound(T1 &t, Value value,dummy<0> = 0) {
1477      for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1478        colUpperBound(*i, value);
1479      }
1480    }
1481    template<class T1>
1482    typename enable_if<typename T1::value_type::second_type::LpCol,
1483                       void>::type
1484    colUpperBound(T1 &t, Value value,dummy<1> = 1) {
1485      for(typename T1::iterator i=t.begin();i!=t.end();++i) {
1486        colUpperBound(i->second, value);
1487      }
1488    }
1489    template<class T1>
1490    typename enable_if<typename T1::MapIt::Value::LpCol,
1491                       void>::type
1492    colUpperBound(T1 &t, Value value,dummy<2> = 2) {
1493      for(typename T1::MapIt i(t); i!=INVALID; ++i){
1494        colUpperBound(*i, value);
1495      }
1496    }
1497#endif
1498
1499    /// Set the lower and the upper bounds of a column (i.e a variable)
1500
1501    /// The lower and the upper bounds of
1502    /// a variable (column) have to be given by an
1503    /// extended number of type Value, i.e. a finite number of type
1504    /// Value, -\ref INF or \ref INF.
1505    void colBounds(Col c, Value lower, Value upper) {
1506      _setColLowerBound(_cols(id(c)),lower);
1507      _setColUpperBound(_cols(id(c)),upper);
1508    }
1509
1510    ///\brief Set the lower and the upper bound of several columns
1511    ///(i.e variables) at once
1512    ///
1513    ///This magic function takes a container as its argument
1514    ///and applies the function on all of its elements.
1515    /// The lower and the upper bounds of
1516    /// a variable (column) have to be given by an
1517    /// extended number of type Value, i.e. a finite number of type
1518    /// Value, -\ref INF or \ref INF.
1519#ifdef DOXYGEN
1520    template<class T>
1521    void colBounds(T &t, Value lower, Value upper) { return 0;}
1522#else
1523    template<class T2>
1524    typename enable_if<typename T2::value_type::LpCol,void>::type
1525    colBounds(T2 &t, Value lower, Value upper,dummy<0> = 0) {
1526      for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1527        colBounds(*i, lower, upper);
1528      }
1529    }
1530    template<class T2>
1531    typename enable_if<typename T2::value_type::second_type::LpCol, void>::type
1532    colBounds(T2 &t, Value lower, Value upper,dummy<1> = 1) {
1533      for(typename T2::iterator i=t.begin();i!=t.end();++i) {
1534        colBounds(i->second, lower, upper);
1535      }
1536    }
1537    template<class T2>
1538    typename enable_if<typename T2::MapIt::Value::LpCol, void>::type
1539    colBounds(T2 &t, Value lower, Value upper,dummy<2> = 2) {
1540      for(typename T2::MapIt i(t); i!=INVALID; ++i){
1541        colBounds(*i, lower, upper);
1542      }
1543    }
1544#endif
1545
1546    /// Set the lower bound of a row (i.e a constraint)
1547
1548    /// The lower bound of a constraint (row) has to be given by an
1549    /// extended number of type Value, i.e. a finite number of type
1550    /// Value or -\ref INF.
1551    void rowLowerBound(Row r, Value value) {
1552      _setRowLowerBound(_rows(id(r)),value);
1553    }
1554
1555    /// Get the lower bound of a row (i.e a constraint)
1556
1557    /// This function returns the lower bound for row (constraint) \c c
1558    /// (this might be -\ref INF as well).
1559    ///\return The lower bound for row \c r
1560    Value rowLowerBound(Row r) const {
1561      return _getRowLowerBound(_rows(id(r)));
1562    }
1563
1564    /// Set the upper bound of a row (i.e a constraint)
1565
1566    /// The upper bound of a constraint (row) has to be given by an
1567    /// extended number of type Value, i.e. a finite number of type
1568    /// Value or -\ref INF.
1569    void rowUpperBound(Row r, Value value) {
1570      _setRowUpperBound(_rows(id(r)),value);
1571    }
1572
1573    /// Get the upper bound of a row (i.e a constraint)
1574
1575    /// This function returns the upper bound for row (constraint) \c c
1576    /// (this might be -\ref INF as well).
1577    ///\return The upper bound for row \c r
1578    Value rowUpperBound(Row r) const {
1579      return _getRowUpperBound(_rows(id(r)));
1580    }
1581
1582    ///Set an element of the objective function
1583    void objCoeff(Col c, Value v) {_setObjCoeff(_cols(id(c)),v); };
1584
1585    ///Get an element of the objective function
1586    Value objCoeff(Col c) const { return _getObjCoeff(_cols(id(c))); };
1587
1588    ///Set the objective function
1589
1590    ///\param e is a linear expression of type \ref Expr.
1591    ///
1592    void obj(const Expr& e) {
1593      _setObjCoeffs(ExprIterator(e.comps.begin(), _cols),
1594                    ExprIterator(e.comps.end(), _cols));
1595      obj_const_comp = *e;
1596    }
1597
1598    ///Get the objective function
1599
1600    ///\return the objective function as a linear expression of type
1601    ///Expr.
1602    Expr obj() const {
1603      Expr e;
1604      _getObjCoeffs(InsertIterator(e.comps, _cols));
1605      *e = obj_const_comp;
1606      return e;
1607    }
1608
1609
1610    ///Set the direction of optimization
1611    void sense(Sense sense) { _setSense(sense); }
1612
1613    ///Query the direction of the optimization
1614    Sense sense() const {return _getSense(); }
1615
1616    ///Set the sense to maximization
1617    void max() { _setSense(MAX); }
1618
1619    ///Set the sense to maximization
1620    void min() { _setSense(MIN); }
1621
1622    ///Clear the problem
1623    void clear() { _clear(); _rows.clear(); _cols.clear(); }
1624
1625    /// Set the message level of the solver
1626    void messageLevel(MessageLevel level) { _messageLevel(level); }
1627
1628    /// Write the problem to a file in the given format
1629
1630    /// This function writes the problem to a file in the given format.
1631    /// Different solver backends may support different formats.
1632    /// Trying to write in an unsupported format will trigger
1633    /// \ref UnsupportedFormatError. For the supported formats,
1634    /// visit the documentation of the base class of the related backends
1635    /// (\ref CplexBase, \ref GlpkBase etc.)
1636    /// \param file The file path
1637    /// \param format The output file format.
1638    void write(std::string file, std::string format = "MPS") const
1639    {
1640      _write(file.c_str(),format.c_str());
1641    }
1642
1643    ///@}
1644
1645  };
1646
1647  /// Addition
1648
1649  ///\relates LpBase::Expr
1650  ///
1651  inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
1652    LpBase::Expr tmp(a);
1653    tmp+=b;
1654    return tmp;
1655  }
1656  ///Substraction
1657
1658  ///\relates LpBase::Expr
1659  ///
1660  inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
1661    LpBase::Expr tmp(a);
1662    tmp-=b;
1663    return tmp;
1664  }
1665  ///Multiply with constant
1666
1667  ///\relates LpBase::Expr
1668  ///
1669  inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
1670    LpBase::Expr tmp(a);
1671    tmp*=b;
1672    return tmp;
1673  }
1674
1675  ///Multiply with constant
1676
1677  ///\relates LpBase::Expr
1678  ///
1679  inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
1680    LpBase::Expr tmp(b);
1681    tmp*=a;
1682    return tmp;
1683  }
1684  ///Divide with constant
1685
1686  ///\relates LpBase::Expr
1687  ///
1688  inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
1689    LpBase::Expr tmp(a);
1690    tmp/=b;
1691    return tmp;
1692  }
1693
1694  ///Create constraint
1695
1696  ///\relates LpBase::Constr
1697  ///
1698  inline LpBase::Constr operator<=(const LpBase::Expr &e,
1699                                   const LpBase::Expr &f) {
1700    return LpBase::Constr(0, f - e, LpBase::NaN);
1701  }
1702
1703  ///Create constraint
1704
1705  ///\relates LpBase::Constr
1706  ///
1707  inline LpBase::Constr operator<=(const LpBase::Value &e,
1708                                   const LpBase::Expr &f) {
1709    return LpBase::Constr(e, f, LpBase::NaN);
1710  }
1711
1712  ///Create constraint
1713
1714  ///\relates LpBase::Constr
1715  ///
1716  inline LpBase::Constr operator<=(const LpBase::Expr &e,
1717                                   const LpBase::Value &f) {
1718    return LpBase::Constr(LpBase::NaN, e, f);
1719  }
1720
1721  ///Create constraint
1722
1723  ///\relates LpBase::Constr
1724  ///
1725  inline LpBase::Constr operator>=(const LpBase::Expr &e,
1726                                   const LpBase::Expr &f) {
1727    return LpBase::Constr(0, e - f, LpBase::NaN);
1728  }
1729
1730
1731  ///Create constraint
1732
1733  ///\relates LpBase::Constr
1734  ///
1735  inline LpBase::Constr operator>=(const LpBase::Value &e,
1736                                   const LpBase::Expr &f) {
1737    return LpBase::Constr(LpBase::NaN, f, e);
1738  }
1739
1740
1741  ///Create constraint
1742
1743  ///\relates LpBase::Constr
1744  ///
1745  inline LpBase::Constr operator>=(const LpBase::Expr &e,
1746                                   const LpBase::Value &f) {
1747    return LpBase::Constr(f, e, LpBase::NaN);
1748  }
1749
1750  ///Create constraint
1751
1752  ///\relates LpBase::Constr
1753  ///
1754  inline LpBase::Constr operator==(const LpBase::Expr &e,
1755                                   const LpBase::Value &f) {
1756    return LpBase::Constr(f, e, f);
1757  }
1758
1759  ///Create constraint
1760
1761  ///\relates LpBase::Constr
1762  ///
1763  inline LpBase::Constr operator==(const LpBase::Expr &e,
1764                                   const LpBase::Expr &f) {
1765    return LpBase::Constr(0, f - e, 0);
1766  }
1767
1768  ///Create constraint
1769
1770  ///\relates LpBase::Constr
1771  ///
1772  inline LpBase::Constr operator<=(const LpBase::Value &n,
1773                                   const LpBase::Constr &c) {
1774    LpBase::Constr tmp(c);
1775    LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1776    tmp.lowerBound()=n;
1777    return tmp;
1778  }
1779  ///Create constraint
1780
1781  ///\relates LpBase::Constr
1782  ///
1783  inline LpBase::Constr operator<=(const LpBase::Constr &c,
1784                                   const LpBase::Value &n)
1785  {
1786    LpBase::Constr tmp(c);
1787    LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1788    tmp.upperBound()=n;
1789    return tmp;
1790  }
1791
1792  ///Create constraint
1793
1794  ///\relates LpBase::Constr
1795  ///
1796  inline LpBase::Constr operator>=(const LpBase::Value &n,
1797                                   const LpBase::Constr &c) {
1798    LpBase::Constr tmp(c);
1799    LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1800    tmp.upperBound()=n;
1801    return tmp;
1802  }
1803  ///Create constraint
1804
1805  ///\relates LpBase::Constr
1806  ///
1807  inline LpBase::Constr operator>=(const LpBase::Constr &c,
1808                                   const LpBase::Value &n)
1809  {
1810    LpBase::Constr tmp(c);
1811    LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1812    tmp.lowerBound()=n;
1813    return tmp;
1814  }
1815
1816  ///Addition
1817
1818  ///\relates LpBase::DualExpr
1819  ///
1820  inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
1821                                    const LpBase::DualExpr &b) {
1822    LpBase::DualExpr tmp(a);
1823    tmp+=b;
1824    return tmp;
1825  }
1826  ///Substraction
1827
1828  ///\relates LpBase::DualExpr
1829  ///
1830  inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
1831                                    const LpBase::DualExpr &b) {
1832    LpBase::DualExpr tmp(a);
1833    tmp-=b;
1834    return tmp;
1835  }
1836  ///Multiply with constant
1837
1838  ///\relates LpBase::DualExpr
1839  ///
1840  inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
1841                                    const LpBase::Value &b) {
1842    LpBase::DualExpr tmp(a);
1843    tmp*=b;
1844    return tmp;
1845  }
1846
1847  ///Multiply with constant
1848
1849  ///\relates LpBase::DualExpr
1850  ///
1851  inline LpBase::DualExpr operator*(const LpBase::Value &a,
1852                                    const LpBase::DualExpr &b) {
1853    LpBase::DualExpr tmp(b);
1854    tmp*=a;
1855    return tmp;
1856  }
1857  ///Divide with constant
1858
1859  ///\relates LpBase::DualExpr
1860  ///
1861  inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
1862                                    const LpBase::Value &b) {
1863    LpBase::DualExpr tmp(a);
1864    tmp/=b;
1865    return tmp;
1866  }
1867
1868  /// \ingroup lp_group
1869  ///
1870  /// \brief Common base class for LP solvers
1871  ///
1872  /// This class is an abstract base class for LP solvers. This class
1873  /// provides a full interface for set and modify an LP problem,
1874  /// solve it and retrieve the solution. You can use one of the
1875  /// descendants as a concrete implementation, or the \c Lp
1876  /// default LP solver. However, if you would like to handle LP
1877  /// solvers as reference or pointer in a generic way, you can use
1878  /// this class directly.
1879  class LpSolver : virtual public LpBase {
1880  public:
1881
1882    /// The problem types for primal and dual problems
1883    enum ProblemType {
1884      /// = 0. Feasible solution hasn't been found (but may exist).
1885      UNDEFINED = 0,
1886      /// = 1. The problem has no feasible solution.
1887      INFEASIBLE = 1,
1888      /// = 2. Feasible solution found.
1889      FEASIBLE = 2,
1890      /// = 3. Optimal solution exists and found.
1891      OPTIMAL = 3,
1892      /// = 4. The cost function is unbounded.
1893      UNBOUNDED = 4
1894    };
1895
1896    ///The basis status of variables
1897    enum VarStatus {
1898      /// The variable is in the basis
1899      BASIC,
1900      /// The variable is free, but not basic
1901      FREE,
1902      /// The variable has active lower bound
1903      LOWER,
1904      /// The variable has active upper bound
1905      UPPER,
1906      /// The variable is non-basic and fixed
1907      FIXED
1908    };
1909
1910  protected:
1911
1912    virtual SolveExitStatus _solve() = 0;
1913
1914    virtual Value _getPrimal(int i) const = 0;
1915    virtual Value _getDual(int i) const = 0;
1916
1917    virtual Value _getPrimalRay(int i) const = 0;
1918    virtual Value _getDualRay(int i) const = 0;
1919
1920    virtual Value _getPrimalValue() const = 0;
1921
1922    virtual VarStatus _getColStatus(int i) const = 0;
1923    virtual VarStatus _getRowStatus(int i) const = 0;
1924
1925    virtual ProblemType _getPrimalType() const = 0;
1926    virtual ProblemType _getDualType() const = 0;
1927
1928  public:
1929
1930    ///Allocate a new LP problem instance
1931    virtual LpSolver* newSolver() const = 0;
1932    ///Make a copy of the LP problem
1933    virtual LpSolver* cloneSolver() const = 0;
1934
1935    ///\name Solve the LP
1936
1937    ///@{
1938
1939    ///\e Solve the LP problem at hand
1940    ///
1941    ///\return The result of the optimization procedure. Possible
1942    ///values and their meanings can be found in the documentation of
1943    ///\ref SolveExitStatus.
1944    SolveExitStatus solve() { return _solve(); }
1945
1946    ///@}
1947
1948    ///\name Obtain the Solution
1949
1950    ///@{
1951
1952    /// The type of the primal problem
1953    ProblemType primalType() const {
1954      return _getPrimalType();
1955    }
1956
1957    /// The type of the dual problem
1958    ProblemType dualType() const {
1959      return _getDualType();
1960    }
1961
1962    /// Return the primal value of the column
1963
1964    /// Return the primal value of the column.
1965    /// \pre The problem is solved.
1966    Value primal(Col c) const { return _getPrimal(_cols(id(c))); }
1967
1968    /// Return the primal value of the expression
1969
1970    /// Return the primal value of the expression, i.e. the dot
1971    /// product of the primal solution and the expression.
1972    /// \pre The problem is solved.
1973    Value primal(const Expr& e) const {
1974      double res = *e;
1975      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1976        res += *c * primal(c);
1977      }
1978      return res;
1979    }
1980    /// Returns a component of the primal ray
1981
1982    /// The primal ray is solution of the modified primal problem,
1983    /// where we change each finite bound to 0, and we looking for a
1984    /// negative objective value in case of minimization, and positive
1985    /// objective value for maximization. If there is such solution,
1986    /// that proofs the unsolvability of the dual problem, and if a
1987    /// feasible primal solution exists, then the unboundness of
1988    /// primal problem.
1989    ///
1990    /// \pre The problem is solved and the dual problem is infeasible.
1991    /// \note Some solvers does not provide primal ray calculation
1992    /// functions.
1993    Value primalRay(Col c) const { return _getPrimalRay(_cols(id(c))); }
1994
1995    /// Return the dual value of the row
1996
1997    /// Return the dual value of the row.
1998    /// \pre The problem is solved.
1999    Value dual(Row r) const { return _getDual(_rows(id(r))); }
2000
2001    /// Return the dual value of the dual expression
2002
2003    /// Return the dual value of the dual expression, i.e. the dot
2004    /// product of the dual solution and the dual expression.
2005    /// \pre The problem is solved.
2006    Value dual(const DualExpr& e) const {
2007      double res = 0.0;
2008      for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
2009        res += *r * dual(r);
2010      }
2011      return res;
2012    }
2013
2014    /// Returns a component of the dual ray
2015
2016    /// The dual ray is solution of the modified primal problem, where
2017    /// we change each finite bound to 0 (i.e. the objective function
2018    /// coefficients in the primal problem), and we looking for a
2019    /// ositive objective value. If there is such solution, that
2020    /// proofs the unsolvability of the primal problem, and if a
2021    /// feasible dual solution exists, then the unboundness of
2022    /// dual problem.
2023    ///
2024    /// \pre The problem is solved and the primal problem is infeasible.
2025    /// \note Some solvers does not provide dual ray calculation
2026    /// functions.
2027    Value dualRay(Row r) const { return _getDualRay(_rows(id(r))); }
2028
2029    /// Return the basis status of the column
2030
2031    /// \see VarStatus
2032    VarStatus colStatus(Col c) const { return _getColStatus(_cols(id(c))); }
2033
2034    /// Return the basis status of the row
2035
2036    /// \see VarStatus
2037    VarStatus rowStatus(Row r) const { return _getRowStatus(_rows(id(r))); }
2038
2039    ///The value of the objective function
2040
2041    ///\return
2042    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
2043    /// of the primal problem, depending on whether we minimize or maximize.
2044    ///- \ref NaN if no primal solution is found.
2045    ///- The (finite) objective value if an optimal solution is found.
2046    Value primal() const { return _getPrimalValue()+obj_const_comp;}
2047    ///@}
2048
2049  protected:
2050
2051  };
2052
2053
2054  /// \ingroup lp_group
2055  ///
2056  /// \brief Common base class for MIP solvers
2057  ///
2058  /// This class is an abstract base class for MIP solvers. This class
2059  /// provides a full interface for set and modify an MIP problem,
2060  /// solve it and retrieve the solution. You can use one of the
2061  /// descendants as a concrete implementation, or the \c Lp
2062  /// default MIP solver. However, if you would like to handle MIP
2063  /// solvers as reference or pointer in a generic way, you can use
2064  /// this class directly.
2065  class MipSolver : virtual public LpBase {
2066  public:
2067
2068    /// The problem types for MIP problems
2069    enum ProblemType {
2070      /// = 0. Feasible solution hasn't been found (but may exist).
2071      UNDEFINED = 0,
2072      /// = 1. The problem has no feasible solution.
2073      INFEASIBLE = 1,
2074      /// = 2. Feasible solution found.
2075      FEASIBLE = 2,
2076      /// = 3. Optimal solution exists and found.
2077      OPTIMAL = 3,
2078      /// = 4. The cost function is unbounded.
2079      ///The Mip or at least the relaxed problem is unbounded.
2080      UNBOUNDED = 4
2081    };
2082
2083    ///Allocate a new MIP problem instance
2084    virtual MipSolver* newSolver() const = 0;
2085    ///Make a copy of the MIP problem
2086    virtual MipSolver* cloneSolver() const = 0;
2087
2088    ///\name Solve the MIP
2089
2090    ///@{
2091
2092    /// Solve the MIP problem at hand
2093    ///
2094    ///\return The result of the optimization procedure. Possible
2095    ///values and their meanings can be found in the documentation of
2096    ///\ref SolveExitStatus.
2097    SolveExitStatus solve() { return _solve(); }
2098
2099    ///@}
2100
2101    ///\name Set Column Type
2102    ///@{
2103
2104    ///Possible variable (column) types (e.g. real, integer, binary etc.)
2105    enum ColTypes {
2106      /// = 0. Continuous variable (default).
2107      REAL = 0,
2108      /// = 1. Integer variable.
2109      INTEGER = 1
2110    };
2111
2112    ///Sets the type of the given column to the given type
2113
2114    ///Sets the type of the given column to the given type.
2115    ///
2116    void colType(Col c, ColTypes col_type) {
2117      _setColType(_cols(id(c)),col_type);
2118    }
2119
2120    ///Gives back the type of the column.
2121
2122    ///Gives back the type of the column.
2123    ///
2124    ColTypes colType(Col c) const {
2125      return _getColType(_cols(id(c)));
2126    }
2127    ///@}
2128
2129    ///\name Obtain the Solution
2130
2131    ///@{
2132
2133    /// The type of the MIP problem
2134    ProblemType type() const {
2135      return _getType();
2136    }
2137
2138    /// Return the value of the row in the solution
2139
2140    ///  Return the value of the row in the solution.
2141    /// \pre The problem is solved.
2142    Value sol(Col c) const { return _getSol(_cols(id(c))); }
2143
2144    /// Return the value of the expression in the solution
2145
2146    /// Return the value of the expression in the solution, i.e. the
2147    /// dot product of the solution and the expression.
2148    /// \pre The problem is solved.
2149    Value sol(const Expr& e) const {
2150      double res = *e;
2151      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
2152        res += *c * sol(c);
2153      }
2154      return res;
2155    }
2156    ///The value of the objective function
2157
2158    ///\return
2159    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
2160    /// of the problem, depending on whether we minimize or maximize.
2161    ///- \ref NaN if no primal solution is found.
2162    ///- The (finite) objective value if an optimal solution is found.
2163    Value solValue() const { return _getSolValue()+obj_const_comp;}
2164    ///@}
2165
2166  protected:
2167
2168    virtual SolveExitStatus _solve() = 0;
2169    virtual ColTypes _getColType(int col) const = 0;
2170    virtual void _setColType(int col, ColTypes col_type) = 0;
2171    virtual ProblemType _getType() const = 0;
2172    virtual Value _getSol(int i) const = 0;
2173    virtual Value _getSolValue() const = 0;
2174
2175  };
2176
2177
2178
2179} //namespace lemon
2180
2181#endif //LEMON_LP_BASE_H
Note: See TracBrowser for help on using the repository browser.