/home/coin/SVN-release/OS-2.0.1/Bonmin/src/Algorithms/OaGenerators/BonOaDecBase.cpp

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines (IBM) 2006
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // P. Bonami, International Business Machines
00007 //
00008 // Date :  12/07/2006
00009 
00010 #include <sstream>
00011 #include <climits>
00012 #include <algorithm>
00013 
00014 #include "BonOaDecBase.hpp"
00015 
00016 
00017 #include "BonminConfig.h"
00018 
00019 #include "OsiClpSolverInterface.hpp"
00020 
00021 #include "CbcModel.hpp"
00022 #include "CbcStrategy.hpp"
00023 #ifdef COIN_HAS_CPX
00024 #include "OsiCpxSolverInterface.hpp"
00025 #endif
00026 #include "OsiAuxInfo.hpp"
00027 
00028 //The following two are to interupt the solution of sub-mip through CTRL-C
00029 extern CbcModel * OAModel;
00030 
00031 namespace Bonmin
00032 {
00033 
00034   OaDecompositionBase::OaDecompositionBase
00035   (OsiTMINLPInterface * nlp,
00036    OsiSolverInterface * si,
00037    CbcStrategy * strategy,
00038    double cbcCutoffIncrement,
00039    double cbcIntegerTolerance,
00040    bool leaveSiUnchanged
00041   )
00042       :
00043       CglCutGenerator(),
00044       nlp_(nlp),
00045       nSolve_(0),
00046       lp_(si),
00047       objects_(NULL),
00048       nObjects_(0),
00049       nLocalSearch_(0),
00050       handler_(NULL),
00051       leaveSiUnchanged_(leaveSiUnchanged),
00052       reassignLpsolver_(false),
00053       timeBegin_(0),
00054       parameters_()
00055   {
00056     handler_ = new CoinMessageHandler();
00057     handler_ -> setLogLevel(2);
00058     messages_ = OaMessages();
00059     if (strategy)
00060       parameters_.setStrategy(*strategy);
00061     timeBegin_ = CoinCpuTime();
00062     parameters_.cbcCutoffIncrement_  = cbcCutoffIncrement;
00063     parameters_.cbcIntegerTolerance_ = cbcIntegerTolerance;
00064   }
00065   OaDecompositionBase::OaDecompositionBase(BabSetupBase &b, bool leaveSiUnchanged,
00066       bool reassignLpsolver):
00067       CglCutGenerator(),
00068       nlp_(b.nonlinearSolver()),
00069       lp_(b.continuousSolver()),
00070       objects_(NULL),
00071       nObjects_(0),
00072       nLocalSearch_(0),
00073       handler_(NULL),
00074       leaveSiUnchanged_(leaveSiUnchanged),
00075       reassignLpsolver_(reassignLpsolver),
00076       timeBegin_(0),
00077       parameters_()
00078   {
00079     handler_ = new CoinMessageHandler();
00080     int logLevel;
00081     b.options()->GetIntegerValue("oa_log_level",logLevel,"bonmin.");
00082     b.options()->GetNumericValue("oa_log_frequency",parameters_.logFrequency_,"bonmin.");
00083 
00084     handler_ -> setLogLevel(logLevel);
00085 
00086     messages_ = OaMessages();
00087     timeBegin_ = CoinCpuTime();
00088     b.options()->GetIntegerValue("milp_log_level",parameters_.subMilpLogLevel_,"bonmin.");
00089     b.options()->GetNumericValue("cutoff_decr",parameters_.cbcCutoffIncrement_,"bonmin.");
00090     b.options()->GetNumericValue("integer_tolerance",parameters_.cbcIntegerTolerance_,"bonmin.");
00091     int ivalue;
00092     b.options()->GetEnumValue("add_only_violated_oa", ivalue,"bonmin.");
00093     parameters_.addOnlyViolated_ = ivalue;
00094     b.options()->GetEnumValue("oa_cuts_scope", ivalue,"bonmin.");
00095     parameters_.global_ = ivalue;
00096   }
00097 
00098   OaDecompositionBase::OaDecompositionBase
00099   (const OaDecompositionBase & other)
00100       :
00101       CglCutGenerator(other),
00102       nlp_(other.nlp_),
00103       lp_(other.lp_),
00104       objects_(other.objects_),
00105       nObjects_(other.nObjects_),
00106       nLocalSearch_(0),
00107       messages_(other.messages_),
00108       leaveSiUnchanged_(other.leaveSiUnchanged_),
00109       reassignLpsolver_(other.reassignLpsolver_),
00110       timeBegin_(other.timeBegin_),
00111       parameters_(other.parameters_)
00112   {
00113     //timeBegin_ = CoinCpuTime();
00114     handler_ = other.handler_->clone();
00115   }
00117   OaDecompositionBase::Parameters::Parameters():
00118       global_(true),
00119       addOnlyViolated_(false),
00120       cbcCutoffIncrement_(1e-06),
00121       cbcIntegerTolerance_(1e-05),
00122       localSearchNodeLimit_(0),
00123       maxLocalSearchPerNode_(0),
00124       maxLocalSearch_(0),
00125       maxLocalSearchTime_(3600),
00126       subMilpLogLevel_(0),
00127       logFrequency_(1000.),
00128       strategy_(NULL)
00129   {}
00130 
00132   OaDecompositionBase::~OaDecompositionBase()
00133   {
00134     delete handler_;
00135   }
00136 
00137 
00139   OaDecompositionBase::Parameters::Parameters(const Parameters & other):
00140       global_(other.global_),
00141       addOnlyViolated_(other.addOnlyViolated_),
00142       cbcCutoffIncrement_(other.cbcCutoffIncrement_),
00143       cbcIntegerTolerance_(other.cbcIntegerTolerance_),
00144       localSearchNodeLimit_(other.localSearchNodeLimit_),
00145       maxLocalSearchPerNode_(other.maxLocalSearchPerNode_),
00146       maxLocalSearch_(other.maxLocalSearch_),
00147       maxLocalSearchTime_(other.maxLocalSearchTime_),
00148       subMilpLogLevel_(other.subMilpLogLevel_),
00149       logFrequency_(other.logFrequency_),
00150       strategy_(NULL)
00151   {
00152     if (other.strategy_)
00153       strategy_ = other.strategy_->clone();
00154   }
00155 
00156 
00158   OaDecompositionBase::SubMipSolver::SubMipSolver(OsiSolverInterface * lp,
00159       const CbcStrategy * strategy):
00160       lp_(lp),
00161       clp_(NULL),
00162       cpx_(NULL),
00163       cbc_(NULL),
00164       lowBound_(-COIN_DBL_MAX),
00165       optimal_(false),
00166       integerSolution_(NULL),
00167       strategy_(NULL)
00168   {
00169     clp_ = (lp_ == NULL)? NULL :
00170         dynamic_cast<OsiClpSolverInterface *>(lp_);
00171 #ifdef COIN_HAS_CPX
00172     cpx_ = (lp_ == NULL)? NULL :
00173         dynamic_cast<OsiCpxSolverInterface *>(lp_);
00174 #endif
00175     if (strategy) strategy_ = strategy->clone();
00176   }
00177   OaDecompositionBase::SubMipSolver::~SubMipSolver()
00178   {
00179     if (strategy_) delete strategy_;
00180     if (integerSolution_) delete [] integerSolution_;
00181     if (cbc_) delete cbc_;
00182   }
00183 
00185   void
00186   OaDecompositionBase::SubMipSolver::setLpSolver(OsiSolverInterface * lp)
00187   {
00188     lp_ = lp;
00189     clp_ = (lp_ == NULL) ? NULL :
00190         dynamic_cast<OsiClpSolverInterface *>(lp_);
00191 #ifdef COIN_HAS_CPX
00192     cpx_ = (lp_ == NULL) ? NULL :
00193         dynamic_cast<OsiCpxSolverInterface *>(lp_);
00194 #endif
00195     lowBound_ = -COIN_DBL_MAX;
00196     optimal_ = false;
00197     if (integerSolution_) {
00198       delete [] integerSolution_;
00199       integerSolution_ = NULL;
00200     }
00201   }
00202 
00203 
00204 
00205   void
00206   OaDecompositionBase::SubMipSolver::performLocalSearch(double cutoff, int loglevel, double maxTime,
00207       int maxNodes)
00208   {
00209     if (clp_) {
00210       if (!strategy_)
00211         strategy_ = new CbcStrategyDefault(1,0,0, loglevel);
00212 
00213       OsiBabSolver empty;
00214       if (cbc_) delete cbc_;
00215       OAModel = cbc_ = new CbcModel(*clp_);
00216       cbc_->solver()->setAuxiliaryInfo(&empty);
00217 
00218       //Change Cbc messages prefixes
00219       strcpy(cbc_->messagesPointer()->source_,"OaCbc");
00220 
00221       cbc_->setLogLevel(loglevel);
00222       cbc_->solver()->messageHandler()->setLogLevel(0);
00223       clp_->resolve();
00224       cbc_->setStrategy(*strategy_);
00225       cbc_->setLogLevel(loglevel);
00226       cbc_->solver()->messageHandler()->setLogLevel(0);
00227       cbc_->setMaximumNodes(maxNodes);
00228       cbc_->setMaximumSeconds(maxTime);
00229       cbc_->setCutoff(cutoff);
00230       printf("Solving local search");
00231       cbc_->branchAndBound();
00232       OAModel = NULL;
00233       lowBound_ = cbc_->getBestPossibleObjValue();
00234 
00235       if (cbc_->isProvenOptimal() || cbc_->isProvenInfeasible())
00236         optimal_ = true;
00237       else optimal_ = false;
00238 
00239       if (cbc_->getSolutionCount()) {
00240         if (!integerSolution_)
00241           integerSolution_ = new double[lp_->getNumCols()];
00242         CoinCopyN(cbc_->bestSolution(), lp_->getNumCols(), integerSolution_);
00243       }
00244       else if (integerSolution_) {
00245         delete [] integerSolution_;
00246         integerSolution_ = NULL;
00247       }
00248       nodeCount_ = cbc_->getNodeCount();
00249       iterationCount_ = cbc_->getIterationCount();
00250     }
00251     else {
00252       lp_->messageHandler()->setLogLevel(loglevel);
00253 #ifdef COIN_HAS_CPX
00254       if (cpx_) {
00255         CPXENVptr env = cpx_->getEnvironmentPtr();
00256         CPXsetintparam(env, CPX_PARAM_NODELIM, maxNodes);
00257         CPXsetdblparam(env, CPX_PARAM_TILIM, maxTime);
00258         CPXsetdblparam(env, CPX_PARAM_CUTUP, cutoff);
00259         //CpxModel = cpx_;
00260       }
00261       else
00262 #endif 
00263      {
00264         throw CoinError("Unsuported solver, for local searches you should use clp or cplex",
00265             "performLocalSearch",
00266             "OaDecompositionBase::SubMipSolver");
00267     }
00268 
00269     lp_->branchAndBound();
00270 
00271    optimal_ = lp_->isProvenOptimal();
00272 #ifdef COIN_HAS_CPX
00273     if (cpx_) {
00274       //CpxModel = NULL;
00275       CPXENVptr env = cpx_->getEnvironmentPtr();
00276       CPXLPptr cpxlp = cpx_->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
00277 
00278        
00279       int status = CPXgetbestobjval(env, cpxlp, &lowBound_);
00280      
00281       int stat = CPXgetstat( env, cpxlp);
00282       optimal_ |= (stat == CPXMIP_INFEASIBLE); 
00283        nodeCount_ = CPXgetnodecnt(env , cpxlp);
00284       iterationCount_ = CPXgetmipitcnt(env , cpxlp);
00285       if (status)
00286         throw CoinError("Error in getting some CPLEX information","OaDecompositionBase::SubMipSolver","performLocalSearch");
00287     }
00288 #endif
00289 
00290     if (lp_->getFractionalIndices().size() == 0) {
00291       if (!integerSolution_)
00292         integerSolution_ = new double[lp_->getNumCols()];
00293       CoinCopyN(lp_->getColSolution(), lp_->getNumCols() , integerSolution_);
00294     }
00295     else if (integerSolution_) {
00296       delete [] integerSolution_;
00297       integerSolution_ = NULL;
00298     }
00299   }
00300 }
00301 
00302 OaDecompositionBase::solverManip::solverManip
00303 (OsiSolverInterface * si,
00304  bool saveNumRows,
00305  bool saveBasis,
00306  bool saveBounds,
00307  bool saveCutoff,
00308  bool resolve):
00309     si_(si),
00310     initialNumberRows_(-1),
00311     colLower_(NULL),
00312     colUpper_(NULL),
00313     warm_(NULL),
00314     cutoff_(COIN_DBL_MAX),
00315     deleteSolver_(false),
00316     objects_(NULL),
00317     nObjects_(0),
00318     integerTolerance_(1e-08)
00319 {
00320   getCached();
00321   if (saveNumRows)
00322     initialNumberRows_ = numrows_;
00323   if (saveBasis)
00324     warm_ = si->getWarmStart();
00325   if (saveBounds) {
00326     colLower_ = new double[numcols_];
00327     colUpper_ = new double[numcols_];
00328     CoinCopyN(si->getColLower(), numcols_ , colLower_);
00329     CoinCopyN(si->getColUpper(), numcols_ , colUpper_);
00330   }
00331   if (saveCutoff)
00332     si->getDblParam(OsiDualObjectiveLimit, cutoff_);
00333   si->messageHandler()->setLogLevel(0);
00334   if (resolve) si->resolve();
00335 }
00336 
00337 
00338 OaDecompositionBase::solverManip::solverManip
00339 (const OsiSolverInterface & si):
00340     si_(NULL),
00341     initialNumberRows_(-1),
00342     colLower_(NULL),
00343     colUpper_(NULL),
00344     warm_(NULL),
00345     cutoff_(COIN_DBL_MAX),
00346     deleteSolver_(true),
00347     objects_(NULL),
00348     nObjects_(0),
00349     integerTolerance_(1e-08)
00350 {
00351   si_ = si.clone();
00352   getCached();
00353 }
00354 
00355 OaDecompositionBase::solverManip::~solverManip()
00356 {
00357   if (warm_) delete warm_;
00358   if (colLower_) delete [] colLower_;
00359   if (colUpper_) delete [] colUpper_;
00360   if (deleteSolver_) delete si_;
00361 }
00362 
00363 void
00364 OaDecompositionBase::solverManip::restore()
00365 {
00366   if (initialNumberRows_ >= 0) {
00367     int nRowsToDelete = numrows_ - initialNumberRows_;
00368     int * rowsToDelete = new int[nRowsToDelete];
00369     for (int i = 0 ; i < nRowsToDelete ; i++) {
00370       rowsToDelete[i] = i + initialNumberRows_;
00371     }
00372     si_->deleteRows(nRowsToDelete, rowsToDelete);
00373     delete [] rowsToDelete;
00374     numrows_ -= nRowsToDelete;
00375   }
00376 
00377   if (colLower_) {
00378     si_->setColLower(colLower_);
00379   }
00380 
00381   if (colUpper_) {
00382     si_->setColUpper(colUpper_);
00383   }
00384 
00385   if (cutoff_<COIN_DBL_MAX) {
00386     si_->setDblParam(OsiDualObjectiveLimit, cutoff_);
00387   }
00388 
00389   if (warm_) {
00390     if (si_->setWarmStart(warm_)==false) {
00391       throw CoinError("Fail restoring the warm start at the end of procedure",
00392           "restore","OaDecompositionBase::SaveSolverState") ;
00393     }
00394   }
00395   getCached();
00396 }
00397 
00398 void
00399 OaDecompositionBase::passInMessageHandler(CoinMessageHandler * handler)
00400 {
00401   int logLevel = handler_->logLevel();
00402   delete handler_;
00403   handler_=handler->clone();
00404   handler_->setLogLevel(logLevel);
00405 }
00406 
00408 bool
00409 OaDecompositionBase::solverManip::integerFeasible(const OsiBranchingInformation& info) const
00410 {
00411   if (objects_) {
00412     int dummy;
00413     for (int i = 0 ; i < nObjects_ ; i++) {
00414       if (objects_[i]->infeasibility(&info, dummy) > 0.0) return false;
00415     }
00416   }
00417   else {
00418     const double * sol = info.solution_;
00419     int numcols = si_->getNumCols();
00420     for (int i = 0 ; i < numcols ; i++) {
00421       if (si_->isInteger(i)) {
00422         if (fabs(sol[i] - floor(sol[i] + 0.5)) >
00423             integerTolerance_) {
00424           return false;
00425         }
00426       }
00427     }
00428   }
00429   return true;
00430 }
00431 
00434 void
00435 OaDecompositionBase::solverManip::fixIntegers(const OsiBranchingInformation& info)
00436 {
00437   if (objects_) {
00438     for (int i = 0 ; i < nObjects_ ; i++) {
00439       if (objects_[i]->feasibleRegion(si_, &info));
00440     }
00441   }
00442   else {
00443     const double * colsol = info.solution_;
00444     for (int i = 0; i < numcols_; i++) {
00445       if (si_->isInteger(i)) {
00446         double  value =  colsol[i];
00447         if (fabs(value - floor(value+0.5)) > integerTolerance_) {
00448           std::stringstream stream;
00449           stream<<"Error not integer valued solution"<<std::endl;
00450           stream<<"---------------- x["<<i<<"] = "<<value<<std::endl;
00451           throw CoinError(stream.str(),"fixIntegers","OaDecompositionBase::solverManip");
00452         }
00453         value = floor(value+0.5);
00454         value = std::max(colLower_[i],value);
00455         value = std::min(value, colUpper_[i]);
00456 
00457         if (fabs(value) > 1e10) {
00458           std::stringstream stream;
00459           stream<<"Can not fix variable in nlp because it has too big a value ("<<value
00460           <<") at optimium of LP relaxation. You should try running the problem with B-BB"<<std::endl;
00461           throw CoinError(stream.str(),
00462               "fixIntegers","OaDecompositionBase::solverManip") ;
00463         }
00464 #ifdef OA_DEBUG
00465         //         printf("xx %d at %g (bounds %g, %g)",i,value,nlp_->getColLower()[i],
00466         //                nlp_->getColUpper()[i]);
00467         std::cout<<(int)value;
00468 #endif
00469         si_->setColLower(i,value);
00470         si_->setColUpper(i,value);
00471       }
00472     }
00473 #ifdef OA_DEBUG
00474     std::cout<<std::endl;
00475 #endif
00476   }
00477 }
00479 bool
00480 OaDecompositionBase::solverManip::isDifferentOnIntegers(const double * colsol)
00481 {
00482   return isDifferentOnIntegers(colsol, si_->getColSolution());
00483 }
00484 
00486 bool
00487 OaDecompositionBase::solverManip::isDifferentOnIntegers(const double * colsol, const double *otherSol)
00488 {
00489   if (objects_) {
00490     for (int i = 0 ; i < nObjects_ ; i++) {
00491       OsiObject * obj = objects_[i];
00492       int colnum = obj->columnNumber();
00493       if (colnum >= 0) {//Variable branching object
00494         if (fabs(otherSol[colnum] - colsol[colnum]) > 100*integerTolerance_) {
00495           return true;
00496         }
00497       }
00498       else {//It is a sos branching object
00499         OsiSOS * sos = dynamic_cast<OsiSOS *>(obj);
00500         assert(sos);
00501         const int * members = sos->members();
00502         int end = sos->numberMembers();
00503         for (int k = 0 ; k < end ; k++) {
00504           if (fabs(otherSol[members[k]] - colsol[members[k]]) > 0.001) {
00505             return true;
00506           }
00507         }
00508       }
00509     }
00510   }
00511   else {
00512     for (int i = 0; i < numcols_ ; i++) {
00513       if (si_->isInteger(i) && fabs(otherSol[i] - colsol[i])>0.001)
00514         return true;
00515     }
00516   }
00517   return false;
00518 }
00519 
00521 void
00522 OaDecompositionBase::solverManip::cloneOther(const OsiSolverInterface &si)
00523 {
00524   //Install current active cuts into local solver
00525   int numberCutsToAdd = si.getNumRows();
00526   numberCutsToAdd -= numrows_;
00527   if (numberCutsToAdd > 0)//Have to install some cuts
00528   {
00529     CoinPackedVector * * addCuts = new CoinPackedVector *[numberCutsToAdd];
00530     for (int i = 0 ; i < numberCutsToAdd ; i++)
00531     {
00532       addCuts[i] = new CoinPackedVector;
00533     }
00534     //Get the current matrix and fill the addCuts
00535     const CoinPackedMatrix * mat = si.getMatrixByCol();
00536     const CoinBigIndex * start = mat->getVectorStarts();
00537     const int * length = mat->getVectorLengths();
00538     const double * elements = mat->getElements();
00539     const int * indices = mat->getIndices();
00540     for (int i = 0 ; i < numcols_ ; i++)
00541       for (int k = start[i] ; k < start[i] + length[i] ; k++)
00542       {
00543         if (indices[k] >= numrows_) {
00544           addCuts[ indices[k] - numrows_ ]->insert(i, elements[k]);
00545         }
00546       }
00547     si_->addRows(numberCutsToAdd, (const CoinPackedVectorBase * const *) addCuts, si.getRowLower() + numrows_,
00548         si.getRowUpper() + numrows_);
00549     for (int i = 0 ; i < numberCutsToAdd ; i++){
00550       delete addCuts[i];
00551     }
00552     delete [] addCuts;
00553   }
00554   else if (numberCutsToAdd < 0) { //Oups some error
00555     throw CoinError("Internal error number of cuts wrong",
00556         "generateCuts","OACutGenerator2");
00557   }
00558 
00559   si_->setColLower(si.getColLower());
00560   si_->setColUpper(si.getColUpper());
00561   //Install basis in problem
00562   CoinWarmStart * warm = si.getWarmStart();
00563   if (si_->setWarmStart(warm)==false) {
00564     delete warm;
00565     throw CoinError("Fail installing the warm start in the subproblem",
00566         "generateCuts","OACutGenerator2") ;
00567   }
00568   delete warm;
00569   //put the cutoff
00570   double cutoff;
00571   si.getDblParam(OsiDualObjectiveLimit, cutoff);
00572   si_->setDblParam(OsiDualObjectiveLimit, cutoff);
00573   si_->messageHandler()->setLogLevel(0);
00574   si_->resolve();
00575 
00576   numrows_ = si.getNumRows();
00577 #ifdef OA_DEBUG
00578 
00579   std::cout<<"Resolve with hotstart :"<<std::endl
00580   <<"number of iterations(should be 0) : "<<lp_->getIterationCount()<<std::endl
00581   <<"Objective value and diff to original : "<<lp_->getObjValue()<<", "
00582   <<fabs(si_->getObjValue() - si.getObjValue())<<std::endl;
00583   for (int i = 0 ; i <= numcols ; i++) {
00584     if (fabs(si.getColSolution()[i]-si_->getColSolution()[i])>1e-08) {
00585       std::cout<<"Diff between solution at node and solution with local solver : "<<fabs(si.getColSolution()[i]-si_->getColSolution()[i])<<std::endl;
00586     }
00587   }
00588 #endif
00589 
00590 }
00591 
00592 
00594 void
00595 OaDecompositionBase::solverManip::installCuts(const OsiCuts& cs, int numberCuts)
00596 {
00597   int numberCutsBefore = cs.sizeRowCuts() - numberCuts;
00598 
00599   CoinWarmStartBasis * basis
00600   = dynamic_cast<CoinWarmStartBasis*>(si_->getWarmStart()) ;
00601   assert(basis != NULL); // make sure not volume
00602   basis->resize(numrows_ + numberCuts,numcols_) ;
00603   for (int i = 0 ; i < numberCuts ; i++) {
00604     basis->setArtifStatus(numrows_ + i,
00605         CoinWarmStartBasis::basic) ;
00606   }
00607 
00608   const OsiRowCut ** addCuts = new const OsiRowCut * [numberCuts] ;
00609   for (int i = 0 ; i < numberCuts ; i++) {
00610     addCuts[i] = &cs.rowCut(i + numberCutsBefore) ;
00611   }
00612   si_->applyRowCuts(numberCuts,addCuts) ;
00613   numrows_ += numberCuts;
00614   delete [] addCuts ;
00615   if (si_->setWarmStart(basis) == false) {
00616     delete basis;
00617     throw CoinError("Fail setWarmStart() after cut installation.",
00618         "generateCuts","OACutGenerator2") ;
00619   }
00620   delete basis;
00621 }
00622 
00624 void
00625 OaDecompositionBase::generateCuts(const OsiSolverInterface &si,  OsiCuts & cs,
00626     const CglTreeInfo info) const
00627 {
00628   if (nlp_ == NULL) {
00629     throw CoinError("Error in cut generator for outer approximation no NLP ipopt assigned", "generateCuts", "OaDecompositionBase");
00630   }
00631 
00632   // babInfo is used to communicate with the b-and-b solver (Cbc or Bcp).
00633   OsiBabSolver * babInfo = dynamic_cast<OsiBabSolver *> (si.getAuxiliaryInfo());
00634   if (babInfo)
00635     if (!babInfo->mipFeasible())
00636       return;
00637 
00638   //Get the continuous solution
00639   const double *colsol = si.getColSolution();
00640 
00641 
00642   solverManip nlpManip(nlp_, false, false, true, false, false);
00643   nlpManip.setIntegerTolerance(parameters_.cbcIntegerTolerance_);
00644   nlpManip.setObjects(objects_, nObjects_);
00645   OsiBranchingInformation brInfo(nlp_, false);
00646   brInfo.solution_ = colsol;
00647   //Check integer infeasibility
00648   bool isInteger = nlpManip.integerFeasible(brInfo);
00649 
00650   SubMipSolver * subMip = NULL;
00651 
00652   if (!isInteger) {
00653     if (doLocalSearch())//create sub mip solver.
00654     {
00655       subMip = new SubMipSolver(lp_, parameters_.strategy());
00656     }
00657     else {
00658       return;
00659     }
00660   }
00661 
00662 
00663   //If we are going to modify things copy current information to restore it in the end
00664 
00665 
00666   //get the current cutoff
00667   double cutoff;
00668   si.getDblParam(OsiDualObjectiveLimit, cutoff);
00669 
00670   // Save solvers state if needed
00671 
00672   solverManip * lpManip = NULL;
00673   if (lp_ != NULL) {
00674     if (lp_!=&si) {
00675       lpManip = new solverManip(lp_, true, false, false, true, true);
00676       lpManip->cloneOther(si);
00677     }
00678     else {
00679 #if 0
00680       throw CoinError("Not allowed to modify si in a cutGenerator",
00681           "OACutGenerator2","generateCuts");
00682 #else
00683       lpManip = new solverManip(lp_, true, leaveSiUnchanged_, true, true);
00684 #endif
00685     }
00686   }
00687   else {
00688     lpManip = new solverManip(si);
00689   }
00690   lpManip->setObjects(objects_, nObjects_);
00691   lpManip->setIntegerTolerance(parameters_.cbcIntegerTolerance_);
00692   double milpBound = performOa(cs, nlpManip, *lpManip, subMip, babInfo, cutoff);
00693 
00694   //Transmit the bound found by the milp
00695   {
00696     if (milpBound>-1e100)
00697     {
00698       // Also store into solver
00699       if (babInfo)
00700         babInfo->setMipBound(milpBound);
00701     }
00702   }  //Clean everything :
00703 
00704   //free subMip
00705   if (subMip!= NULL) {
00706     delete subMip;
00707     subMip = NULL;
00708   }
00709 
00710   //  Reset the two solvers
00711   if (leaveSiUnchanged_)
00712     lpManip->restore();
00713   delete lpManip;
00714   nlpManip.restore();
00715   return;
00716 }
00717 
00718 void
00719 OaDecompositionBase::solverManip::getCached()
00720 {
00721   numrows_ = si_->getNumRows();
00722   numcols_ = si_->getNumCols();
00723   siColLower_ = si_->getColLower();
00724   siColUpper_ = si_->getColUpper();
00725 }
00726 
00727 
00729 bool
00730 OaDecompositionBase::solveNlp(OsiBabSolver * babInfo, double cutoff) const
00731 {
00732   nSolve_++;
00733   nlp_->resolve();
00734   bool return_value = false;
00735   if (nlp_->isProvenOptimal()) {
00736     handler_->message(FEASIBLE_NLP, messages_)
00737     <<nlp_->getIterationCount()
00738     <<nlp_->getObjValue()<<CoinMessageEol;
00739 
00740 #ifdef OA_DEBUG
00741     const double * colsol2 = nlp_->getColSolution();
00742     debug_.checkInteger(*nlp_,std::cerr);
00743 #endif
00744 
00745     if ((nlp_->getObjValue() < cutoff) ) {
00746       handler_->message(UPDATE_UB, messages_)
00747       <<nlp_->getObjValue()
00748       <<CoinCpuTime()-timeBegin_
00749       <<CoinMessageEol;
00750 
00751       return_value = true;
00752       // Also pass it to solver
00753       if (babInfo) {
00754         int numcols = nlp_->getNumCols();
00755         double * lpSolution = new double[numcols + 1];
00756         CoinCopyN(nlp_->getColSolution(), numcols, lpSolution);
00757         lpSolution[numcols] = nlp_->getObjValue();
00758         babInfo->setSolution(lpSolution,
00759             numcols + 1, lpSolution[numcols]);
00760         delete [] lpSolution;
00761       }
00762       else {
00763         printf("No auxiliary info in nlp solve!\n");
00764         throw -1;
00765       }
00766     }
00767   }
00768   else if (nlp_->isAbandoned() || nlp_->isIterationLimitReached()) {
00769     (*handler_)<<"Unsolved NLP... exit"<<CoinMessageEol;
00770   }
00771   else {
00772     handler_->message(INFEASIBLE_NLP, messages_)
00773     <<nlp_->getIterationCount()
00774     <<CoinMessageEol;
00775   }
00776   return return_value;
00777 }
00778 
00779 
00780 
00781 #ifdef OA_DEBUG
00782 bool
00783 OaDecompositionBase::OaDebug::checkInteger(const OsiSolverInterface &nlp, std::ostream & os) const
00784 {
00785    const double * colsol = nlp.getColSolution();
00786    int numcols = nlp.getNumCols();
00787   for (int i = 0 ; i < numcols ; i++) {
00788     if (nlp.isInteger(i)) {
00789       if (fabs(colsol[i]) - floor(colsol[i] + 0.5) >
00790           1e-07) {
00791         std::cerr<<"Integer infeasible point (should not be), integer infeasibility for variable "<<i
00792         <<" is, "<<fabs(colsol[i] - floor(colsol[i] + 0.5))<<std::endl;
00793       }
00794     }
00795     return true;
00796   }
00797 
00798 }
00799 
00800 void
00801 OaDecompositionBase::OaDebug::printEndOfProcedureDebugMessage(const OsiCuts &cs,
00802     bool foundSolution,
00803     double solValue,
00804     double milpBound,
00805     bool isInteger,
00806     bool feasible,
00807     std::ostream & os) const
00808 {
00809   std::cout<<"------------------------------------------------------------------"
00810   <<std::endl;
00811   std::cout<<"OA procedure finished"<<std::endl;
00812   std::cout<<"Generated "<<cs.sizeRowCuts()<<std::endl;
00813   if (foundSolution)
00814     std::cout <<"Found NLP-integer feasible solution of  value : "<<solValue<<std::endl;
00815   std::cout<<"Current MILP lower bound is : "<<milpBound<<std::endl;
00816   std::cout<<"-------------------------------------------------------------------"<<std::endl;
00817   std::cout<<"Stopped because : isInteger "<<isInteger<<", feasible "<<feasible<<std::endl<<std::endl;
00818 
00819 }
00820 
00821 
00822 
00823 #endif
00824 }/* End namespace Bonmin. */
00825 

Generated on Thu Oct 8 03:02:54 2009 by  doxygen 1.4.7