COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/lp/lp_base.h @ 1275:16980bf77bd3

Last change on this file since 1275:16980bf77bd3 was 1275:16980bf77bd3, checked in by Alpar Juttner, 15 years ago

Minor improvements

File size: 20.0 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      ///\e
302      bool lowerBounded() const { return std::isfinite(_lb); }
303      ///\e
304      bool upperBounded() const { return std::isfinite(_ub); }
305    };
306   
307
308  protected:
309    _FixId rows;
310    _FixId cols;
311
312    /// \e
313    virtual int _addCol() = 0;
314    /// \e
315    virtual int _addRow() = 0;
316    /// \e
317
318    /// \warning Arrays are indexed from 1 (datum at index 0 is ignored)
319    ///
320    virtual void _setRowCoeffs(int i,
321                               int length,
322                               int  const * indices,
323                               Value  const * values ) = 0;
324    /// \e
325
326    /// \warning Arrays are indexed from 1 (datum at index 0 is ignored)
327    ///
328    virtual void _setColCoeffs(int i,
329                               int length,
330                               int  const * indices,
331                               Value  const * values ) = 0;
332   
333    /// \e
334
335    /// The lower bound of a variable (column) have to be given by an
336    /// extended number of type Value, i.e. a finite number of type
337    /// Value or -\ref INF.
338    virtual void _setColLowerBound(int i, Value value) = 0;
339    /// \e
340
341    /// The upper bound of a variable (column) have to be given by an
342    /// extended number of type Value, i.e. a finite number of type
343    /// Value or \ref INF.
344    virtual void _setColUpperBound(int i, Value value) = 0;
345    /// \e
346
347    /// The lower bound of a linear expression (row) have to be given by an
348    /// extended number of type Value, i.e. a finite number of type
349    /// Value or -\ref INF.
350    virtual void _setRowLowerBound(int i, Value value) = 0;
351    /// \e
352
353    /// The upper bound of a linear expression (row) have to be given by an
354    /// extended number of type Value, i.e. a finite number of type
355    /// Value or \ref INF.
356    virtual void _setRowUpperBound(int i, Value value) = 0;
357
358    /// \e
359    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
360
361    ///\e
362   
363    ///\bug Wrong interface
364    ///
365    virtual SolutionType _solve() = 0;
366
367    ///\e
368
369    ///\bug Wrong interface
370    ///
371    virtual Value _getSolution(int i) = 0;
372    ///\e
373
374    ///\bug unimplemented!!!!
375    void clearObj() {}
376  public:
377
378
379    ///\e
380    virtual ~LpSolverBase() {}
381
382    ///\name Building up and modification of the LP
383
384    ///@{
385
386    ///Add a new empty column (i.e a new variable) to the LP
387    Col addCol() { Col c; c.id=cols.insert(_addCol()); return c;}
388
389    ///\brief Fill the elements of a container with newly created columns
390    ///(i.e a new variables)
391    ///
392    ///This magic function takes a container as its argument
393    ///and fills its elements
394    ///with new columns (i.e. variables)
395    ///\param t can be
396    ///- a standard STL compatible iterable container with
397    ///\ref Col as its \c values_type
398    ///like
399    ///\code
400    ///std::vector<LpSolverBase::Col>
401    ///std::list<LpSolverBase::Col>
402    ///\endcode
403    ///- a standard STL compatible iterable container with
404    ///\ref Col as its \c mapped_type
405    ///like
406    ///\code
407    ///std::map<AnyType,LpSolverBase::Col>
408    ///\endcode
409    ///- an iterable lemon \ref concept::WriteMap "write map" like
410    ///\code
411    ///ListGraph::NodeMap<LpSolverBase::Col>
412    ///ListGraph::EdgeMap<LpSolverBase::Col>
413    ///\endcode
414    ///\return The number of the created column.
415    ///\bug Iterable nodemap hasn't been implemented yet.
416#ifdef DOXYGEN
417    template<class T>
418    int addColSet(T &t) { return 0;}
419#else
420    template<class T>
421    typename enable_if<typename T::value_type::LpSolverCol,int>::type
422    addColSet(T &t,dummy<0> = 0) {
423      int s=0;
424      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
425      return s;
426    }
427    template<class T>
428    typename enable_if<typename T::value_type::second_type::LpSolverCol,
429                       int>::type
430    addColSet(T &t,dummy<1> = 1) {
431      int s=0;
432      for(typename T::iterator i=t.begin();i!=t.end();++i) {
433        i->second=addCol();
434        s++;
435      }
436      return s;
437    }
438    template<class T>
439    typename enable_if<typename T::ValueSet::value_type::LpSolverCol,
440                       int>::type
441    addColSet(T &t,dummy<2> = 2) {
442      ///\bug <tt>return addColSet(t.valueSet());</tt> should also work.
443      int s=0;
444      for(typename T::ValueSet::iterator i=t.valueSet().begin();
445          i!=t.valueSet().end();
446          ++i)
447        {
448          *i=addCol();
449          s++;
450        }
451      return s;
452    }
453#endif
454
455    ///Add a new empty row (i.e a new constaint) to the LP
456
457    ///This function adds a new empty row (i.e a new constaint) to the LP.
458    ///\return The created row
459    Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
460
461    ///Set a row (i.e a constaint) of the LP
462
463    ///\param r is the row to be modified
464    ///\param l is lower bound (-\ref INF means no bound)
465    ///\param e is a linear expression (see \ref Expr)
466    ///\param u is the upper bound (\ref INF means no bound)
467    ///\bug This is a temportary function. The interface will change to
468    ///a better one.
469    void setRow(Row r, Value l,const Expr &e, Value u) {
470      std::vector<int> indices;
471      std::vector<Value> values;
472      indices.push_back(0);
473      values.push_back(0);
474      for(Expr::const_iterator i=e.begin(); i!=e.end(); ++i)
475        if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
476          indices.push_back(cols.floatingId((*i).first.id));
477          values.push_back((*i).second);
478        }
479      _setRowCoeffs(rows.floatingId(r.id),indices.size()-1,
480                    &indices[0],&values[0]);
481      _setRowLowerBound(rows.floatingId(r.id),l-e.constComp());
482      _setRowUpperBound(rows.floatingId(r.id),u-e.constComp());
483    }
484
485    ///Set a row (i.e a constaint) of the LP
486
487    ///\param r is the row to be modified
488    ///\param c is a linear expression (see \ref Constr)
489    void setRow(Row r, const Constr &c) {
490      setRow(r,
491             c.lowerBounded()?c.lowerBound():-INF,
492             c.expr(),
493             c.upperBounded()?c.upperBound():INF);
494    }
495
496    ///Add a new row (i.e a new constaint) to the LP
497
498    ///\param l is the lower bound (-\ref INF means no bound)
499    ///\param e is a linear expression (see \ref Expr)
500    ///\param u is the upper bound (\ref INF means no bound)
501    ///\return The created row.
502    ///\bug This is a temportary function. The interface will change to
503    ///a better one.
504    Row addRow(Value l,const Expr &e, Value u) {
505      Row r=addRow();
506      setRow(r,l,e,u);
507      return r;
508    }
509
510    ///Add a new row (i.e a new constaint) to the LP
511
512    ///\param c is a linear expression (see \ref Constr)
513    ///\return The created row.
514    Row addRow(const Constr &c) {
515      Row r=addRow();
516      setRow(r,c);
517      return r;
518    }
519
520    /// Set the lower bound of a column (i.e a variable)
521
522    /// The upper bound of a variable (column) have to be given by an
523    /// extended number of type Value, i.e. a finite number of type
524    /// Value or -\ref INF.
525    virtual void setColLowerBound(Col c, Value value) {
526      _setColLowerBound(cols.floatingId(c.id),value);
527    }
528    /// Set the upper bound of a column (i.e a variable)
529
530    /// The upper bound of a variable (column) have to be given by an
531    /// extended number of type Value, i.e. a finite number of type
532    /// Value or \ref INF.
533    virtual void setColUpperBound(Col c, Value value) {
534      _setColUpperBound(cols.floatingId(c.id),value);
535    };
536    /// Set the lower bound of a row (i.e a constraint)
537
538    /// The lower bound of a linear expression (row) have to be given by an
539    /// extended number of type Value, i.e. a finite number of type
540    /// Value or -\ref INF.
541    virtual void setRowLowerBound(Row r, Value value) {
542      _setRowLowerBound(rows.floatingId(r.id),value);
543    };
544    /// Set the upper bound of a row (i.e a constraint)
545
546    /// The upper bound of a linear expression (row) have to be given by an
547    /// extended number of type Value, i.e. a finite number of type
548    /// Value or \ref INF.
549    virtual void setRowUpperBound(Row r, Value value) {
550      _setRowUpperBound(rows.floatingId(r.id),value);
551    };
552    ///Set an element of the objective function
553    void setObjCoeff(Col c, Value v) {_setObjCoeff(cols.floatingId(c.id),v); };
554    ///Set the objective function
555   
556    ///\param e is a linear expression of type \ref Expr.
557    ///\todo What to do with the constant component?
558    void setObj(Expr e) {
559      clearObj();
560      for (Expr::iterator i=e.begin(); i!=e.end(); ++i)
561        setObjCoeff((*i).first,(*i).second);
562    }
563
564    ///@}
565
566
567    ///\name Solving the LP
568
569    ///@{
570
571    ///\e
572    SolutionType solve() { return _solve(); }
573   
574    ///@}
575   
576    ///\name Obtaining the solution LP
577
578    ///@{
579
580    ///\e
581    Value solution(Col c) { return _getSolution(cols.floatingId(c.id)); }
582
583    ///@}
584   
585  }; 
586
587  ///\e
588 
589  ///\relates LpSolverBase::Expr
590  ///
591  inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,
592                                      const LpSolverBase::Expr &b)
593  {
594    LpSolverBase::Expr tmp(a);
595    tmp+=b; ///\todo Don't STL have some special 'merge' algorithm?
596    return tmp;
597  }
598  ///\e
599 
600  ///\relates LpSolverBase::Expr
601  ///
602  inline LpSolverBase::Expr operator-(const LpSolverBase::Expr &a,
603                                      const LpSolverBase::Expr &b)
604  {
605    LpSolverBase::Expr tmp(a);
606    tmp-=b; ///\todo Don't STL have some special 'merge' algorithm?
607    return tmp;
608  }
609  ///\e
610 
611  ///\relates LpSolverBase::Expr
612  ///
613  inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,
614                                      const LpSolverBase::Value &b)
615  {
616    LpSolverBase::Expr tmp(a);
617    tmp*=b; ///\todo Don't STL have some special 'merge' algorithm?
618    return tmp;
619  }
620 
621  ///\e
622 
623  ///\relates LpSolverBase::Expr
624  ///
625  inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,
626                                      const LpSolverBase::Expr &b)
627  {
628    LpSolverBase::Expr tmp(b);
629    tmp*=a; ///\todo Don't STL have some special 'merge' algorithm?
630    return tmp;
631  }
632  ///\e
633 
634  ///\relates LpSolverBase::Expr
635  ///
636  inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,
637                                      const LpSolverBase::Value &b)
638  {
639    LpSolverBase::Expr tmp(a);
640    tmp/=b; ///\todo Don't STL have some special 'merge' algorithm?
641    return tmp;
642  }
643 
644  ///\e
645 
646  ///\relates LpSolverBase::Constr
647  ///
648  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
649                                         const LpSolverBase::Expr &f)
650  {
651    return LpSolverBase::Constr(-LpSolverBase::INF,e-f,0);
652  }
653
654  ///\e
655 
656  ///\relates LpSolverBase::Constr
657  ///
658  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,
659                                         const LpSolverBase::Expr &f)
660  {
661    return LpSolverBase::Constr(e,f);
662  }
663
664  ///\e
665 
666  ///\relates LpSolverBase::Constr
667  ///
668  inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,
669                                         const LpSolverBase::Value &f)
670  {
671    return LpSolverBase::Constr(e,f);
672  }
673
674  ///\e
675 
676  ///\relates LpSolverBase::Constr
677  ///
678  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
679                                         const LpSolverBase::Expr &f)
680  {
681    return LpSolverBase::Constr(-LpSolverBase::INF,f-e,0);
682  }
683
684
685  ///\e
686 
687  ///\relates LpSolverBase::Constr
688  ///
689  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,
690                                         const LpSolverBase::Expr &f)
691  {
692    return LpSolverBase::Constr(f,e);
693  }
694
695
696  ///\e
697 
698  ///\relates LpSolverBase::Constr
699  ///
700  inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,
701                                         const LpSolverBase::Value &f)
702  {
703    return LpSolverBase::Constr(f,e);
704  }
705
706  ///\e
707 
708  ///\relates LpSolverBase::Constr
709  ///
710  inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,
711                                         const LpSolverBase::Expr &f)
712  {
713    return LpSolverBase::Constr(0,e-f,0);
714  }
715
716  ///\e
717 
718  ///\relates LpSolverBase::Constr
719  ///
720  inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,
721                                         const LpSolverBase::Constr&c)
722  {
723    LpSolverBase::Constr tmp(c);
724    ///\todo Create an own exception type.
725    if(!isnan(tmp.lowerBound())) throw LogicError();
726    else tmp.lowerBound()=n;
727    return tmp;
728  }
729  ///\e
730 
731  ///\relates LpSolverBase::Constr
732  ///
733  inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,
734                                         const LpSolverBase::Value &n)
735  {
736    LpSolverBase::Constr tmp(c);
737    ///\todo Create an own exception type.
738    if(!isnan(tmp.upperBound())) throw LogicError();
739    else tmp.upperBound()=n;
740    return tmp;
741  }
742
743  ///\e
744 
745  ///\relates LpSolverBase::Constr
746  ///
747  inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,
748                                         const LpSolverBase::Constr&c)
749  {
750    LpSolverBase::Constr tmp(c);
751    ///\todo Create an own exception type.
752    if(!isnan(tmp.upperBound())) throw LogicError();
753    else tmp.upperBound()=n;
754    return tmp;
755  }
756  ///\e
757 
758  ///\relates LpSolverBase::Constr
759  ///
760  inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,
761                                         const LpSolverBase::Value &n)
762  {
763    LpSolverBase::Constr tmp(c);
764    ///\todo Create an own exception type.
765    if(!isnan(tmp.lowerBound())) throw LogicError();
766    else tmp.lowerBound()=n;
767    return tmp;
768  }
769
770
771} //namespace lemon
772
773#endif //LEMON_LP_BASE_H
Note: See TracBrowser for help on using the repository browser.