Changeset 264 for trunk


Ignore:
Timestamp:
Mar 6, 2006 9:37:11 AM (14 years ago)
Author:
forrest
Message:

for nonlinear

Location:
trunk
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/CbcHeuristic.cpp

    r246 r264  
    1717#include "CbcStrategy.hpp"
    1818#include "CglPreProcess.hpp"
     19#include "OsiAuxInfo.hpp"
    1920
    2021// Default Constructor
     
    860861}
    861862
     863// Default Constructor
     864CbcSerendipity::CbcSerendipity()
     865  :CbcHeuristic()
     866{
     867}
     868
     869// Constructor from model
     870CbcSerendipity::CbcSerendipity(CbcModel & model)
     871  :CbcHeuristic(model)
     872{
     873}
     874
     875// Destructor
     876CbcSerendipity::~CbcSerendipity ()
     877{
     878}
     879
     880// Clone
     881CbcHeuristic *
     882CbcSerendipity::clone() const
     883{
     884  return new CbcSerendipity(*this);
     885}
     886
     887// Copy constructor
     888CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
     889:
     890  CbcHeuristic(rhs)
     891{
     892}
     893
     894// Returns 1 if solution, 0 if not
     895int
     896CbcSerendipity::solution(double & solutionValue,
     897                         double * betterSolution)
     898{
     899  if (!model_)
     900    return 0;
     901  // get information on solver type
     902  OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
     903  OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
     904  if (auxiliaryInfo)
     905    return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
     906  else
     907    return 0;
     908}
     909// update model
     910void CbcSerendipity::setModel(CbcModel * model)
     911{
     912  model_ = model;
     913}
     914// Resets stuff if model changes
     915void
     916CbcSerendipity::resetModel(CbcModel * model)
     917{
     918  model_ = model;
     919}
    862920 
  • trunk/CbcMessage.cpp

    r262 r264  
    2222  {CBC_EVENT,27,1,"Exiting on user event"},
    2323  {CBC_SOLUTION,4,1,"Integer solution of %g found after %d iterations and %d nodes"},
     24  {CBC_SOLUTION2,33,1,"Integer solution of %g found (by alternate solver) after %d iterations and %d nodes"},
    2425  {CBC_END,5,1,"Partial search - best objective %g (best possible %g), took %d iterations and %d nodes"},
    2526  {CBC_INFEAS,6,1,"The LP relaxation is infeasible or too expensive"},
  • trunk/CbcModel.cpp

    r262 r264  
    99//#define CHECK_CUT_COUNTS
    1010//#define CHECK_NODE_FULL
     11//#define NODE_LOG
    1112#include <cassert>
    1213#include <cmath>
     
    2021
    2122#include "OsiSolverInterface.hpp"
     23#include "OsiAuxInfo.hpp"
    2224#include "OsiSolverBranch.hpp"
    2325#include "CoinWarmStartBasis.hpp"
     
    422424#ifdef COIN_USE_CLP
    423425  ClpEventHandler * eventHandler=NULL;
    424  {
    425    OsiClpSolverInterface * clpSolver
    426      = dynamic_cast<OsiClpSolverInterface *> (solver_);
    427    if (clpSolver) {
    428      ClpSimplex * clpSimplex = clpSolver->getModelPtr();
    429      eventHandler = clpSimplex->eventHandler();
    430    }
    431  }
     426  {
     427    OsiClpSolverInterface * clpSolver
     428      = dynamic_cast<OsiClpSolverInterface *> (solver_);
     429    if (clpSolver) {
     430      ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     431      eventHandler = clpSimplex->eventHandler();
     432    }
     433  }
    432434#endif
     435 
    433436  if (!nodeCompare_)
    434437    nodeCompare_=new CbcCompareDefault();;
     
    503506*/
    504507  synchronizeModel() ;
     508  assert (!solverCharacteristics_);
     509  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
     510  if (solverCharacteristics) {
     511    solverCharacteristics_ = solverCharacteristics;
     512  } else {
     513    // replace in solver
     514    OsiBabSolver defaultC;
     515    solver_->setAuxiliaryInfo(&defaultC);
     516    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
     517  }
     518  solverCharacteristics_->setSolver(solver_);
    505519  // Set so we can tell we are in initial phase in resolve
    506520  continuousObjective_ = -COIN_DBL_MAX ;
     
    512526  continuous relaxation.
    513527*/
    514   bool feasible = resolve() ;
     528  bool feasible;
     529  // If NLP then we assume already solved outside branchAndbound
     530  if (!solverCharacteristics_->solverType()) {
     531    feasible=resolve() ;
     532  } else {
     533    // pick up given status
     534    feasible = (solver_->isProvenOptimal() &&
     535                !solver_->isDualObjectiveLimitReached()) ;
     536  }
    515537  if (problemFeasibility_->feasible(this,0)<0) {
    516538    feasible=false; // pretend infeasible
     
    602624  installed.
    603625*/
    604   analyzeObjective() ;
     626  if(solverCharacteristics_->reducedCostsAccurate())
     627    analyzeObjective() ;
    605628/*
    606629  Set up for cut generation. addedCuts_ holds the cuts which are relevant for
     
    690713          object_[iObject]->infeasibility(preferredWay) ;
    691714      if (infeasibility ) numberUnsatisfied++ ; }
    692     if (numberUnsatisfied) {
    693       feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,NULL);
    694     }
     715    // replace solverType
     716    if(solverCharacteristics_->tryCuts())  {
     717      if (numberUnsatisfied)   {
     718        feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
     719                                 NULL);
     720      } else if (solverCharacteristics_->solutionAddsCuts()) {
     721        // may generate cuts and turn the solution
     722        //to an infeasible one
     723        feasible = solveWithCuts(cuts, 1,
     724                                 NULL);
     725#if 0
     726        currentNumberCuts_ = cuts.sizeRowCuts();
     727        if (currentNumberCuts_ >= maximumNumberCuts_) {
     728          maximumNumberCuts_ = currentNumberCuts;
     729          delete [] addedCuts_;
     730          addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
     731        }
     732#endif
     733      }
     734    }
     735    // check extra info on feasibility
     736    if (!solverCharacteristics_->mipFeasible())
     737      feasible = false;
    695738  }
    696739  // make cut generators less aggressive
     
    720763  numberFixedNow_=0;
    721764  int numberIterationsAtContinuous = numberIterations_;
    722   if (feasible)
    723   { newNode = new CbcNode ;
    724     newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
     765  //solverCharacteristics_->setSolver(solver_);
     766  if (feasible) {
     767    newNode = new CbcNode ;
     768    double newObjValue = direction*solver_->getObjValue();
     769    if (newObjValue!=solverCharacteristics_->mipBound()) {
     770      newObjValue = CoinMin(newObjValue,solverCharacteristics_->mipBound());
     771      solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
     772    }
     773    newNode->setObjectiveValue(newObjValue);
    725774    anyAction = -1 ;
    726775    // To make depth available we may need a fake node
     
    767816        feasible=false; // pretend infeasible
    768817      }
    769       //int nfix=
    770818      reducedCostFix() ;
    771       //if (nfix)
    772       //printf("%d fixed after resolve root\n",nfix);
    773         resolved = true ;
     819      resolved = true ;
    774820#       ifdef CBC_DEBUG
    775         printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
    776                   solver_->getObjValue(),
    777                   solver_->getNumRows()) ;
     821      printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
     822             solver_->getObjValue(),
     823             solver_->getNumRows()) ;
    778824#       endif
    779         if (!feasible) anyAction = -2 ; }
     825      if (!feasible) anyAction = -2 ; }
    780826      if (anyAction == -2||newNode->objectiveValue() >= cutoff)
    781827      { delete newNode ;
     
    11111157      currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_ ;
    11121158      int saveNumber = numberIterations_;
    1113       feasible = solveWithCuts(cuts,maximumCutPasses_,node);
     1159      if(solverCharacteristics_->solutionAddsCuts()) {
     1160        feasible=resolve();
     1161        if (feasible) {
     1162          int iObject ;
     1163          int preferredWay ;
     1164          int numberUnsatisfied = 0 ;
     1165          memcpy(currentSolution_,solver_->getColSolution(),
     1166                 numberColumns*sizeof(double)) ;
     1167         
     1168          for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
     1169            double infeasibility =
     1170              object_[iObject]->infeasibility(preferredWay) ;
     1171            if (infeasibility ) numberUnsatisfied++ ;
     1172          }
     1173          if (numberUnsatisfied)   {
     1174            feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
     1175                                     NULL);
     1176          } else {
     1177            // may generate cuts and turn the solution
     1178            //to an infeasible one
     1179            feasible = solveWithCuts(cuts, 1,
     1180                                     NULL);
     1181#if 0
     1182            currentNumberCuts_ = cuts.sizeRowCuts();
     1183            if (currentNumberCuts_ >= maximumNumberCuts_) {
     1184              maximumNumberCuts_ = currentNumberCuts;
     1185              delete [] addedCuts_;
     1186              addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
     1187            }
     1188#endif
     1189          }
     1190          // check extra info on feasibility
     1191          if (!solverCharacteristics_->mipFeasible())
     1192            feasible = false;
     1193        }
     1194      } else {
     1195        // normal
     1196        feasible = solveWithCuts(cuts,maximumCutPasses_,node);
     1197      }
    11141198      if ((specialOptions_&1)!=0&&onOptimalPath) {
    11151199        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
     
    11621246        }
    11631247        bool checkingNode=false;
    1164         if (feasible)
    1165         { newNode = new CbcNode ;
    1166           newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
     1248        if (feasible) {
     1249          newNode = new CbcNode ;
     1250          double newObjValue = direction*solver_->getObjValue();
     1251          if (newObjValue!=solverCharacteristics_->mipBound()) {
     1252            newObjValue = CoinMin(newObjValue,solverCharacteristics_->mipBound());
     1253            solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
     1254          }
     1255          newNode->setObjectiveValue(newObjValue);
    11671256          anyAction =-1 ;
    11681257          resolved = false ;
     
    12061295              { newNode->setObjectiveValue(direction*
    12071296                                           solver_->getObjValue()) ;
    1208               //int nfix=
    1209               reducedCostFix() ;
    1210               //if (nfix)
    1211               //printf("%d fixed after resolve\n",nfix);
     1297                reducedCostFix() ;
    12121298              if (newNode->objectiveValue() >= getCutoff())
    12131299                anyAction=-2;
     
    16481734  public methods.
    16491735*/
    1650   if (bestSolution_)
     1736  if (bestSolution_&&solverCharacteristics_->solverType()<2)
    16511737  { setCutoff(1.0e50) ; // As best solution should be worse than cutoff
    16521738    phase_=5;
     
    16721758  delete [] addedCuts_ ;
    16731759  addedCuts_ = NULL ;
     1760  // Get rid of characteristics
     1761  solverCharacteristics_=NULL;
    16741762  if (continuousSolver_)
    16751763  { delete continuousSolver_ ;
     
    17891877*/
    17901878  synchronizeModel() ;
     1879  assert (!solverCharacteristics_);
     1880  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
     1881  if (solverCharacteristics) {
     1882    solverCharacteristics_ = solverCharacteristics;
     1883  } else {
     1884    // replace in solver
     1885    OsiBabSolver defaultC;
     1886    solver_->setAuxiliaryInfo(&defaultC);
     1887    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
     1888  }
     1889  solverCharacteristics_->setSolver(solver_);
    17911890  // Set so we can tell we are in initial phase in resolve
    17921891  continuousObjective_ = -COIN_DBL_MAX ;
     
    18881987  installed.
    18891988*/
    1890   analyzeObjective() ;
     1989  if(solverCharacteristics_->reducedCostsAccurate())
     1990    analyzeObjective() ;
    18911991/*
    18921992  Set up for cut generation. addedCuts_ holds the cuts which are relevant for
     
    20052105  numberFixedNow_=0;
    20062106  int numberIterationsAtContinuous = numberIterations_;
    2007   if (feasible)
    2008   { newNode = new CbcNode ;
    2009     newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
     2107  if (feasible) {
     2108    newNode = new CbcNode ;
     2109    double newObjValue = direction*solver_->getObjValue();
     2110    if (newObjValue!=solverCharacteristics_->mipBound()) {
     2111      newObjValue = CoinMin(newObjValue,solverCharacteristics_->mipBound());
     2112      solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
     2113    }
     2114    newNode->setObjectiveValue(newObjValue);
    20102115    anyAction = -1 ;
    20112116    // To make depth available we may need a fake node
     
    20522157        feasible=false; // pretend infeasible
    20532158      }
    2054       //int nfix=
    20552159      reducedCostFix() ;
    2056       //if (nfix)
    2057       //printf("%d fixed after resolve root\n",nfix);
    2058         resolved = true ;
     2160      resolved = true ;
    20592161#       ifdef CBC_DEBUG
    2060         printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
    2061                   solver_->getObjValue(),
    2062                   solver_->getNumRows()) ;
     2162      printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
     2163             solver_->getObjValue(),
     2164             solver_->getNumRows()) ;
    20632165#       endif
    2064         if (!feasible) anyAction = -2 ; }
     2166      if (!feasible) anyAction = -2 ; }
    20652167      if (anyAction == -2||newNode->objectiveValue() >= cutoff)
    20662168      { delete newNode ;
     
    25852687  delete [] addedCuts_ ;
    25862688  addedCuts_ = NULL ;
     2689  // Get rid of characteristics
     2690  solverCharacteristics_=NULL;
    25872691  if (continuousSolver_) {
    25882692    delete continuousSolver_ ;
     
    27292833    OsiCuts cuts = OsiCuts() ;
    27302834    int saveNumber = numberIterations_;
    2731     bool feasible = solveWithCuts(cuts,maximumCutPasses_,node);
     2835    if(solverCharacteristics_->solutionAddsCuts()) {
     2836      feasible=resolve();
     2837      if (feasible) {
     2838        int iObject ;
     2839        int preferredWay ;
     2840        int numberUnsatisfied = 0 ;
     2841        memcpy(currentSolution_,solver_->getColSolution(),
     2842               numberColumns*sizeof(double)) ;
     2843       
     2844        for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
     2845          double infeasibility =
     2846            object_[iObject]->infeasibility(preferredWay) ;
     2847          if (infeasibility ) numberUnsatisfied++ ;
     2848        }
     2849        if (numberUnsatisfied)   {
     2850          feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
     2851                                   NULL);
     2852        } else {
     2853          // may generate cuts and turn the solution
     2854          //to an infeasible one
     2855          feasible = solveWithCuts(cuts, 1,
     2856                                   NULL);
     2857          currentNumberCuts_ = cuts.sizeRowCuts();
     2858            if (currentNumberCuts_ >= maximumNumberCuts_) {
     2859              maximumNumberCuts_ = currentNumberCuts;
     2860              delete [] addedCuts_;
     2861              addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
     2862            }
     2863        }
     2864        // check extra info on feasibility
     2865        if (!solverCharacteristics_->mipFeasible())
     2866          feasible = false;
     2867      }
     2868    } else {
     2869      // normal
     2870      feasible = solveWithCuts(cuts,maximumCutPasses_,node);
     2871    }
    27322872    if (statistics_) {
    27332873      assert (numberNodes2_);
     
    27642904      if (feasible) {
    27652905        newNode = new CbcNode ;
    2766         newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
     2906        double newObjValue = direction*solver_->getObjValue();
     2907        if (newObjValue!=solverCharacteristics_->mipBound()) {
     2908          newObjValue = CoinMin(newObjValue,solverCharacteristics_->mipBound());
     2909          solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
     2910        }
     2911        newNode->setObjectiveValue(newObjValue);
    27672912        anyAction =-1 ;
    27682913        resolved = false ;
     
    28032948            if (feasible) {
    28042949              newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
    2805               //int nfix=
    28062950              reducedCostFix() ;
    2807               //if (nfix)
    2808               //printf("%d fixed after resolve\n",nfix);
    28092951              if (newNode->objectiveValue() >= getCutoff())
    28102952                anyAction=-2;
     
    32853427{
    32863428  assert (solver_);
     3429  assert (!solverCharacteristics_);
     3430  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
     3431  if (solverCharacteristics) {
     3432    solverCharacteristics_ = solverCharacteristics;
     3433  } else {
     3434    // replace in solver
     3435    OsiBabSolver defaultC;
     3436    solver_->setAuxiliaryInfo(&defaultC);
     3437    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
     3438  }
     3439  solverCharacteristics_->setSolver(solver_);
    32873440  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
    32883441  solver_->initialSolve();
     
    32963449                                             solver_->getNumCols());
    32973450  setPointers(solver_);
     3451  solverCharacteristics_ = NULL;
    32983452}
    32993453
     
    34443598  strongInfo_[1]=0;
    34453599  strongInfo_[2]=0;
     3600  solverCharacteristics_ = NULL;
    34463601  nodeCompare_=new CbcCompareDefault();;
    34473602  problemFeasibility_=new CbcFeasibilityBase();
     
    35643719  strongInfo_[1]=0;
    35653720  strongInfo_[2]=0;
     3721  solverCharacteristics_ = NULL;
    35663722
    35673723  nodeCompare_=new CbcCompareDefault();;
     
    35763732  messages_ = CbcMessage();
    35773733  solver_ = rhs.clone();
    3578   referenceSolver_ = rhs.clone();
     3734  referenceSolver_ = solver_->clone();
    35793735  ourSolver_ = true ;
    35803736  cbcColLower_ = NULL;
     
    37573913  strongInfo_[1]=rhs.strongInfo_[1];
    37583914  strongInfo_[2]=rhs.strongInfo_[2];
     3915  solverCharacteristics_ = NULL;
    37593916  if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
    37603917  if (defaultHandler_) {
     
    40474204    strongInfo_[1]=rhs.strongInfo_[1];
    40484205    strongInfo_[2]=rhs.strongInfo_[2];
     4206    solverCharacteristics_ = NULL;
    40494207    lastHeuristic_ = NULL;
    40504208    numberCutGenerators_ = rhs.numberCutGenerators_;
     
    42194377  delete [] continuousSolution_;
    42204378  continuousSolution_=NULL;
     4379  solverCharacteristics_=NULL;
    42214380  delete [] usedInSolution_;
    42224381  usedInSolution_ = NULL;
     
    46434802    Update lastws
    46444803*/
     4804    int nxrow = lastws->getNumArtificial();
    46454805    for (i=0;i<currentNumberCuts;i++) {
     4806      assert (i+numberRowsAtContinuous_<nxrow);
    46464807      CoinWarmStartBasis::Status status =
    46474808        lastws->getArtifStatus(i+numberRowsAtContinuous_);
     
    47164877int CbcModel::reducedCostFix ()
    47174878
    4718 { double cutoff = getCutoff() ;
     4879{
     4880  if(!solverCharacteristics_->reducedCostsAccurate())
     4881    return 0; //NLP
     4882  double cutoff = getCutoff() ;
    47194883  double direction = solver_->getObjSense() ;
    47204884  double gap = cutoff - solver_->getObjValue()*direction ;
     
    48004964                       
    48014965
    4802 { bool feasible = true ;
     4966{
     4967  bool feasible = true ;
    48034968  int lastNumberCuts = 0 ;
    48044969  double lastObjective = -1.0e100 ;
     
    49895154*/
    49905155    // This should only happen after heuristic solution
    4991     if (!solver_->optimalBasisIsAvailable()) {
     5156    if (solverCharacteristics_->warmStart()&&
     5157        !solver_->optimalBasisIsAvailable()) {
    49925158      //printf("XXXXYY no opt basis\n");
    49935159      resolve();
     
    54525618    }
    54535619  } while (numberTries>0) ;
     5620  //check feasibility.
     5621  //If solution seems to be integer feasible calling setBestSolution
     5622  //will eventually add extra global cuts which we need to install at
     5623  //the nodes
     5624
     5625
     5626  if(feasible&&solverCharacteristics_->solutionAddsCuts())  //check integer feasibility
     5627  {
     5628    bool integerFeasible = true;
     5629    const double * save = testSolution_;
     5630    testSolution_ = solver_->getColSolution();
     5631    for (int i=0;i<numberObjects_ && integerFeasible;i++)
     5632    {
     5633      int preferredWay;
     5634      double infeasibility = object_[i]->infeasibility(preferredWay);
     5635      if(infeasibility)
     5636        integerFeasible = false;
     5637    }
     5638    testSolution_ = save;
     5639    if(integerFeasible)//update
     5640    {
     5641      double objValue = solver_->getObjValue();
     5642      //int numberGlobalBefore = globalCuts_.sizeRowCuts();
     5643      // SOLUTION2 so won't up cutoff or print message
     5644      setBestSolution(CBC_SOLUTION2, objValue,
     5645                      solver_->getColSolution(),0);
     5646#if 0
     5647      int numberGlobalAfter = globalCuts_.sizeRowCuts();
     5648      int numberToAdd = numberGlobalAfter - numberGlobalBefore;
     5649      if(numberToAdd > 0)
     5650        //We have added some cuts say they are tight at that node
     5651        //Basis and lp should already have been updated
     5652      {
     5653        feasible = (solver_->isProvenOptimal() &&
     5654                    !solver_->isDualObjectiveLimitReached()) ;
     5655        if(feasible)
     5656        {
     5657          int numberCuts = numberNewCuts_ = cuts.sizeRowCuts();
     5658      // possibly extend whichGenerator
     5659          resizeWhichGenerator(numberCuts, numberToAdd+numberCuts);
     5660         
     5661          for (int i = numberGlobalBefore ; i < numberGlobalAfter ; i++)
     5662            {       
     5663              whichGenerator_[numberNewCuts_++]=-1;
     5664              cuts.insert(globalCuts_.rowCut(i)) ;
     5665            }
     5666          numberNewCuts_ = lastNumberCuts+numberToAdd;
     5667          //now take off the cuts which are not tight anymore
     5668          takeOffCuts(cuts,resolveAfterTakeOffCuts_, NULL) ;
     5669          if (solver_->isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_)
     5670            { feasible = false ;}
     5671        }
     5672        if(!feasible) //node will be fathomed
     5673          {
     5674          for (int i = 0;i<currentNumberCuts_;i++)
     5675            {
     5676            // take off node
     5677            if (addedCuts_[i])
     5678              {
     5679                if (!addedCuts_[i]->decrement())
     5680                  delete addedCuts_[i] ;
     5681                addedCuts_[i] = NULL ;
     5682              }
     5683            }
     5684        }
     5685      }
     5686#endif
     5687    }
     5688  }
    54545689  // Reduced cost fix at end
    54555690  if (solver_->isProvenOptimal())
     
    55295764  class.
    55305765*/
     5766#ifdef NODE_LOG
     5767  int fatherNum = (node == NULL) ? -1 : node->nodeNumber();
     5768  double value =  (node == NULL) ? -1 : node->branchingObject()->value();
     5769  string bigOne = (solver_->getIterationCount() > 30) ? "*******" : "";
     5770  string way = (node == NULL) ? "" : (node->branchingObject()->way())==1 ? "Down" : "Up";
     5771  std::cout<<"Node "<<numberNodes_<<", father "<<fatherNum<<", #iterations "<<solver_->getIterationCount()<<", sol value : "<<solver_->getObjValue()<<std::endl;
     5772#endif
    55315773  if (fullScan&&numberCutGenerators_) {
    55325774    /* If cuts just at root node then it will probably be faster to
     
    55365778    // get sizes
    55375779    int numberRowsAdded = solver_->getNumRows()-numberRowsAtStart;
    5538     CoinBigIndex numberElementsAdded = solver_->getNumElements()-numberElementsAtStart;
     5780    CoinBigIndex numberElementsAdded =  solver_->getNumElements()-numberElementsAtStart ;
    55395781    double densityOld = ((double) numberElementsAtStart)/((double) numberRowsAtStart);
    55405782    double densityNew = numberRowsAdded ? ((double) (numberElementsAdded))/((double) numberRowsAdded)
     
    59466188    }
    59476189  }
     6190#if 0
     6191  {
     6192    int iColumn;
     6193    int numberColumns = solver_->getNumCols();
     6194    const double * columnLower = solver_->getColLower();
     6195    const double * columnUpper = solver_->getColUpper();
     6196    const double * sol = solver_->getColSolution();
     6197    bool feasible=true;
     6198    double xx[37];
     6199    memset(xx,0,sizeof(xx));
     6200    xx[6]=xx[7]=xx[9]=xx[18]=xx[26]=xx[36]=1.0;
     6201    for (iColumn= 0;iColumn<numberColumns;iColumn++) {
     6202      if (solver_->isInteger(iColumn)) {
     6203        //printf("yy %d at %g (bounds %g, %g)\n",iColumn,
     6204        //     sol[iColumn],columnLower[iColumn],columnUpper[iColumn]);
     6205        if (columnLower[iColumn]>xx[iColumn]||
     6206            columnUpper[iColumn]<xx[iColumn])
     6207          feasible=false;
     6208        //solver_->setColLower(iColumn,CoinMin(xx[iColumn],columnUpper[iColumn]));
     6209        //solver_->setColUpper(iColumn,CoinMax(xx[iColumn],columnLower[iColumn]));
     6210      }
     6211    }
     6212    if (feasible)
     6213      printf("On Feasible path\n");
     6214  }
     6215#endif
    59486216/*
    59496217  Reoptimize. Consider the possibility that we should fathom on bounds. But be
     
    66866954                         bool fixVariables)
    66876955
    6688 { int numberColumns = solver_->getNumCols();
    6689 
    6690   /*
    6691     Grab the continuous solver (the pristine copy of the problem, made before
    6692     starting to work on the root node). Save the bounds on the variables.
    6693     Install the solution passed as a parameter, and copy it to the model's
    6694     currentSolution_.
     6956{
     6957  if (!solverCharacteristics_->solutionAddsCuts()) {
     6958    // Can trust solution
     6959    int numberColumns = solver_->getNumCols();
    66956960   
    6696     TODO: This is a belt-and-suspenders approach. Once the code has settled
    6697           a bit, we can cast a critical eye here.
    6698   */
    6699   OsiSolverInterface * saveSolver = solver_;
    6700   if (continuousSolver_)
    6701     solver_ = continuousSolver_;
    6702   // move solution to continuous copy
    6703   solver_->setColSolution(solution);
    6704   // Put current solution in safe place
    6705   // Point to current solution
    6706   const double * save = testSolution_;
    6707   // Safe as will be const inside infeasibility()
    6708   testSolution_ = solver_->getColSolution();
    6709   //memcpy(currentSolution_,solver_->getColSolution(),
    6710   // numberColumns*sizeof(double));
    6711   //solver_->messageHandler()->setLogLevel(4);
    6712 
    6713   // save original bounds
    6714   double * saveUpper = new double[numberColumns];
    6715   double * saveLower = new double[numberColumns];
    6716   memcpy(saveUpper,getColUpper(),numberColumns*sizeof(double));
    6717   memcpy(saveLower,getColLower(),numberColumns*sizeof(double));
    6718 
    6719   /*
    6720     Run through the objects and use feasibleRegion() to set variable bounds
    6721     so as to fix the variables specified in the objects at their value in this
    6722     solution. Since the object list contains (at least) one object for every
    6723     integer variable, this has the effect of fixing all integer variables.
    6724   */
    6725   int i;
    6726   for (i=0;i<numberObjects_;i++)
    6727     object_[i]->feasibleRegion();
    6728   // We can switch off check
    6729   if ((specialOptions_&4)==0) {
    6730     if ((specialOptions_&2)==0) {
    6731       /*
    6732         Remove any existing warm start information to be sure there is no
    6733         residual influence on initialSolve().
    6734       */
    6735       CoinWarmStartBasis *slack =
    6736         dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
    6737       solver_->setWarmStart(slack);
    6738       delete slack ;
    6739     }
    6740     // Give a hint to do dual
    6741     bool saveTakeHint;
    6742     OsiHintStrength saveStrength;
    6743     bool gotHint = (solver_->getHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength));
    6744     assert (gotHint);
    6745     solver_->setHintParam(OsiDoDualInInitial,true,OsiHintTry);
    6746     solver_->initialSolve();
    6747     if (!solver_->isProvenOptimal())
    6748       { //printf("checkSolution infeas! Retrying wihout scaling.\n");
     6961    /*
     6962      Grab the continuous solver (the pristine copy of the problem, made before
     6963      starting to work on the root node). Save the bounds on the variables.
     6964      Install the solution passed as a parameter, and copy it to the model's
     6965      currentSolution_.
     6966     
     6967      TODO: This is a belt-and-suspenders approach. Once the code has settled
     6968      a bit, we can cast a critical eye here.
     6969    */
     6970    OsiSolverInterface * saveSolver = solver_;
     6971    if (continuousSolver_)
     6972      solver_ = continuousSolver_;
     6973    // move solution to continuous copy
     6974    solver_->setColSolution(solution);
     6975    // Put current solution in safe place
     6976    // Point to current solution
     6977    const double * save = testSolution_;
     6978    // Safe as will be const inside infeasibility()
     6979    testSolution_ = solver_->getColSolution();
     6980    //memcpy(currentSolution_,solver_->getColSolution(),
     6981    // numberColumns*sizeof(double));
     6982    //solver_->messageHandler()->setLogLevel(4);
     6983   
     6984    // save original bounds
     6985    double * saveUpper = new double[numberColumns];
     6986    double * saveLower = new double[numberColumns];
     6987    memcpy(saveUpper,getColUpper(),numberColumns*sizeof(double));
     6988    memcpy(saveLower,getColLower(),numberColumns*sizeof(double));
     6989   
     6990    /*
     6991      Run through the objects and use feasibleRegion() to set variable bounds
     6992      so as to fix the variables specified in the objects at their value in this
     6993      solution. Since the object list contains (at least) one object for every
     6994      integer variable, this has the effect of fixing all integer variables.
     6995    */
     6996    int i;
     6997    for (i=0;i<numberObjects_;i++)
     6998      object_[i]->feasibleRegion();
     6999    // We can switch off check
     7000    if ((specialOptions_&4)==0) {
     7001      if ((specialOptions_&2)==0) {
     7002        /*
     7003          Remove any existing warm start information to be sure there is no
     7004          residual influence on initialSolve().
     7005        */
     7006        CoinWarmStartBasis *slack =
     7007          dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
     7008        solver_->setWarmStart(slack);
     7009        delete slack ;
     7010      }
     7011      // Give a hint to do dual
    67497012      bool saveTakeHint;
    67507013      OsiHintStrength saveStrength;
    6751       bool savePrintHint;
    6752       //solver_->writeMps("infeas");
    6753       bool gotHint = (solver_->getHintParam(OsiDoReducePrint,savePrintHint,saveStrength));
    6754       gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
    6755       solver_->setHintParam(OsiDoScale,false,OsiHintTry);
    6756       //solver_->setHintParam(OsiDoReducePrint,false,OsiHintTry) ;
    6757       solver_->setHintParam(OsiDoDualInInitial,false,OsiHintTry);
     7014      bool gotHint = (solver_->getHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength));
     7015      assert (gotHint);
     7016      solver_->setHintParam(OsiDoDualInInitial,true,OsiHintTry);
    67587017      solver_->initialSolve();
    6759       solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
    6760       //solver_->setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ;
    6761       }
    6762     //assert(solver_->isProvenOptimal());
    6763     solver_->setHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength);
    6764   }
    6765   double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
    6766 
    6767   /*
    6768     Check that the solution still beats the objective cutoff.
    6769 
    6770     If it passes, make a copy of the primal variable values and do some
    6771     cleanup and checks:
    6772     + Values of all variables are are within original bounds and values of
     7018      if (!solver_->isProvenOptimal())
     7019        { //printf("checkSolution infeas! Retrying wihout scaling.\n");
     7020          bool saveTakeHint;
     7021          OsiHintStrength saveStrength;
     7022          bool savePrintHint;
     7023          //solver_->writeMps("infeas");
     7024          bool gotHint = (solver_->getHintParam(OsiDoReducePrint,savePrintHint,saveStrength));
     7025          gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
     7026          solver_->setHintParam(OsiDoScale,false,OsiHintTry);
     7027          //solver_->setHintParam(OsiDoReducePrint,false,OsiHintTry) ;
     7028          solver_->setHintParam(OsiDoDualInInitial,false,OsiHintTry);
     7029          solver_->initialSolve();
     7030          solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
     7031          //solver_->setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ;
     7032        }
     7033      //assert(solver_->isProvenOptimal());
     7034      solver_->setHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength);
     7035    }
     7036    double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
     7037   
     7038    /*
     7039      Check that the solution still beats the objective cutoff.
     7040     
     7041      If it passes, make a copy of the primal variable values and do some
     7042      cleanup and checks:
     7043      + Values of all variables are are within original bounds and values of
    67737044      all integer variables are within tolerance of integral.
    6774     + There are no constraint violations.
    6775     There really should be no need for the check against original bounds.
    6776     Perhaps an opportunity for a sanity check?
    6777   */
    6778   if ((solver_->isProvenOptimal()||(specialOptions_&4)!=0) && objectiveValue <= cutoff)
    6779   {
    6780     double * solution = new double[numberColumns];
    6781     memcpy(solution ,solver_->getColSolution(),numberColumns*sizeof(double)) ;
    6782 
    6783     const double * rowLower = solver_->getRowLower() ;
    6784     const double * rowUpper = solver_->getRowUpper() ;
    6785     int numberRows = solver_->getNumRows() ;
    6786     double *rowActivity = new double[numberRows] ;
    6787     memset(rowActivity,0,numberRows*sizeof(double)) ;
    6788 
     7045      + There are no constraint violations.
     7046      There really should be no need for the check against original bounds.
     7047      Perhaps an opportunity for a sanity check?
     7048    */
     7049    if ((solver_->isProvenOptimal()||(specialOptions_&4)!=0) && objectiveValue <= cutoff)
     7050      {
     7051        double * solution = new double[numberColumns];
     7052        memcpy(solution ,solver_->getColSolution(),numberColumns*sizeof(double)) ;
     7053       
     7054        const double * rowLower = solver_->getRowLower() ;
     7055        const double * rowUpper = solver_->getRowUpper() ;
     7056        int numberRows = solver_->getNumRows() ;
     7057        double *rowActivity = new double[numberRows] ;
     7058        memset(rowActivity,0,numberRows*sizeof(double)) ;
     7059       
    67897060#ifndef NDEBUG
    6790     double integerTolerance = getIntegerTolerance() ;
     7061        double integerTolerance = getIntegerTolerance() ;
    67917062#endif
    6792     int iColumn;
    6793     for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
    6794     { double value = solution[iColumn] ;
    6795       value = CoinMax(value, saveLower[iColumn]) ;
    6796       value = CoinMin(value, saveUpper[iColumn]) ;
    6797       if (solver_->isInteger(iColumn))
    6798         assert(fabs(value-solution[iColumn]) <= integerTolerance) ;
    6799       solution[iColumn] = value ; }
     7063        int iColumn;
     7064        for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
     7065          { double value = solution[iColumn] ;
     7066          value = CoinMax(value, saveLower[iColumn]) ;
     7067          value = CoinMin(value, saveUpper[iColumn]) ;
     7068          if (solver_->isInteger(iColumn))
     7069            assert(fabs(value-solution[iColumn]) <= integerTolerance) ;
     7070          solution[iColumn] = value ; }
     7071       
     7072        solver_->getMatrixByCol()->times(solution,rowActivity) ;
     7073        delete [] solution;
     7074        double primalTolerance ;
     7075        solver_->getDblParam(OsiPrimalTolerance,primalTolerance) ;
     7076        double largestInfeasibility =0.0;
     7077        for (i=0 ; i < numberRows ; i++) {
     7078          largestInfeasibility = CoinMax(largestInfeasibility,
     7079                                         rowLower[i]-rowActivity[i]);
     7080          largestInfeasibility = CoinMax(largestInfeasibility,
     7081                                         rowActivity[i]-rowUpper[i]);
     7082        }
     7083        if (largestInfeasibility>100.0*primalTolerance)
     7084          handler_->message(CBC_NOTFEAS3, messages_)
     7085            << largestInfeasibility << CoinMessageEol ;
     7086       
     7087        delete [] rowActivity ; }
     7088    else
     7089      { objectiveValue=1.0e50 ; }
     7090    /*
     7091      Regardless of what we think of the solution, we may need to restore the
     7092      original bounds of the continuous solver. Unfortunately, const'ness
     7093      prevents us from simply reversing the memcpy used to make these snapshots.
     7094    */
     7095    if (!fixVariables)
     7096      { for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++)
     7097        { solver_->setColLower(iColumn,saveLower[iColumn]) ;
     7098        solver_->setColUpper(iColumn,saveUpper[iColumn]) ; } }
     7099    delete [] saveLower;
     7100    delete [] saveUpper;
    68007101   
    6801     solver_->getMatrixByCol()->times(solution,rowActivity) ;
    6802     delete [] solution;
    6803     double primalTolerance ;
    6804     solver_->getDblParam(OsiPrimalTolerance,primalTolerance) ;
    6805     double largestInfeasibility =0.0;
    6806     for (i=0 ; i < numberRows ; i++) {
    6807       largestInfeasibility = CoinMax(largestInfeasibility,
    6808                                  rowLower[i]-rowActivity[i]);
    6809       largestInfeasibility = CoinMax(largestInfeasibility,
    6810                                  rowActivity[i]-rowUpper[i]);
    6811     }
    6812     if (largestInfeasibility>100.0*primalTolerance)
    6813       handler_->message(CBC_NOTFEAS3, messages_)
    6814         << largestInfeasibility << CoinMessageEol ;
    6815 
    6816     delete [] rowActivity ; }
    6817   else
    6818   { objectiveValue=1.0e50 ; }
    6819   /*
    6820     Regardless of what we think of the solution, we may need to restore the
    6821     original bounds of the continuous solver. Unfortunately, const'ness
    6822     prevents us from simply reversing the memcpy used to make these snapshots.
    6823   */
    6824   if (!fixVariables)
    6825   { for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++)
    6826     { solver_->setColLower(iColumn,saveLower[iColumn]) ;
    6827       solver_->setColUpper(iColumn,saveUpper[iColumn]) ; } }
    6828   delete [] saveLower;
    6829   delete [] saveUpper;
     7102    /*
     7103      Restore the usual solver.
     7104    */
     7105    solver_=saveSolver;
     7106    testSolution_ = save;
     7107    return objectiveValue;
     7108  } else {
     7109    // Outer approximation or similar
     7110    //If this is true then the solution comes from the nlp we don't need to resolve the same nlp with ipopt
     7111    //solverCharacteristics_->setSolver(solver_);
     7112    bool solutionComesFromNlp = solverCharacteristics_->bestObjectiveValue()<cutoff;
     7113    double objectiveValue;
     7114    int numberColumns = solver_->getNumCols();
     7115    double *saveLower = NULL;
     7116    double * saveUpper = NULL;
    68307117   
    6831   /*
    6832     Restore the usual solver.
    6833   */
    6834   solver_=saveSolver;
    6835   testSolution_ = save;
    6836   return objectiveValue;
     7118    if(! solutionComesFromNlp)//Otherwise solution already comes from ipopt and cuts are known
     7119      {
     7120        if(fixVariables)//Will temporarily fix all integer valued var
     7121          {
     7122            // save original bounds
     7123            saveUpper = new double[numberColumns];
     7124            saveLower = new double[numberColumns];
     7125            memcpy(saveUpper,solver_->getColUpper(),numberColumns*sizeof(double));
     7126            memcpy(saveLower,solver_->getColLower(),numberColumns*sizeof(double));
     7127            //in any case solution should be already loaded into solver_
     7128            /*
     7129              Run through the objects and use feasibleRegion() to set variable bounds
     7130              so as to fix the variables specified in the objects at their value in this
     7131              solution. Since the object list contains (at least) one object for every
     7132              integer variable, this has the effect of fixing all integer variables.
     7133            */
     7134            const double * save = testSolution_;
     7135            testSolution_ = solution;
     7136            for (int i=0;i<numberObjects_;i++)
     7137              object_[i]->feasibleRegion();
     7138            testSolution_ = save;
     7139            solver_->resolve();
     7140          }
     7141       
     7142        /*
     7143          Now step through the cut generators and see if any of them are flagged to
     7144          run when a new solution is discovered. Only global cuts are useful.
     7145          (The solution being evaluated may not correspond to the current location in the
     7146          search tree --- discovered by heuristic, for example.)
     7147        */
     7148        OsiCuts theseCuts;
     7149        int i;
     7150        int lastNumberCuts=0;
     7151        for (i=0;i<numberCutGenerators_;i++)
     7152          {
     7153            if (generator_[i]->atSolution())
     7154              {
     7155                generator_[i]->generateCuts(theseCuts,true,NULL);
     7156                int numberCuts = theseCuts.sizeRowCuts();
     7157                for (int j=lastNumberCuts;j<numberCuts;j++)
     7158                  {
     7159                    const OsiRowCut * thisCut = theseCuts.rowCutPtr(j);
     7160                    if (thisCut->globallyValid())
     7161                      {
     7162                        //           if ((specialOptions_&1)!=0)
     7163                        //           {
     7164                        //             /* As these are global cuts -
     7165                        //             a) Always get debugger object
     7166                        // b) Not fatal error to cutoff optimal (if we have just got optimal)
     7167                        // */
     7168                        //             const OsiRowCutDebugger *debugger = solver_->getRowCutDebuggerAlways() ;
     7169                        //             if (debugger)
     7170                        //             {
     7171                        //               if(debugger->invalidCut(*thisCut))
     7172                        //                 printf("ZZZZ Global cut - cuts off optimal solution!\n");
     7173                        //             }
     7174                        //           }
     7175                        // add to global list
     7176                        globalCuts_.insert(*thisCut);
     7177                        generator_[i]->incrementNumberCutsInTotal();
     7178                      }
     7179                  }
     7180              }
     7181          }
     7182        //   int numberCuts = theseCuts.sizeColCuts();
     7183        //   for (i=0;i<numberCuts;i++) {
     7184        //     const OsiColCut * thisCut = theseCuts.colCutPtr(i);
     7185        //     if (thisCut->globallyValid()) {
     7186        //       // add to global list
     7187        //       globalCuts_.insert(*thisCut);
     7188        //     }
     7189        //   }
     7190        //have to retrieve the solution and its value from the nlp
     7191      }
     7192    double newObjectiveValue=cutoff;
     7193    if(solverCharacteristics_->solution(newObjectiveValue,
     7194                                        const_cast<double *> (solution),
     7195                                        numberColumns))
     7196      {
     7197        objectiveValue = newObjectiveValue;
     7198      }
     7199    else
     7200      {
     7201        objectiveValue = 2e50;
     7202      }
     7203    if (!solutionComesFromNlp && fixVariables)
     7204      {
     7205        for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++)
     7206          {
     7207            solver_->setColLower(iColumn,saveLower[iColumn]) ;
     7208            solver_->setColUpper(iColumn,saveUpper[iColumn]) ; 
     7209          }
     7210        delete [] saveLower;
     7211        delete [] saveUpper;
     7212      }
     7213    return objectiveValue;
     7214  }
    68377215}
    68387216
     
    68607238                           bool fixVariables)
    68617239
    6862 { double cutoff = getCutoff() ;
    6863 
    6864 /*
    6865   Double check the solution to catch pretenders.
    6866 */
    6867   objectiveValue = checkSolution(cutoff,solution,fixVariables);
    6868   if (objectiveValue > cutoff)
    6869   { if (objectiveValue>1.0e30)
    6870       handler_->message(CBC_NOTFEAS1, messages_) << CoinMessageEol ;
    6871     else
    6872       handler_->message(CBC_NOTFEAS2, messages_)
    6873         << objectiveValue << cutoff << CoinMessageEol ; }
    6874 /*
    6875   We have a winner. Install it as the new incumbent.
    6876   Bump the objective cutoff value and solution counts. Give the user the
    6877   good news.
    6878 */
    6879   else
    6880   { bestObjective_ = objectiveValue;
    6881     int numberColumns = solver_->getNumCols();
    6882     if (!bestSolution_)
    6883       bestSolution_ = new double[numberColumns];
    6884     CoinCopyN(solution,numberColumns,bestSolution_);
    6885 
    6886     cutoff = bestObjective_-dblParam_[CbcCutoffIncrement];
    6887     // This is not correct - that way cutoff can go up if maximization
    6888     //double direction = solver_->getObjSense();
    6889     //setCutoff(cutoff*direction);
    6890     setCutoff(cutoff);
    6891 
    6892     if (how==CBC_ROUNDING)
    6893       numberHeuristicSolutions_++;
    6894     numberSolutions_++;
    6895     if (numberHeuristicSolutions_==numberSolutions_)
    6896       stateOfSearch_ = 1;
    6897     else
    6898       stateOfSearch_ = 2;
    6899 
    6900     handler_->message(how,messages_)
    6901       <<bestObjective_<<numberIterations_
    6902       <<numberNodes_
    6903       <<CoinMessageEol;
    6904 /*
    6905   Now step through the cut generators and see if any of them are flagged to
    6906   run when a new solution is discovered. Only global cuts are useful. (The
    6907   solution being evaluated may not correspond to the current location in the
    6908   search tree --- discovered by heuristic, for example.)
    6909 */
    6910     OsiCuts theseCuts;
    6911     int i;
    6912     int lastNumberCuts=0;
    6913     for (i=0;i<numberCutGenerators_;i++) {
    6914       bool generate = generator_[i]->atSolution();
    6915       // skip if not optimal and should be (maybe a cut generator has fixed variables)
    6916       if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
    6917         generate=false;
    6918       if (generate) {
    6919         generator_[i]->generateCuts(theseCuts,true,NULL);
    6920         int numberCuts = theseCuts.sizeRowCuts();
    6921         for (int j=lastNumberCuts;j<numberCuts;j++) {
    6922           const OsiRowCut * thisCut = theseCuts.rowCutPtr(j);
    6923           if (thisCut->globallyValid()) {
    6924             if ((specialOptions_&1)!=0) {
    6925               /* As these are global cuts -
    6926                  a) Always get debugger object
    6927                  b) Not fatal error to cutoff optimal (if we have just got optimal)
    6928               */
    6929               const OsiRowCutDebugger *debugger = solver_->getRowCutDebuggerAlways() ;
    6930               if (debugger) {
    6931                 if(debugger->invalidCut(*thisCut))
    6932                   printf("ZZZZ Global cut - cuts off optimal solution!\n");
     7240{
     7241  if (!solverCharacteristics_->solutionAddsCuts()) {
     7242    // Can trust solution
     7243    double cutoff = CoinMin(getCutoff(),bestObjective_) ;
     7244   
     7245    /*
     7246      Double check the solution to catch pretenders.
     7247    */
     7248    objectiveValue = checkSolution(cutoff,solution,fixVariables);
     7249    if (objectiveValue > cutoff) {
     7250      if (objectiveValue>1.0e30)
     7251        handler_->message(CBC_NOTFEAS1, messages_) << CoinMessageEol ;
     7252      else
     7253        handler_->message(CBC_NOTFEAS2, messages_)
     7254          << objectiveValue << cutoff << CoinMessageEol ;
     7255    } else {
     7256      /*
     7257        We have a winner. Install it as the new incumbent.
     7258        Bump the objective cutoff value and solution counts. Give the user the
     7259        good news.
     7260      */
     7261      bestObjective_ = objectiveValue;
     7262      int numberColumns = solver_->getNumCols();
     7263      if (!bestSolution_)
     7264        bestSolution_ = new double[numberColumns];
     7265      CoinCopyN(solution,numberColumns,bestSolution_);
     7266     
     7267      cutoff = bestObjective_-dblParam_[CbcCutoffIncrement];
     7268      // This is not correct - that way cutoff can go up if maximization
     7269      //double direction = solver_->getObjSense();
     7270      //setCutoff(cutoff*direction);
     7271      setCutoff(cutoff);
     7272     
     7273      if (how==CBC_ROUNDING)
     7274        numberHeuristicSolutions_++;
     7275      numberSolutions_++;
     7276      if (numberHeuristicSolutions_==numberSolutions_)
     7277        stateOfSearch_ = 1;
     7278      else
     7279        stateOfSearch_ = 2;
     7280     
     7281      handler_->message(how,messages_)
     7282        <<bestObjective_<<numberIterations_
     7283        <<numberNodes_
     7284        <<CoinMessageEol;
     7285      /*
     7286        Now step through the cut generators and see if any of them are flagged to
     7287        run when a new solution is discovered. Only global cuts are useful. (The
     7288        solution being evaluated may not correspond to the current location in the
     7289        search tree --- discovered by heuristic, for example.)
     7290      */
     7291      OsiCuts theseCuts;
     7292      int i;
     7293      int lastNumberCuts=0;
     7294      for (i=0;i<numberCutGenerators_;i++) {
     7295        bool generate = generator_[i]->atSolution();
     7296        // skip if not optimal and should be (maybe a cut generator has fixed variables)
     7297        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
     7298          generate=false;
     7299        if (generate) {
     7300          generator_[i]->generateCuts(theseCuts,true,NULL);
     7301          int numberCuts = theseCuts.sizeRowCuts();
     7302          for (int j=lastNumberCuts;j<numberCuts;j++) {
     7303            const OsiRowCut * thisCut = theseCuts.rowCutPtr(j);
     7304            if (thisCut->globallyValid()) {
     7305              if ((specialOptions_&1)!=0) {
     7306                /* As these are global cuts -
     7307                   a) Always get debugger object
     7308                   b) Not fatal error to cutoff optimal (if we have just got optimal)
     7309                */
     7310                const OsiRowCutDebugger *debugger = solver_->getRowCutDebuggerAlways() ;
     7311                if (debugger) {
     7312                  if(debugger->invalidCut(*thisCut))
     7313                    printf("ZZZZ Global cut - cuts off optimal solution!\n");
     7314                }
    69337315              }
     7316              // add to global list
     7317              globalCuts_.insert(*thisCut);
     7318              generator_[i]->incrementNumberCutsInTotal();
    69347319            }
    6935             // add to global list
    6936             globalCuts_.insert(*thisCut);
    6937             generator_[i]->incrementNumberCutsInTotal();
    6938           }
    6939         }
    6940       }
    6941     }
    6942     int numberCuts = theseCuts.sizeColCuts();
    6943     for (i=0;i<numberCuts;i++) {
    6944       const OsiColCut * thisCut = theseCuts.colCutPtr(i);
    6945       if (thisCut->globallyValid()) {
    6946         // add to global list
    6947         globalCuts_.insert(*thisCut);
    6948       }
    6949     }
    6950   }
    6951 
    6952   return ; }
     7320          }
     7321        }
     7322      }
     7323      int numberCuts = theseCuts.sizeColCuts();
     7324      for (i=0;i<numberCuts;i++) {
     7325        const OsiColCut * thisCut = theseCuts.colCutPtr(i);
     7326        if (thisCut->globallyValid()) {
     7327          // add to global list
     7328          globalCuts_.insert(*thisCut);
     7329        }
     7330      }
     7331    }
     7332  } else {
     7333    // Outer approximation or similar
     7334    double cutoff = getCutoff() ;
     7335   
     7336    /*
     7337      Double check the solution to catch pretenders.
     7338    */
     7339   
     7340    int numberRowBefore = solver_->getNumRows();
     7341    int numberColBefore = solver_->getNumCols();
     7342    double *saveColSol;
     7343   
     7344    CoinWarmStart * saveWs;
     7345    // if(how!=CBC_SOLUTION) return;
     7346    if(how==CBC_ROUNDING||true)//We don't want to make any change to solver_
     7347      //take a snapshot of current state
     7348      {     
     7349        //save solution
     7350        saveColSol = new double[numberColBefore];
     7351        CoinCopyN(solver_->getColSolution(), numberColBefore, saveColSol);
     7352        //save warm start
     7353        saveWs = solver_->getWarmStart();     
     7354      }
     7355   
     7356    //run check solution this will eventually generate cuts
     7357    //if in strongBranching or heuristic will do only one cut generation iteration
     7358    // by fixing variables.
     7359    fixVariables = fixVariables || (how==CBC_ROUNDING) || (how==CBC_STRONGSOL);
     7360    double * candidate = new double[numberColBefore];
     7361    CoinCopyN(solution, numberColBefore, candidate);
     7362    objectiveValue = checkSolution(cutoff,candidate,fixVariables);
     7363   
     7364    //If it was an heuristic solution we have to clean up the solver
     7365    if (how==CBC_ROUNDING||true)
     7366      {
     7367        //delete the cuts
     7368        int currentNumberRowCuts = solver_->getNumRows() - numberRowBefore;
     7369        int currentNumberColCuts = solver_->getNumCols() - numberColBefore;
     7370        if(CoinMax(currentNumberColCuts, currentNumberRowCuts)>0)
     7371          {
     7372            int *which = new int[CoinMax(currentNumberColCuts, currentNumberRowCuts)];
     7373            if(currentNumberRowCuts)
     7374              {
     7375                for (int i = 0 ; i < currentNumberRowCuts ; i++)
     7376                  which[i] = i + numberRowBefore;
     7377               
     7378                solver_->deleteRows(currentNumberRowCuts,which);
     7379              }
     7380            if(currentNumberColCuts)
     7381              {
     7382                for (int i = 0 ; i < currentNumberColCuts ; i++)
     7383                  which[i] = i+numberColBefore;
     7384                solver_->deleteCols(currentNumberColCuts,which);
     7385              }
     7386            delete [] which;
     7387          }
     7388        // Reset solution and warm start info
     7389        solver_->setColSolution(saveColSol);
     7390        solver_->setWarmStart(saveWs);
     7391        delete [] saveColSol;
     7392        delete saveWs;
     7393      }
     7394   
     7395    if (objectiveValue > cutoff) {
     7396      if (!solverCharacteristics_->solutionAddsCuts()) {
     7397        if (objectiveValue>1.0e30)
     7398          handler_->message(CBC_NOTFEAS1, messages_) << CoinMessageEol ;
     7399        else
     7400          handler_->message(CBC_NOTFEAS2, messages_)
     7401            << objectiveValue << cutoff << CoinMessageEol ;
     7402      }
     7403    } else {
     7404      /*
     7405        We have a winner. Install it as the new incumbent.
     7406        Bump the objective cutoff value and solution counts. Give the user the
     7407        good news.
     7408        NB - Not all of this if from solve with cuts
     7409      */
     7410      bestObjective_ = objectiveValue;
     7411      int numberColumns = solver_->getNumCols();
     7412      if (!bestSolution_)
     7413        bestSolution_ = new double[numberColumns];
     7414      CoinCopyN(candidate,numberColumns,bestSolution_);
     7415
     7416      // don't update if from solveWithCuts
     7417      if (how!=CBC_SOLUTION2) {
     7418        if (how==CBC_ROUNDING)
     7419          numberHeuristicSolutions_++;
     7420        cutoff = bestObjective_-dblParam_[CbcCutoffIncrement];
     7421        // This is not correct - that way cutoff can go up if maximization
     7422        //double direction = solver_->getObjSense();
     7423        //setCutoff(cutoff*direction);
     7424        setCutoff(cutoff);
     7425     
     7426        numberSolutions_++;
     7427        if (numberNodes_== 0 || numberHeuristicSolutions_==numberSolutions_)
     7428          stateOfSearch_ = 1;
     7429        else
     7430          stateOfSearch_ = 2;
     7431       
     7432        handler_->message(how,messages_)
     7433          <<bestObjective_<<numberIterations_
     7434          <<numberNodes_
     7435          <<CoinMessageEol;
     7436      }
     7437    }
     7438    delete [] candidate;
     7439  }
     7440  return ;
     7441}
    69537442
    69547443
     
    83958884  cbcRowPrice_ = solver_->getRowPrice();
    83968885  /// Get a pointer to array[getNumCols()] (for speed) of reduced costs
    8397   cbcReducedCost_ = solver_->getReducedCost();
     8886  if(solverCharacteristics_->reducedCostsAccurate())
     8887    cbcReducedCost_ = solver_->getReducedCost();
     8888  else
     8889    cbcReducedCost_ = NULL;
    83988890  /// Pointer to array[getNumRows()] (for speed) of row activity levels.
    83998891  cbcRowActivity_ = solver_->getRowActivity();
  • trunk/CbcNode.cpp

    r262 r264  
    1313#define CUTS
    1414#include "OsiSolverInterface.hpp"
     15#include "OsiAuxInfo.hpp"
    1516#include "OsiSolverBranch.hpp"
    1617#include "CoinWarmStartBasis.hpp"
     
    491492  int numberColumns = model->getNumCols();
    492493  // move basis - but make sure size stays
    493   int numberRows = basis->getNumArtificial();
     494  // for bon-min - should not be needed int numberRows = model->getNumRows();
     495  int numberRows=basis->getNumArtificial();
    494496  delete basis ;
    495   basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
    496   basis->resize(numberRows,numberColumns);
     497  if (basis_) {
     498    basis = dynamic_cast<CoinWarmStartBasis *>(basis_->clone()) ;
     499    basis->resize(numberRows,numberColumns);
     500  } else {
     501    // We have a solver without a basis
     502    basis=NULL;
     503  }
    497504  for (i=0;i<numberCuts_;i++)
    498505    addCuts[currentNumberCuts+i]= cuts_[i];
     
    734741    //int maxBasisLength = ((iFull+15)>>4)+((numberColumns+15)>>4);
    735742    //printf("l %d full %d\n",maxBasisLength,iFull);
    736     expanded->resize(iFull,numberColumns);
     743    if (expanded)
     744      expanded->resize(iFull,numberColumns);
    737745#ifdef FULL_DEBUG
    738746    printf("Before expansion: orig %d, old %d, new %d, current %d\n",
     
    14021410                                     newObjectiveValue,
    14031411                                     solver->getColSolution()) ;
     1412              // only needed for odd solvers
     1413              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
     1414              objectiveChange = newObjectiveValue-objectiveValue_ ;
    14041415              model->setLastHeuristic(NULL);
    14051416              model->incrementUsed(solver->getColSolution());
     
    15031514                                     newObjectiveValue,
    15041515                                     solver->getColSolution()) ;
     1516              // only needed for odd solvers
     1517              newObjectiveValue = solver->getObjSense()*solver->getObjValue();
     1518              objectiveChange = newObjectiveValue-objectiveValue_ ;
    15051519              model->setLastHeuristic(NULL);
    15061520              model->incrementUsed(solver->getColSolution());
     
    18341848  branch_=NULL;
    18351849  OsiSolverInterface * solver = model->solver();
     1850  // get information on solver type
     1851  const OsiAuxInfo * auxInfo = solver->getAuxiliaryInfo();
     1852  const OsiBabSolver * auxiliaryInfo = dynamic_cast<const OsiBabSolver *> (auxInfo);
    18361853  //assert(objectiveValue_ == solver->getObjSense()*solver->getObjValue());
    18371854  double cutoff =model->getCutoff();
     
    18481865  int saveStateOfSearch = model->stateOfSearch();
    18491866  int numberStrong=model->numberStrong();
     1867  if (!auxiliaryInfo->warmStart())
     1868    numberStrong=0;
    18501869  // But make more likely to get out after some times
    18511870  int changeStrategy=numberStrong;
     
    23172336          solver->markHotStart();
    23182337        }
     2338        assert (auxiliaryInfo->warmStart());
    23192339        doneHotStart=true;
    23202340        xMark++;
     
    23662386          solver->markHotStart();
    23672387          doneHotStart=true;
     2388          assert (auxiliaryInfo->warmStart());
    23682389          xMark++;
    23692390          //if (solver->isProvenPrimalInfeasible())
     
    23782399        // Mark hot start
    23792400        doneHotStart=true;
     2401        assert (auxiliaryInfo->warmStart());
    23802402        solver->markHotStart();
    23812403        xMark++;
     
    25902612          choice.possibleBranch->branch() ;
    25912613          solver->solveFromHotStart() ;
     2614          bool needHotStartUpdate=false;
    25922615          numberStrongDone++;
    25932616          numberStrongIterations += solver->getIterationCount();
     
    26202643                                          choice.numObjInfeasDown)
    26212644                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
     2645                if (auxiliaryInfo->solutionAddsCuts()) {
     2646                  needHotStartUpdate=true;
     2647                  solver->unmarkHotStart();
     2648                }
    26222649                model->setBestSolution(CBC_STRONGSOL,
    26232650                                       newObjectiveValue,
    26242651                                       solver->getColSolution()) ;
     2652                if (needHotStartUpdate) {
     2653                  solver->resolve();
     2654                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
     2655                  objectiveChange = newObjectiveValue  - objectiveValue_;
     2656                  model->feasibleSolution(choice.numIntInfeasDown,
     2657                                          choice.numObjInfeasDown);
     2658                }
    26252659                model->setLastHeuristic(NULL);
    26262660                model->incrementUsed(solver->getColSolution());
     
    26492683              solver->setColUpper(j,saveUpper[j]);
    26502684          }
     2685          if(needHotStartUpdate) {
     2686            needHotStartUpdate = false;
     2687            solver->resolve();
     2688            //we may again have an integer feasible solution
     2689            int numberIntegerInfeasibilities;
     2690            int numberObjectInfeasibilities;
     2691            if (model->feasibleSolution(
     2692                                        numberIntegerInfeasibilities,
     2693                                        numberObjectInfeasibilities)) {
     2694              double objValue = solver->getObjValue();
     2695              model->setBestSolution(CBC_STRONGSOL,
     2696                                     objValue,
     2697                                     solver->getColSolution()) ;
     2698              solver->resolve();
     2699              cutoff =model->getCutoff();
     2700            }
     2701            solver->markHotStart();
     2702          }
     2703          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
    26512704          //printf("Down on %d, status is %d, obj %g its %d cost %g finished %d inf %d infobj %d\n",
    26522705          //     choice.objectNumber,iStatus,newObjectiveValue,choice.numItersDown,
     
    26882741                                          choice.numObjInfeasUp)
    26892742                  &&model->problemFeasibility()->feasible(model,-1)>=0) {
     2743                if (auxiliaryInfo->solutionAddsCuts()) {
     2744                  needHotStartUpdate=true;
     2745                  solver->unmarkHotStart();
     2746                }
    26902747                model->setBestSolution(CBC_STRONGSOL,
    26912748                                       newObjectiveValue,
    26922749                                       solver->getColSolution()) ;
     2750                if (needHotStartUpdate) {
     2751                  solver->resolve();
     2752                  newObjectiveValue = solver->getObjSense()*solver->getObjValue();
     2753                  objectiveChange = newObjectiveValue  - objectiveValue_;
     2754                  model->feasibleSolution(choice.numIntInfeasDown,
     2755                                          choice.numObjInfeasDown);
     2756                }
    26932757                model->setLastHeuristic(NULL);
    26942758                model->incrementUsed(solver->getColSolution());
     
    27162780            if (saveUpper[j] != upper[j])
    27172781              solver->setColUpper(j,saveUpper[j]);
     2782          }
     2783          if(needHotStartUpdate) {
     2784            needHotStartUpdate = false;
     2785            solver->resolve();
     2786            //we may again have an integer feasible solution
     2787            int numberIntegerInfeasibilities;
     2788            int numberObjectInfeasibilities;
     2789            if (model->feasibleSolution(
     2790                                        numberIntegerInfeasibilities,
     2791                                        numberObjectInfeasibilities)) {
     2792              double objValue = solver->getObjValue();
     2793              model->setBestSolution(CBC_STRONGSOL,
     2794                                     objValue,
     2795                                     solver->getColSolution()) ;
     2796              solver->resolve();
     2797              cutoff =model->getCutoff();
     2798            }
     2799            solver->markHotStart();
    27182800          }
    27192801         
  • trunk/include/CbcHeuristic.hpp

    r197 r264  
    144144};
    145145
     146/** heuristic - just picks up any good solution
     147    found by solver - see OsiBabSolver
     148 */
     149
     150class CbcSerendipity : public CbcHeuristic {
     151public:
     152
     153  // Default Constructor
     154  CbcSerendipity ();
     155
     156  /* Constructor with model
     157  */
     158  CbcSerendipity (CbcModel & model);
     159 
     160  // Copy constructor
     161  CbcSerendipity ( const CbcSerendipity &);
     162   
     163  // Destructor
     164  ~CbcSerendipity ();
     165 
     166  /// Clone
     167  virtual CbcHeuristic * clone() const;
     168
     169  /// update model
     170  virtual void setModel(CbcModel * model);
     171 
     172  /** returns 0 if no solution, 1 if valid solution.
     173      Sets solution values if good, sets objective value (only if good)
     174      We leave all variables which are at one at this node of the
     175      tree to that value and will
     176      initially set all others to zero.  We then sort all variables in order of their cost
     177      divided by the number of entries in rows which are not yet covered.  We randomize that
     178      value a bit so that ties will be broken in different ways on different runs of the heuristic.
     179      We then choose the best one and set it to one and repeat the exercise. 
     180
     181  */
     182  virtual int solution(double & objectiveValue,
     183                       double * newSolution);
     184  /// Resets stuff if model changes
     185  virtual void resetModel(CbcModel * model);
     186
     187protected:
     188private:
     189  /// Illegal Assignment operator
     190  CbcSerendipity & operator=(const CbcSerendipity& rhs);
     191};
    146192
    147193#endif
  • trunk/include/CbcMessage.hpp

    r262 r264  
    2727  CBC_EVENT,
    2828  CBC_SOLUTION,
     29  CBC_SOLUTION2,
    2930  CBC_END,
    3031  CBC_INFEAS,
  • trunk/include/CbcModel.hpp

    r262 r264  
    2020class CbcCutGenerator;
    2121class OsiRowCut;
     22class OsiBabSolver;
    2223class OsiRowCutDebugger;
    2324class CglCutGenerator;
     
    12181219    /// Get application data
    12191220    void * getApplicationData() const;
     1221  /**
     1222      For advanced applications you may wish to modify the behavior of Cbc
     1223      e.g. if the solver is a NLP solver then you may not have an exact
     1224      optimum solution at each step.  Information could be built into
     1225      OsiSolverInterface but this is an alternative so that that interface
     1226      does not have to be changed.  If something similar is useful to
     1227      enough solvers then it could be migrated
     1228  */
     1229  void passInSolverCharacteristics(OsiBabSolver * solverCharacteristics);
    12201230  //@}
    12211231 
     
    17621772  /** 0 - number times strong branching done, 1 - number fixed, 2 - number infeasible */
    17631773  int strongInfo_[3];
     1774  /**
     1775      For advanced applications you may wish to modify the behavior of Cbc
     1776      e.g. if the solver is a NLP solver then you may not have an exact
     1777      optimum solution at each step.  This gives characteristics - just for one BAB.
     1778      For actually saving/restoring a solution you need the actual solver one.
     1779  */
     1780  OsiBabSolver * solverCharacteristics_;
    17641781  /// Whether to force a resolve after takeOffCuts
    17651782  bool resolveAfterTakeOffCuts_;
Note: See TracChangeset for help on using the changeset viewer.