COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/lp/lp_base.h @ 1273:2b2ffa625775

Last change on this file since 1273:2b2ffa625775 was 1273:2b2ffa625775, checked in by Alpar Juttner, 19 years ago
  • Better (but still incomplete) doc
  • lp_test runs correctly
File size: 19.9 KB
Line 
1/* -*- C++ -*-
2 * src/lemon/lp_base.h - Part of LEMON, a generic C++ optimization library
3 *
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5 * (Egervary Combinatorial Optimization Research Group, EGRES).
6 *
7 * Permission to use, modify and distribute this software is granted
8 * provided that this copyright notice appears in all copies. For
9 * precise terms see the accompanying LICENSE file.
10 *
11 * This software is provided "AS IS" with no warranty of any kind,
12 * express or implied, and with no claim as to its suitability for any
13 * purpose.
14 *
15 */
16
17#ifndef LEMON_LP_BASE_H
18#define LEMON_LP_BASE_H
19
20#include<vector>
21#include<map>
22#include<limits>
23#include<math.h>
24
25#include<lemon/utility.h>
26#include<lemon/error.h>
27#include<lemon/invalid.h>
28
29//#include"lin_expr.h"
30
31///\file
32///\brief The interface of the LP solver interface.
33namespace lemon {
34 
35  ///Internal data structure to convert floating id's to fix one's
36   
37  ///\todo This might by implemented to be usable in other places.
38  class _FixId
39  {
40    std::vector<int> index;
41    std::vector<int> cross;
42    int first_free;
43  public:
44    _FixId() : first_free(-1) {};
45    ///Convert a floating id to a fix one
46
47    ///\param n is a floating id
48    ///\return the corresponding fix id
49    int fixId(int n) {return cross[n];}
50    ///Convert a fix id to a floating one
51
52    ///\param n is a fix id
53    ///\return the corresponding floating id
54    int floatingId(int n) { return index[n];}
55    ///Add a new floating id.
56
57    ///\param n is a floating id
58    ///\return the fix id of the new value
59    ///\todo Multiple additions should also be handled.
60    int insert(int n)
61    {
62      if(n>=int(cross.size())) {
63        cross.resize(n+1);
64        if(first_free==-1) {
65          cross[n]=index.size();
66          index.push_back(n);
67        }
68        else {
69          cross[n]=first_free;
70          int next=index[first_free];
71          index[first_free]=n;
72          first_free=next;
73        }
74        return cross[n];
75      }
76      ///\todo Create an own exception type.
77      else throw LogicError(); //floatingId-s must form a continuous range;
78    }
79    ///Remove a fix id.
80
81    ///\param n is a fix id
82    ///
83    void erase(int n)
84    {
85      int fl=index[n];
86      index[n]=first_free;
87      first_free=n;
88      for(int i=fl+1;i<int(cross.size());++i) {
89        cross[i-1]=cross[i];
90        index[cross[i]]--;
91      }
92      cross.pop_back();
93    }
94    ///An upper bound on the largest fix id.
95
96    ///\todo Do we need this?
97    ///
98    std::size_t maxFixId() { return cross.size()-1; }
99 
100  };
101   
102  ///Common base class for LP solvers
103  class LpSolverBase {
104   
105  public:
106
107    ///\e
108    enum SolutionType {
109      ///\e
110      INFEASIBLE = 0,
111      ///\e
112      UNBOUNDED = 1,
113      ///\e
114      OPTIMAL = 2,
115      ///\e
116      FEASIBLE = 3,
117    };
118     
119    ///The floating point type used by the solver
120    typedef double Value;
121    ///The infinity constant
122    static const Value INF;
123    ///The not a number constant
124    static const Value NaN;
125   
126    ///Refer to a column of the LP.
127
128    ///This type is used to refer to a column of the LP.
129    ///
130    ///Its value remains valid and correct even after the addition or erase of
131    ///other columns.
132    ///
133    ///\todo Document what can one do with a Col (INVALID, comparing,
134    ///it is similar to Node/Edge)
135    class Col {
136    protected:
137      int id;
138      friend class LpSolverBase;
139    public:
140      typedef Value ExprValue;
141      typedef True LpSolverCol;
142      Col() {}
143      Col(const Invalid&) : id(-1) {}
144      bool operator<(Col c) const  {return id<c.id;}
145      bool operator==(Col c) const  {return id==c.id;}
146      bool operator!=(Col c) const  {return id==c.id;}
147    };
148
149    ///Refer to a row of the LP.
150
151    ///This type is used to refer to a row of the LP.
152    ///
153    ///Its value remains valid and correct even after the addition or erase of
154    ///other rows.
155    ///
156    ///\todo Document what can one do with a Row (INVALID, comparing,
157    ///it is similar to Node/Edge)
158    class Row {
159    protected:
160      int id;
161      friend class LpSolverBase;
162    public:
163      typedef Value ExprValue;
164      typedef True LpSolverRow;
165      Row() {}
166      Row(const Invalid&) : id(-1) {}
167      typedef True LpSolverRow;
168      bool operator<(Row c) const  {return id<c.id;}
169      bool operator==(Row c) const  {return id==c.id;}
170      bool operator!=(Row c) const  {return id==c.id;}
171   };
172   
173    ///Linear expression
174    //    typedef SparseLinExpr<Col> Expr;
175    class Expr : public std::map<Col,Value>
176    {
177    public:
178      typedef LpSolverBase::Col Key;
179      typedef LpSolverBase::Value Value;
180     
181    protected:
182      typedef std::map<Col,Value> Base;
183     
184      Value const_comp;
185  public:
186      typedef True IsLinExpression;
187      ///\e
188      Expr() : Base(), const_comp(0) { }
189      ///\e
190      Expr(const Key &v) : const_comp(0) {
191        Base::insert(std::make_pair(v, 1));
192      }
193      ///\e
194      Expr(const Value &v) : const_comp(v) {}
195      ///\e
196      void set(const Key &v,const Value &c) {
197        Base::insert(std::make_pair(v, c));
198      }
199      ///\e
200      Value &constComp() { return const_comp; }
201      ///\e
202      const Value &constComp() const { return const_comp; }
203     
204      ///Removes the components with zero coefficient.
205      void simplify() {
206        for (Base::iterator i=Base::begin(); i!=Base::end();) {
207          Base::iterator j=i;
208          ++j;
209          if ((*i).second==0) Base::erase(i);
210          j=i;
211        }
212      }
213
214      ///Sets all coefficients and the constant component to 0.
215      void clear() {
216        Base::clear();
217        const_comp=0;
218      }
219
220      ///\e
221      Expr &operator+=(const Expr &e) {
222        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
223          (*this)[j->first]+=j->second;
224        ///\todo it might be speeded up using "hints"
225        const_comp+=e.const_comp;
226        return *this;
227      }
228      ///\e
229      Expr &operator-=(const Expr &e) {
230        for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
231          (*this)[j->first]-=j->second;
232        const_comp-=e.const_comp;
233        return *this;
234      }
235      ///\e
236      Expr &operator*=(const Value &c) {
237        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
238          j->second*=c;
239        const_comp*=c;
240        return *this;
241      }
242      ///\e
243      Expr &operator/=(const Value &c) {
244        for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
245          j->second/=c;
246        const_comp/=c;
247        return *this;
248      }
249    };
250   
251    ///Linear constraint
252    //typedef LinConstr<Expr> Constr;
253    class Constr
254    {
255    public:
256      typedef LpSolverBase::Expr Expr;
257      typedef Expr::Key Key;
258      typedef Expr::Value Value;
259     
260      static const Value INF;
261      static const Value NaN;
262      //     static const Value INF=0;
263      //     static const Value NaN=1;
264     
265    protected:
266      Expr _expr;
267      Value _lb,_ub;
268    public:
269      ///\e
270      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
271      ///\e
272      Constr(Value lb,const Expr &e,Value ub) :
273        _expr(e), _lb(lb), _ub(ub) {}
274      ///\e
275      Constr(const Expr &e,Value ub) :
276        _expr(e), _lb(NaN), _ub(ub) {}
277      ///\e
278      Constr(Value lb,const Expr &e) :
279        _expr(e), _lb(lb), _ub(NaN) {}
280      ///\e
281      Constr(const Expr &e) :
282        _expr(e), _lb(NaN), _ub(NaN) {}
283      ///\e
284      void clear()
285      {
286        _expr.clear();
287        _lb=_ub=NaN;
288      }
289      ///\e
290      Expr &expr() { return _expr; }
291      ///\e
292      const Expr &expr() const { return _expr; }
293      ///\e
294      Value &lowerBound() { return _lb; }
295      ///\e
296      const Value &lowerBound() const { return _lb; }
297      ///\e
298      Value &upperBound() { return _ub; }
299      ///\e
300      const Value &upperBound() const { return _ub; }
301    };
302   
303
304  protected:
305    _FixId rows;
306    _FixId cols;
307
308    /// \e
309    virtual int _addCol() = 0;
310    /// \e
311    virtual int _addRow() = 0;
312    /// \e
313
314    /// \warning Arrays are indexed from 1 (datum at index 0 is ignored)
315    ///
316    virtual void _setRowCoeffs(int i,
317                               int length,
318                               int  const * indices,
319                               Value  const * values ) = 0;
320    /// \e
321
322    /// \warning Arrays are indexed from 1 (datum at index 0 is ignored)
323    ///
324    virtual void _setColCoeffs(int i,
325                               int length,
326                               int  const * indices,
327                               Value  const * values ) = 0;
328   
329    /// \e
330
331    /// The lower bound of a variable (column) have to be given by an
332    /// extended number of type Value, i.e. a finite number of type
333    /// Value or -\ref INF.
334    virtual void _setColLowerBound(int i, Value value) = 0;
335    /// \e
336
337    /// The upper bound of a variable (column) have to be given by an
338    /// extended number of type Value, i.e. a finite number of type
339    /// Value or \ref INF.
340    virtual void _setColUpperBound(int i, Value value) = 0;
341    /// \e
342
343    /// The lower bound of a linear expression (row) have to be given by an
344    /// extended number of type Value, i.e. a finite number of type
345    /// Value or -\ref INF.
346    virtual void _setRowLowerBound(int i, Value value) = 0;
347    /// \e
348
349    /// The upper bound of a linear expression (row) have to be given by an
350    /// extended number of type Value, i.e. a finite number of type
351    /// Value or \ref INF.
352    virtual void _setRowUpperBound(int i, Value value) = 0;
353
354    /// \e
355    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
356
357    ///\e
358   
359    ///\bug Wrong interface
360    ///
361    virtual SolutionType _solve() = 0;
362
363    ///\e
364
365    ///\bug Wrong interface
366    ///
367    virtual Value _getSolution(int i) = 0;
368    ///\e
369
370    ///\bug unimplemented!!!!
371    void clearObj() {}
372  public:
373
374
375    ///\e
376    virtual ~LpSolverBase() {}
377
378    ///\name Building up and modification of the LP
379
380    ///@{
381
382    ///Add a new empty column (i.e a new variable) to the LP
383    Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
384
385    ///\brief Fill the elements of a container with newly created columns
386    ///(i.e a new variables)
387    ///
388    ///This magic function takes a container as its argument
389    ///and fills its elements
390    ///with new columns (i.e. variables)
391    ///\param t can be
392    ///- a standard STL compatible iterable container with
393    ///\ref Col as its \c values_type
394    ///like
395    ///\code
396    ///std::vector<LpSolverBase::Col>
397    ///std::list<LpSolverBase::Col>
398    ///\endcode
399    ///- a standard STL compatible iterable container with
400    ///\ref Col as its \c mapped_type
401    ///like
402    ///\code
403    ///std::map<AnyType,LpSolverBase::Col>
404    ///\endcode
405    ///- an iterable lemon \ref concept::WriteMap "write map" like
406    ///\code
407    ///ListGraph::NodeMap<LpSolverBase::Col>
408    ///ListGraph::EdgeMap<LpSolverBase::Col>
409    ///\endcode
410    ///\return The number of the created column.
411    ///\bug Iterable nodemap hasn't been implemented yet.
412#ifdef DOXYGEN
413    template<class T>
414    int addColSet(T &t) { return 0;}
415#else
416    template<class T>
417    typename enable_if<typename T::value_type::LpSolverCol,int>::type
418    addColSet(T &t,dummy<0> = 0) {
419      int s=0;
420      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
421      return s;
422    }
423    template<class T>
424    typename enable_if<typename T::value_type::second_type::LpSolverCol,
425                       int>::type
426    addColSet(T &t,dummy<1> = 1) {
427      int s=0;
428      for(typename T::iterator i=t.begin();i!=t.end();++i) {
429        i->second=addCol();
430        s++;
431      }
432      return s;
433    }
434    template<class T>
435    typename enable_if<typename T::ValueSet::value_type::LpSolverCol,
436                       int>::type
437    addColSet(T &t,dummy<2> = 2) {
438      ///\bug <tt>return addColSet(t.valueSet());</tt> should also work.
439      int s=0;
440      for(typename T::ValueSet::iterator i=t.valueSet().begin();
441          i!=t.valueSet().end();
442          ++i)
443        {
444          *i=addCol();
445          s++;
446        }
447      return s;
448    }
449#endif
450
451    ///Add a new empty row (i.e a new constaint) to the LP
452
453    ///This function adds a new empty row (i.e a new constaint) to the LP.
454    ///\return The created row
455    Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
456
457    ///Set a row (i.e a constaint) of the LP
458
459    ///\param r is the row to be modified
460    ///\param l is lower bound (-\ref INF means no bound)
461    ///\param e is a linear expression (see \ref Expr)
462    ///\param u is the upper bound (\ref INF means no bound)
463    ///\bug This is a temportary function. The interface will change to
464    ///a better one.
465    void setRow(Row r, Value l,const Expr &e, Value u) {
466      std::vector<int> indices;
467      std::vector<Value> values;
468      indices.push_back(0);
469      values.push_back(0);
470      for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
471        if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
472          indices.push_back(cols.floatingId((*i).first.id));
473          values.push_back((*i).second);
474        }
475      _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
476                    &indices[0],&values[0]);
477      _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
478      _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
479    }
480
481    ///Set a row (i.e a constaint) of the LP
482
483    ///\param r is the row to be modified
484    ///\param c is a linear expression (see \ref Constr)
485    void setRow(Row r, const Constr &c) {
486      setRow(r,
487             isnan(c.lowerBound())?-INF:c.lowerBound(),
488             c.expr(),
489             isnan(c.upperBound())?INF:c.upperBound());
490    }
491
492    ///Add a new row (i.e a new constaint) to the LP
493
494    ///\param l is the lower bound (-\ref INF means no bound)
495    ///\param e is a linear expression (see \ref Expr)
496    ///\param u is the upper bound (\ref INF means no bound)
497    ///\return The created row.
498    ///\bug This is a temportary function. The interface will change to
499    ///a better one.
500    Row addRow(Value l,const Expr &e, Value u) {
501      Row r=addRow();
502      setRow(r,l,e,u);
503      return r;
504    }
505
506    ///Add a new row (i.e a new constaint) to the LP
507
508    ///\param c is a linear expression (see \ref Constr)
509    ///\return The created row.
510    Row addRow(const Constr &c) {
511      Row r=addRow();
512      setRow(r,c);
513      return r;
514    }
515
516    /// Set the lower bound of a column (i.e a variable)
517
518    /// The upper bound of a variable (column) have to be given by an
519    /// extended number of type Value, i.e. a finite number of type
520    /// Value or -\ref INF.
521    virtual void setColLowerBound(Col c, Value value) {
522      _setColLowerBound(cols.floatingId(c.id),value);
523    }
524    /// Set the upper bound of a column (i.e a variable)
525
526    /// The upper bound of a variable (column) have to be given by an
527    /// extended number of type Value, i.e. a finite number of type
528    /// Value or \ref INF.
529    virtual void setColUpperBound(Col c, Value value) {
530      _setColUpperBound(cols.floatingId(c.id),value);
531    };
532    /// Set the lower bound of a row (i.e a constraint)
533
534    /// The lower bound of a linear expression (row) have to be given by an
535    /// extended number of type Value, i.e. a finite number of type
536    /// Value or -\ref INF.
537    virtual void setRowLowerBound(Row r, Value value) {
538      _setRowLowerBound(rows.floatingId(r.id),value);
539    };
540    /// Set the upper bound of a row (i.e a constraint)
541
542    /// The upper bound of a linear expression (row) have to be given by an
543    /// extended number of type Value, i.e. a finite number of type
544    /// Value or \ref INF.
545    virtual void setRowUpperBound(Row r, Value value) {
546      _setRowUpperBound(rows.floatingId(r.id),value);
547    };
548    ///Set an element of the objective function
549    void setObjCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
550    ///Set the objective function
551   
552    ///\param e is a linear expression of type \ref Expr.
553    ///\todo What to do with the constant component?
554    void setObj(Expr e) {
555      clearObj();
556      for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
557        setObjCoeff((*i).first,(*i).second);
558    }
559
560    ///@}
561
562
563    ///\name Solving the LP
564
565    ///@{
566
567    ///\e
568    SolutionType solve() { return _solve(); }
569   
570    ///@}
571   
572    ///\name Obtaining the solution LP
573
574    ///@{
575
576    ///\e
577    Value solution(Col c) { return _getSolution(cols.floatingId(c.id)); }
578
579    ///@}
580   
581  }; 
582
583  ///\e
584 
585  ///\relates LpSolverBase::Expr
586  ///
587  inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
588                                      const LpSolverBase::Expr &b)
589  {
590    LpSolverBase::Expr tmp(a);
591    tmp+=b; ///\todo Don't STL have some special 'merge' algorithm?
592    return tmp;
593  }
594  ///\e
595 
596  ///\relates LpSolverBase::Expr
597  ///
598  inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
599                                      const LpSolverBase::Expr &b)
600  {
601    LpSolverBase::Expr tmp(a);
602    tmp-=b; ///\todo Don't STL have some special 'merge' algorithm?
603    return tmp;
604  }
605  ///\e
606 
607  ///\relates LpSolverBase::Expr
608  ///
609  inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
610                                      const LpSolverBase::Value &b)
611  {
612    LpSolverBase::Expr tmp(a);
613    tmp*=b; ///\todo Don't STL have some special 'merge' algorithm?
614    return tmp;
615  }
616 
617  ///\e
618 
619  ///\relates LpSolverBase::Expr
620  ///
621  inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
622                                      const LpSolverBase::Expr &b)
623  {
624    LpSolverBase::Expr tmp(b);
625    tmp*=a; ///\todo Don't STL have some special 'merge' algorithm?
626    return tmp;
627  }
628  ///\e
629 
630  ///\relates LpSolverBase::Expr
631  ///
632  inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
633                                      const LpSolverBase::Value &b)
634  {
635    LpSolverBase::Expr tmp(a);
636    tmp/=b; ///\todo Don't STL have some special 'merge' algorithm?
637    return tmp;
638  }
639 
640  ///\e
641 
642  ///\relates LpSolverBase::Constr
643  ///
644  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
645                                         const LpSolverBase::Expr &f)
646  {
647    return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
648  }
649
650  ///\e
651 
652  ///\relates LpSolverBase::Constr
653  ///
654  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
655                                         const LpSolverBase::Expr &f)
656  {
657    return LpSolverBase::Constr(e,f);
658  }
659
660  ///\e
661 
662  ///\relates LpSolverBase::Constr
663  ///
664  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
665                                         const LpSolverBase::Value &f)
666  {
667    return LpSolverBase::Constr(e,f);
668  }
669
670  ///\e
671 
672  ///\relates LpSolverBase::Constr
673  ///
674  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
675                                         const LpSolverBase::Expr &f)
676  {
677    return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
678  }
679
680
681  ///\e
682 
683  ///\relates LpSolverBase::Constr
684  ///
685  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
686                                         const LpSolverBase::Expr &f)
687  {
688    return LpSolverBase::Constr(f,e);
689  }
690
691
692  ///\e
693 
694  ///\relates LpSolverBase::Constr
695  ///
696  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
697                                         const LpSolverBase::Value &f)
698  {
699    return LpSolverBase::Constr(f,e);
700  }
701
702  ///\e
703 
704  ///\relates LpSolverBase::Constr
705  ///
706  inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
707                                         const LpSolverBase::Expr &f)
708  {
709    return LpSolverBase::Constr(0,e-f,0);
710  }
711
712  ///\e
713 
714  ///\relates LpSolverBase::Constr
715  ///
716  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
717                                         const LpSolverBase::Constr&c)
718  {
719    LpSolverBase::Constr tmp(c);
720    ///\todo Create an own exception type.
721    if(!isnan(tmp.lowerBound())) throw LogicError();
722    else tmp.lowerBound()=n;
723    return tmp;
724  }
725  ///\e
726 
727  ///\relates LpSolverBase::Constr
728  ///
729  inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
730                                         const LpSolverBase::Value &n)
731  {
732    LpSolverBase::Constr tmp(c);
733    ///\todo Create an own exception type.
734    if(!isnan(tmp.upperBound())) throw LogicError();
735    else tmp.upperBound()=n;
736    return tmp;
737  }
738
739  ///\e
740 
741  ///\relates LpSolverBase::Constr
742  ///
743  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
744                                         const LpSolverBase::Constr&c)
745  {
746    LpSolverBase::Constr tmp(c);
747    ///\todo Create an own exception type.
748    if(!isnan(tmp.upperBound())) throw LogicError();
749    else tmp.upperBound()=n;
750    return tmp;
751  }
752  ///\e
753 
754  ///\relates LpSolverBase::Constr
755  ///
756  inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
757                                         const LpSolverBase::Value &n)
758  {
759    LpSolverBase::Constr tmp(c);
760    ///\todo Create an own exception type.
761    if(!isnan(tmp.lowerBound())) throw LogicError();
762    else tmp.lowerBound()=n;
763    return tmp;
764  }
765
766
767} //namespace lemon
768
769#endif //LEMON_LP_BASE_H
Note: See TracBrowser for help on using the repository browser.