Ignore:
Timestamp:
Sep 29, 2006 4:48:36 PM (13 years ago)
Author:
pbonami
Message:

astyled the devel branch

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Bonmin/src/OaInterface/BonOACutGenerator2.cpp

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