COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/lp/lp_base.h @ 1272:17be4c5bc6c6

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