Changeset 208


Ignore:
Timestamp:
Nov 11, 2005 3:00:06 PM (14 years ago)
Author:
forrest
Message:

large number of changes

Location:
trunk
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • trunk/CbcBranchDynamic.cpp

    r202 r208  
    274274    preferredWay=-1;
    275275  // weight at 1.0 is max min
    276 #define WEIGHT_AFTER 0.8
     276#define WEIGHT_AFTER 0.9
    277277#define WEIGHT_BEFORE 0.3
    278278  if (fabs(value-nearest)<=integerTolerance) {
     
    388388#endif
    389389  } else {
    390     OsiSolverInterface * solver = model_->solver();
    391390    const double * upper = model_->getCbcColUpper();
    392391    double integerTolerance =
     
    398397      below = above -1;
    399398    }
    400     double objectiveValue = solver->getObjSense()*solver->getObjValue();
     399    double objectiveValue = model_->getCurrentMinimizationObjValue();
    401400    double distanceToCutoff =  model_->getCutoff() - objectiveValue;
    402401    if (distanceToCutoff<1.0e20)
     
    682681  int betterWay=0;
    683682  double value=0.0;
    684   if (stateOfSearch<=1||thisOne->model()->currentNode()->depth()<=10) {
     683  if (stateOfSearch<=1&&thisOne->model()->currentNode()->depth()>10) {
    685684#if 0
    686685    if (!bestObject_) {
  • trunk/CbcCompareActual.cpp

    r202 r208  
    192192  }
    193193#else
    194   if (weight_==-1.0) {
     194  if (weight_==-1.0&&(y->depth()>7||x->depth()>7)) {
    195195    // before solution
    196196    /* printf("x %d %d %g, y %d %d %g\n",
  • trunk/CbcCutGenerator.cpp

    r202 r208  
    159159  int pass=model_->getCurrentPassNumber()-1;
    160160  bool doThis=(model_->getNodeCount()%howOften)==0;
    161   if (depthCutGenerator_>0)
     161  if (depthCutGenerator_>0) {
    162162    doThis = (depth % depthCutGenerator_) ==0;
     163    if (depth<depthCutGenerator_)
     164      doThis=true; // and also at top of tree
     165  }
    163166  // But turn off if 100
    164167  if (howOften==100)
     
    216219    // switch off if first time and no good
    217220    if (node==NULL&&!pass) {
    218       if (cs.sizeCuts()-cutsBefore<switchOffIfLessThan_) {
     221      if (cs.sizeCuts()-cutsBefore<CoinAbs(switchOffIfLessThan_)) {
    219222        whenCutGenerator_=-99;
    220223        whenCutGeneratorInSub_ = -200;
  • trunk/CbcHeuristic.cpp

    r206 r208  
    103103          << model.getMaximumNodes()
    104104          <<CoinMessageEol;
     105      // probably faster to use a basis to get integer solutions
     106      model.setSpecialOptions(2);
    105107      model.branchAndBound();
    106108      if (logLevel>1)
  • trunk/CbcMessage.cpp

    r197 r208  
    4848  {CBC_END_SUB,29,1,"Ending sub-tree for %s"},
    4949  {CBC_HEURISTIC_SOLUTION,30,1,"solution of %g found by %s"},
     50  {CBC_CUTS_STATS,31,1,"%d added rows had average density of %g"},
    5051  {CBC_DUMMY_END,999999,0,""}
    5152};
  • trunk/CbcModel.cpp

    r202 r208  
    1212#include <cmath>
    1313#include <cfloat>
     14#define COIN_USE_CLP
     15#ifdef COIN_USE_CLP
     16// include Presolve from Clp
     17#include "ClpPresolve.hpp"
     18#include "OsiClpSolverInterface.hpp"
     19#include "ClpEventHandler.hpp"
     20#endif
    1421
    1522#include "OsiSolverInterface.hpp"
     
    3542#include "CglProbing.hpp"
    3643
    37 #define COIN_USE_CLP
    38 #ifdef COIN_USE_CLP
    39 // include Presolve from Clp
    40 #include "ClpPresolve.hpp"
    41 #include "OsiClpSolverInterface.hpp"
    42 #include "ClpEventHandler.hpp"
    43 #endif
    4444
    4545#include "CoinTime.hpp"
     
    532532  continuousSolver_ = solver_->clone() ;
    533533#ifdef COIN_USE_CLP
    534   {
    535     OsiClpSolverInterface * clpSolver
    536       = dynamic_cast<OsiClpSolverInterface *> (solver_);
    537     if (clpSolver) {
    538       ClpSimplex * clpSimplex = clpSolver->getModelPtr();
    539       // take off names
    540       clpSimplex->dropNames();
    541     }
     534  OsiClpSolverInterface * clpSolver
     535    = dynamic_cast<OsiClpSolverInterface *> (solver_);
     536  if (clpSolver) {
     537    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     538    // take off names
     539    clpSimplex->dropNames();
    542540  }
    543541#endif
     
    664662  bool resolved = false ;
    665663  CbcNode *newNode = NULL ;
     664  int numberFixedAtRoot=0;
     665  int numberFixedNow=0;
    666666  if (feasible)
    667667  { newNode = new CbcNode ;
     
    675675    }
    676676    phase_=3;
    677     // only allow twenty passes
    678     int numberPassesLeft=20;
     677    // only allow 1000 passes
     678    int numberPassesLeft=1000;
     679    // This is first crude step
     680    if (numberAnalyzeIterations_) {
     681      delete [] analyzeResults_;
     682      analyzeResults_ = new double [4*numberIntegers_];
     683      numberFixedAtRoot=newNode->analyze(this,analyzeResults_);
     684      if (numberFixedAtRoot>0) {
     685        printf("%d fixed by analysis\n",numberFixedAtRoot);
     686        setPointers(solver_);
     687        numberFixedNow = numberFixedAtRoot;
     688      } else if (numberFixedAtRoot<0) {
     689        printf("analysis found to be infeasible\n");
     690        anyAction=-2;
     691        delete newNode ;
     692        newNode = NULL ;
     693        feasible = false ;
     694      }
     695    }
    679696    while (anyAction == -1)
    680697    {
     
    692709        feasible=false; // pretend infeasible
    693710      }
     711      //int nfix=
     712      reducedCostFix() ;
     713      //if (nfix)
     714      //printf("%d fixed after resolve root\n",nfix);
    694715        resolved = true ;
    695716#       ifdef CBC_DEBUG
     
    831852  than the current objective cutoff.
    832853*/
    833   while (!tree_->empty())
    834   { if (cutoff > getCutoff()) {
    835     if (eventHandler) {
    836       if (!eventHandler->event(ClpEventHandler::solution)) {
    837         eventHappened=true; // exit
    838       }
    839     }
     854  while (!tree_->empty()) {
     855    if (cutoff > getCutoff()) {
     856      double newCutoff = getCutoff();
     857      if (analyzeResults_) {
     858        // see if we could fix any (more)
     859        int n=0;
     860        double * newLower = analyzeResults_;
     861        double * objLower = newLower+numberIntegers_;
     862        double * newUpper = objLower+numberIntegers_;
     863        double * objUpper = newUpper+numberIntegers_;
     864        for (int i=0;i<numberIntegers_;i++) {
     865          if (objLower[i]>newCutoff) {
     866            n++;
     867            if (objUpper[i]>newCutoff) {
     868              newCutoff = -COIN_DBL_MAX;
     869              break;
     870            }
     871          } else if (objUpper[i]>newCutoff) {
     872            n++;
     873          }
     874        }
     875        if (newCutoff==-COIN_DBL_MAX) {
     876          printf("Root analysis says finished\n");
     877        } else if (n>numberFixedNow) {
     878          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow,n);
     879          numberFixedNow=n;
     880        }
     881      }
     882      if (eventHandler) {
     883        if (!eventHandler->event(ClpEventHandler::solution)) {
     884          eventHappened=true; // exit
     885        }
     886      }
    840887      // Do from deepest
    841       tree_->cleanTree(this, getCutoff(),bestPossibleObjective_) ;
     888      tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ;
    842889      nodeCompare_->newSolution(this) ;
    843890      nodeCompare_->newSolution(this,continuousObjective_,
     
    946993    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
    947994    newNode = NULL ;
    948     if (!addCuts(node,lastws))
     995    if (!addCuts(node,lastws,numberFixedNow>numberFixedAtRoot))
    949996    { int i ;
    950997      const double * lower = getColLower() ;
     
    10441091        { newNode = new CbcNode ;
    10451092          newNode->setObjectiveValue(direction*solver_->getObjValue()) ;
    1046           if (newNode->objectiveValue() >= getCutoff())
    1047             anyAction=-2;
    10481093          anyAction =-1 ;
    10491094          resolved = false ;
     
    10841129              { newNode->setObjectiveValue(direction*
    10851130                                           solver_->getObjValue()) ;
     1131              //int nfix=
     1132              reducedCostFix() ;
     1133              //if (nfix)
     1134              //printf("%d fixed after resolve\n",nfix);
    10861135              if (newNode->objectiveValue() >= getCutoff())
    10871136                anyAction=-2;
     
    16491698  numberPenalties_(20),
    16501699  penaltyScaleFactor_(3.0),
     1700  numberAnalyzeIterations_(0),
     1701  analyzeResults_(NULL),
    16511702  numberInfeasibleNodes_(0),
    16521703  problemType_(0),
     
    17481799  numberPenalties_(20),
    17491800  penaltyScaleFactor_(3.0),
     1801  numberAnalyzeIterations_(0),
     1802  analyzeResults_(NULL),
    17501803  numberInfeasibleNodes_(0),
    17511804  problemType_(0),
     
    19331986  numberPenalties_(rhs.numberPenalties_),
    19341987  penaltyScaleFactor_(rhs.penaltyScaleFactor_),
     1988  numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
     1989  analyzeResults_(NULL),
    19351990  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
    19361991  problemType_(rhs.problemType_),
     
    21762231    numberBeforeTrust_ = rhs.numberBeforeTrust_;
    21772232    numberPenalties_ = rhs.numberPenalties_;
    2178     penaltyScaleFactor_ = penaltyScaleFactor_;
     2233    penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
     2234    numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
     2235    delete [] analyzeResults_;
     2236    analyzeResults_ = NULL;
    21792237    numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
    21802238    problemType_ = rhs.problemType_;
     
    26172675  and pass it back out if required.
    26182676*/
    2619 int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws)
     2677int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws,bool canFix)
    26202678{
    26212679/*
     
    26292687  double cutoff = getCutoff() ;
    26302688  int currentNumberCuts=currentNumberCuts_;
     2689  if (canFix) {
     2690    bool feasible=true;
     2691    const double *lower = solver_->getColLower() ;
     2692    const double *upper = solver_->getColUpper() ;
     2693    double * newLower = analyzeResults_;
     2694    double * objLower = newLower+numberIntegers_;
     2695    double * newUpper = objLower+numberIntegers_;
     2696    double * objUpper = newUpper+numberIntegers_;
     2697    int n=0;
     2698    for (i=0;i<numberIntegers_;i++) {
     2699      int iColumn = integerVariable_[i];
     2700      bool changed=false;
     2701      double lo = 0.0;
     2702      double up = 0.0;
     2703      if (objLower[i]>cutoff) {
     2704        lo = lower[iColumn];
     2705        up = upper[iColumn];
     2706        if (lo<newLower[i]) {
     2707          lo = newLower[i];
     2708          solver_->setColLower(iColumn,lo);
     2709          changed=true;
     2710          n++;
     2711        }
     2712        if (objUpper[i]>cutoff) {
     2713          if (up>newUpper[i]) {
     2714            up = newUpper[i];
     2715            solver_->setColUpper(iColumn,up);
     2716            changed=true;
     2717            n++;
     2718          }
     2719        }
     2720      } else if (objUpper[i]>cutoff) {
     2721        lo = lower[iColumn];
     2722        up = upper[iColumn];
     2723        if (up>newUpper[i]) {
     2724          up = newUpper[i];
     2725          solver_->setColUpper(iColumn,up);
     2726          changed=true;
     2727          n++;
     2728        }
     2729      }
     2730      if (changed&&lo>up) {
     2731        feasible=false;
     2732        break;
     2733      }
     2734    }
     2735    if (!feasible) {
     2736      printf("analysis says node infeas\n");
     2737      cutoff=-COIN_DBL_MAX;
     2738    }
     2739  }
    26312740/*
    26322741  If the node can't be fathomed by bound, reinstall tight cuts in the
     
    27402849*/
    27412850
    2742 void CbcModel::reducedCostFix ()
     2851int CbcModel::reducedCostFix ()
    27432852
    27442853{ double cutoff = getCutoff() ;
    27452854  double direction = solver_->getObjSense() ;
    27462855  double gap = cutoff - solver_->getObjValue()*direction ;
     2856  double tolerance;
     2857  solver_->getDblParam(OsiDualTolerance,tolerance) ;
     2858  gap += tolerance;
     2859  if (gap<=0.0)
     2860    return 0;
    27472861  double integerTolerance = getDblParam(CbcIntegerTolerance) ;
    27482862
     
    27532867
    27542868  int numberFixed = 0 ;
     2869#ifdef COIN_USE_CLP
     2870  OsiClpSolverInterface * clpSolver
     2871    = dynamic_cast<OsiClpSolverInterface *> (solver_);
     2872  ClpSimplex * clpSimplex=NULL;
     2873  if (clpSolver)
     2874    clpSimplex = clpSolver->getModelPtr();
     2875#endif
    27552876  for (int i = 0 ; i < numberIntegers_ ; i++)
    27562877  { int iColumn = integerVariable_[i] ;
     
    27592880    { if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap)
    27602881      { solver_->setColUpper(iColumn,lower[iColumn]) ;
     2882#ifdef COIN_USE_CLP
     2883      if (clpSimplex)
     2884        assert (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound);
     2885#endif
    27612886        numberFixed++ ; }
    27622887      else
    27632888      if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap)
    27642889      { solver_->setColLower(iColumn,upper[iColumn]) ;
     2890#ifdef COIN_USE_CLP
     2891      if (clpSimplex)
     2892        assert (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound);
     2893#endif
    27652894        numberFixed++ ; } } }
    27662895 
    2767   return ; }
     2896  return numberFixed; }
    27682897
    27692898// Collect coding to replace whichGenerator
     
    28232952  int numberRowsAtStart = solver_->getNumRows() ;
    28242953  int numberColumns = solver_->getNumCols() ;
     2954  CoinBigIndex numberElementsAtStart = solver_->getNumElements();
    28252955
    28262956  numberOldActiveCuts = numberRowsAtStart-numberRowsAtContinuous_ ;
     
    30963226        }
    30973227      }
    3098       // Add in any violated saved cuts
    3099       if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
    3100         int numberCuts = slackCuts.sizeRowCuts() ;
    3101         int i;
    3102         // possibly extend whichGenerator
    3103         whichGenerator = newWhichGenerator(numberViolated, numberViolated+numberCuts,
    3104                                            maximumWhich,  whichGenerator);
    3105         for ( i = 0;i<numberCuts;i++) {
    3106           const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
    3107           if (thisCut->violated(solution)>100.0*primalTolerance) {
    3108             if (messageHandler()->logLevel()>2)
    3109               printf("Old cut added - violation %g\n",
    3110                      thisCut->violated(solution)) ;
    3111             whichGenerator[numberViolated++]=-1;
    3112             theseCuts.insert(*thisCut) ;
    3113           }
    3114         }
    3115       }
    31163228      int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
    31173229      int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
     
    31963308      numberSolutions_ = saveNumberSolutions;
    31973309      numberHeuristicSolutions_ = saveNumberHeuristicSolutions;
     3310    }
     3311    // Add in any violated saved cuts
     3312    if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
     3313      int numberOld = theseCuts.sizeRowCuts();
     3314      int numberCuts = slackCuts.sizeRowCuts() ;
     3315      int i;
     3316      // possibly extend whichGenerator
     3317      whichGenerator = newWhichGenerator(numberOld, numberOld+numberCuts,
     3318                                         maximumWhich,  whichGenerator);
     3319      for ( i = 0;i<numberCuts;i++) {
     3320        const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
     3321        if (thisCut->violated(solution)>100.0*primalTolerance) {
     3322          if (messageHandler()->logLevel()>2)
     3323            printf("Old cut added - violation %g\n",
     3324                   thisCut->violated(solution)) ;
     3325          whichGenerator[numberOld++]=-1;
     3326          theseCuts.insert(*thisCut) ;
     3327        }
     3328      }
    31983329    }
    31993330/*
     
    34673598  } while (numberTries>0) ;
    34683599  // Reduced cost fix at end
    3469   //reducedCostFix();
     3600  if (solver_->isProvenOptimal())
     3601    reducedCostFix();
    34703602  // If at root node do heuristics
    34713603  if (!numberNodes_) {
     3604    // First see if any cuts are slack
     3605    int numberRows = solver_->getNumRows();
     3606    int numberAdded = numberRows-numberRowsAtContinuous_;
     3607    if (numberAdded) {
     3608      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
     3609      assert(basis != NULL);
     3610      int * added = new int[numberAdded];
     3611      int nDelete = 0;
     3612      for (int j=numberRowsAtContinuous_;j<numberRows;j++) {
     3613        if (basis->getArtifStatus(j)==CoinWarmStartBasis::basic) {
     3614          //printf("%d slack!\n",j);
     3615          added[nDelete++]=j;
     3616        }
     3617      }
     3618      if (nDelete) {
     3619        solver_->deleteRows(nDelete,added);
     3620        delete basis_ ;
     3621        basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
     3622      }
     3623      delete [] added;
     3624      delete basis ;
     3625    }
    34723626    // mark so heuristics can tell
    34733627    int savePass=currentPassNumber_;
     
    35273681    int willBeCutsInTree=0;
    35283682    double thisObjective = solver_->getObjValue()*direction ;
    3529     if (thisObjective-startObjective<1.0e-5)
    3530       willBeCutsInTree=-1;
     3683    // get sizes
     3684    int numberRowsAdded = solver_->getNumRows()-numberRowsAtStart;
     3685    CoinBigIndex numberElementsAdded = solver_->getNumElements()-numberElementsAtStart;
     3686    double densityOld = ((double) numberElementsAtStart)/((double) numberRowsAtStart);
     3687    double densityNew = numberRowsAdded ? ((double) (numberElementsAdded))/((double) numberRowsAdded)
     3688      : 0.0;
     3689    if (!numberNodes_) {
     3690      handler_->message(CBC_CUTS_STATS,messages_)
     3691        <<numberRowsAdded
     3692        <<densityNew
     3693        <<CoinMessageEol ;
     3694      if (thisObjective-startObjective<1.0e-5&&numberElementsAdded>0.2*numberElementsAtStart)
     3695        willBeCutsInTree=-1;
     3696    }
     3697    if ((numberRowsAdded>100+0.5*numberRowsAtStart
     3698         ||numberElementsAdded>0.5*numberElementsAtStart)
     3699        &&(densityNew>200.0&&numberRowsAdded>100&&densityNew>2.0*densityOld)) {
     3700      // much bigger
     3701      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
     3702      //willBeCutsInTree=-1;
     3703      //printf("Cuts will be taken off , %d rows added with density %g\n",
     3704      //     numberRowsAdded,densityNew);
     3705    }
     3706    if (densityNew>100.0&&numberRowsAdded>2&&densityNew>2.0*densityOld) {
     3707      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
     3708      //willBeCutsInTree=-2;
     3709      //printf("Density says no cuts ? , %d rows added with density %g\n",
     3710      //     numberRowsAdded,densityNew);
     3711    }
    35313712    // Root node or every so often - see what to turn off
    35323713    int i ;
     3714    for (i = 0;i<numberCutGenerators_;i++) {
     3715      int howOften = generator_[i]->howOften() ;
     3716      if (howOften>-90)
     3717        willBeCutsInTree=0;
     3718    }
    35333719    double totalCuts = 0.0 ;
    35343720    for (i = 0;i<numberCutGenerators_;i++)
     
    35443730    for (i = 0;i<numberNewCuts;i++) {
    35453731      int iGenerator = whichGenerator[i];
    3546       if (iGenerator>=0)
     3732      if (iGenerator>=0&&iGenerator<numberCutGenerators_)
    35473733        count[iGenerator]++ ;
    35483734    }
     
    35513737    for (i = 0;i<numberCutGenerators_;i++) {
    35523738      int howOften = generator_[i]->howOften() ;
    3553       if (howOften==-98) {
     3739      if (willBeCutsInTree<0&&howOften==-98)
     3740        howOften =-99;
     3741      if (howOften==-98&&generator_[i]->switchOffIfLessThan()>=0) {
    35543742        if (thisObjective-startObjective<0.005*fabs(startObjective)+1.0e-5)
    35553743          howOften=-99; // switch off
    3556         if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5&&solver_->getNumRows()<500)
     3744        if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5
     3745            &&5*solver_->getNumRows()<solver_->getNumCols())
    35573746          howOften=-99; // switch off
    35583747      }
     
    35833772        if (numberNodes_>=100000&&sumChangeObjective1_>2.0e2*(sumChangeObjective2_+1.0e-12)) {
    35843773          howOften = 1000000+SCANCUTS; // wait until next time
    3585           printf("switch off cut %d due to lack of use\n",i);
     3774          //printf("switch off cut %d due to lack of use\n",i);
    35863775        }
    35873776      }
     
    36193808    delete [] count ;
    36203809    if( !numberNodes_) {
     3810      if (willBeCutsInTree==-2)
     3811        willBeCutsInTree=0;
    36213812      if( willBeCutsInTree<=0) {
    36223813        // Take off cuts
     
    45904781      delete slack ;
    45914782    }
    4592     // Give a hint not to do scaling
    4593     //bool saveTakeHint;
    4594     //OsiHintStrength saveStrength;
    4595     //bool gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
    4596     //assert (gotHint);
    4597     //solver_->setHintParam(OsiDoScale,false,OsiHintTry);
     4783    // Give a hint to do dual
     4784    bool saveTakeHint;
     4785    OsiHintStrength saveStrength;
     4786    bool gotHint = (solver_->getHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength));
     4787    assert (gotHint);
     4788    solver_->setHintParam(OsiDoDualInInitial,true,OsiHintTry);
    45984789    solver_->initialSolve();
    4599     //solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
    46004790    if (!solver_->isProvenOptimal())
    46014791      { //printf("checkSolution infeas! Retrying wihout scaling.\n");
     
    46034793      OsiHintStrength saveStrength;
    46044794      bool savePrintHint;
    4605       solver_->writeMps("infeas");
     4795      //solver_->writeMps("infeas");
    46064796      bool gotHint = (solver_->getHintParam(OsiDoReducePrint,savePrintHint,saveStrength));
    46074797      gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
    46084798      solver_->setHintParam(OsiDoScale,false,OsiHintTry);
    4609       solver_->setHintParam(OsiDoReducePrint,false,OsiHintTry) ;
     4799      //solver_->setHintParam(OsiDoReducePrint,false,OsiHintTry) ;
     4800      solver_->setHintParam(OsiDoDualInInitial,false,OsiHintTry);
    46104801      solver_->initialSolve();
    46114802      solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
    4612       solver_->setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ;
     4803      //solver_->setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ;
    46134804      }
    46144805    //assert(solver_->isProvenOptimal());
     4806    solver_->setHintParam(OsiDoDualInInitial,saveTakeHint,saveStrength);
    46154807  }
    46164808  double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
     
    59516143{
    59526144/*
    5953   Assume we're done, and see if we're proven wrong.
    5954 */
     6145  Switch off heuristics
     6146*/
     6147  int saveNumberHeuristics=numberHeuristics_;
     6148  numberHeuristics_=0;
    59556149/*
    59566150  Scan the variables, noting the integer variables. Create an
     
    60906284*/
    60916285  globalCuts_= OsiCuts() ;
     6286  numberHeuristics_ = saveNumberHeuristics;
    60926287 
    60936288  return newSolver;
     
    62466441    dblParam_[CbcOptimizationDirection];
    62476442}
     6443#ifdef COIN_USE_CLP
     6444// Pass in Event handler (cloned and deleted at end)
     6445void
     6446CbcModel::passInEventHandler(const ClpEventHandler * eventHandler)
     6447{
     6448  OsiClpSolverInterface * clpSolver
     6449    = dynamic_cast<OsiClpSolverInterface *> (solver_);
     6450  if (clpSolver) {
     6451    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     6452    clpSimplex->passInEventHandler(eventHandler);
     6453  }
     6454}
     6455// Event handler
     6456ClpEventHandler *
     6457CbcModel::eventHandler() const
     6458{
     6459  OsiClpSolverInterface * clpSolver
     6460    = dynamic_cast<OsiClpSolverInterface *> (solver_);
     6461  if (clpSolver) {
     6462    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
     6463    return clpSimplex->eventHandler();
     6464  } else {
     6465    return NULL;
     6466  }
     6467}
     6468#endif
  • trunk/CbcNode.cpp

    r202 r208  
    16871687        if (whichObject>=0) {
    16881688          branch_ = objects[whichObject];
     1689          if (model->messageHandler()->logLevel()>3)
     1690            printf("Choosing column %d\n",choice[whichObject].possibleBranch->variable()) ;
    16891691          choice[whichObject].possibleBranch=NULL;
    16901692        }
     
    17761778  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
    17771779  double cutoff =model->getCutoff();
    1778   double distanceToCutoff = cutoff-objectiveValue_;
    17791780  const double * lower = solver->getColLower();
    17801781  const double * upper = solver->getColUpper();
     
    17861787  }
    17871788  int i;
    1788   int stateOfSearch = model->stateOfSearch();
     1789  int saveStateOfSearch = model->stateOfSearch();
    17891790  int numberStrong=model->numberStrong();
    17901791  // But make more likely to get out after some times
     
    18081809  double * saveUpper = new double[numberColumns];
    18091810  double * saveLower = new double[numberColumns];
    1810   // Ratio to cutoff for pseudo costs
    1811   double tryStrongPseudo = 0.8*distanceToCutoff;
    1812   // Ratio to cutoff for penalties
    1813   //double tryStrongPenalty = 0.5*distanceToCutoff;
    18141811
    18151812  // Save solution in case heuristics need good solution later
     
    18491846    numberBeforeTrust=0;
    18501847  }
    1851   double scaleFactor = model->penaltyScaleFactor();
    18521848  // May go round twice if strong branching fixes all local candidates
    18531849  bool finished=false;
     1850  int numberToFix=0;
     1851  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
     1852  int saveClpOptions=0;
     1853  if (osiclp) {
     1854    // for faster hot start
     1855    saveClpOptions = osiclp->specialOptions();
     1856    osiclp->setSpecialOptions(saveClpOptions|1024);
     1857  }
    18541858  while(!finished) {
    18551859    finished=true;
     
    18571861    // Some objects may compute an estimate of best solution from here
    18581862    estimatedDegradation=0.0;
    1859     int numberToFix=0;
     1863    numberToFix=0;
    18601864    int numberIntegerInfeasibilities=0; // without odd ones
    18611865    int numberToDo=0;
     
    18681872    double best=0.0;
    18691873    int numberNotTrusted=0;
    1870    
     1874    int * which = objectMark+numberObjects+1;
     1875    int neededPenalties;
    18711876    // We may go round this loop twice (only if we think we have solution)
    18721877    for (int iPass=0;iPass<2;iPass++) {
     
    18931898      */
    18941899      numberToDo=0;
     1900      neededPenalties=0;
    18951901      averageDown=0.0;
    18961902      numberDown=0;
     
    19101916        double infeasibility = object->infeasibility(preferredWay);
    19111917        int priorityLevel = object->priority();
    1912         bool gotDown=false;
    1913         int numberThisDown = dynamicObject->numberTimesDown();
    1914         if (numberThisDown) {
    1915           averageDown += dynamicObject->downDynamicPseudoCost();
    1916           numberDown++;
    1917           if (numberThisDown>=numberBeforeTrust)
    1918             gotDown=true;
    1919         }
    1920         bool gotUp=false;
    1921         int numberThisUp = dynamicObject->numberTimesUp();
    1922         if (numberThisUp) {
    1923           averageUp += dynamicObject->upDynamicPseudoCost();
    1924           numberUp++;
    1925           if (numberThisUp>=numberBeforeTrust)
    1926             gotUp=true;
    1927         }
    19281918        if (infeasibility) {
     1919          bool gotDown=false;
     1920          int numberThisDown = dynamicObject->numberTimesDown();
     1921          if (numberThisDown) {
     1922            averageDown += dynamicObject->downDynamicPseudoCost();
     1923            numberDown++;
     1924            if (numberThisDown>=numberBeforeTrust)
     1925              gotDown=true;
     1926          }
     1927          bool gotUp=false;
     1928          int numberThisUp = dynamicObject->numberTimesUp();
     1929          if (numberThisUp) {
     1930            averageUp += dynamicObject->upDynamicPseudoCost();
     1931            numberUp++;
     1932            if (numberThisUp>=numberBeforeTrust)
     1933              gotUp=true;
     1934          }
    19291935          if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
    19301936            printf("%d down %d %g up %d %g - infeas %g\n",
     
    19551961              iBestGot=numberToDo;
    19561962            }
    1957             if (infeasibility<tryStrongPseudo)
    1958               objectMark[i]=2*numberBeforeTrust; // range info not wanted
    1959             else
    1960               objectMark[i]=numberBeforeTrust; // keep as possible cutoff
    19611963          } else {
     1964            objectMark[neededPenalties]=numberToDo;
     1965            which[neededPenalties++]=dynamicObject->columnNumber();
    19621966            int iColumn = dynamicObject->columnNumber();
    19631967            double part =saveSolution[iColumn]-floor(saveSolution[iColumn]);
    1964             sort[numberToDo]=part;
     1968            sort[numberToDo]=-infeasibility;
    19651969            if (1.0-fabs(part-0.5)>bestNot) {
    19661970              iBestNot=numberToDo;
    19671971              bestNot = 1.0-fabs(part-0.5);
    19681972            }
    1969             objectMark[i]=numberThisDown+numberThisUp;
    19701973          }
    19711974          whichObject[numberToDo++]=i;
     
    20432046      ws = solver->getWarmStart();
    20442047      solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
    2045       if (!stateOfSearch&&saveLimit<100)
     2048      if (!saveStateOfSearch&&saveLimit<100)
    20462049        solver->setIntParam(OsiMaxNumIterationHotStart,100);
    20472050    }
     
    20732076      if (!skipAll)
    20742077        solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
    2075       if (numberDown)
    2076         averageDown /= (double) numberDown;
    2077       else
    2078         averageDown=1.0;
    2079       if (numberUp)
    2080         averageUp /= (double) numberUp;
    2081       else
    2082         averageUp=1.0;
    2083       double average = 0.5*(averageUp+averageDown);
    20842078      int iDo;
    2085       // Bias by trust
    2086       int iObj;
    2087       for (iObj=0;iObj<numberToDo;iObj++) {
    2088         int iObject=whichObject[iObj];
    2089         double add = objectMark[iObject];
    2090         if (sort[iObj]<0.0)
    2091           sort[iObj] += add*average;
    2092         else
    2093           sort[iObj] = -(sort[iObj]*averageDown+(1.0-sort[iObj])*averageUp);
    2094       }
    2095       // Sort
    2096       CoinSort_2(sort,sort+numberToDo,whichObject);
    2097       int needed2=numberToDo;
    20982079#define RANGING
    20992080#ifdef RANGING
    2100       OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
    21012081      if (skipAll&&numberBeforeTrust)
    21022082        numberPenalties=0;
    2103       if (osiclp&&numberPenalties) {
    2104         // just get those not touched and best and not trusted
    2105         int n=CoinMin(numberPenalties,numberToDo);
    2106         int * which = objectMark+numberObjects+1;
    2107         int i;
    2108         int needed=0;
    2109         for ( i=0;i<n;i++) {
    2110           int iObject=whichObject[i];
    2111           CbcObject * object = model->modifiableObject(iObject);
    2112           CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
    2113             dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
    2114           which[needed++]=dynamicObject->columnNumber();
    2115         }
     2083      {
     2084        // off penalties if too much
     2085        double needed = neededPenalties;
     2086        needed *= numberRows;
     2087        if (needed>1.0e6) {
     2088          numberPenalties=0;
     2089          neededPenalties=0;
     2090        }
     2091      }
     2092      if (osiclp&&numberPenalties&&neededPenalties) {
    21162093        which--;
    2117         which[0]=needed;
    2118         assert (needed);
     2094        which[0]=neededPenalties;
    21192095        osiclp->passInRanges(which);
    21202096        // Mark hot start and get ranges
     
    21302106        kPass++;
    21312107        osiclp->passInRanges(NULL);
    2132         needed2=0;
    2133         int put=0;
    21342108        const double * downCost=osiclp->upRange();
    21352109        const double * upCost=osiclp->downRange();
    2136         // Bug - so switch off for now
    2137         double distanceToCutoff=COIN_DBL_MAX;
    21382110        //printf("numberTodo %d needed %d numberpenalties %d\n",numberToDo,needed,numberPenalties);
    2139         for ( i=0;i<numberToDo;i++) {
    2140           int iObject = whichObject[i];
     2111        double invTrust = 1.0/((double) numberBeforeTrust);
     2112        for (int i=0;i<neededPenalties;i++) {
     2113          int j = objectMark[i];
     2114          int iObject = whichObject[j];
    21412115          CbcObject * object = model->modifiableObject(iObject);
    21422116          CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
    21432117            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
    21442118          int iSequence=dynamicObject->columnNumber();
    2145           double estimate = sort[i];
    2146           if (i<needed) {
    2147             // We have looked at penalties
    2148             int iAction=0;
    2149             double value = saveSolution[iSequence];
    2150             value -= floor(value);
    2151             double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
    2152             double downPenalty = CoinMin(downCost[i],1.0e110)*value;
    2153             if (!numberBeforeTrust) {
    2154               // override
    2155               downEstimate[iObject]=downPenalty;
    2156               upEstimate[iObject]=upPenalty;
     2119          double value = saveSolution[iSequence];
     2120          value -= floor(value);
     2121          double upPenalty = CoinMin(upCost[i],1.0e110)*(1.0-value);
     2122          double downPenalty = CoinMin(downCost[i],1.0e110)*value;
     2123          if (!numberBeforeTrust) {
     2124            // override
     2125            downEstimate[iObject]=downPenalty;
     2126            upEstimate[iObject]=upPenalty;
     2127          } else {
     2128            int numberThisDown = dynamicObject->numberTimesDown();
     2129            if (numberThisDown<numberBeforeTrust) {
     2130              double fraction = ((double) numberThisDown)*invTrust;
     2131              downEstimate[iObject] = fraction*downEstimate[iObject]+(1.0-fraction)*downPenalty;
    21572132            }
    2158             if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
    2159               printf("%d pen down ps %g -> %g up ps %g -> %g\n",
    2160                      iObject,downCost[i],downPenalty,upCost[i],upPenalty);
    2161             if (upPenalty>distanceToCutoff) {
    2162               if(downPenalty>distanceToCutoff) {
    2163                 //printf("%d infeas both penalty %g %g\n",iObject,upPenalty,downPenalty);
    2164                 iAction=3;
    2165               } else {
    2166                 //printf("%d infeas up penalty %g %g\n",iObject,upPenalty,downPenalty);
    2167                 iAction=2;
    2168               }
    2169             } else if(downPenalty>distanceToCutoff) {
    2170               //printf("%d infeas down penalty %g %g\n",iObject,upPenalty,downPenalty);
    2171               iAction=1;
     2133            int numberThisUp = dynamicObject->numberTimesUp();
     2134            if (numberThisUp<numberBeforeTrust) {
     2135              double fraction = ((double) numberThisUp)*invTrust;
     2136              upEstimate[iObject] = fraction*upEstimate[iObject]+(1.0-fraction)*upPenalty;
    21722137            }
    2173             if (iAction) {
    2174               if (iAction==3) {
    2175                 //printf("%d infeas both penalty %g %g\n",iObject,upPenalty,downPenalty);
    2176                 anyAction=-2;
    2177                 delete ws;
    2178                 ws=NULL;
    2179                 break;
    2180               } else if (iAction==2) {
    2181                 //printf("%d infeas up penalty %g %g\n",iObject,upPenalty,downPenalty);
    2182                 CbcStrongInfo choice;
    2183                 choice.objectNumber=iObject;
    2184                 choice.fix=-1;
    2185                 choice.possibleBranch=object->createBranch(-1);
    2186                 fixObject[numberToFix++]=choice;
    2187                 if (!anyAction)
    2188                   anyAction=-1;
    2189               } else {
    2190                 //printf("%d infeas down penalty %g %g\n",iObject,upPenalty,downPenalty);
    2191                 CbcStrongInfo choice;
    2192                 choice.objectNumber=iObject;
    2193                 choice.fix=1;
    2194                 choice.possibleBranch=object->createBranch(1);
    2195                 fixObject[numberToFix++]=choice;
    2196                 if (!anyAction)
    2197                   anyAction=-1;
    2198               }
    2199             } else {
    2200               double add = objectMark[iObject];
    2201               double trueEstimate = sort[i] - add*average;
    2202               trueEstimate=0.0; // temp
    2203               estimate = trueEstimate-scaleFactor*(upPenalty+downPenalty);
    2204               sort[put]=estimate;
    2205               whichObject[put++]=iObject;
    2206               needed2=put;
    2207             }
    2208           } else {
    2209             // just estimate
    2210             double add = objectMark[iObject];
    2211             double trueEstimate = sort[i] - add*average;
    2212             sort[put]=trueEstimate;
    2213             whichObject[put++]=iObject;
    2214           }
    2215         }
    2216         numberToDo=put;
     2138          }
     2139          sort[j] = - CoinMin(downEstimate[iObject],upEstimate[iObject]);
     2140          if ((numberNodes%PRINT_STUFF)==0&&PRINT_STUFF>0)
     2141            printf("%d pen down ps %g -> %g up ps %g -> %g\n",
     2142                   iObject,downCost[i],downPenalty,upCost[i],upPenalty);
     2143        }
    22172144      } else {
    22182145        if (!skipAll) {
    22192146          // Mark hot start
    22202147          solver->markHotStart();
    2221           if (solver->isProvenPrimalInfeasible())
    2222             printf("**** Hot start says node infeasible\n");
     2148          //if (solver->isProvenPrimalInfeasible())
     2149          //printf("**** Hot start says node infeasible\n");
    22232150        }
    22242151        // make sure best will be first
     
    22392166        numberToDo=0; // skip as we will be trying again
    22402167#endif
    2241       // Sort (but just those re-computed)
    2242       CoinSort_2(sort,sort+needed2,whichObject);
     2168      // Sort
     2169      CoinSort_2(sort,sort+numberToDo,whichObject);
     2170      // Change in objective opposite infeasible
     2171      double worstFeasible=0.0;
    22432172      // Just first if strong off
    22442173      if (!numberStrong)
     
    22682197      }
    22692198#endif
     2199      if (depth_<10&&numberStrong) {
     2200        doQuickly=false;
     2201        if (depth_<7)
     2202          numberStrong *=3;
     2203        if (!depth_)
     2204          numberStrong=CoinMin(6*numberStrong,numberToDo);
     2205        numberTest=numberStrong;
     2206        model->setStateOfSearch(2); // use min min
     2207      }
     2208      // could adjust using average iterations per branch
     2209      // double average = ((double)model->getIterationCount())/
     2210      //((double) model->getNodeCount()+1.0);
    22702211      // if too many and big then just do 10 its
    2271       if (!skipAll&&stateOfSearch) {
    2272         if (numberNotTrusted>3*numberStrong&&numberRows>250&&numberColumns>1000)
    2273           solver->setIntParam(OsiMaxNumIterationHotStart,10);
     2212      if (!skipAll&&saveStateOfSearch) {
     2213        //if (numberNotTrusted>3*numberStrong&&numberRows>250&&numberColumns>1000)
     2214          // off solver->setIntParam(OsiMaxNumIterationHotStart,10);
    22742215      }
    22752216      double distanceToCutoff=model->getCutoff()-objectiveValue_;
    2276       // larger and make negative for test
    2277       distanceToCutoff *= -100.0;
     2217      // make negative for test
     2218      distanceToCutoff = - distanceToCutoff;
     2219      if (numberObjects>-100) {
     2220        // larger
     2221        distanceToCutoff *= 100.0;
     2222      }
    22782223      if (skipAll)
    22792224        distanceToCutoff = -COIN_DBL_MAX;
     2225      // Do at least 5 strong
     2226      if (numberColumns<1000&&(depth_<15||numberNodes<1000000))
     2227        numberTest = CoinMax(numberTest,5);
    22802228      for ( iDo=0;iDo<numberToDo;iDo++) {
    22812229        CbcStrongInfo choice;
     
    25242472            // up feasible, down infeasible
    25252473            anyAction=-1;
     2474            worstFeasible = CoinMax(worstFeasible,choice.upMovement);
    25262475            //printf("Down infeasible for choice %d sequence %d\n",i,
    25272476            // model->object(choice.objectNumber)->columnNumber());
     
    25482497            // down feasible, up infeasible
    25492498            anyAction=-1;
     2499            worstFeasible = CoinMax(worstFeasible,choice.downMovement);
    25502500            //printf("Up infeasible for choice %d sequence %d\n",i,
    25512501            // model->object(choice.objectNumber)->columnNumber());
     
    25942544        }
    25952545      }
     2546      double averageChange = model->sumChangeObjective()/((double) model->getNodeCount());
     2547      if (depth_<10||worstFeasible>0.2*averageChange)
     2548        solveAll=false;
    25962549      if (model->messageHandler()->logLevel()>3||false) {
    25972550        if (anyAction==-2)
     
    26342587          }
    26352588          bool feasible=true;
    2636 #if ACTION <2
    2637           // can do quick optimality check
    2638           int easy=2;
    2639           solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
    2640           solver->resolve() ;
    2641           solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
    2642           feasible = solver->isProvenOptimal();
    2643           if (feasible) {
    2644             anyAction=0;
    2645             memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
    2646             model->reserveCurrentSolution(saveSolution);
    2647             memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
    2648             memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
    2649             // See if candidate still possible
    2650             if (branch_) {
    2651               const CbcObject * object = model->object(bestChoice);
    2652               int preferredWay;
    2653               double infeasibility = object->infeasibility(preferredWay);
    2654               if (!infeasibility) {
    2655                 // take out
    2656                 delete branch_;
    2657                 branch_=NULL;
    2658               } else {
    2659                 branch_->way(preferredWay);
     2589#if ACTION <2
     2590          if (solveAll) {
     2591            // can do quick optimality check
     2592            int easy=2;
     2593            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
     2594            solver->resolve() ;
     2595            solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
     2596            feasible = solver->isProvenOptimal();
     2597            if (feasible) {
     2598              anyAction=0;
     2599              memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
     2600              model->reserveCurrentSolution(saveSolution);
     2601              memcpy(saveLower,solver->getColLower(),numberColumns*sizeof(double));
     2602              memcpy(saveUpper,solver->getColUpper(),numberColumns*sizeof(double));
     2603              model->setPointers(solver);
     2604              // See if candidate still possible
     2605              if (branch_) {
     2606                const CbcObject * object = model->object(bestChoice);
     2607                int preferredWay;
     2608                double infeasibility = object->infeasibility(preferredWay);
     2609                if (!infeasibility) {
     2610                  // take out
     2611                  delete branch_;
     2612                  branch_=NULL;
     2613                } else {
     2614                  branch_->way(preferredWay);
     2615                }
    26602616              }
     2617            } else {
     2618              anyAction=-2;
     2619              finished=true;
    26612620            }
    2662           } else {
    2663             anyAction=-2;
    2664             finished=true;
    26652621          }
    26662622#endif
     
    26802636    }
    26812637  }
    2682  
     2638  //if (numberToFix&&depth_<5)
     2639  //printf("%d fixed by strong at depth %d\n",numberToFix,depth_);
    26832640  // Set guessed solution value
    26842641  guessedObjectiveValue_ = objectiveValue_+estimatedDegradation;
     
    26972654  delete [] upEstimate;
    26982655  delete [] downEstimate;
     2656  if (osiclp)
     2657    osiclp->setSpecialOptions(saveClpOptions);
    26992658 
    27002659  // restore solution
     
    27022661  model->reserveCurrentSolution(saveSolution);
    27032662  delete [] saveSolution;
     2663  model->setStateOfSearch(saveStateOfSearch);
    27042664  return anyAction;
     2665}
     2666int CbcNode::analyze (CbcModel *model, double * results)
     2667{
     2668  int i;
     2669  int numberIterationsAllowed = model->numberAnalyzeIterations();
     2670  OsiSolverInterface * solver = model->solver();
     2671  objectiveValue_ = solver->getObjSense()*solver->getObjValue();
     2672  double cutoff =model->getCutoff();
     2673  const double * lower = solver->getColLower();
     2674  const double * upper = solver->getColUpper();
     2675  const double * dj = solver->getReducedCost();
     2676  int numberObjects = model->numberObjects();
     2677  int numberColumns = model->getNumCols();
     2678  // Initialize arrays
     2679  int numberIntegers = model->numberIntegers();
     2680  int * back = new int[numberColumns];
     2681  const int * integerVariable = model->integerVariable();
     2682  for (i=0;i<numberColumns;i++)
     2683    back[i]=-1;
     2684  // What results is
     2685  double * newLower = results;
     2686  double * objLower = newLower+numberIntegers;
     2687  double * newUpper = objLower+numberIntegers;
     2688  double * objUpper = newUpper+numberIntegers;
     2689  for (i=0;i<numberIntegers;i++) {
     2690    int iColumn = integerVariable[i];
     2691    back[iColumn]=i;
     2692    newLower[i]=0.0;
     2693    objLower[i]=COIN_DBL_MAX;
     2694    newUpper[i]=0.0;
     2695    objUpper[i]=COIN_DBL_MAX;
     2696  }
     2697  double * saveUpper = new double[numberColumns];
     2698  double * saveLower = new double[numberColumns];
     2699  int anyAction=0;
     2700  // Save solution in case heuristics need good solution later
     2701 
     2702  double * saveSolution = new double[numberColumns];
     2703  memcpy(saveSolution,solver->getColSolution(),numberColumns*sizeof(double));
     2704  model->reserveCurrentSolution(saveSolution);
     2705  for (i=0;i<numberColumns;i++) {
     2706    saveLower[i] = lower[i];
     2707    saveUpper[i] = upper[i];
     2708  }
     2709  // Get arrays to sort
     2710  double * sort = new double[numberObjects];
     2711  int * whichObject = new int[numberObjects];
     2712  int numberToFix=0;
     2713  int numberToDo=0;
     2714     
     2715  // compute current state
     2716  int numberObjectInfeasibilities; // just odd ones
     2717  int numberIntegerInfeasibilities;
     2718  model->feasibleSolution(
     2719                          numberIntegerInfeasibilities,
     2720                          numberObjectInfeasibilities);
     2721     
     2722  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
     2723  int saveClpOptions=0;
     2724  if (osiclp) {
     2725    // for faster hot start
     2726    saveClpOptions = osiclp->specialOptions();
     2727    osiclp->setSpecialOptions(saveClpOptions|1024);
     2728  }
     2729  /*
     2730    Scan for branching objects that indicate infeasibility. Choose candidates
     2731    using priority as the first criteria, then integer infeasibility.
     2732   
     2733    The algorithm is to fill the array with a set of good candidates (by
     2734    infeasibility) with priority bestPriority.  Finding a candidate with
     2735    priority better (less) than bestPriority flushes the choice array. (This
     2736    serves as initialization when the first candidate is found.)
     2737   
     2738  */
     2739  numberToDo=0;
     2740  for (i=0;i<numberObjects;i++) {
     2741    CbcObject * object = model->modifiableObject(i);
     2742    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
     2743      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
     2744    if(!dynamicObject)
     2745      continue;
     2746    int preferredWay;
     2747    double infeasibility = object->infeasibility(preferredWay);
     2748    int iColumn = dynamicObject->columnNumber();
     2749    if (saveUpper[iColumn]==saveLower[iColumn])
     2750      continue;
     2751    if (infeasibility)
     2752      sort[numberToDo]=-1.0e10-infeasibility;
     2753    else
     2754      sort[numberToDo]=-fabs(dj[iColumn]);
     2755    whichObject[numberToDo++]=i;
     2756  }
     2757  // Save basis
     2758  CoinWarmStart * ws = solver->getWarmStart();
     2759  int saveLimit;
     2760  solver->getIntParam(OsiMaxNumIterationHotStart,saveLimit);
     2761  if (saveLimit<200)
     2762    solver->setIntParam(OsiMaxNumIterationHotStart,200);
     2763  // Mark hot start
     2764  solver->markHotStart();
     2765  // Sort
     2766  CoinSort_2(sort,sort+numberToDo,whichObject);
     2767  //double distanceToCutoff=model->getCutoff()-objectiveValue_;
     2768  double * currentSolution = model->currentSolution();
     2769  double integerTolerance =
     2770    model->getDblParam(CbcModel::CbcIntegerTolerance);
     2771  double objMin = 1.0e50;
     2772  double objMax = -1.0e50;
     2773  bool needResolve=false;
     2774  for (int iDo=0;iDo<numberToDo;iDo++) {
     2775    CbcStrongInfo choice;
     2776    int iObject = whichObject[iDo];
     2777    CbcObject * object = model->modifiableObject(iObject);
     2778    CbcSimpleIntegerDynamicPseudoCost * dynamicObject =
     2779      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object) ;
     2780    int iColumn = dynamicObject->columnNumber();
     2781    int preferredWay;
     2782    object->infeasibility(preferredWay);
     2783    double value = currentSolution[iColumn];
     2784    double nearest = floor(value+0.5);
     2785    double lowerValue = floor(value);
     2786    bool satisfied=false;
     2787    if (fabs(value-nearest)<=integerTolerance) {
     2788      satisfied=true;
     2789      double newValue;
     2790      if (nearest<saveUpper[iColumn]) {
     2791        newValue = nearest + 1.0001*integerTolerance;
     2792        lowerValue = nearest;
     2793      } else {
     2794        newValue = nearest - 1.0001*integerTolerance;
     2795        lowerValue = nearest-1;
     2796      }
     2797      currentSolution[iColumn]=newValue;
     2798    }
     2799    double upperValue = lowerValue+1.0;
     2800    choice.possibleBranch=object->createBranch(preferredWay);
     2801    currentSolution[iColumn]=value;
     2802    // Save which object it was
     2803    choice.objectNumber=iObject;
     2804    choice.numIntInfeasUp = numberUnsatisfied_;
     2805    choice.numIntInfeasDown = numberUnsatisfied_;
     2806    choice.downMovement = 0.0;
     2807    choice.upMovement = 0.0;
     2808    choice.numItersDown = 0;
     2809    choice.numItersUp = 0;
     2810    choice.fix=0; // say not fixed
     2811    double objectiveChange ;
     2812    double newObjectiveValue=1.0e100;
     2813    int j;
     2814    // status is 0 finished, 1 infeasible and other
     2815    int iStatus;
     2816    /*
     2817      Try the down direction first. (Specify the initial branching alternative as
     2818      down with a call to way(-1). Each subsequent call to branch() performs the
     2819      specified branch and advances the branch object state to the next branch
     2820      alternative.)
     2821    */
     2822    choice.possibleBranch->way(-1) ;
     2823    choice.possibleBranch->branch() ;
     2824    if (fabs(value-lowerValue)>integerTolerance) {
     2825      solver->solveFromHotStart() ;
     2826      /*
     2827        We now have an estimate of objective degradation that we can use for strong
     2828        branching. If we're over the cutoff, the variable is monotone up.
     2829        If we actually made it to optimality, check for a solution, and if we have
     2830        a good one, call setBestSolution to process it. Note that this may reduce the
     2831        cutoff, so we check again to see if we can declare this variable monotone.
     2832      */
     2833      if (solver->isProvenOptimal())
     2834        iStatus=0; // optimal
     2835      else if (solver->isIterationLimitReached()
     2836               &&!solver->isDualObjectiveLimitReached())
     2837        iStatus=2; // unknown
     2838      else
     2839        iStatus=1; // infeasible
     2840      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
     2841      choice.numItersDown = solver->getIterationCount();
     2842      numberIterationsAllowed -= choice.numItersDown;
     2843      objectiveChange = newObjectiveValue  - objectiveValue_;
     2844      if (!iStatus) {
     2845        choice.finishedDown = true ;
     2846        if (newObjectiveValue>=cutoff) {
     2847          objectiveChange = 1.0e100; // say infeasible
     2848        } else {
     2849          // See if integer solution
     2850          if (model->feasibleSolution(choice.numIntInfeasDown,
     2851                                      choice.numObjInfeasDown)
     2852              &&model->problemFeasibility()->feasible(model,-1)>=0) {
     2853            model->setBestSolution(CBC_STRONGSOL,
     2854                                   newObjectiveValue,
     2855                                   solver->getColSolution()) ;
     2856            model->setLastHeuristic(NULL);
     2857            model->incrementUsed(solver->getColSolution());
     2858            cutoff =model->getCutoff();
     2859            if (newObjectiveValue >= cutoff)    //  *new* cutoff
     2860              objectiveChange = 1.0e100 ;
     2861          }
     2862        }
     2863      } else if (iStatus==1) {
     2864        objectiveChange = 1.0e100 ;
     2865      } else {
     2866        // Can't say much as we did not finish
     2867        choice.finishedDown = false ;
     2868      }
     2869      choice.downMovement = objectiveChange ;
     2870    }
     2871    // restore bounds
     2872    for ( j=0;j<numberColumns;j++) {
     2873      if (saveLower[j] != lower[j])
     2874        solver->setColLower(j,saveLower[j]);
     2875      if (saveUpper[j] != upper[j])
     2876        solver->setColUpper(j,saveUpper[j]);
     2877    }
     2878    // repeat the whole exercise, forcing the variable up
     2879    choice.possibleBranch->branch();
     2880    if (fabs(value-upperValue)>integerTolerance) {
     2881      solver->solveFromHotStart() ;
     2882      /*
     2883        We now have an estimate of objective degradation that we can use for strong
     2884        branching. If we're over the cutoff, the variable is monotone up.
     2885        If we actually made it to optimality, check for a solution, and if we have
     2886        a good one, call setBestSolution to process it. Note that this may reduce the
     2887        cutoff, so we check again to see if we can declare this variable monotone.
     2888      */
     2889      if (solver->isProvenOptimal())
     2890        iStatus=0; // optimal
     2891      else if (solver->isIterationLimitReached()
     2892               &&!solver->isDualObjectiveLimitReached())
     2893        iStatus=2; // unknown
     2894      else
     2895        iStatus=1; // infeasible
     2896      newObjectiveValue = solver->getObjSense()*solver->getObjValue();
     2897      choice.numItersUp = solver->getIterationCount();
     2898      numberIterationsAllowed -= choice.numItersUp;
     2899      objectiveChange = newObjectiveValue  - objectiveValue_;
     2900      if (!iStatus) {
     2901        choice.finishedUp = true ;
     2902        if (newObjectiveValue>=cutoff) {
     2903          objectiveChange = 1.0e100; // say infeasible
     2904        } else {
     2905          // See if integer solution
     2906          if (model->feasibleSolution(choice.numIntInfeasUp,
     2907                                      choice.numObjInfeasUp)
     2908              &&model->problemFeasibility()->feasible(model,-1)>=0) {
     2909            model->setBestSolution(CBC_STRONGSOL,
     2910                                   newObjectiveValue,
     2911                                   solver->getColSolution()) ;
     2912            model->setLastHeuristic(NULL);
     2913            model->incrementUsed(solver->getColSolution());
     2914            cutoff =model->getCutoff();
     2915            if (newObjectiveValue >= cutoff)    //  *new* cutoff
     2916              objectiveChange = 1.0e100 ;
     2917          }
     2918        }
     2919      } else if (iStatus==1) {
     2920        objectiveChange = 1.0e100 ;
     2921      } else {
     2922        // Can't say much as we did not finish
     2923        choice.finishedUp = false ;
     2924      }
     2925      choice.upMovement = objectiveChange ;
     2926     
     2927      // restore bounds
     2928      for ( j=0;j<numberColumns;j++) {
     2929        if (saveLower[j] != lower[j])
     2930          solver->setColLower(j,saveLower[j]);
     2931        if (saveUpper[j] != upper[j])
     2932          solver->setColUpper(j,saveUpper[j]);
     2933      }
     2934    }
     2935    // If objective goes above certain amount we can set bound
     2936    int jInt = back[iColumn];
     2937    newLower[jInt]=upperValue;
     2938    objLower[jInt]=choice.downMovement+objectiveValue_;
     2939    newUpper[jInt]=lowerValue;
     2940    objUpper[jInt]=choice.upMovement+objectiveValue_;
     2941    objMin = CoinMin(CoinMin(objLower[jInt],objUpper[jInt]),objMin);
     2942    /*
     2943      End of evaluation for this candidate variable. Possibilities are:
     2944      * Both sides below cutoff; this variable is a candidate for branching.
     2945      * Both sides infeasible or above the objective cutoff: no further action
     2946      here. Break from the evaluation loop and assume the node will be purged
     2947      by the caller.
     2948      * One side below cutoff: Install the branch (i.e., fix the variable). Break
     2949      from the evaluation loop and assume the node will be reoptimised by the
     2950      caller.
     2951    */
     2952    if (choice.upMovement<1.0e100) {
     2953      if(choice.downMovement<1.0e100) {
     2954        objMax = CoinMax(CoinMax(objLower[jInt],objUpper[jInt]),objMax);
     2955        // In case solution coming in was odd
     2956        choice.upMovement = CoinMax(0.0,choice.upMovement);
     2957        choice.downMovement = CoinMax(0.0,choice.downMovement);
     2958        // feasible -
     2959        model->messageHandler()->message(CBC_STRONG,model->messages())
     2960          << iObject << iColumn
     2961          <<choice.downMovement<<choice.numIntInfeasDown
     2962          <<choice.upMovement<<choice.numIntInfeasUp
     2963          <<value
     2964          <<CoinMessageEol;
     2965      } else {
     2966        // up feasible, down infeasible
     2967        anyAction=-1;
     2968        if (!satisfied)
     2969          needResolve=true;
     2970        choice.fix=1;
     2971        numberToFix++;
     2972        saveLower[iColumn]=upperValue;
     2973        solver->setColLower(iColumn,upperValue);
     2974      }
     2975    } else {
     2976      if(choice.downMovement<1.0e100) {
     2977        // down feasible, up infeasible
     2978        anyAction=-1;
     2979        if (!satisfied)
     2980          needResolve=true;
     2981        choice.fix=-1;
     2982        numberToFix++;
     2983        saveUpper[iColumn]=lowerValue;
     2984        solver->setColUpper(iColumn,lowerValue);
     2985      } else {
     2986        // neither side feasible
     2987        anyAction=-2;
     2988        delete choice.possibleBranch;
     2989        printf("Both infeasible for choice %d sequence %d\n",i,
     2990               model->object(choice.objectNumber)->columnNumber());
     2991        delete ws;
     2992        ws=NULL;
     2993        solver->writeMps("bad");
     2994        numberToFix=-1;
     2995        break;
     2996      }
     2997    }
     2998    if (numberIterationsAllowed<=0)
     2999      break;
     3000    //printf("obj %d, col %d, down %g up %g value %g\n",iObject,iColumn,
     3001    //     choice.downMovement,choice.upMovement,value);
     3002  }
     3003  printf("Best possible solution %g, can fix more if solution of %g found\n",
     3004         objMin,objMax);
     3005  model->setNumberAnalyzeIterations(numberIterationsAllowed);
     3006  // Delete the snapshot
     3007  solver->unmarkHotStart();
     3008  // back to normal
     3009  solver->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
     3010  solver->setIntParam(OsiMaxNumIterationHotStart,saveLimit);
     3011  // restore basis
     3012  solver->setWarmStart(ws);
     3013  delete ws;
     3014   
     3015  delete [] sort;
     3016  delete [] whichObject;
     3017  delete [] saveLower;
     3018  delete [] saveUpper;
     3019  delete [] back;
     3020  // restore solution
     3021  solver->setColSolution(saveSolution);
     3022  if (osiclp)
     3023    osiclp->setSpecialOptions(saveClpOptions);
     3024  model->reserveCurrentSolution(saveSolution);
     3025  delete [] saveSolution;
     3026  if (needResolve)
     3027    solver->resolve();
     3028  return numberToFix;
    27053029}
    27063030
  • trunk/include/CbcMessage.hpp

    r181 r208  
    5353  CBC_END_SUB,
    5454  CBC_HEURISTIC_SOLUTION,
     55  CBC_CUTS_STATS,
    5556  CBC_DUMMY_END
    5657};
  • trunk/include/CbcModel.hpp

    r202 r208  
    33#ifndef CbcModel_H
    44#define CbcModel_H
     5#ifndef COIN_USE_CLP
     6#define COIN_USE_CLP
     7#endif
    58
    69#include <string>
    710#include <vector>
    8 
    911#include "CoinFinite.hpp"
    1012#include "CoinMessageHandler.hpp"
     
    1416#include "CbcCompareBase.hpp"
    1517#include "CbcMessage.hpp"
     18#include "ClpEventHandler.hpp"
    1619
    1720//class OsiSolverInterface;
     
    590593  inline int numberPenalties() const
    591594  { return numberPenalties_;};
     595  /// Number of analyze iterations to do
     596  inline void setNumberAnalyzeIterations(int number)
     597  { numberAnalyzeIterations_=number;};
     598  inline int numberAnalyzeIterations() const
     599  { return numberAnalyzeIterations_;};
    592600  /** Get scale factor to make penalties match strong.
    593601      Should/will be computed */
     
    9931001  inline double rootObjectiveAfterCuts() const
    9941002  { return continuousObjective_;};
     1003  /// Sum of Changes to objective by first solve
     1004  inline double sumChangeObjective() const
     1005  { return sumChangeObjective1_;};
    9951006  /** Number of times global cuts violated.  When global cut pool then this
    9961007      should be kept for each cut and type of cut */
     
    10831094
    10841095    Fixes integer variables at their current value based on reduced cost
    1085     penalties.
    1086   */
    1087   void reducedCostFix() ;
     1096    penalties.  Returns number fixed
     1097  */
     1098  int reducedCostFix() ;
    10881099
    10891100  /** Return an empty basis object of the specified size
     
    11181129    If it turns out that the node should really be fathomed by bound,
    11191130    addCuts() simply treats all the cuts as loose as it does the bookkeeping.
    1120   */
    1121   int addCuts(CbcNode * node, CoinWarmStartBasis *&lastws);
     1131
     1132    canFix true if extra information being passed
     1133  */
     1134  int addCuts(CbcNode * node, CoinWarmStartBasis *&lastws,bool canFix);
    11221135
    11231136  /** Traverse the tree from node to root and prep the model
     
    11591172  inline int stateOfSearch() const
    11601173  { return stateOfSearch_;};
     1174  inline void setStateOfSearch(int state)
     1175  { stateOfSearch_=state;};
    11611176
    11621177  /// Get the number of cut generators
     
    12881303  inline int logLevel() const
    12891304  { return handler_->logLevel();};
     1305   /// Pass in Event handler (cloned and deleted at end)
     1306   void passInEventHandler(const ClpEventHandler * eventHandler);
     1307   /// Event handler
     1308  ClpEventHandler * eventHandler() const;
    12901309  //@}
    12911310  //---------------------------------------------------------------------------
     
    15841603      Should/will be computed */
    15851604  double penaltyScaleFactor_;
     1605  /// Number of analyze iterations to do
     1606  int numberAnalyzeIterations_;
     1607  /// Arrays with analysis results
     1608  double * analyzeResults_;
    15861609  /// Number of nodes infeasible by normal branching (before cuts)
    15871610  int numberInfeasibleNodes_;
  • trunk/include/CbcNode.hpp

    r146 r208  
    473473                    CbcNode * lastNode,
    474474                    int numberPassesLeft);
    475  
     475  int analyze(CbcModel * model,double * results);
    476476  /// Decrement active cut counts
    477477  void decrementCuts(int change=1);
Note: See TracChangeset for help on using the changeset viewer.