lemon/clp.cc
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 03 Apr 2009 18:59:15 +0200
changeset 600 6ac5d9ae1d3d
parent 461 08d495d48089
child 528 9db62975c32b
permissions -rw-r--r--
Support real types + numerical stability fix in NS (#254)

- Real types are supported by appropriate inicialization.
- A feature of the XTI spanning tree structure is removed to ensure
numerical stability (could cause problems using integer types).
The node potentials are updated always on the lower subtree,
in order to prevent overflow problems.
The former method isn't notably faster during to our tests.
     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_NO_OUTPUT);
    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_NO_OUTPUT);
    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 
    82   void ClpLp::_eraseCol(int c) {
    83     _col_names_ref.erase(_prob->getColumnName(c));
    84     _prob->deleteColumns(1, &c);
    85   }
    86 
    87   void ClpLp::_eraseRow(int r) {
    88     _row_names_ref.erase(_prob->getRowName(r));
    89     _prob->deleteRows(1, &r);
    90   }
    91 
    92   void ClpLp::_eraseColId(int i) {
    93     cols.eraseIndex(i);
    94     cols.shiftIndices(i);
    95   }
    96 
    97   void ClpLp::_eraseRowId(int i) {
    98     rows.eraseIndex(i);
    99     rows.shiftIndices(i);
   100   }
   101 
   102   void ClpLp::_getColName(int c, std::string& name) const {
   103     name = _prob->getColumnName(c);
   104   }
   105 
   106   void ClpLp::_setColName(int c, const std::string& name) {
   107     _prob->setColumnName(c, const_cast<std::string&>(name));
   108     _col_names_ref[name] = c;
   109   }
   110 
   111   int ClpLp::_colByName(const std::string& name) const {
   112     std::map<std::string, int>::const_iterator it = _col_names_ref.find(name);
   113     return it != _col_names_ref.end() ? it->second : -1;
   114   }
   115 
   116   void ClpLp::_getRowName(int r, std::string& name) const {
   117     name = _prob->getRowName(r);
   118   }
   119 
   120   void ClpLp::_setRowName(int r, const std::string& name) {
   121     _prob->setRowName(r, const_cast<std::string&>(name));
   122     _row_names_ref[name] = r;
   123   }
   124 
   125   int ClpLp::_rowByName(const std::string& name) const {
   126     std::map<std::string, int>::const_iterator it = _row_names_ref.find(name);
   127     return it != _row_names_ref.end() ? it->second : -1;
   128   }
   129 
   130 
   131   void ClpLp::_setRowCoeffs(int ix, ExprIterator b, ExprIterator e) {
   132     std::map<int, Value> coeffs;
   133 
   134     int n = _prob->clpMatrix()->getNumCols();
   135 
   136     const int* indices = _prob->clpMatrix()->getIndices();
   137     const double* elements = _prob->clpMatrix()->getElements();
   138 
   139     for (int i = 0; i < n; ++i) {
   140       CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
   141       CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
   142 
   143       const int* it = std::lower_bound(indices + begin, indices + end, ix);
   144       if (it != indices + end && *it == ix && elements[it - indices] != 0.0) {
   145         coeffs[i] = 0.0;
   146       }
   147     }
   148 
   149     for (ExprIterator it = b; it != e; ++it) {
   150       coeffs[it->first] = it->second;
   151     }
   152 
   153     for (std::map<int, Value>::iterator it = coeffs.begin();
   154          it != coeffs.end(); ++it) {
   155       _prob->modifyCoefficient(ix, it->first, it->second);
   156     }
   157   }
   158 
   159   void ClpLp::_getRowCoeffs(int ix, InsertIterator b) const {
   160     int n = _prob->clpMatrix()->getNumCols();
   161 
   162     const int* indices = _prob->clpMatrix()->getIndices();
   163     const double* elements = _prob->clpMatrix()->getElements();
   164 
   165     for (int i = 0; i < n; ++i) {
   166       CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[i];
   167       CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[i];
   168 
   169       const int* it = std::lower_bound(indices + begin, indices + end, ix);
   170       if (it != indices + end && *it == ix) {
   171         *b = std::make_pair(i, elements[it - indices]);
   172       }
   173     }
   174   }
   175 
   176   void ClpLp::_setColCoeffs(int ix, ExprIterator b, ExprIterator e) {
   177     std::map<int, Value> coeffs;
   178 
   179     CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
   180     CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
   181 
   182     const int* indices = _prob->clpMatrix()->getIndices();
   183     const double* elements = _prob->clpMatrix()->getElements();
   184 
   185     for (CoinBigIndex i = begin; i != end; ++i) {
   186       if (elements[i] != 0.0) {
   187         coeffs[indices[i]] = 0.0;
   188       }
   189     }
   190     for (ExprIterator it = b; it != e; ++it) {
   191       coeffs[it->first] = it->second;
   192     }
   193     for (std::map<int, Value>::iterator it = coeffs.begin();
   194          it != coeffs.end(); ++it) {
   195       _prob->modifyCoefficient(it->first, ix, it->second);
   196     }
   197   }
   198 
   199   void ClpLp::_getColCoeffs(int ix, InsertIterator b) const {
   200     CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
   201     CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
   202 
   203     const int* indices = _prob->clpMatrix()->getIndices();
   204     const double* elements = _prob->clpMatrix()->getElements();
   205 
   206     for (CoinBigIndex i = begin; i != end; ++i) {
   207       *b = std::make_pair(indices[i], elements[i]);
   208       ++b;
   209     }
   210   }
   211 
   212   void ClpLp::_setCoeff(int ix, int jx, Value value) {
   213     _prob->modifyCoefficient(ix, jx, value);
   214   }
   215 
   216   ClpLp::Value ClpLp::_getCoeff(int ix, int jx) const {
   217     CoinBigIndex begin = _prob->clpMatrix()->getVectorStarts()[ix];
   218     CoinBigIndex end = begin + _prob->clpMatrix()->getVectorLengths()[ix];
   219 
   220     const int* indices = _prob->clpMatrix()->getIndices();
   221     const double* elements = _prob->clpMatrix()->getElements();
   222 
   223     const int* it = std::lower_bound(indices + begin, indices + end, jx);
   224     if (it != indices + end && *it == jx) {
   225       return elements[it - indices];
   226     } else {
   227       return 0.0;
   228     }
   229   }
   230 
   231   void ClpLp::_setColLowerBound(int i, Value lo) {
   232     _prob->setColumnLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
   233   }
   234 
   235   ClpLp::Value ClpLp::_getColLowerBound(int i) const {
   236     double val = _prob->getColLower()[i];
   237     return val == - COIN_DBL_MAX ? - INF : val;
   238   }
   239 
   240   void ClpLp::_setColUpperBound(int i, Value up) {
   241     _prob->setColumnUpper(i, up == INF ? COIN_DBL_MAX : up);
   242   }
   243 
   244   ClpLp::Value ClpLp::_getColUpperBound(int i) const {
   245     double val = _prob->getColUpper()[i];
   246     return val == COIN_DBL_MAX ? INF : val;
   247   }
   248 
   249   void ClpLp::_setRowLowerBound(int i, Value lo) {
   250     _prob->setRowLower(i, lo == - INF ? - COIN_DBL_MAX : lo);
   251   }
   252 
   253   ClpLp::Value ClpLp::_getRowLowerBound(int i) const {
   254     double val = _prob->getRowLower()[i];
   255     return val == - COIN_DBL_MAX ? - INF : val;
   256   }
   257 
   258   void ClpLp::_setRowUpperBound(int i, Value up) {
   259     _prob->setRowUpper(i, up == INF ? COIN_DBL_MAX : up);
   260   }
   261 
   262   ClpLp::Value ClpLp::_getRowUpperBound(int i) const {
   263     double val = _prob->getRowUpper()[i];
   264     return val == COIN_DBL_MAX ? INF : val;
   265   }
   266 
   267   void ClpLp::_setObjCoeffs(ExprIterator b, ExprIterator e) {
   268     int num = _prob->clpMatrix()->getNumCols();
   269     for (int i = 0; i < num; ++i) {
   270       _prob->setObjectiveCoefficient(i, 0.0);
   271     }
   272     for (ExprIterator it = b; it != e; ++it) {
   273       _prob->setObjectiveCoefficient(it->first, it->second);
   274     }
   275   }
   276 
   277   void ClpLp::_getObjCoeffs(InsertIterator b) const {
   278     int num = _prob->clpMatrix()->getNumCols();
   279     for (int i = 0; i < num; ++i) {
   280       Value coef = _prob->getObjCoefficients()[i];
   281       if (coef != 0.0) {
   282         *b = std::make_pair(i, coef);
   283         ++b;
   284       }
   285     }
   286   }
   287 
   288   void ClpLp::_setObjCoeff(int i, Value obj_coef) {
   289     _prob->setObjectiveCoefficient(i, obj_coef);
   290   }
   291 
   292   ClpLp::Value ClpLp::_getObjCoeff(int i) const {
   293     return _prob->getObjCoefficients()[i];
   294   }
   295 
   296   ClpLp::SolveExitStatus ClpLp::_solve() {
   297     return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
   298   }
   299 
   300   ClpLp::SolveExitStatus ClpLp::solvePrimal() {
   301     return _prob->primal() >= 0 ? SOLVED : UNSOLVED;
   302   }
   303 
   304   ClpLp::SolveExitStatus ClpLp::solveDual() {
   305     return _prob->dual() >= 0 ? SOLVED : UNSOLVED;
   306   }
   307 
   308   ClpLp::SolveExitStatus ClpLp::solveBarrier() {
   309     return _prob->barrier() >= 0 ? SOLVED : UNSOLVED;
   310   }
   311 
   312   ClpLp::Value ClpLp::_getPrimal(int i) const {
   313     return _prob->primalColumnSolution()[i];
   314   }
   315   ClpLp::Value ClpLp::_getPrimalValue() const {
   316     return _prob->objectiveValue();
   317   }
   318 
   319   ClpLp::Value ClpLp::_getDual(int i) const {
   320     return _prob->dualRowSolution()[i];
   321   }
   322 
   323   ClpLp::Value ClpLp::_getPrimalRay(int i) const {
   324     if (!_primal_ray) {
   325       _primal_ray = _prob->unboundedRay();
   326       LEMON_ASSERT(_primal_ray != 0, "Primal ray is not provided");
   327     }
   328     return _primal_ray[i];
   329   }
   330 
   331   ClpLp::Value ClpLp::_getDualRay(int i) const {
   332     if (!_dual_ray) {
   333       _dual_ray = _prob->infeasibilityRay();
   334       LEMON_ASSERT(_dual_ray != 0, "Dual ray is not provided");
   335     }
   336     return _dual_ray[i];
   337   }
   338 
   339   ClpLp::VarStatus ClpLp::_getColStatus(int i) const {
   340     switch (_prob->getColumnStatus(i)) {
   341     case ClpSimplex::basic:
   342       return BASIC;
   343     case ClpSimplex::isFree:
   344       return FREE;
   345     case ClpSimplex::atUpperBound:
   346       return UPPER;
   347     case ClpSimplex::atLowerBound:
   348       return LOWER;
   349     case ClpSimplex::isFixed:
   350       return FIXED;
   351     case ClpSimplex::superBasic:
   352       return FREE;
   353     default:
   354       LEMON_ASSERT(false, "Wrong column status");
   355       return VarStatus();
   356     }
   357   }
   358 
   359   ClpLp::VarStatus ClpLp::_getRowStatus(int i) const {
   360     switch (_prob->getColumnStatus(i)) {
   361     case ClpSimplex::basic:
   362       return BASIC;
   363     case ClpSimplex::isFree:
   364       return FREE;
   365     case ClpSimplex::atUpperBound:
   366       return UPPER;
   367     case ClpSimplex::atLowerBound:
   368       return LOWER;
   369     case ClpSimplex::isFixed:
   370       return FIXED;
   371     case ClpSimplex::superBasic:
   372       return FREE;
   373     default:
   374       LEMON_ASSERT(false, "Wrong row status");
   375       return VarStatus();
   376     }
   377   }
   378 
   379 
   380   ClpLp::ProblemType ClpLp::_getPrimalType() const {
   381     if (_prob->isProvenOptimal()) {
   382       return OPTIMAL;
   383     } else if (_prob->isProvenPrimalInfeasible()) {
   384       return INFEASIBLE;
   385     } else if (_prob->isProvenDualInfeasible()) {
   386       return UNBOUNDED;
   387     } else {
   388       return UNDEFINED;
   389     }
   390   }
   391 
   392   ClpLp::ProblemType ClpLp::_getDualType() const {
   393     if (_prob->isProvenOptimal()) {
   394       return OPTIMAL;
   395     } else if (_prob->isProvenDualInfeasible()) {
   396       return INFEASIBLE;
   397     } else if (_prob->isProvenPrimalInfeasible()) {
   398       return INFEASIBLE;
   399     } else {
   400       return UNDEFINED;
   401     }
   402   }
   403 
   404   void ClpLp::_setSense(ClpLp::Sense sense) {
   405     switch (sense) {
   406     case MIN:
   407       _prob->setOptimizationDirection(1);
   408       break;
   409     case MAX:
   410       _prob->setOptimizationDirection(-1);
   411       break;
   412     }
   413   }
   414 
   415   ClpLp::Sense ClpLp::_getSense() const {
   416     double dir = _prob->optimizationDirection();
   417     if (dir > 0.0) {
   418       return MIN;
   419     } else {
   420       return MAX;
   421     }
   422   }
   423 
   424   void ClpLp::_clear() {
   425     delete _prob;
   426     _prob = new ClpSimplex();
   427     rows.clear();
   428     cols.clear();
   429     _col_names_ref.clear();
   430     _clear_temporals();
   431   }
   432 
   433   void ClpLp::messageLevel(MessageLevel m) {
   434     _prob->setLogLevel(static_cast<int>(m));
   435   }
   436 
   437 } //END OF NAMESPACE LEMON