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

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

astyled the devel branch

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