COIN-OR::LEMON - Graph Library

source: lemon/lemon/lp_base.h @ 1336:0759d974de81

Last change on this file since 1336:0759d974de81 was 1336:0759d974de81, checked in by Gabor Gevay <ggab90@…>, 6 years ago

STL style iterators (#325)

For

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