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

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

Change CPLEX code for Hassan

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