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