Changeset 459:ed54c0d13df0 in lemon1.2 for lemon/lp_base.h
 Timestamp:
 12/02/08 22:48:28 (12 years ago)
 Branch:
 default
 Children:
 460:76ec7bd57026, 502:17cabb114d52
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/lp_base.h
r458 r459 26 26 #include<lemon/math.h> 27 27 28 #include<lemon/error.h> 29 #include<lemon/assert.h> 30 28 31 #include<lemon/core.h> 29 #include<lemon/bits/ lp_id.h>32 #include<lemon/bits/solver_bits.h> 30 33 31 34 ///\file … … 34 37 namespace lemon { 35 38 36 /// Function to decide whether a floating point value is finite or not. 37 38 /// Retruns true if the argument is not infinity, minus infinity or NaN. 39 /// It does the same as the isfinite() function defined by C99. 40 template <typename T> 41 bool isFinite(T value) 42 { 43 typedef std::numeric_limits<T> Lim; 44 if ((Lim::has_infinity && (value == Lim::infinity()  value == 45 Lim::infinity()))  46 ((Lim::has_quiet_NaN  Lim::has_signaling_NaN) && value != value)) 47 { 48 return false; 49 } 50 return true; 51 } 52 53 ///Common base class for LP solvers 54 55 ///\todo Much more docs 39 ///Common base class for LP and MIP solvers 40 41 ///Usually this class is not used directly, please use one of the concrete 42 ///implementations of the solver interface. 56 43 ///\ingroup lp_group 57 class Lp SolverBase {44 class LpBase { 58 45 59 46 protected: 60 47 61 _ lp_bits::LpIdrows;62 _ lp_bits::LpIdcols;48 _solver_bits::VarIndex rows; 49 _solver_bits::VarIndex cols; 63 50 64 51 public: … … 75 62 }; 76 63 77 ///\e 78 enum SolutionStatus { 79 ///Feasible solution hasn't been found (but may exist). 80 81 ///\todo NOTFOUND might be a better name. 82 /// 83 UNDEFINED = 0, 84 ///The problem has no feasible solution 85 INFEASIBLE = 1, 86 ///Feasible solution found 87 FEASIBLE = 2, 88 ///Optimal solution exists and found 89 OPTIMAL = 3, 90 ///The cost function is unbounded 91 92 ///\todo Give a feasible solution and an infinite ray (and the 93 ///corresponding bases) 94 INFINITE = 4 95 }; 96 97 ///\e The type of the investigated LP problem 98 enum ProblemTypes { 99 ///Primaldual feasible 100 PRIMAL_DUAL_FEASIBLE = 0, 101 ///Primal feasible dual infeasible 102 PRIMAL_FEASIBLE_DUAL_INFEASIBLE = 1, 103 ///Primal infeasible dual feasible 104 PRIMAL_INFEASIBLE_DUAL_FEASIBLE = 2, 105 ///Primaldual infeasible 106 PRIMAL_DUAL_INFEASIBLE = 3, 107 ///Could not determine so far 108 UNKNOWN = 4 64 ///Direction of the optimization 65 enum Sense { 66 /// Minimization 67 MIN, 68 /// Maximization 69 MAX 109 70 }; 110 71 … … 116 77 static const Value NaN; 117 78 118 static inline bool isNaN(const Value& v) { return v!=v; }119 120 79 friend class Col; 121 80 friend class ColIt; 122 81 friend class Row; 82 friend class RowIt; 123 83 124 84 ///Refer to a column of the LP. … … 129 89 ///other columns. 130 90 /// 131 ///\ todo Document what can one do with a Col (INVALID, comparing,132 /// it is similar to Node/Edge)91 ///\note This class is similar to other Item types in LEMON, like 92 ///Node and Arc types in digraph. 133 93 class Col { 94 friend class LpBase; 134 95 protected: 135 int id; 136 friend class LpSolverBase; 137 friend class MipSolverBase; 138 explicit Col(int _id) : id(_id) {} 96 int _id; 97 explicit Col(int id) : _id(id) {} 139 98 public: 140 99 typedef Value ExprValue; 141 typedef True LpSolverCol; 100 typedef True LpCol; 101 /// Default constructor 102 103 /// \warning The default constructor sets the Col to an 104 /// undefined value. 142 105 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 bool operator!=(Col c) const {return id!=c.id;} 106 /// Invalid constructor \& conversion. 107 108 /// This constructor initializes the Col to be invalid. 109 /// \sa Invalid for more details. 110 Col(const Invalid&) : _id(1) {} 111 /// Equality operator 112 113 /// Two \ref Col "Col"s are equal if and only if they point to 114 /// the same LP column or both are invalid. 115 bool operator==(Col c) const {return _id == c._id;} 116 /// Inequality operator 117 118 /// \sa operator==(Col c) 119 /// 120 bool operator!=(Col c) const {return _id != c._id;} 121 /// Artificial ordering operator. 122 123 /// To allow the use of this object in std::map or similar 124 /// associative container we require this. 125 /// 126 /// \note This operator only have to define some strict ordering of 127 /// the items; this order has nothing to do with the iteration 128 /// ordering of the items. 129 bool operator<(Col c) const {return _id < c._id;} 148 130 }; 149 131 132 ///Iterator for iterate over the columns of an LP problem 133 134 /// Its usage is quite simple, for example you can count the number 135 /// of columns in an LP \c lp: 136 ///\code 137 /// int count=0; 138 /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count; 139 ///\endcode 150 140 class ColIt : public Col { 151 const Lp SolverBase *_lp;141 const LpBase *_solver; 152 142 public: 143 /// Default constructor 144 145 /// \warning The default constructor sets the iterator 146 /// to an undefined value. 153 147 ColIt() {} 154 ColIt(const LpSolverBase &lp) : _lp(&lp) 148 /// Sets the iterator to the first Col 149 150 /// Sets the iterator to the first Col. 151 /// 152 ColIt(const LpBase &solver) : _solver(&solver) 155 153 { 156 _lp>cols.firstFix(id); 157 } 154 _solver>cols.firstItem(_id); 155 } 156 /// Invalid constructor \& conversion 157 158 /// Initialize the iterator to be invalid. 159 /// \sa Invalid for more details. 158 160 ColIt(const Invalid&) : Col(INVALID) {} 161 /// Next column 162 163 /// Assign the iterator to the next column. 164 /// 159 165 ColIt &operator++() 160 166 { 161 _ lp>cols.nextFix(id);167 _solver>cols.nextItem(_id); 162 168 return *this; 163 169 } 164 170 }; 165 171 166 static int id(const Col& col) { return col.id; } 167 172 /// \brief Returns the ID of the column. 173 static int id(const Col& col) { return col._id; } 174 /// \brief Returns the column with the given ID. 175 /// 176 /// \pre The argument should be a valid column ID in the LP problem. 177 static Col colFromId(int id) { return Col(id); } 168 178 169 179 ///Refer to a row of the LP. … … 174 184 ///other rows. 175 185 /// 176 ///\ todo Document what can one do with a Row (INVALID, comparing,177 /// it is similar to Node/Edge)186 ///\note This class is similar to other Item types in LEMON, like 187 ///Node and Arc types in digraph. 178 188 class Row { 189 friend class LpBase; 179 190 protected: 180 int id; 181 friend class LpSolverBase; 182 explicit Row(int _id) : id(_id) {} 191 int _id; 192 explicit Row(int id) : _id(id) {} 183 193 public: 184 194 typedef Value ExprValue; 185 typedef True LpSolverRow; 195 typedef True LpRow; 196 /// Default constructor 197 198 /// \warning The default constructor sets the Row to an 199 /// undefined value. 186 200 Row() {} 187 Row(const Invalid&) : id(1) {} 188 189 bool operator< (Row c) const {return id< c.id;} 190 bool operator> (Row c) const {return id> c.id;} 191 bool operator==(Row c) const {return id==c.id;} 192 bool operator!=(Row c) const {return id!=c.id;} 201 /// Invalid constructor \& conversion. 202 203 /// This constructor initializes the Row to be invalid. 204 /// \sa Invalid for more details. 205 Row(const Invalid&) : _id(1) {} 206 /// Equality operator 207 208 /// Two \ref Row "Row"s are equal if and only if they point to 209 /// the same LP row or both are invalid. 210 bool operator==(Row r) const {return _id == r._id;} 211 /// Inequality operator 212 213 /// \sa operator==(Row r) 214 /// 215 bool operator!=(Row r) const {return _id != r._id;} 216 /// Artificial ordering operator. 217 218 /// To allow the use of this object in std::map or similar 219 /// associative container we require this. 220 /// 221 /// \note This operator only have to define some strict ordering of 222 /// the items; this order has nothing to do with the iteration 223 /// ordering of the items. 224 bool operator<(Row r) const {return _id < r._id;} 193 225 }; 194 226 227 ///Iterator for iterate over the rows of an LP problem 228 229 /// Its usage is quite simple, for example you can count the number 230 /// of rows in an LP \c lp: 231 ///\code 232 /// int count=0; 233 /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count; 234 ///\endcode 195 235 class RowIt : public Row { 196 const Lp SolverBase *_lp;236 const LpBase *_solver; 197 237 public: 238 /// Default constructor 239 240 /// \warning The default constructor sets the iterator 241 /// to an undefined value. 198 242 RowIt() {} 199 RowIt(const LpSolverBase &lp) : _lp(&lp) 243 /// Sets the iterator to the first Row 244 245 /// Sets the iterator to the first Row. 246 /// 247 RowIt(const LpBase &solver) : _solver(&solver) 200 248 { 201 _lp>rows.firstFix(id); 202 } 249 _solver>rows.firstItem(_id); 250 } 251 /// Invalid constructor \& conversion 252 253 /// Initialize the iterator to be invalid. 254 /// \sa Invalid for more details. 203 255 RowIt(const Invalid&) : Row(INVALID) {} 256 /// Next row 257 258 /// Assign the iterator to the next row. 259 /// 204 260 RowIt &operator++() 205 261 { 206 _ lp>rows.nextFix(id);262 _solver>rows.nextItem(_id); 207 263 return *this; 208 264 } 209 265 }; 210 266 211 static int id(const Row& row) { return row.id; } 212 213 protected: 214 215 int _lpId(const Col& c) const { 216 return cols.floatingId(id(c)); 217 } 218 219 int _lpId(const Row& r) const { 220 return rows.floatingId(id(r)); 221 } 222 223 Col _item(int i, Col) const { 224 return Col(cols.fixId(i)); 225 } 226 227 Row _item(int i, Row) const { 228 return Row(rows.fixId(i)); 229 } 230 267 /// \brief Returns the ID of the row. 268 static int id(const Row& row) { return row._id; } 269 /// \brief Returns the row with the given ID. 270 /// 271 /// \pre The argument should be a valid row ID in the LP problem. 272 static Row rowFromId(int id) { return Row(id); } 231 273 232 274 public: … … 239 281 ///There are several ways to access and modify the contents of this 240 282 ///container. 241 /// Its it fully compatible with \c std::map<Col,double>, so for expamle242 ///if \c e is an Expr and \c v and \c w are of type \ref Col, then you can243 ///read and modify the coefficients like244 ///these.245 283 ///\code 246 284 ///e[v]=5; … … 251 289 ///\code 252 290 ///double s=0; 253 ///for(Lp SolverBase::Expr::iterator i=e.begin();i!=e.end();++i)254 /// s+= i>second;291 ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i) 292 /// s+=*i * primal(i); 255 293 ///\endcode 256 ///(This code computes the sum of all coefficients).294 ///(This code computes the primal value of the expression). 257 295 /// Numbers (<tt>double</tt>'s) 258 296 ///and variables (\ref Col "Col"s) directly convert to an … … 263 301 ///v*2.1+(3*v+(v*12+w+6)*3)/2 264 302 ///\endcode 265 ///are valid \ref Expr "Expr"essions.303 ///are valid expressions. 266 304 ///The usual assignment operations are also defined. 267 305 ///\code … … 271 309 ///e/=5; 272 310 ///\endcode 273 /// The constant member can be set and read by \ref constComp() 311 /// The constant member can be set and read by dereference 312 /// operator (unary *) 313 /// 274 314 ///\code 275 /// e.constComp()=12;276 ///double c= e.constComp();315 ///*e=12; 316 ///double c=*e; 277 317 ///\endcode 278 318 /// 279 ///\note \ref clear() not only sets all coefficients to 0 but also280 ///clears the constant components.281 ///282 319 ///\sa Constr 283 /// 284 class Expr : public std::map<Col,Value> 285 { 320 class Expr { 321 friend class LpBase; 286 322 public: 287 typedef LpSolverBase::Col Key; 288 typedef LpSolverBase::Value Value; 323 /// The key type of the expression 324 typedef LpBase::Col Key; 325 /// The value type of the expression 326 typedef LpBase::Value Value; 289 327 290 328 protected: 291 typedef std::map<Col,Value> Base;292 293 329 Value const_comp; 330 std::map<int, Value> comps; 331 294 332 public: 295 typedef True IsLinExpression; 296 ///\e 297 Expr() : Base(), const_comp(0) { } 298 ///\e 299 Expr(const Key &v) : const_comp(0) { 300 Base::insert(std::make_pair(v, 1)); 301 } 302 ///\e 333 typedef True SolverExpr; 334 /// Default constructor 335 336 /// Construct an empty expression, the coefficients and 337 /// the constant component are initialized to zero. 338 Expr() : const_comp(0) {} 339 /// Construct an expression from a column 340 341 /// Construct an expression, which has a term with \c c variable 342 /// and 1.0 coefficient. 343 Expr(const Col &c) : const_comp(0) { 344 typedef std::map<int, Value>::value_type pair_type; 345 comps.insert(pair_type(id(c), 1)); 346 } 347 /// Construct an expression from a constant 348 349 /// Construct an expression, which's constant component is \c v. 350 /// 303 351 Expr(const Value &v) : const_comp(v) {} 304 ///\e 305 void set(const Key &v,const Value &c) { 306 Base::insert(std::make_pair(v, c)); 307 } 308 ///\e 309 Value &constComp() { return const_comp; } 310 ///\e 311 const Value &constComp() const { return const_comp; } 312 313 ///Removes the components with zero coefficient. 314 void simplify() { 315 for (Base::iterator i=Base::begin(); i!=Base::end();) { 316 Base::iterator j=i; 317 ++j; 318 if ((*i).second==0) Base::erase(i); 319 i=j; 352 /// Returns the coefficient of the column 353 Value operator[](const Col& c) const { 354 std::map<int, Value>::const_iterator it=comps.find(id(c)); 355 if (it != comps.end()) { 356 return it>second; 357 } else { 358 return 0; 320 359 } 321 360 } 322 323 void simplify() const { 324 const_cast<Expr*>(this)>simplify(); 325 } 326 327 ///Removes the coefficients closer to zero than \c tolerance. 328 void simplify(double &tolerance) { 329 for (Base::iterator i=Base::begin(); i!=Base::end();) { 330 Base::iterator j=i; 331 ++j; 332 if (std::fabs((*i).second)<tolerance) Base::erase(i); 333 i=j; 361 /// Returns the coefficient of the column 362 Value& operator[](const Col& c) { 363 return comps[id(c)]; 364 } 365 /// Sets the coefficient of the column 366 void set(const Col &c, const Value &v) { 367 if (v != 0.0) { 368 typedef std::map<int, Value>::value_type pair_type; 369 comps.insert(pair_type(id(c), v)); 370 } else { 371 comps.erase(id(c)); 334 372 } 373 } 374 /// Returns the constant component of the expression 375 Value& operator*() { return const_comp; } 376 /// Returns the constant component of the expression 377 const Value& operator*() const { return const_comp; } 378 /// \brief Removes the coefficients which's absolute value does 379 /// not exceed \c epsilon. It also sets to zero the constant 380 /// component, if it does not exceed epsilon in absolute value. 381 void simplify(Value epsilon = 0.0) { 382 std::map<int, Value>::iterator it=comps.begin(); 383 while (it != comps.end()) { 384 std::map<int, Value>::iterator jt=it; 385 ++jt; 386 if (std::fabs((*it).second) <= epsilon) comps.erase(it); 387 it=jt; 388 } 389 if (std::fabs(const_comp) <= epsilon) const_comp = 0; 390 } 391 392 void simplify(Value epsilon = 0.0) const { 393 const_cast<Expr*>(this)>simplify(epsilon); 335 394 } 336 395 337 396 ///Sets all coefficients and the constant component to 0. 338 397 void clear() { 339 Base::clear();398 comps.clear(); 340 399 const_comp=0; 341 400 } 342 401 343 /// \e402 ///Compound assignment 344 403 Expr &operator+=(const Expr &e) { 345 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) 346 (*this)[j>first]+=j>second; 404 for (std::map<int, Value>::const_iterator it=e.comps.begin(); 405 it!=e.comps.end(); ++it) 406 comps[it>first]+=it>second; 347 407 const_comp+=e.const_comp; 348 408 return *this; 349 409 } 350 /// \e410 ///Compound assignment 351 411 Expr &operator=(const Expr &e) { 352 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) 353 (*this)[j>first]=j>second; 412 for (std::map<int, Value>::const_iterator it=e.comps.begin(); 413 it!=e.comps.end(); ++it) 414 comps[it>first]=it>second; 354 415 const_comp=e.const_comp; 355 416 return *this; 356 417 } 357 ///\e 358 Expr &operator*=(const Value &c) { 359 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) 360 j>second*=c; 361 const_comp*=c; 418 ///Multiply with a constant 419 Expr &operator*=(const Value &v) { 420 for (std::map<int, Value>::iterator it=comps.begin(); 421 it!=comps.end(); ++it) 422 it>second*=v; 423 const_comp*=v; 362 424 return *this; 363 425 } 364 /// \e426 ///Division with a constant 365 427 Expr &operator/=(const Value &c) { 366 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) 367 j>second/=c; 428 for (std::map<int, Value>::iterator it=comps.begin(); 429 it!=comps.end(); ++it) 430 it>second/=c; 368 431 const_comp/=c; 369 432 return *this; 370 433 } 434 435 ///Iterator over the expression 436 437 ///The iterator iterates over the terms of the expression. 438 /// 439 ///\code 440 ///double s=0; 441 ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i) 442 /// s+= *i * primal(i); 443 ///\endcode 444 class CoeffIt { 445 private: 446 447 std::map<int, Value>::iterator _it, _end; 448 449 public: 450 451 /// Sets the iterator to the first term 452 453 /// Sets the iterator to the first term of the expression. 454 /// 455 CoeffIt(Expr& e) 456 : _it(e.comps.begin()), _end(e.comps.end()){} 457 458 /// Convert the iterator to the column of the term 459 operator Col() const { 460 return colFromId(_it>first); 461 } 462 463 /// Returns the coefficient of the term 464 Value& operator*() { return _it>second; } 465 466 /// Returns the coefficient of the term 467 const Value& operator*() const { return _it>second; } 468 /// Next term 469 470 /// Assign the iterator to the next term. 471 /// 472 CoeffIt& operator++() { ++_it; return *this; } 473 474 /// Equality operator 475 bool operator==(Invalid) const { return _it == _end; } 476 /// Inequality operator 477 bool operator!=(Invalid) const { return _it != _end; } 478 }; 479 480 /// Const iterator over the expression 481 482 ///The iterator iterates over the terms of the expression. 483 /// 484 ///\code 485 ///double s=0; 486 ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i) 487 /// s+=*i * primal(i); 488 ///\endcode 489 class ConstCoeffIt { 490 private: 491 492 std::map<int, Value>::const_iterator _it, _end; 493 494 public: 495 496 /// Sets the iterator to the first term 497 498 /// Sets the iterator to the first term of the expression. 499 /// 500 ConstCoeffIt(const Expr& e) 501 : _it(e.comps.begin()), _end(e.comps.end()){} 502 503 /// Convert the iterator to the column of the term 504 operator Col() const { 505 return colFromId(_it>first); 506 } 507 508 /// Returns the coefficient of the term 509 const Value& operator*() const { return _it>second; } 510 511 /// Next term 512 513 /// Assign the iterator to the next term. 514 /// 515 ConstCoeffIt& operator++() { ++_it; return *this; } 516 517 /// Equality operator 518 bool operator==(Invalid) const { return _it == _end; } 519 /// Inequality operator 520 bool operator!=(Invalid) const { return _it != _end; } 521 }; 371 522 372 523 }; … … 395 546 /// e>=t 396 547 ///\endcode 397 ///\warning The validity of a constraint is checked only at run time, so398 /// e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will compile, but will throw399 /// an assertion.548 ///\warning The validity of a constraint is checked only at run 549 ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will 550 ///compile, but will fail an assertion. 400 551 class Constr 401 552 { 402 553 public: 403 typedef Lp SolverBase::Expr Expr;554 typedef LpBase::Expr Expr; 404 555 typedef Expr::Key Key; 405 556 typedef Expr::Value Value; … … 412 563 Constr() : _expr(), _lb(NaN), _ub(NaN) {} 413 564 ///\e 414 Constr(Value lb, const Expr &e,Value ub) :565 Constr(Value lb, const Expr &e, Value ub) : 415 566 _expr(e), _lb(lb), _ub(ub) {} 416 ///\e417 Constr(const Expr &e,Value ub) :418 _expr(e), _lb(NaN), _ub(ub) {}419 ///\e420 Constr(Value lb,const Expr &e) :421 _expr(e), _lb(lb), _ub(NaN) {}422 ///\e423 567 Constr(const Expr &e) : 424 568 _expr(e), _lb(NaN), _ub(NaN) {} … … 454 598 ///Is the constraint lower bounded? 455 599 bool lowerBounded() const { 456 return isFinite(_lb);600 return _lb != INF && !std::isnan(_lb); 457 601 } 458 602 ///Is the constraint upper bounded? 459 603 bool upperBounded() const { 460 return isFinite(_ub);604 return _ub != INF && !std::isnan(_ub); 461 605 } 462 606 … … 471 615 ///There are several ways to access and modify the contents of this 472 616 ///container. 473 /// Its it fully compatible with \c std::map<Row,double>, so for expamle474 ///if \c e is an DualExpr and \c v475 ///and \c w are of type \ref Row, then you can476 ///read and modify the coefficients like477 ///these.478 617 ///\code 479 618 ///e[v]=5; … … 484 623 ///\code 485 624 ///double s=0; 486 ///for(Lp SolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)487 /// s+= i>second;625 ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i) 626 /// s+=*i; 488 627 ///\endcode 489 628 ///(This code computes the sum of all coefficients). … … 496 635 ///v*2.1+(3*v+(v*12+w)*3)/2 497 636 ///\endcode 498 ///are valid \ref DualExpr "DualExpr"essions.637 ///are valid \ref DualExpr dual expressions. 499 638 ///The usual assignment operations are also defined. 500 639 ///\code … … 506 645 /// 507 646 ///\sa Expr 508 /// 509 class DualExpr : public std::map<Row,Value> 510 { 647 class DualExpr { 648 friend class LpBase; 511 649 public: 512 typedef LpSolverBase::Row Key; 513 typedef LpSolverBase::Value Value; 650 /// The key type of the expression 651 typedef LpBase::Row Key; 652 /// The value type of the expression 653 typedef LpBase::Value Value; 514 654 515 655 protected: 516 typedef std::map<Row,Value> Base;656 std::map<int, Value> comps; 517 657 518 658 public: 519 typedef True IsLinExpression; 520 ///\e 521 DualExpr() : Base() { } 522 ///\e 523 DualExpr(const Key &v) { 524 Base::insert(std::make_pair(v, 1)); 525 } 526 ///\e 527 void set(const Key &v,const Value &c) { 528 Base::insert(std::make_pair(v, c)); 529 } 530 531 ///Removes the components with zero coefficient. 532 void simplify() { 533 for (Base::iterator i=Base::begin(); i!=Base::end();) { 534 Base::iterator j=i; 535 ++j; 536 if ((*i).second==0) Base::erase(i); 537 i=j; 659 typedef True SolverExpr; 660 /// Default constructor 661 662 /// Construct an empty expression, the coefficients are 663 /// initialized to zero. 664 DualExpr() {} 665 /// Construct an expression from a row 666 667 /// Construct an expression, which has a term with \c r dual 668 /// variable and 1.0 coefficient. 669 DualExpr(const Row &r) { 670 typedef std::map<int, Value>::value_type pair_type; 671 comps.insert(pair_type(id(r), 1)); 672 } 673 /// Returns the coefficient of the row 674 Value operator[](const Row& r) const { 675 std::map<int, Value>::const_iterator it = comps.find(id(r)); 676 if (it != comps.end()) { 677 return it>second; 678 } else { 679 return 0; 538 680 } 539 681 } 540 541 void simplify() const { 542 const_cast<DualExpr*>(this)>simplify(); 543 } 544 545 ///Removes the coefficients closer to zero than \c tolerance. 546 void simplify(double &tolerance) { 547 for (Base::iterator i=Base::begin(); i!=Base::end();) { 548 Base::iterator j=i; 549 ++j; 550 if (std::fabs((*i).second)<tolerance) Base::erase(i); 551 i=j; 682 /// Returns the coefficient of the row 683 Value& operator[](const Row& r) { 684 return comps[id(r)]; 685 } 686 /// Sets the coefficient of the row 687 void set(const Row &r, const Value &v) { 688 if (v != 0.0) { 689 typedef std::map<int, Value>::value_type pair_type; 690 comps.insert(pair_type(id(r), v)); 691 } else { 692 comps.erase(id(r)); 552 693 } 694 } 695 /// \brief Removes the coefficients which's absolute value does 696 /// not exceed \c epsilon. 697 void simplify(Value epsilon = 0.0) { 698 std::map<int, Value>::iterator it=comps.begin(); 699 while (it != comps.end()) { 700 std::map<int, Value>::iterator jt=it; 701 ++jt; 702 if (std::fabs((*it).second) <= epsilon) comps.erase(it); 703 it=jt; 704 } 705 } 706 707 void simplify(Value epsilon = 0.0) const { 708 const_cast<DualExpr*>(this)>simplify(epsilon); 553 709 } 554 710 555 711 ///Sets all coefficients to 0. 556 712 void clear() { 557 Base::clear(); 558 } 559 560 ///\e 713 comps.clear(); 714 } 715 ///Compound assignment 561 716 DualExpr &operator+=(const DualExpr &e) { 562 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) 563 (*this)[j>first]+=j>second; 717 for (std::map<int, Value>::const_iterator it=e.comps.begin(); 718 it!=e.comps.end(); ++it) 719 comps[it>first]+=it>second; 564 720 return *this; 565 721 } 566 /// \e722 ///Compound assignment 567 723 DualExpr &operator=(const DualExpr &e) { 568 for (Base::const_iterator j=e.begin(); j!=e.end(); ++j) 569 (*this)[j>first]=j>second; 724 for (std::map<int, Value>::const_iterator it=e.comps.begin(); 725 it!=e.comps.end(); ++it) 726 comps[it>first]=it>second; 570 727 return *this; 571 728 } 572 ///\e 573 DualExpr &operator*=(const Value &c) { 574 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) 575 j>second*=c; 729 ///Multiply with a constant 730 DualExpr &operator*=(const Value &v) { 731 for (std::map<int, Value>::iterator it=comps.begin(); 732 it!=comps.end(); ++it) 733 it>second*=v; 576 734 return *this; 577 735 } 578 ///\e 579 DualExpr &operator/=(const Value &c) { 580 for (Base::iterator j=Base::begin(); j!=Base::end(); ++j) 581 j>second/=c; 736 ///Division with a constant 737 DualExpr &operator/=(const Value &v) { 738 for (std::map<int, Value>::iterator it=comps.begin(); 739 it!=comps.end(); ++it) 740 it>second/=v; 582 741 return *this; 583 742 } 743 744 ///Iterator over the expression 745 746 ///The iterator iterates over the terms of the expression. 747 /// 748 ///\code 749 ///double s=0; 750 ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i) 751 /// s+= *i * dual(i); 752 ///\endcode 753 class CoeffIt { 754 private: 755 756 std::map<int, Value>::iterator _it, _end; 757 758 public: 759 760 /// Sets the iterator to the first term 761 762 /// Sets the iterator to the first term of the expression. 763 /// 764 CoeffIt(DualExpr& e) 765 : _it(e.comps.begin()), _end(e.comps.end()){} 766 767 /// Convert the iterator to the row of the term 768 operator Row() const { 769 return rowFromId(_it>first); 770 } 771 772 /// Returns the coefficient of the term 773 Value& operator*() { return _it>second; } 774 775 /// Returns the coefficient of the term 776 const Value& operator*() const { return _it>second; } 777 778 /// Next term 779 780 /// Assign the iterator to the next term. 781 /// 782 CoeffIt& operator++() { ++_it; return *this; } 783 784 /// Equality operator 785 bool operator==(Invalid) const { return _it == _end; } 786 /// Inequality operator 787 bool operator!=(Invalid) const { return _it != _end; } 788 }; 789 790 ///Iterator over the expression 791 792 ///The iterator iterates over the terms of the expression. 793 /// 794 ///\code 795 ///double s=0; 796 ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i) 797 /// s+= *i * dual(i); 798 ///\endcode 799 class ConstCoeffIt { 800 private: 801 802 std::map<int, Value>::const_iterator _it, _end; 803 804 public: 805 806 /// Sets the iterator to the first term 807 808 /// Sets the iterator to the first term of the expression. 809 /// 810 ConstCoeffIt(const DualExpr& e) 811 : _it(e.comps.begin()), _end(e.comps.end()){} 812 813 /// Convert the iterator to the row of the term 814 operator Row() const { 815 return rowFromId(_it>first); 816 } 817 818 /// Returns the coefficient of the term 819 const Value& operator*() const { return _it>second; } 820 821 /// Next term 822 823 /// Assign the iterator to the next term. 824 /// 825 ConstCoeffIt& operator++() { ++_it; return *this; } 826 827 /// Equality operator 828 bool operator==(Invalid) const { return _it == _end; } 829 /// Inequality operator 830 bool operator!=(Invalid) const { return _it != _end; } 831 }; 584 832 }; 585 833 586 834 587 private: 588 589 template <typename _Expr> 590 class MappedOutputIterator { 835 protected: 836 837 class InsertIterator { 838 private: 839 840 std::map<int, Value>& _host; 841 const _solver_bits::VarIndex& _index; 842 591 843 public: 592 593 typedef std::insert_iterator<_Expr> Base;594 844 595 845 typedef std::output_iterator_tag iterator_category; … … 599 849 typedef void pointer; 600 850 601 MappedOutputIterator(const Base& _base, const LpSolverBase& _lp) 602 : base(_base), lp(_lp) {} 603 604 MappedOutputIterator& operator*() { 851 InsertIterator(std::map<int, Value>& host, 852 const _solver_bits::VarIndex& index) 853 : _host(host), _index(index) {} 854 855 InsertIterator& operator=(const std::pair<int, Value>& value) { 856 typedef std::map<int, Value>::value_type pair_type; 857 _host.insert(pair_type(_index[value.first], value.second)); 605 858 return *this; 606 859 } 607 860 608 MappedOutputIterator& operator=(const std::pair<int, Value>& value) { 609 *base = std::make_pair(lp._item(value.first, typename _Expr::Key()), 610 value.second); 611 return *this; 612 } 613 614 MappedOutputIterator& operator++() { 615 ++base; 616 return *this; 617 } 618 619 MappedOutputIterator operator++(int) { 620 MappedOutputIterator tmp(*this); 621 ++base; 622 return tmp; 623 } 624 625 bool operator==(const MappedOutputIterator& it) const { 626 return base == it.base; 627 } 628 629 bool operator!=(const MappedOutputIterator& it) const { 630 return base != it.base; 631 } 632 861 InsertIterator& operator*() { return *this; } 862 InsertIterator& operator++() { return *this; } 863 InsertIterator operator++(int) { return *this; } 864 865 }; 866 867 class ExprIterator { 633 868 private: 634 Base base; 635 const LpSolverBase& lp; 636 }; 637 638 template <typename Expr> 639 class MappedInputIterator { 869 std::map<int, Value>::const_iterator _host_it; 870 const _solver_bits::VarIndex& _index; 640 871 public: 641 872 642 typedef typename Expr::const_iterator Base; 643 644 typedef typename Base::iterator_category iterator_category; 645 typedef typename Base::difference_type difference_type; 873 typedef std::bidirectional_iterator_tag iterator_category; 874 typedef std::ptrdiff_t difference_type; 646 875 typedef const std::pair<int, Value> value_type; 647 876 typedef value_type reference; 877 648 878 class pointer { 649 879 public: … … 654 884 }; 655 885 656 MappedInputIterator(const Base& _base, const LpSolverBase& _lp) 657 : base(_base), lp(_lp) {} 886 ExprIterator(const std::map<int, Value>::const_iterator& host_it, 887 const _solver_bits::VarIndex& index) 888 : _host_it(host_it), _index(index) {} 658 889 659 890 reference operator*() { 660 return std::make_pair( lp._lpId(base>first), base>second);891 return std::make_pair(_index(_host_it>first), _host_it>second); 661 892 } 662 893 … … 665 896 } 666 897 667 MappedInputIterator& operator++() { 668 ++base; 669 return *this; 670 } 671 672 MappedInputIterator operator++(int) { 673 MappedInputIterator tmp(*this); 674 ++base; 675 return tmp; 676 } 677 678 bool operator==(const MappedInputIterator& it) const { 679 return base == it.base; 680 } 681 682 bool operator!=(const MappedInputIterator& it) const { 683 return base != it.base; 684 } 685 686 private: 687 Base base; 688 const LpSolverBase& lp; 898 ExprIterator& operator++() { ++_host_it; return *this; } 899 ExprIterator operator++(int) { 900 ExprIterator tmp(*this); ++_host_it; return tmp; 901 } 902 903 ExprIterator& operator() { _host_it; return *this; } 904 ExprIterator operator(int) { 905 ExprIterator tmp(*this); _host_it; return tmp; 906 } 907 908 bool operator==(const ExprIterator& it) const { 909 return _host_it == it._host_it; 910 } 911 912 bool operator!=(const ExprIterator& it) const { 913 return _host_it != it._host_it; 914 } 915 689 916 }; 690 917 691 918 protected: 692 919 693 /// STL compatible iterator for lp col694 typedef MappedInputIterator<Expr> ConstRowIterator;695 /// STL compatible iterator for lp row696 typedef MappedInputIterator<DualExpr> ConstColIterator;697 698 /// STL compatible iterator for lp col699 typedef MappedOutputIterator<Expr> RowIterator;700 /// STL compatible iterator for lp row701 typedef MappedOutputIterator<DualExpr> ColIterator;702 703 920 //Abstract virtual functions 704 virtual LpSolverBase* _newLp() = 0; 705 virtual LpSolverBase* _copyLp(){ 706 LpSolverBase* newlp = _newLp(); 707 708 std::map<Col, Col> ref; 709 for (LpSolverBase::ColIt it(*this); it != INVALID; ++it) { 710 Col ccol = newlp>addCol(); 711 ref[it] = ccol; 712 newlp>colName(ccol, colName(it)); 713 newlp>colLowerBound(ccol, colLowerBound(it)); 714 newlp>colUpperBound(ccol, colUpperBound(it)); 715 } 716 717 for (LpSolverBase::RowIt it(*this); it != INVALID; ++it) { 718 Expr e = row(it), ce; 719 for (Expr::iterator jt = e.begin(); jt != e.end(); ++jt) { 720 ce[ref[jt>first]] = jt>second; 721 } 722 ce += e.constComp(); 723 Row r = newlp>addRow(ce); 724 725 double lower, upper; 726 getRowBounds(it, lower, upper); 727 newlp>rowBounds(r, lower, upper); 728 } 729 730 return newlp; 731 }; 921 virtual LpBase* _newSolver() const = 0; 922 virtual LpBase* _cloneSolver() const = 0; 923 924 virtual int _addColId(int col) { return cols.addIndex(col); } 925 virtual int _addRowId(int row) { return rows.addIndex(row); } 926 927 virtual void _eraseColId(int col) { cols.eraseIndex(col); } 928 virtual void _eraseRowId(int row) { rows.eraseIndex(row); } 732 929 733 930 virtual int _addCol() = 0; … … 737 934 virtual void _eraseRow(int row) = 0; 738 935 739 virtual void _getColName(int col, std::string 740 virtual void _setColName(int col, const std::string 936 virtual void _getColName(int col, std::string& name) const = 0; 937 virtual void _setColName(int col, const std::string& name) = 0; 741 938 virtual int _colByName(const std::string& name) const = 0; 742 939 743 virtual void _setRowCoeffs(int i, ConstRowIterator b, 744 ConstRowIterator e) = 0; 745 virtual void _getRowCoeffs(int i, RowIterator b) const = 0; 746 virtual void _setColCoeffs(int i, ConstColIterator b, 747 ConstColIterator e) = 0; 748 virtual void _getColCoeffs(int i, ColIterator b) const = 0; 940 virtual void _getRowName(int row, std::string& name) const = 0; 941 virtual void _setRowName(int row, const std::string& name) = 0; 942 virtual int _rowByName(const std::string& name) const = 0; 943 944 virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0; 945 virtual void _getRowCoeffs(int i, InsertIterator b) const = 0; 946 947 virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0; 948 virtual void _getColCoeffs(int i, InsertIterator b) const = 0; 949 749 950 virtual void _setCoeff(int row, int col, Value value) = 0; 750 951 virtual Value _getCoeff(int row, int col) const = 0; 952 751 953 virtual void _setColLowerBound(int i, Value value) = 0; 752 954 virtual Value _getColLowerBound(int i) const = 0; 955 753 956 virtual void _setColUpperBound(int i, Value value) = 0; 754 957 virtual Value _getColUpperBound(int i) const = 0; 755 virtual void _setRowBounds(int i, Value lower, Value upper) = 0; 756 virtual void _getRowBounds(int i, Value &lower, Value &upper) const = 0; 958 959 virtual void _setRowLowerBound(int i, Value value) = 0; 960 virtual Value _getRowLowerBound(int i) const = 0; 961 962 virtual void _setRowUpperBound(int i, Value value) = 0; 963 virtual Value _getRowUpperBound(int i) const = 0; 964 965 virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0; 966 virtual void _getObjCoeffs(InsertIterator b) const = 0; 757 967 758 968 virtual void _setObjCoeff(int i, Value obj_coef) = 0; 759 969 virtual Value _getObjCoeff(int i) const = 0; 760 virtual void _clearObj()=0; 761 762 virtual SolveExitStatus _solve() = 0; 763 virtual Value _getPrimal(int i) const = 0; 764 virtual Value _getDual(int i) const = 0; 765 virtual Value _getPrimalValue() const = 0; 766 virtual bool _isBasicCol(int i) const = 0; 767 virtual SolutionStatus _getPrimalStatus() const = 0; 768 virtual SolutionStatus _getDualStatus() const = 0; 769 virtual ProblemTypes _getProblemType() const = 0; 770 771 virtual void _setMax() = 0; 772 virtual void _setMin() = 0; 773 774 775 virtual bool _isMax() const = 0; 970 971 virtual void _setSense(Sense) = 0; 972 virtual Sense _getSense() const = 0; 973 974 virtual void _clear() = 0; 975 976 virtual const char* _solverName() const = 0; 776 977 777 978 //Own protected stuff … … 780 981 Value obj_const_comp; 781 982 983 LpBase() : rows(), cols(), obj_const_comp(0) {} 984 782 985 public: 783 986 784 ///\e 785 LpSolverBase() : obj_const_comp(0) {} 786 787 ///\e 788 virtual ~LpSolverBase() {} 987 /// Virtual destructor 988 virtual ~LpBase() {} 789 989 790 990 ///Creates a new LP problem 791 Lp SolverBase* newLp() {return _newLp();}991 LpBase* newSolver() {return _newSolver();} 792 992 ///Makes a copy of the LP problem 793 LpSolverBase* copyLp() {return _copyLp();} 993 LpBase* cloneSolver() {return _cloneSolver();} 994 995 ///Gives back the name of the solver. 996 const char* solverName() const {return _solverName();} 794 997 795 998 ///\name Build up and modify the LP … … 798 1001 799 1002 ///Add a new empty column (i.e a new variable) to the LP 800 Col addCol() { Col c; _addCol(); c.id = cols.addId(); return c;} 801 802 ///\brief Adds several new columns 803 ///(i.e a variables) at once 804 /// 805 ///This magic function takes a container as its argument 806 ///and fills its elements 807 ///with new columns (i.e. variables) 1003 Col addCol() { Col c; c._id = _addColId(_addCol()); return c;} 1004 1005 ///\brief Adds several new columns (i.e variables) at once 1006 /// 1007 ///This magic function takes a container as its argument and fills 1008 ///its elements with new columns (i.e. variables) 808 1009 ///\param t can be 809 1010 /// a standard STL compatible iterable container with 810 ///\ref Col as its \c values_type 811 ///like 1011 ///\ref Col as its \c values_type like 812 1012 ///\code 813 ///std::vector<Lp SolverBase::Col>814 ///std::list<Lp SolverBase::Col>1013 ///std::vector<LpBase::Col> 1014 ///std::list<LpBase::Col> 815 1015 ///\endcode 816 1016 /// a standard STL compatible iterable container with 817 ///\ref Col as its \c mapped_type 818 ///like 1017 ///\ref Col as its \c mapped_type like 819 1018 ///\code 820 ///std::map<AnyType,Lp SolverBase::Col>1019 ///std::map<AnyType,LpBase::Col> 821 1020 ///\endcode 822 1021 /// an iterable lemon \ref concepts::WriteMap "write map" like 823 1022 ///\code 824 ///ListGraph::NodeMap<Lp SolverBase::Col>825 ///ListGraph:: EdgeMap<LpSolverBase::Col>1023 ///ListGraph::NodeMap<LpBase::Col> 1024 ///ListGraph::ArcMap<LpBase::Col> 826 1025 ///\endcode 827 1026 ///\return The number of the created column. … … 831 1030 #else 832 1031 template<class T> 833 typename enable_if<typename T::value_type::Lp SolverCol,int>::type1032 typename enable_if<typename T::value_type::LpCol,int>::type 834 1033 addColSet(T &t,dummy<0> = 0) { 835 1034 int s=0; … … 838 1037 } 839 1038 template<class T> 840 typename enable_if<typename T::value_type::second_type::Lp SolverCol,1039 typename enable_if<typename T::value_type::second_type::LpCol, 841 1040 int>::type 842 1041 addColSet(T &t,dummy<1> = 1) { … … 849 1048 } 850 1049 template<class T> 851 typename enable_if<typename T::MapIt::Value::Lp SolverCol,1050 typename enable_if<typename T::MapIt::Value::LpCol, 852 1051 int>::type 853 1052 addColSet(T &t,dummy<2> = 2) { … … 867 1066 ///\param e is a dual linear expression (see \ref DualExpr) 868 1067 ///a better one. 869 void col(Col c, const DualExpr &e) {1068 void col(Col c, const DualExpr &e) { 870 1069 e.simplify(); 871 _setColCoeffs( _lpId(c), ConstColIterator(e.begin(), *this),872 ConstColIterator(e.end(), *this));1070 _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), cols), 1071 ExprIterator(e.comps.end(), cols)); 873 1072 } 874 1073 875 1074 ///Get a column (i.e a dual constraint) of the LP 876 1075 877 ///\param ris the column to get1076 ///\param c is the column to get 878 1077 ///\return the dual expression associated to the column 879 1078 DualExpr col(Col c) const { 880 1079 DualExpr e; 881 _getColCoeffs( _lpId(c), ColIterator(std::inserter(e, e.end()), *this));1080 _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows)); 882 1081 return e; 883 1082 } … … 886 1085 887 1086 ///\param e is a dual linear expression (see \ref DualExpr) 888 ///\param o bjis the corresponding component of the objective1087 ///\param o is the corresponding component of the objective 889 1088 ///function. It is 0 by default. 890 1089 ///\return The created column. … … 900 1099 ///This function adds a new empty row (i.e a new constraint) to the LP. 901 1100 ///\return The created row 902 Row addRow() { Row r; _addRow(); r.id = rows.addId(); return r;} 903 904 ///\brief Add several new rows 905 ///(i.e a constraints) at once 906 /// 907 ///This magic function takes a container as its argument 908 ///and fills its elements 909 ///with new row (i.e. variables) 1101 Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;} 1102 1103 ///\brief Add several new rows (i.e constraints) at once 1104 /// 1105 ///This magic function takes a container as its argument and fills 1106 ///its elements with new row (i.e. variables) 910 1107 ///\param t can be 911 1108 /// a standard STL compatible iterable container with 912 ///\ref Row as its \c values_type 913 ///like 1109 ///\ref Row as its \c values_type like 914 1110 ///\code 915 ///std::vector<Lp SolverBase::Row>916 ///std::list<Lp SolverBase::Row>1111 ///std::vector<LpBase::Row> 1112 ///std::list<LpBase::Row> 917 1113 ///\endcode 918 1114 /// a standard STL compatible iterable container with 919 ///\ref Row as its \c mapped_type 920 ///like 1115 ///\ref Row as its \c mapped_type like 921 1116 ///\code 922 ///std::map<AnyType,Lp SolverBase::Row>1117 ///std::map<AnyType,LpBase::Row> 923 1118 ///\endcode 924 1119 /// an iterable lemon \ref concepts::WriteMap "write map" like 925 1120 ///\code 926 ///ListGraph::NodeMap<Lp SolverBase::Row>927 ///ListGraph:: EdgeMap<LpSolverBase::Row>1121 ///ListGraph::NodeMap<LpBase::Row> 1122 ///ListGraph::ArcMap<LpBase::Row> 928 1123 ///\endcode 929 1124 ///\return The number of rows created. … … 933 1128 #else 934 1129 template<class T> 935 typename enable_if<typename T::value_type::Lp SolverRow,int>::type936 addRowSet(T &t, dummy<0> = 0) {1130 typename enable_if<typename T::value_type::LpRow,int>::type 1131 addRowSet(T &t, dummy<0> = 0) { 937 1132 int s=0; 938 1133 for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;} … … 940 1135 } 941 1136 template<class T> 942 typename enable_if<typename T::value_type::second_type::LpSolverRow, 943 int>::type 944 addRowSet(T &t,dummy<1> = 1) { 1137 typename enable_if<typename T::value_type::second_type::LpRow, int>::type 1138 addRowSet(T &t, dummy<1> = 1) { 945 1139 int s=0; 946 1140 for(typename T::iterator i=t.begin();i!=t.end();++i) { … … 951 1145 } 952 1146 template<class T> 953 typename enable_if<typename T::MapIt::Value::LpSolverRow, 954 int>::type 955 addRowSet(T &t,dummy<2> = 2) { 1147 typename enable_if<typename T::MapIt::Value::LpRow, int>::type 1148 addRowSet(T &t, dummy<2> = 2) { 956 1149 int s=0; 957 1150 for(typename T::MapIt i(t); i!=INVALID; ++i) … … 970 1163 ///\param e is a linear expression (see \ref Expr) 971 1164 ///\param u is the upper bound (\ref INF means no bound) 972 ///\bug This is a temporary function. The interface will change to973 ///a better one.974 ///\todo Option to control whether a constraint with a single variable is975 ///added or not.976 1165 void row(Row r, Value l, const Expr &e, Value u) { 977 1166 e.simplify(); 978 _setRowCoeffs(_lpId(r), ConstRowIterator(e.begin(), *this), 979 ConstRowIterator(e.end(), *this)); 980 _setRowBounds(_lpId(r),le.constComp(),ue.constComp()); 1167 _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols), 1168 ExprIterator(e.comps.end(), cols)); 1169 _setRowLowerBound(rows(id(r)),l  *e); 1170 _setRowUpperBound(rows(id(r)),u  *e); 981 1171 } 982 1172 … … 997 1187 Expr row(Row r) const { 998 1188 Expr e; 999 _getRowCoeffs( _lpId(r), RowIterator(std::inserter(e, e.end()), *this));1189 _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols)); 1000 1190 return e; 1001 1191 } … … 1007 1197 ///\param u is the upper bound (\ref INF means no bound) 1008 1198 ///\return The created row. 1009 ///\bug This is a temporary function. The interface will change to1010 ///a better one.1011 1199 Row addRow(Value l,const Expr &e, Value u) { 1012 1200 Row r=addRow(); … … 1024 1212 return r; 1025 1213 } 1026 ///Erase a coloumn (i.e a variable) from the LP 1027 1028 ///\param c is the coloumn to be deleted 1029 ///\todo Please check this 1030 void eraseCol(Col c) { 1031 _eraseCol(_lpId(c)); 1032 cols.eraseId(c.id); 1033 } 1034 ///Erase a row (i.e a constraint) from the LP 1214 ///Erase a column (i.e a variable) from the LP 1215 1216 ///\param c is the column to be deleted 1217 void erase(Col c) { 1218 _eraseCol(cols(id(c))); 1219 _eraseColId(cols(id(c))); 1220 } 1221 ///Erase a row (i.e a constraint) from the LP 1035 1222 1036 1223 ///\param r is the row to be deleted 1037 ///\todo Please check this 1038 void eraseRow(Row r) { 1039 _eraseRow(_lpId(r)); 1040 rows.eraseId(r.id); 1224 void erase(Row r) { 1225 _eraseRow(rows(id(r))); 1226 _eraseRowId(rows(id(r))); 1041 1227 } 1042 1228 1043 1229 /// Get the name of a column 1044 1230 1045 ///\param c is the coresponding col oumn1231 ///\param c is the coresponding column 1046 1232 ///\return The name of the colunm 1047 1233 std::string colName(Col c) const { 1048 1234 std::string name; 1049 _getColName( _lpId(c), name);1235 _getColName(cols(id(c)), name); 1050 1236 return name; 1051 1237 } … … 1053 1239 /// Set the name of a column 1054 1240 1055 ///\param c is the coresponding col oumn1241 ///\param c is the coresponding column 1056 1242 ///\param name The name to be given 1057 1243 void colName(Col c, const std::string& name) { 1058 _setColName( _lpId(c), name);1244 _setColName(cols(id(c)), name); 1059 1245 } 1060 1246 … … 1065 1251 Col colByName(const std::string& name) const { 1066 1252 int k = _colByName(name); 1067 return k != 1 ? Col(cols.fixId(k)) : Col(INVALID); 1253 return k != 1 ? Col(cols[k]) : Col(INVALID); 1254 } 1255 1256 /// Get the name of a row 1257 1258 ///\param r is the coresponding row 1259 ///\return The name of the row 1260 std::string rowName(Row r) const { 1261 std::string name; 1262 _getRowName(rows(id(r)), name); 1263 return name; 1264 } 1265 1266 /// Set the name of a row 1267 1268 ///\param r is the coresponding row 1269 ///\param name The name to be given 1270 void rowName(Row r, const std::string& name) { 1271 _setRowName(rows(id(r)), name); 1272 } 1273 1274 /// Get the row by its name 1275 1276 ///\param name The name of the row 1277 ///\return the proper row or \c INVALID 1278 Row rowByName(const std::string& name) const { 1279 int k = _rowByName(name); 1280 return k != 1 ? Row(rows[k]) : Row(INVALID); 1068 1281 } 1069 1282 … … 1071 1284 1072 1285 ///\param r is the row of the element to be modified 1073 ///\param c is the col oumn of the element to be modified1286 ///\param c is the column of the element to be modified 1074 1287 ///\param val is the new value of the coefficient 1075 1076 1288 void coeff(Row r, Col c, Value val) { 1077 _setCoeff( _lpId(r),_lpId(c), val);1289 _setCoeff(rows(id(r)),cols(id(c)), val); 1078 1290 } 1079 1291 1080 1292 /// Get an element of the coefficient matrix of the LP 1081 1293 1082 ///\param r is the row of the element in question1083 ///\param c is the col oumn of the element in question1294 ///\param r is the row of the element 1295 ///\param c is the column of the element 1084 1296 ///\return the corresponding coefficient 1085 1086 1297 Value coeff(Row r, Col c) const { 1087 return _getCoeff( _lpId(r),_lpId(c));1298 return _getCoeff(rows(id(r)),cols(id(c))); 1088 1299 } 1089 1300 … … 1094 1305 /// Value or \ref INF. 1095 1306 void colLowerBound(Col c, Value value) { 1096 _setColLowerBound( _lpId(c),value);1307 _setColLowerBound(cols(id(c)),value); 1097 1308 } 1098 1309 1099 1310 /// Get the lower bound of a column (i.e a variable) 1100 1311 1101 /// This function returns the lower bound for column (variable) \ tc1312 /// This function returns the lower bound for column (variable) \c c 1102 1313 /// (this might be \ref INF as well). 1103 ///\return The lower bound for col oumn \tc1314 ///\return The lower bound for column \c c 1104 1315 Value colLowerBound(Col c) const { 1105 return _getColLowerBound( _lpId(c));1316 return _getColLowerBound(cols(id(c))); 1106 1317 } 1107 1318 1108 1319 ///\brief Set the lower bound of several columns 1109 ///(i.e avariables) at once1320 ///(i.e variables) at once 1110 1321 /// 1111 1322 ///This magic function takes a container as its argument 1112 1323 ///and applies the function on all of its elements. 1113 /// 1114 /// 1115 /// 1324 ///The lower bound of a variable (column) has to be given by an 1325 ///extended number of type Value, i.e. a finite number of type 1326 ///Value or \ref INF. 1116 1327 #ifdef DOXYGEN 1117 1328 template<class T> … … 1119 1330 #else 1120 1331 template<class T> 1121 typename enable_if<typename T::value_type::Lp SolverCol,void>::type1332 typename enable_if<typename T::value_type::LpCol,void>::type 1122 1333 colLowerBound(T &t, Value value,dummy<0> = 0) { 1123 1334 for(typename T::iterator i=t.begin();i!=t.end();++i) { … … 1126 1337 } 1127 1338 template<class T> 1128 typename enable_if<typename T::value_type::second_type::Lp SolverCol,1339 typename enable_if<typename T::value_type::second_type::LpCol, 1129 1340 void>::type 1130 1341 colLowerBound(T &t, Value value,dummy<1> = 1) { … … 1134 1345 } 1135 1346 template<class T> 1136 typename enable_if<typename T::MapIt::Value::Lp SolverCol,1347 typename enable_if<typename T::MapIt::Value::LpCol, 1137 1348 void>::type 1138 1349 colLowerBound(T &t, Value value,dummy<2> = 2) { … … 1149 1360 /// Value or \ref INF. 1150 1361 void colUpperBound(Col c, Value value) { 1151 _setColUpperBound( _lpId(c),value);1362 _setColUpperBound(cols(id(c)),value); 1152 1363 }; 1153 1364 1154 1365 /// Get the upper bound of a column (i.e a variable) 1155 1366 1156 /// This function returns the upper bound for column (variable) \ tc1367 /// This function returns the upper bound for column (variable) \c c 1157 1368 /// (this might be \ref INF as well). 1158 /// \return The upper bound for coloumn \tc1369 /// \return The upper bound for column \c c 1159 1370 Value colUpperBound(Col c) const { 1160 return _getColUpperBound( _lpId(c));1371 return _getColUpperBound(cols(id(c))); 1161 1372 } 1162 1373 1163 1374 ///\brief Set the upper bound of several columns 1164 ///(i.e avariables) at once1375 ///(i.e variables) at once 1165 1376 /// 1166 1377 ///This magic function takes a container as its argument 1167 1378 ///and applies the function on all of its elements. 1168 /// 1169 /// 1170 /// 1379 ///The upper bound of a variable (column) has to be given by an 1380 ///extended number of type Value, i.e. a finite number of type 1381 ///Value or \ref INF. 1171 1382 #ifdef DOXYGEN 1172 1383 template<class T> … … 1174 1385 #else 1175 1386 template<class T> 1176 typename enable_if<typename T::value_type::Lp SolverCol,void>::type1387 typename enable_if<typename T::value_type::LpCol,void>::type 1177 1388 colUpperBound(T &t, Value value,dummy<0> = 0) { 1178 1389 for(typename T::iterator i=t.begin();i!=t.end();++i) { … … 1181 1392 } 1182 1393 template<class T> 1183 typename enable_if<typename T::value_type::second_type::Lp SolverCol,1394 typename enable_if<typename T::value_type::second_type::LpCol, 1184 1395 void>::type 1185 1396 colUpperBound(T &t, Value value,dummy<1> = 1) { … … 1189 1400 } 1190 1401 template<class T> 1191 typename enable_if<typename T::MapIt::Value::Lp SolverCol,1402 typename enable_if<typename T::MapIt::Value::LpCol, 1192 1403 void>::type 1193 1404 colUpperBound(T &t, Value value,dummy<2> = 2) { … … 1205 1416 /// Value, \ref INF or \ref INF. 1206 1417 void colBounds(Col c, Value lower, Value upper) { 1207 _setColLowerBound( _lpId(c),lower);1208 _setColUpperBound( _lpId(c),upper);1418 _setColLowerBound(cols(id(c)),lower); 1419 _setColUpperBound(cols(id(c)),upper); 1209 1420 } 1210 1421 1211 1422 ///\brief Set the lower and the upper bound of several columns 1212 ///(i.e avariables) at once1423 ///(i.e variables) at once 1213 1424 /// 1214 1425 ///This magic function takes a container as its argument … … 1223 1434 #else 1224 1435 template<class T> 1225 typename enable_if<typename T::value_type::Lp SolverCol,void>::type1436 typename enable_if<typename T::value_type::LpCol,void>::type 1226 1437 colBounds(T &t, Value lower, Value upper,dummy<0> = 0) { 1227 1438 for(typename T::iterator i=t.begin();i!=t.end();++i) { … … 1230 1441 } 1231 1442 template<class T> 1232 typename enable_if<typename T::value_type::second_type::LpSolverCol, 1233 void>::type 1443 typename enable_if<typename T::value_type::second_type::LpCol, void>::type 1234 1444 colBounds(T &t, Value lower, Value upper,dummy<1> = 1) { 1235 1445 for(typename T::iterator i=t.begin();i!=t.end();++i) { … … 1238 1448 } 1239 1449 template<class T> 1240 typename enable_if<typename T::MapIt::Value::LpSolverCol, 1241 void>::type 1450 typename enable_if<typename T::MapIt::Value::LpCol, void>::type 1242 1451 colBounds(T &t, Value lower, Value upper,dummy<2> = 2) { 1243 1452 for(typename T::MapIt i(t); i!=INVALID; ++i){ … … 1247 1456 #endif 1248 1457 1249 1250 /// Set the lower and the upper bounds of a row (i.e a constraint) 1251 1252 /// The lower and the upper bound of a constraint (row) have to be 1253 /// given by an extended number of type Value, i.e. a finite 1254 /// number of type Value, \ref INF or \ref INF. There is no 1255 /// separate function for the lower and the upper bound because 1256 /// that would have been hard to implement for CPLEX. 1257 void rowBounds(Row c, Value lower, Value upper) { 1258 _setRowBounds(_lpId(c),lower, upper); 1259 } 1260 1261 /// Get the lower and the upper bounds of a row (i.e a constraint) 1262 1263 /// The lower and the upper bound of 1264 /// a constraint (row) are 1265 /// extended numbers of type Value, i.e. finite numbers of type 1266 /// Value, \ref INF or \ref INF. 1267 /// \todo There is no separate function for the 1268 /// lower and the upper bound because we had problems with the 1269 /// implementation of the setting functions for CPLEX: 1270 /// check out whether this can be done for these functions. 1271 void getRowBounds(Row c, Value &lower, Value &upper) const { 1272 _getRowBounds(_lpId(c),lower, upper); 1458 /// Set the lower bound of a row (i.e a constraint) 1459 1460 /// The lower bound of a constraint (row) has to be given by an 1461 /// extended number of type Value, i.e. a finite number of type 1462 /// Value or \ref INF. 1463 void rowLowerBound(Row r, Value value) { 1464 _setRowLowerBound(rows(id(r)),value); 1465 } 1466 1467 /// Get the lower bound of a row (i.e a constraint) 1468 1469 /// This function returns the lower bound for row (constraint) \c c 1470 /// (this might be \ref INF as well). 1471 ///\return The lower bound for row \c r 1472 Value rowLowerBound(Row r) const { 1473 return _getRowLowerBound(rows(id(r))); 1474 } 1475 1476 /// Set the upper bound of a row (i.e a constraint) 1477 1478 /// The upper bound of a constraint (row) has to be given by an 1479 /// extended number of type Value, i.e. a finite number of type 1480 /// Value or \ref INF. 1481 void rowUpperBound(Row r, Value value) { 1482 _setRowUpperBound(rows(id(r)),value); 1483 } 1484 1485 /// Get the upper bound of a row (i.e a constraint) 1486 1487 /// This function returns the upper bound for row (constraint) \c c 1488 /// (this might be \ref INF as well). 1489 ///\return The upper bound for row \c r 1490 Value rowUpperBound(Row r) const { 1491 return _getRowUpperBound(rows(id(r))); 1273 1492 } 1274 1493 1275 1494 ///Set an element of the objective function 1276 void objCoeff(Col c, Value v) {_setObjCoeff( _lpId(c),v); };1495 void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); }; 1277 1496 1278 1497 ///Get an element of the objective function 1279 Value objCoeff(Col c) const { return _getObjCoeff( _lpId(c)); };1498 Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); }; 1280 1499 1281 1500 ///Set the objective function 1282 1501 1283 1502 ///\param e is a linear expression of type \ref Expr. 1284 void obj(Expr e) {1285 _clearObj();1286 for (Expr::iterator i=e.begin(); i!=e.end(); ++i)1287 objCoeff((*i).first,(*i).second);1288 obj_const_comp =e.constComp();1503 /// 1504 void obj(const Expr& e) { 1505 _setObjCoeffs(ExprIterator(e.comps.begin(), cols), 1506 ExprIterator(e.comps.end(), cols)); 1507 obj_const_comp = *e; 1289 1508 } 1290 1509 1291 1510 ///Get the objective function 1292 1511 1293 ///\return the objective function as a linear expression of type \ref Expr. 1512 ///\return the objective function as a linear expression of type 1513 ///Expr. 1294 1514 Expr obj() const { 1295 1515 Expr e; 1296 for (ColIt it(*this); it != INVALID; ++it) { 1297 double c = objCoeff(it); 1298 if (c != 0.0) { 1299 e.insert(std::make_pair(it, c)); 1300 } 1301 } 1516 _getObjCoeffs(InsertIterator(e.comps, cols)); 1517 *e = obj_const_comp; 1302 1518 return e; 1303 1519 } 1304 1520 1305 1521 1306 ///Maximize 1307 void max() { _setMax(); } 1308 ///Minimize 1309 void min() { _setMin(); } 1310 1311 ///Query function: is this a maximization problem? 1312 bool isMax() const {return _isMax(); } 1313 1314 ///Query function: is this a minimization problem? 1315 bool isMin() const {return !isMax(); } 1522 ///Set the direction of optimization 1523 void sense(Sense sense) { _setSense(sense); } 1524 1525 ///Query the direction of the optimization 1526 Sense sense() const {return _getSense(); } 1527 1528 ///Set the sense to maximization 1529 void max() { _setSense(MAX); } 1530 1531 ///Set the sense to maximization 1532 void min() { _setSense(MIN); } 1533 1534 ///Clears the problem 1535 void clear() { _clear(); } 1316 1536 1317 1537 ///@} 1318 1538 1539 }; 1540 1541 /// Addition 1542 1543 ///\relates LpBase::Expr 1544 /// 1545 inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) { 1546 LpBase::Expr tmp(a); 1547 tmp+=b; 1548 return tmp; 1549 } 1550 ///Substraction 1551 1552 ///\relates LpBase::Expr 1553 /// 1554 inline LpBase::Expr operator(const LpBase::Expr &a, const LpBase::Expr &b) { 1555 LpBase::Expr tmp(a); 1556 tmp=b; 1557 return tmp; 1558 } 1559 ///Multiply with constant 1560 1561 ///\relates LpBase::Expr 1562 /// 1563 inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) { 1564 LpBase::Expr tmp(a); 1565 tmp*=b; 1566 return tmp; 1567 } 1568 1569 ///Multiply with constant 1570 1571 ///\relates LpBase::Expr 1572 /// 1573 inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) { 1574 LpBase::Expr tmp(b); 1575 tmp*=a; 1576 return tmp; 1577 } 1578 ///Divide with constant 1579 1580 ///\relates LpBase::Expr 1581 /// 1582 inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) { 1583 LpBase::Expr tmp(a); 1584 tmp/=b; 1585 return tmp; 1586 } 1587 1588 ///Create constraint 1589 1590 ///\relates LpBase::Constr 1591 /// 1592 inline LpBase::Constr operator<=(const LpBase::Expr &e, 1593 const LpBase::Expr &f) { 1594 return LpBase::Constr(0, f  e, LpBase::INF); 1595 } 1596 1597 ///Create constraint 1598 1599 ///\relates LpBase::Constr 1600 /// 1601 inline LpBase::Constr operator<=(const LpBase::Value &e, 1602 const LpBase::Expr &f) { 1603 return LpBase::Constr(e, f, LpBase::NaN); 1604 } 1605 1606 ///Create constraint 1607 1608 ///\relates LpBase::Constr 1609 /// 1610 inline LpBase::Constr operator<=(const LpBase::Expr &e, 1611 const LpBase::Value &f) { 1612 return LpBase::Constr( LpBase::INF, e, f); 1613 } 1614 1615 ///Create constraint 1616 1617 ///\relates LpBase::Constr 1618 /// 1619 inline LpBase::Constr operator>=(const LpBase::Expr &e, 1620 const LpBase::Expr &f) { 1621 return LpBase::Constr(0, e  f, LpBase::INF); 1622 } 1623 1624 1625 ///Create constraint 1626 1627 ///\relates LpBase::Constr 1628 /// 1629 inline LpBase::Constr operator>=(const LpBase::Value &e, 1630 const LpBase::Expr &f) { 1631 return LpBase::Constr(LpBase::NaN, f, e); 1632 } 1633 1634 1635 ///Create constraint 1636 1637 ///\relates LpBase::Constr 1638 /// 1639 inline LpBase::Constr operator>=(const LpBase::Expr &e, 1640 const LpBase::Value &f) { 1641 return LpBase::Constr(f, e, LpBase::INF); 1642 } 1643 1644 ///Create constraint 1645 1646 ///\relates LpBase::Constr 1647 /// 1648 inline LpBase::Constr operator==(const LpBase::Expr &e, 1649 const LpBase::Value &f) { 1650 return LpBase::Constr(f, e, f); 1651 } 1652 1653 ///Create constraint 1654 1655 ///\relates LpBase::Constr 1656 /// 1657 inline LpBase::Constr operator==(const LpBase::Expr &e, 1658 const LpBase::Expr &f) { 1659 return LpBase::Constr(0, f  e, 0); 1660 } 1661 1662 ///Create constraint 1663 1664 ///\relates LpBase::Constr 1665 /// 1666 inline LpBase::Constr operator<=(const LpBase::Value &n, 1667 const LpBase::Constr &c) { 1668 LpBase::Constr tmp(c); 1669 LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint"); 1670 tmp.lowerBound()=n; 1671 return tmp; 1672 } 1673 ///Create constraint 1674 1675 ///\relates LpBase::Constr 1676 /// 1677 inline LpBase::Constr operator<=(const LpBase::Constr &c, 1678 const LpBase::Value &n) 1679 { 1680 LpBase::Constr tmp(c); 1681 LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint"); 1682 tmp.upperBound()=n; 1683 return tmp; 1684 } 1685 1686 ///Create constraint 1687 1688 ///\relates LpBase::Constr 1689 /// 1690 inline LpBase::Constr operator>=(const LpBase::Value &n, 1691 const LpBase::Constr &c) { 1692 LpBase::Constr tmp(c); 1693 LEMON_ASSERT(std::isnan(tmp.upperBound()), "Wrong LP constraint"); 1694 tmp.upperBound()=n; 1695 return tmp; 1696 } 1697 ///Create constraint 1698 1699 ///\relates LpBase::Constr 1700 /// 1701 inline LpBase::Constr operator>=(const LpBase::Constr &c, 1702 const LpBase::Value &n) 1703 { 1704 LpBase::Constr tmp(c); 1705 LEMON_ASSERT(std::isnan(tmp.lowerBound()), "Wrong LP constraint"); 1706 tmp.lowerBound()=n; 1707 return tmp; 1708 } 1709 1710 ///Addition 1711 1712 ///\relates LpBase::DualExpr 1713 /// 1714 inline LpBase::DualExpr operator+(const LpBase::DualExpr &a, 1715 const LpBase::DualExpr &b) { 1716 LpBase::DualExpr tmp(a); 1717 tmp+=b; 1718 return tmp; 1719 } 1720 ///Substraction 1721 1722 ///\relates LpBase::DualExpr 1723 /// 1724 inline LpBase::DualExpr operator(const LpBase::DualExpr &a, 1725 const LpBase::DualExpr &b) { 1726 LpBase::DualExpr tmp(a); 1727 tmp=b; 1728 return tmp; 1729 } 1730 ///Multiply with constant 1731 1732 ///\relates LpBase::DualExpr 1733 /// 1734 inline LpBase::DualExpr operator*(const LpBase::DualExpr &a, 1735 const LpBase::Value &b) { 1736 LpBase::DualExpr tmp(a); 1737 tmp*=b; 1738 return tmp; 1739 } 1740 1741 ///Multiply with constant 1742 1743 ///\relates LpBase::DualExpr 1744 /// 1745 inline LpBase::DualExpr operator*(const LpBase::Value &a, 1746 const LpBase::DualExpr &b) { 1747 LpBase::DualExpr tmp(b); 1748 tmp*=a; 1749 return tmp; 1750 } 1751 ///Divide with constant 1752 1753 ///\relates LpBase::DualExpr 1754 /// 1755 inline LpBase::DualExpr operator/(const LpBase::DualExpr &a, 1756 const LpBase::Value &b) { 1757 LpBase::DualExpr tmp(a); 1758 tmp/=b; 1759 return tmp; 1760 } 1761 1762 /// \ingroup lp_group 1763 /// 1764 /// \brief Common base class for LP solvers 1765 /// 1766 /// This class is an abstract base class for LP solvers. This class 1767 /// provides a full interface for set and modify an LP problem, 1768 /// solve it and retrieve the solution. You can use one of the 1769 /// descendants as a concrete implementation, or the \c Lp 1770 /// default LP solver. However, if you would like to handle LP 1771 /// solvers as reference or pointer in a generic way, you can use 1772 /// this class directly. 1773 class LpSolver : virtual public LpBase { 1774 public: 1775 1776 /// The problem types for primal and dual problems 1777 enum ProblemType { 1778 ///Feasible solution hasn't been found (but may exist). 1779 UNDEFINED = 0, 1780 ///The problem has no feasible solution 1781 INFEASIBLE = 1, 1782 ///Feasible solution found 1783 FEASIBLE = 2, 1784 ///Optimal solution exists and found 1785 OPTIMAL = 3, 1786 ///The cost function is unbounded 1787 UNBOUNDED = 4 1788 }; 1789 1790 ///The basis status of variables 1791 enum VarStatus { 1792 /// The variable is in the basis 1793 BASIC, 1794 /// The variable is free, but not basic 1795 FREE, 1796 /// The variable has active lower bound 1797 LOWER, 1798 /// The variable has active upper bound 1799 UPPER, 1800 /// The variable is nonbasic and fixed 1801 FIXED 1802 }; 1803 1804 protected: 1805 1806 virtual SolveExitStatus _solve() = 0; 1807 1808 virtual Value _getPrimal(int i) const = 0; 1809 virtual Value _getDual(int i) const = 0; 1810 1811 virtual Value _getPrimalRay(int i) const = 0; 1812 virtual Value _getDualRay(int i) const = 0; 1813 1814 virtual Value _getPrimalValue() const = 0; 1815 1816 virtual VarStatus _getColStatus(int i) const = 0; 1817 virtual VarStatus _getRowStatus(int i) const = 0; 1818 1819 virtual ProblemType _getPrimalType() const = 0; 1820 virtual ProblemType _getDualType() const = 0; 1821 1822 public: 1319 1823 1320 1824 ///\name Solve the LP … … 1327 1831 ///values and their meanings can be found in the documentation of 1328 1832 ///\ref SolveExitStatus. 1329 ///1330 ///\todo Which method is used to solve the problem1331 1833 SolveExitStatus solve() { return _solve(); } 1332 1834 … … 1337 1839 ///@{ 1338 1840 1339 /// The status of the primal problem (the original LP problem) 1340 SolutionStatus primalStatus() const { 1341 return _getPrimalStatus(); 1342 } 1343 1344 /// The status of the dual (of the original LP) problem 1345 SolutionStatus dualStatus() const { 1346 return _getDualStatus(); 1347 } 1348 1349 ///The type of the original LP problem 1350 ProblemTypes problemType() const { 1351 return _getProblemType(); 1352 } 1353 1354 ///\e 1355 Value primal(Col c) const { return _getPrimal(_lpId(c)); } 1356 ///\e 1841 /// The type of the primal problem 1842 ProblemType primalType() const { 1843 return _getPrimalType(); 1844 } 1845 1846 /// The type of the dual problem 1847 ProblemType dualType() const { 1848 return _getDualType(); 1849 } 1850 1851 /// Return the primal value of the column 1852 1853 /// Return the primal value of the column. 1854 /// \pre The problem is solved. 1855 Value primal(Col c) const { return _getPrimal(cols(id(c))); } 1856 1857 /// Return the primal value of the expression 1858 1859 /// Return the primal value of the expression, i.e. the dot 1860 /// product of the primal solution and the expression. 1861 /// \pre The problem is solved. 1357 1862 Value primal(const Expr& e) const { 1358 double res = e.constComp(); 1359 for (std::map<Col, double>::const_iterator it = e.begin(); 1360 it != e.end(); ++it) { 1361 res += _getPrimal(_lpId(it>first)) * it>second; 1863 double res = *e; 1864 for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) { 1865 res += *c * primal(c); 1362 1866 } 1363 1867 return res; 1364 1868 } 1365 1366 ///\e 1367 Value dual(Row r) const { return _getDual(_lpId(r)); } 1368 ///\e 1869 /// Returns a component of the primal ray 1870 1871 /// The primal ray is solution of the modified primal problem, 1872 /// where we change each finite bound to 0, and we looking for a 1873 /// negative objective value in case of minimization, and positive 1874 /// objective value for maximization. If there is such solution, 1875 /// that proofs the unsolvability of the dual problem, and if a 1876 /// feasible primal solution exists, then the unboundness of 1877 /// primal problem. 1878 /// 1879 /// \pre The problem is solved and the dual problem is infeasible. 1880 /// \note Some solvers does not provide primal ray calculation 1881 /// functions. 1882 Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); } 1883 1884 /// Return the dual value of the row 1885 1886 /// Return the dual value of the row. 1887 /// \pre The problem is solved. 1888 Value dual(Row r) const { return _getDual(rows(id(r))); } 1889 1890 /// Return the dual value of the dual expression 1891 1892 /// Return the dual value of the dual expression, i.e. the dot 1893 /// product of the dual solution and the dual expression. 1894 /// \pre The problem is solved. 1369 1895 Value dual(const DualExpr& e) const { 1370 1896 double res = 0.0; 1371 for (std::map<Row, double>::const_iterator it = e.begin(); 1372 it != e.end(); ++it) { 1373 res += _getPrimal(_lpId(it>first)) * it>second; 1897 for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) { 1898 res += *r * dual(r); 1374 1899 } 1375 1900 return res; 1376 1901 } 1377 1902 1378 ///\e 1379 bool isBasicCol(Col c) const { return _isBasicCol(_lpId(c)); } 1380 1381 ///\e 1903 /// Returns a component of the dual ray 1904 1905 /// The dual ray is solution of the modified primal problem, where 1906 /// we change each finite bound to 0 (i.e. the objective function 1907 /// coefficients in the primal problem), and we looking for a 1908 /// ositive objective value. If there is such solution, that 1909 /// proofs the unsolvability of the primal problem, and if a 1910 /// feasible dual solution exists, then the unboundness of 1911 /// dual problem. 1912 /// 1913 /// \pre The problem is solved and the primal problem is infeasible. 1914 /// \note Some solvers does not provide dual ray calculation 1915 /// functions. 1916 Value dualRay(Row r) const { return _getDualRay(rows(id(r))); } 1917 1918 /// Return the basis status of the column 1919 1920 /// \see VarStatus 1921 VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); } 1922 1923 /// Return the basis status of the row 1924 1925 /// \see VarStatus 1926 VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); } 1927 1928 ///The value of the objective function 1382 1929 1383 1930 ///\return … … 1386 1933 /// \ref NaN if no primal solution is found. 1387 1934 /// The (finite) objective value if an optimal solution is found. 1388 Value primal Value() const { return _getPrimalValue()+obj_const_comp;}1935 Value primal() const { return _getPrimalValue()+obj_const_comp;} 1389 1936 ///@} 1390 1937 1938 LpSolver* newSolver() {return _newSolver();} 1939 LpSolver* cloneSolver() {return _cloneSolver();} 1940 1941 protected: 1942 1943 virtual LpSolver* _newSolver() const = 0; 1944 virtual LpSolver* _cloneSolver() const = 0; 1391 1945 }; 1392 1946 … … 1395 1949 /// 1396 1950 /// \brief Common base class for MIP solvers 1397 /// \todo Much more docs 1398 class MipSolverBase : virtual public LpSolverBase{ 1951 /// 1952 /// This class is an abstract base class for MIP solvers. This class 1953 /// provides a full interface for set and modify an MIP problem, 1954 /// solve it and retrieve the solution. You can use one of the 1955 /// descendants as a concrete implementation, or the \c Lp 1956 /// default MIP solver. However, if you would like to handle MIP 1957 /// solvers as reference or pointer in a generic way, you can use 1958 /// this class directly. 1959 class MipSolver : virtual public LpBase { 1399 1960 public: 1400 1961 1401 ///Possible variable (coloumn) types (e.g. real, integer, binary etc.) 1962 /// The problem types for MIP problems 1963 enum ProblemType { 1964 ///Feasible solution hasn't been found (but may exist). 1965 UNDEFINED = 0, 1966 ///The problem has no feasible solution 1967 INFEASIBLE = 1, 1968 ///Feasible solution found 1969 FEASIBLE = 2, 1970 ///Optimal solution exists and found 1971 OPTIMAL = 3, 1972 ///The cost function is unbounded 1973 /// 1974 ///The Mip or at least the relaxed problem is unbounded 1975 UNBOUNDED = 4 1976 }; 1977 1978 ///\name Solve the MIP 1979 1980 ///@{ 1981 1982 /// Solve the MIP problem at hand 1983 /// 1984 ///\return The result of the optimization procedure. Possible 1985 ///values and their meanings can be found in the documentation of 1986 ///\ref SolveExitStatus. 1987 SolveExitStatus solve() { return _solve(); } 1988 1989 ///@} 1990 1991 ///\name Setting column type 1992 ///@{ 1993 1994 ///Possible variable (column) types (e.g. real, integer, binary etc.) 1402 1995 enum ColTypes { 1403 ///Continuous variable 1996 ///Continuous variable (default) 1404 1997 REAL = 0, 1405 1998 ///Integer variable 1406 1407 ///Unfortunately, cplex 7.5 somewhere writes something like 1408 ///#define INTEGER 'I' 1409 INT = 1 1410 ///\todo No support for other types yet. 1999 INTEGER = 1 1411 2000 }; 1412 2001 1413 ///Sets the type of the given coloumn to the given type 1414 /// 1415 ///Sets the type of the given coloumn to the given type. 2002 ///Sets the type of the given column to the given type 2003 2004 ///Sets the type of the given column to the given type. 2005 /// 1416 2006 void colType(Col c, ColTypes col_type) { 1417 _ colType(_lpId(c),col_type);2007 _setColType(cols(id(c)),col_type); 1418 2008 } 1419 2009 1420 2010 ///Gives back the type of the column. 1421 /// 2011 1422 2012 ///Gives back the type of the column. 2013 /// 1423 2014 ColTypes colType(Col c) const { 1424 return _colType(_lpId(c)); 1425 } 1426 1427 ///Sets the type of the given Col to integer or remove that property. 1428 /// 1429 ///Sets the type of the given Col to integer or remove that property. 1430 void integer(Col c, bool enable) { 1431 if (enable) 1432 colType(c,INT); 1433 else 1434 colType(c,REAL); 1435 } 1436 1437 ///Gives back whether the type of the column is integer or not. 1438 /// 1439 ///Gives back the type of the column. 1440 ///\return true if the column has integer type and false if not. 1441 bool integer(Col c) const { 1442 return (colType(c)==INT); 1443 } 1444 1445 /// The status of the MIP problem 1446 SolutionStatus mipStatus() const { 1447 return _getMipStatus(); 1448 } 2015 return _getColType(cols(id(c))); 2016 } 2017 ///@} 2018 2019 ///\name Obtain the solution 2020 2021 ///@{ 2022 2023 /// The type of the MIP problem 2024 ProblemType type() const { 2025 return _getType(); 2026 } 2027 2028 /// Return the value of the row in the solution 2029 2030 /// Return the value of the row in the solution. 2031 /// \pre The problem is solved. 2032 Value sol(Col c) const { return _getSol(cols(id(c))); } 2033 2034 /// Return the value of the expression in the solution 2035 2036 /// Return the value of the expression in the solution, i.e. the 2037 /// dot product of the solution and the expression. 2038 /// \pre The problem is solved. 2039 Value sol(const Expr& e) const { 2040 double res = *e; 2041 for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) { 2042 res += *c * sol(c); 2043 } 2044 return res; 2045 } 2046 ///The value of the objective function 2047 2048 ///\return 2049 /// \ref INF or \ref INF means either infeasibility or unboundedness 2050 /// of the problem, depending on whether we minimize or maximize. 2051 /// \ref NaN if no primal solution is found. 2052 /// The (finite) objective value if an optimal solution is found. 2053 Value solValue() const { return _getSolValue()+obj_const_comp;} 2054 ///@} 1449 2055 1450 2056 protected: 1451 2057 1452 virtual ColTypes _colType(int col) const = 0; 1453 virtual void _colType(int col, ColTypes col_type) = 0; 1454 virtual SolutionStatus _getMipStatus() const = 0; 1455 2058 virtual SolveExitStatus _solve() = 0; 2059 virtual ColTypes _getColType(int col) const = 0; 2060 virtual void _setColType(int col, ColTypes col_type) = 0; 2061 virtual ProblemType _getType() const = 0; 2062 virtual Value _getSol(int i) const = 0; 2063 virtual Value _getSolValue() const = 0; 2064 2065 public: 2066 2067 MipSolver* newSolver() {return _newSolver();} 2068 MipSolver* cloneSolver() {return _cloneSolver();} 2069 2070 protected: 2071 2072 virtual MipSolver* _newSolver() const = 0; 2073 virtual MipSolver* _cloneSolver() const = 0; 1456 2074 }; 1457 2075 1458 ///\relates LpSolverBase::Expr1459 ///1460 inline LpSolverBase::Expr operator+(const LpSolverBase::Expr &a,1461 const LpSolverBase::Expr &b)1462 {1463 LpSolverBase::Expr tmp(a);1464 tmp+=b;1465 return tmp;1466 }1467 ///\e1468 1469 ///\relates LpSolverBase::Expr1470 ///1471 inline LpSolverBase::Expr operator(const LpSolverBase::Expr &a,1472 const LpSolverBase::Expr &b)1473 {1474 LpSolverBase::Expr tmp(a);1475 tmp=b;1476 return tmp;1477 }1478 ///\e1479 1480 ///\relates LpSolverBase::Expr1481 ///1482 inline LpSolverBase::Expr operator*(const LpSolverBase::Expr &a,1483 const LpSolverBase::Value &b)1484 {1485 LpSolverBase::Expr tmp(a);1486 tmp*=b;1487 return tmp;1488 }1489 1490 ///\e1491 1492 ///\relates LpSolverBase::Expr1493 ///1494 inline LpSolverBase::Expr operator*(const LpSolverBase::Value &a,1495 const LpSolverBase::Expr &b)1496 {1497 LpSolverBase::Expr tmp(b);1498 tmp*=a;1499 return tmp;1500 }1501 ///\e1502 1503 ///\relates LpSolverBase::Expr1504 ///1505 inline LpSolverBase::Expr operator/(const LpSolverBase::Expr &a,1506 const LpSolverBase::Value &b)1507 {1508 LpSolverBase::Expr tmp(a);1509 tmp/=b;1510 return tmp;1511 }1512 1513 ///\e1514 1515 ///\relates LpSolverBase::Constr1516 ///1517 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,1518 const LpSolverBase::Expr &f)1519 {1520 return LpSolverBase::Constr(LpSolverBase::INF,ef,0);1521 }1522 1523 ///\e1524 1525 ///\relates LpSolverBase::Constr1526 ///1527 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &e,1528 const LpSolverBase::Expr &f)1529 {1530 return LpSolverBase::Constr(e,f);1531 }1532 1533 ///\e1534 1535 ///\relates LpSolverBase::Constr1536 ///1537 inline LpSolverBase::Constr operator<=(const LpSolverBase::Expr &e,1538 const LpSolverBase::Value &f)1539 {1540 return LpSolverBase::Constr(LpSolverBase::INF,e,f);1541 }1542 1543 ///\e1544 1545 ///\relates LpSolverBase::Constr1546 ///1547 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,1548 const LpSolverBase::Expr &f)1549 {1550 return LpSolverBase::Constr(LpSolverBase::INF,fe,0);1551 }1552 1553 1554 ///\e1555 1556 ///\relates LpSolverBase::Constr1557 ///1558 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &e,1559 const LpSolverBase::Expr &f)1560 {1561 return LpSolverBase::Constr(f,e);1562 }1563 1564 1565 ///\e1566 1567 ///\relates LpSolverBase::Constr1568 ///1569 inline LpSolverBase::Constr operator>=(const LpSolverBase::Expr &e,1570 const LpSolverBase::Value &f)1571 {1572 return LpSolverBase::Constr(f,e,LpSolverBase::INF);1573 }1574 1575 ///\e1576 1577 ///\relates LpSolverBase::Constr1578 ///1579 inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,1580 const LpSolverBase::Value &f)1581 {1582 return LpSolverBase::Constr(f,e,f);1583 }1584 1585 ///\e1586 1587 ///\relates LpSolverBase::Constr1588 ///1589 inline LpSolverBase::Constr operator==(const LpSolverBase::Expr &e,1590 const LpSolverBase::Expr &f)1591 {1592 return LpSolverBase::Constr(0,ef,0);1593 }1594 1595 ///\e1596 1597 ///\relates LpSolverBase::Constr1598 ///1599 inline LpSolverBase::Constr operator<=(const LpSolverBase::Value &n,1600 const LpSolverBase::Constr&c)1601 {1602 LpSolverBase::Constr tmp(c);1603 LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");1604 tmp.lowerBound()=n;1605 return tmp;1606 }1607 ///\e1608 1609 ///\relates LpSolverBase::Constr1610 ///1611 inline LpSolverBase::Constr operator<=(const LpSolverBase::Constr& c,1612 const LpSolverBase::Value &n)1613 {1614 LpSolverBase::Constr tmp(c);1615 LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");1616 tmp.upperBound()=n;1617 return tmp;1618 }1619 1620 ///\e1621 1622 ///\relates LpSolverBase::Constr1623 ///1624 inline LpSolverBase::Constr operator>=(const LpSolverBase::Value &n,1625 const LpSolverBase::Constr&c)1626 {1627 LpSolverBase::Constr tmp(c);1628 LEMON_ASSERT(LpSolverBase::isNaN(tmp.upperBound()), "Wrong LP constraint");1629 tmp.upperBound()=n;1630 return tmp;1631 }1632 ///\e1633 1634 ///\relates LpSolverBase::Constr1635 ///1636 inline LpSolverBase::Constr operator>=(const LpSolverBase::Constr& c,1637 const LpSolverBase::Value &n)1638 {1639 LpSolverBase::Constr tmp(c);1640 LEMON_ASSERT(LpSolverBase::isNaN(tmp.lowerBound()), "Wrong LP constraint");1641 tmp.lowerBound()=n;1642 return tmp;1643 }1644 1645 ///\e1646 1647 ///\relates LpSolverBase::DualExpr1648 ///1649 inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,1650 const LpSolverBase::DualExpr &b)1651 {1652 LpSolverBase::DualExpr tmp(a);1653 tmp+=b;1654 return tmp;1655 }1656 ///\e1657 1658 ///\relates LpSolverBase::DualExpr1659 ///1660 inline LpSolverBase::DualExpr operator(const LpSolverBase::DualExpr &a,1661 const LpSolverBase::DualExpr &b)1662 {1663 LpSolverBase::DualExpr tmp(a);1664 tmp=b;1665 return tmp;1666 }1667 ///\e1668 1669 ///\relates LpSolverBase::DualExpr1670 ///1671 inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,1672 const LpSolverBase::Value &b)1673 {1674 LpSolverBase::DualExpr tmp(a);1675 tmp*=b;1676 return tmp;1677 }1678 1679 ///\e1680 1681 ///\relates LpSolverBase::DualExpr1682 ///1683 inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,1684 const LpSolverBase::DualExpr &b)1685 {1686 LpSolverBase::DualExpr tmp(b);1687 tmp*=a;1688 return tmp;1689 }1690 ///\e1691 1692 ///\relates LpSolverBase::DualExpr1693 ///1694 inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,1695 const LpSolverBase::Value &b)1696 {1697 LpSolverBase::DualExpr tmp(a);1698 tmp/=b;1699 return tmp;1700 }1701 2076 1702 2077
Note: See TracChangeset
for help on using the changeset viewer.