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