source: branches/devel/Bonmin/src/OaInterface/BonOACutGenerator2.cpp @ 57

Last change on this file since 57 was 57, checked in by pbonami, 13 years ago

Rename file with Bon prefix part II

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 27.1 KB
Line 
1// (C) Copyright Carnegie Mellon University 2005
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// Authors :
6// P. Bonami, Carnegie Mellon University
7//
8// Date : 05/26/2005
9//#define OA_DEBUG
10
11#include "BonminConfig.h"
12
13#include "BonOACutGenerator2.hpp"
14#include "OsiClpSolverInterface.hpp"
15
16#include "CbcModel.hpp"
17#include "CbcStrategy.hpp"
18#ifdef COIN_HAS_CPX
19#include "OsiCpxSolverInterface.hpp"
20#endif
21#include "OsiAuxInfo.hpp"
22
23
24extern CbcModel * OAModel;
25namespace Bonmin{
26/// Default constructor
27OACutGenerator2::OACutGenerator2():
28    CglCutGenerator(),
29    nlp_(NULL),
30    nSolve_(0),
31    si_(NULL),
32    cbcCutoffIncrement_(1e-06),
33    cbcIntegerTolerance_(1e-05),
34    localSearchNodeLimit_(0),
35    maxLocalSearchPerNode_(0),
36    maxLocalSearch_(0),
37    maxLocalSearchTime_(3600),
38    nLocalSearch_(0),
39    solveAuxiliaryProblem_(0),
40    handler_(NULL),
41    subMilpLogLevel_(0),
42    leaveSiUnchanged_(0),
43    strategy_(NULL),
44    timeBegin_(0.),
45    logFrequency_(1000.)
46{
47  handler_ = new CoinMessageHandler();
48  handler_ -> setLogLevel(2);
49  messages_ = OaMessages();
50  timeBegin_ = CoinCpuTime();
51}
52
53
54
55OACutGenerator2::OACutGenerator2
56(OsiTMINLPInterface * nlp,
57 OsiSolverInterface * si,
58 CbcStrategy * strategy,
59 double cbcCutoffIncrement,
60 double cbcIntegerTolerance,
61 bool solveAuxiliaryProblem,
62 bool leaveSiUnchanged
63)
64    :
65    CglCutGenerator(),
66    nlp_(nlp),
67    nSolve_(0),
68    si_(si),
69    cbcCutoffIncrement_(cbcCutoffIncrement),
70    cbcIntegerTolerance_(cbcIntegerTolerance),
71    localSearchNodeLimit_(0),
72    maxLocalSearchPerNode_(0),
73    maxLocalSearch_(0),
74    maxLocalSearchTime_(3600),
75    nLocalSearch_(0),
76    solveAuxiliaryProblem_(solveAuxiliaryProblem),
77    handler_(NULL),
78    subMilpLogLevel_(0),
79    leaveSiUnchanged_(leaveSiUnchanged),
80    strategy_(NULL),
81    timeBegin_(0),
82    logFrequency_(1000.)
83{
84  handler_ = new CoinMessageHandler();
85  handler_ -> setLogLevel(2);
86  messages_ = OaMessages();
87  if(strategy)
88    strategy_ = strategy->clone();
89  timeBegin_ = CoinCpuTime();
90}
91
92OACutGenerator2::~OACutGenerator2()
93{
94  delete handler_;
95  if(strategy_)
96    delete strategy_;
97}
98/// Assign an OsiTMINLPInterface
99void
100OACutGenerator2::assignNlpInterface(OsiTMINLPInterface * nlp)
101{
102  nlp_ = nlp;
103}
104
105/// Assign an OsiTMINLPInterface
106void
107OACutGenerator2::assignLpInterface(OsiSolverInterface * si)
108{
109  si_ = si;
110  if(maxLocalSearch_>0) {
111    setTheNodeLimit();
112  }
113}
114
115double OACutGenerator2::siBestObj(CbcModel * model) const
116{
117  if(model == NULL) {
118    //Check solver name to see if local searches can be performed
119    std::string solverName;
120    si_->getStrParam(OsiSolverName,solverName);
121    //for cbc needs to get the first three letter of solver name
122    std::string shortSolverName(solverName,0,3);
123    if(shortSolverName == "cbc") {
124      throw CoinError("OsiCbc is not supported for doing local searches use OsiClpSolverInterface instead",
125          "OACutGenerator2","setTheNodeLimit");
126    }
127    else if(solverName == "cplex") {
128#ifdef COIN_HAS_CPX
129      OsiCpxSolverInterface * cpx = dynamic_cast<OsiCpxSolverInterface *>(si_);
130      double value;
131      int status = CPXgetbestobjval(cpx->getEnvironmentPtr(),cpx->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL), &value);
132      if(status)
133        throw CoinError("Error in getting CPLEX best bound","OACutGenerator2","siBestObj");
134      return value;
135#else
136
137      throw CoinError("You have to define COIN_HAS_CPX at compilation to be able to use cplex for local searches",
138          "OACutGenerator2","siBestObj");
139#endif
140
141    }
142    else {
143      std::string mesg="Unsuported solver";
144      mesg += solverName;
145      mesg +=" for local searches, you should use Cbc or Cplex";
146      throw CoinError(mesg,
147          "OACutGenerator2","assignLpInterface");
148    }
149  }
150  else
151    return model->getBestPossibleObjValue();
152}
153void OACutGenerator2::setTheNodeLimit()
154{
155  //Check solver name to see if local searches can be performed
156  std::string solverName;
157  si_->getStrParam(OsiSolverName,solverName);
158
159  //for cbc needs to get the first three letter of solver name
160  std::string shortSolverName(solverName,0,3);
161  if(shortSolverName == "cbc") {
162    throw CoinError("OsiCbc is not supported for doing local searches use OsiClpSolverInterface instead",
163        "OACutGenerator2","setTheNodeLimit");
164  }
165  else if(solverName == "cplex") {
166#ifdef COIN_HAS_CPX
167    OsiCpxSolverInterface * cpx = dynamic_cast<OsiCpxSolverInterface *>(si_);
168    CPXsetintparam(cpx->getEnvironmentPtr(),CPX_PARAM_NODELIM, localSearchNodeLimit_);
169#else
170
171    throw CoinError("You have to define COIN_HAS_CPX at compilation to be able to use cplex for local searches",
172        "OACutGenerator2","setTheNodeLimit");
173#endif
174
175  }
176  else if(solverName == "clp") {
177    //do nothing will set up node limit when creating a CbcModel
178  }
179  else {
180    std::string mesg="Unsuported solver";
181    mesg += solverName;
182    mesg +=" for local searches, you should use Cbc or Cplex";
183    throw CoinError(mesg,
184        "OACutGenerator2","setTheNodeLimit");
185  }
186}
187
188
189void OACutGenerator2::setTimeLimit(double time) const
190{
191  //Check solver name to see if local searches can be performed
192  std::string solverName;
193  si_->getStrParam(OsiSolverName,solverName);
194  //for cbc needs to get the first three letter of solver name
195  std::string shortSolverName(solverName,0,3);
196  if(shortSolverName == "cbc") {
197    throw CoinError("OsiCbc is not supported for doing local searches use OsiClpSolverInterface instead",
198        "OACutGenerator2","setTheNodeLimit");
199  }
200  else if(solverName == "cplex") {
201#ifdef COIN_HAS_CPX
202    OsiCpxSolverInterface * cpx = dynamic_cast<OsiCpxSolverInterface *>(si_);
203    CPXsetdblparam(cpx->getEnvironmentPtr(),CPX_PARAM_TILIM, time);
204#else
205
206    throw CoinError("You have to define COIN_HAS_CPX at compilation to be able to use cplex for local searches",
207        "OACutGenerator2","setTheNodeLimit");
208#endif
209
210  }
211  else if(solverName == "clp") {
212    //do nothing will set up node limit when creating a CbcModel
213  }
214  else {
215    std::string mesg="Unsuported solver";
216    mesg += solverName;
217    mesg +=" for local searches, you should use Cbc or Cplex";
218    throw CoinError(mesg,
219        "OACutGenerator2","setTheNodeLimit");
220  }
221}
222
223void OACutGenerator2::setCutoff(double bestKnown) const
224{
225  //Check solver name to see if local searches can be performed
226  std::string solverName;
227  si_->getStrParam(OsiSolverName,solverName);
228
229  //for cbc needs to get the first three letter of solver name
230  std::string shortSolverName(solverName,0,3);
231  if(shortSolverName == "cbc") {
232    throw CoinError("OsiCbc is not supported for doing local searches use OsiClpSolverInterface instead",
233        "OACutGenerator2","setTheNodeLimit");
234  }
235  else if(solverName == "cplex") {
236#ifdef COIN_HAS_CPX
237    OsiCpxSolverInterface * cpx = dynamic_cast<OsiCpxSolverInterface *>(si_);
238    CPXsetdblparam(cpx->getEnvironmentPtr(),CPX_PARAM_CUTUP, bestKnown);
239#else
240
241    throw CoinError("You have to define COIN_HAS_CPX at compilation to be able to use cplex for local searches",
242        "OACutGenerator2","setTheNodeLimit");
243#endif
244
245  }
246  else if(solverName == "clp") {
247    //do nothing will set up node limit when creating a CbcModel
248  }
249  else {
250    std::string mesg="Unsuported solver";
251    mesg += solverName;
252    mesg +=" for local searches, you should use Cbc or Cplex";
253    throw CoinError(mesg,
254        "OACutGenerator2","setTheNodeLimit");
255  }
256}/// cut generation method
257void
258OACutGenerator2::generateCuts( const OsiSolverInterface & si, OsiCuts & cs,
259    const CglTreeInfo info) const
260{
261  CbcStrategy * strategy = strategy_;
262  double lastPeriodicLog= CoinCpuTime();
263
264  if(nlp_ == NULL) {
265    std::cerr<<"Error in cut generator for outer approximation no NLP ipopt assigned"<<std::endl;
266    throw -1;
267  }
268
269  //This node may already be fathomed by a previous call to this
270  // this information is stored in babInfo
271  OsiBabSolver * babInfo = dynamic_cast<OsiBabSolver *> (si.getAuxiliaryInfo());
272#if 1
273  //if algorithms would converge this should never happen and it seems very dangerous if somebody forgot to reset the bound
274  if(babInfo)
275    if(!babInfo->mipFeasible())
276      return;
277
278#endif
279
280  const int numcols = nlp_->getNumCols();
281
282  //Get the continuous solution
283  const double *colsol = si.getColSolution();
284
285  bool doLocalSearch = 0;
286  bool isInteger = 1;
287  //   nlp_->setMilpBound(-1e100);
288  //   nlp_->setMilpFeasible(true);
289
290
291  //Check integer infeasibility
292  for(int i = 0 ; i < numcols ; i++) {
293    if(nlp_->isInteger(i)) {
294      if(fabs(colsol[i] - floor(colsol[i] + 0.5) ) >
295          cbcIntegerTolerance_) {
296        isInteger = 0;
297        break;
298      }
299    }
300  }
301
302  if(!isInteger) {
303    if(nLocalSearch_<maxLocalSearch_ &&
304        localSearchNodeLimit_ > 0 &&
305        CoinCpuTime() - timeBegin_ < maxLocalSearchTime_)//do a local search
306    {
307      doLocalSearch = 1;
308    }
309    else {
310      //          nlp_->setMilpBound(-1e100);
311      //          nlp_->setMilpFeasible(true);
312      if(strategy && ! strategy_)
313        delete strategy;
314      return;
315    }
316  }
317
318
319  //We are going to generate some cuts copy LP information
320
321
322  //get the current cutoff
323  double cutoff;
324  si.getDblParam(OsiDualObjectiveLimit, cutoff);
325  double saveCutoff=cutoff;
326
327  // Save column bounds and basis
328  double * saveColLb = new double[numcols];
329  double * saveColUb = new double[numcols];
330  CoinCopyN(nlp_->getColLower(), numcols , saveColLb);
331  CoinCopyN(nlp_->getColUpper(), numcols , saveColUb);
332  int originalRowNumber = (si_!=NULL) ?si_->getNumRows() : -1;
333  CoinWarmStart * saveWarmStart = NULL;
334  bool milpOptimal = 1;
335
336#ifdef NO_NULL_SI
337  if(si_ == NULL) {
338    std::cerr<<"Error in cut generator for outer approximation no lp solver interface assigned"<<std::endl;
339    throw -1;
340  }
341#else
342
343  bool deleteSi = false;
344  if(si_ == NULL) {
345    si_ = si.clone();
346    deleteSi = true;
347  }
348  else
349#endif
350  if(si_!=&si)//We may have the same lp as the one passed by the model or a local copy
351  {
352
353    //Install current active cuts into local solver
354    int numberCutsToAdd = si.getNumRows();
355    numberCutsToAdd -= si_->getNumRows();
356    if(numberCutsToAdd > 0)//Have to install some cuts
357    {
358      CoinPackedVector * * addCuts = new CoinPackedVector *[numberCutsToAdd];
359      for(int i = 0 ; i < numberCutsToAdd ; i++)
360      {
361        addCuts[i] = new CoinPackedVector;
362      }
363      //Get the current matrix and fill the addCuts
364      const CoinPackedMatrix * mat = si.getMatrixByCol();
365      const CoinBigIndex * start = mat->getVectorStarts();
366      const int * length = mat->getVectorLengths();
367      const double * elements = mat->getElements();
368      const int * indices = mat->getIndices();
369      for(int i = 0 ; i <= numcols ; i++)
370        for(int k = start[i] ; k < start[i] + length[i] ; k++)
371        {
372          if(indices[k] >= si_->getNumRows()) {
373            addCuts[ indices[k] - si_->getNumRows() ]->insert(i, elements[k]);
374          }
375        }
376      si_->addRows(numberCutsToAdd, (const CoinPackedVectorBase * const *) addCuts, &si.getRowLower()[si_->getNumRows()],
377          &si.getRowUpper()[si_->getNumRows()]);
378    }
379    else if (numberCutsToAdd < 0)//Oups some error
380    {
381      std::cerr<<"Internal error in OACutGenerator2 : number of cuts wrong"<<std::endl;
382    }
383
384    //Set the bounds on columns
385    for(int i = 0 ; i <= numcols ; i++)
386    {
387      si_->setColBounds(i, si.getColLower()[i], si.getColUpper()[i]);
388    }
389    //Install basis in problem
390    CoinWarmStart * warm = si.getWarmStart();
391    if(si_->setWarmStart(warm)==false)
392    {
393      delete warm;
394      throw CoinError("Fail installing the warm start in the subproblem",
395          "generateCuts","OACutGenerator2") ;
396    }
397    delete warm;
398    //put the cutoff
399    si_->setDblParam(OsiDualObjectiveLimit, cutoff);
400    si_->resolve();
401
402#ifdef OA_DEBUG
403
404    std::cout<<"Resolve with hotstart :"<<std::endl
405    <<"number of iterations(should be 0) : "<<si_->getIterationCount()<<std::endl
406    <<"Objective value and diff to original : "<<si_->getObjValue()<<", "
407    <<fabs(si_->getObjValue() - si.getObjValue())<<std::endl;
408    for(int i = 0 ; i <= numcols ; i++)
409    {
410      if(fabs(si.getColSolution()[i]-si_->getColSolution()[i])>1e-08) {
411        std::cout<<"Diff between solution at node and solution with local solver : "<<fabs(si.getColSolution()[i]-si_->getColSolution()[i])<<std::endl;
412      }
413    }
414#endif
415  }
416  else {
417#ifdef NO_NULL_SI
418    throw CoinError("Not allowed to modify si in a cutGenerator",
419        "OACutGenerator2","generateCuts");
420#else //Seems that nobody wants to allow me to do this
421    if(leaveSiUnchanged_)
422      saveWarmStart = si.getWarmStart();
423#endif
424  }
425
426  OsiClpSolverInterface * clp = dynamic_cast<OsiClpSolverInterface *>(si_);
427  CbcModel * model = NULL;//We will need the model through all the function
428  double milpBound = -DBL_MAX;
429  bool milpFeasible = 1;
430  //double milpBound=si_->getObjValue();
431  bool feasible = 1;
432
433  if(feasible && doLocalSearch)//Perform a local search
434  {
435    if(clp)//Clp does not do branch-and-bound have to create a CbcModel
436    {
437      OsiBabSolver empty;
438      model = new CbcModel(*clp); // which clones
439      OAModel = model;
440      model->solver()->setAuxiliaryInfo(&empty);
441
442      //Change Cbc messages prefixes
443      strcpy(model->messagesPointer()->source_,"OaCbc");
444
445      if(!strategy)
446        strategy = new CbcStrategyDefault(1,0,0,subMilpLogLevel_);
447
448      clp->resolve();
449      model->setLogLevel(subMilpLogLevel_);
450      model->solver()->messageHandler()->setLogLevel(0);
451      model->setStrategy(*strategy);
452      model->setMaximumNodes(localSearchNodeLimit_);
453      model->setMaximumSeconds(maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
454      model->setCutoff(cutoff);
455      model->branchAndBound();
456      milpBound = siBestObj(model);
457      if(model->isProvenOptimal())
458      {
459        milpOptimal = 1;
460      }
461      else
462        milpOptimal = 0;
463      feasible = milpBound < cutoff;
464      milpFeasible = feasible;
465      isInteger = model->getSolutionCount();
466      nLocalSearch_++;
467
468      if(milpOptimal)
469        handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<model->getNodeCount()<<model->getIterationCount()<<CoinMessageEol;
470      else
471      {
472        handler_->message(LOCAL_SEARCH_ABORT, messages_)<<model->getNodeCount()<<model->getIterationCount()<<CoinMessageEol;
473        if(isInteger)
474          std::cout<<"Integer feasible solution found"<<std::endl;
475      }
476    }
477    else//use OsiSolverInterface::branchAndBound
478    {
479      setTimeLimit(maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
480      si_->messageHandler()->setLogLevel(subMilpLogLevel_);
481      si_->branchAndBound();
482      nLocalSearch_++;
483      //Did we find any integer feasible solution
484      isInteger = si_->getFractionalIndices().size() == 0;
485      milpBound = siBestObj();
486      feasible = milpBound < cutoff;
487      milpFeasible = feasible;
488    }
489  }
490  int numberPasses = 0;
491  bool foundSolution = 0;
492  while(isInteger && feasible ) {
493    numberPasses++;
494
495    //eventually give some information to user
496    double time = CoinCpuTime();
497    if(time - lastPeriodicLog > logFrequency_) {
498      double lb = (model == NULL) ?si_->getObjValue():model->getBestPossibleObjValue();
499      handler_->message(PERIODIC_MSG,messages_)
500      <<time - timeBegin_<<cutoff
501      <<lb
502      <<CoinMessageEol;
503      lastPeriodicLog = CoinCpuTime();
504    }
505
506
507    //setup the nlp
508    int numberCutsBefore = cs.sizeRowCuts();
509
510    bool fixed = true;
511
512#ifdef OA_DEBUG
513  std::cout<<"FX_P   : ";
514#endif
515  //Fix the variable which have to be fixed, after having saved the bounds
516  for(int i = 0; i < numcols ; i++) {
517      if(nlp_->isInteger(i)) {
518        double value =  (model == NULL) ? si_->getColSolution()[i] : model->bestSolution()[i];
519        if(saveColUb[i] < saveColLb[i] + cbcIntegerTolerance_)
520          fixed = false;
521      value = floor(value+0.5);
522        value = max(saveColLb[i],value);
523        value = min(value, saveColUb[i]);
524        if(fabs(value) > 1e10) { std::cerr<<"FATAL ERROR: Variable taking a big value ("<<value
525                                 <<") at optimium of LP relaxation, can not construct outer approximation. You should try running the problem with B-BB"<<std::endl;
526         throw -1;
527        }
528#ifdef OA_DEBUG
529        //         printf("xx %d at %g (bounds %g, %g)",i,value,nlp_->getColLower()[i],
530        //                nlp_->getColUpper()[i]);
531        std::cout<<(int)value;
532#endif
533        nlp_->setColLower(i,value);
534        nlp_->setColUpper(i,value);
535      }
536    }
537#ifdef OA_DEBUG
538    std::cout<<std::endl;
539#endif
540    //Now solve the NLP get the cuts, and intall them in the local LP
541
542    //  nlp_->turnOnIpoptOutput();
543    nSolve_++;
544    nlp_->resolve();
545    if(nlp_->isProvenOptimal()) {
546      nlp_->getOuterApproximation(cs);
547      handler_->message(FEASIBLE_NLP, messages_)
548      <<nlp_->getIterationCount()
549      <<nlp_->getObjValue()<<CoinMessageEol;
550
551
552#ifdef OA_DEBUG
553
554      const double * colsol2 = nlp_->getColSolution();
555      for(int i = 0 ; i < numcols; i++) {
556        if(nlp_->isInteger(i)) {
557          if(fabs(colsol2[i] - floor(colsol2[i] + 0.5) ) >
558              cbcIntegerTolerance_) {
559            std::cerr<<"Integer infeasible point (should not be), integer infeasibility for variable "<<i
560            <<" is, "<<fabs(colsol2[i] - floor(colsol2[i] + 0.5))<<std::endl;
561            //              integerFeasible = 0;
562          }
563        }
564      }
565#endif
566
567      if((nlp_->getObjValue() < cutoff) ) {
568        handler_->message(UPDATE_UB, messages_)
569        <<nlp_->getObjValue()
570        <<CoinCpuTime()-timeBegin_
571        <<CoinMessageEol;
572
573        foundSolution = 1;
574#if 1
575        // Also store into solver
576        if(babInfo) {
577          double * lpSolution = new double[numcols + 1];
578          CoinCopyN(nlp_->getColSolution(), numcols, lpSolution);
579          lpSolution[numcols] = nlp_->getObjValue();
580          babInfo->setSolution(lpSolution,
581              numcols + 1, lpSolution[numcols]);
582          delete [] lpSolution;
583        }
584        else {
585          printf("No auxiliary info in nlp solve!\n");
586          throw -1;
587        }
588#endif
589        // Update the cutoff
590        cutoff = nlp_->getObjValue() *(1 - cbcCutoffIncrement_);
591        // Update the lp solver cutoff
592        si_->setDblParam(OsiDualObjectiveLimit, cutoff);
593      }
594    }
595    else if(nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
596      std::cerr<<"Unsolved NLP... exit"<<std::endl;
597    }
598    else {
599      handler_->message(INFEASIBLE_NLP, messages_)
600      <<nlp_->getIterationCount()
601      <<CoinMessageEol;
602      if(solveAuxiliaryProblem_)  //solve feasibility NLP and add corresponding outer approximation constraints
603      {
604        double *x = new double[numcols];
605        int * ind = new int[numcols];
606        int nInd = 0;
607        for(int i = 0; i < numcols ; i++)
608        {
609          if(nlp_->isInteger(i)) {
610            ind[nInd] = i;
611            x[nInd++] = (model == NULL) ? si_->getColSolution()[i] : model->bestSolution()[i];
612            //reset the bounds
613            nlp_->setColBounds(i,saveColLb[i],saveColUb[i]);
614          }
615        }
616        nlp_->getFeasibilityOuterApproximation(nInd, x, ind,cs);
617        if(nlp_->isProvenOptimal())
618        {
619          assert((nlp_->getObjValue() > cbcIntegerTolerance_) );
620          std::cout<<"Solved auxiliary infeasibility problem"<<std::endl;
621        }
622        else
623        {
624          std::cout<<"Failed to solve auxiliary infeasibility problem"<<std::endl;
625          bool * s = const_cast<bool *> (&solveAuxiliaryProblem_);
626          *s =0;
627        }
628        delete []x;
629        delete[]ind;
630      }
631      else
632        nlp_->getOuterApproximation(cs);
633    }
634
635    //install the cuts in si_ reoptimize and check for integer feasibility
636    // Code taken and adapted from CbcModel::solveWithCuts
637    int numberCuts = cs.sizeRowCuts() - numberCutsBefore;
638    if (numberCuts > 0) {
639      CoinWarmStartBasis * basis
640      = dynamic_cast<CoinWarmStartBasis*>(si_->getWarmStart()) ;
641      assert(basis != NULL); // make sure not volume
642      basis->resize(si_->getNumRows()+numberCuts,numcols + 1) ;
643      for (int i = 0 ; i < numberCuts ; i++) {
644        basis->setArtifStatus(si_->getNumRows()+i,
645            CoinWarmStartBasis::basic) ;
646      }
647
648      const OsiRowCut ** addCuts = new const OsiRowCut * [numberCuts] ;
649      for (int i = 0 ; i < numberCuts ; i++) {
650        addCuts[i] = &cs.rowCut(i + numberCutsBefore) ;
651      }
652      si_->applyRowCuts(numberCuts,addCuts) ;
653
654      delete [] addCuts ;
655      if (si_->setWarmStart(basis) == false) {
656        delete basis;
657        throw CoinError("Fail setWarmStart() after cut installation.",
658            "generateCuts","OACutGenerator2") ;
659      }
660      delete basis;
661      si_->resolve();
662      double objvalue = si_->getObjValue();
663      //milpBound = max(milpBound, si_->getObjValue());
664      feasible = (si_->isProvenOptimal() &&
665          !si_->isDualObjectiveLimitReached() && (objvalue<cutoff)) ;
666      //if value of integers are unchanged then we have to get out
667      bool changed = //!nlp_->isProvenOptimal()//if nlp_ is infeasible solution has necessarily changed
668          !feasible;//if lp is infeasible we don't have to check anything
669      const double * nlpSol = nlp_->getColSolution();
670      const double * lpSol = si_->getColSolution();
671      for(int i = 0; i < numcols && !changed; i++) {
672        if(nlp_->isInteger(i) && fabs(nlpSol[i] - lpSol[i])>0.001)
673          changed = 1;
674      }
675      if(changed) {
676        isInteger = 1;
677        for(int i = 0 ; i < numcols ; i++) {
678          if(nlp_->isInteger(i)) {
679            if(fabs(lpSol[i] - floor(lpSol[i] + 0.5) ) > cbcIntegerTolerance_ ) {
680              isInteger = 0;
681              break;
682            }
683          }
684        }
685      }
686      else {
687        isInteger = 0;
688        //        if(!fixed)//fathom on bounds
689        milpBound = 1e200;
690      }
691#ifdef OA_DEBUG
692
693      printf("Obj value after cuts %g %d rows\n",si_->getObjValue(),
694          numberCuts) ;
695#endif
696      //do we perform a new local search ?
697      if(milpOptimal && feasible && !isInteger &&
698          nLocalSearch_ < maxLocalSearch_ &&
699          numberPasses < maxLocalSearchPerNode_ &&
700          localSearchNodeLimit_ > 0 &&
701          CoinCpuTime() - timeBegin_ < maxLocalSearchTime_) {
702         std::cout<<"Perform new local search"<<std::endl;
703        if(clp==NULL) {
704          nLocalSearch_++;
705
706          setTimeLimit(maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
707          setCutoff(cutoff);
708          si_->branchAndBound();
709          //Did we find any integer feasible solution
710          milpBound = siBestObj();
711          //If integer solution is the same as nlp solution problem is solved
712          changed = !nlp_->isProvenOptimal();
713          for(int i = 0; i < numcols && !changed; i++) {
714            if(nlp_->isInteger(i) && fabs(nlp_->getColSolution()[i] - si_->getColSolution()[i])>cbcIntegerTolerance_)
715              changed = 1;
716          }
717          if(changed) {
718            feasible = (milpBound < cutoff);
719            std::cout<<"milp bound "<<milpBound<<" cutoff "<<cutoff<<std::endl;
720          }
721          else
722           {
723            feasible = 0;
724            milpBound = 1e50;
725           }
726          isInteger = si_->getFractionalIndices().size() == 0;
727        }/** endif solved by branchAndBound capable OsiInterface*/
728        else {
729          if(model)
730            {
731              OAModel = NULL;       
732              delete model;
733            }
734          model = new CbcModel(*clp); // which clones
735          OAModel = model;
736          OsiBabSolver empty;
737          model->solver()->setAuxiliaryInfo(&empty);
738          //Change Cbc messages prefixes
739          strcpy(model->messagesPointer()->source_,"OaCbc");
740          model->solver()->messageHandler()->setLogLevel(0);
741
742          if(!strategy)
743            strategy = new CbcStrategyDefault(1,0,0, subMilpLogLevel_);
744          model->setLogLevel(subMilpLogLevel_);
745          model->setStrategy(*strategy);
746          model->setMaximumNodes(localSearchNodeLimit_);
747          model->setMaximumSeconds(maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
748          model->setCutoff(cutoff);
749          model->branchAndBound();
750          nLocalSearch_++;
751          milpBound = siBestObj(model);
752          handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<model->getNodeCount()<<model->getIterationCount()<<CoinMessageEol;
753          feasible =  (milpBound < cutoff);
754          isInteger = model->getSolutionCount();
755          if(model->isProvenOptimal() || model->isProvenInfeasible()) {
756            bool changed = 0;             //If integer solution is the same as nlp solution problem is solved
757            for(int i = 0; feasible && isInteger && i < numcols && !changed; i++) {
758              if(nlp_->isInteger(i) && fabs(nlp_->getColSolution()[i] - model->bestSolution()[i])>cbcIntegerTolerance_)
759                changed = 1;
760            }
761          if(!changed) {
762            feasible = 0;
763            milpBound = 1e50;
764           }
765            milpFeasible = feasible;
766          }
767          else {
768            std::cout<<"Exiting on time limit"<<std::endl;
769            milpOptimal = 0;
770            //feasible = 1;
771            break;
772          }
773
774        }/** endif solved by clp/cbc*/
775
776        if(milpBound < cutoff)
777          handler_->message(UPDATE_LB, messages_)
778          <<milpBound<<CoinCpuTime() - timeBegin_
779          <<CoinMessageEol;
780        else
781        {
782          milpBound = 1e50;
783          feasible = 0;
784          handler_->message(OASUCCESS, messages_)
785          <<CoinCpuTime() - timeBegin_ <<CoinMessageEol;
786        }
787      }/** endif localSearch*/
788      else if(model!=NULL)/** have to delete model for next iteration*/
789      { 
790        delete model;
791        OAModel = NULL;
792        model=NULL;
793      }
794
795    }
796
797
798
799  }
800#ifdef OA_DEBUG
801    std::cout<<"-------------------------------------------------------------------------------------------------------------------------"<<std::endl;
802    std::cout<<"OA cut generation finished"<<std::endl;
803    std::cout<<"Generated "<<cs.sizeRowCuts()<<std::endl;
804    if(foundSolution)
805      std::cout <<"Found NLP-integer feasible solution of  value : "<<cutoff<<std::endl;
806    std::cout<<"Current MILP lower bound is : "<<milpBound<<std::endl;
807    std::cout<<"-------------------------------------------------------------------------------------------------------------------------"<<std::endl;
808    std::cout<<"Stopped because : isInteger "<<isInteger<<", feasible "<<feasible<<std::endl<<std::endl;
809#endif
810    //Transmit the bound found by the milp
811    {
812      if(milpBound>-1e100)
813      {
814#if 1
815        // Also store into solver
816        if (babInfo)
817          babInfo->setMipBound(milpBound);
818#endif
819
820      }
821    }  //Clean everything :
822
823  //free model
824  if(model!= NULL) {
825    delete model;
826  }
827
828
829
830  //  Reset bounds in the NLP
831
832  for(int i = 0; i < numcols ; i++) {
833    if(nlp_->isInteger(i)) {
834      nlp_->setColBounds(i,saveColLb[i],saveColUb[i]);
835    }
836  }
837#ifndef NO_NULL_SI
838  if(deleteSi) {
839    delete si_;
840    si_ = NULL;
841  }
842  else
843#endif
844   {
845    if(leaveSiUnchanged_) {
846      int nRowsToDelete = si_->getNumRows() - originalRowNumber;
847      int * rowsToDelete = new int[nRowsToDelete];
848      for(int i = 0 ; i < nRowsToDelete ; i++) {
849        rowsToDelete[i] = i + originalRowNumber;
850      }
851      si_->deleteRows(nRowsToDelete, rowsToDelete);
852      delete [] rowsToDelete;
853      if(si_==&si) {
854        si_->setDblParam(OsiDualObjectiveLimit, saveCutoff);
855        if(si_->setWarmStart(saveWarmStart)==false) {
856          throw CoinError("Fail restoring the warm start at the end of procedure",
857              "generateCuts","OACutGenerator2") ;
858        }
859        delete saveWarmStart;
860      }
861    }
862}
863
864delete [] saveColLb;
865delete [] saveColUb;
866if(strategy && ! strategy_)
867  delete strategy;
868}
869}
Note: See TracBrowser for help on using the repository browser.