CBCSolver.cc 9.39 KB
Newer Older
1 2
// (C) Copyright 2015 by Autodesk, Inc.

3 4 5 6 7 8 9 10 11 12 13 14 15 16
//=============================================================================
//
//  CLASS CBCSolver - IMPLEMENTATION
//
//=============================================================================

//== INCLUDES =================================================================

//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_CBC_AVAILABLE

//=============================================================================
#include "CBCSolver.hh"
17
#include <CoMISo/Utils/CoMISoError.hh>
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

#include <Base/Debug/DebTime.hh>

// For Branch and bound
#include "OsiSolverInterface.hpp"
#include "CbcModel.hpp"
#include "CbcCutGenerator.hpp"
#include "CbcStrategy.hpp"
#include "CbcHeuristicLocal.hpp"
#include "OsiClpSolverInterface.hpp"

// Cuts
#include <CbcModel.hpp>
#include <OsiClpSolverInterface.hpp>

#include "CglGomory.hpp"
#include "CglProbing.hpp"
#include "CglKnapsackCover.hpp"
#include "CglRedSplit.hpp"
#include "CglClique.hpp"
#include "CglFlowCover.hpp"
#include "CglMixedIntegerRounding2.hpp"

// Preprocessing
#include "CglPreProcess.hpp"

// Heuristics
#include "CbcHeuristic.hpp"
46
#include "CbcCompareDepth.hpp"
47 48 49 50 51 52 53 54 55 56 57

#include <stdexcept>

#define CBC_INFINITY COIN_DBL_MAX

//== NAMESPACES ===============================================================

