/home/coin/SVN-release/OS-2.0.1/Bonmin/experimental/FP/OAFeasibilityPump.cpp

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines Corporation and Carnegie Mellon University 2006 
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // This code is published under the Common Public License.
00006 //
00007 // Authors :
00008 // Pierre Bonami, Carnegie Mellon University,
00009 //
00010 // Date : 06/18/2005
00011 // The ideas developped in this code are co-authored with: G. Cornuejols, A. Lodi, F. Margot and can be found 
00012 // in:
00013 // "A Feasibility Pump for Mixed Integer Nonlinear Programming" 
00014 //  IBM Research Report RC 23682 February 2006
00015 //
00016 
00017 #include "BonminConfig.h"
00018 
00019 #if defined(_MSC_VER)
00020 // Turn off compiler warning about long names
00021 #  pragma warning(disable:4786)
00022 #endif
00023 
00024 #include <cassert>
00025 #include <iomanip>
00026 #include <sstream>
00027 
00028 //#define EXTRACT_LIN_REL
00029 
00030 // For Branch and bound
00031 #include "OsiSolverInterface.hpp"
00032 #include "OsiClpSolverInterface.hpp"
00033 #include "CbcModel.hpp"
00034 #include "CbcBranchUser.hpp"
00035 #include "CbcCompareUser.hpp"
00036 #include "CbcCompareActual.hpp"
00037 #include "CbcCutGenerator.hpp"
00038 //#include "CbcHeuristicUser.hpp"
00039 #include "BonAmplInterface.hpp"
00040 #include "BonDummyHeuristic.hpp"
00041 #include "BonOACutGenerator2.hpp"
00042 
00043 //Use heuristics
00044 #include "CbcHeuristicFPump.hpp"
00045 #include "CbcHeuristicGreedy.hpp"
00046 
00047 #include "BonAmplTMINLP.hpp"
00048 
00049 #include "CglGomory.hpp"
00050 //#include "CglProbing.hpp"
00051 
00052 //#include "ClpQuadInterface.hpp"
00053 
00054 // Heuristics would need adapting
00055 
00056 //#include "CbcHeuristic.hpp"
00057 //#define COIN_HAS_CPX
00058 #ifdef COIN_HAS_CPX
00059 #include "OsiCpxSolverInterface.hpp"
00060 #else
00061 #include "OsiCbcSolverInterface.hpp"
00062 #endif
00063 
00064 // Time
00065 #include "CoinTime.hpp"
00066 #include "FP.hpp"
00067 
00068 using namespace Bonmin;
00069 
00070 OptParam params;
00071 static double BeginTimeGLOB;
00072 
00077 #ifdef COIN_HAS_CPX
00078 int findGoodSolution(OsiCpxSolverInterface &mip, double maxTime,
00079                      int& nTotalNodes)
00080 {
00081   double beginTime = CoinCpuTime();
00082   int saveIntSolLim;
00083   int saveNodLim;
00084   double saveTiLim;
00085   CPXCENVptr env2 = mip.getEnvironmentPtr();
00086   CPXCLPptr lp = mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL);
00087   CPXgetintparam(env2,CPX_PARAM_INTSOLLIM, &saveIntSolLim);
00088   CPXgetintparam(env2,CPX_PARAM_NODELIM, &saveNodLim);
00089   CPXgetdblparam(env2,CPX_PARAM_TILIM, &saveTiLim);
00090 
00091   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 10000);
00092   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_NODELIM, params.minNodes_);
00093   CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, maxTime);
00094 
00095   mip.branchAndBound();
00096   nTotalNodes += CPXgetnodecnt(env2,lp);
00097   int mipstat = CPXgetstat(env2,lp);
00098   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 1);
00099   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_NODELIM, params.nodeInterval_);
00100   while (mipstat == CPXMIP_NODE_LIM_INFEAS) {
00101     CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, maxTime - CoinCpuTime() + beginTime);
00102     mip.branchAndBound();
00103     nTotalNodes += CPXgetnodecnt(env2,lp);
00104     mipstat = CPXgetstat(env2,lp);
00105   }
00106   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, saveIntSolLim);
00107   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_NODELIM, saveNodLim);
00108   CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, saveTiLim);
00109   return mipstat;
00110 }
00111 #else
00112 int findGoodSolution(OsiSolverInterface &mip, int&nodeNumber)
00113 {
00114   std::cerr<<"Does not work without CPX"<<std::endl;
00115   throw -1;
00116 }
00117 
00118 #endif
00119 struct ResolutionInformation
00120 {
00121   double time;
00122   double nlp_time;
00123   double mip_time;
00124   int n_iterations;
00125   int n_nlp_iterations;
00126   int n_mip_nodes;
00127   ResolutionInformation():time(0.), nlp_time(0.), mip_time(0.), n_iterations(0), n_nlp_iterations(0), n_mip_nodes(0)
00128   {}
00129   ResolutionInformation& operator +=(ResolutionInformation &other)
00130   {
00131     time +=other.time;
00132     nlp_time+=(other.nlp_time);
00133     mip_time+=(other.mip_time);
00134     n_iterations+=(other.n_iterations);
00135     n_nlp_iterations+=other.n_nlp_iterations;
00136     n_mip_nodes+=other.n_mip_nodes;
00137     return *this;
00138   }
00139 };
00140 
00141 #if 0
00142 void
00143 writeBoundFiles(AmplInterface& nlp, const double * originalLower, const double * originalUpper)
00144 {
00145   const double * currentLower = nlp.getColLower();
00146   const double * currentUpper = nlp.getColUpper();
00147 
00148   CoinRelFltEq eq;
00149   std::string fBoundsName;
00150   nlp.getStrParam(OsiProbName,fBoundsName);
00151   fBoundsName+="_bounds_FP";
00152   std::string fModName = fBoundsName + ".mod";
00153   std::ofstream fBounds;
00154   std::ofstream fMod;
00155   bool hasVarNames = 0;
00156   if(nlp.getVarNames()!=NULL )
00157     hasVarNames=1;
00158   if(hasVarNames)
00159     fMod.open(fModName.c_str());
00160   fBounds.open(fBoundsName.c_str());
00161 
00162   for(int i = 0 ; i < nlp.getNumCols() ; i++) {
00163 
00164     if(nlp.isInteger(i) && !eq(currentLower[i],originalLower[i])) {
00165       if(hasVarNames)
00166         fMod<<"bounds"<<i<<": "
00167         <<nlp.getVarNames()[i]<<" >= "
00168         <<currentLower[i]<<";\n";
00169       fBounds<<"LO"<<"\t"<<i<<"\t"<<currentLower[i]<<std::endl;
00170     }
00171     if(nlp.isInteger(i) && !eq(currentUpper[i],originalUpper[i])) {
00172       if(hasVarNames)
00173         fMod<<"bounds"<<i<<": "
00174         <<nlp.getVarNames()[i]<<" <= "
00175         <<currentUpper[i]<<";\n";
00176       fBounds<<"UP"<<"\t"<<i<<"\t"<<currentUpper[i]<<std::endl;
00177     }
00178   }
00179 }
00180 #endif
00181 
00182 double FP(AmplInterface &nlp,
00183           OsiSolverInterface &linearModel,
00184           int numIntCols, int * inds, double * vals,
00185           double maxTime, int maxIter, ResolutionInformation& info,
00186           double ub, int &provenInfeas, double *&solution)
00187 {
00188   provenInfeas=0;
00189 #ifdef COIN_HAS_CPX
00190 
00191   OsiCpxSolverInterface * cpxSi = dynamic_cast<OsiCpxSolverInterface *>
00192                                   (&linearModel);
00193   OsiCpxSolverInterface &mip = *cpxSi;
00194 #endif
00195 
00196   int &nMajorIt = info.n_iterations;
00197   int &nTotalNodes = info.n_mip_nodes;
00198   int &nNlpIterations = info.n_nlp_iterations;
00199   double& bbTime = info.mip_time;
00200   double& nlpTime = info.nlp_time;
00201   info.time = -CoinCpuTime();
00202   double objValue = DBL_MAX;
00203   {
00204     OsiCuts cs;
00205     //First get the closest point to the current integer (NLP-infeasible) optimum
00206     nlpTime -= CoinCpuTime();
00207     double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
00208     nlpTime += CoinCpuTime();
00209     if(dist < 1e-06)//do integer infeasibility check on variables
00210     {
00211       //Something wrong shouldn't be
00212       std::cout<<"Feasibility subproblem has objective 0 while problem was claimed infeasible before"<<std::endl
00213       <<"I am confused, exiting with error"<<std::endl;
00214       throw -1;
00215     }
00216     OsiSolverInterface::ApplyCutsReturnCode acRc;
00217 
00218     //       linearModel.writeMps("test1");
00219     //       acRc = linearModel.applyCuts(cs);
00220     //       linearModel.writeMps("test2");
00221 
00222     // Print applyCuts return code
00223     std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00224     std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00225     std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
00226     <<" were inconsistent for this problem" <<std::endl;
00227     std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00228     std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00229     std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00230     std::cout << std::endl << std::endl;
00231   }
00232   //Modify the linearModel for feasibility pump
00233   //save objective
00234   int numcols = linearModel.getNumCols();
00235   double * saveObj = new double[numcols];
00236   CoinCopyN(linearModel.getObjCoefficients(), numcols, saveObj);
00237   //double * obj = new double[2*numIntCols];
00238   //int k = 0;
00239 
00240   //if ub is given add a "cutoff" row
00241   if(ub < 1e100) {
00242     //change bound on alpha
00243     linearModel.setColUpper(numcols-1,ub);
00244     //     addCutOff=1;
00245     //     CoinPackedVector * rowToAdd = (CoinPackedVector *) emptyCols[k];
00246     //     rowToAdd->insert(origNumCols-1,1.);
00247     //     colUb[k] = ub;
00248     //     colLb[k] = -linearModel.getInfinity();
00249   }
00250   for(int i = 0;i< numcols ;i++)
00251     if(!linearModel.isInteger(i))
00252       linearModel.setObjCoeff(i,0.);
00253 
00254   //done
00255   int nAddedCuts = 0;
00256 
00257   bool solved = 0;
00258   bool numericFailure = 0;
00259   int numScaleIncrease = 0;
00260   while((nMajorIt < maxIter) && (!solved && ( CoinCpuTime() + info.time < maxTime) && !numericFailure && numScaleIncrease < 5)) {
00261     //Change the objective of the MIP to get the closest point to rounding of NLP optimum
00262     for(int i = 0; i < numIntCols; i++) {
00263       linearModel.setObjCoeff(inds[i],1 - 2* nlp.getColSolution()[inds[i]]);
00264     }
00265 #ifdef COIN_HAS_CPX
00266     //        CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_NODELIM, nMaxNodes - nTotalNodes);
00267     //      CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 1);
00268     //      CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_HEURFREQ, 10);
00269     //      CPXsetdblparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_TILIM, maxTime - CoinCpuTime() + time);
00270 
00271 #else
00272     linearModel.initialSolve();
00273     CbcModel mip(linearModel);
00274     //        mip.writeMps("oa");
00275     CbcStrategyDefault defaultStrategy;
00276     mip.setStrategy(defaultStrategy);
00277     mip.solver()->messageHandler()->setLogLevel(0);
00278 
00279     //Add some heuristics to get feasible solutions
00280     mip.setMaximumSeconds(maxTime - CoinCpuTime() - info.time);
00281     mip.setMaximumSolutions(1);
00282 #endif
00283 
00284     bbTime -= CoinCpuTime();
00285 #ifdef COIN_HAS_CPX
00286 
00287     int mipstat = findGoodSolution(mip, maxTime - CoinCpuTime() - info.time, nTotalNodes);
00288 #else
00289     //  mip.branchAndBound();
00290 #endif
00291 
00292     bbTime += CoinCpuTime();
00293     //      mip.writeMps("oa1");
00294     OsiCuts cs;
00295     nMajorIt++;
00296 
00297 #ifdef COIN_HAS_CPX
00298 
00299     const double * colsol = mip.getColSolution();
00300     //int nNodes = 0;
00301     nTotalNodes += CPXgetnodecnt(cpxSi->getEnvironmentPtr(),cpxSi->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00302     //    int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00303     //   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_HEURFREQ, 0);
00304     if(mipstat == CPXMIP_INFEASIBLE) {
00305       provenInfeas=1;
00306       info.time += CoinCpuTime();
00307       for(int i = 0 ; i < numcols; i++)
00308         linearModel.setObjCoeff(i, saveObj[i]);
00309       return 1e75;
00310     }
00311     else if(mipstat == CPXMIP_TIME_LIM_FEAS || mipstat == CPXMIP_TIME_LIM_INFEAS) {
00312       provenInfeas= -1;
00313       info.time += CoinCpuTime();
00314       for(int i = 0 ; i < numcols; i++)
00315         linearModel.setObjCoeff(i, saveObj[i]);
00316       return 1e75;
00317     }
00318 #else
00319     if(mip.bestSolution() == NULL)
00320       break;
00321     if(mip.isNodeLimitReached()) {
00322       break;
00323     }
00324     //      if(!mip.isSolutionLimitReached())
00325     //          break;
00326     nTotalNodes += mip.getNodeCount();
00327     const double * colsol = mip.bestSolution();
00328 #endif
00329 
00330     for(int i = 0 ; i < numIntCols ; i++) {
00331       vals[i] = floor(colsol[inds[i]] + 0.5);
00332       vals[i] = max(mip.getColLower()[inds[i]],vals[i]);
00333       vals[i] = min(mip.getColUpper()[inds[i]],vals[i]);
00334     }
00335     nlpTime -= CoinCpuTime();
00336     double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
00337     nlpTime += CoinCpuTime();
00338     nNlpIterations += nlp.getIterationCount();
00339     if(!nlp.isProvenOptimal()) {
00340       std::cout<<"?????"<<std::endl;
00341       throw -1;
00342     }
00343     nlpTime += CoinCpuTime();
00344     if(dist < 1e-04)//do integer infeasibility check on variables
00345     {
00346       solved = 1;
00347       double norm_inf = DBL_MAX;
00348       for(int i = 0 ; i < nlp.getNumCols() && solved; i++)
00349       {
00350         if(nlp.isInteger(i)) {
00351           const double &value = nlp.getColSolution()[i];
00352           double IIf = fabs( value- floor(value + 0.5));
00353           norm_inf = min(norm_inf, IIf);
00354           if(fabs( value- floor(value + 0.5)) > 1e-06) {
00355             std::cout<<"Variable "<<i<<" has integer infeasibility : "<<IIf<<std::endl;
00356             solved=0;
00357             numericFailure = 1;
00358 #if 0
00359 
00360             numScaleIncrease++;
00361             nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
00362             for(int i = 0 ; i < numIntCols ; i++) {
00363               nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
00364               nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
00365             }
00366 #endif
00367 
00368           }
00369         }
00370       }
00371       std::cout<<"Found a solution with maximal integer infeasibility of "<<norm_inf<<std::endl;
00372     }
00373 
00374     OsiSolverInterface::ApplyCutsReturnCode acRc;
00375     std::cout<<"Cut violation :"<<cs.rowCut(0).violated(colsol)
00376     <<std::endl;
00377     //       cs.printCuts();
00378 
00379     linearModel.writeMps("test1");
00380     acRc = linearModel.applyCuts(cs);
00381     linearModel.writeMps("test2");
00382     nAddedCuts += cs.sizeRowCuts();
00383     // Print applyCuts return code
00384     std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00385     std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00386     std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
00387     <<" were inconsistent for this problem" <<std::endl;
00388     std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00389     std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00390     std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00391     std::cout << std::endl << std::endl;
00392 
00393 
00394     if(solved) {
00395       //Set warm start point to the last point found (which is feasible for this relaxation)
00396       nlp.setColSolution(nlp.getColSolution());
00397       nlp.setRowPrice(nlp.getRowPrice());
00398       nlp.solver()->enableWarmStart();
00399       //Resolve the NLP with fixed variables and original objective function
00400       for(int i = 0; i < numIntCols; i++) {
00401         nlp.setColLower(inds[i], vals[i]);
00402         nlp.setColUpper(inds[i], vals[i]);
00403       }
00404       //nlp.turnOnSolverOutput();
00405       nlp.initialSolve();
00406       if(nlp.isProvenOptimal()) {
00407         OsiCuts cs;
00408         nlp.getOuterApproximation(cs,1, NULL, true);
00409         linearModel.applyCuts(cs);
00410         nAddedCuts += cs.sizeRowCuts();
00411         std::cout<<" FP found easible solution of value "<<nlp.getObjValue()<<" in "<< CoinCpuTime() - BeginTimeGLOB<<" seconds, "<<nMajorIt<<" major iterations, took"
00412         <<nTotalNodes<<" nodes."<<std::endl;
00413         objValue = nlp.getObjValue();
00414         if (solution==NULL) solution = new double[numcols];
00415         CoinCopyN(nlp.getColSolution(), numcols, solution);
00416       }
00417       else {
00418         std::cout<<" FP converged in "<<CoinCpuTime()- BeginTimeGLOB<<" seconds, "<<nMajorIt<<" major iterations, took"
00419         <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
00420         std::cout<<"Increasing scaling of objective and restarting"<<std::endl;
00421         solved = 0;
00422         numScaleIncrease++;
00423         numericFailure = 1;
00424 #if 0
00425 
00426         nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
00427         for(int i = 0 ; i < numIntCols ; i++) {
00428           nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
00429           nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
00430         }
00431 #endif
00432 
00433       }
00434 
00435     }
00436 
00437   }
00438 
00439   if(!solved) {
00440     if(numericFailure)
00441       std::cout<<"FP aborted because of a numberical failure : "<<CoinCpuTime() - BeginTimeGLOB<<",  "<<nMajorIt<<" major iterations, took"
00442       <<nTotalNodes<<" nodes."<<std::endl;
00443     else
00444       std::cout<<"FP aborted on time limit elapsed time : "<<CoinCpuTime() - BeginTimeGLOB<<",  "<<nMajorIt<<" major iterations, took"
00445       <<nTotalNodes<<" nodes."<<std::endl;
00446     provenInfeas = -1;
00447 
00448   }
00449   std::cout<<"Nlp solve time : "<<nlpTime<<" B&B solve time : "<<bbTime<<std::endl;
00450   //    restore objective
00451   for(int i = 0 ; i < numcols; i++)
00452     linearModel.setObjCoeff(i, saveObj[i]);
00453   delete [] saveObj;
00454   info.time += CoinCpuTime();
00455   return objValue;
00456 }
00457 
00458 
00459 #if 0 // Old enhanced OA code without the iFP at beginning
00460 int main3 (int argc, char *argv[])
00461 {
00462   std::cout.precision(11);
00463   using namespace Ipopt;
00464 
00465   // Read in model using argv[1]
00466   char * pbName = new char[strlen(argv[1])+1];
00467   strcpy(pbName, argv[1]);
00468   BonminAmplInterface solver1(argv);
00469   solver1.turnOnSolverOutput();
00470   bool doFp=true;
00471 
00472 #ifdef COIN_HAS_CPX
00473 
00474   OsiCpxSolverInterface solver;
00475   OsiCpxSolverInterface &mip = solver;
00476 #else
00477 
00478   OsiClpSolverInterface solver;
00479 #endif
00480 
00481   solver.messageHandler()->setLogLevel(0);
00482 
00483 
00484   //Setup timers since an nlp is solved to extract the linear relaxation
00485   BeginTimeGLOB= CoinCpuTime();
00486   double bbTime = 0;
00487   double nlpTime = - BeginTimeGLOB;
00488 
00489   ResolutionInformation FP_infos;
00490   solver1.extractLinearRelaxation(solver, 1);
00491   nlpTime += CoinCpuTime();
00492 
00493 #ifdef EXTRACT_LIN_REL
00494 
00495   std::string linearName="Lin";
00496   linearName += pbName;
00497   solver.writeMpsNative(linearName.c_str(),NULL,NULL);
00498 #endif
00499 
00500   int numIntCols = 0;
00501   int * inds = new int[solver1.getNumCols()]; //indices of integer variables
00502   double * x = new double[solver1.getNumCols()]; //to store values of integer variables later on
00503   //int origNumCols = solver.getNumCols();
00504   //int origNumRows = solver.getNumRows();
00505   for(int i = 0; i < solver1.getNumCols(); i++) {
00506     if(solver1.isInteger(i)) {
00507       inds[numIntCols++] = i;
00508     }
00509   }
00510 
00511   double lb = solver1.getObjValue();
00512   double ub = DBL_MAX;
00513   //Limits of the procedure
00514   //  int nMaxNodes = 100000;
00515   params.maxTime_ = 2.*3600.;
00516 
00517   //counters for iterations, nodes,...
00518   int nMajorIt = 0;
00519   int nTotalNodes = 0;
00520   int nTimesFPCalled =0;
00521   double precision = 1e-02;
00522   //bool solved = 0;
00523   double firstIterationTime;
00524   solver.messageHandler()->setLogLevel(0);
00525   int numNotFound = 0;
00526   bool minutePassed = 0;
00527   while( (ub - lb) > precision && (CoinCpuTime() - BeginTimeGLOB < params.maxTime_) ) {
00528     solver.resolve();
00529 #ifndef COIN_HAS_CPX
00530 
00531     CbcModel mip(solver);
00532     //  solver.writeMps("oa");
00533     CbcStrategyDefault defaultStrategy(1,8,4);
00534     mip.setStrategy(defaultStrategy);
00535 
00536     mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
00537     mip.setMaximumSolutions(3);
00538     mip.setCutoff(ub);
00539     mip.setLogLevel(1);
00540     mip.solver()->messageHandler()->setLogLevel(0);
00541     //      mip.solver()->writeMps("test");
00542 #else
00543     //        CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_NODELIM, nMaxNodes - nTotalNodes);
00544     //   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 1);
00545     CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
00546     CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, ub);
00547 
00548 #endif
00549 
00550     bbTime -= CoinCpuTime();
00551 
00552     //      mip.solver()->writeMps("oa");
00553 
00554     //int mipstat = findGoodSolution(mip, minNodes,
00555     //                             nodeInterval, maxTime - CoinCpuTime() + time, nTotalNodes)
00556 
00557     mip.branchAndBound();
00558     bbTime += CoinCpuTime();
00559     if(nMajorIt==0) {
00560       firstIterationTime = bbTime;
00561     }
00562 #ifdef EXTRACT_LIN_REL
00563     std::string oa_name=linearName;
00564     std::ostringstream t;
00565     t<<nMajorIt;
00566     oa_name+=t.str();
00567 #endif
00568 
00569     OsiCuts cs;
00570 #ifdef COIN_HAS_CPX
00571 
00572     const double *colsol=mip.getColSolution();
00573     nTotalNodes += CPXgetnodecnt(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00574     double new_lb = 0;
00575     int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00576     //    if(mipstat ==
00577     if((mipstat != CPXMIP_OPTIMAL && mipstat != CPXMIP_OPTIMAL_TOL)) {
00578       int status = CPXgetbestobjval(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL), &new_lb);
00579       lb=max(lb,new_lb);
00580       if(status)
00581         throw CoinError("Error in getting CPLEX best bound","IpCbcOACutGenerator2","siBestObj");
00582     }
00583     else
00584       lb = mip.getObjValue();
00585 
00586     CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, mip.getInfinity());
00587 
00588 #else
00589     //      if(!mip.isSolutionLimitReached())
00590     //          break;
00591     nTotalNodes += mip.getNodeCount();
00592     const double * colsol = mip.bestSolution();
00593     lb = mip.getBestPossibleObjValue();
00594     if(mip.bestSolution() == NULL)
00595       break;
00596 #endif
00597 
00598     nMajorIt++;
00599 
00600     std::cout<<"Found new lower bound of : "<<lb<<std::endl;
00601     if ((ub - lb) < precision)
00602       break;
00603     for(int i = 0 ; i < numIntCols ; i++) {
00604       x[i] = max(solver.getColLower()[inds[i]],colsol[inds[i]]);
00605       x[i] = min(solver.getColUpper()[inds[i]],x[i]);
00606       solver1.setColLower(inds[i], x[i]);
00607       solver1.setColUpper(inds[i], x[i]);
00608     }
00609     solver1.turnOnSolverOutput();
00610     solver1.initialSolve();
00611     if(solver1.isProvenOptimal()) {
00612       std::cout<<pbName<<" OA found easible solution of value "
00613       <<solver1.getObjValue()<<" in "
00614       <<CoinCpuTime() -BeginTimeGLOB<<" seconds, "
00615       <<nMajorIt<<" major iterations, took"
00616       <<nTotalNodes<<" nodes."<<std::endl;
00617       if(ub > 1e100)
00618         std::cout<<"First solution "<< solver1.getObjValue()<<" after "
00619         <<CoinCpuTime() -BeginTimeGLOB<<" secs."<<std::endl;
00620       if(!minutePassed && (CoinCpuTime() -BeginTimeGLOB) > 61.) {
00621         minutePassed = true;
00622         std::cout<<"Solution at minute "<< ub<<" after "
00623         <<CoinCpuTime() -BeginTimeGLOB<<" secs."<<std::endl;
00624       }
00625 
00626       OsiCuts cs;
00627       solver1.getOuterApproximation(cs,1);
00628       OsiSolverInterface::ApplyCutsReturnCode acRc;
00629       ub = min (solver1.getObjValue(), ub);
00630       acRc = solver.applyCuts(cs);
00631       // Print applyCuts return code
00632       std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00633       std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00634       std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
00635       <<" were inconsistent for this problem" <<std::endl;
00636       std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00637       std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00638       std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00639       std::cout << std::endl << std::endl;
00640 
00641     }
00642     else if(solver1.isAbandoned() || solver1.isIterationLimitReached()) {
00643       TNLPSolver::UnsolvedError * E=solver1.newUnsolvedError(0,argv[1],"not solved");
00644       E->writeDiffFiles();
00645       solver1.turnOnSolverOutput();
00646       solver1.initialSolve();
00647 
00648       std::cerr<<"Error"<<std::endl;
00649       throw -1;
00650     }
00651 
00652     else {
00653       if (doFp)//Launch FP
00654       {
00655         //Reset NLP bounds
00656         for(int i = 0 ; i < numIntCols ; i++)
00657         {
00658           solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
00659           solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
00660         }
00661         ResolutionInformation infos;
00662         nTimesFPCalled++;
00663         int provenInfeas =0;
00664 
00665         //compute FP maxTime
00666         //            firstIterationTime = 1.;
00667         double fpTime = min( params.maxTime_ - CoinCpuTime() + BeginTimeGLOB,120.);
00668         int fpMaxIter = 30;
00669         ub = min(ub, FP(solver1, solver,numIntCols, inds, x,fpTime, fpMaxIter,  infos, ub*(1-1e-03), provenInfeas, solution));
00670         if(provenInfeas == 1)
00671           lb = mip.getInfinity();
00672         if(provenInfeas == -1)
00673         {
00674           numNotFound++;
00675           if(numNotFound==2)
00676             doFp = 0;
00677         }
00678         //            if(infos.n_iterations >= 10)
00679         //              doFp = 0;
00680         FP_infos+=infos;
00681       }
00682       else {
00683         std::cout<<"Adding feasibility cuts based on 1-norm of"
00684         <<"constraint satisfaction"<<std::endl;
00685         // solver1.getOuterApproximation(cs,1);
00686         solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs);
00687         OsiSolverInterface::ApplyCutsReturnCode acRc;
00688         acRc = solver.applyCuts(cs);
00689         // Print applyCuts return code
00690         std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00691         std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00692         std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
00693         <<" were inconsistent for this problem" <<std::endl;
00694         std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00695         std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00696         std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00697         std::cout << std::endl << std::endl;
00698       }
00699     }
00700 
00701   }
00702   double time = CoinCpuTime() - BeginTimeGLOB;
00703   std::cout<<pbName<<" FP driven OA found optimal solution of value "<<ub<<" (lower bound is "<<lb<<") in "<<time<<" seconds, "<<nMajorIt<<" major iterations, "
00704   <<nTimesFPCalled<<" calls to FP ( "
00705   <<FP_infos.n_iterations <<" it, and "<< FP_infos.time
00706   <<"sec) overall took "
00707   <<nTotalNodes + FP_infos.n_mip_nodes<<" nodes."<<std::endl;
00708   std::cout<<"Nlp solve time : "<<nlpTime + FP_infos.nlp_time<<" B&B solve time : "<<bbTime + FP_infos.mip_time<<std::endl;
00709   delete[] inds;
00710   delete[] x;
00711   delete [] pbName;
00712   return 0;
00713 }
00714 #endif
00715 
00716 
00717 
00718 
00720 int enhancedOA(AmplInterface & solver1, bool doFp,
00721                double *& solution)
00722 {
00723   bool nonConvex = 0;
00724   std::string pbName;
00725   solver1.getStrParam(OsiProbName, pbName);
00726 #ifdef COIN_HAS_CPX
00727 
00728   OsiCpxSolverInterface solver;
00729   OsiCpxSolverInterface &mip = solver;
00730 #else
00731 
00732   OsiClpSolverInterface solver;
00733 #endif
00734 
00735   solver.messageHandler()->setLogLevel(0);
00736 
00737 
00738   //Setup timers since an nlp is solved to extract the linear relaxation
00739   BeginTimeGLOB= CoinCpuTime();
00740   double bbTime = 0;
00741   double nlpTime = - BeginTimeGLOB;
00742 
00743   ResolutionInformation FP_infos;
00744   solver1.extractLinearRelaxation(solver, !nonConvex);
00745   //solver1.fpnlp()->setObjectiveScaling(1000* solver1.fpnlp()->getObjectiveScaling());
00746   nlpTime += CoinCpuTime();
00747   int numcols = solver.getNumCols();
00748   double * saveObj = new double[numcols];
00749   CoinCopyN(solver.getObjCoefficients(), numcols, saveObj);
00750 
00751 #ifdef EXTRACT_LIN_REL
00752   {
00753     std::string linearName="Lin";
00754     linearName += basename(pbName);
00755     solver.writeMpsNative(linearName.c_str(),NULL,NULL);
00756   }
00757 #endif
00758 
00759   int numIntCols = 0;
00760   int * inds = new int[solver1.getNumCols()]; //indices of integer variables
00761   double * x = new double[solver1.getNumCols()]; //to store values of integer variables later on
00762   //int origNumCols = solver.getNumCols();
00763   //int origNumRows = solver.getNumRows();
00764   for(int i = 0; i < solver1.getNumCols(); i++) {
00765     if(solver1.isInteger(i)) {
00766       inds[numIntCols++] = i;
00767       x[i]=0.5;
00768     }
00769   }
00770 
00771   double lb = solver1.getObjValue();
00772   double ub = DBL_MAX;
00773   //Limits of the procedure
00774   //  int nMaxNodes = 100000;
00775   double maxFpTime = 60.;//params.maxTime_+1;
00776 
00777 
00778   //counters for iterations, nodes,...
00779   int nMajorIt = 0;
00780   int nTotalNodes = 0;
00781   int nTimesFPCalled =0;
00782   int nTimesInitFPCalled =0;
00783   int nFPIterations =0;
00784   double precision = 1e-04;
00785   bool solved = 0;
00786   double firstIterationTime;
00787   solver.messageHandler()->setLogLevel(0);
00788   int numNotFound = 0;
00789 
00790   bool feasible = 1;
00791 
00792   double FPTime = 0.;
00793   
00794   if(doFp)
00795   {
00796   
00797   while(ub-lb> precision && feasible && ( CoinCpuTime() - BeginTimeGLOB < maxFpTime)) {
00798     nTimesInitFPCalled++;
00799     int nIt = 0;
00800     solved = 0;
00801     while(feasible & !solved &&
00802           (CoinCpuTime() - BeginTimeGLOB < maxFpTime)
00803           //   && nIt < 20
00804          ) {
00805       nIt++;
00806       nFPIterations++;
00807       for(int i = 0 ; i < numIntCols; i++) {
00808         solver.setObjCoeff(inds[i],1 - 2* solver1.getColSolution()[inds[i]]);
00809       }
00810       //change bound on alpha
00811       if(!nonConvex) {
00812         int sign = ub > 0 ? -1 : 1;
00813         solver.setColUpper(solver.getNumCols()-1,ub*(1 + sign *1e-03));
00814       }
00815 #ifndef COIN_HAS_CPX
00816       solver.initialSolve();
00817       solver.messageHandler()->setLogLevel(0);
00818 
00819       CbcStrategyDefault defaultStrategy;
00820       CbcModel mip(solver);
00821       mip.setStrategy(defaultStrategy);
00822       mip.solver()->messageHandler()->setLogLevel(0);
00823       mip.setLogLevel(0);
00824       mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() - BeginTimeGLOB);
00825       mip.setMaximumSolutions(1);
00826 #else
00827       //              CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_NODELIM, nMaxNodes - nTotalNodes);
00828     //  CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 3);
00829       CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM,
00830                      maxFpTime - CoinCpuTime() + BeginTimeGLOB);
00831       
00832 #endif
00833 
00834       bbTime -= CoinCpuTime();
00835 
00836 #ifdef COIN_HAS_CPX
00837 
00838       int mipstat = findGoodSolution(mip, maxFpTime - CoinCpuTime() - BeginTimeGLOB, nTotalNodes);
00839 #else
00840 
00841       mip.branchAndBound();
00842 #endif
00843 
00844       bbTime += CoinCpuTime();
00845       //      mip.solver()->writeMps("oa");
00846 
00847       OsiCuts cs;
00848 #ifdef COIN_HAS_CPX
00849 
00850       const double *colsol=mip.getColSolution();
00851       //          nTotalNodes += CPXgetnodecnt(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00852       //int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
00853       if(mipstat == CPXMIP_INFEASIBLE || mipstat == CPXMIP_INForUNBD) {
00854         std::cout<<"I found an Infeasible mip"<<std::endl;
00855         feasible = 0;
00856         break;
00857       }
00858 
00859 #else
00860       if(mip.bestSolution() == NULL)
00861         break;
00862       if(mip.isNodeLimitReached()) {
00863         break;
00864       }
00865       //                if(!mip.isSolutionLimitReached())
00866       //                        break;
00867       nTotalNodes += mip.getNodeCount();
00868       const double * colsol = mip.bestSolution();
00869 #endif
00870 
00871       nMajorIt++;
00872       for(int i = 0 ; i < numIntCols ; i++) {
00873         x[i] = max(solver1.getColLower()[inds[i]],colsol[inds[i]]);
00874         x[i] = min(solver1.getColUpper()[inds[i]],x[i]);
00875         //                              std::cout<<"Var "<<ind[i]<<" value "<<x[i]<<std::endl;
00876       }
00877       nlpTime -= CoinCpuTime();
00878       double dist = solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs, false, true);
00879       nlpTime += CoinCpuTime();
00880       std::cout<<"Dist : "<<dist<<std::endl;
00881       if(dist < 1e-04)//do integer infeasibility check on variables
00882       {
00883         solved = 1;
00884         for(int i = 0 ; i < solver1.getNumCols() && solved; i++)
00885         {
00886           if(solver1.isInteger(i)) {
00887             const double &value = solver1.getColSolution()[i];
00888             if(fabs( value- floor(value + 0.5)) > 1e-04) {
00889               solved = 0;
00890               std::cout<<"Variable "<<i<<" has integer infeasibility : "<<fabs( value- floor(value + 0.5))<<std::endl;
00891             }
00892           }
00893         }
00894         //          feasible = !solved;
00895       }
00896 
00897       OsiSolverInterface::ApplyCutsReturnCode acRc;
00898       acRc = solver.applyCuts(cs);
00899       // Print applyCuts return code
00900       std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00901       std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00902       std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
00903       <<" were inconsistent for this problem" <<std::endl;
00904       std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00905       std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00906       std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00907       std::cout << std::endl << std::endl;
00908     }
00909     double time = CoinCpuTime() - BeginTimeGLOB;
00910     if(solved) {
00911       //Resolve the NLP with fixed variables and original objective function
00912       for(int i = 0; i < numIntCols; i++) {
00913         solver1.setColLower(inds[i], x[i]);
00914         solver1.setColUpper(inds[i], x[i]);
00915       }
00916       solver1.turnOnSolverOutput();
00917       solver1.initialSolve();
00918       if(solver1.isProvenOptimal()) {
00919         ub = min(solver1.getObjValue(), ub);
00920         if (solution==NULL) solution = new double[numcols];
00921         CoinCopyN(solver1.getColSolution(), numcols, solution);
00922         std::cout<<pbName<<" FP found easible solution of value "
00923         <<solver1.getObjValue()<<" in "
00924         <<CoinCpuTime() - BeginTimeGLOB<<" seconds, "
00925         <<nMajorIt<<" major iterations, took"
00926         <<nTotalNodes<<" nodes."<<std::endl;
00927         if(nTimesInitFPCalled == 1)
00928           std::cout<<"FIRST sol found in : "<<time<<" seconds value is "<<solver1.getObjValue()
00929           <<" iterations "<<nFPIterations<<std::endl;
00930 
00931         solved = 1;
00932         OsiCuts cs;
00933         solver1.getOuterApproximation(cs,1, NULL, true);
00934         OsiSolverInterface::ApplyCutsReturnCode acRc;
00935         acRc = solver.applyCuts(cs);
00936         for(int i = 0 ; i < numIntCols ; i++) {
00937           solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
00938           solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
00939         }
00940         solver1.initialSolve();
00941       }
00942       else// may have cq problems add a feasibility cut
00943       {
00944         std::cout<<pbName<<" FP converged in "<<time<<" seconds, "<<nMajorIt<<" major iterations, took"
00945         <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
00946         CoinPackedVector v;
00947         double rhs = 1.;
00948         for(int i = 0 ; i < numIntCols ; i++) {
00949           v.insert(inds[i], -(2*x[i] - 1));
00950           rhs -= x[i];
00951         }
00952         OsiCuts cs;
00953         OsiRowCut cut;
00954         cut.setRow(v);
00955         cut.setLb(rhs);
00956         cut.setUb(1e300);
00957         cut.setGloballyValid();
00958         cs.insert(cut);
00959 
00960         OsiSolverInterface::ApplyCutsReturnCode acRc;
00961         acRc = solver.applyCuts(cs);
00962         // Print applyCuts return code
00963         std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
00964         std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
00965         std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
00966         <<" were inconsistent for this problem" <<std::endl;
00967         std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
00968         std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
00969         std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
00970         std::cout << std::endl << std::endl;
00971 
00972       }
00973     }
00974     else if (feasible) {
00975       //restart from NLP-relaxation optimum
00976       solver1.initialSolve();
00977       nIt = 0;
00978     }
00979     if(nonConvex)
00980       break;
00981   }
00982 
00983   //revert mip to original objective function
00984   for(int i = 0 ; i < numcols; i++)
00985     solver.setObjCoeff(i, saveObj[i]);
00986 
00987   solver.setColUpper(solver.getNumCols()-1,ub);
00988 
00989   FPTime = CoinCpuTime() - BeginTimeGLOB;
00990   std::cout<<"Best known sol after : "<<FPTime<<" is "<<ub<<std::endl;
00991 
00992 #ifdef EXTRACT_LIN_REL
00993 
00994   {
00995     std::string linearName="Lin1Min";
00996     linearName += basename(pbName);
00997     solver.writeMpsNative(linearName.c_str(),NULL,NULL);
00998   }
00999 #endif
01000 
01001   solver.setColUpper(solver.getNumCols()-1,ub);
01002   
01003 }
01004   while( (ub - lb) > precision && (CoinCpuTime() - BeginTimeGLOB < params.maxTime_) ) {
01005 #ifdef EXTRACT_LIN_REL
01006     {
01007       if(CoinCpuTime() - BeginTimeGLOB > 5. * 60.)
01008       {
01009         std::string linearName="Lin5min";
01010         linearName += basename(pbName);
01011         solver.writeMpsNative(linearName.c_str(),NULL,NULL);
01012         return 0;
01013       }
01014     }
01015 #endif
01016 
01017     solver.resolve();
01018 #ifndef COIN_HAS_CPX
01019 
01020     solver.messageHandler()->setLogLevel(0);
01021     CbcModel mip(solver);
01022     //  solver.writeMps("oa");
01023     CbcStrategyDefault defaultStrategy(1,8,4);
01024     mip.setStrategy(defaultStrategy);
01025 
01026     mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
01027     mip.setMaximumSolutions(3);
01028     mip.setCutoff(ub);
01029     mip.setLogLevel(1);
01030     mip.solver()->messageHandler()->setLogLevel(0);
01031     //      mip.solver()->writeMps("test");
01032 #else
01033     //        CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_NODELIM, nMaxNodes - nTotalNodes);
01034     //   CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 1);
01035     CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM, params.maxTime_ - CoinCpuTime() + BeginTimeGLOB);
01036     CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, ub);
01037 
01038 #endif
01039 
01040     bbTime -= CoinCpuTime();
01041 
01042     //      mip.solver()->writeMps("oa");
01043     //int mipstat = findGoodSolution(mip, minNodes,
01044     //                             nodeInterval, maxTime - CoinCpuTime() + time, nTotalNodes)
01045 
01046     mip.branchAndBound();
01047     bbTime += CoinCpuTime();
01048     if(nMajorIt==0) {
01049       firstIterationTime = bbTime;
01050     }
01051 
01052     OsiCuts cs;
01053 #ifdef COIN_HAS_CPX
01054 
01055     const double *colsol=mip.getColSolution();
01056     nTotalNodes += CPXgetnodecnt(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01057     double new_lb = 0;
01058     int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01059     //    if(mipstat ==
01060     if((mipstat != CPXMIP_OPTIMAL && mipstat != CPXMIP_OPTIMAL_TOL)) {
01061       int status = CPXgetbestobjval(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL), &new_lb);
01062       lb=max(lb,new_lb);
01063       if(status)
01064         throw CoinError("Error in getting CPLEX best bound","IpCbcOACutGenerator2","siBestObj");
01065     }
01066     else
01067       lb = mip.getObjValue();
01068 
01069     CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_CUTUP, mip.getInfinity());
01070 
01071 #else
01072     //      if(!mip.isSolutionLimitReached())
01073     //          break;
01074     nTotalNodes += mip.getNodeCount();
01075     const double * colsol = mip.bestSolution();
01076     lb = mip.getBestPossibleObjValue();
01077     if(mip.bestSolution() == NULL)
01078       break;
01079 #endif
01080 
01081     nMajorIt++;
01082 
01083     std::cout<<"Found new lower bound of : "<<lb<<std::endl;
01084     if ((ub - lb) < precision)
01085       break;
01086     bool same = 1;
01087     for(int i = 0 ; i < numIntCols ; i++) {
01088       double value;
01089       value = max(solver.getColLower()[inds[i]],colsol[inds[i]]);
01090       value = min(solver.getColUpper()[inds[i]],value);
01091       if(fabs(value - x[i])>0.0001) {
01092         same = 0;
01093       }
01094       x[i] = value;
01095       solver1.setColLower(inds[i], x[i]);
01096       solver1.setColUpper(inds[i], x[i]);
01097     }
01098     if(same) {
01099       std::cout<<"Converged on same solution"<<std::endl;
01100       lb += 1e10;
01101       break;
01102     }
01103     solver1.turnOnSolverOutput();
01104     solver1.initialSolve();
01105     if(solver1.isProvenOptimal()) {
01106       if (solution==NULL) solution = new double[numcols];
01107       CoinCopyN(solver1.getColSolution(), numcols, solution);
01108       std::cout<<pbName<<" OA found easible solution of value "
01109       <<solver1.getObjValue()<<" in "
01110       <<CoinCpuTime() -BeginTimeGLOB<<" seconds, "
01111       <<nMajorIt<<" major iterations, took"
01112       <<nTotalNodes<<" nodes."<<std::endl;
01113 
01114       OsiCuts cs;
01115       solver1.getOuterApproximation(cs,1, false, true);
01116       OsiSolverInterface::ApplyCutsReturnCode acRc;
01117       ub = min (solver1.getObjValue(), ub);
01118       acRc = solver.applyCuts(cs);
01119       // Print applyCuts return code
01120       std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01121       std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01122       std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
01123       <<" were inconsistent for this problem" <<std::endl;
01124       std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01125       std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01126       std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01127       std::cout << std::endl << std::endl;
01128 
01129     }
01130 
01131 
01132     else {
01133       if (doFp)//Launch FP
01134       {
01135         //Reset NLP bounds
01136         for(int i = 0 ; i < numIntCols ; i++)
01137         {
01138           solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
01139           solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
01140         }
01141         ResolutionInformation infos;
01142         nTimesFPCalled++;
01143         int provenInfeas =0;
01144 
01145         //compute FP maxTime
01146         //            firstIterationTime = 1.;
01147         double fpTime = min( params.maxTime_ - CoinCpuTime() + BeginTimeGLOB,120.);
01148         int fpMaxIter = 5;
01149         ub = min(ub, FP(solver1, solver,numIntCols, inds, x,fpTime, fpMaxIter,  infos, ub*(1-1e-03), provenInfeas, solution));
01150         if(provenInfeas == 1)
01151           lb = mip.getInfinity();
01152         if(provenInfeas == -1)
01153         {
01154           numNotFound++;
01155           //              if(numNotFound==10)
01156           //                doFp = 0;
01157         }
01158         //            if(infos.n_iterations >= 10)
01159         //              doFp = 0;
01160         FP_infos+=infos;
01161       }
01162       else {
01163         std::cout<<"Adding feasibility cuts based on 1-norm of"
01164         <<"constraint satisfaction"<<std::endl;
01165         //            solver1.getOuterApproximation(cs,1);
01166         solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs, false, true);
01167         OsiSolverInterface::ApplyCutsReturnCode acRc;
01168         acRc = solver.applyCuts(cs);
01169         // Print applyCuts return code
01170         std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01171         std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01172         std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
01173         <<" were inconsistent for this problem" <<std::endl;
01174         std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01175         std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01176         std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01177         std::cout << std::endl << std::endl;
01178       }
01179     }
01180 
01181   }
01182   std::string algoName=(doFp)?" enhanced OA":" OA";
01183   double time = CoinCpuTime() - BeginTimeGLOB;
01184   std::cout<<pbName<<algoName<<" found optimal solution of value "<<ub<<" (lower bound is "<<lb<<") in "<<time<<" seconds, "<<nMajorIt<<" major iterations, "<<std::endl;
01185   if(doFp)
01186   {
01187   std::cout<<nTimesInitFPCalled<<" calls to initial FP ( "
01188   <<nFPIterations <<" it, and "<< FPTime
01189   <<"sec)"
01190   <<std::endl
01191   <<nTimesFPCalled<<"subsequent calls to FP ( "
01192   <<FP_infos.n_iterations <<" it, and "<< FP_infos.time
01193   <<"sec)"<<std::endl;
01194   }
01195   std::cout<<"Nlp solve time : "<<nlpTime + FP_infos.nlp_time<<" B&B solve time : "<<bbTime + FP_infos.mip_time<<std::endl;
01196   delete[] inds;
01197   delete[] x;
01198   return 0;
01199 }
01200 
01201 
01202 
01204 int iteratedFP (AmplInterface& solver1, bool standAlone, 
01205                 double * &solution)
01206 {
01207   // Define a Solver which inherits from OsiClpsolverInterface -> OsiSolverInterface
01208 
01209   using namespace Ipopt;
01210   std::string pbName;
01211   solver1.getStrParam(OsiProbName,pbName);
01212 
01213 #ifdef COIN_HAS_CPX
01214 
01215   OsiCpxSolverInterface solver;
01216   OsiCpxSolverInterface &mip = solver;
01217 #else
01218 
01219   OsiClpSolverInterface solver;
01220 #endif
01221 
01222   solver.messageHandler()->setLogLevel(0);
01223 
01224   //Setup timers since an nlp is solved to extract the linear relaxation
01225   double beginTime= CoinCpuTime();
01226   double bbTime = 0;
01227   double nlpTime = - beginTime;
01228 
01229   solver1.extractLinearRelaxation(solver, 1);
01230   nlpTime += CoinCpuTime();
01231   //Add extra variables for linearizing norm-1
01232   //get the number of 0-1 variables
01233   int numIntCols = 0;
01234   int numcols = solver1.getNumCols();
01235   int * inds = new int[numcols]; //indices of integer variables
01236   double * x = new double[numcols]; //to store values of integer variables later on
01237 
01238   for(int i = 0; i < numcols; i++) {
01239     if(solver1.isInteger(i)) {
01240       inds[numIntCols++] = i;
01241     }
01242   }
01243   //int k = 0;
01244 
01245 
01246   for(int i = 0;i<solver.getNumCols();i++)
01247     solver.setObjCoeff(i,0.);
01248 
01249   //Limits of the procedure
01250   //int nMaxNodes = 100000;
01251 
01252   //counters for iterations, nodes,...
01253   int nMajorIt = 0;
01254   int nTotalNodes = 0;
01255 
01256 
01257   double lb=-DBL_MAX;
01258   double ub=DBL_MAX;
01259   bool solved = 0;
01260   while(ub-lb>1e-05) {
01261     double ub = 1e100;
01262     bool feasible = 1;
01263     while(feasible && ( CoinCpuTime() - beginTime < params.maxTime_) ) {
01264       int nIt = 0;
01265       while(feasible & !solved && (CoinCpuTime() - beginTime < params.maxTime_)) {
01266         nIt++;
01267         for(int i = 0 ; i < numIntCols; i++) {
01268           solver.setObjCoeff(inds[i],1 - 2* solver1.getColSolution()[inds[i]]);
01269         }
01270         //change bound on alpha
01271         solver.setColUpper(solver.getNumCols()-1,ub);
01272 
01273 #ifndef COIN_HAS_CPX
01274   
01275         solver.initialSolve();
01276         CbcModel mip(solver);
01277         CbcStrategyDefault defaultStrategy;
01278         mip.setStrategy(defaultStrategy);
01279         mip.solver()->messageHandler()->setLogLevel(0);
01280 
01281         //Add some heuristics to get feasible solutions
01282         mip.setMaximumSeconds(params.maxTime_ - CoinCpuTime() +beginTime);
01283         mip.setMaximumSolutions(1);
01284 #else
01285         //            CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_NODELIM, nMaxNodes - nTotalNodes);
01286         CPXsetintparam(mip.getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 3);
01287         CPXsetdblparam(mip.getEnvironmentPtr(),CPX_PARAM_TILIM,
01288                        params.maxTime_ - CoinCpuTime() + beginTime);
01289 
01290 #endif
01291 
01292         bbTime -= CoinCpuTime();
01293         mip.branchAndBound();
01294         bbTime += CoinCpuTime();
01295         //      mip.solver()->writeMps("oa");
01296 
01297         OsiCuts cs;
01298 #ifdef COIN_HAS_CPX
01299 
01300         const double *colsol=mip.getColSolution();
01301         nTotalNodes += CPXgetnodecnt(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01302         int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01303         if(mipstat == CPXMIP_INFEASIBLE || mipstat == CPXMIP_INForUNBD) {
01304           std::cout<<"Infeasible mip"<<std::endl;
01305           feasible = 0;
01306           break;
01307         }
01308 
01309 #else
01310         if(mip.bestSolution() == NULL)
01311           break;
01312         if(mip.isNodeLimitReached()) {
01313           break;
01314         }
01315         //              if(!mip.isSolutionLimitReached())
01316         //                      break;
01317         nTotalNodes += mip.getNodeCount();
01318         const double * colsol = mip.bestSolution();
01319 #endif
01320 
01321         nMajorIt++;
01322         for(int i = 0 ; i < numIntCols ; i++) {
01323           x[i] = max(solver1.getColLower()[inds[i]],colsol[inds[i]]);
01324           x[i] = min(solver1.getColUpper()[inds[i]],x[i]);
01325           //                            std::cout<<"Var "<<ind[i]<<" value "<<x[i]<<std::endl;
01326         }
01327         nlpTime -= CoinCpuTime();
01328         double dist = solver1.getFeasibilityOuterApproximation( numIntCols, x, inds, cs, false, true);
01329         nlpTime += CoinCpuTime();
01330         std::cout<<"Dist : "<<dist<<std::endl;
01331         if(dist < 1e-05)//do integer infeasibility check on variables
01332         {
01333           solved = 1;
01334           for(int i = 0 ; i < numcols && solved; i++)
01335           {
01336             if(solver1.isInteger(i)) {
01337               const double &value = solver1.getColSolution()[i];
01338               if(fabs( value- floor(value + 0.5)) > 1e-04) {
01339                 solved = 0;
01340                 std::cout<<"Variable "<<i<<" has integer infeasibility : "<<fabs( value- floor(value + 0.5))<<std::endl;
01341               }
01342             }
01343           }
01344           feasible = !solved;
01345         }
01346 
01347         OsiSolverInterface::ApplyCutsReturnCode acRc;
01348         acRc = solver.applyCuts(cs);
01349         // Print applyCuts return code
01350         std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01351         std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01352         std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
01353         <<" were inconsistent for this problem" <<std::endl;
01354         std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01355         std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01356         std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01357         std::cout << std::endl << std::endl;
01358       }
01359       double time = CoinCpuTime() - beginTime;
01360       if(solved) {
01361         //Resolve the NLP with fixed variables and original objective function
01362         for(int i = 0; i < numIntCols; i++) {
01363           solver1.setColLower(inds[i], x[i]);
01364           solver1.setColUpper(inds[i], x[i]);
01365         }
01366         solver1.turnOnSolverOutput();
01367         solver1.initialSolve();
01368         if(solver1.isProvenOptimal()) {
01369           ub = min(solver1.getObjValue() * (1-1e-4), ub);
01370           if (solution==NULL) solution = new double[numcols];
01371           CoinCopyN(solver1.getColSolution(), numcols, solution);
01372           std::cout<<pbName<<" FP found easible solution of value "
01373           <<solver1.getObjValue()<<" in "
01374           <<CoinCpuTime() - beginTime<<" seconds, "
01375           <<nMajorIt<<" major iterations, took"
01376           <<nTotalNodes<<" nodes."<<std::endl;
01377           solved = 1;
01378           OsiCuts cs;
01379           solver1.getOuterApproximation(cs,1, NULL, true);
01380           OsiSolverInterface::ApplyCutsReturnCode acRc;
01381           acRc = solver.applyCuts(cs);
01382           for(int i = 0 ; i < numIntCols ; i++) {
01383             solver1.setColLower(inds[i], solver.getColLower()[inds[i]]);
01384             solver1.setColUpper(inds[i], solver.getColUpper()[inds[i]]);
01385           }
01386           //  solver1.initialSolve();
01387         }
01388         else
01389           std::cout<<pbName<<" FP converged in "<<time<<" seconds, "<<nMajorIt<<" major iterations, took"
01390           <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
01391 
01392       }
01393       else if (feasible) {
01394         //restart from NLP-relaxation optimum
01395         solver1.initialSolve();
01396         nIt = 0;
01397       }
01398     }
01399       if(standAlone)
01400         break;
01401   }
01402   int returncode = 0;
01403   if(!standAlone)
01404   {
01405   //double time = CoinCpuTime() - beginTime;
01406   if(solved) {
01407     std::cout<<"iterated FP finished in : "<<CoinCpuTime() - beginTime
01408     <<", value of optimum "<<ub
01409     <<", "<<nMajorIt<<" major iterations, took"
01410     <<nTotalNodes<<" nodes."<<std::endl;
01411     returncode = 1;
01412   }
01413   else
01414     std::cout<<"iterated FP aborted on time limit elapsed time : "
01415     <<CoinCpuTime() - beginTime<<",  "<<nMajorIt
01416     <<" major iterations, took"
01417     <<nTotalNodes<<" nodes."<<std::endl;
01418     returncode = 0;
01419   }
01420   else returncode = 0;
01421   std::cout<<"Nlp solve time : "<<nlpTime
01422   <<" B&B solve time : "<<bbTime<<std::endl;
01423   delete[] inds;
01424   delete[] x;
01425   return returncode;
01426 }
01427 
01428 
01429 
01430 double FPGeneralIntegers(AmplInterface &nlp, OsiSolverInterface &linearModel,
01431                          int numIntCols, int * inds, double * vals, double maxTime,
01432                          ResolutionInformation& info, double ub, bool &provenInfeas)
01433 
01434 {
01435   provenInfeas=0;
01436 #ifdef COIN_HAS_CPX
01437 
01438   OsiCpxSolverInterface * cpxSi = dynamic_cast<OsiCpxSolverInterface *>
01439                                   (&linearModel);
01440   OsiCpxSolverInterface &mip = *cpxSi;
01441 #endif
01442 
01443   int &nMajorIt = info.n_iterations;
01444   int &nTotalNodes = info.n_mip_nodes;
01445   int &nNlpIterations = info.n_nlp_iterations;
01446   double& bbTime = info.mip_time;
01447   double& nlpTime = info.nlp_time;
01448   double &time = info.time = -CoinCpuTime();
01449   double objValue = DBL_MAX;
01450   {
01451     OsiCuts cs;
01452     //First get the closest point to the current integer (NLP-infeasible) optimum
01453     nlpTime -= CoinCpuTime();
01454     double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
01455     nlpTime += CoinCpuTime();
01456     if(dist < 1e-06)//do integer infeasibility check on variables
01457     {
01458       //Something wrong shouldn't be
01459       std::cout<<"Feasibility subproblem has objective 0 while problem was claimed infeasible before"<<std::endl
01460       <<"I am confused, exiting with error"<<std::endl;
01461       throw -1;
01462     }
01463 
01464     OsiSolverInterface::ApplyCutsReturnCode acRc;
01465 
01466     //       linearModel.writeMps("test1");
01467     //       acRc = linearModel.applyCuts(cs);
01468     //       linearModel.writeMps("test2");
01469 
01470     // Print applyCuts return code
01471     std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01472     std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01473     std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
01474     <<" were inconsistent for this problem" <<std::endl;
01475     std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01476     std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01477     std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01478     std::cout << std::endl << std::endl;
01479   }
01480   //Modify the linearModel for feasibility pump
01481   //save objective
01482   double * saveObj = new double[linearModel.getNumCols()];
01483   CoinCopyN(linearModel.getObjCoefficients(), linearModel.getNumCols(), saveObj);
01484   //Add extra variables for linearizing norm-1
01485   //get the number of 0-1 variables
01486 
01487   int origNumCols = linearModel.getNumCols();
01488   int origNumRows = linearModel.getNumRows();
01489 
01490   double * colLb = new double[2*numIntCols+1];
01491   double * colUb = new double[2*numIntCols+1];
01492   double * obj = new double[2*numIntCols];
01493   int k = 0;
01494   CoinPackedVectorBase ** emptyCols = new CoinPackedVectorBase*[2*numIntCols+1];
01495   for(int i = 0; i < nlp.getNumCols(); i++) {
01496     if(nlp.isInteger(i)) {
01497       colLb[k] = colLb[k+1] = 0.;
01498       colUb[k] = colUb[k+1] = nlp.getColUpper()[i] - nlp.getColLower()[i];
01499       obj[k] = obj[k+1] = 1.;
01500       emptyCols[k] = new CoinPackedVector;
01501       emptyCols[k + 1] = new CoinPackedVector;
01502       k+=2;
01503     }
01504   }
01505 
01506   linearModel.addCols(2*numIntCols, emptyCols, colLb, colUb, obj);
01507   k = 0;
01508   //add the constraint x_i - x^+_i + x^-_i = 0
01509   for(int i = 0; i < nlp.getNumCols(); i++) {
01510     if(nlp.isInteger(i)) {
01511       CoinPackedVector * rowToAdd = (CoinPackedVector *) emptyCols[k];
01512       rowToAdd->insert(i,1.);
01513       rowToAdd->insert(origNumCols + 2*k ,-1.);
01514       rowToAdd->insert(origNumCols + 2*k + 1,1.);
01515       colUb[k] = 0.;
01516       k++;
01517     }
01518   }
01519   //int addCutOff = 0;
01520   //if ub is given add a "cutoff" row
01521   if(ub < 1e100) {
01522     //change bound on alpha
01523     linearModel.setColUpper(origNumCols-1,ub);
01524     //     addCutOff=1;
01525     //     CoinPackedVector * rowToAdd = (CoinPackedVector *) emptyCols[k];
01526     //     rowToAdd->insert(origNumCols-1,1.);
01527     //     colUb[k] = ub;
01528     //     colLb[k] = -linearModel.getInfinity();
01529   }
01530   for(int i = 0;i<origNumCols;i++)
01531     linearModel.setObjCoeff(i,0.);
01532   linearModel.addRows(numIntCols, emptyCols, colLb, colUb);
01533   //done
01534   int nAddedCuts = 0;
01535 
01536   bool solved = 0;
01537   bool numericFailure = 0;
01538   int numScaleIncrease = 0;
01539   while(!solved && (time + CoinCpuTime() < maxTime) && !numericFailure && numScaleIncrease < 5 ) {
01540     for(int i = 0; i < numIntCols; i++) {
01541       linearModel.setRowBounds(origNumRows + i, nlp.getColSolution()[inds[i]], nlp.getColSolution()[inds[i]]);
01542     }
01543 #ifdef COIN_HAS_CPX
01544     //        CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_NODELIM, nMaxNodes - nTotalNodes);
01545     CPXsetintparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_INTSOLLIM, 1);
01546     CPXsetdblparam(cpxSi->getEnvironmentPtr(),CPX_PARAM_TILIM, maxTime - CoinCpuTime() + time);
01547 
01548 #else
01549 
01550     linearModel.initialSolve();
01551     CbcModel mip(linearModel);
01552     //        mip.writeMps("oa");
01553     CbcStrategyDefault defaultStrategy;
01554     mip.setStrategy(defaultStrategy);
01555     mip.solver()->messageHandler()->setLogLevel(0);
01556 
01557     mip.setMaximumSeconds(maxTime - CoinCpuTime() + time);
01558     mip.setMaximumSolutions(1);
01559 #endif
01560 
01561     bbTime -= CoinCpuTime();
01562     mip.branchAndBound();
01563     bbTime += CoinCpuTime();
01564     //      mip.writeMps("oa1");
01565     OsiCuts cs;
01566     nMajorIt++;
01567 
01568 #ifdef COIN_HAS_CPX
01569 
01570     const double * colsol = mip.getColSolution();
01571     //int nNodes = 0;
01572     nTotalNodes += CPXgetnodecnt(cpxSi->getEnvironmentPtr(),cpxSi->getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01573     int mipstat = CPXgetstat(mip.getEnvironmentPtr(),mip.getLpPtr(OsiCpxSolverInterface::KEEPCACHED_ALL));
01574     if(mipstat == CPXMIP_INFEASIBLE) {
01575       provenInfeas=1;
01576       return 1e75;
01577     }
01578 #else
01579     if(mip.bestSolution() == NULL)
01580       break;
01581     if(mip.isNodeLimitReached()) {
01582       break;
01583     }
01584     //      if(!mip.isSolutionLimitReached())
01585     //          break;
01586     nTotalNodes += mip.getNodeCount();
01587     const double * colsol = mip.bestSolution();
01588 #endif
01589 
01590     for(int i = 0 ; i < numIntCols ; i++) {
01591       vals[i] = floor(colsol[inds[i]] + 0.5);
01592       vals[i] = max(mip.getColLower()[inds[i]],vals[i]);
01593       vals[i] = min(mip.getColUpper()[inds[i]],vals[i]);
01594     }
01595     nlpTime -= CoinCpuTime();
01596     double dist = nlp.getFeasibilityOuterApproximation( numIntCols, vals, inds, cs, false, true);
01597     nlpTime += CoinCpuTime();
01598     nNlpIterations += nlp.getIterationCount();
01599     if(!nlp.isProvenOptimal()) {
01600       std::cout<<"?????"<<std::endl;
01601       throw -1;
01602     }
01603     nlpTime += CoinCpuTime();
01604     if(dist < 1e-04)//do integer infeasibility check on variables
01605     {
01606       solved = 1;
01607       double norm_inf = DBL_MAX;
01608       for(int i = 0 ; i < nlp.getNumCols() && solved; i++)
01609       {
01610         if(nlp.isInteger(i)) {
01611           const double &value = nlp.getColSolution()[i];
01612           double IIf = fabs( value- floor(value + 0.5));
01613           norm_inf = min(norm_inf, IIf);
01614           if(fabs( value- floor(value + 0.5)) > 1e-04) {
01615             std::cout<<"Variable "<<i<<" has integer infeasibility : "<<IIf<<std::endl;
01616             solved=0;
01617             numericFailure = 1;
01618 #if 0
01619 
01620             numScaleIncrease++;
01621             nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
01622             for(int i = 0 ; i < numIntCols ; i++) {
01623               nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
01624               nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
01625             }
01626 #endif
01627 
01628           }
01629         }
01630       }
01631       std::cout<<"Found a solution with maximal integer infeasibility of "<<norm_inf<<std::endl;
01632     }
01633 
01634     OsiSolverInterface::ApplyCutsReturnCode acRc;
01635     std::cout<<"Cut violation :"<<cs.rowCut(0).violated(colsol)
01636     <<std::endl;
01637     //       cs.printCuts();
01638 
01639     linearModel.writeMps("test1");
01640     acRc = linearModel.applyCuts(cs);
01641     linearModel.writeMps("test2");
01642     nAddedCuts += cs.sizeRowCuts();
01643     // Print applyCuts return code
01644     std::cout <<cs.sizeCuts() <<" cuts were generated" <<std::endl;
01645     std::cout <<"  " <<acRc.getNumInconsistent() <<" were inconsistent" <<std::endl;
01646     std::cout <<"  " <<acRc.getNumInconsistentWrtIntegerModel()
01647     <<" were inconsistent for this problem" <<std::endl;
01648     std::cout <<"  " <<acRc.getNumInfeasible() <<" were infeasible" <<std::endl;
01649     std::cout <<"  " <<acRc.getNumIneffective() <<" were ineffective" <<std::endl;
01650     std::cout <<"  " <<acRc.getNumApplied() <<" were applied" <<std::endl;
01651     std::cout << std::endl << std::endl;
01652 
01653 
01654     if(solved) {
01655       //Set warm start point to the last point found (which is feasible for this relaxation)
01656       nlp.setColSolution(nlp.getColSolution());
01657       nlp.setRowPrice(nlp.getRowPrice());
01658       nlp.solver()->enableWarmStart();
01659       //Resolve the NLP with fixed variables and original objective function
01660       for(int i = 0; i < numIntCols; i++) {
01661         nlp.setColLower(inds[i], vals[i]);
01662         nlp.setColUpper(inds[i], vals[i]);
01663       }
01664       //nlp.turnOnSolverOutput();
01665       nlp.initialSolve();
01666       if(nlp.isProvenOptimal()) {
01667         OsiCuts cs;
01668         nlp.getOuterApproximation(cs,1, NULL, true);
01669         linearModel.applyCuts(cs);
01670         nAddedCuts += cs.sizeRowCuts();
01671         std::cout<<" FP found easible solution of value "<<nlp.getObjValue()<<" in "<<time + CoinCpuTime()<<" seconds, "<<nMajorIt<<" major iterations, took"
01672         <<nTotalNodes<<" nodes."<<std::endl;
01673         objValue = nlp.getObjValue();
01674       }
01675       else {
01676         std::cout<<" FP converged in "<<time<<" seconds, "<<nMajorIt<<" major iterations, took"
01677         <<nTotalNodes<<" nodes, but solution seems infeasible"<<std::endl;
01678         std::cout<<"Increasing scaling of objective and restarting"<<std::endl;
01679         solved = 0;
01680         numScaleIncrease++;
01681         numericFailure = 1;
01682 #if 0
01683 
01684         nlp.fpnlp()->setObjectiveScaling(10* nlp.fpnlp()->getObjectiveScaling());
01685         for(int i = 0 ; i < numIntCols ; i++) {
01686           nlp.setColLower(inds[i], linearModel.getColLower()[inds[i]]);
01687           nlp.setColUpper(inds[i], linearModel.getColUpper()[inds[i]]);
01688         }
01689 #endif
01690 
01691       }
01692 
01693     }
01694 
01695   }
01696   time+=CoinCpuTime();
01697 
01698   if(!solved) {
01699     if(numericFailure)
01700       std::cout<<"FP aborted because of a numberical failure : "<<time<<",  "<<nMajorIt<<" major iterations, took"
01701       <<nTotalNodes<<" nodes."<<std::endl;
01702     else
01703       std::cout<<"FP aborted on time limit elapsed time : "<<time<<",  "<<nMajorIt<<" major iterations, took"
01704       <<nTotalNodes<<" nodes."<<std::endl;
01705 
01706   }
01707   std::cout<<"Nlp solve time : "<<nlpTime<<" B&B solve time : "<<bbTime<<std::endl;
01708 
01709 
01710   //erase the extra columns and extra rows added to linearModel and resore the objective function
01711   int * colToDelete = new int[2*numIntCols];
01712   int * rowToDelete = new int[numIntCols];
01713   for(int i = 0 ; i < numIntCols ; i++) {
01714     colToDelete[2*i] = origNumCols + 2*i;
01715     colToDelete[2*i +1] = origNumCols + 2*i + 1;
01716     rowToDelete[i] = origNumRows + i;
01717   }
01718   //   if(addCutOff)
01719   //     rowToDelete[numIntCols] = linearModel.getNumCols();
01720   linearModel.deleteCols(2*numIntCols, colToDelete);
01721   linearModel.deleteRows(numIntCols, rowToDelete);
01722   //Check that size are correct
01723   //  mip.writeMps("oa2");
01724   assert(linearModel.getNumCols()==origNumCols);
01725   assert(linearModel.getNumRows()==origNumRows + nAddedCuts);
01726   //    restore objective
01727   for(int i = 0 ; i < origNumCols; i++)
01728     linearModel.setObjCoeff(i, saveObj[i]);
01729   delete [] saveObj;
01730   return objValue;
01731 }

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