Changeset 1240


Ignore:
Timestamp:
Oct 2, 2009 2:41:44 PM (10 years ago)
Author:
forrest
Message:

adding stuff for heuristics

Location:
trunk/Cbc/src
Files:
15 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/CbcBranchDynamic.cpp

    r1232 r1240  
    541541//#define FUNNY_BRANCHING 
    542542double
    543 CbcSimpleIntegerDynamicPseudoCost::infeasibility(const OsiBranchingInformation * /*info*/,
     543CbcSimpleIntegerDynamicPseudoCost::infeasibility(const OsiBranchingInformation * info,
    544544                               int &preferredWay) const
    545545{
     
    699699  if (preferredWay_)
    700700    preferredWay=preferredWay_;
     701  if (info->hotstartSolution_) {
     702    double targetValue = info->hotstartSolution_[columnNumber_];
     703    if (value>targetValue)
     704      preferredWay=-1;
     705    else
     706      preferredWay=1;
     707  }
    701708  // weight at 1.0 is max min
    702709#define WEIGHT_AFTER 0.8
  • trunk/Cbc/src/CbcCompareActual.cpp

    r1219 r1240  
    311311    double testX = x->guessedObjectiveValue();
    312312    double testY = y->guessedObjectiveValue();
    313 #elif THY_THIS==3
     313#elif TRY_THIS==3
    314314#define THRESH 0.95
     315    // Use estimates
     316    double testX = x->guessedObjectiveValue();
     317    double testY = y->guessedObjectiveValue();
    315318    if (x->objectiveValue()-bestPossible_>THRESH*(cutoff_-bestPossible_))
    316319      testX *= 2.0; // make worse
  • trunk/Cbc/src/CbcHeuristicDive.cpp

    r1224 r1240  
    179179  return shouldHeurRun_randomChoice();
    180180}
    181 
    182 struct PseudoReducedCost {
    183   int var;
    184   double pseudoRedCost;
    185 };
    186181
    187182inline bool compareBinaryVars(const PseudoReducedCost obj1,
     
    453448#endif
    454449    if (reducedCost&&true) {
     450#if 1
     451      cnt=fixOtherVariables(solver,solution,candidate,random);
     452#else
    455453#ifdef GAP
    456454      double cutoff = model_->getCutoff() ;
     
    501499      printf("cutoff %g obj %g - %d free, %d fixed\n",
    502500             model_->getCutoff(),solver->getObjValue(),numberFree,numberFixed);
     501#endif
    503502#endif
    504503#endif
     
    1000999  return numberFixed;
    10011000}
     1001// Fix other variables at bounds
     1002int
     1003CbcHeuristicDive::fixOtherVariables(OsiSolverInterface * solver,
     1004                                    const double * solution,
     1005                                    PseudoReducedCost * candidate,
     1006                                    const double * random)
     1007{
     1008  const double * lower = solver->getColLower();
     1009  const double * upper = solver->getColUpper();
     1010  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1011  double primalTolerance;
     1012  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
     1013
     1014  int numberIntegers = model_->numberIntegers();
     1015  const int * integerVariable = model_->integerVariable();
     1016  const double* reducedCost = solver->getReducedCost();
     1017  // fix other integer variables that are at their bounds
     1018  int cnt=0;
     1019#ifdef GAP
     1020  double direction = solver->getObjSense(); // 1 for min, -1 for max
     1021  double gap=1.0e30;
     1022#endif
     1023#ifdef GAP
     1024  double cutoff = model_->getCutoff() ;
     1025  if (cutoff<1.0e20&&false) {
     1026    double direction = solver->getObjSense() ;
     1027    gap = cutoff - solver->getObjValue()*direction ;
     1028    gap *= 0.1; // Fix more if plausible
     1029    double tolerance;
     1030    solver->getDblParam(OsiDualTolerance,tolerance) ;
     1031    if (gap<=0.0)
     1032      gap = tolerance;
     1033    gap += 100.0*tolerance;
     1034  }
     1035  int nOverGap=0;
     1036#endif
     1037  int numberFree=0;
     1038  int numberFixedAlready=0;
     1039  for (int i=0; i<numberIntegers; i++) {
     1040    int iColumn = integerVariable[i];
     1041    if (upper[iColumn]>lower[iColumn]) {
     1042      numberFree++;
     1043      double value=solution[iColumn];
     1044      if(fabs(floor(value+0.5)-value)<=integerTolerance) {
     1045        candidate[cnt].var = iColumn;
     1046        candidate[cnt++].pseudoRedCost =
     1047          fabs(reducedCost[iColumn]*random[i]);
     1048#ifdef GAP
     1049        if (fabs(reducedCost[iColumn])>gap)
     1050          nOverGap++;
     1051#endif
     1052      }
     1053    } else {
     1054      numberFixedAlready++;
     1055    }
     1056  }
     1057#ifdef GAP
     1058  int nLeft = maxNumberToFix-numberFixedAlready;
     1059#ifdef CLP_INVESTIGATE
     1060  printf("cutoff %g obj %g nover %d - %d free, %d fixed\n",
     1061         cutoff,solver->getObjValue(),nOverGap,numberFree,
     1062         numberFixedAlready);
     1063#endif
     1064  if (nOverGap>nLeft&&true) {
     1065    nOverGap = CoinMin(nOverGap,nLeft+maxNumberToFix/2);
     1066    maxNumberToFix += nOverGap-nLeft;
     1067  }
     1068#else
     1069#ifdef CLP_INVESTIGATE
     1070  printf("cutoff %g obj %g - %d free, %d fixed\n",
     1071         model_->getCutoff(),solver->getObjValue(),numberFree,
     1072         numberFixedAlready);
     1073#endif
     1074#endif
     1075  return cnt;
     1076}
  • trunk/Cbc/src/CbcHeuristicDive.hpp

    r1224 r1240  
    66
    77#include "CbcHeuristic.hpp"
     8struct PseudoReducedCost {
     9  int var;
     10  double pseudoRedCost;
     11};
     12
    813
    914/** Dive class
     
    98103  /// Perform reduced cost fixing on integer variables
    99104  int reducedCostFix (OsiSolverInterface* solver);
     105  /// Fix other variables at bounds
     106  virtual int fixOtherVariables(OsiSolverInterface * solver,
     107                                const double * solution,
     108                                PseudoReducedCost * candidate,
     109                                const double * random);
    100110
    101111protected:
  • trunk/Cbc/src/CbcHeuristicDivePseudoCost.cpp

    r1230 r1240  
    165165      if (obj1) {
    166166        //int iColumn = obj1->columnNumber();
    167         double downPseudoCost = obj1->downDynamicPseudoCost();
     167        double downPseudoCost = 1.0e-2*obj1->downDynamicPseudoCost();
    168168        double downShadow = obj1->downShadowPrice();
    169         double upPseudoCost = obj1->upDynamicPseudoCost();
     169        double upPseudoCost = 1.0e-2*obj1->upDynamicPseudoCost();
    170170        double upShadow = obj1->upShadowPrice();
    171171        downPseudoCost = CoinMax(downPseudoCost,downShadow);
     
    179179  }
    180180}
     181// Fix other variables at bounds
     182int
     183CbcHeuristicDivePseudoCost::fixOtherVariables(OsiSolverInterface * solver,
     184                                              const double * solution,
     185                                              PseudoReducedCost * candidate,
     186                                              const double * random)
     187{
     188  const double * lower = solver->getColLower();
     189  const double * upper = solver->getColUpper();
     190  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
     191  double primalTolerance;
     192  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
     193
     194  int numberIntegers = model_->numberIntegers();
     195  const int * integerVariable = model_->integerVariable();
     196  const double* reducedCost = solver->getReducedCost();
     197  // fix other integer variables that are at their bounds
     198  int cnt=0;
     199  int numberFree=0;
     200  int numberFixedAlready=0;
     201  for (int i=0; i<numberIntegers; i++) {
     202    int iColumn = integerVariable[i];
     203    if (upper[iColumn]>lower[iColumn]) {
     204      numberFree++;
     205      double value=solution[iColumn];
     206      if(value-lower[iColumn]<=integerTolerance) {
     207        candidate[cnt].var = iColumn;
     208        candidate[cnt++].pseudoRedCost = CoinMax(1.0e-2*reducedCost[iColumn],
     209                                                 downArray_[i])*random[i];
     210      } else if(upper[iColumn]-value<=integerTolerance) {
     211        candidate[cnt].var = iColumn;
     212        candidate[cnt++].pseudoRedCost = CoinMax(-1.0e-2*reducedCost[iColumn],
     213                                                 downArray_[i])*random[i];
     214      }
     215    } else {
     216      numberFixedAlready++;
     217    }
     218  }
     219#ifdef CLP_INVESTIGATE
     220  printf("cutoff %g obj %g - %d free, %d fixed\n",
     221         model_->getCutoff(),solver->getObjValue(),numberFree,
     222         numberFixedAlready);
     223#endif
     224  return cnt;
     225  //return CbcHeuristicDive::fixOtherVariables(solver, solution,
     226  //                                 candidate, random);
     227}
  • trunk/Cbc/src/CbcHeuristicDivePseudoCost.hpp

    r1224 r1240  
    4747      in selectVariableToBranch */
    4848  virtual void initializeData() ;
    49 
     49  /// Fix other variables at bounds
     50  virtual int fixOtherVariables(OsiSolverInterface * solver,
     51                                const double * solution,
     52                                PseudoReducedCost * candidate,
     53                                const double * random);
    5054
    5155};
  • trunk/Cbc/src/CbcHeuristicFPump.cpp

    r1232 r1240  
    1919#include "CbcHeuristicFPump.hpp"
    2020#include "CbcBranchActual.hpp"
     21#include "CbcBranchDynamic.hpp"
    2122#include "CoinHelperFunctions.hpp"
    2223#include "CoinWarmStartBasis.hpp"
     
    328329  int * closestSolution = general ? NULL : new int[numberIntegers];
    329330  double closestObjectiveValue = COIN_DBL_MAX;
    330  
     331
    331332  int numberIntegersOrig = numberIntegers;
    332333  numberIntegers = j;
     
    574575    double primalTolerance;
    575576    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
    576    
    577577   
    578578    // 2 space for last rounded solutions
     
    13251325#endif
    13261326              //roundingObjective = newSolutionValue;
    1327             } else {
    1328               roundingObjective = COIN_DBL_MAX;
     1327              //} else {
     1328              //roundingObjective = COIN_DBL_MAX;
    13291329            }
    13301330            model_->swapSolver(saveSolver);
     
    19401940      if (absoluteIncrement_>0.0||relativeIncrement_>0.0) {
    19411941        double gap = relativeIncrement_*fabs(solutionValue);
    1942         cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
     1942        double change = CoinMax(gap,absoluteIncrement_);
     1943        cutoff = CoinMin(cutoff,solutionValue-change);
    19431944      } else {
    19441945        //double weights[10]={0.1,0.1,0.2,0.2,0.2,0.3,0.3,0.3,0.4,0.5};
  • trunk/Cbc/src/CbcHeuristicRINS.cpp

    r1221 r1240  
    777777}
    778778
     779// Default Constructor
     780CbcHeuristicVND::CbcHeuristicVND()
     781  :CbcHeuristic()
     782{
     783  numberSolutions_=0;
     784  numberSuccesses_=0;
     785  numberTries_=0;
     786  lastNode_=-999999;
     787  howOften_=100;
     788  decayFactor_ = 0.5;
     789  baseSolution_=NULL;
     790  whereFrom_=1+8+255*256;
     791  stepSize_=0;
     792  k_=0;
     793  kmax_=0;
     794  nDifferent_=0;
     795}
     796
     797// Constructor with model - assumed before cuts
     798
     799CbcHeuristicVND::CbcHeuristicVND(CbcModel & model)
     800  :CbcHeuristic(model)
     801{
     802  numberSolutions_=0;
     803  numberSuccesses_=0;
     804  numberTries_=0;
     805  lastNode_=-999999;
     806  howOften_=100;
     807  decayFactor_ = 0.5;
     808  assert(model.solver());
     809  int numberColumns = model.solver()->getNumCols();
     810  baseSolution_ = new double [numberColumns];
     811  memset(baseSolution_,0,numberColumns*sizeof(double));
     812  whereFrom_=1+8+255*256;
     813  stepSize_=0;
     814  k_=0;
     815  kmax_=0;
     816  nDifferent_=0;
     817}
     818
     819// Destructor
     820CbcHeuristicVND::~CbcHeuristicVND ()
     821{
     822  delete [] baseSolution_;
     823}
     824
     825// Clone
     826CbcHeuristic *
     827CbcHeuristicVND::clone() const
     828{
     829  return new CbcHeuristicVND(*this);
     830}
     831
     832// Assignment operator
     833CbcHeuristicVND &
     834CbcHeuristicVND::operator=( const CbcHeuristicVND& rhs)
     835{
     836  if (this!=&rhs) {
     837    CbcHeuristic::operator=(rhs);
     838    numberSolutions_ = rhs.numberSolutions_;
     839    howOften_ = rhs.howOften_;
     840    numberSuccesses_ = rhs.numberSuccesses_;
     841    numberTries_ = rhs.numberTries_;
     842    lastNode_=rhs.lastNode_;
     843    delete [] baseSolution_;
     844    if (model_&&rhs.baseSolution_) {
     845      int numberColumns = model_->solver()->getNumCols();
     846      baseSolution_ = new double [numberColumns];
     847      memcpy(baseSolution_,rhs.baseSolution_,numberColumns*sizeof(double));
     848    } else {
     849      baseSolution_=NULL;
     850    }
     851    stepSize_=rhs.stepSize_;
     852    k_=rhs.k_;
     853    kmax_=rhs.kmax_;
     854    nDifferent_=rhs.nDifferent_;
     855  }
     856  return *this;
     857}
     858
     859// Create C++ lines to get to current state
     860void
     861CbcHeuristicVND::generateCpp( FILE * fp)
     862{
     863  CbcHeuristicVND other;
     864  fprintf(fp,"0#include \"CbcHeuristicVND.hpp\"\n");
     865  fprintf(fp,"3  CbcHeuristicVND heuristicVND(*cbcModel);\n");
     866  CbcHeuristic::generateCpp(fp,"heuristicVND");
     867  if (howOften_!=other.howOften_)
     868    fprintf(fp,"3  heuristicVND.setHowOften(%d);\n",howOften_);
     869  else
     870    fprintf(fp,"4  heuristicVND.setHowOften(%d);\n",howOften_);
     871  fprintf(fp,"3  cbcModel->addHeuristic(&heuristicVND);\n");
     872}
     873
     874// Copy constructor
     875CbcHeuristicVND::CbcHeuristicVND(const CbcHeuristicVND & rhs)
     876:
     877  CbcHeuristic(rhs),
     878  numberSolutions_(rhs.numberSolutions_),
     879  howOften_(rhs.howOften_),
     880  numberSuccesses_(rhs.numberSuccesses_),
     881  numberTries_(rhs.numberTries_),
     882  lastNode_(rhs.lastNode_)
     883{
     884  if (model_&&rhs.baseSolution_) {
     885    int numberColumns = model_->solver()->getNumCols();
     886    baseSolution_ = new double [numberColumns];
     887    memcpy(baseSolution_,rhs.baseSolution_,numberColumns*sizeof(double));
     888  } else {
     889    baseSolution_=NULL;
     890  }
     891  stepSize_=rhs.stepSize_;
     892  k_=rhs.k_;
     893  kmax_=rhs.kmax_;
     894  nDifferent_=rhs.nDifferent_;
     895}
     896// Resets stuff if model changes
     897void
     898CbcHeuristicVND::resetModel(CbcModel * /*model*/)
     899{
     900  //CbcHeuristic::resetModel(model);
     901  delete [] baseSolution_;
     902  if (model_&&baseSolution_) {
     903    int numberColumns = model_->solver()->getNumCols();
     904    baseSolution_ = new double [numberColumns];
     905    memset(baseSolution_,0,numberColumns*sizeof(double));
     906  } else {
     907    baseSolution_=NULL;
     908  }
     909}
     910/*
     911  First tries setting a variable to better value.  If feasible then
     912  tries setting others.  If not feasible then tries swaps
     913  Returns 1 if solution, 0 if not */
     914int
     915CbcHeuristicVND::solution(double & solutionValue,
     916                         double * betterSolution)
     917{
     918  numCouldRun_++;
     919  int returnCode=0;
     920  const double * bestSolution = model_->bestSolution();
     921  if (!bestSolution)
     922    return 0; // No solution found yet
     923  if (numberSolutions_<model_->getSolutionCount()) {
     924    // new solution - add info
     925    numberSolutions_=model_->getSolutionCount();
     926
     927    int numberIntegers = model_->numberIntegers();
     928    const int * integerVariable = model_->integerVariable();
    779929 
     930    int i;
     931    for (i=0;i<numberIntegers;i++) {
     932      int iColumn = integerVariable[i];
     933      const OsiObject * object = model_->object(i);
     934      // get original bounds
     935      double originalLower;
     936      double originalUpper;
     937      getIntegerInformation( object,originalLower, originalUpper);
     938      double value=bestSolution[iColumn];
     939      if (value<originalLower) {
     940        value=originalLower;
     941      } else if (value>originalUpper) {
     942        value=originalUpper;
     943      }
     944    }
     945  }
     946  int numberNodes=model_->getNodeCount();
     947  if (howOften_==100) {
     948    if (numberNodes<lastNode_+12)
     949      return 0;
     950    // Do at 50 and 100
     951    if ((numberNodes>40&&numberNodes<=50)||(numberNodes>90&&numberNodes<100))
     952      numberNodes=howOften_;
     953  }
     954  if ((numberNodes%howOften_)==0&&(model_->getCurrentPassNumber()==1||
     955                                   model_->getCurrentPassNumber()==999999)) {
     956    lastNode_=model_->getNodeCount();
     957    OsiSolverInterface * solver = model_->solver();
     958
     959    int numberIntegers = model_->numberIntegers();
     960    const int * integerVariable = model_->integerVariable();
     961 
     962    const double * currentSolution = solver->getColSolution();
     963    OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
     964    //const double * colLower = newSolver->getColLower();
     965    //const double * colUpper = newSolver->getColUpper();
     966
     967    double primalTolerance;
     968    solver->getDblParam(OsiPrimalTolerance,primalTolerance);
     969
     970    // Sort on distance
     971    double * distance = new double [numberIntegers];
     972    int * which = new int [numberIntegers];
     973   
     974    int i;
     975    int nFix=0;
     976    double tolerance = 10.0*primalTolerance;
     977    for (i=0;i<numberIntegers;i++) {
     978      int iColumn=integerVariable[i];
     979      const OsiObject * object = model_->object(i);
     980      // get original bounds
     981      double originalLower;
     982      double originalUpper;
     983      getIntegerInformation( object,originalLower, originalUpper);
     984      double valueInt=bestSolution[iColumn];
     985      if (valueInt<originalLower) {
     986        valueInt=originalLower;
     987      } else if (valueInt>originalUpper) {
     988        valueInt=originalUpper;
     989      }
     990      baseSolution_[iColumn]=currentSolution[iColumn];
     991      distance[i]=fabs(currentSolution[iColumn]-valueInt);
     992      which[i]=i;
     993      if (fabs(currentSolution[iColumn]-valueInt)<tolerance)
     994        nFix++;
     995    }
     996    CoinSort_2(distance,distance+numberIntegers,which);
     997    nDifferent_=numberIntegers-nFix;
     998    stepSize_=nDifferent_/10;
     999    k_ = stepSize_;
     1000    //nFix = numberIntegers-stepSize_;
     1001    for (i=0;i<nFix;i++) {
     1002      int j=which[i];
     1003      int iColumn=integerVariable[j];
     1004      const OsiObject * object = model_->object(i);
     1005      // get original bounds
     1006      double originalLower;
     1007      double originalUpper;
     1008      getIntegerInformation( object,originalLower, originalUpper);
     1009      double valueInt=bestSolution[iColumn];
     1010      if (valueInt<originalLower) {
     1011        valueInt=originalLower;
     1012      } else if (valueInt>originalUpper) {
     1013        valueInt=originalUpper;
     1014      }
     1015      double nearest=floor(valueInt+0.5);
     1016      newSolver->setColLower(iColumn,nearest);
     1017      newSolver->setColUpper(iColumn,nearest);
     1018    }
     1019    delete [] distance;
     1020    delete [] which;
     1021    if (nFix>numberIntegers/5) {
     1022      //printf("%d integers have samish value\n",nFix);
     1023      returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
     1024                                         model_->getCutoff(),"CbcHeuristicVND");
     1025      if (returnCode<0)
     1026        returnCode=0; // returned on size
     1027      else
     1028        numRuns_++;
     1029      if ((returnCode&1)!=0)
     1030        numberSuccesses_++;
     1031      //printf("return code %d",returnCode);
     1032      if ((returnCode&2)!=0) {
     1033        // could add cut
     1034        returnCode &= ~2;
     1035        //printf("could add cut with %d elements (if all 0-1)\n",nFix);
     1036      } else {
     1037        //printf("\n");
     1038      }
     1039      numberTries_++;
     1040      if ((numberTries_%10)==0&&numberSuccesses_*3<numberTries_)
     1041        howOften_ += static_cast<int> (howOften_*decayFactor_);
     1042    }
     1043
     1044    delete newSolver;
     1045  }
     1046  return returnCode;
     1047}
     1048// update model
     1049void CbcHeuristicVND::setModel(CbcModel * model)
     1050{
     1051  model_ = model;
     1052  // Get a copy of original matrix
     1053  assert(model_->solver());
     1054  delete [] baseSolution_;
     1055  int numberColumns = model->solver()->getNumCols();
     1056  baseSolution_ = new double [numberColumns];
     1057  memset(baseSolution_,0,numberColumns*sizeof(double));
     1058}
     1059
     1060 
  • trunk/Cbc/src/CbcHeuristicRINS.hpp

    r1221 r1240  
    204204};
    205205
     206/** LocalSearch class
     207 */
     208
     209class CbcHeuristicVND : public CbcHeuristic {
     210public:
     211
     212  // Default Constructor
     213  CbcHeuristicVND ();
     214
     215  /* Constructor with model - assumed before cuts
     216     Initial version does not do Lps
     217  */
     218  CbcHeuristicVND (CbcModel & model);
     219 
     220  // Copy constructor
     221  CbcHeuristicVND ( const CbcHeuristicVND &);
     222   
     223  // Destructor
     224  ~CbcHeuristicVND ();
     225 
     226  /// Clone
     227  virtual CbcHeuristic * clone() const;
     228
     229
     230  /// Assignment operator
     231  CbcHeuristicVND & operator=(const CbcHeuristicVND& rhs);
     232
     233  /// Create C++ lines to get to current state
     234  virtual void generateCpp( FILE * fp) ;
     235
     236  /// Resets stuff if model changes
     237  virtual void resetModel(CbcModel * model);
     238
     239  /// update model (This is needed if cliques update matrix etc)
     240  virtual void setModel(CbcModel * model);
     241 
     242  using CbcHeuristic::solution ;
     243  /** returns 0 if no solution, 1 if valid solution.
     244      Sets solution values if good, sets objective value (only if good)
     245      This does Relaxation Induced Neighborhood Search
     246  */
     247  virtual int solution(double & objectiveValue,
     248                       double * newSolution);
     249  /// This version fixes stuff and does IP
     250  int solutionFix(double & objectiveValue,
     251                  double * newSolution,
     252                  const int * keep);
     253
     254  /// Sets how often to do it
     255  inline void setHowOften(int value)
     256  { howOften_=value;}
     257  /// base solution array so we can set
     258  inline double * baseSolution() const
     259  { return baseSolution_;}
     260
     261protected:
     262  // Data
     263
     264  /// Number of solutions so we can do something at solution
     265  int numberSolutions_;
     266  /// How often to do (code can change)
     267  int howOften_;
     268  /// Number of successes
     269  int numberSuccesses_;
     270  /// Number of tries
     271  int numberTries_;
     272  /// Node when last done
     273  int lastNode_;
     274  /// Step size for decomposition
     275  int stepSize_;
     276  int k_;
     277  int kmax_;
     278  int nDifferent_;
     279  /// Base solution
     280  double * baseSolution_;
     281};
     282
    206283#endif
  • trunk/Cbc/src/CbcModel.cpp

    r1239 r1240  
    609609    } else if (largestObj<smallestObj*5.0&&!parentModel_&&
    610610               !numberContinuousObj&&
     611               !numberGeneralIntegerObj&&
    611612               numberIntegerObj*2<numberColumns) {
    612613      // up priorities on costed
     
    13011302          iColumn = originalObject[iObject]->columnNumber();
    13021303          if (iColumn<0) {
     1304            // already has column numbers changed
    13031305            object_[numberObjects_] = originalObject[iObject]->clone();
     1306#if 0
    13041307            // redo ids etc
    13051308            CbcObject * obj =
     
    13071310            assert (obj);
    13081311            obj->redoSequenceEtc(this,numberColumns,originalColumns);
     1312#endif
    13091313            numberObjects_++;
    13101314          }
     
    25062510  //solverCharacteristics_->setSolver(solver_);
    25072511  if (feasible) {
    2508     //#define HOTSTART 1
    2509 #if HOTSTART
     2512    //#define HOTSTART -1
     2513#if HOTSTART<0
     2514    if (bestSolution_&&!parentModel_&&!hotstartSolution_&&
     2515        (moreSpecialOptions_&1024)!=0) {
     2516      // Set priorities so only branch on ones we need to
     2517      // use djs and see if only few branches needed
     2518#ifndef NDEBUG
     2519      double integerTolerance = getIntegerTolerance() ;
     2520#endif
     2521      bool possible=true;
     2522      const double * saveLower = continuousSolver_->getColLower();
     2523      const double * saveUpper = continuousSolver_->getColUpper();
     2524      for (int i=0;i<numberObjects_;i++) {
     2525        const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object_[i]);
     2526        if (thisOne) {
     2527          int iColumn = thisOne->columnNumber();
     2528          if (saveUpper[iColumn]>saveLower[iColumn]+1.5) {
     2529            possible=false;
     2530            break;
     2531          }
     2532        } else {
     2533          possible=false;
     2534          break;
     2535        }
     2536      }
     2537      if (possible) {
     2538        OsiSolverInterface * solver = continuousSolver_->clone();
     2539        int numberColumns = solver->getNumCols();
     2540        for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
     2541          double value = bestSolution_[iColumn] ;
     2542          value = CoinMax(value, saveLower[iColumn]) ;
     2543          value = CoinMin(value, saveUpper[iColumn]) ;
     2544          value = floor(value+0.5);
     2545          if (solver->isInteger(iColumn)) {
     2546            solver->setColLower(iColumn,value);
     2547            solver->setColUpper(iColumn,value);
     2548          }
     2549        }
     2550        solver->setHintParam(OsiDoDualInResolve,false,OsiHintTry);
     2551        // objlim and all slack
     2552        double direction = solver->getObjSense();
     2553        solver->setDblParam(OsiDualObjectiveLimit,1.0e50*direction);
     2554        CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis *> (solver->getEmptyWarmStart());
     2555        solver->setWarmStart(basis);
     2556        delete basis;
     2557        bool changed=true;
     2558        hotstartPriorities_ = new int [numberColumns];
     2559        for (int iColumn=0;iColumn<numberColumns;iColumn++)
     2560          hotstartPriorities_[iColumn]=1;
     2561        while (changed) {
     2562          changed=false;
     2563          solver->resolve();
     2564          if (!solver->isProvenOptimal()) {
     2565            possible=false;
     2566            break;
     2567          }
     2568          const double * dj = solver->getReducedCost();
     2569          const double * colLower = solver->getColLower();
     2570          const double * colUpper = solver->getColUpper();
     2571          const double * solution = solver->getColSolution();
     2572          int nAtLbNatural=0;
     2573          int nAtUbNatural=0;
     2574          int nZeroDj=0;
     2575          int nForced=0;
     2576          for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++) {
     2577            double value = solution[iColumn] ;
     2578            value = CoinMax(value, saveLower[iColumn]) ;
     2579            value = CoinMin(value, saveUpper[iColumn]) ;
     2580            if (solver->isInteger(iColumn)) {
     2581              assert(fabs(value-solution[iColumn]) <= integerTolerance) ;
     2582              if (hotstartPriorities_[iColumn]==1) {
     2583                if (dj[iColumn]<-1.0e-6) {
     2584                  // negative dj
     2585                  if (saveUpper[iColumn]==colUpper[iColumn]) {
     2586                    nAtUbNatural++;
     2587                    hotstartPriorities_[iColumn]=2;
     2588                    solver->setColLower(iColumn,saveLower[iColumn]);
     2589                    solver->setColUpper(iColumn,saveUpper[iColumn]);
     2590                  } else {
     2591                    nForced++;
     2592                  }
     2593                } else if (dj[iColumn]>1.0e-6) {
     2594                  // positive dj
     2595                  if (saveLower[iColumn]==colLower[iColumn]) {
     2596                    nAtLbNatural++;
     2597                    hotstartPriorities_[iColumn]=2;
     2598                    solver->setColLower(iColumn,saveLower[iColumn]);
     2599                    solver->setColUpper(iColumn,saveUpper[iColumn]);
     2600                  } else {
     2601                    nForced++;
     2602                  }
     2603                } else {
     2604                  // zero dj
     2605                  nZeroDj++;
     2606                }
     2607              }
     2608            }
     2609          }
     2610#ifdef CLP_INVESTIGATE
     2611          printf("%d forced, %d naturally at lower, %d at upper - %d zero dj\n",
     2612                 nForced,nAtLbNatural,nAtUbNatural,nZeroDj);
     2613#endif
     2614          if (nAtLbNatural||nAtUbNatural) {
     2615            changed=true;
     2616          } else {
     2617            if (nForced+nZeroDj>50||
     2618                (nForced+nZeroDj)*4>numberIntegers_)
     2619              possible=false;
     2620          }
     2621        }
     2622        delete solver;
     2623      }
     2624      if (possible) {
     2625        setHotstartSolution(bestSolution_);
     2626        if (!saveCompare) {
     2627          // create depth first comparison
     2628          saveCompare = nodeCompare_;
     2629          // depth first
     2630          nodeCompare_ = new CbcCompareDepth();
     2631          tree_->setComparison(*nodeCompare_) ;
     2632        }
     2633      } else {
     2634        delete [] hotstartPriorities_;
     2635        hotstartPriorities_=NULL;
     2636      }
     2637    }
     2638#endif
     2639#if HOTSTART>0
    25102640    if (hotstartSolution_&&!hotstartPriorities_) {
    25112641      // Set up hot start
     
    94339563          dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
    94349564        if (obj1) {
    9435           assert (obj1->downShadowPrice()>0.0);
     9565          //assert (obj1->downShadowPrice()>0.0);
    94369566#define P_FACTOR 1.0
    94379567#if 1
     
    97259855    double smallDown = 0.0001*(downSum/static_cast<double> (numberIntegers));
    97269856    double smallUp = 0.0001*(upSum/static_cast<double> (numberIntegers));
    9727 #define PSEUDO_FACTOR 1.0e-1
     9857#define PSEUDO_FACTOR 5.0e-1
     9858    double pseudoFactor=PSEUDO_FACTOR;
     9859    //if (!numberNodes_)
     9860    //pseudoFactor=0.0;
    97289861    for (int i=0;i<numberObjects_;i++) {
    97299862      CbcSimpleIntegerDynamicPseudoCost * obj1 =
     
    97339866        double upPseudoCost = obj1->upDynamicPseudoCost();
    97349867        double saveUp = upPseudoCost;
    9735         upPseudoCost = CoinMax(PSEUDO_FACTOR*upPseudoCost,smallUp);
     9868        upPseudoCost = CoinMax(pseudoFactor*upPseudoCost,smallUp);
    97369869        upPseudoCost = CoinMax(upPseudoCost,up[iColumn]);
    97379870        upPseudoCost = CoinMax(upPseudoCost,0.001*down[iColumn]);
     
    97429875        double downPseudoCost = obj1->downDynamicPseudoCost();
    97439876        double saveDown = downPseudoCost;
    9744         downPseudoCost = CoinMax(PSEUDO_FACTOR*downPseudoCost,smallDown);
     9877        downPseudoCost = CoinMax(pseudoFactor*downPseudoCost,smallDown);
    97459878        downPseudoCost = CoinMax(downPseudoCost,down[iColumn]);
    97469879        downPseudoCost = CoinMax(downPseudoCost,0.001*up[iColumn]);
     
    1050510638    if ((solver_->isProvenOptimal()||(specialOptions_&4)!=0) && objectiveValue <= cutoff) {
    1050610639      memcpy(solution ,solver_->getColSolution(),numberColumns*sizeof(double)) ;
    10507      
    1050810640      int iColumn;
    1050910641#ifndef NDEBUG
     
    1327113403        if (node->depth()>=go_fathom &&(specialOptions_&2048)==0
    1327213404            //if (node->depth()>=FATHOM_BIAS-fastNodeDepth_&&!parentModel_
    13273             &&numberNodes_>=numberNodesBeforeFathom) {
     13405            &&numberNodes_>=numberNodesBeforeFathom&&!hotstartSolution_) {
    1327413406#ifndef COIN_HAS_CPX
    1327513407          specialOptions_ &= ~16384;
     
    1330413436            fillPseudoCosts(down,up,priority,numberDown,numberUp,
    1330513437                            numberDownInfeasible,numberUpInfeasible);
     13438            // See if all priorities same
     13439            bool allSame=true;
     13440            int kPriority=priority[0];
     13441            for (int i=1;i<numberIntegers_;i++) {
     13442              if (kPriority!=priority[i]) {
     13443                allSame=false;
     13444                break;
     13445              }
     13446            }
     13447            ClpSimplex * simplex = clpSolver->getModelPtr();
     13448            if (allSame&&false) {
     13449              // change priorities on general
     13450              const double * lower = simplex->columnLower();
     13451              const double * upper = simplex->columnUpper();
     13452              for (int i=0;i<numberIntegers_;i++) {
     13453                int iColumn = integerVariable_[i];
     13454                if (upper[iColumn]>lower[iColumn]+1.1)
     13455                  priority[i]=kPriority+1;
     13456              }
     13457            }
    1330613458            info->fillPseudoCosts(down,up,priority,numberDown,numberUp,
    1330713459                                  numberDownInfeasible,
     
    1332013472            OsiHintStrength strength;
    1332113473            solver_->getHintParam(OsiDoReducePrint,takeHint,strength);
    13322             ClpSimplex * simplex = clpSolver->getModelPtr();
    1332313474            //printf("mod cutoff %g solver %g offset %g\n",
    1332413475            //   getCutoff(),simplex->dualObjectiveLimit(),simplex->objectiveOffset());
     
    1494215093  }
    1494315094#ifdef CLP_INVESTIGATE
    14944   printf("Before fathom %d not trusted out of %d\n",
    14945          numberNot,numberIntegers_);
     15095  if (priority)
     15096    printf("Before fathom %d not trusted out of %d\n",
     15097           numberNot,numberIntegers_);
    1494615098#endif
    1494715099  delete [] back;
  • trunk/Cbc/src/CbcModel.hpp

    r1232 r1240  
    15931593  { return (specialOptions_&16)==0;}
    15941594  /** Set more special options
    1595       at present bottom 3 bits used for shadow price mode
     1595      at present bottom 6 bits used for shadow price mode
     1596      1024 for experimental hotstart
    15961597  */
    15971598  inline void setMoreSpecialOptions(int value)
  • trunk/Cbc/src/CbcNode.cpp

    r1224 r1240  
    14061406  int numberStrong=model->numberStrong();
    14071407  // switch off strong if hotstart
    1408   if (model->hotstartSolution())
    1409     numberStrong=0;
    1410   int numberStrongDone=0;
    1411   int numberUnfinished=0;
    1412   int numberStrongInfeasible=0;
    1413   int numberStrongIterations=0;
    1414   int saveNumberStrong=numberStrong;
     1408  const double * hotstartSolution = model->hotstartSolution();
     1409  const int * hotstartPriorities = model->hotstartPriorities();
    14151410  int numberObjects = model->numberObjects();
    1416   bool checkFeasibility = numberObjects>model->numberIntegers();
    1417   int maximumStrong = CoinMax(CoinMin(numberStrong,numberObjects),1);
    14181411  int numberColumns = model->getNumCols();
    14191412  double * saveUpper = new double[numberColumns];
    14201413  double * saveLower = new double[numberColumns];
     1414  for (i=0;i<numberColumns;i++) {
     1415    saveLower[i] = lower[i];
     1416    saveUpper[i] = upper[i];
     1417  }
    14211418 
    14221419  // Save solution in case heuristics need good solution later
     
    14251422  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
    14261423  model->reserveCurrentSolution(saveSolution);
     1424  if (hotstartSolution) {
     1425    numberStrong=0;
     1426    if((model->moreSpecialOptions()&1024)!=0) {
     1427      int nBad=0;
     1428      int nUnsat=0;
     1429      int nDiff=0;
     1430      for (int i=0;i<numberObjects;i++) {
     1431        OsiObject * object = model->modifiableObject(i);
     1432        const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
     1433        if (thisOne) {
     1434          int iColumn = thisOne->columnNumber();
     1435          double targetValue = hotstartSolution[iColumn];
     1436          double value = saveSolution[iColumn];
     1437          if (fabs(value-floor(value+0.5))>1.0e-6) {
     1438            nUnsat++;
     1439#ifdef CLP_INVESTIGATE
     1440            printf("H %d is %g target %g\n",iColumn,value,targetValue);
     1441#endif
     1442          } else if (fabs(targetValue-value)>1.0e-6) {
     1443            nDiff++;
     1444          }
     1445          if (targetValue<saveLower[iColumn]||
     1446              targetValue>saveUpper[iColumn]) {
     1447#ifdef CLP_INVESTIGATE
     1448            printf("%d has target %g and current bounds %g and %g\n",
     1449                   iColumn,targetValue,saveLower[iColumn],saveUpper[iColumn]);
     1450#endif
     1451            nBad++;
     1452          }
     1453        }
     1454      }
     1455#ifdef CLP_INVESTIGATE
     1456      printf("Hot %d unsatisfied, %d outside limits, %d different\n",
     1457             nUnsat,nBad,nDiff);
     1458#endif
     1459      if (nBad) {
     1460        // switch off as not possible
     1461          hotstartSolution=NULL;
     1462          model->setHotstartSolution(NULL,NULL);
     1463          usefulInfo.hotstartSolution_=NULL;
     1464      }
     1465    }
     1466  }
     1467  int numberStrongDone=0;
     1468  int numberUnfinished=0;
     1469  int numberStrongInfeasible=0;
     1470  int numberStrongIterations=0;
     1471  int saveNumberStrong=numberStrong;
     1472  bool checkFeasibility = numberObjects>model->numberIntegers();
     1473  int maximumStrong = CoinMax(CoinMin(numberStrong,numberObjects),1);
    14271474  /*
    14281475    Get a branching decision object. Use the default decision criteria unless
     
    14361483  decision->initialize(model);
    14371484  CbcStrongInfo * choice = new CbcStrongInfo[maximumStrong];
    1438   for (i=0;i<numberColumns;i++) {
    1439     saveLower[i] = lower[i];
    1440     saveUpper[i] = upper[i];
    1441   }
    14421485  // May go round twice if strong branching fixes all local candidates
    14431486  bool finished=false;
     
    14611504      //                      numberIntegerInfeasibilities,
    14621505      //                      numberObjectInfeasibilities);
    1463       const double * hotstartSolution = model->hotstartSolution();
    1464       const int * hotstartPriorities = model->hotstartPriorities();
    1465      
    14661506      // Some objects may compute an estimate of best solution from here
    14671507      estimatedDegradation=0.0;
     
    15171557                */
    15181558                if (fabs(value-targetValue)>integerTolerance) {
    1519                   infeasibility = fabs(value-targetValue);
     1559                  //if (infeasibility>0.01)
     1560                  //infeasibility = fabs(1.0e6-fabs(value-targetValue));
     1561                  //else
     1562                  infeasibility = fabs(value-targetValue);
    15201563                  //if (targetValue==1.0)
    15211564                  //infeasibility += 1.0;
     
    16031646        hotstartSolution=NULL;
    16041647        model->setHotstartSolution(NULL,NULL);
     1648        usefulInfo.hotstartSolution_=NULL;
    16051649      }
    16061650      if (numberUnsatisfied_) {
    16071651        // some infeasibilities - go to next steps
     1652#ifdef CLP_INVESTIGATE
     1653        if (hotstartSolution) {
     1654          int k=choice[0].objectNumber;
     1655          OsiObject * object = model->modifiableObject(k);
     1656          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
     1657          assert (thisOne);
     1658          int iColumn = thisOne->columnNumber();
     1659          double targetValue = hotstartSolution[iColumn];
     1660          double value = saveSolution[iColumn];
     1661          printf("Branch on %d has target %g (value %g) and current bounds %g and %g\n",
     1662                 iColumn,targetValue,value,saveLower[iColumn],saveUpper[iColumn]);
     1663        }
     1664#endif
    16081665        break;
    16091666      } else if (!iPass) {
     
    23922449    depth_ = 0;
    23932450  // Go to other choose if hot start
    2394   if (model->hotstartSolution())
     2451  if (model->hotstartSolution()&&
     2452      (((model->moreSpecialOptions()&1024)==0)||false))
    23952453    return -3;
    23962454  delete branch_;
     
    24912549  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
    24922550  model->reserveCurrentSolution(saveSolution);
     2551  const double * hotstartSolution = model->hotstartSolution();
     2552  const int * hotstartPriorities = model->hotstartPriorities();
     2553  double integerTolerance =
     2554    model->getDblParam(CbcModel::CbcIntegerTolerance);
     2555  if (hotstartSolution) {
     2556    if((model->moreSpecialOptions()&1024)!=0) {
     2557      int nBad=0;
     2558      int nUnsat=0;
     2559      int nDiff=0;
     2560      for (int i=0;i<numberObjects;i++) {
     2561        OsiObject * object = model->modifiableObject(i);
     2562        const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
     2563        if (thisOne) {
     2564          int iColumn = thisOne->columnNumber();
     2565          double targetValue = hotstartSolution[iColumn];
     2566          double value = saveSolution[iColumn];
     2567          if (fabs(value-floor(value+0.5))>1.0e-6) {
     2568            nUnsat++;
     2569#ifdef CLP_INVESTIGATE
     2570            printf("H %d is %g target %g\n",iColumn,value,targetValue);
     2571#endif
     2572          } else if (fabs(targetValue-value)>1.0e-6) {
     2573            nDiff++;
     2574          }
     2575          if (targetValue<saveLower[iColumn]||
     2576              targetValue>saveUpper[iColumn]) {
     2577#ifdef CLP_INVESTIGATE
     2578            printf("%d has target %g and current bounds %g and %g\n",
     2579                   iColumn,targetValue,saveLower[iColumn],saveUpper[iColumn]);
     2580#endif
     2581            nBad++;
     2582          }
     2583        }
     2584      }
     2585#ifdef CLP_INVESTIGATE
     2586      printf("Hot %d unsatisfied, %d outside limits, %d different\n",
     2587             nUnsat,nBad,nDiff);
     2588#endif
     2589      if (nBad) {
     2590        // switch off as not possible
     2591          hotstartSolution=NULL;
     2592          model->setHotstartSolution(NULL,NULL);
     2593          usefulInfo.hotstartSolution_=NULL;
     2594      }
     2595    }
     2596  }
    24932597  /*
    24942598    Get a branching decision object. Use the default dynamic decision criteria unless
     
    25822686#endif
    25832687  // Get average up and down costs
    2584   double averageUp=0.0;
    2585   double averageDown=0.0;
    25862688  {
     2689    double averageUp=0.0;
     2690    double averageDown=0.0;
    25872691    int numberUp=0;
    25882692    int numberDown=0;
     
    26512755  if (useShadow) {
    26522756    // pseudo shadow prices
    2653     model->pseudoShadow(model->moreSpecialOptions()>>3);
     2757    model->pseudoShadow((model->moreSpecialOptions()>>3)&63);
    26542758  }
    26552759#ifdef DEPRECATED_STRATEGY
     
    27622866      */
    27632867      int problemType = model->problemType();
     2868      bool canDoOneHot=false;
    27642869      for (i=0;i<numberObjects;i++) {
    27652870        OsiObject * object = model->modifiableObject(i);
     
    27692874        double infeasibility = object->infeasibility(&usefulInfo,preferredWay);
    27702875        int priorityLevel = object->priority();
     2876        if (hotstartSolution) {
     2877          // we are doing hot start
     2878          const CbcSimpleInteger * thisOne = dynamic_cast <const CbcSimpleInteger *> (object);
     2879          if (thisOne) {
     2880            int iColumn = thisOne->columnNumber();
     2881            bool canDoThisHot=true;
     2882            double targetValue = hotstartSolution[iColumn];
     2883            if (saveUpper[iColumn]>saveLower[iColumn]) {
     2884              double value = saveSolution[iColumn];
     2885              if (hotstartPriorities)
     2886                priorityLevel=hotstartPriorities[iColumn];
     2887              //double originalLower = thisOne->originalLower();
     2888              //double originalUpper = thisOne->originalUpper();
     2889              // switch off if not possible
     2890              if (targetValue>=saveLower[iColumn]&&targetValue<=saveUpper[iColumn]) {
     2891                /* priority outranks rest always if negative
     2892                   otherwise can be downgraded if at correct level.
     2893                   Infeasibility may be increased to choose 1.0 values first.
     2894                   choose one near wanted value
     2895                */
     2896                if (fabs(value-targetValue)>integerTolerance) {
     2897                  //if (infeasibility>0.01)
     2898                  //infeasibility = fabs(1.0e6-fabs(value-targetValue));
     2899                  //else
     2900                  infeasibility = fabs(value-targetValue);
     2901                  //if (targetValue==1.0)
     2902                  //infeasibility += 1.0;
     2903                  if (value>targetValue) {
     2904                    preferredWay=-1;
     2905                  } else {
     2906                    preferredWay=1;
     2907                  }
     2908                  priorityLevel = CoinAbs(priorityLevel);
     2909                } else if (priorityLevel<0) {
     2910                  priorityLevel = CoinAbs(priorityLevel);
     2911                  if (targetValue==saveLower[iColumn]) {
     2912                    infeasibility = integerTolerance+1.0e-12;
     2913                    preferredWay=-1;
     2914                  } else if (targetValue==saveUpper[iColumn]) {
     2915                    infeasibility = integerTolerance+1.0e-12;
     2916                    preferredWay=1;
     2917                  } else {
     2918                    // can't
     2919                    priorityLevel += 10000000;
     2920                    canDoThisHot=false;
     2921                  }
     2922                } else {
     2923                  priorityLevel += 10000000;
     2924                  canDoThisHot=false;
     2925                }
     2926              } else {
     2927                // switch off if not possible
     2928                canDoThisHot=false;
     2929              }
     2930              if (canDoThisHot)
     2931                canDoOneHot=true;
     2932            } else if (targetValue<saveLower[iColumn]||targetValue>saveUpper[iColumn]) {
     2933            }
     2934          } else {
     2935            priorityLevel += 10000000;
     2936          }
     2937        }
    27712938#define ZERO_ONE 0
    27722939#define ZERO_FAKE 1.0e20;
     
    29203087                 numberUnsatisProbed,numberUnsatisNotProbed);
    29213088        // some infeasibilities - go to next steps
     3089        if (!canDoOneHot&&hotstartSolution) {
     3090          // switch off as not possible
     3091          hotstartSolution=NULL;
     3092          model->setHotstartSolution(NULL,NULL);
     3093          usefulInfo.hotstartSolution_=NULL;
     3094        }
    29223095        break;
    29233096      } else if (!iPass) {
  • trunk/Cbc/src/CbcSolver.cpp

    r1236 r1240  
    34143414        }
    34153415      }
     3416      if (accumulate) {
     3417        heuristic4.setAccumulate(accumulate);
     3418        if (printStuff) {
     3419          if (accumulate) {
     3420            sprintf(generalPrint,"Accumulate of %d",accumulate);
     3421            generalMessageHandler->message(CBC_GENERAL,generalMessages)
     3422              << generalPrint
     3423              <<CoinMessageEol;
     3424          }
     3425        }
     3426      }
    34163427      if (r) {
    34173428        // also set increment
     
    34223433          increment = fakeIncrement;
    34233434        heuristic4.setAbsoluteIncrement(increment);
    3424         heuristic4.setAccumulate(accumulate);
    34253435        heuristic4.setMaximumRetries(r+1);
    34263436        if (printStuff) {
    34273437          if (increment) {
    34283438            sprintf(generalPrint,"Increment of %g",increment);
    3429             generalMessageHandler->message(CBC_GENERAL,generalMessages)
    3430               << generalPrint
    3431               <<CoinMessageEol;
    3432           }
    3433           if (accumulate) {
    3434             sprintf(generalPrint,"Accumulate of %d",accumulate);
    34353439            generalMessageHandler->message(CBC_GENERAL,generalMessages)
    34363440              << generalPrint
     
    34963500  }
    34973501  if (useRENS>=kType&&useRENS<=kType+1) {
     3502#if 1
    34983503    CbcHeuristicRENS heuristic6(*model);
    34993504    heuristic6.setHeuristicName("RENS");
     
    35023507    int nodes []={-2,50,50,50,200,1000,10000};
    35033508    heuristic6.setNumberNodes(nodes[useRENS]);
     3509#else
     3510    CbcHeuristicVND heuristic6(*model);
     3511    heuristic6.setHeuristicName("VND");
     3512    heuristic5.setFractionSmall(0.5);
     3513    heuristic5.setDecayFactor(5.0);
     3514#endif
    35043515    model->addHeuristic(&heuristic6) ;
    35053516    anyToDo=true;
     
    54675478                  iStat = 3;
    54685479                  pos += sprintf(buf+pos,"stopped on ctrl-c,");
     5480                } else if (iStat==6) {
     5481                  // bab infeasible
     5482                  pos += sprintf(buf+pos,"integer infeasible,");
     5483                  iStat=1;
    54695484                } else {
    54705485                  pos += sprintf(buf+pos,"status unknown,");
     
    88408855                    iStat = 3;
    88418856                    pos += sprintf(buf+pos,"stopped on ctrl-c,");
     8857                  } else if (iStat==6) {
     8858                    // bab infeasible
     8859                    pos += sprintf(buf+pos,"integer infeasible,");
     8860                    iStat=1;
    88428861                  } else {
    88438862                    pos += sprintf(buf+pos,"status unknown,");
     
    1019210211                  iStat = 3;
    1019310212                  pos += sprintf(buf+pos,"stopped on ctrl-c,");
     10213                } else if (iStat==6) {
     10214                  // bab infeasible
     10215                  pos += sprintf(buf+pos,"integer infeasible,");
     10216                  iStat=1;
    1019410217                } else {
    1019510218                  pos += sprintf(buf+pos,"status unknown,");
  • trunk/Cbc/src/CbcTreeLocal.cpp

    r1173 r1240  
    163163      if (newSolutionValue<bestCutoff_) {
    164164        model_->setBestSolution(CBC_ROUNDING,newSolutionValue,solution);
    165         bestCutoff_=newSolutionValue;
     165        bestCutoff_=model_->getCutoff();
    166166        // save as best solution
    167167        memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
     
    324324  int goodSolution=createCut(solution,cut_);
    325325  if (goodSolution>=0) {
    326     bestCutoff_=solutionValue;
     326    bestCutoff_=CoinMin(solutionValue,model_->getCutoff());
    327327  } else {
    328328    model_=NULL;
     
    475475        if (nextStrong_) {
    476476          diversification_++;
    477           model_->setCutoff(1.0e50);
     477          // cut is valid so don't model_->setCutoff(1.0e50);
    478478          searchType_=0;
    479479        }
     
    880880
    881881
    882 
     882CbcTreeVariable::CbcTreeVariable()
     883  : localNode_(NULL),
     884    bestSolution_(NULL),
     885    savedSolution_(NULL),
     886    saveNumberSolutions_(0),
     887    model_(NULL),
     888    originalLower_(NULL),
     889    originalUpper_(NULL),
     890    range_(0),
     891    typeCuts_(-1),
     892    maxDiversification_(0),
     893    diversification_(0),
     894    nextStrong_(false),
     895    rhs_(0.0),
     896    savedGap_(0.0),
     897    bestCutoff_(0.0),
     898    timeLimit_(0),
     899    startTime_(0),
     900    nodeLimit_(0),
     901    startNode_(-1),
     902    searchType_(-1),
     903    refine_(false)
     904{
     905 
     906}
     907/* Constructor with solution.
     908   range is upper bound on difference from given solution.
     909   maxDiversification is maximum number of diversifications to try
     910   timeLimit is seconds in subTree
     911   nodeLimit is nodes in subTree
     912*/
     913CbcTreeVariable::CbcTreeVariable(CbcModel * model,const double * solution ,
     914                                  int range,int typeCuts,int maxDiversification,
     915                                  int timeLimit, int nodeLimit,bool refine)
     916  : localNode_(NULL),
     917    bestSolution_(NULL),
     918    savedSolution_(NULL),
     919    saveNumberSolutions_(0),
     920    model_(model),
     921    originalLower_(NULL),
     922    originalUpper_(NULL),
     923    range_(range),
     924    typeCuts_(typeCuts),
     925    maxDiversification_(maxDiversification),
     926    diversification_(0),
     927    nextStrong_(false),
     928    rhs_(0.0),
     929    savedGap_(0.0),
     930    bestCutoff_(0.0),
     931    timeLimit_(timeLimit),
     932    startTime_(0),
     933    nodeLimit_(nodeLimit),
     934    startNode_(-1),
     935    searchType_(-1),
     936    refine_(refine)
     937{
     938
     939  OsiSolverInterface * solver = model_->solver();
     940  const double * lower = solver->getColLower();
     941  const double * upper = solver->getColUpper();
     942  //const double * solution = solver->getColSolution();
     943  //const double * objective = solver->getObjCoefficients();
     944  double primalTolerance;
     945  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
     946
     947  // Get increment
     948  model_->analyzeObjective();
     949
     950  {
     951    // needed to sync cutoffs
     952    double value ;
     953    solver->getDblParam(OsiDualObjectiveLimit,value) ;
     954    model_->setCutoff(value * solver->getObjSense());
     955  }
     956  bestCutoff_ = model_->getCutoff();
     957  // save current gap
     958  savedGap_ = model_->getDblParam(CbcModel::CbcAllowableGap);
     959
     960  // make sure integers found
     961  model_->findIntegers(false);
     962  int numberIntegers = model_->numberIntegers();
     963  const int * integerVariable = model_->integerVariable();
     964  int i;
     965  double direction = solver->getObjSense();
     966  double newSolutionValue=1.0e50;
     967  if (solution) {
     968    // copy solution
     969    solver->setColSolution(solution);
     970    newSolutionValue = direction*solver->getObjValue();
     971  }
     972  originalLower_=new double [numberIntegers];
     973  originalUpper_=new double [numberIntegers];
     974  bool all01=true;
     975  int number01=0;
     976  for (i=0;i<numberIntegers;i++) {
     977    int iColumn = integerVariable[i];
     978    originalLower_[i]=lower[iColumn];
     979    originalUpper_[i]=upper[iColumn];
     980    if (upper[iColumn]-lower[iColumn]>1.5)
     981      all01=false;
     982    else if (upper[iColumn]-lower[iColumn]==1.0)
     983      number01++;
     984  }
     985  if (all01 && !typeCuts_)
     986    typeCuts_=1; // may as well so we don't have to deal with refine
     987  if (!number01&&!typeCuts_) {
     988    if (model_->messageHandler()->logLevel()>0)
     989      printf("** No 0-1 variables and local search only on 0-1 - switching off\n");
     990    typeCuts_=-1;
     991  } else {
     992    if (model_->messageHandler()->logLevel()>0) {
     993      std::string type;
     994      if (all01) {
     995        printf("%d 0-1 variables normal local  cuts\n",
     996               number01);
     997      } else if (typeCuts_) {
     998        printf("%d 0-1 variables, %d other - general integer local cuts\n",
     999               number01,numberIntegers-number01);
     1000      } else {
     1001        printf("%d 0-1 variables, %d other - local cuts but just on 0-1 variables\n",
     1002               number01,numberIntegers-number01);
     1003      }
     1004      printf("maximum diversifications %d, initial cutspace %d, max time %d seconds, max nodes %d\n",
     1005             maxDiversification_,range_,timeLimit_,nodeLimit_);
     1006    }
     1007  }
     1008  int numberColumns = model_->getNumCols();
     1009  savedSolution_ = new double [numberColumns];
     1010  memset(savedSolution_,0,numberColumns*sizeof(double));
     1011  if (solution) {
     1012    rhs_=range_;
     1013    // Check feasible
     1014    int goodSolution=createCut(solution,cut_);
     1015    if (goodSolution>=0) {
     1016      for (i=0;i<numberIntegers;i++) {
     1017        int iColumn = integerVariable[i];
     1018        double value = floor(solution[iColumn]+0.5);
     1019        // fix so setBestSolution will work
     1020        solver->setColLower(iColumn,value);
     1021        solver->setColUpper(iColumn,value);
     1022      }
     1023      model_->reserveCurrentSolution();
     1024      // Create cut and get total gap
     1025      if (newSolutionValue<bestCutoff_) {
     1026        model_->setBestSolution(CBC_ROUNDING,newSolutionValue,solution);
     1027        bestCutoff_=model_->getCutoff();
     1028        // save as best solution
     1029        memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
     1030      }
     1031      for (i=0;i<numberIntegers;i++) {
     1032        int iColumn = integerVariable[i];
     1033        // restore bounds
     1034        solver->setColLower(iColumn,originalLower_[i]);
     1035        solver->setColUpper(iColumn,originalUpper_[i]);
     1036      }
     1037      // make sure can't stop on gap
     1038      model_->setDblParam(CbcModel::CbcAllowableGap,-1.0e50);
     1039    } else {
     1040      model_=NULL;
     1041    }
     1042  } else {
     1043    // no solution
     1044    rhs_=1.0e50;
     1045    // make sure can't stop on gap
     1046    model_->setDblParam(CbcModel::CbcAllowableGap,-1.0e50);
     1047  }
     1048}
     1049CbcTreeVariable::~CbcTreeVariable()
     1050{
     1051  delete [] originalLower_;
     1052  delete [] originalUpper_;
     1053  delete [] bestSolution_;
     1054  delete [] savedSolution_;
     1055  delete localNode_;
     1056}
     1057// Copy constructor
     1058CbcTreeVariable::CbcTreeVariable ( const CbcTreeVariable & rhs)
     1059  :CbcTree(rhs),
     1060   saveNumberSolutions_(rhs.saveNumberSolutions_),
     1061   model_(rhs.model_),
     1062   range_(rhs.range_),
     1063   typeCuts_(rhs.typeCuts_),
     1064   maxDiversification_(rhs.maxDiversification_),
     1065   diversification_(rhs.diversification_),
     1066   nextStrong_(rhs.nextStrong_),
     1067   rhs_(rhs.rhs_),
     1068   savedGap_(rhs.savedGap_),
     1069   bestCutoff_(rhs.bestCutoff_),
     1070   timeLimit_(rhs.timeLimit_),
     1071   startTime_(rhs.startTime_),
     1072   nodeLimit_(rhs.nodeLimit_),
     1073   startNode_(rhs.startNode_),
     1074   searchType_(rhs.searchType_),
     1075   refine_(rhs.refine_)
     1076
     1077  cut_=rhs.cut_;
     1078  fixedCut_=rhs.fixedCut_;
     1079  if (rhs.localNode_)
     1080    localNode_=new CbcNode(*rhs.localNode_);
     1081  else
     1082    localNode_=NULL;
     1083  if (rhs.originalLower_) {
     1084    int numberIntegers = model_->numberIntegers();
     1085    originalLower_=new double [numberIntegers];
     1086    memcpy(originalLower_,rhs.originalLower_,numberIntegers*sizeof(double));
     1087    originalUpper_=new double [numberIntegers];
     1088    memcpy(originalUpper_,rhs.originalUpper_,numberIntegers*sizeof(double));
     1089  } else {
     1090    originalLower_=NULL;
     1091    originalUpper_=NULL;
     1092  }
     1093  if (rhs.bestSolution_) {
     1094    int numberColumns = model_->getNumCols();
     1095    bestSolution_ = new double [numberColumns];
     1096    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
     1097  } else {
     1098    bestSolution_=NULL;
     1099  }
     1100  if (rhs.savedSolution_) {
     1101    int numberColumns = model_->getNumCols();
     1102    savedSolution_ = new double [numberColumns];
     1103    memcpy(savedSolution_,rhs.savedSolution_,numberColumns*sizeof(double));
     1104  } else {
     1105    savedSolution_=NULL;
     1106  }
     1107}
     1108//----------------------------------------------------------------
     1109// Assignment operator
     1110//-------------------------------------------------------------------
     1111CbcTreeVariable &
     1112CbcTreeVariable::operator=(const CbcTreeVariable& rhs)
     1113{
     1114  if (this != &rhs) {
     1115    CbcTree::operator=(rhs);
     1116    saveNumberSolutions_ = rhs.saveNumberSolutions_;
     1117    cut_=rhs.cut_;
     1118    fixedCut_=rhs.fixedCut_;
     1119    delete localNode_;
     1120    if (rhs.localNode_)
     1121      localNode_=new CbcNode(*rhs.localNode_);
     1122    else
     1123      localNode_=NULL;
     1124    model_ = rhs.model_;
     1125    range_ = rhs.range_;
     1126    typeCuts_ = rhs.typeCuts_;
     1127    maxDiversification_ = rhs.maxDiversification_;
     1128    diversification_ = rhs.diversification_;
     1129    nextStrong_ = rhs.nextStrong_;
     1130    rhs_ = rhs.rhs_;
     1131    savedGap_ = rhs.savedGap_;
     1132    bestCutoff_ = rhs.bestCutoff_;
     1133    timeLimit_ = rhs.timeLimit_;
     1134    startTime_ = rhs.startTime_;
     1135    nodeLimit_ = rhs.nodeLimit_;
     1136    startNode_ = rhs.startNode_;
     1137    searchType_ = rhs.searchType_;
     1138    refine_ = rhs.refine_;
     1139    delete [] originalLower_;
     1140    delete [] originalUpper_;
     1141    if (rhs.originalLower_) {
     1142      int numberIntegers = model_->numberIntegers();
     1143      originalLower_=new double [numberIntegers];
     1144      memcpy(originalLower_,rhs.originalLower_,numberIntegers*sizeof(double));
     1145      originalUpper_=new double [numberIntegers];
     1146      memcpy(originalUpper_,rhs.originalUpper_,numberIntegers*sizeof(double));
     1147    } else {
     1148      originalLower_=NULL;
     1149      originalUpper_=NULL;
     1150    }
     1151    delete [] bestSolution_;
     1152    if (rhs.bestSolution_) {
     1153      int numberColumns = model_->getNumCols();
     1154      bestSolution_ = new double [numberColumns];
     1155      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
     1156    } else {
     1157      bestSolution_=NULL;
     1158    }
     1159    delete [] savedSolution_;
     1160    if (rhs.savedSolution_) {
     1161      int numberColumns = model_->getNumCols();
     1162      savedSolution_ = new double [numberColumns];
     1163      memcpy(savedSolution_,rhs.savedSolution_,numberColumns*sizeof(double));
     1164    } else {
     1165      savedSolution_=NULL;
     1166    }
     1167  }
     1168  return *this;
     1169}
     1170// Clone
     1171CbcTree *
     1172CbcTreeVariable::clone() const
     1173{
     1174  return new CbcTreeVariable(*this);
     1175}
     1176// Pass in solution (so can be used after heuristic)
     1177void
     1178CbcTreeVariable::passInSolution(const double * solution, double solutionValue)
     1179{
     1180  int numberColumns = model_->getNumCols();
     1181  delete [] savedSolution_;
     1182  savedSolution_ = new double [numberColumns];
     1183  memcpy(savedSolution_,solution,numberColumns*sizeof(double));
     1184  rhs_=range_;
     1185  // Check feasible
     1186  int goodSolution=createCut(solution,cut_);
     1187  if (goodSolution>=0) {
     1188    bestCutoff_=CoinMin(solutionValue,model_->getCutoff());
     1189  } else {
     1190    model_=NULL;
     1191  }
     1192}
     1193// Return the top node of the heap
     1194CbcNode *
     1195CbcTreeVariable::top() const{
     1196#ifdef CBC_DEBUG
     1197  int smallest=9999999;
     1198  int largest=-1;
     1199  double smallestD=1.0e30;
     1200  double largestD=-1.0e30;
     1201  int n=nodes_.size();
     1202  for (int i=0;i<n;i++) {
     1203    int nn=nodes_[i]->nodeInfo()->nodeNumber();
     1204    double dd = nodes_[i]->objectiveValue();
     1205    largest=CoinMax(largest,nn);
     1206    smallest=CoinMin(smallest,nn);
     1207    largestD=CoinMax(largestD,dd);
     1208    smallestD=CoinMin(smallestD,dd);
     1209  }
     1210  if (model_->messageHandler()->logLevel()>0) {
     1211    printf("smallest %d, largest %d, top %d\n",smallest,largest,
     1212           nodes_.front()->nodeInfo()->nodeNumber());
     1213    printf("smallestD %g, largestD %g, top %g\n",smallestD,largestD,nodes_.front()->objectiveValue());
     1214  }
     1215#endif
     1216  return nodes_.front();
     1217}
     1218
     1219// Add a node to the heap
     1220void
     1221CbcTreeVariable::push(CbcNode * x) {
     1222  if (typeCuts_>=0&&!nodes_.size()&&searchType_<0) {
     1223    startNode_=model_->getNodeCount();
     1224    // save copy of node
     1225    localNode_ = new CbcNode(*x);
     1226
     1227    if (cut_.row().getNumElements()) {
     1228      // Add to global cuts
     1229      // we came in with solution
     1230      model_->globalCuts()->insert(cut_);
     1231      if (model_->messageHandler()->logLevel()>0)
     1232        printf("initial cut - rhs %g %g\n",
     1233               cut_.lb(),cut_.ub());
     1234      searchType_=1;
     1235    } else {
     1236      // stop on first solution
     1237      searchType_=0;
     1238    }
     1239    startTime_ = static_cast<int> (CoinCpuTime());
     1240    saveNumberSolutions_ = model_->getSolutionCount();
     1241  }
     1242  nodes_.push_back(x);
     1243#ifdef CBC_DEBUG
     1244  if (model_->messageHandler()->logLevel()>0)
     1245    printf("pushing node onto heap %d %x %x\n",
     1246           x->nodeInfo()->nodeNumber(),x,x->nodeInfo());
     1247#endif
     1248  std::push_heap(nodes_.begin(), nodes_.end(), comparison_);
     1249}
     1250
     1251// Remove the top node from the heap
     1252void
     1253CbcTreeVariable::pop() {
     1254  std::pop_heap(nodes_.begin(), nodes_.end(), comparison_);
     1255  nodes_.pop_back();
     1256}
     1257// Test if empty - does work if so
     1258bool
     1259CbcTreeVariable::empty()
     1260{
     1261  if (typeCuts_<0)
     1262    return !nodes_.size();
     1263  /* state -
     1264     0 iterating
     1265     1 subtree finished optimal solution for subtree found
     1266     2 subtree finished and no solution found
     1267     3 subtree exiting and solution found
     1268     4 subtree exiting and no solution found
     1269  */
     1270  int state=0;
     1271  assert (searchType_!=2);
     1272  if (searchType_) {
     1273    if (CoinCpuTime()-startTime_>timeLimit_||model_->getNodeCount()-startNode_>=nodeLimit_) {
     1274      state=4;
     1275    }
     1276  } else {
     1277    if (model_->getSolutionCount()>saveNumberSolutions_) {
     1278      state=4;
     1279    }
     1280  }
     1281  if (!nodes_.size())
     1282    state=2;
     1283  if (!state) {
     1284    return false;
     1285  }
     1286  // Finished this phase
     1287  int numberColumns = model_->getNumCols();
     1288  if (model_->getSolutionCount()>saveNumberSolutions_) {
     1289    if(model_->getCutoff()<bestCutoff_) {
     1290      // Save solution
     1291      if (!bestSolution_)
     1292        bestSolution_ = new double [numberColumns];
     1293      memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
     1294      bestCutoff_=model_->getCutoff();
     1295    }
     1296    state--;
     1297  }
     1298  // get rid of all nodes (safe even if already done)
     1299  double bestPossibleObjective;
     1300  cleanTree(model_,-COIN_DBL_MAX,bestPossibleObjective);
     1301
     1302  double increment = model_->getDblParam(CbcModel::CbcCutoffIncrement) ;
     1303  if (model_->messageHandler()->logLevel()>0)
     1304    printf("local state %d after %d nodes and %d seconds, new solution %g, best solution %g, k was %g\n",
     1305           state,
     1306           model_->getNodeCount()-startNode_,
     1307           static_cast<int> (CoinCpuTime())-startTime_,
     1308           model_->getCutoff()+increment,bestCutoff_+increment,rhs_);
     1309  saveNumberSolutions_ = model_->getSolutionCount();
     1310  bool finished=false;
     1311  bool lastTry=false;
     1312  switch (state) {
     1313  case 1:
     1314    // solution found and subtree exhausted
     1315    if (rhs_>1.0e30) {
     1316      finished=true;
     1317    } else {
     1318      // find global cut and reverse
     1319      reverseCut(1);
     1320      searchType_=1; // first false
     1321      rhs_ = range_; // reset range
     1322      nextStrong_=false;
     1323
     1324      // save best solution in this subtree
     1325      memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
     1326    }
     1327    break;
     1328  case 2:
     1329    // solution not found and subtree exhausted
     1330    if (rhs_>1.0e30) {
     1331      finished=true;
     1332    } else {
     1333      // find global cut and reverse
     1334      reverseCut(2);
     1335      searchType_=1; // first false
     1336      if (diversification_<maxDiversification_) {
     1337        if (nextStrong_) {
     1338          diversification_++;
     1339          // cut is valid so don't model_->setCutoff(1.0e50);
     1340          searchType_=0;
     1341        }
     1342        nextStrong_=true;
     1343        rhs_ += range_/2;
     1344      } else {
     1345        // This will be last try (may hit max time)
     1346        lastTry=true;
     1347        if (!maxDiversification_)
     1348          typeCuts_=-1; // make sure can't start again
     1349        model_->setCutoff(bestCutoff_);
     1350        if (model_->messageHandler()->logLevel()>0)
     1351          printf("Exiting local search with current set of cuts\n");
     1352        rhs_=1.0e100;
     1353        // Can now stop on gap
     1354        model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
     1355      }
     1356    }
     1357    break;
     1358  case 3:
     1359    // solution found and subtree not exhausted
     1360    if (rhs_<1.0e30) {
     1361      if (searchType_) {
     1362        if (!typeCuts_&&refine_&&searchType_==1) {
     1363          // We need to check we have best solution given these 0-1 values
     1364          OsiSolverInterface * subSolver = model_->continuousSolver()->clone();
     1365          CbcModel * subModel = model_->subTreeModel(subSolver);
     1366          CbcTree normalTree;
     1367          subModel->passInTreeHandler(normalTree);
     1368          int numberIntegers = model_->numberIntegers();
     1369          const int * integerVariable = model_->integerVariable();
     1370          const double * solution = model_->bestSolution();
     1371          int i;
     1372          int numberColumns = model_->getNumCols();
     1373          for (i=0;i<numberIntegers;i++) {
     1374            int iColumn = integerVariable[i];
     1375            double value = floor(solution[iColumn]+0.5);
     1376            if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
     1377              continue; // skip as not 0-1
     1378            if (originalLower_[i]==originalUpper_[i])
     1379              continue;
     1380            subSolver->setColLower(iColumn,value);
     1381            subSolver->setColUpper(iColumn,value);
     1382          }
     1383          subSolver->initialSolve();
     1384          // We can copy cutoff
     1385          // But adjust
     1386          subModel->setCutoff(model_->getCutoff()+model_->getDblParam(CbcModel::CbcCutoffIncrement)+1.0e-6);
     1387          subModel->setSolutionCount(0);
     1388          assert (subModel->isProvenOptimal());
     1389          if (!subModel->typePresolve()) {
     1390            subModel->branchAndBound();
     1391            if (subModel->status()) {
     1392              model_->incrementSubTreeStopped();
     1393            }
     1394            //printf("%g %g %g %g\n",subModel->getCutoff(),model_->getCutoff(),
     1395            //   subModel->getMinimizationObjValue(),model_->getMinimizationObjValue());
     1396            double newCutoff = subModel->getMinimizationObjValue()-
     1397              subModel->getDblParam(CbcModel::CbcCutoffIncrement) ;
     1398            if (subModel->getSolutionCount()) {
     1399              if (!subModel->status())
     1400                assert (subModel->isProvenOptimal());
     1401              memcpy(model_->bestSolution(),subModel->bestSolution(),
     1402                     numberColumns*sizeof(double));
     1403              model_->setCutoff(newCutoff);
     1404            }
     1405          } else if (subModel->typePresolve()==1) {
     1406            CbcModel * model2 = subModel->integerPresolve(true);
     1407            if (model2) {
     1408              // Do complete search
     1409              model2->branchAndBound();
     1410              // get back solution
     1411              subModel->originalModel(model2,false);
     1412              if (model2->status()) {
     1413                model_->incrementSubTreeStopped();
     1414              }
     1415              double newCutoff = model2->getMinimizationObjValue()-
     1416                model2->getDblParam(CbcModel::CbcCutoffIncrement) ;
     1417              if (model2->getSolutionCount()) {
     1418                if (!model2->status())
     1419                  assert (model2->isProvenOptimal());
     1420                memcpy(model_->bestSolution(),subModel->bestSolution(),
     1421                       numberColumns*sizeof(double));
     1422                model_->setCutoff(newCutoff);
     1423              }
     1424              delete model2;
     1425            } else {
     1426              // infeasible - could just be - due to cutoff
     1427            }
     1428          } else {
     1429            // too dangerous at present
     1430            assert (subModel->typePresolve()!=2);
     1431          }
     1432          if(model_->getCutoff()<bestCutoff_) {
     1433            // Save solution
     1434            if (!bestSolution_)
     1435              bestSolution_ = new double [numberColumns];
     1436            memcpy(bestSolution_,model_->bestSolution(),numberColumns*sizeof(double));
     1437            bestCutoff_=model_->getCutoff();
     1438          }
     1439          delete subModel;
     1440        }
     1441        // we have done search to make sure best general solution
     1442        searchType_=1;
     1443        // Reverse cut weakly
     1444        reverseCut(3,rhs_);
     1445      } else {
     1446        searchType_=1;
     1447        // delete last cut
     1448        deleteCut(cut_);
     1449      }
     1450    } else {
     1451      searchType_=1;
     1452    }
     1453    // save best solution in this subtree
     1454    memcpy(savedSolution_,model_->bestSolution(),numberColumns*sizeof(double));
     1455    nextStrong_=false;
     1456    rhs_=range_;
     1457    break;
     1458  case 4:
     1459    // solution not found and subtree not exhausted
     1460    if (maxDiversification_) {
     1461      if (nextStrong_) {
     1462        // Reverse cut weakly
     1463        reverseCut(4,rhs_);
     1464        model_->setCutoff(1.0e50);
     1465        diversification_++;
     1466        searchType_=0;
     1467      } else {
     1468        // delete last cut
     1469        deleteCut(cut_);
     1470        searchType_=1;
     1471      }
     1472      nextStrong_=true;
     1473      rhs_ += range_/2;
     1474    } else {
     1475      // special case when using as heuristic
     1476      // Reverse cut weakly if lb -infinity
     1477      reverseCut(4,rhs_);
     1478      // This will be last try (may hit max time0
     1479      lastTry=true;
     1480      model_->setCutoff(bestCutoff_);
     1481      if (model_->messageHandler()->logLevel()>0)
     1482        printf("Exiting local search with current set of cuts\n");
     1483      rhs_=1.0e100;
     1484      // Can now stop on gap
     1485      model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
     1486      typeCuts_ =-1;
     1487    }
     1488    break;
     1489  }
     1490  if (rhs_<1.0e30||lastTry) {
     1491    int goodSolution=createCut(savedSolution_,cut_);
     1492    if (goodSolution>=0) {
     1493      // Add to global cuts
     1494      model_->globalCuts()->insert(cut_);
     1495      OsiCuts * global = model_->globalCuts();
     1496      int n = global->sizeRowCuts();
     1497      OsiRowCut * rowCut = global->rowCutPtr(n-1);
     1498      if (model_->messageHandler()->logLevel()>0)
     1499        printf("inserting cut - now %d cuts, rhs %g %g, cutspace %g, diversification %d\n",
     1500               n,rowCut->lb(),rowCut->ub(),rhs_,diversification_);
     1501      const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
     1502      if (debugger) {
     1503        if(debugger->invalidCut(*rowCut))
     1504          printf("ZZZZTree Global cut - cuts off optimal solution!\n");
     1505      }
     1506      for (int i=0;i<n;i++) {
     1507        rowCut = global->rowCutPtr(i);
     1508        if (model_->messageHandler()->logLevel()>0)
     1509          printf("%d - rhs %g %g\n",
     1510                 i,rowCut->lb(),rowCut->ub());
     1511      }
     1512    }
     1513    // put back node
     1514    startTime_ = static_cast<int> (CoinCpuTime());
     1515    startNode_=model_->getNodeCount();
     1516    if (localNode_) {
     1517      // save copy of node
     1518      CbcNode * localNode2 = new CbcNode(*localNode_);
     1519      // But localNode2 now owns cuts so swap
     1520      //printf("pushing local node2 onto heap %d %x %x\n",localNode_->nodeNumber(),
     1521      //   localNode_,localNode_->nodeInfo());
     1522      nodes_.push_back(localNode_);
     1523      localNode_=localNode2;
     1524      std::make_heap(nodes_.begin(), nodes_.end(), comparison_);
     1525    }
     1526  }
     1527  return finished;
     1528}
     1529// We may have got an intelligent tree so give it one more chance
     1530void
     1531CbcTreeVariable::endSearch()
     1532{
     1533  if (typeCuts_>=0) {
     1534    // copy best solution to model
     1535    int numberColumns = model_->getNumCols();
     1536    if (bestSolution_&&bestCutoff_<model_->getCutoff()) {
     1537      memcpy(model_->bestSolution(),bestSolution_,numberColumns*sizeof(double));
     1538      model_->setCutoff(bestCutoff_);
     1539      // recompute objective value
     1540      const double * objCoef = model_->getObjCoefficients();
     1541      double objOffset=0.0;
     1542      model_->continuousSolver()->getDblParam(OsiObjOffset,objOffset);
     1543 
     1544      // Compute dot product of objCoef and colSol and then adjust by offset
     1545      double objValue = -objOffset;
     1546      for ( int i=0 ; i<numberColumns ; i++ )
     1547        objValue += objCoef[i]*bestSolution_[i];
     1548      model_->setMinimizationObjValue(objValue);
     1549    }
     1550    // Can now stop on gap
     1551    model_->setDblParam(CbcModel::CbcAllowableGap,savedGap_);
     1552  }
     1553}
     1554// Create cut
     1555int
     1556CbcTreeVariable::createCut(const double * solution,OsiRowCut & rowCut)
     1557{
     1558  if (rhs_>1.0e20)
     1559    return -1;
     1560  OsiSolverInterface * solver = model_->solver();
     1561  const double * rowLower = solver->getRowLower();
     1562  const double * rowUpper = solver->getRowUpper();
     1563  //const double * solution = solver->getColSolution();
     1564  //const double * objective = solver->getObjCoefficients();
     1565  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
     1566  double primalTolerance;
     1567  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
     1568  // relax
     1569  primalTolerance *= 1000.0;
     1570
     1571  int numberRows = model_->getNumRows();
     1572
     1573  int numberIntegers = model_->numberIntegers();
     1574  const int * integerVariable = model_->integerVariable();
     1575  int i;
     1576
     1577  // Check feasible
     1578
     1579  double * rowActivity = new double[numberRows];
     1580  memset(rowActivity,0,numberRows*sizeof(double));
     1581  solver->getMatrixByCol()->times(solution,rowActivity) ;
     1582  int goodSolution=0;
     1583  // check was feasible
     1584  for (i=0;i<numberRows;i++) {
     1585    if(rowActivity[i]<rowLower[i]-primalTolerance) {
     1586      goodSolution=-1;
     1587    } else if(rowActivity[i]>rowUpper[i]+primalTolerance) {
     1588      goodSolution=-1;
     1589    }
     1590  }
     1591  delete [] rowActivity;
     1592  for (i=0;i<numberIntegers;i++) {
     1593    int iColumn = integerVariable[i];
     1594    double value=solution[iColumn];
     1595    if (fabs(floor(value+0.5)-value)>integerTolerance) {
     1596      goodSolution=-1;
     1597    }
     1598  }
     1599  // zap cut
     1600  if (goodSolution==0) {
     1601    // Create cut and get total gap
     1602    CoinPackedVector cut;
     1603    double rhs = rhs_;
     1604    double maxValue=0.0;
     1605    for (i=0;i<numberIntegers;i++) {
     1606      int iColumn = integerVariable[i];
     1607      double value = floor(solution[iColumn]+0.5);
     1608      if (!typeCuts_&&originalUpper_[i]-originalLower_[i]>1.0)
     1609        continue; // skip as not 0-1
     1610      if (originalLower_[i]==originalUpper_[i])
     1611        continue;
     1612      double mu =1.0/(originalUpper_[i]-originalLower_[i]);
     1613      if (value==originalLower_[i]) {
     1614        rhs += mu*originalLower_[i];
     1615        cut.insert(iColumn,1.0);
     1616        maxValue += originalUpper_[i];
     1617      } else if (value==originalUpper_[i]) {
     1618        rhs -= mu*originalUpper_[i];
     1619        cut.insert(iColumn,-1.0);
     1620        maxValue += originalLower_[i];
     1621      }
     1622    }
     1623    if (maxValue<rhs-primalTolerance) {
     1624      if (model_->messageHandler()->logLevel()>0)
     1625        printf("slack cut\n");
     1626      goodSolution=1;
     1627    }
     1628    rowCut.setRow(cut);
     1629    rowCut.setLb(-COIN_DBL_MAX);
     1630    rowCut.setUb(rhs);   
     1631    rowCut.setGloballyValid();
     1632    if (model_->messageHandler()->logLevel()>0)
     1633    printf("Cut size: %i Cut rhs: %g\n",cut.getNumElements(), rhs);
     1634#ifdef CBC_DEBUG
     1635    if (model_->messageHandler()->logLevel()>0) {
     1636      int k;
     1637      for (k=0; k<cut.getNumElements(); k++){
     1638        printf("%i    %g ",cut.getIndices()[k], cut.getElements()[k]);
     1639        if ((k+1)%5==0)
     1640          printf("\n");
     1641      }
     1642      if (k%5!=0)
     1643        printf("\n");
     1644    }
     1645#endif
     1646    return goodSolution;
     1647  } else {
     1648    if (model_->messageHandler()->logLevel()>0)
     1649      printf("Not a good solution\n");
     1650    return -1;
     1651  }
     1652}
     1653// Other side of last cut branch
     1654void
     1655CbcTreeVariable::reverseCut(int state, double bias)
     1656{
     1657  // find global cut
     1658  OsiCuts * global = model_->globalCuts();
     1659  int n = global->sizeRowCuts();
     1660  int i;
     1661  OsiRowCut * rowCut = NULL;
     1662  for ( i=0;i<n;i++) {
     1663    rowCut = global->rowCutPtr(i);
     1664    if (cut_==*rowCut) {
     1665      break;
     1666    }
     1667  }
     1668  if (!rowCut) {
     1669    // must have got here in odd way e.g. strong branching
     1670    return;
     1671  }
     1672  if (rowCut->lb()>-1.0e10)
     1673    return;
     1674  // get smallest element
     1675  double smallest=COIN_DBL_MAX;
     1676  CoinPackedVector row = cut_.row();
     1677  for (int k=0; k<row.getNumElements(); k++)
     1678    smallest = CoinMin(smallest,fabs(row.getElements()[k]));
     1679  if (!typeCuts_&&!refine_) {
     1680    // Reverse cut very very weakly
     1681    if (state>2)
     1682      smallest=0.0;
     1683  }
     1684  // replace by other way
     1685  if (model_->messageHandler()->logLevel()>0)
     1686    printf("reverseCut - changing cut %d out of %d, old rhs %g %g ",
     1687           i,n,rowCut->lb(),rowCut->ub());
     1688  rowCut->setLb(rowCut->ub()+smallest-bias);
     1689  rowCut->setUb(COIN_DBL_MAX);
     1690  if (model_->messageHandler()->logLevel()>0)
     1691    printf("new rhs %g %g, bias %g smallest %g ",
     1692           rowCut->lb(),rowCut->ub(),bias,smallest);
     1693  const OsiRowCutDebugger *debugger = model_->solver()->getRowCutDebuggerAlways() ;
     1694  if (debugger) {
     1695    if(debugger->invalidCut(*rowCut))
     1696      printf("ZZZZTree Global cut - cuts off optimal solution!\n");
     1697  }
     1698}
     1699// Delete last cut branch
     1700void
     1701CbcTreeVariable::deleteCut(OsiRowCut & cut)
     1702{
     1703  // find global cut
     1704  OsiCuts * global = model_->globalCuts();
     1705  int n = global->sizeRowCuts();
     1706  int i;
     1707  OsiRowCut * rowCut = NULL;
     1708  for ( i=0;i<n;i++) {
     1709    rowCut = global->rowCutPtr(i);
     1710    if (cut==*rowCut) {
     1711      break;
     1712    }
     1713  }
     1714  assert (i<n);
     1715  // delete last cut
     1716  if (model_->messageHandler()->logLevel()>0)
     1717    printf("deleteCut - deleting cut %d out of %d, rhs %g %g\n",
     1718           i,n,rowCut->lb(),rowCut->ub());
     1719  global->eraseRowCut(i);
     1720}
     1721// Create C++ lines to get to current state
     1722void
     1723CbcTreeVariable::generateCpp( FILE * fp)
     1724{
     1725  CbcTreeVariable other;
     1726  fprintf(fp,"0#include \"CbcTreeVariable.hpp\"\n");
     1727  fprintf(fp,"5  CbcTreeVariable variableTree(cbcModel,NULL);\n");
     1728  if (range_!=other.range_)
     1729    fprintf(fp,"5  variableTree.setRange(%d);\n",range_);
     1730  if (typeCuts_!=other.typeCuts_)
     1731    fprintf(fp,"5  variableTree.setTypeCuts(%d);\n",typeCuts_);
     1732  if (maxDiversification_!=other.maxDiversification_)
     1733    fprintf(fp,"5  variableTree.setMaxDiversification(%d);\n",maxDiversification_);
     1734  if (timeLimit_!=other.timeLimit_)
     1735    fprintf(fp,"5  variableTree.setTimeLimit(%d);\n",timeLimit_);
     1736  if (nodeLimit_!=other.nodeLimit_)
     1737    fprintf(fp,"5  variableTree.setNodeLimit(%d);\n",nodeLimit_);
     1738  if (refine_!=other.refine_)
     1739    fprintf(fp,"5  variableTree.setRefine(%s);\n",refine_ ? "true" : "false");
     1740  fprintf(fp,"5  cbcModel->passInTreeHandler(variableTree);\n");
     1741}
     1742
     1743
     1744
  • trunk/Cbc/src/CbcTreeLocal.hpp

    r1173 r1240  
    189189
    190190};
     191
     192class CbcTreeVariable : public CbcTree {
     193
     194public:
     195
     196  // Default Constructor
     197  CbcTreeVariable ();
     198
     199  /* Constructor with solution.
     200     If solution NULL no solution, otherwise must be integer
     201     range is initial upper bound (k) on difference from given solution.
     202     typeCuts -
     203              0 means just 0-1 cuts and will need to refine 0-1 solution
     204              1 uses weaker cuts on all integer variables
     205     maxDiversification is maximum number of range widenings to try
     206     timeLimit is seconds in subTree
     207     nodeLimit is nodes in subTree
     208     refine is whether to see if we can prove current solution is optimal
     209     when we fix all 0-1 (in case typeCuts==0 and there are general integer variables)
     210     if false then no refinement but reverse cuts weaker
     211  */
     212  CbcTreeVariable (CbcModel * model,const double * solution ,int range=10,
     213                   int typeCuts=0,int maxDiversification=0,
     214                   int timeLimit=1000000, int nodeLimit=1000000,bool refine=true);
     215  // Copy constructor
     216  CbcTreeVariable ( const CbcTreeVariable & rhs);
     217
     218  // = operator
     219  CbcTreeVariable & operator=(const CbcTreeVariable & rhs);
     220   
     221  virtual ~CbcTreeVariable();
     222
     223  /// Clone
     224  virtual CbcTree * clone() const;
     225  /// Create C++ lines to get to current state
     226  virtual void generateCpp( FILE * fp) ;
     227
     228/*! \name Heap access and maintenance methods */
     229//@{
     230
     231  /// Return the top node of the heap
     232  virtual CbcNode * top() const;
     233
     234  /// Add a node to the heap
     235  virtual void push(CbcNode * x);
     236
     237  /// Remove the top node from the heap
     238  virtual void pop() ;
     239
     240//@}
     241/*! \name Other stuff */
     242//@{
     243
     244  /// Create cut - return -1 if bad, 0 if okay and 1 if cut is everything
     245  int createCut(const double * solution, OsiRowCut & cut);
     246
     247  /// Test if empty *** note may be overridden
     248  virtual bool empty() ;
     249
     250  /// We may have got an intelligent tree so give it one more chance
     251  virtual void endSearch() ;
     252  /// Other side of last cut branch (if bias==rhs_ will be weakest possible)
     253  void reverseCut(int state, double bias=0.0);
     254  /// Delete last cut branch
     255  void deleteCut(OsiRowCut & cut);
     256  /// Pass in solution (so can be used after heuristic)
     257  void passInSolution(const double * solution, double solutionValue);
     258  // range i.e. k
     259  inline int range() const
     260  { return range_;}
     261  // setrange i.e. k
     262  inline void setRange(int value)
     263  { range_ = value;}
     264  // Type of cuts - 0=just 0-1, 1=all
     265  inline int typeCuts() const
     266  { return typeCuts_;}
     267  // Type of cuts - 0=just 0-1, 1=all
     268  inline void setTypeCuts(int value)
     269  { typeCuts_ = value;}
     270  // maximum number of diversifications
     271  inline int maxDiversification() const
     272  { return maxDiversification_;}
     273  // maximum number of diversifications
     274  inline void setMaxDiversification(int value)
     275  { maxDiversification_ = value;}
     276  // time limit per subtree
     277  inline int timeLimit() const
     278  { return timeLimit_;}
     279  // time limit per subtree
     280  inline void setTimeLimit(int value)
     281  { timeLimit_ = value;}
     282  // node limit for subtree
     283  inline int nodeLimit() const
     284  { return nodeLimit_;}
     285  // node limit for subtree
     286  inline void setNodeLimit(int value)
     287  { nodeLimit_ = value;}
     288  // Whether to do refinement step
     289  inline bool refine() const
     290  { return refine_;}
     291  // Whether to do refinement step
     292  inline void setRefine(bool yesNo)
     293    { refine_ = yesNo;}
     294
     295//@}
     296private:
     297  // Node for local cuts
     298  CbcNode * localNode_;
     299  // best solution
     300  double * bestSolution_;
     301  // saved solution
     302  double * savedSolution_;
     303  // solution number at start of pass
     304  int saveNumberSolutions_;
     305  /* Cut.  If zero size then no solution yet.  Otherwise is left hand branch */
     306  OsiRowCut cut_;
     307  // This cut fixes all 0-1 variables
     308  OsiRowCut fixedCut_;
     309  // Model
     310  CbcModel * model_;
     311  // Original lower bounds
     312  double * originalLower_;
     313  // Original upper bounds
     314  double * originalUpper_;
     315  // range i.e. k
     316  int range_;
     317  // Type of cuts - 0=just 0-1, 1=all
     318  int typeCuts_;
     319  // maximum number of diversifications
     320  int maxDiversification_;
     321  // current diversification
     322  int diversification_;
     323  // Whether next will be strong diversification
     324  bool nextStrong_;
     325  // Current rhs
     326  double rhs_;
     327  // Save allowable gap
     328  double savedGap_;
     329  // Best solution
     330  double bestCutoff_;
     331  // time limit per subtree
     332  int timeLimit_;
     333  // time when subtree started
     334  int startTime_;
     335  // node limit for subtree
     336  int nodeLimit_;
     337  // node count when subtree started
     338  int startNode_;
     339  // -1 not started, 0 == stop on first solution, 1 don't stop on first, 2 refinement step
     340  int searchType_;
     341  // Whether to do refinement step
     342  bool refine_;
     343
     344};
    191345#endif
    192346
Note: See TracChangeset for help on using the changeset viewer.