source: stable/1.1/Bonmin/src/Algorithms/OaGenerators/BonOaDecBase.cpp @ 1524

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

Remove some unused option in OA add Claudi d'Ambrosio code

File size: 24.0 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#include <algorithm>
13
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 "OsiAuxInfo.hpp"
30
31//The following two are to interupt the solution of sub-mip through CTRL-C
32extern CbcModel * OAModel;
33
34namespace Bonmin
35{
36
37  OaDecompositionBase::OaDecompositionBase
38  (OsiTMINLPInterface * nlp,
39   OsiSolverInterface * si,
40   CbcStrategy * strategy,
41   double cbcCutoffIncrement,
42   double cbcIntegerTolerance,
43   bool leaveSiUnchanged
44  )
45      :
46      CglCutGenerator(),
47      nlp_(nlp),
48      nSolve_(0),
49      lp_(si),
50      objects_(NULL),
51      nObjects_(0),
52      nLocalSearch_(0),
53      handler_(NULL),
54      leaveSiUnchanged_(leaveSiUnchanged),
55      reassignLpsolver_(false),
56      timeBegin_(0),
57      parameters_()
58  {
59    handler_ = new CoinMessageHandler();
60    handler_ -> setLogLevel(2);
61    messages_ = OaMessages();
62    if (strategy)
63      parameters_.setStrategy(*strategy);
64    timeBegin_ = CoinCpuTime();
65    parameters_.cbcCutoffIncrement_  = cbcCutoffIncrement;
66    parameters_.cbcIntegerTolerance_ = cbcIntegerTolerance;
67  }
68  OaDecompositionBase::OaDecompositionBase(BabSetupBase &b, bool leaveSiUnchanged,
69      bool reassignLpsolver):
70      CglCutGenerator(),
71      nlp_(b.nonlinearSolver()),
72      lp_(b.continuousSolver()),
73      objects_(NULL),
74      nObjects_(0),
75      nLocalSearch_(0),
76      handler_(NULL),
77      leaveSiUnchanged_(leaveSiUnchanged),
78      reassignLpsolver_(reassignLpsolver),
79      timeBegin_(0),
80      parameters_()
81  {
82    handler_ = new CoinMessageHandler();
83    int logLevel;
84    b.options()->GetIntegerValue("oa_log_level",logLevel,"bonmin.");
85    b.options()->GetNumericValue("oa_log_frequency",parameters_.logFrequency_,"bonmin.");
86
87    handler_ -> setLogLevel(logLevel);
88
89    messages_ = OaMessages();
90    timeBegin_ = CoinCpuTime();
91    b.options()->GetIntegerValue("milp_log_level",parameters_.subMilpLogLevel_,"bonmin.");
92    b.options()->GetNumericValue("cutoff_decr",parameters_.cbcCutoffIncrement_,"bonmin.");
93    b.options()->GetNumericValue("integer_tolerance",parameters_.cbcIntegerTolerance_,"bonmin.");
94    int ivalue;
95    b.options()->GetEnumValue("add_only_violated_oa", ivalue,"bonmin.");
96    parameters_.addOnlyViolated_ = ivalue;
97    b.options()->GetEnumValue("oa_cuts_scope", ivalue,"bonmin.");
98    parameters_.global_ = ivalue;
99  }
100
101  OaDecompositionBase::OaDecompositionBase
102  (const OaDecompositionBase & other)
103      :
104      CglCutGenerator(other),
105      nlp_(other.nlp_),
106      lp_(other.lp_),
107      objects_(other.objects_),
108      nObjects_(other.nObjects_),
109      nLocalSearch_(0),
110      messages_(other.messages_),
111      leaveSiUnchanged_(other.leaveSiUnchanged_),
112      reassignLpsolver_(other.reassignLpsolver_),
113      timeBegin_(other.timeBegin_),
114      parameters_(other.parameters_)
115  {
116    //timeBegin_ = CoinCpuTime();
117    handler_ = other.handler_->clone();
118  }
119/// Constructor with default values for parameters
120  OaDecompositionBase::Parameters::Parameters():
121      global_(true),
122      addOnlyViolated_(false),
123      cbcCutoffIncrement_(1e-06),
124      cbcIntegerTolerance_(1e-05),
125      maxLocalSearch_(0),
126      maxLocalSearchTime_(3600),
127      subMilpLogLevel_(0),
128      logFrequency_(1000.),
129      strategy_(NULL)
130  {}
131
132  /** Destructor.*/
133  OaDecompositionBase::~OaDecompositionBase()
134  {
135    delete handler_;
136  }
137
138
139/// Constructor with default values for parameters
140  OaDecompositionBase::Parameters::Parameters(const Parameters & other):
141      global_(other.global_),
142      addOnlyViolated_(other.addOnlyViolated_),
143      cbcCutoffIncrement_(other.cbcCutoffIncrement_),
144      cbcIntegerTolerance_(other.cbcIntegerTolerance_),
145      maxLocalSearch_(other.maxLocalSearch_),
146      maxLocalSearchTime_(other.maxLocalSearchTime_),
147      subMilpLogLevel_(other.subMilpLogLevel_),
148      logFrequency_(other.logFrequency_),
149      strategy_(NULL)
150  {
151    if (other.strategy_)
152      strategy_ = other.strategy_->clone();
153  }
154
155
156  /** Constructor */
157  OaDecompositionBase::SubMipSolver::SubMipSolver(OsiSolverInterface * lp,
158      const CbcStrategy * strategy):
159      lp_(lp),
160      clp_(NULL),
161      cpx_(NULL),
162      cbc_(NULL),
163      lowBound_(-COIN_DBL_MAX),
164      optimal_(false),
165      integerSolution_(NULL),
166      strategy_(NULL)
167  {
168    clp_ = (lp_ == NULL)? NULL :
169        dynamic_cast<OsiClpSolverInterface *>(lp_);
170#ifdef COIN_HAS_CPX
171    cpx_ = (lp_ == NULL)? NULL :
172        dynamic_cast<OsiCpxSolverInterface *>(lp_);
173#endif
174    if (strategy) strategy_ = strategy->clone();
175  }
176  OaDecompositionBase::SubMipSolver::~SubMipSolver()
177  {
178    if (strategy_) delete strategy_;
179    if (integerSolution_) delete [] integerSolution_;
180    if (cbc_) delete cbc_;
181  }
182
183  /** Assign lp solver. */
184  void
185  OaDecompositionBase::SubMipSolver::setLpSolver(OsiSolverInterface * lp)
186  {
187    lp_ = lp;
188    clp_ = (lp_ == NULL) ? NULL :
189        dynamic_cast<OsiClpSolverInterface *>(lp_);
190#ifdef COIN_HAS_CPX
191    cpx_ = (lp_ == NULL) ? NULL :
192        dynamic_cast<OsiCpxSolverInterface *>(lp_);
193#endif
194    lowBound_ = -COIN_DBL_MAX;
195    optimal_ = false;
196    if (integerSolution_) {
197      delete [] integerSolution_;
198      integerSolution_ = NULL;
199    }
200  }
201
202
203
204  void
205  OaDecompositionBase::SubMipSolver::performLocalSearch(double cutoff, int loglevel, double maxTime)
206  {
207    if (clp_) {
208      if (!strategy_)
209        strategy_ = new CbcStrategyDefault(1,0,0, loglevel);
210
211      OsiBabSolver empty;
212      if (cbc_) delete cbc_;
213      OAModel = cbc_ = new CbcModel(*clp_);
214      cbc_->solver()->setAuxiliaryInfo(&empty);
215
216      //Change Cbc messages prefixes
217      strcpy(cbc_->messagesPointer()->source_,"OaCbc");
218
219      cbc_->setLogLevel(loglevel);
220      cbc_->solver()->messageHandler()->setLogLevel(0);
221      clp_->resolve();
222      cbc_->setStrategy(*strategy_);
223      cbc_->setLogLevel(loglevel);
224      cbc_->solver()->messageHandler()->setLogLevel(0);
225      cbc_->setMaximumSeconds(maxTime);
226      cbc_->setCutoff(cutoff);
227      cbc_->branchAndBound();
228      OAModel = NULL;
229      lowBound_ = cbc_->getBestPossibleObjValue();
230
231      if (cbc_->isProvenOptimal() || cbc_->isProvenInfeasible())
232        optimal_ = true;
233      else optimal_ = false;
234
235      if (cbc_->getSolutionCount()) {
236        if (!integerSolution_)
237          integerSolution_ = new double[lp_->getNumCols()];
238        CoinCopyN(cbc_->bestSolution(), lp_->getNumCols(), integerSolution_);
239      }
240      else if (integerSolution_) {
241        delete [] integerSolution_;
242        integerSolution_ = NULL;
243      }
244      nodeCount_ = cbc_->getNodeCount();
245      iterationCount_ = cbc_->getIterationCount();
246    }
247    else {
248      lp_->messageHandler()->setLogLevel(loglevel);
249#ifdef COIN_HAS_CPX
250      if (cpx_) {
251        CPXENVptr env = cpx_->getEnvironmentPtr();
252        CPXsetdblparam(env, CPX_PARAM_TILIM, maxTime);
253        CPXsetdblparam(env, CPX_PARAM_CUTUP, cutoff);
254        //CpxModel = cpx_;
255      }
256      else
257#endif
258     {
259        throw CoinError("Unsuported solver, for local searches you should use clp or cplex",
260            "performLocalSearch",
261            "OaDecompositionBase::SubMipSolver");
262    }
263
264#ifdef COIN_HAS_CPX
265    if (cpx_) {
266      cpx_->switchToMIP();
267      //CpxModel = NULL;
268      CPXENVptr env = cpx_->getEnvironmentPtr();
269      CPXLPptr cpxlp = cpx_->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
270
271      int status = CPXmipopt(env,cpxlp);
272      CHECK_CPX_STAT("mipopt",status)
273
274      int stat = CPXgetstat( env, cpxlp);
275//#define FOR_HASSAN
276#ifdef FOR_HASSAN
277      std::cout << "STATUS = "<<stat<<"\n";
278      if(stat==119){
279          cpx_->writeMpsNative("OA.mps", NULL, NULL, 1);
280          throw CoinError("Ok program stoped by me !","OaDecompositionBase::SubMipSolver","performLocalSearch");
281      }
282#endif
283      bool infeasible = (stat == CPXMIP_INFEASIBLE) || (stat == CPXMIP_ABORT_INFEAS) || (stat == CPXMIP_TIME_LIM_INFEAS) || (stat == CPXMIP_NODE_LIM_INFEAS) || (stat == CPXMIP_FAIL_INFEAS)
284                        || (stat == CPXMIP_MEM_LIM_INFEAS) || (stat == CPXMIP_INForUNBD);
285      optimal_ = (stat == CPXMIP_INFEASIBLE) || (stat == CPXMIP_OPTIMAL) || (stat == CPXMIP_OPTIMAL_TOL); 
286      nodeCount_ = CPXgetnodecnt(env , cpxlp);
287      iterationCount_ = CPXgetmipitcnt(env , cpxlp);
288      status = CPXgetbestobjval(env, cpxlp, &lowBound_);
289      CHECK_CPX_STAT("getbestobjval",status)
290       
291      if(!infeasible){
292         if(!integerSolution_){
293           integerSolution_ = new double[lp_->getNumCols()];
294         }
295         CPXgetmipx(env, cpxlp, integerSolution_, 0, lp_->getNumCols() -1);
296         CHECK_CPX_STAT("getmipx",status)
297      }
298      else {
299        if (integerSolution_) {
300          delete [] integerSolution_;
301          integerSolution_ = NULL;
302        }
303      }
304    }
305    else {
306#endif
307    lp_->branchAndBound();
308
309   optimal_ = lp_->isProvenOptimal();
310
311    if (lp_->getFractionalIndices().size() == 0) {
312      if (!integerSolution_)
313        integerSolution_ = new double[lp_->getNumCols()];
314      CoinCopyN(lp_->getColSolution(), lp_->getNumCols() , integerSolution_);
315    }
316    else if (integerSolution_) {
317      delete [] integerSolution_;
318      integerSolution_ = NULL;
319    }
320  }
321#ifdef COIN_HAS_CPX
322  }
323#endif
324}
325
326OaDecompositionBase::solverManip::solverManip
327(OsiSolverInterface * si,
328 bool saveNumRows,
329 bool saveBasis,
330 bool saveBounds,
331 bool saveCutoff,
332 bool resolve):
333    si_(si),
334    initialNumberRows_(-1),
335    colLower_(NULL),
336    colUpper_(NULL),
337    warm_(NULL),
338    cutoff_(COIN_DBL_MAX),
339    deleteSolver_(false),
340    objects_(NULL),
341    nObjects_(0),
342    integerTolerance_(1e-08)
343{
344  getCached();
345  if (saveNumRows)
346    initialNumberRows_ = numrows_;
347  if (saveBasis)
348    warm_ = si->getWarmStart();
349  if (saveBounds) {
350    colLower_ = new double[numcols_];
351    colUpper_ = new double[numcols_];
352    CoinCopyN(si->getColLower(), numcols_ , colLower_);
353    CoinCopyN(si->getColUpper(), numcols_ , colUpper_);
354  }
355  if (saveCutoff)
356    si->getDblParam(OsiDualObjectiveLimit, cutoff_);
357  si->messageHandler()->setLogLevel(0);
358  if (resolve) si->resolve();
359}
360
361
362OaDecompositionBase::solverManip::solverManip
363(const OsiSolverInterface & si):
364    si_(NULL),
365    initialNumberRows_(-1),
366    colLower_(NULL),
367    colUpper_(NULL),
368    warm_(NULL),
369    cutoff_(COIN_DBL_MAX),
370    deleteSolver_(true),
371    objects_(NULL),
372    nObjects_(0),
373    integerTolerance_(1e-08)
374{
375  si_ = si.clone();
376  getCached();
377}
378
379OaDecompositionBase::solverManip::~solverManip()
380{
381  if (warm_) delete warm_;
382  if (colLower_) delete [] colLower_;
383  if (colUpper_) delete [] colUpper_;
384  if (deleteSolver_) delete si_;
385}
386
387void
388OaDecompositionBase::solverManip::restore()
389{
390  if (initialNumberRows_ >= 0) {
391    int nRowsToDelete = numrows_ - initialNumberRows_;
392    int * rowsToDelete = new int[nRowsToDelete];
393    for (int i = 0 ; i < nRowsToDelete ; i++) {
394      rowsToDelete[i] = i + initialNumberRows_;
395    }
396    si_->deleteRows(nRowsToDelete, rowsToDelete);
397    delete [] rowsToDelete;
398    numrows_ -= nRowsToDelete;
399  }
400
401  if (colLower_) {
402    si_->setColLower(colLower_);
403  }
404
405  if (colUpper_) {
406    si_->setColUpper(colUpper_);
407  }
408
409  if (cutoff_<COIN_DBL_MAX) {
410    si_->setDblParam(OsiDualObjectiveLimit, cutoff_);
411  }
412
413  if (warm_) {
414    if (si_->setWarmStart(warm_)==false) {
415      throw CoinError("Fail restoring the warm start at the end of procedure",
416          "restore","OaDecompositionBase::SaveSolverState") ;
417    }
418  }
419  getCached();
420}
421
422void
423OaDecompositionBase::passInMessageHandler(CoinMessageHandler * handler)
424{
425  int logLevel = handler_->logLevel();
426  delete handler_;
427  handler_=handler->clone();
428  handler_->setLogLevel(logLevel);
429}
430
431/// Check for integer feasibility of a solution return 1 if it is
432bool
433OaDecompositionBase::solverManip::integerFeasible(const OsiBranchingInformation& info) const
434{
435  if (objects_) {
436    int dummy;
437    for (int i = 0 ; i < nObjects_ ; i++) {
438      if (objects_[i]->infeasibility(&info, dummy) > 0.0) return false;
439    }
440  }
441  else {
442    const double * sol = info.solution_;
443    int numcols = si_->getNumCols();
444    for (int i = 0 ; i < numcols ; i++) {
445      if (si_->isInteger(i)) {
446        if (fabs(sol[i] - floor(sol[i] + 0.5)) >
447            integerTolerance_) {
448          return false;
449        }
450      }
451    }
452  }
453  return true;
454}
455
456/** Fix integer variables in si to their values in colsol.
457*/
458void
459OaDecompositionBase::solverManip::fixIntegers(const OsiBranchingInformation& info)
460{
461  if (objects_) {
462    for (int i = 0 ; i < nObjects_ ; i++) {
463      if (objects_[i]->feasibleRegion(si_, &info));
464    }
465  }
466  else {
467    const double * colsol = info.solution_;
468    for (int i = 0; i < numcols_; i++) {
469      if (si_->isInteger(i)) {
470        double  value =  colsol[i];
471        if (fabs(value - floor(value+0.5)) > integerTolerance_) {
472          std::stringstream stream;
473          stream<<"Error not integer valued solution"<<std::endl;
474          stream<<"---------------- x["<<i<<"] = "<<value<<std::endl;
475          throw CoinError(stream.str(),"fixIntegers","OaDecompositionBase::solverManip");
476        }
477        value = floor(value+0.5);
478        value = std::max(colLower_[i],value);
479        value = std::min(value, colUpper_[i]);
480
481        if (fabs(value) > 1e10) {
482          std::stringstream stream;
483          stream<<"Can not fix variable in nlp because it has too big a value ("<<value
484          <<") at optimium of LP relaxation. You should try running the problem with B-BB"<<std::endl;
485          throw CoinError(stream.str(),
486              "fixIntegers","OaDecompositionBase::solverManip") ;
487        }
488#ifdef OA_DEBUG
489        //         printf("xx %d at %g (bounds %g, %g)",i,value,nlp_->getColLower()[i],
490        //                nlp_->getColUpper()[i]);
491        std::cout<<(int)value;
492#endif
493        si_->setColLower(i,value);
494        si_->setColUpper(i,value);
495      }
496    }
497#ifdef OA_DEBUG
498    std::cout<<std::endl;
499#endif
500  }
501}
502/** Check if solution in solver is the same as colsol on integer variables. */
503bool
504OaDecompositionBase::solverManip::isDifferentOnIntegers(const double * colsol)
505{
506  return isDifferentOnIntegers(colsol, si_->getColSolution());
507}
508
509/** Check if two solutions are the same on integer variables. */
510bool
511OaDecompositionBase::solverManip::isDifferentOnIntegers(const double * colsol, const double *otherSol)
512{
513  if (objects_) {
514    for (int i = 0 ; i < nObjects_ ; i++) {
515      OsiObject * obj = objects_[i];
516      int colnum = obj->columnNumber();
517      if (colnum >= 0) {//Variable branching object
518        if (fabs(otherSol[colnum] - colsol[colnum]) > 100*integerTolerance_) {
519          return true;
520        }
521      }
522      else {//It is a sos branching object
523        OsiSOS * sos = dynamic_cast<OsiSOS *>(obj);
524        assert(sos);
525        const int * members = sos->members();
526        int end = sos->numberMembers();
527        for (int k = 0 ; k < end ; k++) {
528          if (fabs(otherSol[members[k]] - colsol[members[k]]) > 0.001) {
529            return true;
530          }
531        }
532      }
533    }
534  }
535  else {
536    for (int i = 0; i < numcols_ ; i++) {
537      if (si_->isInteger(i) && fabs(otherSol[i] - colsol[i])>0.001)
538        return true;
539    }
540  }
541  return false;
542}
543
544/** Clone the state of another solver (bounds, cutoff, basis).*/
545void
546OaDecompositionBase::solverManip::cloneOther(const OsiSolverInterface &si)
547{
548  //Install current active cuts into local solver
549  int numberCutsToAdd = si.getNumRows();
550  numberCutsToAdd -= numrows_;
551  if (numberCutsToAdd > 0)//Have to install some cuts
552  {
553    CoinPackedVector * * addCuts = new CoinPackedVector *[numberCutsToAdd];
554    for (int i = 0 ; i < numberCutsToAdd ; i++)
555    {
556      addCuts[i] = new CoinPackedVector;
557    }
558    //Get the current matrix and fill the addCuts
559    const CoinPackedMatrix * mat = si.getMatrixByCol();
560    const CoinBigIndex * start = mat->getVectorStarts();
561    const int * length = mat->getVectorLengths();
562    const double * elements = mat->getElements();
563    const int * indices = mat->getIndices();
564    for (int i = 0 ; i < numcols_ ; i++)
565      for (int k = start[i] ; k < start[i] + length[i] ; k++)
566      {
567        if (indices[k] >= numrows_) {
568          addCuts[ indices[k] - numrows_ ]->insert(i, elements[k]);
569        }
570      }
571    si_->addRows(numberCutsToAdd, (const CoinPackedVectorBase * const *) addCuts, si.getRowLower() + numrows_,
572        si.getRowUpper() + numrows_);
573    for (int i = 0 ; i < numberCutsToAdd ; i++){
574      delete addCuts[i];
575    }
576    delete [] addCuts;
577  }
578  else if (numberCutsToAdd < 0) { //Oups some error
579    throw CoinError("Internal error number of cuts wrong",
580        "generateCuts","OACutGenerator2");
581  }
582
583  si_->setColLower(si.getColLower());
584  si_->setColUpper(si.getColUpper());
585  //Install basis in problem
586  CoinWarmStart * warm = si.getWarmStart();
587  if (si_->setWarmStart(warm)==false) {
588    delete warm;
589    throw CoinError("Fail installing the warm start in the subproblem",
590        "generateCuts","OACutGenerator2") ;
591  }
592  delete warm;
593  //put the cutoff
594  double cutoff;
595  si.getDblParam(OsiDualObjectiveLimit, cutoff);
596  si_->setDblParam(OsiDualObjectiveLimit, cutoff);
597  si_->messageHandler()->setLogLevel(0);
598  si_->resolve();
599
600  numrows_ = si.getNumRows();
601#ifdef OA_DEBUG
602
603  std::cout<<"Resolve with hotstart :"<<std::endl
604  <<"number of iterations(should be 0) : "<<lp_->getIterationCount()<<std::endl
605  <<"Objective value and diff to original : "<<lp_->getObjValue()<<", "
606  <<fabs(si_->getObjValue() - si.getObjValue())<<std::endl;
607  for (int i = 0 ; i <= numcols ; i++) {
608    if (fabs(si.getColSolution()[i]-si_->getColSolution()[i])>1e-08) {
609      std::cout<<"Diff between solution at node and solution with local solver : "<<fabs(si.getColSolution()[i]-si_->getColSolution()[i])<<std::endl;
610    }
611  }
612#endif
613
614}
615
616
617/** Install cuts in solver. */
618void
619OaDecompositionBase::solverManip::installCuts(const OsiCuts& cs, int numberCuts)
620{
621  int numberCutsBefore = cs.sizeRowCuts() - numberCuts;
622
623  CoinWarmStartBasis * basis
624  = dynamic_cast<CoinWarmStartBasis*>(si_->getWarmStart()) ;
625  assert(basis != NULL); // make sure not volume
626  basis->resize(numrows_ + numberCuts,numcols_) ;
627  for (int i = 0 ; i < numberCuts ; i++) {
628    basis->setArtifStatus(numrows_ + i,
629        CoinWarmStartBasis::basic) ;
630  }
631
632  const OsiRowCut ** addCuts = new const OsiRowCut * [numberCuts] ;
633  for (int i = 0 ; i < numberCuts ; i++) {
634    addCuts[i] = &cs.rowCut(i + numberCutsBefore) ;
635  }
636  si_->applyRowCuts(numberCuts,addCuts) ;
637  numrows_ += numberCuts;
638  delete [] addCuts ;
639  if (si_->setWarmStart(basis) == false) {
640    delete basis;
641    throw CoinError("Fail setWarmStart() after cut installation.",
642        "generateCuts","OACutGenerator2") ;
643  }
644  delete basis;
645}
646
647/** Standard cut generation methods. */
648void
649OaDecompositionBase::generateCuts(const OsiSolverInterface &si,  OsiCuts & cs,
650    const CglTreeInfo info) const
651{
652  if (nlp_ == NULL) {
653    throw CoinError("Error in cut generator for outer approximation no NLP ipopt assigned", "generateCuts", "OaDecompositionBase");
654  }
655
656  // babInfo is used to communicate with the b-and-b solver (Cbc or Bcp).
657  OsiBabSolver * babInfo = dynamic_cast<OsiBabSolver *> (si.getAuxiliaryInfo());
658  if (babInfo)
659    if (!babInfo->mipFeasible())
660      return;
661
662  //Get the continuous solution
663  const double *colsol = si.getColSolution();
664
665
666  solverManip nlpManip(nlp_, false, false, true, false, false);
667  nlpManip.setIntegerTolerance(parameters_.cbcIntegerTolerance_);
668  nlpManip.setObjects(objects_, nObjects_);
669  OsiBranchingInformation brInfo(nlp_, false);
670  brInfo.solution_ = colsol;
671  //Check integer infeasibility
672  bool isInteger = nlpManip.integerFeasible(brInfo);
673
674  SubMipSolver * subMip = NULL;
675
676  if (!isInteger) {
677    if (doLocalSearch())//create sub mip solver.
678    {
679      subMip = new SubMipSolver(lp_, parameters_.strategy());
680    }
681    else {
682      return;
683    }
684  }
685
686
687  //If we are going to modify things copy current information to restore it in the end
688
689
690  //get the current cutoff
691  double cutoff;
692  si.getDblParam(OsiDualObjectiveLimit, cutoff);
693
694  // Save solvers state if needed
695
696  solverManip * lpManip = NULL;
697  if (lp_ != NULL) {
698    if (lp_!=&si) {
699      lpManip = new solverManip(lp_, true, false, false, true, true);
700      lpManip->cloneOther(si);
701    }
702    else {
703#if 0
704      throw CoinError("Not allowed to modify si in a cutGenerator",
705          "OACutGenerator2","generateCuts");
706#else
707      lpManip = new solverManip(lp_, true, leaveSiUnchanged_, true, true);
708#endif
709    }
710  }
711  else {
712    lpManip = new solverManip(si);
713  }
714  lpManip->setObjects(objects_, nObjects_);
715  lpManip->setIntegerTolerance(parameters_.cbcIntegerTolerance_);
716  double milpBound = performOa(cs, nlpManip, *lpManip, subMip, babInfo, cutoff);
717
718  //Transmit the bound found by the milp
719  {
720    if (milpBound>-1e100)
721    {
722      // Also store into solver
723      if (babInfo)
724        babInfo->setMipBound(milpBound);
725    }
726  }  //Clean everything :
727
728  //free subMip
729  if (subMip!= NULL) {
730    delete subMip;
731    subMip = NULL;
732  }
733
734  //  Reset the two solvers
735  if (leaveSiUnchanged_)
736    lpManip->restore();
737  delete lpManip;
738  nlpManip.restore();
739  return;
740}
741
742void
743OaDecompositionBase::solverManip::getCached()
744{
745  numrows_ = si_->getNumRows();
746  numcols_ = si_->getNumCols();
747  siColLower_ = si_->getColLower();
748  siColUpper_ = si_->getColUpper();
749}
750
751
752/** Solve the nlp and do output return true if feasible*/
753bool
754OaDecompositionBase::solveNlp(OsiBabSolver * babInfo, double cutoff) const
755{
756  nSolve_++;
757  nlp_->resolve();
758  bool return_value = false;
759  if (nlp_->isProvenOptimal()) {
760    handler_->message(FEASIBLE_NLP, messages_)
761    <<nlp_->getIterationCount()
762    <<nlp_->getObjValue()<<CoinMessageEol;
763
764#ifdef OA_DEBUG
765    const double * colsol2 = nlp_->getColSolution();
766    debug_.checkInteger(*nlp_,std::cerr);
767#endif
768
769    if ((nlp_->getObjValue() < cutoff) ) {
770      handler_->message(UPDATE_UB, messages_)
771      <<nlp_->getObjValue()
772      <<CoinCpuTime()-timeBegin_
773      <<CoinMessageEol;
774
775      return_value = true;
776      // Also pass it to solver
777      if (babInfo) {
778        int numcols = nlp_->getNumCols();
779        double * lpSolution = new double[numcols + 1];
780        CoinCopyN(nlp_->getColSolution(), numcols, lpSolution);
781        lpSolution[numcols] = nlp_->getObjValue();
782        babInfo->setSolution(lpSolution,
783            numcols + 1, lpSolution[numcols]);
784        delete [] lpSolution;
785      }
786      else {
787        printf("No auxiliary info in nlp solve!\n");
788        throw -1;
789      }
790    }
791  }
792  else if (nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
793    (*handler_)<<"Unsolved NLP... exit"<<CoinMessageEol;
794  }
795  else {
796    handler_->message(INFEASIBLE_NLP, messages_)
797    <<nlp_->getIterationCount()
798    <<CoinMessageEol;
799  }
800  return return_value;
801}
802
803
804
805#ifdef OA_DEBUG
806bool
807OaDecompositionBase::OaDebug::checkInteger(const OsiSolverInterface &nlp, std::ostream & os) const
808{
809   const double * colsol = nlp.getColSolution();
810   int numcols = nlp.getNumCols();
811  for (int i = 0 ; i < numcols ; i++) {
812    if (nlp.isInteger(i)) {
813      if (fabs(colsol[i]) - floor(colsol[i] + 0.5) >
814          1e-07) {
815        std::cerr<<"Integer infeasible point (should not be), integer infeasibility for variable "<<i
816        <<" is, "<<fabs(colsol[i] - floor(colsol[i] + 0.5))<<std::endl;
817      }
818    }
819    return true;
820  }
821
822}
823
824void
825OaDecompositionBase::OaDebug::printEndOfProcedureDebugMessage(const OsiCuts &cs,
826    bool foundSolution,
827    double solValue,
828    double milpBound,
829    bool isInteger,
830    bool feasible,
831    std::ostream & os) const
832{
833  std::cout<<"------------------------------------------------------------------"
834  <<std::endl;
835  std::cout<<"OA procedure finished"<<std::endl;
836  std::cout<<"Generated "<<cs.sizeRowCuts()<<std::endl;
837  if (foundSolution)
838    std::cout <<"Found NLP-integer feasible solution of  value : "<<solValue<<std::endl;
839  std::cout<<"Current MILP lower bound is : "<<milpBound<<std::endl;
840  std::cout<<"-------------------------------------------------------------------"<<std::endl;
841  std::cout<<"Stopped because : isInteger "<<isInteger<<", feasible "<<feasible<<std::endl<<std::endl;
842
843}
844
845
846
847#endif
848}/* End namespace Bonmin. */
849
Note: See TracBrowser for help on using the repository browser.