namespace COMISO {

//== IMPLEMENTATION ==========================================================

58 59
namespace {

60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
// These are some "tweaks" for the CbcModel, copied from some sample Cbc code
// I have no idea what any of these do!
void model_tweak(CbcModel& model)
{
  CglProbing generator1;
  generator1.setUsingObjective(true);
  generator1.setMaxPass(1);
  generator1.setMaxPassRoot(5);
  // Number of unsatisfied variables to look at
  generator1.setMaxProbe(10);
  generator1.setMaxProbeRoot(1000);
  // How far to follow the consequences
  generator1.setMaxLook(50);
  generator1.setMaxLookRoot(500);
  // Only look at rows with fewer than this number of elements
  generator1.setMaxElements(200);
  generator1.setRowCuts(3);

  CglGomory generator2;
  // try larger limit
  generator2.setLimit(300);

  CglKnapsackCover generator3;

  CglRedSplit generator4;
  // try larger limit
  generator4.setLimit(200);

  CglClique generator5;
  generator5.setStarCliqueReport(false);
  generator5.setRowCliqueReport(false);

  CglMixedIntegerRounding2 mixedGen;
  CglFlowCover flowGen;

  // Add in generators
  // Experiment with -1 and -99 etc
  //model.addCutGenerator(&generator1, -1, "Probing");
  model.addCutGenerator(&generator2, -1, "Gomory");
  //model.addCutGenerator(&generator3, -1, "Knapsack");
  //model.addCutGenerator(&generator4,-1,"RedSplit");
  //model.addCutGenerator(&generator5, -1, "Clique");
  //model.addCutGenerator(&flowGen, -1, "FlowCover");
  //model.addCutGenerator(&mixedGen, -1, "MixedIntegerRounding");
  // Say we want timings
  int numberGenerators = model.numberCutGenerators();
  int iGenerator;
  for (iGenerator = 0; iGenerator < numberGenerators; iGenerator++)
  {
    CbcCutGenerator* generator = model.cutGenerator(iGenerator);
    generator->setTiming(true);
  }
  //auto osiclp = dynamic_cast<OsiClpSolverInterface*>(model.solver());
  //// go faster stripes
  //if (osiclp)
  //{
  //  // Turn this off if you get problems
  //  // Used to be automatically set
  //  osiclp->setSpecialOptions(128);
  //  if (osiclp->getNumRows() < 300 && osiclp->getNumCols() < 500)
  //  {
  //    //osiclp->setupForRepeatedUse(2,0);
  //    osiclp->setupForRepeatedUse(0, 0);
  //  }
  //}
  // Uncommenting this should switch off all CBC messages
  // model.messagesPointer()->setDetailMessages(10,10000,NULL);
  // Allow rounding heuristic

  CbcRounding heuristic1(model);
  model.addHeuristic(&heuristic1);

  // And local search when new solution found
  //CbcHeuristicLocal heuristic2(model);
  //model.addHeuristic(&heuristic2);
}

#define TRACE_CBT(DESCR, EXPR) DEB_line(7, DESCR << ": " << EXPR)
#define P(X) ((X).data())

140
bool solve_impl(
141 142 143
  NProblemInterface*                        _problem,
  const std::vector<NConstraintInterface*>& _constraints,
  const std::vector<PairIndexVtype>&        _discrete_constraints,
144
  const double                              _time_limit
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
)
{
  DEB_enter_func;
  DEB_warning_if(!_problem->constant_hessian(), 1,
                 "CBCSolver received a problem with non-constant hessian!");
  DEB_warning_if(!_problem->constant_gradient(), 1,
                 "CBCSolver received a problem with non-constant gradient!");

  const int n_rows = _constraints.size(); // Constraints #
  const int n_cols = _problem->n_unknowns(); // Unknowns #

  // expand the variable types from discrete mtrx array
  std::vector<VariableType> var_type(n_cols, Real);
  for (size_t i = 0; i < _discrete_constraints.size(); ++i)
    var_type[_discrete_constraints[i].first] = _discrete_constraints[i].second;

  //----------------------------------------------
  // set up constraints
  //----------------------------------------------
  std::vector<double> x(n_cols, 0.0); // solution

  CoinPackedMatrix mtrx(false, 0, 0);// make a row-ordered, empty mtrx
  mtrx.setDimensions(0, n_cols);

  std::vector<double> row_lb(n_rows, -CBC_INFINITY);
  std::vector<double> row_ub(n_rows,  CBC_INFINITY);

  for (size_t i = 0; i < _constraints.size();  ++i)
  {
    DEB_error_if(!_constraints[i]->is_linear(), "constraint " << i <<
                 " is non-linear");

    NConstraintInterface::SVectorNC gc;
    _constraints[i]->eval_gradient(P(x), gc);

    // the setup below appears inefficient, the interface is rather restrictive
    CoinPackedVector lin_expr;
    lin_expr.reserve(gc.size());
    for (NConstraintInterface::SVectorNC::InnerIterator v_it(gc); v_it; ++v_it)
      lin_expr.insert(v_it.index(), v_it.value());
    mtrx.appendRow(lin_expr);

    const double b = -_constraints[i]->eval_constraint(P(x));
    switch (_constraints[i]->constraint_type())
    {
    case NConstraintInterface::NC_EQUAL         :
      row_lb[i] = row_ub[i] = b;
      break;
    case NConstraintInterface::NC_LESS_EQUAL    :
      row_ub[i] = b;
      break;
    case NConstraintInterface::NC_GREATER_EQUAL :
      row_lb[i] = b;
      break;
    }

    TRACE_CBT("Constraint " << i << " is of type " <<
              _constraints[i]->constraint_type() << "; RHS val", b);
  }

  //----------------------------------------------
  // setup energy
  //----------------------------------------------

  std::vector<double> objective(n_cols);
  _problem->eval_gradient(P(x), P(objective));
  TRACE_CBT("Objective linear term", objective);

  const double c = _problem->eval_f(P(x));
  // ICGB: Where does objective constant term go in CBC?
  // MCM: I could not find this either: It is possible that the entire model
  // can be translated to accomodate the constant (if != 0). Ask DB!
  DEB_error_if(c > FLT_EPSILON, "Ignoring a non-zero constant objective term: "
               << c);
  TRACE_CBT("Objective constant term", c);

  // CBC Problem initialize
  OsiClpSolverInterface si;
  si.setHintParam(OsiDoReducePrint, true, OsiHintTry);

  const std::vector<double> col_lb(n_cols, -CBC_INFINITY);
  const std::vector<double> col_ub(n_cols,  CBC_INFINITY);
  si.loadProblem(mtrx, P(col_lb), P(col_ub), P(objective), P(row_lb), P(row_ub));

  // set variables
  for (int i = 0; i < n_cols; ++i)
  {
    switch (var_type[i])
    {
    case Real:
      si.setContinuous(i);
      break;
    case Integer:
      si.setInteger(i);
      break;
    case Binary :
      si.setInteger(i);
      si.setColBounds(i, 0.0, 1.0);
      break;
    }
  }
  si.initialSolve();

  // TODO: make this accessible through the DEB system instead
  volatile static bool dump_problem = false; // change on run-time if necessary
  if (dump_problem)
251 252
    si.writeMps("CBC_problem_dump"); //output problem as .MPS

253 254 255
  // Pass the OsiSolver with the problem to be solved to CbcModel
  CbcModel model(si);
  model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
256 257 258 259
//  model.setMaximumSolutions(4);
//  TRACE_CBT("CbcModel::getMaximumSolutions()", model.getMaximumSolutions());
  model.setMaximumSeconds(_time_limit);
//  TRACE_CBT("CbcModel::getMaximumSeconds()", model.getMaximumSeconds());
260 261

  model_tweak(model);
262 263 264 265

  CbcCompareDepth compare_depth;
  model.setNodeComparison(&compare_depth);

266 267
  model.branchAndBound();

268 269 270 271 272 273 274 275 276 277 278
  const double* solution = model.bestSolution();
  //COMISO_THROW_if(solution == nullptr, UNSPECIFIED_CBC_EXCEPTION);

  // store if solution available
  if(solution != 0)
  {
    _problem->store_result(solution);
    return true;
  }
  else
    return false;
279 280 281 282
}

#undef P

283 284
}//namespace 

285
bool CBCSolver::solve(
286 287 288 289 290 291 292
  NProblemInterface*                        _problem,
  const std::vector<NConstraintInterface*>& _constraints,
  const std::vector<PairIndexVtype>&        _discrete_constraints,
  const double                              _time_limit
)
{
  DEB_enter_func;
293
  bool valid_solution = false;
294 295
  try
  {
296
    valid_solution = solve_impl(_problem, _constraints, _discrete_constraints, _time_limit);
297 298 299 300
  }
  catch (CoinError& ce)
  {
    DEB_warning(1, "CoinError code = " << ce.message() << "]\n");
301
    //COMISO_THROW(UNSPECIFIED_CBC_EXCEPTION);
302
  }
303
  return valid_solution;
304 305 306 307 308 309 310 311 312 313
}


//=============================================================================
} // namespace COMISO
//=============================================================================

#endif // COMISO_CBC_AVAILABLE
//=============================================================================