source: trunk/Bonmin/src/Algorithms/OaGenerators/BonOaDecBase.cpp @ 1536

Last change on this file since 1536 was 1536, checked in by pbonami, 10 years ago

Port Claudia's code to trunk

File size: 15.8 KB
Line 
1// (C) Copyright International Business Machines (IBM) 2006
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// P. Bonami, International Business Machines
7//
8// Date :  12/07/2006
9
10#include <sstream>
11#include <climits>
12
13#include <algorithm>
14#include "BonOaDecBase.hpp"
15
16
17#include "BonminConfig.h"
18
19#include "OsiClpSolverInterface.hpp"
20
21#include "CbcModel.hpp"
22#include "CbcStrategy.hpp"
23#ifdef COIN_HAS_CPX
24#include "OsiCpxSolverInterface.hpp"
25
26#define CHECK_CPX_STAT(a,b) if(b) throw CoinError("Error in CPLEX call",__FILE__,a);
27
28#endif
29#include "BonCbc.hpp"
30#include "BonSolverHelp.hpp"
31//The following two are to interupt the solution of sub-mip through CTRL-C
32extern CbcModel * OAModel;
33
34namespace Bonmin {
35
36  OaDecompositionBase::OaDecompositionBase
37  (OsiTMINLPInterface * nlp)
38      :
39      CglCutGenerator(),
40      nlp_(nlp),
41      s_(NULL),
42      nSolve_(0),
43      lp_(NULL),
44      objects_(NULL),
45      nObjects_(0),
46      nLocalSearch_(0),
47      handler_(NULL),
48      leaveSiUnchanged_(true),
49      reassignLpsolver_(false),
50      timeBegin_(0),
51      parameters_(),
52      currentNodeNumber_(-1)
53  {
54    handler_ = new CoinMessageHandler();
55    handler_ -> setLogLevel(2);
56    messages_ = OaMessages();
57    timeBegin_ = CoinCpuTime();
58  }
59
60  OaDecompositionBase::OaDecompositionBase(BabSetupBase &b, bool leaveSiUnchanged,
61      bool reassignLpsolver):
62      CglCutGenerator(),
63      nlp_(b.nonlinearSolver()),
64      lp_(NULL),
65      objects_(NULL),
66      nObjects_(0),
67      nLocalSearch_(0),
68      handler_(NULL),
69      leaveSiUnchanged_(leaveSiUnchanged),
70      reassignLpsolver_(reassignLpsolver),
71      timeBegin_(0),
72      numSols_(0),
73      parameters_(),
74      currentNodeNumber_(-1)
75  {
76    handler_ = new CoinMessageHandler();
77    int logLevel;
78    b.options()->GetIntegerValue("oa_log_level",logLevel,b.prefix());
79    b.options()->GetNumericValue("oa_log_frequency",parameters_.logFrequency_,b.prefix());
80
81    handler_ -> setLogLevel(logLevel);
82    b.options()->GetIntegerValue("solution_limit", parameters_.maxSols_,b.prefix());
83
84    messages_ = OaMessages();
85    timeBegin_ = CoinCpuTime();
86    b.options()->GetIntegerValue("milp_log_level",parameters_.subMilpLogLevel_,b.prefix());
87    b.options()->GetNumericValue("cutoff_decr",parameters_.cbcCutoffIncrement_,b.prefix());
88    b.options()->GetNumericValue("integer_tolerance",parameters_.cbcIntegerTolerance_,b.prefix());
89    int ivalue;
90    b.options()->GetEnumValue("add_only_violated_oa", ivalue,b.prefix());
91    parameters_.addOnlyViolated_ = ivalue;
92    b.options()->GetEnumValue("oa_cuts_scope", ivalue,b.prefix());
93    parameters_.global_ = ivalue;
94  }
95
96  OaDecompositionBase::OaDecompositionBase
97  (const OaDecompositionBase & other)
98      :
99      CglCutGenerator(other),
100      nlp_(other.nlp_),
101      lp_(other.lp_),
102      objects_(other.objects_),
103      nObjects_(other.nObjects_),
104      nLocalSearch_(0),
105      messages_(other.messages_),
106      leaveSiUnchanged_(other.leaveSiUnchanged_),
107      reassignLpsolver_(other.reassignLpsolver_),
108      timeBegin_(0),
109      numSols_(other.numSols_),
110      parameters_(other.parameters_),
111      currentNodeNumber_(other.currentNodeNumber_)
112  {
113    timeBegin_ = CoinCpuTime();
114    handler_ = other.handler_->clone();
115  }
116/// Constructor with default values for parameters
117  OaDecompositionBase::Parameters::Parameters():
118      global_(true),
119      addOnlyViolated_(false),
120      cbcCutoffIncrement_(1e-06),
121      cbcIntegerTolerance_(1e-05),
122      maxLocalSearch_(0),
123      maxLocalSearchTime_(3600),
124      subMilpLogLevel_(0),
125      maxSols_(INT_MAX),
126      logFrequency_(1000.),
127      strategy_(NULL)
128  {}
129
130  /** Destructor.*/
131  OaDecompositionBase::~OaDecompositionBase()
132  {
133    delete handler_;
134  }
135
136
137/// Constructor with default values for parameters
138  OaDecompositionBase::Parameters::Parameters(const Parameters & other):
139      global_(other.global_),
140      addOnlyViolated_(other.addOnlyViolated_),
141      cbcCutoffIncrement_(other.cbcCutoffIncrement_),
142      cbcIntegerTolerance_(other.cbcIntegerTolerance_),
143      maxLocalSearch_(other.maxLocalSearch_),
144      maxLocalSearchTime_(other.maxLocalSearchTime_),
145      subMilpLogLevel_(other.subMilpLogLevel_),
146      maxSols_(other.maxSols_),
147      logFrequency_(other.logFrequency_),
148      strategy_(NULL)
149  {
150    if (other.strategy_)
151      strategy_ = other.strategy_->clone();
152  }
153
154
155
156OaDecompositionBase::solverManip::solverManip
157(OsiSolverInterface * si,
158 bool saveNumRows,
159 bool saveBasis,
160 bool saveBounds,
161 bool saveCutoff,
162 bool resolve):
163    si_(si),
164    initialNumberRows_(-1),
165    colLower_(NULL),
166    colUpper_(NULL),
167    warm_(NULL),
168    cutoff_(COIN_DBL_MAX),
169    deleteSolver_(false),
170    objects_(NULL),
171    nObjects_(0)
172{
173  getCached();
174  if (saveNumRows)
175    initialNumberRows_ = numrows_;
176  if (saveBasis)
177    warm_ = si->getWarmStart();
178  if (saveBounds) {
179    colLower_ = new double[numcols_];
180    colUpper_ = new double[numcols_];
181    CoinCopyN(si->getColLower(), numcols_ , colLower_);
182    CoinCopyN(si->getColUpper(), numcols_ , colUpper_);
183  }
184  if (saveCutoff)
185    si->getDblParam(OsiDualObjectiveLimit, cutoff_);
186  si->messageHandler()->setLogLevel(0);
187  if (resolve) si->resolve();
188}
189
190
191OaDecompositionBase::solverManip::solverManip
192(const OsiSolverInterface & si):
193    si_(NULL),
194    initialNumberRows_(-1),
195    colLower_(NULL),
196    colUpper_(NULL),
197    warm_(NULL),
198    cutoff_(COIN_DBL_MAX),
199    deleteSolver_(true),
200    objects_(NULL),
201    nObjects_(0)
202{
203  si_ = si.clone();
204  getCached();
205}
206
207OaDecompositionBase::solverManip::~solverManip()
208{
209  if (warm_) delete warm_;
210  if (colLower_) delete [] colLower_;
211  if (colUpper_) delete [] colUpper_;
212  if (deleteSolver_) delete si_;
213}
214
215void
216OaDecompositionBase::solverManip::restore()
217{
218  if (initialNumberRows_ >= 0) {
219    int nRowsToDelete = numrows_ - initialNumberRows_;
220    int * rowsToDelete = new int[nRowsToDelete];
221    for (int i = 0 ; i < nRowsToDelete ; i++) {
222      rowsToDelete[i] = i + initialNumberRows_;
223    }
224    si_->deleteRows(nRowsToDelete, rowsToDelete);
225    delete [] rowsToDelete;
226    numrows_ -= nRowsToDelete;
227  }
228
229  if (colLower_) {
230    si_->setColLower(colLower_);
231  }
232
233  if (colUpper_) {
234    si_->setColUpper(colUpper_);
235  }
236
237  if (cutoff_<COIN_DBL_MAX) {
238    si_->setDblParam(OsiDualObjectiveLimit, cutoff_);
239  }
240
241  if (warm_) {
242    if (si_->setWarmStart(warm_)==false) {
243      throw CoinError("Fail restoring the warm start at the end of procedure",
244          "restore","OaDecompositionBase::SaveSolverState") ;
245    }
246  }
247  getCached();
248}
249
250void
251OaDecompositionBase::passInMessageHandler(CoinMessageHandler * handler)
252{
253  int logLevel = handler_->logLevel();
254  delete handler_;
255  handler_=handler->clone();
256  handler_->setLogLevel(logLevel);
257}
258
259
260
261#if 0
262/** Clone the state of another solver (bounds, cutoff, basis).*/
263void
264OaDecompositionBase::solverManip::cloneOther(const OsiSolverInterface &si){
265  //Install current active cuts into local solver
266  int numberCutsToAdd = si.getNumRows();
267  numberCutsToAdd -= numrows_;
268  if (numberCutsToAdd > 0)//Have to install some cuts
269  {
270    CoinPackedVector * * addCuts = new CoinPackedVector *[numberCutsToAdd];
271    for (int i = 0 ; i < numberCutsToAdd ; i++)
272    {
273      addCuts[i] = new CoinPackedVector;
274    }
275    //Get the current matrix and fill the addCuts
276    const CoinPackedMatrix * mat = si.getMatrixByCol();
277    const CoinBigIndex * start = mat->getVectorStarts();
278    const int * length = mat->getVectorLengths();
279    const double * elements = mat->getElements();
280    const int * indices = mat->getIndices();
281    for (int i = 0 ; i < numcols_ ; i++)
282      for (int k = start[i] ; k < start[i] + length[i] ; k++)
283      {
284        if (indices[k] >= numrows_) {
285          addCuts[ indices[k] - numrows_ ]->insert(i, elements[k]);
286        }
287      }
288    si_->addRows(numberCutsToAdd, (const CoinPackedVectorBase * const *) addCuts, si.getRowLower() + numrows_,
289        si.getRowUpper() + numrows_);
290    for (int i = 0 ; i < numberCutsToAdd ; i++){
291      delete addCuts[i];
292    }
293    delete [] addCuts;
294  }
295  else if (numberCutsToAdd < 0) { //Oups some error
296    throw CoinError("Internal error number of cuts wrong",
297        "generateCuts","OACutGenerator2");
298  }
299
300  si_->setColLower(si.getColLower());
301  si_->setColUpper(si.getColUpper());
302  //Install basis in problem
303  CoinWarmStart * warm = si.getWarmStart();
304  if (si_->setWarmStart(warm)==false) {
305    delete warm;
306    throw CoinError("Fail installing the warm start in the subproblem",
307        "generateCuts","OACutGenerator2") ;
308  }
309  delete warm;
310  //put the cutoff
311  double cutoff;
312  si.getDblParam(OsiDualObjectiveLimit, cutoff);
313  si_->setDblParam(OsiDualObjectiveLimit, cutoff);
314  si_->messageHandler()->setLogLevel(0);
315  si_->resolve();
316
317  numrows_ = si.getNumRows();
318#ifdef OA_DEBUG
319
320  std::cout<<"Resolve with hotstart :"<<std::endl
321  <<"number of iterations(should be 0) : "<<si_->getIterationCount()<<std::endl
322  <<"Objective value and diff to original : "<<si_->getObjValue()<<", "
323  <<fabs(si_->getObjValue() - si.getObjValue())<<std::endl;
324  for (int i = 0 ; i <= numcols_ ; i++) {
325    if (fabs(si.getColSolution()[i]-si_->getColSolution()[i])>1e-08) {
326      std::cout<<"Diff between solution at node and solution with local solver : "<<fabs(si.getColSolution()[i]-si_->getColSolution()[i])<<std::endl;
327    }
328  }
329#endif
330
331}
332#endif
333
334
335/** Standard cut generation methods. */
336void
337OaDecompositionBase::generateCuts(const OsiSolverInterface &si,  OsiCuts & cs,
338    const CglTreeInfo info) const{
339  if (nlp_ == NULL) {
340    throw CoinError("Error in cut generator for outer approximation no NLP ipopt assigned", "generateCuts", "OaDecompositionBase");
341  }
342
343  // babInfo is used to communicate with the b-and-b solver (Cbc or Bcp).
344  BabInfo * babInfo = dynamic_cast<BabInfo *> (si.getAuxiliaryInfo());
345  assert(babInfo);
346  assert(babInfo->babPtr());
347  numSols_ = babInfo->babPtr()->model().getSolutionCount ();
348  CglTreeInfo info_copy = info;
349  const CbcNode * node = babInfo->babPtr()->model().currentNode();
350  info_copy.level = (node == NULL) ? 0 : babInfo->babPtr()->model().currentNode()->depth();
351  if(babInfo->hasSolution()) numSols_ ++;
352  if (babInfo)
353    if (!babInfo->mipFeasible())
354      return;
355
356  //Get the continuous solution
357  const double *colsol = si.getColSolution();
358
359
360  vector<double> savedColLower(nlp_->getNumCols());
361  CoinCopyN(nlp_->getColLower(), nlp_->getNumCols(), savedColLower());
362  vector<double> savedColUpper(nlp_->getNumCols());
363  CoinCopyN(nlp_->getColUpper(), nlp_->getNumCols(), savedColUpper());
364
365
366  OsiBranchingInformation brInfo(nlp_, false);
367  brInfo.solution_ = colsol;
368  //Check integer infeasibility
369  bool isInteger = integerFeasible(*nlp_, brInfo, parameters_.cbcIntegerTolerance_,
370                              objects_, nObjects_);
371
372  SubMipSolver * subMip = NULL;
373
374  //Check nodeNumber if it did not change scan savedCuts_ if one is violated force it and exit
375  int nodeNumber = babInfo->babPtr()->model().getNodeCount();
376  if(nodeNumber == currentNodeNumber_){
377#ifdef OA_DEBUG
378    printf("OA decomposition recalled from the same node!\n");
379#endif
380    int numCuts = savedCuts_.sizeRowCuts();
381    for(int i = 0 ; i < numCuts ; i++){
382       //Check if cuts off solution
383       if(savedCuts_.rowCut(i).violated(colsol) > 0.){
384#ifdef OA_DEBUG
385         printf("A violated cut has been found\n");
386#endif
387         savedCuts_.rowCut(i).setEffectiveness(9.99e99);
388         cs.insert(savedCuts_.rowCut(i));
389         savedCuts_.eraseRowCut(i);
390         return;
391         i--; numCuts--;
392       }
393    }
394  }
395  else {
396    currentNodeNumber_ = nodeNumber;
397    savedCuts_.dumpCuts();
398  } 
399         
400  if (!isInteger) {
401    if (!doLocalSearch(babInfo))//create sub mip solver.
402      return;
403  }
404
405  //get the current cutoff
406  double cutoff;
407  si.getDblParam(OsiDualObjectiveLimit, cutoff);
408
409  // Save solvers state if needed
410
411  solverManip * lpManip = NULL;
412  if (lp_ != NULL) {
413    if (lp_!=&si) {
414      throw;
415#if 0
416      lpManip = new solverManip(lp_, true, false, false, true, true);
417      lpManip->cloneOther(si);
418#endif
419    }
420    else {
421#if 0
422      throw CoinError("Not allowed to modify si in a cutGenerator",
423          "OACutGenerator2","generateCuts");
424#else
425      lpManip = new solverManip(lp_, true, leaveSiUnchanged_, true, true);
426#endif
427    }
428  }
429  else {
430    lpManip = new solverManip(si);
431  }
432  lpManip->setObjects(objects_, nObjects_);
433  if (!isInteger) {
434      subMip = new SubMipSolver(lpManip->si(), parameters_.strategy());
435  }
436
437  double milpBound = performOa(cs, *lpManip, subMip, babInfo, cutoff, info_copy);
438
439  if(babInfo->hasSolution()){
440     babInfo->babPtr()->model().setSolutionCount (numSols_ - 1);
441  }
442
443  //Transmit the bound found by the milp
444  {
445    if (milpBound>-1e100)
446    {
447      // Also store into solver
448      if (babInfo)
449        babInfo->setMipBound(milpBound);
450    }
451  }  //Clean everything :
452
453  //free subMip
454  if (subMip!= NULL) {
455    delete subMip;
456    subMip = NULL;
457  }
458
459  //  Reset the two solvers
460  if (leaveSiUnchanged_)
461    lpManip->restore();
462  delete lpManip;
463
464  nlp_->setColLower(savedColLower());
465  nlp_->setColUpper(savedColUpper());
466
467  return;
468}
469
470void
471OaDecompositionBase::solverManip::getCached(){
472  numrows_ = si_->getNumRows();
473  numcols_ = si_->getNumCols();
474  siColLower_ = si_->getColLower();
475  siColUpper_ = si_->getColUpper();
476}
477
478
479/** Do update after an nlp has been solved*/
480bool
481OaDecompositionBase::post_nlp_solve(BabInfo * babInfo, double cutoff) const{
482  nSolve_++;
483  bool return_value = false;
484  if (nlp_->isProvenOptimal()) {
485    handler_->message(FEASIBLE_NLP, messages_)
486    <<nlp_->getIterationCount()
487    <<nlp_->getObjValue()<<CoinMessageEol;
488
489#ifdef OA_DEBUG
490    const double * colsol2 = nlp_->getColSolution();
491    debug_.checkInteger(*nlp_,std::cerr);
492#endif
493
494    if ((nlp_->getObjValue() < cutoff) ) {
495      handler_->message(UPDATE_UB, messages_)
496      <<nlp_->getObjValue()
497      <<CoinCpuTime()-timeBegin_
498      <<CoinMessageEol;
499
500      return_value = true;
501      // Also pass it to solver
502      assert(babInfo);
503      if (babInfo) {
504        int numcols = nlp_->getNumCols();
505        double * lpSolution = new double[numcols + 1];
506        CoinCopyN(nlp_->getColSolution(), numcols, lpSolution);
507        lpSolution[numcols] = nlp_->getObjValue();
508        babInfo->setSolution(lpSolution,
509            numcols + 1, lpSolution[numcols]);
510        delete [] lpSolution;
511      }
512    }
513  }
514  else if (nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
515    (*handler_)<<"Unsolved NLP... exit"<<CoinMessageEol;
516  }
517  else {
518    handler_->message(INFEASIBLE_NLP, messages_)
519    <<nlp_->getIterationCount()
520    <<CoinMessageEol;
521  }
522  return return_value;
523}
524
525
526
527#ifdef OA_DEBUG
528bool
529OaDecompositionBase::OaDebug::checkInteger(const OsiSolverInterface &nlp, 
530                                           std::ostream & os) const {
531   const double * colsol = nlp.getColSolution();
532   int numcols = nlp.getNumCols();
533  for (int i = 0 ; i < numcols ; i++) {
534    if (nlp.isInteger(i)) {
535      if (fabs(colsol[i]) - floor(colsol[i] + 0.5) >
536          1e-07) {
537        std::cerr<<"Integer infeasible point (should not be), integer infeasibility for variable "<<i
538        <<" is, "<<fabs(colsol[i] - floor(colsol[i] + 0.5))<<std::endl;
539      }
540    }
541    return true;
542  }
543
544}
545
546void
547OaDecompositionBase::OaDebug::printEndOfProcedureDebugMessage(const OsiCuts &cs,
548    bool foundSolution,
549    double solValue,
550    double milpBound,
551    bool isInteger,
552    bool feasible,
553    std::ostream & os) const{
554  std::cout<<"------------------------------------------------------------------"
555  <<std::endl;
556  std::cout<<"OA procedure finished"<<std::endl;
557  std::cout<<"Generated "<<cs.sizeRowCuts()<<std::endl;
558  if (foundSolution)
559    std::cout <<"Found NLP-integer feasible solution of  value : "<<solValue<<std::endl;
560  std::cout<<"Current MILP lower bound is : "<<milpBound<<std::endl;
561  std::cout<<"-------------------------------------------------------------------"<<std::endl;
562  std::cout<<"Stopped because : isInteger "<<isInteger<<", feasible "<<feasible<<std::endl<<std::endl;
563
564}
565
566
567
568#endif
569}/* End namespace Bonmin. */
570
Note: See TracBrowser for help on using the repository browser.