lemon/clp.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Sat, 20 Feb 2010 18:39:03 +0100
changeset 839 f3bc4e9b5f3a
parent 576 745e182d0139
child 877 141f9c0db4a3
permissions -rw-r--r--
New heuristics for MCF algorithms (#340)
and some implementation improvements.

- A useful heuristic is added to NetworkSimplex to make the
initial pivots faster.
- A powerful global update heuristic is added to CostScaling
and the implementation is reworked with various improvements.
- Better relabeling in CostScaling to improve numerical stability
and make the code faster.
- A small improvement is made in CapacityScaling for better
delta computation.
- Add notes to the classes about the usage of vector<char> instead
of vector<bool> for efficiency reasons.
     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-2008
     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 #include <lemon/clp.h>
    20 #include <coin/ClpSimplex.hpp>
    21 
    22 namespace lemon {
    23 
    24   ClpLp::ClpLp() {
    25     _prob = new ClpSimplex();
    26     _init_temporals();
    27     messageLevel(MESSAGE_NOTHING);
    28   }
    29 
    30   ClpLp::ClpLp(const ClpLp& other) {
    31     _prob = new ClpSimplex(*other._prob);
    32     rows = other.rows;
    33     cols = other.cols;
    34     _init_temporals();
    35     messageLevel(MESSAGE_NOTHING);
    36   }
    37 
    38   ClpLp::~ClpLp() {
    39     delete _prob;
    40     _clear_temporals();
    41   }
    42 
    43   void ClpLp::_init_temporals() {
    44     _primal_ray = 0;
    45     _dual_ray = 0;
    46   }
    47 
    48   void ClpLp::_clear_temporals() {
    49     if (_primal_ray) {
    50       delete[] _primal_ray;
    51       _primal_ray = 0;
    52     }
    53     if (_dual_ray) {
    54       delete[] _dual_ray;
    55       _dual_ray = 0;
    56     }
    57   }
    58 
    59   ClpLp* ClpLp::newSolver() const {
    60     ClpLp* newlp = new ClpLp;
    61     return newlp;
    62   }
    63 
    64   ClpLp* ClpLp::cloneSolver() const {
    65     ClpLp* copylp = new ClpLp(*this);
    66     return copylp;
    67   }
    68 
    69   const char* ClpLp::_solverName() const { return "ClpLp"; }
    70 
    71   int ClpLp::_addCol() {
    72     _prob->addColumn(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX, 0.0);
    73     return _prob->numberColumns() - 1;
    74   }
    75 
    76   int ClpLp::_addRow() {
    77     _prob->addRow(0, 0, 0, -COIN_DBL_MAX, COIN_DBL_MAX);
    78     return _prob->numberRows() - 1;
    79   }
    80 
    81   int ClpLp::_addRow(Value l, ExprIterator b, ExprIterator e, Value u) {
    82     std::vector<int> indexes;
    83     std::vector<Value> values;
    84 
    85     for(ExprIterator it = b; it != e; ++it) {
    86       indexes.push_back(it->first);
    87       values.push_back(it->second);
    88     }
    89 
    90     _prob->addRow(values.size(), &indexes.front(), &values.front(), l, u);
    91     return _prob->numberRows() - 1;
    92   }
    93 
    94 
    95   void ClpLp::_eraseCol(int c) {
    96     _col_names_ref.erase(_prob->getColumnName(c));
    97     _prob->deleteColumns(1, &c);
    98   }
    99 
   100   void ClpLp::_eraseRow(int r) {
   101     _row_names_ref.erase(_prob->getRowName(r));
   102     _prob->deleteRows(1, &r);
   103   }
   104 
   105   void ClpLp::_eraseColId(int i) {
   106     cols.eraseIndex(i);
   107     cols.shiftIndices(i);
   108   }
   109 
   110   void ClpLp::_eraseRowId(int i) {
   111     rows.eraseIndex(i);
   112     rows.shiftIndices(i);
   113   }
   114 
   115   void ClpLp::_getColName(int c, std::string& name) const {
   116     name = _prob->getColumnName(c);
   117   }
   118 
   119   void ClpLp::_setColName(int c, const std::string& name) {
   120     _prob->setColumnName(c, const_cast<std::string&>(name));
   121     _col_names_ref[name] = c;
   122   }
   123 
   124   int ClpLp::_colByName(const std::string& name) const {
   125     std::map<std::string, int>::const_iterator it = _col_names_ref.find(name);
   126     return it != _col_names_ref.end() ? it->second : -1;
   127   }
   128 
   129   void ClpLp::_getRowName(int r, std::string& name) const {
   130     name = _prob->getRowName(r);
   131   }
   132 
   133   void ClpLp::_setRowName(int r, const std::string& name) {
   134     _prob->setRowName(r, const_cast<std::string&>(name));
   135     _row_names_ref[name] = r;
   136   }
   137 
   138   int ClpLp::_rowByName(const std::string& name) const {
   139     std::map<std::string, int>::const_iterator it = _row_names_ref.find(name);
   140     return it != _row_names_ref.end() ? it->second : -1;
   141   }
   142 
   143 
   144   void ClpLp::_setRowCoeffs(int ix, ExprIterator b, ExprIterator e) {
   145     std::map<int, Value> coeffs;
   146 
   147     int n = _prob->clpMatrix()->getNumCols();
   148 
   149     const int* indices = _prob->clpMatrix()->getIndices();
   150     const double* elements = _prob->clpMatrix()->getElements();
   151 
   152     for (int i = 0; i < n; ++i) {
   153       CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
   154       CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
   155 
   156       const int* it = std::lower_bound(indices + begin, indices + end, ix);
   157       if (it != indices + end && *it == ix && elements[it - indices] != 0.0) {
   158         coeffs[i] = 0.0;
   159       }
   160     }
   161 
   162     for (ExprIterator it = b; it != e; ++it) {
   163       coeffs[it->first] = it->second;
   164     }
   165 
   166     for (std::map<int, Value>::iterator it = coeffs.begin();
   167          it != coeffs.end(); ++it) {
   168       _prob->modifyCoefficient(ix, it->first, it->second);
   169     }
   170   }
   171 
   172   void ClpLp::_getRowCoeffs(int ix, InsertIterator b) const {
   173     int n = _prob->clpMatrix()->getNumCols();
   174 
   175     const int* indices = _prob->clpMatrix()->getIndices();
   176     const double* elements = _prob->clpMatrix()->getElements();
   177 
   178     for (int i = 0; i < n; ++i) {
   179       CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
   180       CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
   181 
   182       const int* it = std::lower_bound(indices + begin, indices + end, ix);
   183       if (it != indices + end && *it == ix) {
   184         *b = std::make_pair(i, elements[it - indices]);
   185       }
   186     }
   187   }
   188 
   189   void ClpLp::_setColCoeffs(int ix, ExprIterator b, ExprIterator e) {
   190     std::map<int, Value> coeffs;
   191 
   192     CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
   193     CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
   194 
   195     const int* indices = _prob->clpMatrix()->getIndices();
   196     const double* elements = _prob->clpMatrix()->getElements();
   197 
   198     for (CoinBigIndex i = begin; i != end; ++i) {
   199       if (elements[i] != 0.0) {
   200         coeffs[indices[i]] = 0.0;
   201       }
   202     }
   203     for (ExprIterator it = b; it != e; ++it) {
   204       coeffs[it->first] = it->second;
   205     }
   206     for (std::map<int, Value>::iterator it = coeffs.begin();
   207          it != coeffs.end(); ++it) {
   208       _prob->modifyCoefficient(it->first, ix, it->second);
   209     }
   210   }
   211 
   212   void ClpLp::_getColCoeffs(int ix, InsertIterator b) const {
   213     CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
   214     CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
   215 
   216     const int* indices = _prob->clpMatrix()->getIndices();
   217     const double* elements = _prob->clpMatrix()->getElements();
   218 
   219     for (CoinBigIndex i = begin; i != end; ++i) {
   220       *b = std::make_pair(indices[i], elements[i]);
   221       ++b;
   222     }
   223   }
   224 
   225   void ClpLp::_setCoeff(int ix, int jx, Value value) {
   226     _prob->modifyCoefficient(ix, jx, value);
   227   }
   228 
   229   ClpLp::Value ClpLp::_getCoeff(int ix, int jx) const {
   230     CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
   231     CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
   232 
   233     const int* indices = _prob->clpMatrix()->getIndices();
   234     const double* elements = _prob->clpMatrix()->getElements();
   235 
   236     const int* it = std::lower_bound(indices + begin, indices + end, jx);
   237     if (it != indices + end && *it == jx) {
   238       return elements[it - indices];
   239     } else {
   240       return 0.0;
   241     }
   242   }
   243 
   244   void ClpLp::_setColLowerBound(int i, Value lo) {
   245     _prob->setColumnLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
   246   }
   247 
   248   ClpLp::Value ClpLp::_getColLowerBound(int i) const {
   249     double val = _prob->getColLower()[i];
   250     return val == - COIN_DBL_MAX ? - INF : val;
   251   }
   252 
   253   void ClpLp::_setColUpperBound(int i, Value up) {
   254     _prob->setColumnUpper(i, up == INF ? COIN_DBL_MAX : up);
   255   }
   256 
   257   ClpLp::Value ClpLp::_getColUpperBound(int i) const {
   258     double val = _prob->getColUpper()[i];
   259     return val == COIN_DBL_MAX ? INF : val;
   260   }
   261 
   262   void ClpLp::_setRowLowerBound(int i, Value lo) {
   263     _prob->setRowLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
   264   }
   265 
   266   ClpLp::Value ClpLp::_getRowLowerBound(int i) const {
   267     double val = _prob->getRowLower()[i];
   268     return val == - COIN_DBL_MAX ? - INF : val;
   269   }
   270 
   271   void ClpLp::_setRowUpperBound(int i, Value up) {
   272     _prob->setRowUpper(i, up == INF ? COIN_DBL_MAX : up);
   273   }
   274 
   275   ClpLp::Value ClpLp::_getRowUpperBound(int i) const {
   276     double val = _prob->getRowUpper()[i];
   277     return val == COIN_DBL_MAX ? INF : val;
   278   }
   279 
   280   void ClpLp::_setObjCoeffs(ExprIterator b, ExprIterator e) {
   281     int num = _prob->clpMatrix()->getNumCols();
   282     for (int i = 0; i < num; ++i) {
   283       _prob->setObjectiveCoefficient(i, 0.0);
   284     }
   285     for (ExprIterator it = b; it != e; ++it) {
   286       _prob->setObjectiveCoefficient(it->first, it->second);
   287     }
   288   }
   289 
   290   void ClpLp::_getObjCoeffs(InsertIterator b) const {
   291     int num = _prob->clpMatrix()->getNumCols();
   292     for (int i = 0; i < num; ++i) {
   293       Value coef = _prob->getObjCoefficients()[i];
   294       if (coef != 0.0) {
   295         *b = std::make_pair(i, coef);
   296         ++b;
   297       }
   298     }
   299   }
   300 
   301   void ClpLp::_setObjCoeff(int i, Value obj_coef) {
   302     _prob->setObjectiveCoefficient(i, obj_coef);
   303   }
   304 
   305   ClpLp::Value ClpLp::_getObjCoeff(int i) const {
   306     return _prob->getObjCoefficients()[i];
   307   }
   308 
   309   ClpLp::SolveExitStatus ClpLp::_solve() {
   310     return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
   311   }
   312 
   313   ClpLp::SolveExitStatus ClpLp::solvePrimal() {
   314     return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
   315   }
   316 
   317   ClpLp::SolveExitStatus ClpLp::solveDual() {
   318     return _prob->dual() >= 0 ? SOLVED : UNSOLVED;
   319   }
   320 
   321   ClpLp::SolveExitStatus ClpLp::solveBarrier() {
   322     return _prob->barrier() >= 0 ? SOLVED : UNSOLVED;
   323   }
   324 
   325   ClpLp::Value ClpLp::_getPrimal(int i) const {
   326     return _prob->primalColumnSolution()[i];
   327   }
   328   ClpLp::Value ClpLp::_getPrimalValue() const {
   329     return _prob->objectiveValue();
   330   }
   331 
   332   ClpLp::Value ClpLp::_getDual(int i) const {
   333     return _prob->dualRowSolution()[i];
   334   }
   335 
   336   ClpLp::Value ClpLp::_getPrimalRay(int i) const {
   337     if (!_primal_ray) {
   338       _primal_ray = _prob->unboundedRay();
   339       LEMON_ASSERT(_primal_ray != 0, "Primal ray is not provided");
   340     }
   341     return _primal_ray[i];
   342   }
   343 
   344   ClpLp::Value ClpLp::_getDualRay(int i) const {
   345     if (!_dual_ray) {
   346       _dual_ray = _prob->infeasibilityRay();
   347       LEMON_ASSERT(_dual_ray != 0, "Dual ray is not provided");
   348     }
   349     return _dual_ray[i];
   350   }
   351 
   352   ClpLp::VarStatus ClpLp::_getColStatus(int i) const {
   353     switch (_prob->getColumnStatus(i)) {
   354     case ClpSimplex::basic:
   355       return BASIC;
   356     case ClpSimplex::isFree:
   357       return FREE;
   358     case ClpSimplex::atUpperBound:
   359       return UPPER;
   360     case ClpSimplex::atLowerBound:
   361       return LOWER;
   362     case ClpSimplex::isFixed:
   363       return FIXED;
   364     case ClpSimplex::superBasic:
   365       return FREE;
   366     default:
   367       LEMON_ASSERT(false, "Wrong column status");
   368       return VarStatus();
   369     }
   370   }
   371 
   372   ClpLp::VarStatus ClpLp::_getRowStatus(int i) const {
   373     switch (_prob->getColumnStatus(i)) {
   374     case ClpSimplex::basic:
   375       return BASIC;
   376     case ClpSimplex::isFree:
   377       return FREE;
   378     case ClpSimplex::atUpperBound:
   379       return UPPER;
   380     case ClpSimplex::atLowerBound:
   381       return LOWER;
   382     case ClpSimplex::isFixed:
   383       return FIXED;
   384     case ClpSimplex::superBasic:
   385       return FREE;
   386     default:
   387       LEMON_ASSERT(false, "Wrong row status");
   388       return VarStatus();
   389     }
   390   }
   391 
   392 
   393   ClpLp::ProblemType ClpLp::_getPrimalType() const {
   394     if (_prob->isProvenOptimal()) {
   395       return OPTIMAL;
   396     } else if (_prob->isProvenPrimalInfeasible()) {
   397       return INFEASIBLE;
   398     } else if (_prob->isProvenDualInfeasible()) {
   399       return UNBOUNDED;
   400     } else {
   401       return UNDEFINED;
   402     }
   403   }
   404 
   405   ClpLp::ProblemType ClpLp::_getDualType() const {
   406     if (_prob->isProvenOptimal()) {
   407       return OPTIMAL;
   408     } else if (_prob->isProvenDualInfeasible()) {
   409       return INFEASIBLE;
   410     } else if (_prob->isProvenPrimalInfeasible()) {
   411       return INFEASIBLE;
   412     } else {
   413       return UNDEFINED;
   414     }
   415   }
   416 
   417   void ClpLp::_setSense(ClpLp::Sense sense) {
   418     switch (sense) {
   419     case MIN:
   420       _prob->setOptimizationDirection(1);
   421       break;
   422     case MAX:
   423       _prob->setOptimizationDirection(-1);
   424       break;
   425     }
   426   }
   427 
   428   ClpLp::Sense ClpLp::_getSense() const {
   429     double dir = _prob->optimizationDirection();
   430     if (dir > 0.0) {
   431       return MIN;
   432     } else {
   433       return MAX;
   434     }
   435   }
   436 
   437   void ClpLp::_clear() {
   438     delete _prob;
   439     _prob = new ClpSimplex();
   440     rows.clear();
   441     cols.clear();
   442     _col_names_ref.clear();
   443     _clear_temporals();
   444   }
   445 
   446   void ClpLp::_messageLevel(MessageLevel level) {
   447     switch (level) {
   448     case MESSAGE_NOTHING:
   449       _prob->setLogLevel(0);
   450       break;
   451     case MESSAGE_ERROR:
   452       _prob->setLogLevel(1);
   453       break;
   454     case MESSAGE_WARNING:
   455       _prob->setLogLevel(2);
   456       break;
   457     case MESSAGE_NORMAL:
   458       _prob->setLogLevel(3);
   459       break;
   460     case MESSAGE_VERBOSE:
   461       _prob->setLogLevel(4);
   462       break;
   463     }
   464   }
   465 
   466 } //END OF NAMESPACE LEMON