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