source: trunk/Bonmin/src/OaInterface/IpCbcOACutGenerator2.cpp @ 1

Last change on this file since 1 was 1, checked in by andreasw, 13 years ago

imported initial code

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 26.9 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 "IpCbcOACutGenerator2.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;
25
26/// Default constructor
27IpCbcOACutGenerator2::IpCbcOACutGenerator2():
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
55IpCbcOACutGenerator2::IpCbcOACutGenerator2
56(IpoptInterface * 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
92IpCbcOACutGenerator2::~IpCbcOACutGenerator2()
93{
94  delete handler_;
95  if(strategy_)
96    delete strategy_;
97}
98/// Assign an IpoptInterface
99void
100IpCbcOACutGenerator2::assignNlpInterface(IpoptInterface * nlp)
101{
102  nlp_ = nlp;
103}
104
105/// Assign an IpoptInterface
106void
107IpCbcOACutGenerator2::assignLpInterface(OsiSolverInterface * si)
108{
109  si_ = si;
110  if(maxLocalSearch_>0) {
111    setTheNodeLimit();
112  }
113}
114
115double IpCbcOACutGenerator2::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          "IpCbcOACutGenerator2","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","IpCbcOACutGenerator2","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          "IpCbcOACutGenerator2","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          "IpCbcOACutGenerator2","assignLpInterface");
148    }
149  }
150  else
151    return model->getBestPossibleObjValue();
152}
153void IpCbcOACutGenerator2::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        "IpCbcOACutGenerator2","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        "IpCbcOACutGenerator2","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        "IpCbcOACutGenerator2","setTheNodeLimit");
185  }
186}
187
188
189void IpCbcOACutGenerator2::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        "IpCbcOACutGenerator2","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        "IpCbcOACutGenerator2","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        "IpCbcOACutGenerator2","setTheNodeLimit");
220  }
221}
222
223void IpCbcOACutGenerator2::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        "IpCbcOACutGenerator2","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        "IpCbcOACutGenerator2","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        "IpCbcOACutGenerator2","setTheNodeLimit");
255  }
256}/// cut generation method
257void
258IpCbcOACutGenerator2::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 IpCbcOACutGenerator2 : 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","IpCbcOACutGenerator2") ;
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        "IpCbcOACutGenerator2","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#ifdef OA_DEBUG
525        //         printf("xx %d at %g (bounds %g, %g)",i,value,nlp_->getColLower()[i],
526        //                nlp_->getColUpper()[i]);
527        std::cout<<(int)value;
528#endif
529        nlp_->setColLower(i,value);
530        nlp_->setColUpper(i,value);
531      }
532    }
533#ifdef OA_DEBUG
534    std::cout<<std::endl;
535#endif
536    //Now solve the NLP get the cuts, and intall them in the local LP
537
538    //  nlp_->turnOnIpoptOutput();
539    nSolve_++;
540    nlp_->resolve();
541    if(nlp_->isProvenOptimal()) {
542      nlp_->getOuterApproximation(cs);
543      handler_->message(FEASIBLE_NLP, messages_)
544      <<nlp_->getIterationCount()
545      <<nlp_->getObjValue()<<CoinMessageEol;
546
547
548#ifdef OA_DEBUG
549
550      const double * colsol2 = nlp_->getColSolution();
551      for(int i = 0 ; i < numcols; i++) {
552        if(nlp_->isInteger(i)) {
553          if(fabs(colsol2[i] - floor(colsol2[i] + 0.5) ) >
554              cbcIntegerTolerance_) {
555            std::cerr<<"Integer infeasible point (should not be), integer infeasibility for variable "<<i
556            <<" is, "<<fabs(colsol2[i] - floor(colsol2[i] + 0.5))<<std::endl;
557            //              integerFeasible = 0;
558          }
559        }
560      }
561#endif
562
563      if((nlp_->getObjValue() < cutoff) ) {
564        handler_->message(UPDATE_UB, messages_)
565        <<nlp_->getObjValue()
566        <<CoinCpuTime()-timeBegin_
567        <<CoinMessageEol;
568
569        foundSolution = 1;
570#if 1
571        // Also store into solver
572        if(babInfo) {
573          double * lpSolution = new double[numcols + 1];
574          CoinCopyN(nlp_->getColSolution(), numcols, lpSolution);
575          lpSolution[numcols] = nlp_->getObjValue();
576          babInfo->setSolution(lpSolution,
577              numcols + 1, lpSolution[numcols]);
578          delete [] lpSolution;
579        }
580        else {
581          printf("No auxiliary info in nlp solve!\n");
582          throw -1;
583        }
584#endif
585        // Update the cutoff
586        cutoff = nlp_->getObjValue() *(1 - cbcCutoffIncrement_);
587        // Update the lp solver cutoff
588        si_->setDblParam(OsiDualObjectiveLimit, cutoff);
589      }
590    }
591    else if(nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
592      std::cerr<<"Unsolved NLP... exit"<<std::endl;
593    }
594    else {
595      handler_->message(INFEASIBLE_NLP, messages_)
596      <<nlp_->getIterationCount()
597      <<CoinMessageEol;
598      if(solveAuxiliaryProblem_)  //solve feasibility NLP and add corresponding outer approximation constraints
599      {
600        double *x = new double[numcols];
601        int * ind = new int[numcols];
602        int nInd = 0;
603        for(int i = 0; i < numcols ; i++)
604        {
605          if(nlp_->isInteger(i)) {
606            ind[nInd] = i;
607            x[nInd++] = (model == NULL) ? si_->getColSolution()[i] : model->bestSolution()[i];
608            //reset the bounds
609            nlp_->setColBounds(i,saveColLb[i],saveColUb[i]);
610          }
611        }
612        nlp_->getFeasibilityOuterApproximation(nInd, x, ind,cs);
613        if(nlp_->isProvenOptimal())
614        {
615          assert((nlp_->getObjValue() > cbcIntegerTolerance_) );
616          std::cout<<"Solved auxiliary infeasibility problem"<<std::endl;
617        }
618        else
619        {
620          std::cout<<"Failed to solve auxiliary infeasibility problem"<<std::endl;
621          bool * s = const_cast<bool *> (&solveAuxiliaryProblem_);
622          *s =0;
623        }
624        delete []x;
625        delete[]ind;
626      }
627      else
628        nlp_->getOuterApproximation(cs);
629    }
630
631    //install the cuts in si_ reoptimize and check for integer feasibility
632    // Code taken and adapted from CbcModel::solveWithCuts
633    int numberCuts = cs.sizeRowCuts() - numberCutsBefore;
634    if (numberCuts > 0) {
635      CoinWarmStartBasis * basis
636      = dynamic_cast<CoinWarmStartBasis*>(si_->getWarmStart()) ;
637      assert(basis != NULL); // make sure not volume
638      basis->resize(si_->getNumRows()+numberCuts,numcols + 1) ;
639      for (int i = 0 ; i < numberCuts ; i++) {
640        basis->setArtifStatus(si_->getNumRows()+i,
641            CoinWarmStartBasis::basic) ;
642      }
643
644      const OsiRowCut ** addCuts = new const OsiRowCut * [numberCuts] ;
645      for (int i = 0 ; i < numberCuts ; i++) {
646        addCuts[i] = &cs.rowCut(i + numberCutsBefore) ;
647      }
648      si_->applyRowCuts(numberCuts,addCuts) ;
649
650      delete [] addCuts ;
651      if (si_->setWarmStart(basis) == false) {
652        delete basis;
653        throw CoinError("Fail setWarmStart() after cut installation.",
654            "generateCuts","IpCbcOACutGenerator2") ;
655      }
656      delete basis;
657      si_->resolve();
658      double objvalue = si_->getObjValue();
659      //milpBound = max(milpBound, si_->getObjValue());
660      feasible = (si_->isProvenOptimal() &&
661          !si_->isDualObjectiveLimitReached() && (objvalue<cutoff)) ;
662      //if value of integers are unchanged then we have to get out
663      bool changed = !nlp_->isProvenOptimal()//if nlp_ is infeasible solution has necessarily changed
664          && feasible;//if lp is infeasible we don't have to check anything
665      const double * nlpSol = nlp_->getColSolution();
666      const double * lpSol = si_->getColSolution();
667      for(int i = 0; i < numcols && !changed; i++) {
668        if(nlp_->isInteger(i) && fabs(nlpSol[i] - lpSol[i])>cbcIntegerTolerance_)
669          changed = 1;
670      }
671      if(changed) {
672        isInteger = 1;
673        for(int i = 0 ; i < numcols ; i++) {
674          if(nlp_->isInteger(i)) {
675            if(fabs(lpSol[i] - floor(lpSol[i] + 0.5) ) > cbcIntegerTolerance_ ) {
676              isInteger = 0;
677              break;
678            }
679          }
680        }
681      }
682      else {
683        isInteger = 0;
684        //        if(!fixed)//fathom on bounds
685        milpBound = 1e200;
686      }
687#ifdef OA_DEBUG
688
689      printf("Obj value after cuts %g %d rows\n",si_->getObjValue(),
690          numberCuts) ;
691#endif
692      //do we perform a new local search ?
693      if(milpOptimal && feasible && !isInteger &&
694          nLocalSearch_ < maxLocalSearch_ &&
695          numberPasses < maxLocalSearchPerNode_ &&
696          localSearchNodeLimit_ > 0 &&
697          CoinCpuTime() - timeBegin_ < maxLocalSearchTime_) {
698        if(clp==NULL) {
699          nLocalSearch_++;
700
701          setTimeLimit(maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
702          setCutoff(cutoff);
703          si_->branchAndBound();
704          //Did we find any integer feasible solution
705          milpBound = siBestObj();
706          //If integer solution is the same as nlp solution problem is solved
707          changed = !nlp_->isProvenOptimal();
708          for(int i = 0; i < numcols && !changed; i++) {
709            if(nlp_->isInteger(i) && fabs(nlp_->getColSolution()[i] - si_->getColSolution()[i])>cbcIntegerTolerance_)
710              changed = 1;
711          }
712          if(changed) {
713            feasible = (milpBound < cutoff);
714          }
715          else
716           {
717            feasible = 0;
718            milpBound = 1e50;
719           }
720          isInteger = si_->getFractionalIndices().size() == 0;
721        }/** endif solved by branchAndBound capable OsiInterface*/
722        else {
723          if(model)
724            {
725              OAModel = NULL;       
726              delete model;
727            }
728          model = new CbcModel(*clp); // which clones
729          OAModel = model;
730          OsiBabSolver empty;
731          model->solver()->setAuxiliaryInfo(&empty);
732          //Change Cbc messages prefixes
733          strcpy(model->messagesPointer()->source_,"OaCbc");
734          model->solver()->messageHandler()->setLogLevel(0);
735
736          if(!strategy)
737            strategy = new CbcStrategyDefault(1,0,0, subMilpLogLevel_);
738          model->setLogLevel(subMilpLogLevel_);
739          model->setStrategy(*strategy);
740          model->setMaximumNodes(localSearchNodeLimit_);
741          model->setMaximumSeconds(maxLocalSearchTime_ + timeBegin_ - CoinCpuTime());
742          model->setCutoff(cutoff);
743          model->branchAndBound();
744          nLocalSearch_++;
745          milpBound = siBestObj(model);
746          handler_->message(SOLVED_LOCAL_SEARCH, messages_)<<model->getNodeCount()<<model->getIterationCount()<<CoinMessageEol;
747
748          feasible =  (milpBound < cutoff);
749          milpFeasible = feasible;
750          isInteger = model->getSolutionCount();
751          if(model->isProvenOptimal() || model->isProvenInfeasible()) {
752            bool feasible = model->isProvenOptimal();
753
754            bool changed = 0;             //If integer solution is the same as nlp solution problem is solved
755            for(int i = 0; feasible && isInteger && i < numcols && !changed; i++) {
756              if(nlp_->isInteger(i) && fabs(nlp_->getColSolution()[i] - model->bestSolution()[i])>cbcIntegerTolerance_)
757                changed = 1;
758            }
759          if(!changed) {
760            feasible = 0;
761            milpBound = 1e50;
762           }
763            milpFeasible = feasible;
764          }
765          else {
766            std::cout<<"Exiting on time limit"<<std::endl;
767            milpOptimal = 0;
768            //feasible = 1;
769            break;
770          }
771
772        }/** endif solved by clp/cbc*/
773        if(milpBound < cutoff)
774          handler_->message(UPDATE_LB, messages_)
775          <<milpBound<<CoinCpuTime() - timeBegin_
776          <<CoinMessageEol;
777        else
778          handler_->message(OASUCCESS, messages_)
779          <<CoinCpuTime() - timeBegin_ <<CoinMessageEol;
780      }/** endif localSearch*/
781      else if(model!=NULL)/** have to delete model for next iteration*/
782      { 
783        delete model;
784        OAModel = NULL;
785        model=NULL;
786      }
787
788    }
789
790
791
792  }
793#ifdef OA_DEBUG
794    std::cout<<"-------------------------------------------------------------------------------------------------------------------------"<<std::endl;
795    std::cout<<"OA cut generation finished"<<std::endl;
796    std::cout<<"Generated "<<cs.sizeRowCuts()<<std::endl;
797    if(foundSolution)
798      std::cout <<"Found NLP-integer feasible solution of  value : "<<cutoff<<std::endl;
799    std::cout<<"Current MILP lower bound is : "<<milpBound<<std::endl;
800    std::cout<<"-------------------------------------------------------------------------------------------------------------------------"<<std::endl;
801    std::cout<<"Stopped because : isInteger "<<isInteger<<", feasible "<<feasible<<std::endl<<std::endl;
802#endif
803    //Transmit the bound found by the milp
804    {
805      if(milpBound>-1e100)
806      {
807#if 1
808        // Also store into solver
809        if (babInfo)
810          babInfo->setMipBound(milpBound);
811#endif
812
813      }
814    }  //Clean everything :
815
816  //free model
817  if(model!= NULL) {
818    delete model;
819  }
820
821
822
823  //  Reset bounds in the NLP
824
825  for(int i = 0; i < numcols ; i++) {
826    if(nlp_->isInteger(i)) {
827      nlp_->setColBounds(i,saveColLb[i],saveColUb[i]);
828    }
829  }
830#ifndef NO_NULL_SI
831  if(deleteSi) {
832    delete si_;
833    si_ = NULL;
834  }
835  else
836#endif
837   {
838    if(leaveSiUnchanged_) {
839      int nRowsToDelete = si_->getNumRows() - originalRowNumber;
840      int * rowsToDelete = new int[nRowsToDelete];
841      for(int i = 0 ; i < nRowsToDelete ; i++) {
842        rowsToDelete[i] = i + originalRowNumber;
843      }
844      si_->deleteRows(nRowsToDelete, rowsToDelete);
845      delete [] rowsToDelete;
846      if(si_==&si) {
847        si_->setDblParam(OsiDualObjectiveLimit, saveCutoff);
848        if(si_->setWarmStart(saveWarmStart)==false) {
849          throw CoinError("Fail restoring the warm start at the end of procedure",
850              "generateCuts","IpCbcOACutGenerator2") ;
851        }
852        delete saveWarmStart;
853      }
854    }
855}
856
857delete [] saveColLb;
858delete [] saveColUb;
859if(strategy && ! strategy_)
860  delete strategy;
861}
Note: See TracBrowser for help on using the repository browser.