source: trunk/CbcModel.cpp @ 267

Last change on this file since 267 was 267, checked in by forrest, 14 years ago

for nonlinear

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 302.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7#include <string>
8//#define CBC_DEBUG 1
9//#define CHECK_CUT_COUNTS
10//#define CHECK_NODE_FULL
11//#define NODE_LOG
12#include <cassert>
13#include <cmath>
14#include <cfloat>
15#ifdef COIN_USE_CLP
16// include Presolve from Clp
17#include "ClpPresolve.hpp"
18#include "OsiClpSolverInterface.hpp"
19#include "ClpEventHandler.hpp"
20#endif
21
22#include "OsiSolverInterface.hpp"
23#include "OsiAuxInfo.hpp"
24#include "OsiSolverBranch.hpp"
25#include "CoinWarmStartBasis.hpp"
26#include "CoinPackedMatrix.hpp"
27#include "CoinHelperFunctions.hpp"
28#include "CbcBranchActual.hpp"
29#include "CbcBranchDynamic.hpp"
30#include "CbcHeuristic.hpp"
31#include "CbcModel.hpp"
32#include "CbcTreeLocal.hpp"
33#include "CbcStatistics.hpp"
34#include "CbcStrategy.hpp"
35#include "CbcMessage.hpp"
36#include "OsiRowCut.hpp"
37#include "OsiColCut.hpp"
38#include "OsiRowCutDebugger.hpp"
39#include "OsiCuts.hpp"
40#include "CbcCountRowCut.hpp"
41#include "CbcCutGenerator.hpp"
42#include "CbcFeasibilityBase.hpp"
43// include Probing
44#include "CglProbing.hpp"
45// include preprocessing
46#include "CglPreProcess.hpp"
47
48#include "CoinTime.hpp"
49
50#include "CbcCompareActual.hpp"
51#include "CbcTree.hpp"
52/* Various functions local to CbcModel.cpp */
53
54namespace {
55
56//-------------------------------------------------------------------
57// Returns the greatest common denominator of two
58// positive integers, a and b, found using Euclid's algorithm
59//-------------------------------------------------------------------
60static int gcd(int a, int b) 
61{
62  int remainder = -1;
63  // make sure a<=b (will always remain so)
64  if(a > b) {
65    // Swap a and b
66    int temp = a;
67    a = b;
68    b = temp;
69  }
70  // if zero then gcd is nonzero (zero may occur in rhs of packed)
71  if (!a) {
72    if (b) {
73      return b;
74    } else {
75      printf("**** gcd given two zeros!!\n");
76      abort();
77    }
78  }
79  while (remainder) {
80    remainder = b % a;
81    b = a;
82    a = remainder;
83  }
84  return b;
85}
86
87
88
89#ifdef CHECK_NODE_FULL
90
91/*
92  Routine to verify that tree linkage is correct. The invariant that is tested
93  is
94
95  reference count = (number of actual references) + (number of branches left)
96
97  The routine builds a set of paired arrays, info and count, by traversing the
98  tree. Each CbcNodeInfo is recorded in info, and the number of times it is
99  referenced (via the parent field) is recorded in count. Then a final check is
100  made to see if the numberPointingToThis_ field agrees.
101*/
102
103void verifyTreeNodes (const CbcTree * branchingTree, const CbcModel &model)
104
105{ printf("*** CHECKING tree after %d nodes\n",model.getNodeCount()) ;
106
107  int j ;
108  int nNodes = branchingTree->size() ;
109# define MAXINFO 1000
110  int *count = new int [MAXINFO] ;
111  CbcNodeInfo **info = new CbcNodeInfo*[MAXINFO] ;
112  int nInfo = 0 ;
113/*
114  Collect all CbcNodeInfo objects in info, by starting from each live node and
115  traversing back to the root. Nodes in the live set should have unexplored
116  branches remaining.
117
118  TODO: The `while (nodeInfo)' loop could be made to break on reaching a
119        common ancester (nodeInfo is found in info[k]). Alternatively, the
120        check could change to signal an error if nodeInfo is not found above a
121        common ancestor.
122*/
123  for (j = 0 ; j < nNodes ; j++)
124  { CbcNode *node = branchingTree->nodePointer(j) ;
125  if (!node)
126    continue;
127    CbcNodeInfo *nodeInfo = node->nodeInfo() ; 
128    int change = node->nodeInfo()->numberBranchesLeft() ;
129    assert(change) ;
130    while (nodeInfo)
131    { int k ;
132      for (k = 0 ; k < nInfo ; k++)
133      { if (nodeInfo == info[k]) break ; }
134      if (k == nInfo)
135      { assert(nInfo < MAXINFO) ;
136        nInfo++ ;
137        info[k] = nodeInfo ;
138        count[k] = 0 ; }
139      nodeInfo = nodeInfo->parent() ; } }
140/*
141  Walk the info array. For each nodeInfo, look up its parent in info and
142  increment the corresponding count.
143*/
144  for (j = 0 ; j < nInfo ; j++)
145  { CbcNodeInfo *nodeInfo = info[j] ;
146    nodeInfo = nodeInfo->parent() ;
147    if (nodeInfo)
148    { int k ;
149      for (k = 0 ; k < nInfo ; k++)
150      { if (nodeInfo == info[k]) break ; }
151      assert (k < nInfo) ;
152      count[k]++ ; } }
153/*
154  Walk the info array one more time and check that the invariant holds. The
155  number of references (numberPointingToThis()) should equal the sum of the
156  number of actual references (held in count[]) plus the number of potential
157  references (unexplored branches, numberBranchesLeft()).
158*/
159  for (j = 0;j < nInfo;j++) {
160    CbcNodeInfo * nodeInfo = info[j] ;
161    if (nodeInfo) {
162      int k ;
163      for (k = 0;k < nInfo;k++)
164        if (nodeInfo == info[k])
165          break ;
166      printf("Nodeinfo %x - %d left, %d count\n",
167             nodeInfo,
168             nodeInfo->numberBranchesLeft(),
169             nodeInfo->numberPointingToThis()) ;
170      assert(nodeInfo->numberPointingToThis() ==
171             count[k]+nodeInfo->numberBranchesLeft()) ; } }
172
173  delete [] count ;
174  delete [] info ;
175 
176  return ; }
177
178#endif  /* CHECK_NODE_FULL */
179
180
181
182#ifdef CHECK_CUT_COUNTS
183
184/*
185  Routine to verify that cut reference counts are correct.
186*/
187void verifyCutCounts (const CbcTree * branchingTree, CbcModel &model)
188
189{ printf("*** CHECKING cuts after %d nodes\n",model.getNodeCount()) ;
190
191  int j ;
192  int nNodes = branchingTree->size() ;
193
194/*
195  cut.tempNumber_ exists for the purpose of doing this verification. Clear it
196  in all cuts. We traverse the tree by starting from each live node and working
197  back to the root. At each CbcNodeInfo, check for cuts.
198*/
199  for (j = 0 ; j < nNodes ; j++)
200  { CbcNode *node = branchingTree->nodePointer(j) ;
201    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
202    assert (node->nodeInfo()->numberBranchesLeft()) ;
203    while (nodeInfo)
204    { int k ;
205      for (k = 0 ; k < nodeInfo->numberCuts() ; k++)
206      { CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
207        if (cut) cut->tempNumber_ = 0; }
208      nodeInfo = nodeInfo->parent() ; } }
209/*
210  Walk the live set again, this time collecting the list of cuts in use at each
211  node. addCuts1 will collect the cuts in model.addedCuts_. Take into account
212  that when we recreate the basis for a node, we compress out the slack cuts.
213*/
214  for (j = 0 ; j < nNodes ; j++)
215  { CoinWarmStartBasis *debugws = model.getEmptyBasis() ;
216    CbcNode *node = branchingTree->nodePointer(j) ;
217    CbcNodeInfo *nodeInfo = node->nodeInfo(); 
218    int change = node->nodeInfo()->numberBranchesLeft() ;
219    printf("Node %d %x (info %x) var %d way %d obj %g",j,node,
220           node->nodeInfo(),node->variable(),node->way(),
221           node->objectiveValue()) ;
222
223    model.addCuts1(node,debugws) ;
224
225    int i ;
226    int numberRowsAtContinuous = model.numberRowsAtContinuous() ;
227    CbcCountRowCut **addedCuts = model.addedCuts() ;
228    for (i = 0 ; i < model.currentNumberCuts() ; i++)
229    { CoinWarmStartBasis::Status status = 
230        debugws->getArtifStatus(i+numberRowsAtContinuous) ;
231      if (status != CoinWarmStartBasis::basic && addedCuts[i])
232      { addedCuts[i]->tempNumber_ += change ; } }
233
234    while (nodeInfo)
235    { nodeInfo = nodeInfo->parent() ;
236      if (nodeInfo) printf(" -> %x",nodeInfo); }
237    printf("\n") ;
238    delete debugws ; }
239/*
240  The moment of truth: We've tallied up the references by direct scan of the
241  search tree. Check for agreement with the count in the cut.
242
243  TODO: Rewrite to check and print mismatch only when tempNumber_ == 0?
244*/
245  for (j = 0 ; j < nNodes ; j++)
246  { CbcNode *node = branchingTree->nodePointer(j) ;
247    CbcNodeInfo *nodeInfo = node->nodeInfo(); 
248    while (nodeInfo)
249    { int k ;
250      for (k = 0 ; k < nodeInfo->numberCuts() ; k++)
251      { CbcCountRowCut *cut = nodeInfo->cuts()[k] ;
252        if (cut && cut->tempNumber_ >= 0)
253        { if (cut->tempNumber_ != cut->numberPointingToThis()) 
254            printf("mismatch %x %d %x %d %d\n",nodeInfo,k,
255                    cut,cut->tempNumber_,cut->numberPointingToThis()) ;
256          else
257            printf("   match %x %d %x %d %d\n", nodeInfo,k,
258                   cut,cut->tempNumber_,cut->numberPointingToThis()) ;
259          cut->tempNumber_ = -1 ; } }
260      nodeInfo = nodeInfo->parent() ; } }
261
262  return ; }
263
264#endif /* CHECK_CUT_COUNTS */
265
266}
267
268 /* End unnamed namespace for CbcModel.cpp */
269
270
271
272void 
273CbcModel::analyzeObjective ()
274/*
275  Try to find a minimum change in the objective function. The first scan
276  checks that there are no continuous variables with non-zero coefficients,
277  and grabs the largest objective coefficient associated with an unfixed
278  integer variable. The second scan attempts to scale up the objective
279  coefficients to a point where they are sufficiently close to integer that
280  we can pretend they are integer, and calculate a gcd over the coefficients
281  of interest. This will be the minimum increment for the scaled coefficients.
282  The final action is to scale the increment back for the original coefficients
283  and install it, if it's better than the existing value.
284
285  John's note: We could do better than this.
286
287  John's second note - apologies for changing s to z
288*/
289{ const double *objective = getObjCoefficients() ;
290  const double *lower = getColLower() ;
291  const double *upper = getColUpper() ;
292/*
293  Take a first scan to see if there are unfixed continuous variables in the
294  objective.  If so, the minimum objective change could be arbitrarily small.
295  Also pick off the maximum coefficient of an unfixed integer variable.
296
297  If the objective is found to contain only integer variables, set the
298  fathoming discipline to strict.
299*/
300  double maximumCost = 0.0 ;
301  bool possibleMultiple = true ;
302  int iColumn ;
303  int numberColumns = getNumCols() ;
304  for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
305  { if (upper[iColumn] > lower[iColumn]+1.0e-8)
306    { if (isInteger(iColumn)) 
307        maximumCost = CoinMax(maximumCost,fabs(objective[iColumn])) ;
308      else if (objective[iColumn]) 
309        possibleMultiple = false ; } }
310  setIntParam(CbcModel::CbcFathomDiscipline,possibleMultiple) ;
311/*
312  If a nontrivial increment is possible, try and figure it out. We're looking
313  for gcd(c<j>) for all c<j> that are coefficients of unfixed integer
314  variables. Since the c<j> might not be integers, try and inflate them
315  sufficiently that they look like integers (and we'll deflate the gcd
316  later).
317
318  2520.0 is used as it is a nice multiple of 2,3,5,7
319*/
320    if (possibleMultiple&&maximumCost)
321    { int increment = 0 ;
322      double multiplier = 2520.0 ;
323      while (10.0*multiplier*maximumCost < 1.0e8)
324        multiplier *= 10.0 ;
325
326      for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
327      { if (upper[iColumn] > lower[iColumn]+1.0e-8)
328        { if (isInteger(iColumn)&&objective[iColumn])
329          { double value = fabs(objective[iColumn])*multiplier ;
330            int nearest = (int) floor(value+0.5) ;
331            if (fabs(value-floor(value+0.5)) > 1.0e-8)
332            { increment = 0 ;
333              break ; }
334            else if (!increment)
335            { increment = nearest ; }
336            else
337            { increment = gcd(increment,nearest) ; } } } }
338/*
339  If the increment beats the current value for objective change, install it.
340*/
341      if (increment)
342      { double value = increment ;
343        double cutoff = getDblParam(CbcModel::CbcCutoffIncrement) ;
344        value /= multiplier ;
345        if (value*0.999 > cutoff)
346        { messageHandler()->message(CBC_INTEGERINCREMENT,
347                                          messages())
348            << value << CoinMessageEol ;
349          setDblParam(CbcModel::CbcCutoffIncrement,value*0.999) ; } } }
350
351  return ; 
352}
353#ifndef CBC_NEXT_VERSION
354
355
356/**
357  \todo
358  Normally, it looks like we enter here from command dispatch in the main
359  routine, after calling the solver for an initial solution
360  (CbcModel::initialSolve, which simply calls the solver's initialSolve
361  routine.) The first thing we do is call resolve. Presumably there are
362  circumstances where this is nontrivial? There's also a call from
363  CbcModel::originalModel (tied up with integer presolve), which should be
364  checked.
365
366*/
367
368/*
369  The overall flow can be divided into three stages:
370    * Prep: Check that the lp relaxation remains feasible at the root. If so,
371      do all the setup for B&C.
372    * Process the root node: Generate cuts, apply heuristics, and in general do
373      the best we can to resolve the problem without B&C.
374    * Do B&C search until we hit a limit or exhaust the search tree.
375 
376  Keep in mind that in general there is no node in the search tree that
377  corresponds to the active subproblem. The active subproblem is represented
378  by the current state of the model,  of the solver, and of the constraint
379  system held by the solver.
380*/
381
382void CbcModel::branchAndBound(int doStatistics) 
383
384{
385/*
386  Capture a time stamp before we start.
387*/
388  dblParam_[CbcStartSeconds] = CoinCpuTime();
389  strongInfo_[0]=0;
390  strongInfo_[1]=0;
391  strongInfo_[2]=0;
392  numberStrongIterations_ = 0;
393  // original solver (only set if pre-processing)
394  OsiSolverInterface * originalSolver=NULL;
395  // Set up strategies
396  if (strategy_) {
397    // May do preprocessing
398    originalSolver = solver_;
399    strategy_->setupOther(*this);
400    if (strategy_->preProcessState()) {
401      // pre-processing done
402      if (strategy_->preProcessState()<0) {
403        // infeasible
404        handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
405        status_ = 0 ;
406        secondaryStatus_ = 1;
407        originalContinuousObjective_ = COIN_DBL_MAX;
408        return ; 
409      }
410    } else {
411      //no preprocessing
412      originalSolver=NULL;
413    }
414    strategy_->setupCutGenerators(*this);
415    strategy_->setupHeuristics(*this);
416    // Set strategy print level to models
417    strategy_->setupPrinting(*this,handler_->logLevel());
418  }
419  eventHappened_=false;
420#ifdef COIN_USE_CLP
421  ClpEventHandler * eventHandler=NULL;
422  {
423    OsiClpSolverInterface * clpSolver
424      = dynamic_cast<OsiClpSolverInterface *> (solver_);
425    if (clpSolver) {
426      ClpSimplex * clpSimplex = clpSolver->getModelPtr();
427      eventHandler = clpSimplex->eventHandler();
428    }
429  }
430#endif
431 
432  if (!nodeCompare_)
433    nodeCompare_=new CbcCompareDefault();;
434  // See if hot start wanted
435  CbcCompareBase * saveCompare = NULL;
436  if (hotstartSolution_) {
437    if (strategy_&&strategy_->preProcessState()>0) {
438      CglPreProcess * process = strategy_->process();
439      assert (process);
440      int n = solver_->getNumCols();
441      const int * originalColumns = process->originalColumns();
442      // columns should be in order ... but
443      double * tempS = new double[n];
444      for (int i=0;i<n;i++) {
445        int iColumn = originalColumns[i];
446        tempS[i]=hotstartSolution_[iColumn];
447      }
448      delete [] hotstartSolution_;
449      hotstartSolution_=tempS;
450      if (hotstartPriorities_) {
451        int * tempP = new int [n];
452        for (int i=0;i<n;i++) {
453          int iColumn = originalColumns[i];
454          tempP[i]=hotstartPriorities_[iColumn];
455        }
456        delete [] hotstartPriorities_;
457        hotstartPriorities_=tempP;
458      }
459    }
460    saveCompare = nodeCompare_;
461    // depth first
462    nodeCompare_ = new CbcCompareDepth();
463  }
464  if (!problemFeasibility_)
465    problemFeasibility_=new CbcFeasibilityBase();
466# ifdef CBC_DEBUG
467  std::string problemName ;
468  solver_->getStrParam(OsiProbName,problemName) ;
469  printf("Problem name - %s\n",problemName.c_str()) ;
470  solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
471# endif
472/*
473  Assume we're done, and see if we're proven wrong.
474*/
475  status_ = 0 ;
476  secondaryStatus_ = 0;
477  phase_=0;
478/*
479  Scan the variables, noting the integer variables. Create an
480  CbcSimpleInteger object for each integer variable.
481*/
482  findIntegers(false) ;
483  // If dynamic pseudo costs then do
484  if (numberBeforeTrust_)
485    convertToDynamic();
486  // Set up char array to say if integer
487  delete [] integerInfo_;
488  {
489    int n = solver_->getNumCols();
490    integerInfo_ = new char [n];
491    for (int i=0;i<n;i++) {
492      if (solver_->isInteger(i))
493        integerInfo_[i]=1;
494      else
495        integerInfo_[i]=0;
496    }
497  }
498   
499/*
500  Ensure that objects on the lists of CbcObjects, heuristics, and cut
501  generators attached to this model all refer to this model.
502*/
503  synchronizeModel() ;
504  assert (!solverCharacteristics_);
505  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
506  if (solverCharacteristics) {
507    solverCharacteristics_ = solverCharacteristics;
508  } else {
509    // replace in solver
510    OsiBabSolver defaultC;
511    solver_->setAuxiliaryInfo(&defaultC);
512    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
513  }
514  solverCharacteristics_->setSolver(solver_);
515  // Set so we can tell we are in initial phase in resolve
516  continuousObjective_ = -COIN_DBL_MAX ;
517/*
518  Solve the relaxation.
519
520  Apparently there are circumstances where this will be non-trivial --- i.e.,
521  we've done something since initialSolve that's trashed the solution to the
522  continuous relaxation.
523*/
524  bool feasible;
525  // If NLP then we assume already solved outside branchAndbound
526  if (!solverCharacteristics_->solverType()) {
527    feasible=resolve() ;
528  } else {
529    // pick up given status
530    feasible = (solver_->isProvenOptimal() &&
531                !solver_->isDualObjectiveLimitReached()) ;
532  }
533  if (problemFeasibility_->feasible(this,0)<0) {
534    feasible=false; // pretend infeasible
535  }
536/*
537  If the linear relaxation of the root is infeasible, bail out now. Otherwise,
538  continue with processing the root node.
539*/
540  if (!feasible)
541  { handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
542    status_ = 0 ;
543    secondaryStatus_ = 1;
544    originalContinuousObjective_ = COIN_DBL_MAX;
545    solverCharacteristics_ = NULL;
546    return ; }
547  // Save objective (just so user can access it)
548  originalContinuousObjective_ = solver_->getObjValue();
549  bestPossibleObjective_=originalContinuousObjective_;
550  sumChangeObjective1_=0.0;
551  sumChangeObjective2_=0.0;
552/*
553  OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
554  Assuming it recognises the problem, when called upon it will check a cut to
555  see if it cuts off the optimal answer.
556*/
557  // If debugger exists set specialOptions_ bit
558  if (solver_->getRowCutDebuggerAlways())
559    specialOptions_ |= 1;
560
561# ifdef CBC_DEBUG
562  if ((specialOptions_&1)==0)
563    solver_->activateRowCutDebugger(problemName.c_str()) ;
564  if (solver_->getRowCutDebuggerAlways())
565    specialOptions_ |= 1;
566# endif
567
568/*
569  Begin setup to process a feasible root node.
570*/
571  bestObjective_ = CoinMin(bestObjective_,1.0e50) ;
572  numberSolutions_ = 0 ;
573  stateOfSearch_=0;
574  numberHeuristicSolutions_ = 0 ;
575  // Everything is minimization
576  { 
577    // needed to sync cutoffs
578    double value ;
579    solver_->getDblParam(OsiDualObjectiveLimit,value) ;
580    dblParam_[CbcCurrentCutoff]= value * solver_->getObjSense();
581  }
582  double cutoff=getCutoff() ;
583  double direction = solver_->getObjSense() ;
584  dblParam_[CbcOptimizationDirection]=direction;
585  if (cutoff < 1.0e20&&direction<0.0)
586    messageHandler()->message(CBC_CUTOFF_WARNING1,
587                                    messages())
588                                      << cutoff << -cutoff << CoinMessageEol ;
589  if (cutoff > bestObjective_)
590    cutoff = bestObjective_ ;
591  setCutoff(cutoff) ;
592/*
593  We probably already have a current solution, but just in case ...
594*/
595  int numberColumns = getNumCols() ;
596  if (!currentSolution_)
597    currentSolution_ = new double[numberColumns] ;
598  testSolution_ = currentSolution_;
599  /* Tell solver we are in Branch and Cut
600     Could use last parameter for subtle differences */
601  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
602/*
603  Create a copy of the solver, thus capturing the original (root node)
604  constraint system (aka the continuous system).
605*/
606  continuousSolver_ = solver_->clone() ;
607#ifdef COIN_USE_CLP
608  OsiClpSolverInterface * clpSolver
609    = dynamic_cast<OsiClpSolverInterface *> (solver_);
610  if (clpSolver) {
611    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
612    // take off names
613    clpSimplex->dropNames();
614  }
615#endif
616
617  numberRowsAtContinuous_ = getNumRows() ;
618/*
619  Check the objective to see if we can deduce a nontrivial increment. If
620  it's better than the current value for CbcCutoffIncrement, it'll be
621  installed.
622*/
623  if(solverCharacteristics_->reducedCostsAccurate())
624    analyzeObjective() ;
625/*
626  Set up for cut generation. addedCuts_ holds the cuts which are relevant for
627  the active subproblem. whichGenerator will be used to record the generator
628  that produced a given cut.
629*/
630  maximumWhich_ = 1000 ;
631  delete [] whichGenerator_;
632  whichGenerator_ = new int[maximumWhich_] ;
633  memset(whichGenerator_,0,maximumWhich_*sizeof(int));
634  int currentNumberCuts = 0 ;
635  maximumNumberCuts_ = 0 ;
636  currentNumberCuts_ = 0 ;
637  delete [] addedCuts_ ;
638  addedCuts_ = NULL ;
639/*
640  Set up an empty heap and associated data structures to hold the live set
641  (problems which require further exploration).
642*/
643  tree_->setComparison(*nodeCompare_) ;
644/*
645  Used to record the path from a node to the root of the search tree, so that
646  we can then traverse from the root to the node when restoring a subproblem.
647*/
648  maximumDepth_ = 10 ;
649  delete [] walkback_ ;
650  walkback_ = new CbcNodeInfo * [maximumDepth_] ;
651/*
652  Used to generate bound edits for CbcPartialNodeInfo.
653*/
654  double * lowerBefore = new double [numberColumns] ;
655  double * upperBefore = new double [numberColumns] ;
656/*
657 
658  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
659  lifting. It will iterate a generate/reoptimise loop (including reduced cost
660  fixing) until no cuts are generated, the change in objective falls off,  or
661  the limit on the number of rounds of cut generation is exceeded.
662
663  At the end of all this, any cuts will be recorded in cuts and also
664  installed in the solver's constraint system. We'll have reoptimised, and
665  removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
666  adjusted accordingly).
667
668  Tell cut generators they can be a bit more aggressive at root node
669
670  TODO: Why don't we make a copy of the solution after solveWithCuts?
671  TODO: If numberUnsatisfied == 0, don't we have a solution?
672*/
673  phase_=1;
674  int iCutGenerator;
675  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
676    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
677    generator->setAggressiveness(generator->getAggressiveness()+100);
678  }
679  OsiCuts cuts ;
680  int anyAction = -1 ;
681  numberOldActiveCuts_ = 0 ;
682  numberNewCuts_ = 0 ;
683  // Array to mark solution
684  delete [] usedInSolution_;
685  usedInSolution_ = new int[numberColumns];
686  CoinZeroN(usedInSolution_,numberColumns);
687/*
688  For printing totals and for CbcNode (numberNodes_)
689*/
690  numberIterations_ = 0 ;
691  numberNodes_ = 0 ;
692  numberNodes2_ = 0 ;
693  maximumStatistics_=0;
694  statistics_ = NULL;
695  // Do on switch
696  if (doStatistics) {
697    maximumStatistics_=10000;
698    statistics_ = new CbcStatistics * [maximumStatistics_];
699    memset(statistics_,0,maximumStatistics_*sizeof(CbcStatistics *));
700  }
701
702  { int iObject ;
703    int preferredWay ;
704    int numberUnsatisfied = 0 ;
705    memcpy(currentSolution_,solver_->getColSolution(),
706           numberColumns*sizeof(double)) ;
707
708    for (iObject = 0 ; iObject < numberObjects_ ; iObject++)
709    { double infeasibility =
710          object_[iObject]->infeasibility(preferredWay) ;
711      if (infeasibility ) numberUnsatisfied++ ; }
712    // replace solverType
713    if(solverCharacteristics_->tryCuts())  {
714      if (numberUnsatisfied)   {
715        feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
716                                 NULL);
717      } else if (solverCharacteristics_->solutionAddsCuts()) {
718        // may generate cuts and turn the solution
719        //to an infeasible one
720        feasible = solveWithCuts(cuts, 1,
721                                 NULL);
722#if 0
723        currentNumberCuts_ = cuts.sizeRowCuts();
724        if (currentNumberCuts_ >= maximumNumberCuts_) {
725          maximumNumberCuts_ = currentNumberCuts;
726          delete [] addedCuts_;
727          addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
728        }
729#endif
730      }
731    }
732    // check extra info on feasibility
733    if (!solverCharacteristics_->mipFeasible())
734      feasible = false;
735  }
736  // make cut generators less aggressive
737  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
738    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
739    generator->setAggressiveness(generator->getAggressiveness()-100);
740  }
741  currentNumberCuts_ = numberNewCuts_ ;
742/*
743  We've taken the continuous relaxation as far as we can. Time to branch.
744  The first order of business is to actually create a node. chooseBranch
745  currently uses strong branching to evaluate branch object candidates,
746  unless forced back to simple branching. If chooseBranch concludes that a
747  branching candidate is monotone (anyAction == -1) or infeasible (anyAction
748  == -2) when forced to integer values, it returns here immediately.
749
750  Monotone variables trigger a call to resolve(). If the problem remains
751  feasible, try again to choose a branching variable. At the end of the loop,
752  resolved == true indicates that some variables were fixed.
753
754  Loss of feasibility will result in the deletion of newNode.
755*/
756
757  bool resolved = false ;
758  CbcNode *newNode = NULL ;
759  numberFixedAtRoot_=0;
760  numberFixedNow_=0;
761  int numberIterationsAtContinuous = numberIterations_;
762  //solverCharacteristics_->setSolver(solver_);
763  if (feasible) {
764    newNode = new CbcNode ;
765    // Set objective value (not so obvious if NLP etc)
766    setObjectiveValue(newNode,NULL);
767    anyAction = -1 ;
768    // To make depth available we may need a fake node
769    CbcNode fakeNode;
770    if (!currentNode_) {
771      // Not true if sub trees assert (!numberNodes_);
772      currentNode_=&fakeNode;
773    }
774    phase_=3;
775    // only allow 1000 passes
776    int numberPassesLeft=1000;
777    // This is first crude step
778    if (numberAnalyzeIterations_) {
779      delete [] analyzeResults_;
780      analyzeResults_ = new double [4*numberIntegers_];
781      numberFixedAtRoot_=newNode->analyze(this,analyzeResults_);
782      if (numberFixedAtRoot_>0) {
783        printf("%d fixed by analysis\n",numberFixedAtRoot_);
784        setPointers(solver_);
785        numberFixedNow_ = numberFixedAtRoot_;
786      } else if (numberFixedAtRoot_<0) {
787        printf("analysis found to be infeasible\n");
788        anyAction=-2;
789        delete newNode ;
790        newNode = NULL ;
791        feasible = false ;
792      }
793    }
794    while (anyAction == -1)
795    {
796      // Set objective value (not so obvious if NLP etc)
797      setObjectiveValue(newNode,NULL);
798      if (numberBeforeTrust_==0 ) {
799        anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ;
800      } else {
801        OsiSolverBranch * branches=NULL;
802        anyAction = newNode->chooseDynamicBranch(this,NULL,branches,numberPassesLeft) ;
803        if (anyAction==-3) 
804          anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ; // dynamic did nothing
805      }
806      numberPassesLeft--;
807      if (anyAction == -1)
808      { feasible = resolve() ;
809      if (problemFeasibility_->feasible(this,0)<0) {
810        feasible=false; // pretend infeasible
811      }
812      reducedCostFix() ;
813      resolved = true ;
814#       ifdef CBC_DEBUG
815      printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
816             solver_->getObjValue(),
817             solver_->getNumRows()) ;
818#       endif
819      if (!feasible) anyAction = -2 ; }
820      if (anyAction == -2||newNode->objectiveValue() >= cutoff)
821      { delete newNode ;
822        newNode = NULL ;
823        feasible = false ; } } }
824/*
825  At this point, the root subproblem is infeasible or fathomed by bound
826  (newNode == NULL), or we're live with an objective value that satisfies the
827  current objective cutoff.
828*/
829  assert (!newNode || newNode->objectiveValue() <= cutoff) ;
830  // Save address of root node as we don't want to delete it
831  CbcNode * rootNode = newNode;
832/*
833  The common case is that the lp relaxation is feasible but doesn't satisfy
834  integrality (i.e., newNode->variable() >= 0, indicating we've been able to
835  select a branching variable). Remove any cuts that have gone slack due to
836  forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
837  basis (via createInfo()) and stash the new cuts in the nodeInfo (via
838  addCuts()). If, by some miracle, we have an integral solution at the root
839  (newNode->variable() < 0), takeOffCuts() will ensure that the solver holds
840  a valid solution for use by setBestSolution().
841*/
842  CoinWarmStartBasis *lastws = 0 ;
843  if (feasible && newNode->variable() >= 0)
844  { if (resolved)
845    { bool needValidSolution = (newNode->variable() < 0) ;
846      takeOffCuts(cuts,needValidSolution,NULL) ;
847#     ifdef CHECK_CUT_COUNTS
848      { printf("Number of rows after chooseBranch fix (root)"
849               "(active only) %d\n",
850                numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ;
851        const CoinWarmStartBasis* debugws =
852          dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
853        debugws->print() ;
854        delete debugws ; }
855#     endif
856    }
857    newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
858    newNode->nodeInfo()->addCuts(cuts,
859                                 newNode->numberBranches(),whichGenerator_) ;
860    if (lastws) delete lastws ;
861    lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
862  }
863/*
864  Continuous data to be used later
865*/
866  continuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
867  continuousInfeasibilities_ = 0 ;
868  if (newNode)
869  { continuousObjective_ = newNode->objectiveValue() ;
870    delete [] continuousSolution_;
871    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
872                                             numberColumns);
873    continuousInfeasibilities_ = newNode->numberUnsatisfied() ; }
874/*
875  Bound may have changed so reset in objects
876*/
877  { int i ;
878    for (i = 0;i < numberObjects_;i++)
879      object_[i]->resetBounds() ; }
880  stoppedOnGap_ = false ;
881/*
882  Feasible? Then we should have either a live node prepped for future
883  expansion (indicated by variable() >= 0), or (miracle of miracles) an
884  integral solution at the root node.
885
886  initializeInfo sets the reference counts in the nodeInfo object.  Since
887  this node is still live, push it onto the heap that holds the live set.
888*/
889  double bestValue = 0.0 ;
890  if (newNode) {
891    bestValue = newNode->objectiveValue();
892    if (newNode->variable() >= 0) {
893      newNode->initializeInfo() ;
894      tree_->push(newNode) ;
895      if (statistics_) {
896        if (numberNodes2_==maximumStatistics_) {
897          maximumStatistics_ = 2*maximumStatistics_;
898          CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
899          memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
900          memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
901          delete [] statistics_;
902          statistics_=temp;
903        }
904        assert (!statistics_[numberNodes2_]);
905        statistics_[numberNodes2_]=new CbcStatistics(newNode);
906      }
907      numberNodes2_++;
908#     ifdef CHECK_NODE
909      printf("Node %x on tree\n",newNode) ;
910#     endif
911    } else {
912      // continuous is integer
913      double objectiveValue = newNode->objectiveValue();
914      setBestSolution(CBC_SOLUTION,objectiveValue,
915                      solver_->getColSolution()) ;
916      delete newNode ;
917      newNode = NULL ;
918    }
919  }
920
921  if (printFrequency_ <= 0) {
922    printFrequency_ = 1000 ;
923    if (getNumCols() > 2000)
924      printFrequency_ = 100 ;
925  }
926  /*
927    It is possible that strong branching fixes one variable and then the code goes round
928    again and again.  This can take too long.  So we need to warn user - just once.
929  */
930  numberLongStrong_=0;
931/*
932  At last, the actual branch-and-cut search loop, which will iterate until
933  the live set is empty or we hit some limit (integrality gap, time, node
934  count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
935  solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
936  add the node to the live set.
937
938  The first action is to winnow the live set to remove nodes which are worse
939  than the current objective cutoff.
940*/
941  while (!tree_->empty()) {
942    if (cutoff > getCutoff()) {
943      double newCutoff = getCutoff();
944      if (analyzeResults_) {
945        // see if we could fix any (more)
946        int n=0;
947        double * newLower = analyzeResults_;
948        double * objLower = newLower+numberIntegers_;
949        double * newUpper = objLower+numberIntegers_;
950        double * objUpper = newUpper+numberIntegers_;
951        for (int i=0;i<numberIntegers_;i++) {
952          if (objLower[i]>newCutoff) {
953            n++;
954            if (objUpper[i]>newCutoff) {
955              newCutoff = -COIN_DBL_MAX;
956              break;
957            }
958          } else if (objUpper[i]>newCutoff) {
959            n++;
960          }
961        }
962        if (newCutoff==-COIN_DBL_MAX) {
963          printf("Root analysis says finished\n");
964        } else if (n>numberFixedNow_) {
965          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
966          numberFixedNow_=n;
967        }
968      }
969#ifdef COIN_USE_CLP
970      if (eventHandler) {
971        if (!eventHandler->event(ClpEventHandler::solution)) {
972          eventHappened_=true; // exit
973        }
974      }
975#endif
976      // Do from deepest
977      tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ;
978      nodeCompare_->newSolution(this) ;
979      nodeCompare_->newSolution(this,continuousObjective_,
980                                continuousInfeasibilities_) ;
981      tree_->setComparison(*nodeCompare_) ;
982      if (tree_->empty())
983        break; // finished
984    }
985    cutoff = getCutoff() ;
986/*
987  Periodic activities: Opportunities to
988    + tweak the nodeCompare criteria,
989    + check if we've closed the integrality gap enough to quit,
990    + print a summary line to let the user know we're working
991*/
992    if ((numberNodes_%1000) == 0) {
993      bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ;
994      // redo tree if wanted
995      if (redoTree)
996        tree_->setComparison(*nodeCompare_) ;
997    }
998    if (saveCompare&&!hotstartSolution_) {
999      // hotstart switched off
1000      delete nodeCompare_; // off depth first
1001      nodeCompare_=saveCompare;
1002      saveCompare=NULL;
1003      // redo tree
1004      tree_->setComparison(*nodeCompare_) ;
1005    }
1006    if ((numberNodes_%printFrequency_) == 0) {
1007      int j ;
1008      int nNodes = tree_->size() ;
1009      bestPossibleObjective_ = 1.0e100 ;
1010      for (j = 0;j < nNodes;j++) {
1011        CbcNode * node = tree_->nodePointer(j) ;
1012        if (node&&node->objectiveValue() < bestPossibleObjective_)
1013          bestPossibleObjective_ = node->objectiveValue() ;
1014      }
1015      messageHandler()->message(CBC_STATUS,messages())
1016        << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
1017        << CoinMessageEol ;
1018#ifdef COIN_USE_CLP
1019      if (eventHandler) {
1020        if (!eventHandler->event(ClpEventHandler::treeStatus)) {
1021          eventHappened_=true; // exit
1022        }
1023      }
1024#endif
1025    }
1026    // If no solution but many nodes - signal change in strategy
1027    if (numberNodes_>2*numberObjects_+1000&&stateOfSearch_!=2)
1028      stateOfSearch_=3;
1029    // See if can stop on gap
1030    double testGap = CoinMax(dblParam_[CbcAllowableGap],
1031                             CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_))
1032                             *dblParam_[CbcAllowableFractionGap]);
1033    if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
1034      stoppedOnGap_ = true ;
1035    }
1036
1037#   ifdef CHECK_NODE_FULL
1038    verifyTreeNodes(tree_,*this) ;
1039#   endif
1040#   ifdef CHECK_CUT_COUNTS
1041    verifyCutCounts(tree_,*this) ;
1042#   endif
1043
1044/*
1045  Now we come to the meat of the loop. To create the active subproblem, we'll
1046  pop the most promising node in the live set, rebuild the subproblem it
1047  represents, and then execute the current arm of the branch to create the
1048  active subproblem.
1049*/
1050    CbcNode *node = tree_->bestNode(cutoff) ;
1051    // Possible one on tree worse than cutoff
1052    if (!node)
1053      continue;
1054    currentNode_=node; // so can be accessed elsewhere
1055#ifdef CBC_DEBUG
1056    printf("%d unsat, way %d, obj %g est %g\n",
1057           node->numberUnsatisfied(),node->way(),node->objectiveValue(),
1058           node->guessedObjectiveValue());
1059#endif
1060    // Save clone in branching decision
1061    if(branchingMethod_)
1062      branchingMethod_->saveBranchingObject(node->modifiableBranchingObject());
1063    bool nodeOnTree=false; // Node has been popped
1064    // Say not on optimal path
1065    bool onOptimalPath=false;
1066#   ifdef CHECK_NODE
1067/*
1068  WARNING: The use of integerVariable_[*] here will break as soon as the
1069           branching object is something other than an integer variable.
1070           This needs some thought.
1071*/
1072    printf("Node %x popped from tree - %d left, %d count\n",node,
1073           node->nodeInfo()->numberBranchesLeft(),
1074           node->nodeInfo()->numberPointingToThis()) ;
1075    printf("\tdepth = %d, z =  %g, unsat = %d, var = %d.\n",
1076           node->depth(),node->objectiveValue(),
1077           node->numberUnsatisfied(),
1078           integerVariable_[node->variable()]) ;
1079#   endif
1080
1081/*
1082  Rebuild the subproblem for this node:  Call addCuts() to adjust the model
1083  to recreate the subproblem for this node (set proper variable bounds, add
1084  cuts, create a basis).  This may result in the problem being fathomed by
1085  bound or infeasibility. Returns 1 if node is fathomed.
1086  Execute the current arm of the branch: If the problem survives, save the
1087  resulting variable bounds and call branch() to modify variable bounds
1088  according to the current arm of the branching object. If we're processing
1089  the final arm of the branching object, flag the node for removal from the
1090  live set.
1091*/
1092    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
1093    newNode = NULL ;
1094    if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_))
1095    { int i ;
1096      const double * lower = getColLower() ;
1097      const double * upper = getColUpper() ;
1098      for (i = 0 ; i < numberColumns ; i++)
1099      { lowerBefore[i]= lower[i] ;
1100        upperBefore[i]= upper[i] ; }
1101      bool deleteNode ;
1102      if (node->branch())
1103      { 
1104        // set nodenumber correctly
1105        node->nodeInfo()->setNodeNumber(numberNodes2_);
1106        tree_->push(node) ;
1107        if (statistics_) {
1108          if (numberNodes2_==maximumStatistics_) {
1109            maximumStatistics_ = 2*maximumStatistics_;
1110            CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
1111            memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
1112            memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
1113            delete [] statistics_;
1114            statistics_=temp;
1115          }
1116          assert (!statistics_[numberNodes2_]);
1117          statistics_[numberNodes2_]=new CbcStatistics(node);
1118        }
1119        numberNodes2_++;
1120        nodeOnTree=true; // back on tree
1121        deleteNode = false ;
1122#       ifdef CHECK_NODE
1123        printf("Node %x pushed back on tree - %d left, %d count\n",node,
1124               nodeInfo->numberBranchesLeft(),
1125               nodeInfo->numberPointingToThis()) ;
1126#       endif
1127      }
1128      else
1129      { deleteNode = true ; }
1130
1131      if ((specialOptions_&1)!=0) {
1132        /*
1133          This doesn't work as intended --- getRowCutDebugger will return null
1134          unless the current feasible solution region includes the optimal solution
1135          that RowCutDebugger knows. There's no way to tell inactive from off the
1136          optimal path.
1137        */
1138        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
1139        if (debugger) {
1140          onOptimalPath=true;
1141          printf("On optimal path\n") ;
1142        }
1143      }
1144/*
1145  Reoptimize, possibly generating cuts and/or using heuristics to find
1146  solutions.  Cut reference counts are unaffected unless we lose feasibility,
1147  in which case solveWithCuts() will make the adjustment.
1148*/
1149      phase_=2;
1150      cuts = OsiCuts() ;
1151      currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_ ;
1152      int saveNumber = numberIterations_;
1153      if(solverCharacteristics_->solutionAddsCuts()) {
1154        feasible=resolve();
1155        if (feasible) {
1156          int iObject ;
1157          int preferredWay ;
1158          int numberUnsatisfied = 0 ;
1159          memcpy(currentSolution_,solver_->getColSolution(),
1160                 numberColumns*sizeof(double)) ;
1161         
1162          for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
1163            double infeasibility =
1164              object_[iObject]->infeasibility(preferredWay) ;
1165            if (infeasibility ) numberUnsatisfied++ ;
1166          }
1167          if (numberUnsatisfied)   {
1168            feasible = solveWithCuts(cuts,maximumCutPasses_,node);
1169          } else {
1170            // may generate cuts and turn the solution
1171            //to an infeasible one
1172            feasible = solveWithCuts(cuts, 1,
1173                                     node);
1174#if 0
1175            currentNumberCuts_ = cuts.sizeRowCuts();
1176            if (currentNumberCuts_ >= maximumNumberCuts_) {
1177              maximumNumberCuts_ = currentNumberCuts;
1178              delete [] addedCuts_;
1179              addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
1180            }
1181#endif
1182          }
1183          // check extra info on feasibility
1184          if (!solverCharacteristics_->mipFeasible())
1185            feasible = false;
1186        }
1187      } else {
1188        // normal
1189        feasible = solveWithCuts(cuts,maximumCutPasses_,node);
1190      }
1191      if ((specialOptions_&1)!=0&&onOptimalPath) {
1192        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
1193        assert (debugger) ;
1194      }
1195      if (statistics_) {
1196        assert (numberNodes2_);
1197        assert (statistics_[numberNodes2_-1]);
1198        assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
1199        statistics_[numberNodes2_-1]->endOfBranch(numberIterations_-saveNumber,
1200                                               feasible ? solver_->getObjValue()
1201                                               : COIN_DBL_MAX);
1202    }
1203/*
1204  Check for abort on limits: node count, solution count, time, integrality gap.
1205*/
1206      double totalTime = CoinCpuTime()-dblParam_[CbcStartSeconds] ;
1207      if (numberNodes_ < intParam_[CbcMaxNumNode] &&
1208          numberSolutions_ < intParam_[CbcMaxNumSol] &&
1209          totalTime < dblParam_[CbcMaximumSeconds] &&
1210          !stoppedOnGap_&&!eventHappened_) 
1211      {
1212/*
1213  Are we still feasible? If so, create a node and do the work to attach a
1214  branching object, reoptimising as needed if chooseBranch() identifies
1215  monotone objects.
1216
1217  Finally, attach a partial nodeInfo object and store away any cuts that we
1218  created back in solveWithCuts. addCuts() will also deal with the cut
1219  reference counts.
1220
1221  TODO: (lh) I'm confused. We create a nodeInfo without checking whether we
1222        have a solution or not. Then we use numberUnsatisfied() to decide
1223        whether to stash the cuts and bump reference counts. Other places we
1224        use variable() (i.e., presence of a branching variable). Equivalent?
1225*/
1226        if (onOptimalPath) {
1227          if (!feasible) {
1228            printf("infeas2\n");
1229            solver_->writeMps("infeas");
1230            CoinWarmStartBasis *slack =
1231              dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
1232            solver_->setWarmStart(slack);
1233            delete slack ;
1234            solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
1235            solver_->initialSolve();
1236            assert (!solver_->isProvenOptimal());
1237          }
1238          assert (feasible);
1239        }
1240        bool checkingNode=false;
1241        if (feasible) {
1242          newNode = new CbcNode ;
1243          // Set objective value (not so obvious if NLP etc)
1244          setObjectiveValue(newNode,node);
1245          anyAction =-1 ;
1246          resolved = false ;
1247          if (newNode->objectiveValue() >= getCutoff()) 
1248            anyAction=-2;
1249          // only allow twenty passes
1250          int numberPassesLeft=20;
1251          checkingNode=true;
1252          while (anyAction == -1)
1253          { 
1254            // Set objective value (not so obvious if NLP etc)
1255            setObjectiveValue(newNode,node);
1256            if (numberBeforeTrust_==0 ) {
1257              anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ;
1258            } else {
1259              OsiSolverBranch * branches=NULL;
1260              anyAction = newNode->chooseDynamicBranch(this,node,branches,numberPassesLeft) ;
1261              if (anyAction==-3) 
1262                anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ; // dynamic did nothing
1263            }
1264            if (onOptimalPath)
1265              assert (anyAction!=-2); // can be useful but gives false positives on strong
1266            numberPassesLeft--;
1267            if (numberPassesLeft<=-1) {
1268              if (!numberLongStrong_)
1269                messageHandler()->message(CBC_WARNING_STRONG,
1270                                          messages()) << CoinMessageEol ;
1271              numberLongStrong_++;
1272            }
1273            if (anyAction == -1)
1274            {
1275              // can do quick optimality check
1276              int easy=2;
1277              solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
1278              feasible = resolve() ;
1279              solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1280              resolved = true ;
1281              if (problemFeasibility_->feasible(this,0)<0) {
1282                feasible=false; // pretend infeasible
1283              }
1284              if (feasible)
1285              { 
1286                // Set objective value (not so obvious if NLP etc)
1287                setObjectiveValue(newNode,node);
1288                reducedCostFix() ;
1289                if (newNode->objectiveValue() >= getCutoff()) 
1290                  anyAction=-2;
1291              }
1292              else
1293              { anyAction = -2 ; } } }
1294          if (anyAction >= 0)
1295          { if (resolved)
1296            { bool needValidSolution = (newNode->variable() < 0) ;
1297              takeOffCuts(cuts,needValidSolution,NULL) ; 
1298#             ifdef CHECK_CUT_COUNTS
1299              { printf("Number of rows after chooseBranch fix (node)"
1300                       "(active only) %d\n",
1301                        numberRowsAtContinuous_+numberNewCuts_+
1302                        numberOldActiveCuts_) ;
1303                const CoinWarmStartBasis* debugws =
1304                  dynamic_cast<const CoinWarmStartBasis*>
1305                    (solver_->getWarmStart()) ;
1306                debugws->print() ;
1307                delete debugws ; }
1308#             endif
1309            }
1310            newNode->createInfo(this,node,lastws,lowerBefore,upperBefore,
1311                                numberOldActiveCuts_,numberNewCuts_) ;
1312            if (newNode->numberUnsatisfied()) {
1313              newNode->initializeInfo() ;
1314              newNode->nodeInfo()->addCuts(cuts,newNode->numberBranches(),
1315                                           whichGenerator_) ; } } }
1316        else
1317        { anyAction = -2 ; }
1318        // May have slipped through i.e. anyAction == 0 and objective above cutoff
1319        if ( anyAction >=0 ) {
1320          assert (newNode);
1321          if (newNode->objectiveValue() >= getCutoff()) 
1322            anyAction = -2; // say bad after all
1323        }
1324/*
1325  If we end up infeasible, we can delete the new node immediately. Since this
1326  node won't be needing the cuts we collected, decrement the reference counts.
1327  If we are feasible, then we'll be placing this node into the live set, so
1328  increment the reference count in the current (parent) nodeInfo.
1329*/
1330        if (anyAction == -2)
1331        { delete newNode ;
1332          newNode = NULL ;
1333          // say strong doing well
1334          if (checkingNode)
1335            setSpecialOptions(specialOptions_|8);
1336          for (i = 0 ; i < currentNumberCuts_ ; i++)
1337          { if (addedCuts_[i])
1338            { if (!addedCuts_[i]->decrement(1))
1339                delete addedCuts_[i] ; } } }
1340        else
1341        { nodeInfo->increment() ;
1342        if ((numberNodes_%20)==0) {
1343          // say strong not doing as well
1344          setSpecialOptions(specialOptions_&~8);
1345        }
1346        }
1347/*
1348  At this point, there are three possibilities:
1349    * We have a live node (variable() >= 0) which will require further
1350      branching to resolve. Before we push it onto the search tree, try for
1351      a heuristic solution.
1352    * We have a solution, in which case newNode is non-null but we have no
1353      branching variable. Decrement the cut counts and save the solution.
1354    * The node was found to be infeasible, in which case it's already been
1355      deleted, and newNode is null.
1356
1357  TODO: (lh) Now I'm more confused. I thought that the call to addCuts() above
1358        took care of incrementing the reference counts for cuts at newNode.
1359        Clearly I need to look more carefully.
1360*/
1361#ifdef COIN_USE_CLP
1362        if (eventHandler) {
1363          if (!eventHandler->event(ClpEventHandler::node)) {
1364            eventHappened_=true; // exit
1365          }
1366        }
1367#endif
1368        assert (!newNode || newNode->objectiveValue() <= getCutoff()) ;
1369        if (statistics_) {
1370          assert (numberNodes2_);
1371          assert (statistics_[numberNodes2_-1]);
1372          assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
1373          if (newNode)
1374            statistics_[numberNodes2_-1]->updateInfeasibility(newNode->numberUnsatisfied());
1375          else
1376            statistics_[numberNodes2_-1]->sayInfeasible();
1377        }
1378        if (newNode)
1379        { if (newNode->variable() >= 0)
1380          { handler_->message(CBC_BRANCH,messages_)
1381               << numberNodes_<< newNode->objectiveValue()
1382               << newNode->numberUnsatisfied()<< newNode->depth()
1383               << CoinMessageEol ;
1384            // Increment cut counts (taking off current)
1385            int numberLeft = newNode->numberBranches() ;
1386            for (i = 0;i < currentNumberCuts_;i++)
1387            { if (addedCuts_[i])
1388              {
1389#               ifdef CHECK_CUT_COUNTS
1390                printf("Count on cut %x increased by %d\n",addedCuts_[i],
1391                        numberLeft-1) ;
1392#               endif
1393                addedCuts_[i]->increment(numberLeft-1) ; } }
1394
1395            double estValue = newNode->guessedObjectiveValue() ;
1396            int found = -1 ;
1397            // no - overhead on small problems solver_->resolve() ;     // double check current optimal
1398            // assert (!solver_->getIterationCount());
1399            double * newSolution = new double [numberColumns] ;
1400            double heurValue = getCutoff() ;
1401            int iHeur ;
1402            for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++)
1403            { double saveValue = heurValue ;
1404              int ifSol = heuristic_[iHeur]->solution(heurValue,newSolution) ;
1405              if (ifSol > 0) {
1406                // new solution found
1407                found = iHeur ;
1408                incrementUsed(newSolution);
1409              }
1410              else
1411              if (ifSol < 0)    // just returning an estimate
1412              { estValue = CoinMin(heurValue,estValue) ;
1413                heurValue = saveValue ; } }
1414            if (found >= 0) {
1415              setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
1416              lastHeuristic_ = heuristic_[found];
1417            }
1418            delete [] newSolution ;
1419            newNode->setGuessedObjectiveValue(estValue) ;
1420            tree_->push(newNode) ;
1421            if (statistics_) {
1422              if (numberNodes2_==maximumStatistics_) {
1423                maximumStatistics_ = 2*maximumStatistics_;
1424                CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
1425                memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
1426                memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
1427                delete [] statistics_;
1428                statistics_=temp;
1429              }
1430              assert (!statistics_[numberNodes2_]);
1431              statistics_[numberNodes2_]=new CbcStatistics(newNode);
1432            }
1433            numberNodes2_++;
1434#           ifdef CHECK_NODE
1435            printf("Node %x pushed on tree c\n",newNode) ;
1436#           endif
1437          }
1438          else
1439          { for (i = 0 ; i < currentNumberCuts_ ; i++)
1440            { if (addedCuts_[i])
1441              { if (!addedCuts_[i]->decrement(1))
1442                  delete addedCuts_[i] ; } }
1443          double objectiveValue = newNode->objectiveValue();
1444            setBestSolution(CBC_SOLUTION,objectiveValue,
1445                            solver_->getColSolution()) ;
1446            lastHeuristic_ = NULL;
1447            incrementUsed(solver_->getColSolution());
1448            //assert(nodeInfo->numberPointingToThis() <= 2) ;
1449            // avoid accidental pruning, if newNode was final branch arm
1450            nodeInfo->increment();
1451            delete newNode ;
1452            nodeInfo->decrement() ; } }
1453/*
1454  This node has been completely expanded and can be removed from the live
1455  set.
1456*/
1457        if (deleteNode) 
1458          delete node ; }
1459/*
1460  End of the non-abort actions. The next block of code is executed if we've
1461  aborted because we hit one of the limits. Clean up by deleting the live set
1462  and break out of the node processing loop.
1463*/
1464      else
1465      { 
1466        tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
1467        delete nextRowCut_;
1468        // We need to get rid of node if is has already been popped from tree
1469        if (!nodeOnTree&&!stoppedOnGap_&&node!=rootNode)
1470          delete node;
1471        if (stoppedOnGap_)
1472        { messageHandler()->message(CBC_GAP,messages())
1473            << bestObjective_-bestPossibleObjective_
1474            << dblParam_[CbcAllowableGap]
1475            << dblParam_[CbcAllowableFractionGap]*100.0
1476            << CoinMessageEol ;
1477        secondaryStatus_ = 2;
1478          status_ = 0 ; }
1479        else
1480        if (isNodeLimitReached())
1481        { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
1482        secondaryStatus_ = 3;
1483          status_ = 1 ; }
1484        else
1485        if (totalTime >= dblParam_[CbcMaximumSeconds])
1486        { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
1487        secondaryStatus_ = 4;
1488          status_ = 1 ; }
1489        else
1490        if (eventHappened_)
1491        { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
1492        secondaryStatus_ = 5;
1493          status_ = 5 ; }
1494        else
1495        { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
1496        secondaryStatus_ = 6;
1497          status_ = 1 ; }
1498        break ; }
1499/*
1500  Delete cuts to get back to the original system.
1501
1502  I'm thinking this is redundant --- the call to addCuts that conditions entry
1503  to this code block also performs this action.
1504*/
1505      int numberToDelete = getNumRows()-numberRowsAtContinuous_ ;
1506      if (numberToDelete)
1507      { int * delRows = new int[numberToDelete] ;
1508        int i ;
1509        for (i = 0 ; i < numberToDelete ; i++)
1510        { delRows[i] = i+numberRowsAtContinuous_ ; }
1511        solver_->deleteRows(numberToDelete,delRows) ;
1512        delete [] delRows ; } }
1513/*
1514  This node fathomed when addCuts atttempted to revive it. Toss it.
1515*/
1516    else
1517      { delete node ; } }
1518/*
1519  That's it, we've exhausted the search tree, or broken out of the loop because
1520  we hit some limit on evaluation.
1521
1522  We may have got an intelligent tree so give it one more chance
1523*/
1524  // Tell solver we are not in Branch and Cut
1525  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
1526  tree_->endSearch();
1527  //  If we did any sub trees - did we give up on any?
1528  if ( numberStoppedSubTrees_)
1529    status_=1;
1530  if (!status_) {
1531    bestPossibleObjective_=bestObjective_;
1532    handler_->message(CBC_END_GOOD,messages_)
1533      << bestObjective_ << numberIterations_ << numberNodes_
1534      << CoinMessageEol ;
1535  } else {
1536    handler_->message(CBC_END,messages_)
1537      << bestObjective_ <<bestPossibleObjective_
1538      << numberIterations_ << numberNodes_
1539      << CoinMessageEol ;
1540  }
1541  if (numberStrongIterations_)
1542    handler_->message(CBC_STRONG_STATS,messages_)
1543      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
1544      << strongInfo_[1] << CoinMessageEol ;
1545  if (statistics_) {
1546    // report in some way
1547    int * lookup = new int[numberObjects_];
1548    int i;
1549    for (i=0;i<numberObjects_;i++) 
1550      lookup[i]=-1;
1551    bool goodIds=true;
1552    for (i=0;i<numberObjects_;i++) {
1553      int id = object_[i]->id();
1554      int iColumn = object_[i]->columnNumber();
1555      if (iColumn<0)
1556        iColumn = id+numberColumns;
1557      if(id>=0&&id<numberObjects_) {
1558        if (lookup[id]==-1) {
1559          lookup[id]=iColumn;
1560        } else {
1561          goodIds=false;
1562          break;
1563        }
1564      } else {
1565        goodIds=false;
1566        break;
1567      }
1568    }
1569    if (!goodIds) {
1570      delete [] lookup;
1571      lookup=NULL;
1572    }
1573    if (doStatistics==3) {
1574      printf("  node parent depth column   value                    obj      inf\n");
1575      for ( i=0;i<numberNodes2_;i++) {
1576        statistics_[i]->print(lookup);
1577      }
1578    }
1579    if (doStatistics>1) {
1580      // Find last solution
1581      int k;
1582      for (k=numberNodes2_-1;k>=0;k--) {
1583        if (statistics_[k]->endingObjective()!=COIN_DBL_MAX&&
1584            !statistics_[k]->endingInfeasibility())
1585          break;
1586      }
1587      if (k>=0) {
1588        int depth=statistics_[k]->depth();
1589        int * which = new int[depth+1];
1590        for (i=depth;i>=0;i--) {
1591          which[i]=k;
1592          k=statistics_[k]->parentNode();
1593        }
1594        printf("  node parent depth column   value                    obj      inf\n");
1595        for (i=0;i<=depth;i++) {
1596          statistics_[which[i]]->print(lookup);
1597        }
1598        delete [] which;
1599      }
1600    }
1601    // now summary
1602    int maxDepth=0;
1603    double averageSolutionDepth=0.0;
1604    int numberSolutions=0;
1605    double averageCutoffDepth=0.0;
1606    double averageSolvedDepth=0.0;
1607    int numberCutoff=0;
1608    int numberDown=0;
1609    int numberFirstDown=0;
1610    double averageInfDown=0.0;
1611    double averageObjDown=0.0;
1612    int numberCutoffDown=0;
1613    int numberUp=0;
1614    int numberFirstUp=0;
1615    double averageInfUp=0.0;
1616    double averageObjUp=0.0;
1617    int numberCutoffUp=0;
1618    double averageNumberIterations1=0.0;
1619    double averageValue=0.0;
1620    for ( i=0;i<numberNodes2_;i++) {
1621      int depth =  statistics_[i]->depth(); 
1622      int way =  statistics_[i]->way(); 
1623      double value = statistics_[i]->value(); 
1624      double startingObjective =  statistics_[i]->startingObjective(); 
1625      int startingInfeasibility = statistics_[i]->startingInfeasibility(); 
1626      double endingObjective = statistics_[i]->endingObjective(); 
1627      int endingInfeasibility = statistics_[i]->endingInfeasibility(); 
1628      maxDepth = CoinMax(depth,maxDepth);
1629      // Only for completed
1630      averageNumberIterations1 += statistics_[i]->numberIterations();
1631      averageValue += value;
1632      if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) {
1633        numberSolutions++;
1634        averageSolutionDepth += depth;
1635      }
1636      if (endingObjective==COIN_DBL_MAX) {
1637        numberCutoff++;
1638        averageCutoffDepth += depth;
1639        if (way<0) {
1640          numberDown++;
1641          numberCutoffDown++;
1642          if (way==-1)
1643            numberFirstDown++;
1644        } else {
1645          numberUp++;
1646          numberCutoffUp++;
1647          if (way==1)
1648            numberFirstUp++;
1649        }
1650      } else {
1651        averageSolvedDepth += depth;
1652        if (way<0) {
1653          numberDown++;
1654          averageInfDown += startingInfeasibility-endingInfeasibility;
1655          averageObjDown += endingObjective-startingObjective;
1656          if (way==-1)
1657            numberFirstDown++;
1658        } else {
1659          numberUp++;
1660          averageInfUp += startingInfeasibility-endingInfeasibility;
1661          averageObjUp += endingObjective-startingObjective;
1662          if (way==1)
1663            numberFirstUp++;
1664        }
1665      }
1666    }
1667    // Now print
1668    if (numberSolutions)
1669      averageSolutionDepth /= (double) numberSolutions;
1670    int numberSolved = numberNodes2_-numberCutoff;
1671    double averageNumberIterations2=numberIterations_-averageNumberIterations1
1672      -numberIterationsAtContinuous;
1673    if(numberCutoff) {
1674      averageCutoffDepth /= (double) numberCutoff;
1675      averageNumberIterations2 /= (double) numberCutoff;
1676    }
1677    if (numberNodes2_) 
1678      averageValue /= (double) numberNodes2_;
1679    if (numberSolved) {
1680      averageNumberIterations1 /= (double) numberSolved;
1681      averageSolvedDepth /= (double) numberSolved;
1682    }
1683    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
1684           numberSolutions,averageSolutionDepth);
1685    printf("average value of variable being branched on was %g\n",
1686           averageValue);
1687    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
1688           numberCutoff,averageCutoffDepth,averageNumberIterations2);
1689    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
1690           numberSolved,averageSolvedDepth,averageNumberIterations1);
1691    if (numberDown) {
1692      averageInfDown /= (double) numberDown;
1693      averageObjDown /= (double) numberDown;
1694    }
1695    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
1696           numberDown,numberFirstDown,numberDown-numberFirstDown,numberCutoffDown,
1697           averageInfDown,averageObjDown);
1698    if (numberUp) {
1699      averageInfUp /= (double) numberUp;
1700      averageObjUp /= (double) numberUp;
1701    }
1702    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
1703           numberUp,numberFirstUp,numberUp-numberFirstUp,numberCutoffUp,
1704           averageInfUp,averageObjUp);
1705    for ( i=0;i<numberNodes2_;i++) 
1706      delete statistics_[i];
1707    delete [] statistics_;
1708    statistics_=NULL;
1709    maximumStatistics_=0;
1710    delete [] lookup;
1711  }
1712/*
1713  If we think we have a solution, restore and confirm it with a call to
1714  setBestSolution().  We need to reset the cutoff value so as not to fathom
1715  the solution on bounds.  Note that calling setBestSolution( ..., true)
1716  leaves the continuousSolver_ bounds vectors fixed at the solution value.
1717
1718  Running resolve() here is a failsafe --- setBestSolution has already
1719  reoptimised using the continuousSolver_. If for some reason we fail to
1720  prove optimality, run the problem again after instructing the solver to
1721  tell us more.
1722
1723  If all looks good, replace solver_ with continuousSolver_, so that the
1724  outside world will be able to obtain information about the solution using
1725  public methods.
1726*/
1727  if (bestSolution_&&solverCharacteristics_->solverType()<2)
1728  { setCutoff(1.0e50) ; // As best solution should be worse than cutoff
1729    phase_=5;
1730    setBestSolution(CBC_SOLUTION,bestObjective_,bestSolution_,true) ;
1731    continuousSolver_->resolve() ;
1732    if (!continuousSolver_->isProvenOptimal())
1733    { continuousSolver_->messageHandler()->setLogLevel(2) ;
1734      continuousSolver_->initialSolve() ; }
1735    delete solver_ ;
1736    // above deletes solverCharacteristics_
1737    solverCharacteristics_ = NULL;
1738    solver_ = continuousSolver_ ;
1739    setPointers(solver_);
1740    continuousSolver_ = NULL ; }
1741/*
1742  Clean up dangling objects. continuousSolver_ may already be toast.
1743*/
1744  delete lastws ;
1745  delete [] whichGenerator_ ;
1746  whichGenerator_=NULL;
1747  delete [] lowerBefore ;
1748  delete [] upperBefore ;
1749  delete [] walkback_ ;
1750  walkback_ = NULL ;
1751  delete [] addedCuts_ ;
1752  addedCuts_ = NULL ;
1753  // Get rid of characteristics
1754  solverCharacteristics_=NULL;
1755  if (continuousSolver_)
1756  { delete continuousSolver_ ;
1757    continuousSolver_ = NULL ; }
1758/*
1759  Destroy global cuts by replacing with an empty OsiCuts object.
1760*/
1761  globalCuts_= OsiCuts() ;
1762  if (strategy_&&strategy_->preProcessState()>0) {
1763    // undo preprocessing
1764    CglPreProcess * process = strategy_->process();
1765    assert (process);
1766    int n = originalSolver->getNumCols();
1767    if (bestSolution_) {
1768      delete [] bestSolution_;
1769      bestSolution_ = new double [n];
1770      process->postProcess(*solver_);
1771    }
1772    strategy_->deletePreProcess();
1773    // Solution now back in originalSolver
1774    delete solver_;
1775    solver_=originalSolver;
1776    if (bestSolution_)
1777      memcpy(bestSolution_,solver_->getColSolution(),n*sizeof(double));
1778  }
1779  return ; }
1780#else
1781void CbcModel::branchAndBound(int doStatistics) 
1782
1783{
1784/*
1785  Capture a time stamp before we start.
1786*/
1787  dblParam_[CbcStartSeconds] = CoinCpuTime();
1788  strongInfo_[0]=0;
1789  strongInfo_[1]=0;
1790  strongInfo_[2]=0;
1791  numberStrongIterations_ = 0;
1792  // original solver (only set if pre-processing)
1793  OsiSolverInterface * originalSolver=NULL;
1794  // Set up strategies
1795  if (strategy_) {
1796    // May do preprocessing
1797    originalSolver = solver_;
1798    strategy_->setupOther(*this);
1799    if (strategy_->preProcessState()) {
1800      // pre-processing done
1801      if (strategy_->preProcessState()<0) {
1802        // infeasible
1803        handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
1804        status_ = 0 ;
1805        secondaryStatus_ = 1;
1806        originalContinuousObjective_ = COIN_DBL_MAX;
1807        return ; 
1808      }
1809    } else {
1810      //no preprocessing
1811      originalSolver=NULL;
1812    }
1813    strategy_->setupCutGenerators(*this);
1814    strategy_->setupHeuristics(*this);
1815    // Set strategy print level to models
1816    strategy_->setupPrinting(*this,handler_->logLevel());
1817  }
1818  eventHappened_=false;
1819#ifdef COIN_USE_CLP
1820  ClpEventHandler * eventHandler=NULL;
1821 {
1822   OsiClpSolverInterface * clpSolver
1823     = dynamic_cast<OsiClpSolverInterface *> (solver_);
1824   if (clpSolver) {
1825     ClpSimplex * clpSimplex = clpSolver->getModelPtr();
1826     eventHandler = clpSimplex->eventHandler();
1827   }
1828 }
1829#endif
1830  if (!nodeCompare_)
1831    nodeCompare_=new CbcCompareDefault();;
1832  if (!problemFeasibility_)
1833    problemFeasibility_=new CbcFeasibilityBase();
1834# ifdef CBC_DEBUG
1835  std::string problemName ;
1836  solver_->getStrParam(OsiProbName,problemName) ;
1837  printf("Problem name - %s\n",problemName.c_str()) ;
1838  solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
1839# endif
1840/*
1841  Assume we're done, and see if we're proven wrong.
1842*/
1843  status_ = 0 ;
1844  secondaryStatus_ = 0;
1845  phase_=0;
1846/*
1847  Scan the variables, noting the integer variables. Create an
1848  CbcSimpleInteger object for each integer variable.
1849*/
1850  findIntegers(false) ;
1851  // If dynamic pseudo costs then do
1852  if (numberBeforeTrust_)
1853    convertToDynamic();
1854  // Set up char array to say if integer
1855  delete [] integerInfo_;
1856  {
1857    int n = solver_->getNumCols();
1858    integerInfo_ = new char [n];
1859    for (int i=0;i<n;i++) {
1860      if (solver_->isInteger(i))
1861        integerInfo_[i]=1;
1862      else
1863        integerInfo_[i]=0;
1864    }
1865  }
1866   
1867/*
1868  Ensure that objects on the lists of CbcObjects, heuristics, and cut
1869  generators attached to this model all refer to this model.
1870*/
1871  synchronizeModel() ;
1872  assert (!solverCharacteristics_);
1873  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1874  if (solverCharacteristics) {
1875    solverCharacteristics_ = solverCharacteristics;
1876  } else {
1877    // replace in solver
1878    OsiBabSolver defaultC;
1879    solver_->setAuxiliaryInfo(&defaultC);
1880    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
1881  }
1882  solverCharacteristics_->setSolver(solver_);
1883  // Set so we can tell we are in initial phase in resolve
1884  continuousObjective_ = -COIN_DBL_MAX ;
1885/*
1886  Solve the relaxation.
1887
1888  Apparently there are circumstances where this will be non-trivial --- i.e.,
1889  we've done something since initialSolve that's trashed the solution to the
1890  continuous relaxation.
1891*/
1892  bool feasible = resolve() ;
1893  if (problemFeasibility_->feasible(this,0)<0) {
1894    feasible=false; // pretend infeasible
1895  }
1896/*
1897  If the linear relaxation of the root is infeasible, bail out now. Otherwise,
1898  continue with processing the root node.
1899*/
1900  if (!feasible)
1901  { handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
1902    status_ = 0 ;
1903    secondaryStatus_ = 1;
1904    originalContinuousObjective_ = COIN_DBL_MAX;
1905    return ; }
1906  // Save objective (just so user can access it)
1907  originalContinuousObjective_ = solver_->getObjValue();
1908  bestPossibleObjective_=originalContinuousObjective_;
1909  sumChangeObjective1_=0.0;
1910  sumChangeObjective2_=0.0;
1911/*
1912  OsiRowCutDebugger knows an optimal answer for a subset of MIP problems.
1913  Assuming it recognises the problem, when called upon it will check a cut to
1914  see if it cuts off the optimal answer.
1915*/
1916  // If debugger exists set specialOptions_ bit
1917  if (solver_->getRowCutDebuggerAlways())
1918    specialOptions_ |= 1;
1919
1920# ifdef CBC_DEBUG
1921  if ((specialOptions_&1)==0)
1922    solver_->activateRowCutDebugger(problemName.c_str()) ;
1923  if (solver_->getRowCutDebuggerAlways())
1924    specialOptions_ |= 1;
1925# endif
1926
1927/*
1928  Begin setup to process a feasible root node.
1929*/
1930  bestObjective_ = CoinMin(bestObjective_,1.0e50) ;
1931  numberSolutions_ = 0 ;
1932  stateOfSearch_=0;
1933  numberHeuristicSolutions_ = 0 ;
1934  // Everything is minimization
1935  { 
1936    // needed to sync cutoffs
1937    double value ;
1938    solver_->getDblParam(OsiDualObjectiveLimit,value) ;
1939    dblParam_[CbcCurrentCutoff]= value * solver_->getObjSense();
1940  }
1941  double cutoff=getCutoff() ;
1942  double direction = solver_->getObjSense() ;
1943  dblParam_[CbcOptimizationDirection]=direction;
1944  if (cutoff < 1.0e20&&direction<0.0)
1945    messageHandler()->message(CBC_CUTOFF_WARNING1,
1946                                    messages())
1947                                      << cutoff << -cutoff << CoinMessageEol ;
1948  if (cutoff > bestObjective_)
1949    cutoff = bestObjective_ ;
1950  setCutoff(cutoff) ;
1951/*
1952  We probably already have a current solution, but just in case ...
1953*/
1954  int numberColumns = getNumCols() ;
1955  if (!currentSolution_)
1956    currentSolution_ = new double[numberColumns] ;
1957  testSolution_ = currentSolution_;
1958  /* Tell solver we are in Branch and Cut
1959     Could use last parameter for subtle differences */
1960  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
1961/*
1962  Create a copy of the solver, thus capturing the original (root node)
1963  constraint system (aka the continuous system).
1964*/
1965  continuousSolver_ = solver_->clone() ;
1966#ifdef COIN_USE_CLP
1967  OsiClpSolverInterface * clpSolver
1968    = dynamic_cast<OsiClpSolverInterface *> (solver_);
1969  if (clpSolver) {
1970    ClpSimplex * clpSimplex = clpSolver->getModelPtr();
1971    // take off names
1972    clpSimplex->dropNames();
1973  }
1974#endif
1975
1976  numberRowsAtContinuous_ = getNumRows() ;
1977/*
1978  Check the objective to see if we can deduce a nontrivial increment. If
1979  it's better than the current value for CbcCutoffIncrement, it'll be
1980  installed.
1981*/
1982  if(solverCharacteristics_->reducedCostsAccurate())
1983    analyzeObjective() ;
1984/*
1985  Set up for cut generation. addedCuts_ holds the cuts which are relevant for
1986  the active subproblem. whichGenerator will be used to record the generator
1987  that produced a given cut.
1988*/
1989  maximumWhich_ = 1000 ;
1990  delete [] whichGenerator_;
1991  whichGenerator_ = new int[maximumWhich_] ;
1992  memset(whichGenerator_,0,maximumWhich_*sizeof(int));
1993  maximumNumberCuts_ = 0 ;
1994  currentNumberCuts_ = 0 ;
1995  delete [] addedCuts_ ;
1996  addedCuts_ = NULL ;
1997/*
1998  Set up an empty heap and associated data structures to hold the live set
1999  (problems which require further exploration).
2000*/
2001  tree_->setComparison(*nodeCompare_) ;
2002/*
2003  Used to record the path from a node to the root of the search tree, so that
2004  we can then traverse from the root to the node when restoring a subproblem.
2005*/
2006  maximumDepth_ = 10 ;
2007  delete [] walkback_ ;
2008  walkback_ = new CbcNodeInfo * [maximumDepth_] ;
2009/*
2010  Used to generate bound edits for CbcPartialNodeInfo.
2011*/
2012  double * lowerBefore = new double [numberColumns] ;
2013  double * upperBefore = new double [numberColumns] ;
2014/*
2015 
2016  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
2017  lifting. It will iterate a generate/reoptimise loop (including reduced cost
2018  fixing) until no cuts are generated, the change in objective falls off,  or
2019  the limit on the number of rounds of cut generation is exceeded.
2020
2021  At the end of all this, any cuts will be recorded in cuts and also
2022  installed in the solver's constraint system. We'll have reoptimised, and
2023  removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
2024  adjusted accordingly).
2025
2026  Tell cut generators they can be a bit more aggressive at root node
2027
2028  TODO: Why don't we make a copy of the solution after solveWithCuts?
2029  TODO: If numberUnsatisfied == 0, don't we have a solution?
2030*/
2031  phase_=1;
2032  int iCutGenerator;
2033  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
2034    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2035    generator->setAggressiveness(generator->getAggressiveness()+100);
2036  }
2037  OsiCuts cuts ;
2038  int anyAction = -1 ;
2039  numberOldActiveCuts_ = 0 ;
2040  numberNewCuts_ = 0 ;
2041  // Array to mark solution
2042  delete [] usedInSolution_;
2043  usedInSolution_ = new int[numberColumns];
2044  CoinZeroN(usedInSolution_,numberColumns);
2045/*
2046  For printing totals and for CbcNode (numberNodes_)
2047*/
2048  numberIterations_ = 0 ;
2049  numberNodes_ = 0 ;
2050  numberNodes2_ = 0 ;
2051  maximumStatistics_=0;
2052  statistics_ = NULL;
2053  // Do on switch
2054  if (doStatistics) {
2055    maximumStatistics_=10000;
2056    statistics_ = new CbcStatistics * [maximumStatistics_];
2057    memset(statistics_,0,maximumStatistics_*sizeof(CbcStatistics *));
2058  }
2059
2060  { int iObject ;
2061    int preferredWay ;
2062    int numberUnsatisfied = 0 ;
2063    memcpy(currentSolution_,solver_->getColSolution(),
2064           numberColumns*sizeof(double)) ;
2065
2066    for (iObject = 0 ; iObject < numberObjects_ ; iObject++)
2067    { double infeasibility =
2068          object_[iObject]->infeasibility(preferredWay) ;
2069      if (infeasibility ) numberUnsatisfied++ ; }
2070    if (numberUnsatisfied) { 
2071      feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,NULL);
2072    }
2073  }
2074  // make cut generators less aggressive
2075  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
2076    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
2077    generator->setAggressiveness(generator->getAggressiveness()-100);
2078  }
2079  currentNumberCuts_ = numberNewCuts_ ;
2080/*
2081  We've taken the continuous relaxation as far as we can. Time to branch.
2082  The first order of business is to actually create a node. chooseBranch
2083  currently uses strong branching to evaluate branch object candidates,
2084  unless forced back to simple branching. If chooseBranch concludes that a
2085  branching candidate is monotone (anyAction == -1) or infeasible (anyAction
2086  == -2) when forced to integer values, it returns here immediately.
2087
2088  Monotone variables trigger a call to resolve(). If the problem remains
2089  feasible, try again to choose a branching variable. At the end of the loop,
2090  resolved == true indicates that some variables were fixed.
2091
2092  Loss of feasibility will result in the deletion of newNode.
2093*/
2094
2095  bool resolved = false ;
2096  CbcNode *newNode = NULL ;
2097  numberFixedAtRoot_=0;
2098  numberFixedNow_=0;
2099  int numberIterationsAtContinuous = numberIterations_;
2100  if (feasible) {
2101    newNode = new CbcNode ;
2102    // Set objective value (not so obvious if NLP etc)
2103    setObjectiveValue(newNode,NULL);
2104    anyAction = -1 ;
2105    // To make depth available we may need a fake node
2106    CbcNode fakeNode;
2107    if (!currentNode_) {
2108      // Not true if sub trees assert (!numberNodes_);
2109      currentNode_=&fakeNode;
2110    }
2111    phase_=3;
2112    // only allow 1000 passes
2113    int numberPassesLeft=1000;
2114    // This is first crude step
2115    if (numberAnalyzeIterations_) {
2116      delete [] analyzeResults_;
2117      analyzeResults_ = new double [4*numberIntegers_];
2118      numberFixedAtRoot_=newNode->analyze(this,analyzeResults_);
2119      if (numberFixedAtRoot_>0) {
2120        printf("%d fixed by analysis\n",numberFixedAtRoot_);
2121        setPointers(solver_);
2122        numberFixedNow_ = numberFixedAtRoot_;
2123      } else if (numberFixedAtRoot_<0) {
2124        printf("analysis found to be infeasible\n");
2125        anyAction=-2;
2126        delete newNode ;
2127        newNode = NULL ;
2128        feasible = false ;
2129      }
2130    }
2131    while (anyAction == -1)
2132    {
2133      // Set objective value (not so obvious if NLP etc)
2134      setObjectiveValue(newNode,NULL);
2135      if (numberBeforeTrust_==0 ) {
2136        anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ;
2137      } else {
2138        OsiSolverBranch * branches=NULL;
2139        anyAction = newNode->chooseDynamicBranch(this,NULL,branches,numberPassesLeft) ;
2140        if (anyAction==-3) 
2141          anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ; // dynamic did nothing
2142      }
2143      numberPassesLeft--;
2144      if (anyAction == -1)
2145      { feasible = resolve() ;
2146      if (problemFeasibility_->feasible(this,0)<0) {
2147        feasible=false; // pretend infeasible
2148      }
2149      reducedCostFix() ;
2150      resolved = true ;
2151#       ifdef CBC_DEBUG
2152      printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
2153             solver_->getObjValue(),
2154             solver_->getNumRows()) ;
2155#       endif
2156      if (!feasible) anyAction = -2 ; }
2157      if (anyAction == -2||newNode->objectiveValue() >= cutoff)
2158      { delete newNode ;
2159        newNode = NULL ;
2160        feasible = false ; } } }
2161/*
2162  At this point, the root subproblem is infeasible or fathomed by bound
2163  (newNode == NULL), or we're live with an objective value that satisfies the
2164  current objective cutoff.
2165*/
2166  assert (!newNode || newNode->objectiveValue() <= cutoff) ;
2167  // Save address of root node as we don't want to delete it
2168  //CbcNode * rootNode = newNode;
2169/*
2170  The common case is that the lp relaxation is feasible but doesn't satisfy
2171  integrality (i.e., newNode->variable() >= 0, indicating we've been able to
2172  select a branching variable). Remove any cuts that have gone slack due to
2173  forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
2174  basis (via createInfo()) and stash the new cuts in the nodeInfo (via
2175  addCuts()). If, by some miracle, we have an integral solution at the root
2176  (newNode->variable() < 0), takeOffCuts() will ensure that the solver holds
2177  a valid solution for use by setBestSolution().
2178*/
2179  if (feasible && newNode->variable() >= 0)
2180  { if (resolved)
2181    { bool needValidSolution = (newNode->variable() < 0) ;
2182      takeOffCuts(cuts,needValidSolution,NULL) ;
2183#     ifdef CHECK_CUT_COUNTS
2184      { printf("Number of rows after chooseBranch fix (root)"
2185               "(active only) %d\n",
2186                numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ;
2187        const CoinWarmStartBasis* debugws =
2188          dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
2189        debugws->print() ;
2190        delete debugws ; }
2191#     endif
2192    }
2193    newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
2194    newNode->nodeInfo()->addCuts(cuts,
2195                                 newNode->numberBranches(),whichGenerator_) ;
2196  }
2197/*
2198  Continuous data to be used later
2199*/
2200  continuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
2201  continuousInfeasibilities_ = 0 ;
2202  if (newNode)
2203  { continuousObjective_ = newNode->objectiveValue() ;
2204    delete [] continuousSolution_;
2205    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
2206                                             numberColumns);
2207    continuousInfeasibilities_ = newNode->numberUnsatisfied() ; }
2208/*
2209  Bound may have changed so reset in objects
2210*/
2211  { int i ;
2212    for (i = 0;i < numberObjects_;i++)
2213      object_[i]->resetBounds() ; }
2214  stoppedOnGap_ = false ;
2215/*
2216  Feasible? Then we should have either a live node prepped for future
2217  expansion (indicated by variable() >= 0), or (miracle of miracles) an
2218  integral solution at the root node.
2219
2220  initializeInfo sets the reference counts in the nodeInfo object.  Since
2221  this node is still live, push it onto the heap that holds the live set.
2222*/
2223  double bestValue = 0.0 ;
2224  if (newNode) {
2225    bestValue = newNode->objectiveValue();
2226    if (newNode->variable() >= 0) {
2227      newNode->initializeInfo() ;
2228      tree_->push(newNode) ;
2229      if (statistics_) {
2230        if (numberNodes2_==maximumStatistics_) {
2231          maximumStatistics_ = 2*maximumStatistics_;
2232          CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
2233          memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
2234          memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
2235          delete [] statistics_;
2236          statistics_=temp;
2237        }
2238        assert (!statistics_[numberNodes2_]);
2239        statistics_[numberNodes2_]=new CbcStatistics(newNode);
2240      }
2241      numberNodes2_++;
2242#     ifdef CHECK_NODE
2243      printf("Node %x on tree\n",newNode) ;
2244#     endif
2245    } else {
2246      // continuous is integer
2247      double objectiveValue = newNode->objectiveValue();
2248      setBestSolution(CBC_SOLUTION,objectiveValue,
2249                      solver_->getColSolution()) ;
2250      delete newNode ;
2251      newNode = NULL ;
2252    }
2253  }
2254
2255  if (printFrequency_ <= 0) {
2256    printFrequency_ = 1000 ;
2257    if (getNumCols() > 2000)
2258      printFrequency_ = 100 ;
2259  }
2260  /*
2261    It is possible that strong branching fixes one variable and then the code goes round
2262    again and again.  This can take too long.  So we need to warn user - just once.
2263  */
2264  numberLongStrong_=0;
2265/*
2266  At last, the actual branch-and-cut search loop, which will iterate until
2267  the live set is empty or we hit some limit (integrality gap, time, node
2268  count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
2269  solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
2270  add the node to the live set.
2271
2272  The first action is to winnow the live set to remove nodes which are worse
2273  than the current objective cutoff.
2274*/
2275  while (!tree_->empty()) {
2276    if (cutoff > getCutoff()) {
2277      double newCutoff = getCutoff();
2278      if (analyzeResults_) {
2279        // see if we could fix any (more)
2280        int n=0;
2281        double * newLower = analyzeResults_;
2282        double * objLower = newLower+numberIntegers_;
2283        double * newUpper = objLower+numberIntegers_;
2284        double * objUpper = newUpper+numberIntegers_;
2285        for (int i=0;i<numberIntegers_;i++) {
2286          if (objLower[i]>newCutoff) {
2287            n++;
2288            if (objUpper[i]>newCutoff) {
2289              newCutoff = -COIN_DBL_MAX;
2290              break;
2291            }
2292          } else if (objUpper[i]>newCutoff) {
2293            n++;
2294          }
2295        }
2296        if (newCutoff==-COIN_DBL_MAX) {
2297          printf("Root analysis says finished\n");
2298        } else if (n>numberFixedNow_) {
2299          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
2300          numberFixedNow_=n;
2301        }
2302      }
2303#ifdef COIN_USE_CLP
2304      if (eventHandler) {
2305        if (!eventHandler->event(ClpEventHandler::solution)) {
2306          eventHappened_=true; // exit
2307        }
2308      }
2309#endif
2310      // Do from deepest
2311      tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ;
2312      nodeCompare_->newSolution(this) ;
2313      nodeCompare_->newSolution(this,continuousObjective_,
2314                                continuousInfeasibilities_) ;
2315      tree_->setComparison(*nodeCompare_) ;
2316      if (tree_->empty())
2317        break; // finished
2318    }
2319    cutoff = getCutoff() ;
2320    /*
2321      Periodic activities: Opportunities to
2322      + tweak the nodeCompare criteria,
2323      + check if we've closed the integrality gap enough to quit,
2324      + print a summary line to let the user know we're working
2325    */
2326    if ((numberNodes_%1000) == 0) {
2327      bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ;
2328      // redo tree if wanted
2329      if (redoTree)
2330        tree_->setComparison(*nodeCompare_) ;
2331    }
2332    if ((numberNodes_%printFrequency_) == 0) {
2333      int j ;
2334      int nNodes = tree_->size() ;
2335      bestPossibleObjective_ = 1.0e100 ;
2336      for (j = 0;j < nNodes;j++) {
2337        CbcNode * node = tree_->nodePointer(j) ;
2338        if (node&&node->objectiveValue() < bestPossibleObjective_)
2339          bestPossibleObjective_ = node->objectiveValue() ;
2340      }
2341      messageHandler()->message(CBC_STATUS,messages())
2342        << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
2343        << CoinMessageEol ;
2344#ifdef COIN_USE_CLP
2345      if (eventHandler) {
2346        if (!eventHandler->event(ClpEventHandler::treeStatus)) {
2347          eventHappened_=true; // exit
2348        }
2349      }
2350#endif
2351    }
2352    // If no solution but many nodes - signal change in strategy
2353    if (numberNodes_>2*numberObjects_+1000&&stateOfSearch_!=2)
2354      stateOfSearch_=3;
2355    // See if can stop on gap
2356    double testGap = CoinMax(dblParam_[CbcAllowableGap],
2357                             CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_))
2358                             *dblParam_[CbcAllowableFractionGap]);
2359    if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
2360      stoppedOnGap_ = true ;
2361    }
2362
2363#   ifdef CHECK_NODE_FULL
2364    verifyTreeNodes(tree_,*this) ;
2365#   endif
2366#   ifdef CHECK_CUT_COUNTS
2367    verifyCutCounts(tree_,*this) ;
2368#   endif
2369   
2370    /*
2371      Now we come to the meat of the loop. To create the active subproblem, we'll
2372      pop the most promising node in the live set, rebuild the subproblem it
2373      represents, and then execute the current arm of the branch to create the
2374      active subproblem.
2375    */
2376    CbcNode *node = tree_->bestNode(cutoff) ;
2377    // Possible one on tree worse than cutoff
2378    if (!node)
2379      continue;
2380    // Solve node
2381    int status;
2382    int numberNewNodes;
2383    CbcNode ** newNodes = solveOneNode(-1,node,numberNewNodes,status);
2384    // put nodes back on tree
2385    for (int iPush = 0; iPush<numberNewNodes;iPush++) {
2386      tree_->push(newNodes[iPush]) ;
2387    }     
2388    delete [] newNodes;
2389    if (status==2) {
2390      /*
2391        The next block of code is executed if we've
2392        aborted because we hit one of the limits. Clean up by deleting the live set
2393        and break out of the node processing loop.
2394      */
2395      tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
2396      delete nextRowCut_;
2397      // We need to get rid of node if is has already been popped from tree
2398      //if (!nodeOnTree&&!stoppedOnGap_&&node!=rootNode)
2399      //delete node;
2400      double totalTime = CoinCpuTime()-dblParam_[CbcStartSeconds] ;
2401      if (stoppedOnGap_) {
2402        messageHandler()->message(CBC_GAP,messages())
2403          << bestObjective_-bestPossibleObjective_
2404          << dblParam_[CbcAllowableGap]
2405          << dblParam_[CbcAllowableFractionGap]*100.0
2406          << CoinMessageEol ;
2407        secondaryStatus_ = 2;
2408        status_ = 0 ;
2409      } else if (isNodeLimitReached()) {
2410        handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
2411        secondaryStatus_ = 3;
2412        status_ = 1 ;
2413      } else if (totalTime >= dblParam_[CbcMaximumSeconds]) {
2414        handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
2415        secondaryStatus_ = 4;
2416        status_ = 1 ;
2417      } else if (eventHappened_) {
2418        handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
2419        secondaryStatus_ = 5;
2420        status_ = 5 ;
2421      } else { 
2422        handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
2423        secondaryStatus_ = 6;
2424        status_ = 1 ;
2425      }
2426      break ;
2427    } else {
2428      /*
2429        Delete cuts to get back to the original system.
2430       
2431        I'm thinking this is redundant --- the call to addCuts that conditions entry
2432        to this code block also performs this action.
2433      */
2434      int numberToDelete = getNumRows()-numberRowsAtContinuous_ ;
2435      if (numberToDelete) {
2436        int * delRows = new int[numberToDelete] ;
2437        int i ;
2438        for (i = 0 ; i < numberToDelete ; i++) 
2439          delRows[i] = i+numberRowsAtContinuous_ ; 
2440        solver_->deleteRows(numberToDelete,delRows) ;
2441        delete [] delRows ;
2442      }
2443    }
2444  }
2445  /*
2446    That's it, we've exhausted the search tree, or broken out of the loop because
2447    we hit some limit on evaluation.
2448   
2449    We may have got an intelligent tree so give it one more chance
2450  */
2451  // Tell solver we are not in Branch and Cut
2452  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
2453  tree_->endSearch();
2454  //  If we did any sub trees - did we give up on any?
2455  if ( numberStoppedSubTrees_)
2456    status_=1;
2457  if (!status_) {
2458    bestPossibleObjective_=bestObjective_;
2459    handler_->message(CBC_END_GOOD,messages_)
2460      << bestObjective_ << numberIterations_ << numberNodes_
2461      << CoinMessageEol ;
2462  } else {
2463    handler_->message(CBC_END,messages_)
2464      << bestObjective_ <<bestPossibleObjective_
2465      << numberIterations_ << numberNodes_
2466      << CoinMessageEol ;
2467  }
2468  if (numberStrongIterations_)
2469    handler_->message(CBC_STRONG_STATS,messages_)
2470      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
2471      << strongInfo_[1] << CoinMessageEol ;
2472  if (statistics_) {
2473    // report in some way
2474    int * lookup = new int[numberObjects_];
2475    int i;
2476    for (i=0;i<numberObjects_;i++) 
2477      lookup[i]=-1;
2478    bool goodIds=true;
2479    for (i=0;i<numberObjects_;i++) {
2480      int id = object_[i]->id();
2481      int iColumn = object_[i]->columnNumber();
2482      if (iColumn<0)
2483        iColumn = id+numberColumns;
2484      if(id>=0&&id<numberObjects_) {
2485        if (lookup[id]==-1) {
2486          lookup[id]=iColumn;
2487        } else {
2488          goodIds=false;
2489          break;
2490        }
2491      } else {
2492        goodIds=false;
2493        break;
2494      }
2495    }
2496    if (!goodIds) {
2497      delete [] lookup;
2498      lookup=NULL;
2499    }
2500    if (doStatistics==3) {
2501      printf("  node parent depth column   value                    obj      inf\n");
2502      for ( i=0;i<numberNodes2_;i++) {
2503        statistics_[i]->print(lookup);
2504      }
2505    }
2506    if (doStatistics>1) {
2507      // Find last solution
2508      int k;
2509      for (k=numberNodes2_-1;k>=0;k--) {
2510        if (statistics_[k]->endingObjective()!=COIN_DBL_MAX&&
2511            !statistics_[k]->endingInfeasibility())
2512          break;
2513      }
2514      if (k>=0) {
2515        int depth=statistics_[k]->depth();
2516        int * which = new int[depth+1];
2517        for (i=depth;i>=0;i--) {
2518          which[i]=k;
2519          k=statistics_[k]->parentNode();
2520        }
2521        printf("  node parent depth column   value                    obj      inf\n");
2522        for (i=0;i<=depth;i++) {
2523          statistics_[which[i]]->print(lookup);
2524        }
2525        delete [] which;
2526      }
2527    }
2528    // now summary
2529    int maxDepth=0;
2530    double averageSolutionDepth=0.0;
2531    int numberSolutions=0;
2532    double averageCutoffDepth=0.0;
2533    double averageSolvedDepth=0.0;
2534    int numberCutoff=0;
2535    int numberDown=0;
2536    int numberFirstDown=0;
2537    double averageInfDown=0.0;
2538    double averageObjDown=0.0;
2539    int numberCutoffDown=0;
2540    int numberUp=0;
2541    int numberFirstUp=0;
2542    double averageInfUp=0.0;
2543    double averageObjUp=0.0;
2544    int numberCutoffUp=0;
2545    double averageNumberIterations1=0.0;
2546    double averageValue=0.0;
2547    for ( i=0;i<numberNodes2_;i++) {
2548      int depth =  statistics_[i]->depth(); 
2549      int way =  statistics_[i]->way(); 
2550      double value = statistics_[i]->value(); 
2551      double startingObjective =  statistics_[i]->startingObjective(); 
2552      int startingInfeasibility = statistics_[i]->startingInfeasibility(); 
2553      double endingObjective = statistics_[i]->endingObjective(); 
2554      int endingInfeasibility = statistics_[i]->endingInfeasibility(); 
2555      maxDepth = CoinMax(depth,maxDepth);
2556      // Only for completed
2557      averageNumberIterations1 += statistics_[i]->numberIterations();
2558      averageValue += value;
2559      if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) {
2560        numberSolutions++;
2561        averageSolutionDepth += depth;
2562      }
2563      if (endingObjective==COIN_DBL_MAX) {
2564        numberCutoff++;
2565        averageCutoffDepth += depth;
2566        if (way<0) {
2567          numberDown++;
2568          numberCutoffDown++;
2569          if (way==-1)
2570            numberFirstDown++;
2571        } else {
2572          numberUp++;
2573          numberCutoffUp++;
2574          if (way==1)
2575            numberFirstUp++;
2576        }
2577      } else {
2578        averageSolvedDepth += depth;
2579        if (way<0) {
2580          numberDown++;
2581          averageInfDown += startingInfeasibility-endingInfeasibility;
2582          averageObjDown += endingObjective-startingObjective;
2583          if (way==-1)
2584            numberFirstDown++;
2585        } else {
2586          numberUp++;
2587          averageInfUp += startingInfeasibility-endingInfeasibility;
2588          averageObjUp += endingObjective-startingObjective;
2589          if (way==1)
2590            numberFirstUp++;
2591        }
2592      }
2593    }
2594    // Now print
2595    if (numberSolutions)
2596      averageSolutionDepth /= (double) numberSolutions;
2597    int numberSolved = numberNodes2_-numberCutoff;
2598    double averageNumberIterations2=numberIterations_-averageNumberIterations1
2599      -numberIterationsAtContinuous;
2600    if(numberCutoff) {
2601      averageCutoffDepth /= (double) numberCutoff;
2602      averageNumberIterations2 /= (double) numberCutoff;
2603    }
2604    if (numberNodes2_) 
2605      averageValue /= (double) numberNodes2_;
2606    if (numberSolved) {
2607      averageNumberIterations1 /= (double) numberSolved;
2608      averageSolvedDepth /= (double) numberSolved;
2609    }
2610    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
2611           numberSolutions,averageSolutionDepth);
2612    printf("average value of variable being branched on was %g\n",
2613           averageValue);
2614    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
2615           numberCutoff,averageCutoffDepth,averageNumberIterations2);
2616    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
2617           numberSolved,averageSolvedDepth,averageNumberIterations1);
2618    if (numberDown) {
2619      averageInfDown /= (double) numberDown;
2620      averageObjDown /= (double) numberDown;
2621    }
2622    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
2623           numberDown,numberFirstDown,numberDown-numberFirstDown,numberCutoffDown,
2624           averageInfDown,averageObjDown);
2625    if (numberUp) {
2626      averageInfUp /= (double) numberUp;
2627      averageObjUp /= (double) numberUp;
2628    }
2629    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
2630           numberUp,numberFirstUp,numberUp-numberFirstUp,numberCutoffUp,
2631           averageInfUp,averageObjUp);
2632    for ( i=0;i<numberNodes2_;i++) 
2633      delete statistics_[i];
2634    delete [] statistics_;
2635    statistics_=NULL;
2636    maximumStatistics_=0;
2637    delete [] lookup;
2638  }
2639/*
2640  If we think we have a solution, restore and confirm it with a call to
2641  setBestSolution().  We need to reset the cutoff value so as not to fathom
2642  the solution on bounds.  Note that calling setBestSolution( ..., true)
2643  leaves the continuousSolver_ bounds vectors fixed at the solution value.
2644
2645  Running resolve() here is a failsafe --- setBestSolution has already
2646  reoptimised using the continuousSolver_. If for some reason we fail to
2647  prove optimality, run the problem again after instructing the solver to
2648  tell us more.
2649
2650  If all looks good, replace solver_ with continuousSolver_, so that the
2651  outside world will be able to obtain information about the solution using
2652  public methods.
2653*/
2654  if (bestSolution_) {
2655    setCutoff(1.0e50) ; // As best solution should be worse than cutoff
2656    phase_=5;
2657    setBestSolution(CBC_SOLUTION,bestObjective_,bestSolution_,true) ;
2658    continuousSolver_->resolve() ;
2659    if (!continuousSolver_->isProvenOptimal()) {
2660      continuousSolver_->messageHandler()->setLogLevel(2) ;
2661      continuousSolver_->initialSolve() ;
2662    }
2663    delete solver_ ;
2664    solver_ = continuousSolver_ ;
2665    setPointers(solver_);
2666    continuousSolver_ = NULL ;
2667  }
2668  /*
2669    Clean up dangling objects. continuousSolver_ may already be toast.
2670  */
2671  delete [] whichGenerator_ ;
2672  whichGenerator_=NULL;
2673  delete [] lowerBefore ;
2674  delete [] upperBefore ;
2675  delete [] walkback_ ;
2676  walkback_ = NULL ;
2677  delete [] addedCuts_ ;
2678  addedCuts_ = NULL ;
2679  // Get rid of characteristics
2680  solverCharacteristics_=NULL;
2681  if (continuousSolver_) {
2682    delete continuousSolver_ ;
2683    continuousSolver_ = NULL ;
2684  }
2685  /*
2686    Destroy global cuts by replacing with an empty OsiCuts object.
2687  */
2688  globalCuts_= OsiCuts() ;
2689  if (strategy_&&strategy_->preProcessState()>0) {
2690    // undo preprocessing
2691    CglPreProcess * process = strategy_->process();
2692    assert (process);
2693    int n = originalSolver->getNumCols();
2694    if (bestSolution_) {
2695      delete [] bestSolution_;
2696      bestSolution_ = new double [n];
2697      process->postProcess(*solver_);
2698    }
2699    strategy_->deletePreProcess();
2700    // Solution now back in originalSolver
2701    delete solver_;
2702    solver_=originalSolver;
2703    if (bestSolution_)
2704      memcpy(bestSolution_,solver_->getColSolution(),n*sizeof(double));
2705  }
2706  return ;
2707}
2708/* Input one node output N nodes to put on tree and optional solution update
2709   This should be able to operate in parallel so is given a solver and is const(ish)
2710   However we will need to keep an array of solver_ and bases and more
2711   status is 0 for normal, 1 if solution, 2 if should stop
2712   Calling code should always push nodes back on tree
2713 */
2714CbcNode ** 
2715CbcModel::solveOneNode(int whichSolver,CbcNode * node, 
2716                    int & numberNodesOutput, int & status) 
2717{
2718  numberNodesOutput=0;
2719  status=0;
2720  // Maximum number of new nodes will be known
2721  // For normal case can be 2 (old + new)
2722  int maximumNewNodes=2+(1<<sizeMiniTree_);
2723  CbcNode ** newNodes = new CbcNode * [maximumNewNodes];
2724 
2725  currentNode_=node; // so can be accessed elsewhere
2726#ifdef CBC_DEBUG
2727  printf("%d unsat, way %d, obj %g est %g\n",
2728         node->numberUnsatisfied(),node->way(),node->objectiveValue(),
2729         node->guessedObjectiveValue());
2730#endif
2731  // Save clone in branching decision
2732  if(branchingMethod_)
2733    branchingMethod_->saveBranchingObject(node->modifiableBranchingObject());
2734  // Say not on optimal path
2735  bool onOptimalPath=false;
2736#   ifdef CHECK_NODE
2737  /*
2738    WARNING: The use of integerVariable_[*] here will break as soon as the
2739    branching object is something other than an integer variable.
2740    This needs some thought.
2741  */
2742  printf("Node %x popped from tree - %d left, %d count\n",node,
2743         node->nodeInfo()->numberBranchesLeft(),
2744         node->nodeInfo()->numberPointingToThis()) ;
2745  printf("\tdepth = %d, z =  %g, unsat = %d, var = %d.\n",
2746         node->depth(),node->objectiveValue(),
2747         node->numberUnsatisfied(),
2748         integerVariable_[node->variable()]) ;
2749#   endif
2750 
2751  /*
2752    Rebuild the subproblem for this node:        Call addCuts() to adjust the model
2753    to recreate the subproblem for this node (set proper variable bounds, add
2754    cuts, create a basis).  This may result in the problem being fathomed by
2755    bound or infeasibility. Returns 1 if node is fathomed.
2756    Execute the current arm of the branch: If the problem survives, save the
2757    resulting variable bounds and call branch() to modify variable bounds
2758    according to the current arm of the branching object. If we're processing
2759    the final arm of the branching object, flag the node for removal from the
2760    live set.
2761  */
2762  CbcNodeInfo * nodeInfo = node->nodeInfo() ;
2763  // empty basis
2764  CoinWarmStartBasis * lastws = new CoinWarmStartBasis();
2765  int numberColumns = getNumCols() ;
2766  double * lowerBefore = new double[numberColumns];
2767  double * upperBefore = new double[numberColumns];
2768  if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_)) {
2769    int i ;
2770    const double * lower = getColLower() ;
2771    const double * upper = getColUpper() ;
2772    for (i = 0 ; i < numberColumns ; i++) {
2773      lowerBefore[i]= lower[i] ;
2774      upperBefore[i]= upper[i] ;
2775    }
2776    bool deleteNode ;
2777    if (node->branch()) { 
2778      // set nodenumber correctly
2779      node->nodeInfo()->setNodeNumber(numberNodes2_);
2780      if (statistics_) {
2781        if (numberNodes2_==maximumStatistics_) {
2782          maximumStatistics_ = 2*maximumStatistics_;
2783          CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
2784          memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
2785          memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
2786          delete [] statistics_;
2787          statistics_=temp;
2788        }
2789        assert (!statistics_[numberNodes2_]);
2790        statistics_[numberNodes2_]=new CbcStatistics(node);
2791      }
2792      numberNodes2_++;
2793      newNodes[numberNodesOutput++]=node;
2794      deleteNode = false ;
2795#       ifdef CHECK_NODE
2796      printf("Node %x pushed back on tree - %d left, %d count\n",node,
2797             nodeInfo->numberBranchesLeft(),
2798             nodeInfo->numberPointingToThis()) ;
2799#       endif
2800    } else { 
2801      deleteNode = true ;
2802    }
2803   
2804    if ((specialOptions_&1)!=0) {
2805      /*
2806        This doesn't work as intended --- getRowCutDebugger will return null
2807        unless the current feasible solution region includes the optimal solution
2808        that RowCutDebugger knows. There's no way to tell inactive from off the
2809        optimal path.
2810      */
2811      const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
2812      if (debugger) {
2813        onOptimalPath=true;
2814        printf("On optimal path\n") ;
2815      }
2816    }
2817    /*
2818      Reoptimize, possibly generating cuts and/or using heuristics to find
2819      solutions.  Cut reference counts are unaffected unless we lose feasibility,
2820      in which case solveWithCuts() will make the adjustment.
2821    */
2822    phase_=2;
2823    OsiCuts cuts = OsiCuts() ;
2824    int saveNumber = numberIterations_;
2825    if(solverCharacteristics_->solutionAddsCuts()) {
2826      feasible=resolve();
2827      if (feasible) {
2828        int iObject ;
2829        int preferredWay ;
2830        int numberUnsatisfied = 0 ;
2831        memcpy(currentSolution_,solver_->getColSolution(),
2832               numberColumns*sizeof(double)) ;
2833       
2834        for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2835          double infeasibility =
2836            object_[iObject]->infeasibility(preferredWay) ;
2837          if (infeasibility ) numberUnsatisfied++ ;
2838        }
2839        if (numberUnsatisfied)   {
2840          feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
2841                                   NULL);
2842        } else {
2843          // may generate cuts and turn the solution
2844          //to an infeasible one
2845          feasible = solveWithCuts(cuts, 1,
2846                                   NULL);
2847          currentNumberCuts_ = cuts.sizeRowCuts();
2848            if (currentNumberCuts_ >= maximumNumberCuts_) {
2849              maximumNumberCuts_ = currentNumberCuts;
2850              delete [] addedCuts_;
2851              addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
2852            }
2853        }
2854        // check extra info on feasibility
2855        if (!solverCharacteristics_->mipFeasible())
2856          feasible = false;
2857      }
2858    } else {
2859      // normal
2860      feasible = solveWithCuts(cuts,maximumCutPasses_,node);
2861    }
2862    if (statistics_) {
2863      assert (numberNodes2_);
2864      assert (statistics_[numberNodes2_-1]);
2865      assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
2866      statistics_[numberNodes2_-1]->endOfBranch(numberIterations_-saveNumber,
2867                                               feasible ? solver_->getObjValue()
2868                                               : COIN_DBL_MAX);
2869    }
2870    /*
2871      Check for abort on limits: node count, solution count, time, integrality gap.
2872    */
2873    double totalTime = CoinCpuTime()-dblParam_[CbcStartSeconds] ;
2874    if (numberNodes_ < intParam_[CbcMaxNumNode] &&
2875        numberSolutions_ < intParam_[CbcMaxNumSol] &&
2876        totalTime < dblParam_[CbcMaximumSeconds] &&
2877        !stoppedOnGap_&&!eventHappened_) {
2878      /*
2879        Are we still feasible? If so, create a node and do the work to attach a
2880        branching object, reoptimising as needed if chooseBranch() identifies
2881        monotone objects.
2882       
2883        Finally, attach a partial nodeInfo object and store away any cuts that we
2884        created back in solveWithCuts. addCuts() will also deal with the cut
2885        reference counts.
2886      */
2887      if (onOptimalPath)
2888        assert (feasible);
2889      bool checkingNode=false;
2890      int anyAction;
2891      double direction =  dblParam_[CbcOptimizationDirection];
2892      bool resolved = false;
2893      CbcNode * newNode=NULL;
2894      if (feasible) {
2895        newNode = new CbcNode ;
2896        // Set objective value (not so obvious if NLP etc)
2897        setObjectiveValue(newNode,node);
2898        anyAction =-1 ;
2899        resolved = false ;
2900        if (newNode->objectiveValue() >= getCutoff()) 
2901          anyAction=-2;
2902        // only allow twenty passes
2903        int numberPassesLeft=20;
2904        checkingNode=true;
2905        OsiSolverBranch * branches=NULL;
2906        while (anyAction == -1) { 
2907          // Set objective value (not so obvious if NLP etc)
2908          setObjectiveValue(newNode,node);
2909          if (numberBeforeTrust_==0 ) {
2910            anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ;
2911          } else {
2912            anyAction = newNode->chooseDynamicBranch(this,node,branches,numberPassesLeft) ;
2913            if (anyAction==-3) 
2914              anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ; // dynamic did nothing
2915          }
2916          if (onOptimalPath)
2917            assert (anyAction!=-2); // can be useful but gives false positives on strong
2918          numberPassesLeft--;
2919          if (numberPassesLeft<=-1) {
2920            if (!numberLongStrong_)
2921              messageHandler()->message(CBC_WARNING_STRONG,
2922                                        messages()) << CoinMessageEol ;
2923            numberLongStrong_++;
2924          }
2925          if (anyAction == -1) {
2926            // can do quick optimality check
2927            int easy=2;
2928            solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
2929            feasible = resolve() ;
2930            solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2931            resolved = true ;
2932            if (problemFeasibility_->feasible(this,0)<0) {
2933              feasible=false; // pretend infeasible
2934            }
2935            if (feasible) {
2936              // Set objective value (not so obvious if NLP etc)
2937              setObjectiveValue(newNode,node);
2938              reducedCostFix() ;
2939              if (newNode->objectiveValue() >= getCutoff()) 
2940                anyAction=-2;
2941            } else {
2942              anyAction = -2 ;
2943            }
2944          }
2945        }
2946        if (anyAction >= 0) {
2947          if (resolved) {
2948            bool needValidSolution = (newNode->variable() < 0) ;
2949            takeOffCuts(cuts,needValidSolution,NULL) ; 
2950#             ifdef CHECK_CUT_COUNTS
2951            { printf("Number of rows after chooseBranch fix (node)"
2952                     "(active only) %d\n",
2953                     numberRowsAtContinuous_+numberNewCuts_+
2954                     numberOldActiveCuts_) ;
2955            const CoinWarmStartBasis* debugws =
2956              dynamic_cast<const CoinWarmStartBasis*>
2957              (solver_->getWarmStart()) ;
2958            debugws->print() ;
2959            delete debugws ; }
2960#             endif
2961          }
2962          if (anyAction==0||!newNode->numberUnsatisfied()) {
2963            delete [] branches;
2964            newNode->createInfo(this,node,lastws,lowerBefore,upperBefore,
2965                                numberOldActiveCuts_,numberNewCuts_) ;
2966            if (newNode->numberUnsatisfied()) {
2967              newNode->initializeInfo() ;
2968              newNode->nodeInfo()->addCuts(cuts,newNode->numberBranches(),
2969                                           whichGenerator_) ;
2970            }
2971          } else {
2972            // mini tree
2973            int numberBranches=anyAction;
2974            int numberResults=1<<numberBranches;
2975            OsiSolverResult * result = new OsiSolverResult[numberResults];
2976            double * lowerBefore2 = new double[numberColumns];
2977            double * upperBefore2 = new double[numberColumns];
2978            int i ;
2979            const double * lower = solver_->getColLower() ;
2980            const double * upper = solver_->getColUpper() ;
2981            for (i = 0 ; i < numberColumns ; i++) {
2982              lowerBefore2[i]= lower[i] ;
2983              upperBefore2[i]= upper[i] ;
2984            }
2985            int numberIterations;
2986            int numberSolves;
2987            int numberGood = solver_->solveBranches(numberBranches,branches,result,
2988                                                    numberSolves,numberIterations);
2989            numberIterations_ += numberIterations;
2990            // Add in solves needed to get to leaves
2991            numberNodes_ += (numberSolves-numberGood);
2992            int numberCutsAdd=-1; // We will be taking off one anyway
2993            int saveMini = sizeMiniTree_;
2994            sizeMiniTree_=0;
2995            int saveNumber = numberNodesOutput;
2996            for (int iResult=0;iResult<numberGood;iResult++) {
2997              if (result[iResult].objectiveValue()<1.0e50) {
2998                // put in correct bounds etc
2999                solver_->setColLower(lowerBefore2);
3000                solver_->setColUpper(upperBefore2);
3001                result[iResult].restoreResult(*solver_);
3002                numberNodes_++;
3003#if 0
3004                // need to tell solver that current factorization is no good
3005                // should be better way
3006#ifdef COIN_USE_CLP
3007                OsiClpSolverInterface * clpSolver
3008                  = dynamic_cast<OsiClpSolverInterface *> (solver_);
3009                if (clpSolver)
3010                  clpSolver->getModelPtr()->setWhatsChanged(clpSolver->getModelPtr()->whatsChanged()&(~512));
3011#endif
3012                solver_->resolve();
3013                assert (!solver_->getIterationCount());
3014                printf("res obj %g was %g\n",solver_->getObjValue(),
3015                       result[iResult].objectiveValue());
3016#endif
3017                bool checkingNode=false;
3018                bool feasible = true;
3019                int anyAction;
3020                CbcNode * newNode2= new CbcNode ;
3021                newNode2->setObjectiveValue(result[iResult].objectiveValue()) ;
3022                anyAction =-1 ;
3023                // only allow twenty passes
3024                int numberPassesLeft=20;
3025                checkingNode=true;
3026                OsiSolverBranch * branches=NULL;
3027                while (anyAction == -1) { 
3028                  if (numberBeforeTrust_==0 ) {
3029                    anyAction = newNode2->chooseBranch(this,node,numberPassesLeft) ;
3030                  } else {
3031                    anyAction = newNode2->chooseDynamicBranch(this,node,branches,numberPassesLeft) ;
3032                    if (anyAction==-3) 
3033                      anyAction = newNode2->chooseBranch(this,node,numberPassesLeft) ; // dynamic did nothing
3034                  }
3035                  numberPassesLeft--;
3036                  if (anyAction == -1) {
3037                    // can do quick optimality check
3038                    int easy=2;
3039                    solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
3040                    feasible = resolve() ;
3041                    solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3042                    if (problemFeasibility_->feasible(this,0)<0) {
3043                      feasible=false; // pretend infeasible
3044                    }
3045                    if (feasible) {
3046                      // Set objective value (not so obvious if NLP etc)
3047                      setObjectiveValue(newNode2,node);
3048                      reducedCostFix() ;
3049                      if (newNode2->objectiveValue() >= getCutoff()) 
3050                        anyAction=-2;
3051                    } else {
3052                      anyAction = -2 ;
3053                    }
3054                  }
3055                }
3056                // May have slipped through i.e. anyAction == 0 and objective above cutoff
3057                if ( anyAction ==0 ) {
3058                  assert (newNode2);
3059                  if (newNode2->objectiveValue() >= getCutoff()) 
3060                    anyAction = -2; // say bad after all
3061                }
3062                if (anyAction >= 0) {
3063                  if (anyAction==0||!newNode2->numberUnsatisfied()) {
3064                    newNode2->createInfo(this,node,lastws,lowerBefore,upperBefore,
3065                                        numberOldActiveCuts_,numberNewCuts_) ;
3066                    if (newNode2->numberUnsatisfied()) {
3067                      newNode2->initializeInfo() ;
3068                      newNode2->nodeInfo()->addCuts(cuts,newNode2->numberBranches(),
3069                                                   whichGenerator_) ;
3070                    }
3071                  }
3072                }
3073                /*
3074                  If we end up infeasible, we can delete the new node immediately. Since this
3075                  node won't be needing the cuts we collected, decrement the reference counts.
3076                  If we are feasible, then we'll be placing this node into the live set, so
3077                  increment the reference count in the current (parent) nodeInfo.
3078                */
3079                if (anyAction == -2) { 
3080                  delete newNode2 ;
3081                  newNode2 = NULL ;
3082#if 0
3083                  // say strong doing well
3084                  if (checkingNode)
3085                    setSpecialOptions(specialOptions_|8);
3086#endif
3087#if 0
3088                  for (i = 0 ; i < currentNumberCuts_ ; i++) {
3089                    if (addedCuts_[i]) {
3090                      if (!addedCuts_[i]->decrement(1)) {
3091                        delete addedCuts_[i] ;
3092                        addedCuts_[i]=NULL;
3093                      }
3094                    }
3095                  }
3096#endif
3097                } else {
3098                  nodeInfo->increment() ;
3099#if 0
3100                  if ((numberNodes_%20)==0) {
3101                    // say strong not doing as well
3102                    setSpecialOptions(specialOptions_&~8);
3103                  }
3104#endif
3105                }
3106                /*
3107                  At this point, there are three possibilities:
3108                  * We have a live node (variable() >= 0) which will require further
3109                  branching to resolve. Before we push it onto the search tree, try for
3110                  a heuristic solution.
3111                  * We have a solution, in which case newNode is non-null but we have no
3112                  branching variable. Decrement the cut counts and save the solution.
3113                  * The node was found to be infeasible, in which case it's already been
3114                  deleted, and newNode is null.
3115                 
3116                */
3117#ifdef COIN_USE_CLP
3118                if (eventHandler()) {
3119                  if (!eventHandler()->event(ClpEventHandler::node)) {
3120                    eventHappened_=true; // exit
3121                  }
3122                }
3123#endif
3124                assert (!newNode2 || newNode2->objectiveValue() <= getCutoff()) ;
3125                if (statistics_) {
3126                  assert (numberNodes2_);
3127                  assert (statistics_[numberNodes2_-1]);
3128                  assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
3129                  if (newNode2)
3130                    statistics_[numberNodes2_-1]->updateInfeasibility(newNode2->numberUnsatisfied());
3131                  else
3132                    statistics_[numberNodes2_-1]->sayInfeasible();
3133                }
3134                if (newNode2) {
3135                  if (newNode2->variable() >= 0) {
3136                    handler_->message(CBC_BRANCH,messages_)
3137                      << numberNodes_<< newNode2->objectiveValue()
3138                      << newNode2->numberUnsatisfied()<< newNode2->depth()
3139                      << CoinMessageEol ;
3140                    // Increment cut counts
3141                    numberCutsAdd += newNode2->numberBranches();
3142                    double estValue = newNode2->guessedObjectiveValue() ;
3143                    int found = -1 ;
3144                    double * newSolution = new double [numberColumns] ;
3145                    double heurValue = getCutoff() ;
3146                    int iHeur ;
3147                    for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++) {
3148                      double saveValue = heurValue ;
3149                      int ifSol = heuristic_[iHeur]->solution(heurValue,newSolution) ;
3150                      if (ifSol > 0) {
3151                        // new solution found
3152                        found = iHeur ;
3153                        incrementUsed(newSolution);
3154                      } else if (ifSol < 0) {   // just returning an estimate
3155                        estValue = CoinMin(heurValue,estValue) ;
3156                        heurValue = saveValue ;
3157                      } 
3158                    }
3159                    if (found >= 0) {
3160                      setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
3161                      lastHeuristic_ = heuristic_[found];
3162                      status=1;
3163                    }
3164                    delete [] newSolution ;
3165                    newNode2->setGuessedObjectiveValue(estValue) ;
3166                    newNodes[numberNodesOutput++] = newNode2 ;
3167                    nodeInfo->increment();
3168                    if (statistics_) {
3169                      if (numberNodes2_==maximumStatistics_) {
3170                        maximumStatistics_ = 2*maximumStatistics_;
3171                        CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
3172                        memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
3173                        memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
3174                        delete [] statistics_;
3175                        statistics_=temp;
3176                      }
3177                      assert (!statistics_[numberNodes2_]);
3178                      statistics_[numberNodes2_]=new CbcStatistics(newNode2);
3179                    }
3180                    numberNodes2_++;
3181#           ifdef CHECK_NODE
3182                    printf("Node %x pushed on tree c\n",newNode2) ;
3183#           endif
3184                  } else {
3185#if 0
3186                    for (i = 0 ; i < currentNumberCuts_ ; i++) {
3187                      if (addedCuts_[i]) {
3188                        if (!addedCuts_[i]->decrement(1)) {
3189                          delete addedCuts_[i] ;
3190                          addedCuts_[i]=NULL;
3191                        }
3192                      }
3193                    }
3194#endif
3195                    double objectiveValue = newNode2->objectiveValue();
3196                    setBestSolution(CBC_SOLUTION,objectiveValue,
3197                                    solver_->getColSolution()) ;
3198                    lastHeuristic_ = NULL;
3199                    incrementUsed(solver_->getColSolution());
3200                    //assert(nodeInfo->numberPointingToThis() <= 2) ;
3201                    // avoid accidental pruning, if newNode was final branch arm
3202                    nodeInfo->increment();
3203                    delete newNode2 ;
3204                    nodeInfo->decrement() ;
3205                    status=1;
3206                  }
3207                }
3208              }
3209            }
3210            if (numberCutsAdd>0) {
3211              for (i = 0;i < currentNumberCuts_;i++) {
3212                if (addedCuts_[i]) {
3213#               ifdef CHECK_CUT_COUNTS
3214                  printf("Count on cut %x increased by %d\n",addedCuts_[i],
3215                         numberCutsAdd) ;
3216#               endif
3217                  addedCuts_[i]->increment(numberCutsAdd) ;
3218                }
3219              }
3220            }     
3221            sizeMiniTree_=saveMini;
3222            delete [] result;
3223            delete [] branches;
3224            delete newNode ;
3225            newNode = NULL ;
3226            delete [] lowerBefore2;
3227            delete [] upperBefore2;
3228            if (saveNumber < numberNodesOutput)
3229              anyAction = -4;
3230            else
3231              anyAction=-2;
3232          }
3233        }
3234      } else {
3235        anyAction = -2 ;
3236      }
3237#if 0
3238      if (numberNodesOutput) {
3239        if (!sizeMiniTree_)
3240          assert(newNode== newNodes[numberNodesOutput-1]);
3241        newNode= newNodes[numberNodesOutput-1];
3242      } else {
3243        assert (!newNode);
3244      }
3245#endif
3246      // May have slipped through i.e. anyAction == 0 and objective above cutoff
3247      if ( anyAction >=0 ) {
3248        assert (newNode);
3249        if (newNode->objectiveValue() >= getCutoff()) 
3250          anyAction = -2; // say bad after all
3251      }
3252      /*
3253        If we end up infeasible, we can delete the new node immediately. Since this
3254        node won't be needing the cuts we collected, decrement the reference counts.
3255        If we are feasible, then we'll be placing this node into the live set, so
3256        increment the reference count in the current (parent) nodeInfo.
3257      */
3258      if (anyAction == -2) { 
3259        delete newNode ;
3260        newNode = NULL ;
3261        // say strong doing well
3262        if (checkingNode)
3263          setSpecialOptions(specialOptions_|8);
3264        for (i = 0 ; i < currentNumberCuts_ ; i++) {
3265          if (addedCuts_[i]) {
3266            if (!addedCuts_[i]->decrement(1)) {
3267              delete addedCuts_[i] ;
3268              addedCuts_[i]=NULL;
3269            }
3270          }
3271        }
3272      } else if (anyAction>=0) {
3273        nodeInfo->increment() ;
3274        if ((numberNodes_%20)==0) {
3275          // say strong not doing as well
3276          setSpecialOptions(specialOptions_&~8);
3277        }
3278      }
3279      /*
3280        At this point, there are three possibilities:
3281        * We have a live node (variable() >= 0) which will require further
3282        branching to resolve. Before we push it onto the search tree, try for
3283        a heuristic solution.
3284        * We have a solution, in which case newNode is non-null but we have no
3285        branching variable. Decrement the cut counts and save the solution.
3286        * The node was found to be infeasible, in which case it's already been
3287        deleted, and newNode is null.
3288       
3289      */
3290#ifdef COIN_USE_CLP
3291      if (eventHandler()) {
3292        if (!eventHandler()->event(ClpEventHandler::node)) {
3293          eventHappened_=true; // exit
3294        }
3295      }
3296#endif
3297      assert (!newNode || newNode->objectiveValue() <= getCutoff()) ;
3298      if (statistics_) {
3299        assert (numberNodes2_);
3300        assert (statistics_[numberNodes2_-1]);
3301        assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
3302        if (newNode)
3303          statistics_[numberNodes2_-1]->updateInfeasibility(newNode->numberUnsatisfied());
3304        else
3305          statistics_[numberNodes2_-1]->sayInfeasible();
3306      }
3307      if (newNode) {
3308        if (newNode->variable() >= 0) {
3309          handler_->message(CBC_BRANCH,messages_)
3310            << numberNodes_<< newNode->objectiveValue()
3311            << newNode->numberUnsatisfied()<< newNode->depth()
3312            << CoinMessageEol ;
3313          // Increment cut counts (taking off current)
3314          int numberLeft = newNode->numberBranches() ;
3315          for (i = 0;i < currentNumberCuts_;i++) {
3316            if (addedCuts_[i]) {
3317#               ifdef CHECK_CUT_COUNTS
3318              printf("Count on cut %x increased by %d\n",addedCuts_[i],
3319                     numberLeft-1) ;
3320#               endif
3321              addedCuts_[i]->increment(numberLeft-1) ;
3322            }
3323          }
3324         
3325          double estValue = newNode->guessedObjectiveValue() ;
3326          int found = -1 ;
3327          double * newSolution = new double [numberColumns] ;
3328          double heurValue = getCutoff() ;
3329          int iHeur ;
3330          for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++) {
3331            double saveValue = heurValue ;
3332            int ifSol = heuristic_[iHeur]->solution(heurValue,newSolution) ;
3333            if (ifSol > 0) {
3334              // new solution found
3335              found = iHeur ;
3336              incrementUsed(newSolution);
3337            } else if (ifSol < 0) {     // just returning an estimate
3338              estValue = CoinMin(heurValue,estValue) ;
3339              heurValue = saveValue ;
3340            } 
3341          }
3342          if (found >= 0) {
3343            setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
3344            lastHeuristic_ = heuristic_[found];
3345            status=1;
3346          }
3347          delete [] newSolution ;
3348          newNode->setGuessedObjectiveValue(estValue) ;
3349          newNodes[numberNodesOutput++] = newNode ;
3350          if (statistics_) {
3351            if (numberNodes2_==maximumStatistics_) {
3352              maximumStatistics_ = 2*maximumStatistics_;
3353              CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
3354              memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
3355              memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
3356              delete [] statistics_;
3357              statistics_=temp;
3358            }
3359            assert (!statistics_[numberNodes2_]);
3360            statistics_[numberNodes2_]=new CbcStatistics(newNode);
3361          }
3362          numberNodes2_++;
3363#           ifdef CHECK_NODE
3364          printf("Node %x pushed on tree c\n",newNode) ;
3365#           endif
3366        } else {
3367          for (i = 0 ; i < currentNumberCuts_ ; i++) {
3368            if (addedCuts_[i]) {
3369              if (!addedCuts_[i]->decrement(1)) {
3370                delete addedCuts_[i] ;
3371                addedCuts_[i]=NULL;
3372              }
3373            }
3374          }
3375          double objectiveValue = newNode->objectiveValue();
3376          setBestSolution(CBC_SOLUTION,objectiveValue,
3377                          solver_->getColSolution()) ;
3378          lastHeuristic_ = NULL;
3379          incrementUsed(solver_->getColSolution());
3380          //assert(nodeInfo->numberPointingToThis() <= 2) ;
3381          // avoid accidental pruning, if newNode was final branch arm
3382          nodeInfo->increment();
3383          delete newNode ;
3384          nodeInfo->decrement() ;
3385          status=1;
3386        }
3387      }
3388      /*
3389        This node has been completely expanded and can be removed from the live
3390        set.
3391      */
3392      if (deleteNode) 
3393        delete node ;
3394    } else {
3395      // max time etc
3396      status=2;
3397      if (deleteNode) 
3398        delete node ;
3399    }
3400  } else {
3401    // node before branching is infeasible (so no nodes output)
3402    delete node;
3403  }
3404  delete lastws;
3405  delete [] lowerBefore;
3406  delete [] upperBefore;
3407  return newNodes;
3408}
3409#endif
3410
3411
3412// Solve the initial LP relaxation
3413void 
3414CbcModel::initialSolve() 
3415{
3416  assert (solver_);
3417  assert (!solverCharacteristics_);
3418  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
3419  if (solverCharacteristics) {
3420    solverCharacteristics_ = solverCharacteristics;
3421  } else {
3422    // replace in solver
3423    OsiBabSolver defaultC;
3424    solver_->setAuxiliaryInfo(&defaultC);
3425    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
3426  }
3427  solverCharacteristics_->setSolver(solver_);
3428  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3429  solver_->initialSolve();
3430  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
3431  // But set up so Jon Lee will be happy
3432  status_=-1;
3433  secondaryStatus_ = -1;
3434  originalContinuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
3435  delete [] continuousSolution_;
3436  continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
3437                                             solver_->getNumCols());
3438  setPointers(solver_);
3439  solverCharacteristics_ = NULL;
3440}
3441
3442/*! \brief Get an empty basis object
3443
3444  Return an empty CoinWarmStartBasis object with the requested capacity,
3445  appropriate for the current solver. The object is cloned from the object
3446  cached as emptyWarmStart_. If there is no cached object, the routine
3447  queries the solver for a warm start object, empties it, and caches the
3448  result.
3449*/
3450
3451CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
3452
3453{ CoinWarmStartBasis *emptyBasis ;
3454/*
3455  Acquire an empty basis object, if we don't yet have one.
3456*/
3457  if (emptyWarmStart_ == 0)
3458  { if (solver_ == 0)
3459    { throw CoinError("Cannot construct basis without solver!",
3460                      "getEmptyBasis","CbcModel") ; }
3461    emptyBasis =
3462        dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
3463    if (emptyBasis == 0)
3464    { throw CoinError(
3465        "Solver does not appear to use a basis-oriented warm start.",
3466        "getEmptyBasis","CbcModel") ; }
3467    emptyBasis->setSize(0,0) ;
3468    emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ; }
3469/*
3470  Clone the empty basis object, resize it as requested, and return.
3471*/
3472  emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
3473  assert(emptyBasis) ;
3474  if (ns != 0 || na != 0) emptyBasis->setSize(ns,na) ;
3475
3476  return (emptyBasis) ; }
3477   
3478
3479/** Default Constructor
3480
3481  Creates an empty model without an associated solver.
3482*/
3483CbcModel::CbcModel() 
3484
3485:
3486  solver_(NULL),
3487  ourSolver_(true),
3488  continuousSolver_(NULL),
3489  referenceSolver_(NULL),
3490  defaultHandler_(true),
3491  emptyWarmStart_(NULL),
3492  bestObjective_(COIN_DBL_MAX),
3493  bestPossibleObjective_(COIN_DBL_MAX),
3494  sumChangeObjective1_(0.0),
3495  sumChangeObjective2_(0.0),
3496  bestSolution_(NULL),
3497  currentSolution_(NULL),
3498  testSolution_(NULL),
3499  minimumDrop_(1.0e-4),
3500  numberSolutions_(0),
3501  stateOfSearch_(0),
3502  hotstartSolution_(NULL),
3503  hotstartPriorities_(NULL),
3504  numberHeuristicSolutions_(0),
3505  numberNodes_(0),
3506  numberNodes2_(0),
3507  numberIterations_(0),
3508  status_(-1),
3509  secondaryStatus_(-1),
3510  numberIntegers_(0),
3511  numberRowsAtContinuous_(0),
3512  maximumNumberCuts_(0),
3513  phase_(0),
3514  currentNumberCuts_(0),
3515  maximumDepth_(0),
3516  walkback_(NULL),
3517  addedCuts_(NULL),
3518  nextRowCut_(NULL),
3519  currentNode_(NULL),
3520  integerVariable_(NULL),
3521  integerInfo_(NULL),
3522  continuousSolution_(NULL),
3523  usedInSolution_(NULL),
3524  specialOptions_(0),
3525  subTreeModel_(NULL),
3526  numberStoppedSubTrees_(0),
3527  presolve_(0),
3528  numberStrong_(5),
3529  numberBeforeTrust_(0),
3530  numberPenalties_(20),
3531  penaltyScaleFactor_(3.0),
3532  numberAnalyzeIterations_(0),
3533  analyzeResults_(NULL),
3534  numberInfeasibleNodes_(0),
3535  problemType_(0),
3536  printFrequency_(0),
3537  numberCutGenerators_(0),
3538  generator_(NULL),
3539  virginGenerator_(NULL),
3540  numberHeuristics_(0),
3541  heuristic_(NULL),
3542  lastHeuristic_(NULL),
3543  numberObjects_(0),
3544  object_(NULL),
3545  originalColumns_(NULL),
3546  howOftenGlobalScan_(1),
3547  numberGlobalViolations_(0),
3548  continuousObjective_(COIN_DBL_MAX),
3549  originalContinuousObjective_(COIN_DBL_MAX),
3550  continuousInfeasibilities_(INT_MAX),
3551  maximumCutPassesAtRoot_(20),
3552  maximumCutPasses_(10),
3553  currentPassNumber_(0),
3554  maximumWhich_(1000),
3555  whichGenerator_(NULL),
3556  maximumStatistics_(0),
3557  statistics_(NULL),
3558  numberFixedAtRoot_(0),
3559  numberFixedNow_(0),
3560  stoppedOnGap_(false),
3561  eventHappened_(false),
3562  numberLongStrong_(0),
3563  numberOldActiveCuts_(0),
3564  numberNewCuts_(0),
3565  sizeMiniTree_(0),
3566  searchStrategy_(-1),
3567  numberStrongIterations_(0),
3568  resolveAfterTakeOffCuts_(true)
3569{
3570  intParam_[CbcMaxNumNode] = 2147483647;
3571  intParam_[CbcMaxNumSol] = 9999999;
3572  intParam_[CbcFathomDiscipline] = 0;
3573
3574  dblParam_[CbcIntegerTolerance] = 1e-6;
3575  dblParam_[CbcInfeasibilityWeight] = 0.0;
3576  dblParam_[CbcCutoffIncrement] = 1e-5;
3577  dblParam_[CbcAllowableGap] = 1.0e-10;
3578  dblParam_[CbcAllowableFractionGap] = 0.0;
3579  dblParam_[CbcMaximumSeconds] = 1.0e100;
3580  dblParam_[CbcCurrentCutoff] = 1.0e100;
3581  dblParam_[CbcOptimizationDirection] = 1.0;
3582  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
3583  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
3584  dblParam_[CbcStartSeconds] = 0.0;
3585  strongInfo_[0]=0;
3586  strongInfo_[1]=0;
3587  strongInfo_[2]=0;
3588  solverCharacteristics_ = NULL;
3589  nodeCompare_=new CbcCompareDefault();;
3590  problemFeasibility_=new CbcFeasibilityBase();
3591  tree_= new CbcTree();
3592  branchingMethod_=NULL;
3593  strategy_=NULL;
3594  parentModel_=NULL;
3595  cbcColLower_ = NULL;
3596  cbcColUpper_ = NULL;
3597  cbcRowLower_ = NULL;
3598  cbcRowUpper_ = NULL;
3599  cbcColSolution_ = NULL;
3600  cbcRowPrice_ = NULL;
3601  cbcReducedCost_ = NULL;
3602  cbcRowActivity_ = NULL;
3603  appData_=NULL;
3604  handler_ = new CoinMessageHandler();
3605  handler_->setLogLevel(2);
3606  messages_ = CbcMessage();
3607}
3608
3609/** Constructor from solver.
3610
3611  Creates a model complete with a clone of the solver passed as a parameter.
3612*/
3613
3614CbcModel::CbcModel(const OsiSolverInterface &rhs)
3615:
3616  continuousSolver_(NULL),
3617  referenceSolver_(NULL),
3618  defaultHandler_(true),
3619  emptyWarmStart_(NULL),
3620  bestObjective_(COIN_DBL_MAX),
3621  bestPossibleObjective_(COIN_DBL_MAX),
3622  sumChangeObjective1_(0.0),
3623  sumChangeObjective2_(0.0),
3624  minimumDrop_(1.0e-4),
3625  numberSolutions_(0),
3626  stateOfSearch_(0),
3627  hotstartSolution_(NULL),
3628  hotstartPriorities_(NULL),
3629  numberHeuristicSolutions_(0),
3630  numberNodes_(0),
3631  numberNodes2_(0),
3632  numberIterations_(0),
3633  status_(-1),
3634  secondaryStatus_(-1),
3635  numberRowsAtContinuous_(0),
3636  maximumNumberCuts_(0),
3637  phase_(0),
3638  currentNumberCuts_(0),
3639  maximumDepth_(0),
3640  walkback_(NULL),
3641  addedCuts_(NULL),
3642  nextRowCut_(NULL),
3643  currentNode_(NULL),
3644  integerInfo_(NULL),
3645  specialOptions_(0),
3646  subTreeModel_(NULL),
3647  numberStoppedSubTrees_(0),
3648  presolve_(0),
3649  numberStrong_(5),
3650  numberBeforeTrust_(0),
3651  numberPenalties_(20),
3652  penaltyScaleFactor_(3.0),
3653  numberAnalyzeIterations_(0),
3654  analyzeResults_(NULL),
3655  numberInfeasibleNodes_(0),
3656  problemType_(0),
3657  printFrequency_(0),
3658  numberCutGenerators_(0),
3659  generator_(NULL),
3660  virginGenerator_(NULL),
3661  numberHeuristics_(0),
3662  heuristic_(NULL),
3663  lastHeuristic_(NULL),
3664  numberObjects_(0),
3665  object_(NULL),
3666  originalColumns_(NULL),
3667  howOftenGlobalScan_(1),
3668  numberGlobalViolations_(0),
3669  continuousObjective_(COIN_DBL_MAX),
3670  originalContinuousObjective_(COIN_DBL_MAX),
3671  continuousInfeasibilities_(INT_MAX),
3672  maximumCutPassesAtRoot_(20),
3673  maximumCutPasses_(10),
3674  currentPassNumber_(0),
3675  maximumWhich_(1000),
3676  whichGenerator_(NULL),
3677  maximumStatistics_(0),
3678  statistics_(NULL),
3679  numberFixedAtRoot_(0),
3680  numberFixedNow_(0),
3681  stoppedOnGap_(false),
3682  eventHappened_(false),
3683  numberLongStrong_(0),
3684  numberOldActiveCuts_(0),
3685  numberNewCuts_(0),
3686  sizeMiniTree_(0),
3687  searchStrategy_(-1),
3688  numberStrongIterations_(0),
3689  resolveAfterTakeOffCuts_(true)
3690{
3691  intParam_[CbcMaxNumNode] = 2147483647;
3692  intParam_[CbcMaxNumSol] = 9999999;
3693  intParam_[CbcFathomDiscipline] = 0;
3694
3695  dblParam_[CbcIntegerTolerance] = 1e-6;
3696  dblParam_[CbcInfeasibilityWeight] = 0.0;
3697  dblParam_[CbcCutoffIncrement] = 1e-5;
3698  dblParam_[CbcAllowableGap] = 1.0e-10;
3699  dblParam_[CbcAllowableFractionGap] = 0.0;
3700  dblParam_[CbcMaximumSeconds] = 1.0e100;
3701  dblParam_[CbcCurrentCutoff] = 1.0e100;
3702  dblParam_[CbcOptimizationDirection] = 1.0;
3703  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
3704  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
3705  dblParam_[CbcStartSeconds] = 0.0;
3706  strongInfo_[0]=0;
3707  strongInfo_[1]=0;
3708  strongInfo_[2]=0;
3709  solverCharacteristics_ = NULL;
3710
3711  nodeCompare_=new CbcCompareDefault();;
3712  problemFeasibility_=new CbcFeasibilityBase();
3713  tree_= new CbcTree();
3714  branchingMethod_=NULL;
3715  strategy_=NULL;
3716  parentModel_=NULL;
3717  appData_=NULL;
3718  handler_ = new CoinMessageHandler();
3719  handler_->setLogLevel(2);
3720  messages_ = CbcMessage();
3721  solver_ = rhs.clone();
3722  referenceSolver_ = solver_->clone();
3723  ourSolver_ = true ;
3724  cbcColLower_ = NULL;
3725  cbcColUpper_ = NULL;
3726  cbcRowLower_ = NULL;
3727  cbcRowUpper_ = NULL;
3728  cbcColSolution_ = NULL;
3729  cbcRowPrice_ = NULL;
3730  cbcReducedCost_ = NULL;
3731  cbcRowActivity_ = NULL;
3732
3733  // Initialize solution and integer variable vectors
3734  bestSolution_ = NULL; // to say no solution found
3735  numberIntegers_=0;
3736  int numberColumns = solver_->getNumCols();
3737  int iColumn;
3738  if (numberColumns) {
3739    // Space for current solution
3740    currentSolution_ = new double[numberColumns];
3741    continuousSolution_ = new double[numberColumns];
3742    usedInSolution_ = new int[numberColumns];
3743    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3744      if( solver_->isInteger(iColumn)) 
3745        numberIntegers_++;
3746    }
3747  } else {
3748    // empty model
3749    currentSolution_=NULL;
3750    continuousSolution_=NULL;
3751    usedInSolution_=NULL;
3752  }
3753  testSolution_=currentSolution_;
3754  if (numberIntegers_) {
3755    integerVariable_ = new int [numberIntegers_];
3756    numberIntegers_=0;
3757    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3758      if( solver_->isInteger(iColumn)) 
3759        integerVariable_[numberIntegers_++]=iColumn;
3760    }
3761  } else {
3762    integerVariable_ = NULL;
3763  }
3764}
3765
3766/*
3767  Assign a solver to the model (model assumes ownership)
3768
3769  The integer variable vector is initialized if it's not already present.
3770  If deleteSolver then current solver deleted (if model owned)
3771
3772  Assuming ownership matches usage in OsiSolverInterface
3773  (cf. assignProblem, loadProblem).
3774
3775  TODO: What to do about solver parameters? A simple copy likely won't do it,
3776        because the SI must push the settings into the underlying solver. In
3777        the context of switching solvers in cbc, this means that command line
3778        settings will get lost. Stash the command line somewhere and reread it
3779        here, maybe?
3780 
3781  TODO: More generally, how much state should be transferred from the old
3782        solver to the new solver? Best perhaps to see how usage develops.
3783        What's done here mimics the CbcModel(OsiSolverInterface) constructor.
3784*/
3785void
3786CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver)
3787
3788{
3789  // Keep the current message level for solver (if solver exists)
3790  if (solver_)
3791    solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
3792
3793  if (ourSolver_&&deleteSolver) delete solver_ ;
3794  solver_ = solver;
3795  solver = NULL ;
3796  ourSolver_ = true ;
3797/*
3798  Basis information is solver-specific.
3799*/
3800  if (emptyWarmStart_)
3801  { delete emptyWarmStart_  ;
3802    emptyWarmStart_ = 0 ; }
3803/*
3804  Initialize integer variable vector.
3805*/
3806  numberIntegers_=0;
3807  int numberColumns = solver_->getNumCols();
3808  int iColumn;
3809  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3810    if( solver_->isInteger(iColumn)) 
3811      numberIntegers_++;
3812  }
3813  if (numberIntegers_) {
3814    delete [] integerVariable_;
3815    integerVariable_ = new int [numberIntegers_];
3816    numberIntegers_=0;
3817    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3818      if( solver_->isInteger(iColumn)) 
3819        integerVariable_[numberIntegers_++]=iColumn;
3820    }
3821  } else {
3822    integerVariable_ = NULL;
3823  }
3824
3825  return ;
3826}
3827
3828// Copy constructor.
3829
3830CbcModel::CbcModel(const CbcModel & rhs, bool noTree)
3831:
3832  continuousSolver_(NULL),
3833  referenceSolver_(NULL),
3834  defaultHandler_(rhs.defaultHandler_),
3835  emptyWarmStart_(NULL),
3836  bestObjective_(rhs.bestObjective_),
3837  bestPossibleObjective_(rhs.bestPossibleObjective_),
3838  sumChangeObjective1_(rhs.sumChangeObjective1_),
3839  sumChangeObjective2_(rhs.sumChangeObjective2_),
3840  minimumDrop_(rhs.minimumDrop_),
3841  numberSolutions_(rhs.numberSolutions_),
3842  stateOfSearch_(rhs.stateOfSearch_),
3843  numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
3844  numberNodes_(rhs.numberNodes_),
3845  numberNodes2_(rhs.numberNodes2_),
3846  numberIterations_(rhs.numberIterations_),
3847  status_(rhs.status_),
3848  secondaryStatus_(rhs.secondaryStatus_),
3849  specialOptions_(rhs.specialOptions_),
3850  subTreeModel_(rhs.subTreeModel_),
3851  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
3852  presolve_(rhs.presolve_),
3853  numberStrong_(rhs.numberStrong_),
3854  numberBeforeTrust_(rhs.numberBeforeTrust_),
3855  numberPenalties_(rhs.numberPenalties_),
3856  penaltyScaleFactor_(rhs.penaltyScaleFactor_),
3857  numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
3858  analyzeResults_(NULL),
3859  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
3860  problemType_(rhs.problemType_),
3861  printFrequency_(rhs.printFrequency_),
3862  howOftenGlobalScan_(rhs.howOftenGlobalScan_),
3863  numberGlobalViolations_(rhs.numberGlobalViolations_),
3864  continuousObjective_(rhs.continuousObjective_),
3865  originalContinuousObjective_(rhs.originalContinuousObjective_),
3866  continuousInfeasibilities_(rhs.continuousInfeasibilities_),
3867  maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
3868  maximumCutPasses_( rhs.maximumCutPasses_),
3869  currentPassNumber_(rhs.currentPassNumber_),
3870  maximumWhich_(rhs.maximumWhich_),
3871  whichGenerator_(NULL),
3872  maximumStatistics_(0),
3873  statistics_(NULL),
3874  numberFixedAtRoot_(rhs.numberFixedAtRoot_),
3875  numberFixedNow_(rhs.numberFixedNow_),
3876  stoppedOnGap_(rhs.stoppedOnGap_),
3877  eventHappened_(rhs.eventHappened_),
3878  numberLongStrong_(rhs.numberLongStrong_),
3879  numberOldActiveCuts_(rhs.numberOldActiveCuts_),
3880  numberNewCuts_(rhs.numberNewCuts_),
3881  sizeMiniTree_(rhs.sizeMiniTree_),
3882  searchStrategy_(rhs.searchStrategy_),
3883  numberStrongIterations_(rhs.numberStrongIterations_),
3884  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_)
3885{
3886  intParam_[CbcMaxNumNode] = rhs.intParam_[CbcMaxNumNode];
3887  intParam_[CbcMaxNumSol] = rhs.intParam_[CbcMaxNumSol];
3888  intParam_[CbcFathomDiscipline] = rhs.intParam_[CbcFathomDiscipline];
3889  dblParam_[CbcIntegerTolerance] = rhs.dblParam_[CbcIntegerTolerance];
3890  dblParam_[CbcInfeasibilityWeight] = rhs.dblParam_[CbcInfeasibilityWeight];
3891  dblParam_[CbcCutoffIncrement] = rhs.dblParam_[CbcCutoffIncrement]; 
3892  dblParam_[CbcAllowableGap] = rhs.dblParam_[CbcAllowableGap]; 
3893  dblParam_[CbcAllowableFractionGap] = rhs.dblParam_[CbcAllowableFractionGap]; 
3894  dblParam_[CbcMaximumSeconds] = rhs.dblParam_[CbcMaximumSeconds];
3895  dblParam_[CbcCurrentCutoff] = rhs.dblParam_[CbcCurrentCutoff];
3896  dblParam_[CbcOptimizationDirection] = rhs.dblParam_[CbcOptimizationDirection];
3897  dblParam_[CbcCurrentObjectiveValue] = rhs.dblParam_[CbcCurrentObjectiveValue];
3898  dblParam_[CbcCurrentMinimizationObjectiveValue] = rhs.dblParam_[CbcCurrentMinimizationObjectiveValue];
3899  dblParam_[CbcStartSeconds] = dblParam_[CbcStartSeconds]; // will be overwritten hopefully
3900  strongInfo_[0]=rhs.strongInfo_[0];
3901  strongInfo_[1]=rhs.strongInfo_[1];
3902  strongInfo_[2]=rhs.strongInfo_[2];
3903  solverCharacteristics_ = NULL;
3904  if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
3905  if (defaultHandler_) {
3906    handler_ = new CoinMessageHandler();
3907    handler_->setLogLevel(2);
3908  } else {
3909    handler_ = rhs.handler_;
3910  }
3911  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3912  numberCutGenerators_ = rhs.numberCutGenerators_;
3913  if (numberCutGenerators_) {
3914    generator_ = new CbcCutGenerator * [numberCutGenerators_];
3915    virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
3916    int i;
3917    for (i=0;i<numberCutGenerators_;i++) {
3918      generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
3919      virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
3920    }
3921  } else {
3922    generator_=NULL;
3923    virginGenerator_=NULL;
3924  }
3925  if (!noTree)
3926    globalCuts_ = rhs.globalCuts_;
3927  numberHeuristics_ = rhs.numberHeuristics_;
3928  if (numberHeuristics_) {
3929    heuristic_ = new CbcHeuristic * [numberHeuristics_];
3930    int i;
3931    for (i=0;i<numberHeuristics_;i++) {
3932      heuristic_[i]=rhs.heuristic_[i]->clone();
3933    }
3934  } else {
3935    heuristic_=NULL;
3936  }
3937  lastHeuristic_ = NULL;
3938  numberObjects_=rhs.numberObjects_;
3939  if (numberObjects_) {
3940    object_ = new CbcObject * [numberObjects_];
3941    int i;
3942    for (i=0;i<numberObjects_;i++) 
3943      object_[i]=(rhs.object_[i])->clone();
3944  } else {
3945    object_=NULL;
3946  }
3947  if (rhs.referenceSolver_)
3948    referenceSolver_ = rhs.referenceSolver_->clone();
3949  else
3950    referenceSolver_=NULL;
3951  if (!noTree||!rhs.continuousSolver_)
3952    solver_ = rhs.solver_->clone();
3953  else
3954    solver_ = rhs.continuousSolver_->clone();
3955  if (rhs.originalColumns_) {
3956    int numberColumns = solver_->getNumCols();
3957    originalColumns_= new int [numberColumns];
3958    memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
3959  } else {
3960    originalColumns_=NULL;
3961  }
3962  if (maximumWhich_&&rhs.whichGenerator_)
3963    whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
3964  nodeCompare_=rhs.nodeCompare_->clone();
3965  problemFeasibility_=rhs.problemFeasibility_->clone();
3966  tree_= rhs.tree_->clone();
3967  branchingMethod_=rhs.branchingMethod_;
3968  cbcColLower_ = NULL;
3969  cbcColUpper_ = NULL;
3970  cbcRowLower_ = NULL;
3971  cbcRowUpper_ = NULL;
3972  cbcColSolution_ = NULL;
3973  cbcRowPrice_ = NULL;
3974  cbcReducedCost_ = NULL;
3975  cbcRowActivity_ = NULL;
3976  if (rhs.strategy_)
3977    strategy_=rhs.strategy_->clone();
3978  else
3979    strategy_=NULL;
3980  parentModel_=rhs.parentModel_;
3981  appData_=rhs.appData_;
3982  messages_ = rhs.messages_;
3983  ourSolver_ = true ;
3984  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3985  numberIntegers_=rhs.numberIntegers_;
3986  if (numberIntegers_) {
3987    integerVariable_ = new int [numberIntegers_];
3988    memcpy(integerVariable_,rhs.integerVariable_,numberIntegers_*sizeof(int));
3989    integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
3990  } else {
3991    integerVariable_ = NULL;
3992    integerInfo_=NULL;
3993  }
3994  if (rhs.hotstartSolution_) {
3995    int numberColumns = solver_->getNumCols();
3996    hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
3997    hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
3998  } else {
3999    hotstartSolution_ = NULL;
4000    hotstartPriorities_ =NULL;
4001  }
4002  if (rhs.bestSolution_&&!noTree) {
4003    int numberColumns = solver_->getNumCols();
4004    bestSolution_ = new double[numberColumns];
4005    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
4006  } else {
4007    bestSolution_=NULL;
4008  }
4009  if (!noTree) {
4010    int numberColumns = solver_->getNumCols();
4011    currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns);
4012    continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns);
4013    usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns);
4014  } else {
4015    currentSolution_=NULL;
4016    continuousSolution_=NULL;
4017    usedInSolution_=NULL;
4018  }
4019  testSolution_=currentSolution_;
4020  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
4021  maximumNumberCuts_=rhs.maximumNumberCuts_;
4022  phase_ = rhs.phase_;
4023  currentNumberCuts_=rhs.currentNumberCuts_;
4024  maximumDepth_= rhs.maximumDepth_;
4025  if (noTree) {
4026    bestObjective_ = COIN_DBL_MAX;
4027    numberSolutions_ =0;
4028    stateOfSearch_= 0;
4029    numberHeuristicSolutions_=0;
4030    numberNodes_=0;
4031    numberNodes2_=0;
4032    numberIterations_=0;
4033    status_=0;
4034    subTreeModel_=NULL;
4035    numberStoppedSubTrees_=0;
4036    continuousObjective_=COIN_DBL_MAX;
4037    originalContinuousObjective_=COIN_DBL_MAX;
4038    continuousInfeasibilities_=INT_MAX;
4039    maximumNumberCuts_=0;
4040    tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_);
4041    bestPossibleObjective_ = COIN_DBL_MAX;
4042  }
4043  // These are only used as temporary arrays so need not be filled
4044  if (maximumNumberCuts_) {
4045    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4046  } else {
4047    addedCuts_ = NULL;
4048  }
4049  nextRowCut_ = NULL;
4050  currentNode_ = NULL;
4051  if (maximumDepth_)
4052    walkback_ = new CbcNodeInfo * [maximumDepth_];
4053  else
4054    walkback_ = NULL;
4055  synchronizeModel();
4056}
4057 
4058// Assignment operator
4059CbcModel & 
4060CbcModel::operator=(const CbcModel& rhs)
4061{
4062  if (this!=&rhs) {
4063    gutsOfDestructor();
4064    if (defaultHandler_)
4065    { delete handler_;
4066      handler_ = NULL; }
4067    defaultHandler_ = rhs.defaultHandler_;
4068    if (defaultHandler_)
4069    { handler_ = new CoinMessageHandler();
4070      handler_->setLogLevel(2); }
4071    else
4072    { handler_ = rhs.handler_; }
4073    messages_ = rhs.messages_;
4074    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4075
4076    delete solver_;
4077    if (rhs.solver_)
4078    { solver_ = rhs.solver_->clone() ; }
4079    else
4080    { solver_ = 0 ; }
4081    ourSolver_ = true ;
4082    delete continuousSolver_ ;
4083    if (rhs.continuousSolver_)
4084    { solver_ = rhs.continuousSolver_->clone() ; }
4085    else
4086    { continuousSolver_ = 0 ; }
4087
4088    if (rhs.referenceSolver_)
4089    { solver_ = rhs.referenceSolver_->clone() ; }
4090    else
4091    { referenceSolver_ = NULL ; }
4092
4093    delete emptyWarmStart_ ;
4094    if (rhs.emptyWarmStart_)
4095    { emptyWarmStart_ = rhs.emptyWarmStart_->clone() ; }
4096    else
4097    { emptyWarmStart_ = 0 ; }
4098
4099    bestObjective_ = rhs.bestObjective_;
4100    bestPossibleObjective_=rhs.bestPossibleObjective_;
4101    sumChangeObjective1_=rhs.sumChangeObjective1_;
4102    sumChangeObjective2_=rhs.sumChangeObjective2_;
4103    delete [] bestSolution_;
4104    if (rhs.bestSolution_) {
4105      int numberColumns = rhs.getNumCols();
4106      bestSolution_ = new double[numberColumns];
4107      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
4108    } else {
4109      bestSolution_=NULL;
4110    }
4111    int numberColumns = solver_->getNumCols();
4112    currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns);
4113    continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns);
4114    usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns);
4115    testSolution_=currentSolution_;
4116    minimumDrop_ = rhs.minimumDrop_;
4117    numberSolutions_=rhs.numberSolutions_;
4118    stateOfSearch_= rhs.stateOfSearch_;
4119    numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
4120    numberNodes_ = rhs.numberNodes_;
4121    numberNodes2_ = rhs.numberNodes2_;
4122    numberIterations_ = rhs.numberIterations_;
4123    status_ = rhs.status_;
4124    secondaryStatus_ = rhs.secondaryStatus_;
4125    specialOptions_ = rhs.specialOptions_;
4126    subTreeModel_ = rhs.subTreeModel_;
4127    numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
4128    presolve_ = rhs.presolve_;
4129    numberStrong_ = rhs.numberStrong_;
4130    numberBeforeTrust_ = rhs.numberBeforeTrust_;
4131    numberPenalties_ = rhs.numberPenalties_;
4132    penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
4133    numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
4134    delete [] analyzeResults_;
4135    analyzeResults_ = NULL;
4136    numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
4137    problemType_ = rhs.problemType_;
4138    printFrequency_ = rhs.printFrequency_;
4139    howOftenGlobalScan_=rhs.howOftenGlobalScan_;
4140    numberGlobalViolations_=rhs.numberGlobalViolations_;
4141    continuousObjective_=rhs.continuousObjective_;
4142    originalContinuousObjective_ = rhs.originalContinuousObjective_;
4143    continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
4144    maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
4145    maximumCutPasses_ = rhs.maximumCutPasses_;
4146    currentPassNumber_ = rhs.currentPassNumber_;
4147    intParam_[CbcMaxNumNode] = rhs.intParam_[CbcMaxNumNode];
4148    intParam_[CbcMaxNumSol] = rhs.intParam_[CbcMaxNumSol];
4149    intParam_[CbcFathomDiscipline] = rhs.intParam_[CbcFathomDiscipline];
4150    dblParam_[CbcIntegerTolerance] = rhs.dblParam_[CbcIntegerTolerance];
4151    dblParam_[CbcInfeasibilityWeight] = rhs.dblParam_[CbcInfeasibilityWeight];
4152    dblParam_[CbcCutoffIncrement] = rhs.dblParam_[CbcCutoffIncrement]; 
4153    dblParam_[CbcAllowableGap] = rhs.dblParam_[CbcAllowableGap]; 
4154    dblParam_[CbcAllowableFractionGap] = rhs.dblParam_[CbcAllowableFractionGap]; 
4155    dblParam_[CbcMaximumSeconds] = rhs.dblParam_[CbcMaximumSeconds];
4156    dblParam_[CbcCurrentCutoff] = rhs.dblParam_[CbcCurrentCutoff];
4157    dblParam_[CbcOptimizationDirection] = rhs.dblParam_[CbcOptimizationDirection];
4158    dblParam_[CbcCurrentObjectiveValue] = rhs.dblParam_[CbcCurrentObjectiveValue];
4159    dblParam_[CbcCurrentMinimizationObjectiveValue] = rhs.dblParam_[CbcCurrentMinimizationObjectiveValue];
4160    dblParam_[CbcStartSeconds] = dblParam_[CbcStartSeconds]; // will be overwritten hopefully
4161    globalCuts_ = rhs.globalCuts_;
4162    int i;
4163    for (i=0;i<numberCutGenerators_;i++) {
4164      delete generator_[i];
4165      delete virginGenerator_[i];
4166    }
4167    delete [] generator_;
4168    delete [] virginGenerator_;
4169    delete [] heuristic_;
4170    maximumWhich_ = rhs.maximumWhich_;
4171    delete [] whichGenerator_;
4172    whichGenerator_ = NULL;
4173    if (maximumWhich_&&rhs.whichGenerator_)
4174      whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
4175    for (i=0;i<maximumStatistics_;i++)
4176      delete statistics_[i];
4177    delete [] statistics_;
4178    maximumStatistics_ = 0;
4179    statistics_ = NULL;
4180    numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
4181    numberFixedNow_ = rhs.numberFixedNow_;
4182    stoppedOnGap_ = rhs.stoppedOnGap_;
4183    eventHappened_ = rhs.eventHappened_;
4184    numberLongStrong_ = rhs.numberLongStrong_;
4185    numberOldActiveCuts_ = rhs.numberOldActiveCuts_;
4186    numberNewCuts_ = rhs.numberNewCuts_;
4187    resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_;
4188    sizeMiniTree_ = rhs.sizeMiniTree_;
4189    searchStrategy_ = rhs.searchStrategy_;
4190    numberStrongIterations_ = rhs.numberStrongIterations_;
4191    strongInfo_[0]=rhs.strongInfo_[0];
4192    strongInfo_[1]=rhs.strongInfo_[1];
4193    strongInfo_[2]=rhs.strongInfo_[2];
4194    solverCharacteristics_ = NULL;
4195    lastHeuristic_ = NULL;
4196    numberCutGenerators_ = rhs.numberCutGenerators_;
4197    if (numberCutGenerators_) {
4198      generator_ = new CbcCutGenerator * [numberCutGenerators_];
4199      virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
4200      int i;
4201      for (i=0;i<numberCutGenerators_;i++) {
4202        generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
4203        virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
4204      }
4205    } else {
4206      generator_=NULL;
4207      virginGenerator_=NULL;
4208    }
4209    numberHeuristics_ = rhs.numberHeuristics_;
4210    if (numberHeuristics_) {
4211      heuristic_ = new CbcHeuristic * [numberHeuristics_];
4212      memcpy(heuristic_,rhs.heuristic_,
4213             numberHeuristics_*sizeof(CbcHeuristic *));
4214    } else {
4215      heuristic_=NULL;
4216    }
4217    lastHeuristic_ = NULL;
4218    for (i=0;i<numberObjects_;i++)
4219      delete object_[i];
4220    delete [] object_;
4221    numberObjects_=rhs.numberObjects_;
4222    if (numberObjects_) {
4223      object_ = new CbcObject * [numberObjects_];
4224      int i;
4225      for (i=0;i<numberObjects_;i++) 
4226        object_[i]=(rhs.object_[i])->clone();
4227    } else {
4228      object_=NULL;
4229    }
4230    delete [] originalColumns_;
4231    if (rhs.originalColumns_) {
4232      int numberColumns = rhs.getNumCols();
4233      originalColumns_= new int [numberColumns];
4234      memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
4235    } else {
4236      originalColumns_=NULL;
4237    }
4238    nodeCompare_=rhs.nodeCompare_->clone();
4239    problemFeasibility_=rhs.problemFeasibility_->clone();
4240    delete tree_;
4241    tree_= rhs.tree_->clone();
4242    branchingMethod_=rhs.branchingMethod_;
4243    delete strategy_;
4244    if (rhs.strategy_)
4245      strategy_=rhs.strategy_->clone();
4246    else
4247      strategy_=NULL;
4248    parentModel_=rhs.parentModel_;
4249    appData_=rhs.appData_;
4250
4251    delete [] integerVariable_;
4252    numberIntegers_=rhs.numberIntegers_;
4253    if (numberIntegers_) {
4254      integerVariable_ = new int [numberIntegers_];
4255      memcpy(integerVariable_,rhs.integerVariable_,
4256             numberIntegers_*sizeof(int));
4257      integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
4258    } else {
4259      integerVariable_ = NULL;
4260      integerInfo_=NULL;
4261    }
4262    if (rhs.hotstartSolution_) {
4263      int numberColumns = solver_->getNumCols();
4264      hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
4265      hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
4266    } else {
4267      hotstartSolution_ = NULL;
4268      hotstartPriorities_ =NULL;
4269    }
4270    numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
4271    maximumNumberCuts_=rhs.maximumNumberCuts_;
4272    phase_ = rhs.phase_;
4273    currentNumberCuts_=rhs.currentNumberCuts_;
4274    maximumDepth_= rhs.maximumDepth_;
4275    delete [] addedCuts_;
4276    delete [] walkback_;
4277    // These are only used as temporary arrays so need not be filled
4278    if (maximumNumberCuts_) {
4279      addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4280    } else {
4281      addedCuts_ = NULL;
4282    }
4283    nextRowCut_ = NULL;
4284    currentNode_ = NULL;
4285    if (maximumDepth_)
4286      walkback_ = new CbcNodeInfo * [maximumDepth_];
4287    else
4288      walkback_ = NULL;
4289    synchronizeModel();
4290    cbcColLower_ = NULL;
4291    cbcColUpper_ = NULL;
4292    cbcRowLower_ = NULL;
4293    cbcRowUpper_ = NULL;
4294    cbcColSolution_ = NULL;
4295    cbcRowPrice_ = NULL;
4296    cbcReducedCost_ = NULL;
4297    cbcRowActivity_ = NULL;
4298  }
4299  return *this;
4300}
4301 
4302// Destructor
4303CbcModel::~CbcModel ()
4304{
4305  if (defaultHandler_) {
4306    delete handler_;
4307    handler_ = NULL;
4308  }
4309  delete tree_;
4310  tree_=NULL;
4311  if (ourSolver_) delete solver_;
4312  gutsOfDestructor();
4313}
4314// Clears out as much as possible (except solver)
4315void 
4316CbcModel::gutsOfDestructor()
4317{
4318  delete referenceSolver_;
4319  referenceSolver_=NULL;
4320  int i;
4321  for (i=0;i<numberCutGenerators_;i++) {
4322    delete generator_[i];
4323    delete virginGenerator_[i];
4324  }
4325  delete [] generator_;
4326  delete [] virginGenerator_;
4327  generator_=NULL;
4328  virginGenerator_=NULL;
4329  for (i=0;i<numberHeuristics_;i++)
4330    delete heuristic_[i];
4331  delete [] heuristic_;
4332  heuristic_=NULL;
4333  delete nodeCompare_;
4334  nodeCompare_=NULL;
4335  delete problemFeasibility_;
4336  problemFeasibility_=NULL;
4337  delete [] originalColumns_;
4338  originalColumns_=NULL;
4339  delete strategy_;
4340  gutsOfDestructor2();
4341}
4342// Clears out enough to reset CbcModel
4343void 
4344CbcModel::gutsOfDestructor2()
4345{
4346  delete [] integerInfo_;
4347  integerInfo_=NULL;
4348  delete [] integerVariable_;
4349  integerVariable_=NULL;
4350  int i;
4351  for (i=0;i<numberObjects_;i++)
4352    delete object_[i];
4353  delete [] object_;
4354  object_=NULL;
4355  numberIntegers_=0;
4356  numberObjects_=0;
4357  delete emptyWarmStart_ ;
4358  emptyWarmStart_ =NULL;
4359  delete continuousSolver_;
4360  continuousSolver_=NULL;
4361  delete [] bestSolution_;
4362  bestSolution_=NULL;
4363  delete [] currentSolution_;
4364  currentSolution_=NULL;
4365  delete [] continuousSolution_;
4366  continuousSolution_=NULL;
4367  solverCharacteristics_=NULL;
4368  delete [] usedInSolution_;
4369  usedInSolution_ = NULL;
4370  testSolution_=NULL;
4371  lastHeuristic_ = NULL;
4372  delete [] addedCuts_;
4373  addedCuts_=NULL;
4374  nextRowCut_ = NULL;
4375  currentNode_ = NULL;
4376  delete [] walkback_;
4377  walkback_=NULL;
4378  delete [] whichGenerator_;
4379  whichGenerator_ = NULL;
4380  for (i=0;i<maximumStatistics_;i++)
4381    delete statistics_[i];
4382  delete [] statistics_;
4383  statistics_=NULL;
4384  maximumStatistics_=0;
4385  delete [] analyzeResults_;
4386  analyzeResults_=NULL;
4387  // Below here is whatever consensus is
4388  ourSolver_=true;
4389  bestObjective_=COIN_DBL_MAX;
4390  bestPossibleObjective_=COIN_DBL_MAX;
4391  sumChangeObjective1_=0.0;
4392  sumChangeObjective2_=0.0;
4393  numberSolutions_=0;
4394  stateOfSearch_=0;
4395  delete [] hotstartSolution_;
4396  hotstartSolution_=NULL;
4397  delete [] hotstartPriorities_;
4398  hotstartPriorities_=NULL;
4399  numberHeuristicSolutions_=0;
4400  numberNodes_=0;
4401  numberNodes2_=0;
4402  numberIterations_=0;
4403  status_=-1;
4404  secondaryStatus_=-1;
4405  maximumNumberCuts_=0;
4406  phase_=0;
4407  currentNumberCuts_=0;
4408  maximumDepth_=0;
4409  nextRowCut_=NULL;
4410  currentNode_=NULL;
4411  // clear out tree
4412  if (tree_&&tree_->size())
4413    tree_->cleanTree(this, -1.0e100,bestPossibleObjective_) ;
4414  subTreeModel_=NULL;
4415  numberStoppedSubTrees_=0;
4416  numberInfeasibleNodes_=0;
4417  numberGlobalViolations_=0;
4418  continuousObjective_=0.0;
4419  originalContinuousObjective_=0.0;
4420  continuousInfeasibilities_=0;
4421  numberFixedAtRoot_=0;
4422  numberFixedNow_=0;
4423  stoppedOnGap_=false;
4424  eventHappened_=false;
4425  numberLongStrong_=0;
4426  numberOldActiveCuts_=0;
4427  numberNewCuts_=0;
4428  searchStrategy_=-1;
4429  numberStrongIterations_=0;
4430  // Parameters which need to be reset
4431  dblParam_[CbcCutoffIncrement] = 1e-5;
4432  dblParam_[CbcCurrentCutoff] = 1.0e100;
4433  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4434  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4435}
4436// Save a copy of the current solver so can be reset to
4437void 
4438CbcModel::saveReferenceSolver()
4439{
4440  delete referenceSolver_;
4441  referenceSolver_= solver_->clone();
4442}
4443
4444// Uses a copy of reference solver to be current solver
4445void 
4446CbcModel::resetToReferenceSolver()
4447{
4448  delete solver_;
4449  solver_ = referenceSolver_->clone();
4450  // clear many things
4451  gutsOfDestructor2();
4452  // Reset cutoff
4453  // Solvers know about direction
4454  double direction = solver_->getObjSense();
4455  double value;
4456  solver_->getDblParam(OsiDualObjectiveLimit,value); 
4457  setCutoff(value*direction);
4458}
4459
4460// Are there a numerical difficulties?
4461bool 
4462CbcModel::isAbandoned() const
4463{
4464  return status_ == 2;
4465}
4466// Is optimality proven?
4467bool 
4468CbcModel::isProvenOptimal() const
4469{
4470  if (!status_ && bestObjective_<1.0e30)
4471    return true;
4472  else
4473    return false;
4474}
4475// Is  infeasiblity proven (or none better than cutoff)?
4476bool 
4477CbcModel::isProvenInfeasible() const
4478{
4479  if (!status_ && bestObjective_>=1.0e30)
4480    return true;
4481  else
4482    return false;
4483}
4484// Node limit reached?
4485bool 
4486CbcModel::isNodeLimitReached() const
4487{
4488  return numberNodes_ >= intParam_[CbcMaxNumNode];
4489}
4490// Time limit reached?
4491bool 
4492CbcModel::isSecondsLimitReached() const
4493{
4494  if (status_==1&&secondaryStatus_==4)
4495    return true;
4496  else
4497    return false;
4498}
4499// Solution limit reached?
4500bool 
4501CbcModel::isSolutionLimitReached() const
4502{
4503  return numberSolutions_ >= intParam_[CbcMaxNumSol];
4504}
4505// Set language
4506void 
4507CbcModel::newLanguage(CoinMessages::Language language)
4508{
4509  messages_ = CbcMessage(language);
4510}
4511void 
4512CbcModel::setNumberStrong(int number)
4513{
4514  if (number<0)
4515    numberStrong_=0;
4516   else
4517    numberStrong_=number;
4518}
4519void 
4520CbcModel::setNumberBeforeTrust(int number)
4521{
4522  if (number<-1) {
4523    numberBeforeTrust_=0;
4524  } else {
4525    numberBeforeTrust_=number;
4526    //numberStrong_ = CoinMax(numberStrong_,1);
4527  }
4528}
4529void 
4530CbcModel::setNumberPenalties(int number)
4531{
4532  if (number<=0) {
4533    numberPenalties_=0;
4534  } else {
4535    numberPenalties_=number;
4536  }
4537}
4538void 
4539CbcModel::setPenaltyScaleFactor(double value)
4540{
4541  if (value<=0) {
4542    penaltyScaleFactor_=3.0;
4543  } else {
4544    penaltyScaleFactor_=value;
4545  }
4546}
4547void 
4548CbcModel::setHowOftenGlobalScan(int number)
4549{
4550  if (number<-1)
4551    howOftenGlobalScan_=0;
4552   else
4553    howOftenGlobalScan_=number;
4554}
4555
4556// Add one generator
4557void 
4558CbcModel::addCutGenerator(CglCutGenerator * generator,
4559                          int howOften, const char * name,
4560                          bool normal, bool atSolution,
4561                          bool whenInfeasible,int howOftenInSub,
4562                          int whatDepth, int whatDepthInSub)
4563{
4564  CbcCutGenerator ** temp = generator_;
4565  generator_ = new CbcCutGenerator * [numberCutGenerators_+1];
4566  memcpy(generator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
4567  delete[] temp ;
4568  generator_[numberCutGenerators_]= 
4569    new CbcCutGenerator(this,generator, howOften, name,
4570                        normal,atSolution,whenInfeasible,howOftenInSub,
4571                        whatDepth, whatDepthInSub);
4572  // and before any cahnges
4573  temp = virginGenerator_;
4574  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_+1];
4575  memcpy(virginGenerator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
4576  delete[] temp ;
4577  virginGenerator_[numberCutGenerators_++]= 
4578    new CbcCutGenerator(this,generator, howOften, name,
4579                        normal,atSolution,whenInfeasible,howOftenInSub,
4580                        whatDepth, whatDepthInSub);
4581                                                         
4582}
4583// Add one heuristic
4584void 
4585CbcModel::addHeuristic(CbcHeuristic * generator)
4586{
4587  CbcHeuristic ** temp = heuristic_;
4588  heuristic_ = new CbcHeuristic * [numberHeuristics_+1];
4589  memcpy(heuristic_,temp,numberHeuristics_*sizeof(CbcHeuristic *));
4590  delete [] temp;
4591  heuristic_[numberHeuristics_++]=generator->clone();
4592}
4593
4594/*
4595  The last subproblem handled by the solver is not necessarily related to the
4596  one being recreated, so the first action is to remove all cuts from the
4597  constraint system.  Next, traverse the tree from node to the root to
4598  determine the basis size required for this subproblem and create an empty
4599  basis with the right capacity.  Finally, traverse the tree from root to
4600  node, adjusting bounds in the constraint system, adjusting the basis, and
4601  collecting the cuts that must be added to the constraint system.
4602  applyToModel does the heavy lifting.
4603
4604  addCuts1 is used in contexts where all that's desired is the list of cuts:
4605  the node is already fathomed, and we're collecting cuts so that we can
4606  adjust reference counts as we prune nodes. Arguably the two functions
4607  should be separated. The culprit is applyToModel, which performs cut
4608  collection and model adjustment.
4609
4610  Certainly in the contexts where all we need is a list of cuts, there's no
4611  point in passing in a valid basis --- an empty basis will do just fine.
4612*/
4613void CbcModel::addCuts1 (CbcNode * node, CoinWarmStartBasis *&lastws)
4614{ int i;
4615  int nNode=0;
4616  int numberColumns = getNumCols();
4617  CbcNodeInfo * nodeInfo = node->nodeInfo();
4618
4619/*
4620  Remove all cuts from the constraint system.
4621  (original comment includes ``see note below for later efficiency'', but
4622  the reference isn't clear to me).
4623*/
4624  int currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_;
4625  int *which = new int[currentNumberCuts];
4626  for (i = 0 ; i < currentNumberCuts ; i++)
4627    which[i] = i+numberRowsAtContinuous_;
4628  solver_->deleteRows(currentNumberCuts,which);
4629  delete [] which;
4630/*
4631  Accumulate the path from node to the root in walkback_, and accumulate a
4632  cut count in currentNumberCuts.
4633
4634  original comment: when working then just unwind until where new node joins
4635  old node (for cuts?)
4636*/
4637  currentNumberCuts = 0;
4638  while (nodeInfo) {
4639    //printf("nNode = %d, nodeInfo = %x\n",nNode,nodeInfo);
4640    walkback_[nNode++]=nodeInfo;
4641    currentNumberCuts += nodeInfo->numberCuts() ;
4642    nodeInfo = nodeInfo->parent() ;
4643    if (nNode==maximumDepth_) {
4644      maximumDepth_ *= 2;
4645      CbcNodeInfo ** temp = new CbcNodeInfo * [maximumDepth_];
4646      for (i=0;i<nNode;i++) 
4647        temp[i] = walkback_[i];
4648      delete [] walkback_;
4649      walkback_ = temp;
4650    }
4651  }
4652/*
4653  Create an empty basis with sufficient capacity for the constraint system
4654  we'll construct: original system plus cuts. Make sure we have capacity to
4655  record those cuts in addedCuts_.
4656
4657  The method of adjusting the basis at a FullNodeInfo object (the root, for
4658  example) is to use a copy constructor to duplicate the basis held in the
4659  nodeInfo, then resize it and return the new basis object. Guaranteed,
4660  lastws will point to a different basis when it returns. We pass in a basis
4661  because we need the parameter to return the allocated basis, and it's an
4662  easy way to pass in the size. But we take a hit for memory allocation.
4663*/
4664  currentNumberCuts_=currentNumberCuts;
4665  if (currentNumberCuts >= maximumNumberCuts_) {
4666    maximumNumberCuts_ = currentNumberCuts;
4667    delete [] addedCuts_;
4668    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4669  }
4670  lastws->setSize(numberColumns,numberRowsAtContinuous_+currentNumberCuts);
4671/*
4672  This last bit of code traverses the path collected in walkback_ from the
4673  root back to node. At the end of the loop,
4674   * lastws will be an appropriate basis for node;
4675   * variable bounds in the constraint system will be set to be correct for
4676     node; and
4677   * addedCuts_ will be set to a list of cuts that need to be added to the
4678     constraint system at node.
4679  applyToModel does all the heavy lifting.
4680*/
4681  currentNumberCuts=0;
4682  while (nNode) {
4683    --nNode;
4684    walkback_[nNode]->applyToModel(this,lastws,addedCuts_,currentNumberCuts);
4685  }
4686}
4687
4688/*
4689  adjustCuts might be a better name: If the node is feasible, we sift through
4690  the cuts we've collected, add the ones that are tight and omit the ones that
4691  are loose. If the node is infeasible, we just adjust the reference counts to
4692  reflect that we're about to prune this node and its descendants.
4693
4694  The reason we need to pass in lastws is that OsiClp automagically corrects
4695  the basis when it deletes constraints. So when all cuts are stripped within
4696  addCuts1, we lose their basis entries, hence the ability to determine if
4697  they are loose or tight. The question is whether we really need to pass in
4698  a basis or if we can capture it here. I'm thinking we can capture it here
4699  and pass it back out if required.
4700*/
4701int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws,bool canFix)
4702{
4703/*
4704  addCuts1 performs step 1 of restoring the subproblem at this node; see the
4705  comments there.
4706*/
4707  addCuts1(node,lastws);
4708  int i;
4709  int numberColumns = getNumCols();
4710  CbcNodeInfo * nodeInfo = node->nodeInfo();
4711  double cutoff = getCutoff() ;
4712  int currentNumberCuts=currentNumberCuts_;
4713  if (canFix) {
4714    bool feasible=true;
4715    const double *lower = solver_->getColLower() ;
4716    const double *upper = solver_->getColUpper() ;
4717    double * newLower = analyzeResults_;
4718    double * objLower = newLower+numberIntegers_;
4719    double * newUpper = objLower+numberIntegers_;
4720    double * objUpper = newUpper+numberIntegers_;
4721    int n=0;
4722    for (i=0;i<numberIntegers_;i++) {
4723      int iColumn = integerVariable_[i];
4724      bool changed=false;
4725      double lo = 0.0;
4726      double up = 0.0;
4727      if (objLower[i]>cutoff) {
4728        lo = lower[iColumn];
4729        up = upper[iColumn];
4730        if (lo<newLower[i]) {
4731          lo = newLower[i];
4732          solver_->setColLower(iColumn,lo);
4733          changed=true;
4734          n++;
4735        }
4736        if (objUpper[i]>cutoff) {
4737          if (up>newUpper[i]) {
4738            up = newUpper[i];
4739            solver_->setColUpper(iColumn,up);
4740            changed=true;
4741            n++;
4742          }
4743        }
4744      } else if (objUpper[i]>cutoff) {
4745        lo = lower[iColumn];
4746        up = upper[iColumn];
4747        if (up>newUpper[i]) {
4748          up = newUpper[i];
4749          solver_->setColUpper(iColumn,up);
4750          changed=true;
4751          n++;
4752        }
4753      }
4754      if (changed&&lo>up) {
4755        feasible=false;
4756        break;
4757      }
4758    }
4759    if (!feasible) {
4760      printf("analysis says node infeas\n");
4761      cutoff=-COIN_DBL_MAX;
4762    }
4763  }
4764/*
4765  If the node can't be fathomed by bound, reinstall tight cuts in the
4766  constraint system.
4767*/
4768  if (node->objectiveValue() < cutoff)
4769  { int numberToAdd = 0;
4770    const OsiRowCut * * addCuts;
4771    if (currentNumberCuts == 0)
4772      addCuts = NULL;
4773    else
4774      addCuts = new const OsiRowCut  * [currentNumberCuts];
4775#   ifdef CHECK_CUT_COUNTS
4776    printf("addCuts: expanded basis; rows %d+%d\n",
4777           numberRowsAtContinuous_,currentNumberCuts);
4778    lastws->print();
4779#   endif
4780/*
4781  Adjust the basis and constraint system so that we retain only active cuts.
4782  There are three steps:
4783    1) Scan the basis. If the logical associated with the cut is basic, it's
4784       loose and we drop it. The status of the logical for tight cuts is
4785       written back into the status array, compressing as we go.
4786    2) Resize the basis to fit the number of active cuts, stash a clone, and
4787       install with a call to setWarmStart().
4788    3) Install the tight cuts into the constraint system (applyRowCuts).
4789
4790    Update lastws
4791*/
4792    int nxrow = lastws->getNumArtificial();
4793    for (i=0;i<currentNumberCuts;i++) {
4794      assert (i+numberRowsAtContinuous_<nxrow);
4795      CoinWarmStartBasis::Status status = 
4796        lastws->getArtifStatus(i+numberRowsAtContinuous_);
4797      if (addedCuts_[i]&&(status != CoinWarmStartBasis::basic||addedCuts_[i]->effectiveness()==COIN_DBL_MAX)) {
4798#       ifdef CHECK_CUT_COUNTS
4799        printf("Using cut %d %x as row %d\n",i,addedCuts_[i],
4800               numberRowsAtContinuous_+numberToAdd);
4801#       endif
4802        lastws->setArtifStatus(numberToAdd+numberRowsAtContinuous_,status);
4803        addCuts[numberToAdd++] = new OsiRowCut(*addedCuts_[i]);
4804      } else {
4805#       ifdef CHECK_CUT_COUNTS
4806        printf("Dropping cut %d %x\n",i,addedCuts_[i]);
4807#       endif
4808        addedCuts_[i]=NULL;
4809      }
4810    }
4811    int numberRowsNow=numberRowsAtContinuous_+numberToAdd;
4812    lastws->resize(numberRowsNow,numberColumns);
4813#ifdef FULL_DEBUG
4814    printf("addCuts: stripped basis; rows %d + %d\n",
4815           numberRowsAtContinuous_,numberToAdd);
4816    lastws->print();
4817#endif
4818/*
4819  Apply the cuts and set the basis in the solver.
4820*/
4821    solver_->applyRowCuts(numberToAdd,addCuts);
4822    solver_->setWarmStart(lastws);
4823
4824#if 0
4825    if ((numberNodes_%printFrequency_)==0) {
4826      printf("Objective %g, depth %d, unsatisfied %d\n",
4827             node->objectiveValue(),
4828             node->depth(),node->numberUnsatisfied());
4829    }
4830#endif
4831/*
4832  Clean up and we're out of here.
4833*/
4834    for (i=0;i<numberToAdd;i++)
4835      delete addCuts[i];
4836    delete [] addCuts;
4837    numberNodes_++;
4838    return 0;
4839  } 
4840/*
4841  This node has been fathomed by bound as we try to revive it out of the live
4842  set. Adjust the cut reference counts to reflect that we no longer need to
4843  explore the remaining branch arms, hence they will no longer reference any
4844  cuts. Cuts whose reference count falls to zero are deleted.
4845*/
4846  else
4847  { int i;
4848    int numberLeft = nodeInfo->numberBranchesLeft();
4849    for (i = 0 ; i < currentNumberCuts ; i++)
4850    { if (addedCuts_[i])
4851      { if (!addedCuts_[i]->decrement(numberLeft))
4852        { delete addedCuts_[i];
4853          addedCuts_[i] = NULL; } } }
4854    return 1 ; }
4855}
4856
4857
4858/*
4859  Perform reduced cost fixing on integer variables.
4860
4861  The variables in question are already nonbasic at bound. We're just nailing
4862  down the current situation.
4863*/
4864
4865int CbcModel::reducedCostFix ()
4866
4867{
4868  if(!solverCharacteristics_->reducedCostsAccurate())
4869    return 0; //NLP
4870  double cutoff = getCutoff() ;
4871  double direction = solver_->getObjSense() ;
4872  double gap = cutoff - solver_->getObjValue()*direction ;
4873  double tolerance;
4874  solver_->getDblParam(OsiDualTolerance,tolerance) ;
4875  if (gap<=0.0)
4876    return 0;
4877  gap += 100.0*tolerance;
4878  double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4879
4880  const double *lower = solver_->getColLower() ;
4881  const double *upper = solver_->getColUpper() ;
4882  const double *solution = solver_->getColSolution() ;
4883  const double *reducedCost = solver_->getReducedCost() ;
4884
4885  int numberFixed = 0 ;
4886#ifdef COIN_USE_CLP
4887  OsiClpSolverInterface * clpSolver
4888    = dynamic_cast<OsiClpSolverInterface *> (solver_);
4889  ClpSimplex * clpSimplex=NULL;
4890  if (clpSolver) 
4891    clpSimplex = clpSolver->getModelPtr();
4892#endif
4893  for (int i = 0 ; i < numberIntegers_ ; i++)
4894  { int iColumn = integerVariable_[i] ;
4895    double djValue = direction*reducedCost[iColumn] ;
4896    if (upper[iColumn]-lower[iColumn] > integerTolerance)
4897    { if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap)
4898      { solver_->setColUpper(iColumn,lower[iColumn]) ;
4899#ifdef COIN_USE_CLP
4900      if (clpSimplex)
4901        assert (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound);
4902#endif
4903        numberFixed++ ; }
4904      else
4905      if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap)
4906      { solver_->setColLower(iColumn,upper[iColumn]) ;
4907#ifdef COIN_USE_CLP
4908      if (clpSimplex)
4909        assert (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound);
4910#endif
4911        numberFixed++ ; } } }
4912 
4913  return numberFixed; }
4914
4915// Collect coding to replace whichGenerator
4916void
4917CbcModel::resizeWhichGenerator(int numberNow, int numberAfter)
4918{
4919  if (numberAfter > maximumWhich_) {
4920    maximumWhich_ = CoinMax(maximumWhich_*2+100,numberAfter) ;
4921    int * temp = new int[2*maximumWhich_] ;
4922    memcpy(temp,whichGenerator_,numberNow*sizeof(int)) ;
4923    delete [] whichGenerator_ ;
4924    whichGenerator_ = temp ;
4925    memset(whichGenerator_+numberNow,0,(maximumWhich_-numberNow)*sizeof(int));
4926  }
4927}
4928
4929/** Solve the model using cuts
4930
4931  This version takes off redundant cuts from node.
4932  Returns true if feasible.
4933
4934  \todo
4935  Why do I need to resolve the problem? What has been done between the last
4936  relaxation and calling solveWithCuts?
4937
4938  If numberTries == 0 then user did not want any cuts.
4939*/
4940
4941bool 
4942CbcModel::solveWithCuts (OsiCuts &cuts, int numberTries, CbcNode *node)
4943/*
4944  Parameters:
4945    numberTries: (i) the maximum number of iterations for this round of cut
4946                     generation; if negative then we don't mind if drop is tiny.
4947   
4948    cuts:       (o) all cuts generated in this round of cut generation
4949
4950    node: (i)     So we can update dynamic pseudo costs
4951*/
4952                       
4953
4954{
4955  bool feasible = true ;
4956  int lastNumberCuts = 0 ;
4957  double lastObjective = -1.0e100 ;
4958  int violated = 0 ;
4959  int numberRowsAtStart = solver_->getNumRows() ;
4960  //printf("solver had %d rows\n",numberRowsAtStart);
4961  int numberColumns = solver_->getNumCols() ;
4962  CoinBigIndex numberElementsAtStart = solver_->getNumElements();
4963
4964  numberOldActiveCuts_ = numberRowsAtStart-numberRowsAtContinuous_ ;
4965  numberNewCuts_ = 0 ;
4966
4967  bool onOptimalPath = false ;
4968  const OsiRowCutDebugger *debugger = NULL;
4969  if ((specialOptions_&1)!=0) {
4970    /*
4971      See OsiRowCutDebugger for details. In a nutshell, make sure that current
4972      variable values do not conflict with a known optimal solution. (Obviously
4973      this can be fooled when there are multiple solutions.)
4974    */
4975    debugger = solver_->getRowCutDebugger() ;
4976    if (debugger) 
4977      onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
4978  }
4979  OsiCuts slackCuts;
4980/*
4981  Resolve the problem. If we've lost feasibility, might as well bail out right
4982  after the debug stuff.
4983*/
4984  double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
4985  if (node)
4986    objectiveValue= node->objectiveValue();
4987  feasible = resolve() ;
4988  if (problemFeasibility_->feasible(this,0)<0) {
4989    feasible=false; // pretend infeasible
4990  }
4991 
4992  // Update branching information if wanted
4993  if(node &&branchingMethod_)
4994    branchingMethod_->updateInformation(solver_,node);
4995
4996#ifdef CBC_DEBUG
4997  if (feasible)
4998  { printf("Obj value %g (%s) %d rows\n",solver_->getObjValue(),
4999           (solver_->isProvenOptimal())?"proven":"unproven",
5000           solver_->getNumRows()) ; }
5001 
5002  else
5003  { printf("Infeasible %d rows\n",solver_->getNumRows()) ; }
5004#endif
5005  if ((specialOptions_&1)!=0) {
5006/*
5007  If the RowCutDebugger said we were compatible with the optimal solution,
5008  and now we're suddenly infeasible, we might be confused. Then again, we
5009  may have fathomed by bound, heading for a rediscovery of an optimal solution.
5010*/
5011    if (onOptimalPath && !solver_->isDualObjectiveLimitReached()) {
5012      if (!feasible) {
5013        solver_->writeMps("infeas");
5014        CoinWarmStartBasis *slack =
5015          dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
5016        solver_->setWarmStart(slack);
5017        delete slack ;
5018        solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
5019        solver_->initialSolve();
5020      }
5021      assert(feasible) ;
5022    }
5023  }
5024
5025  if (!feasible) {
5026    numberInfeasibleNodes_++;
5027    return (false) ;
5028  }
5029  sumChangeObjective1_ += solver_->getObjValue()*solver_->getObjSense()
5030    - objectiveValue ;
5031  if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
5032    numberTries=0; // exit
5033  //if ((numberNodes_%100)==0)
5034  //printf("XXa sum obj changed by %g\n",sumChangeObjective1_);
5035  objectiveValue = solver_->getObjValue()*solver_->getObjSense();
5036  // Return at once if numberTries zero
5037  if (!numberTries) {
5038    cuts=OsiCuts();
5039    numberNewCuts_=0;
5040    return true;
5041  }
5042/*
5043  Do reduced cost fixing, and then grab the primal solution and bounds vectors.
5044*/
5045  reducedCostFix() ;
5046  const double *lower = solver_->getColLower() ;
5047  const double *upper = solver_->getColUpper() ;
5048  const double *solution = solver_->getColSolution() ;
5049/*
5050  Set up for at most numberTries rounds of cut generation. If numberTries is
5051  negative, we'll ignore the minimumDrop_ cutoff and keep generating cuts for
5052  the specified number of rounds.
5053*/
5054  double minimumDrop = minimumDrop_ ;
5055  if (numberTries<0)
5056  { numberTries = -numberTries ;
5057    minimumDrop = -1.0 ; }
5058/*
5059  Is it time to scan the cuts in order to remove redundant cuts? If so, set
5060  up to do it.
5061*/
5062# define SCANCUTS 100 
5063  int *countColumnCuts = NULL ;
5064  // Always accumulate row cut counts
5065  int * countRowCuts =new int[numberCutGenerators_+numberHeuristics_] ;
5066  memset(countRowCuts,0,
5067         (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ;
5068  bool fullScan = false ;
5069  if ((numberNodes_%SCANCUTS) == 0)
5070  { fullScan = true ;
5071    countColumnCuts = new int[numberCutGenerators_+numberHeuristics_] ;
5072    memset(countColumnCuts,0,
5073           (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ; }
5074
5075  double direction = solver_->getObjSense() ;
5076  double startObjective = solver_->getObjValue()*direction ;
5077
5078  currentPassNumber_ = 0 ;
5079  double primalTolerance = 1.0e-7 ;
5080/*
5081  Begin cut generation loop. Cuts generated during each iteration are
5082  collected in theseCuts. The loop can be divided into four phases:
5083   1) Prep: Fix variables using reduced cost. In the first iteration only,
5084      consider scanning globalCuts_ and activating any applicable cuts.
5085   2) Cut Generation: Call each generator and heuristic registered in the
5086      generator_ and heuristic_ arrays. Newly generated global cuts are
5087      copied to globalCuts_ at this time.
5088   3) Cut Installation and Reoptimisation: Install column and row cuts in
5089      the solver. Copy row cuts to cuts (parameter). Reoptimise.
5090   4) Cut Purging: takeOffCuts() removes inactive cuts from the solver, and
5091      does the necessary bookkeeping in the model.
5092*/
5093  do
5094  { currentPassNumber_++ ;
5095    numberTries-- ;
5096    OsiCuts theseCuts ;
5097/*
5098  Scan previously generated global column and row cuts to see if any are
5099  useful.
5100  I can't see why this code
5101  needs its own copy of the primal solution. Removed the dec'l.
5102*/
5103    int numberViolated=0;
5104    if (currentPassNumber_ == 1 && howOftenGlobalScan_ > 0 &&
5105        (numberNodes_%howOftenGlobalScan_) == 0)
5106    { int numberCuts = globalCuts_.sizeColCuts() ;
5107      int i;
5108      // possibly extend whichGenerator
5109      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
5110      for ( i = 0 ; i < numberCuts ; i++)
5111      { const OsiColCut *thisCut = globalCuts_.colCutPtr(i) ;
5112        if (thisCut->violated(solution)>primalTolerance) {
5113          printf("Global cut added - violation %g\n",
5114                 thisCut->violated(solution)) ;
5115          whichGenerator_[numberViolated++]=-1;
5116          theseCuts.insert(*thisCut) ;
5117        }
5118      }
5119      numberCuts = globalCuts_.sizeRowCuts() ;
5120      // possibly extend whichGenerator
5121      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
5122      for ( i = 0;i<numberCuts;i++) {
5123        const OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ;
5124        if (thisCut->violated(solution)>primalTolerance) {
5125          //printf("Global cut added - violation %g\n",
5126          // thisCut->violated(solution)) ;
5127          whichGenerator_[numberViolated++]=-1;
5128          theseCuts.insert(*thisCut) ;
5129        }
5130      }
5131      numberGlobalViolations_+=numberViolated;
5132    }
5133/*
5134  Generate new cuts (global and/or local) and/or apply heuristics.  If
5135  CglProbing is used, then it should be first as it can fix continuous
5136  variables.
5137
5138  At present, CglProbing is the only case where generateCuts will return
5139  true. generateCuts actually modifies variable bounds in the solver when
5140  CglProbing indicates that it can fix a variable. Reoptimisation is required
5141  to take full advantage.
5142*/
5143    // This should only happen after heuristic solution
5144    if (solverCharacteristics_->warmStart()&&
5145        !solver_->optimalBasisIsAvailable()) {
5146      //printf("XXXXYY no opt basis\n");
5147      resolve();
5148    }
5149    if (nextRowCut_) {
5150      // branch was a cut - add it
5151      theseCuts.insert(*nextRowCut_);
5152      if (handler_->logLevel()>1)
5153        nextRowCut_->print();
5154      const OsiRowCut * cut=nextRowCut_;
5155      const double * solution = solver_->getColSolution();
5156      double lb = cut->lb();
5157      double ub = cut->ub();
5158      int n=cut->row().getNumElements();
5159      const int * column = cut->row().getIndices();
5160      const double * element = cut->row().getElements();
5161      double sum=0.0;
5162      for (int i=0;i<n;i++) {
5163        int iColumn = column[i];
5164        double value = element[i];
5165        //if (solution[iColumn]>1.0e-7)
5166        //printf("value of %d is %g\n",iColumn,solution[iColumn]);
5167        sum += value * solution[iColumn];
5168      }
5169      delete nextRowCut_;
5170      nextRowCut_=NULL;
5171      if (handler_->logLevel()>1)
5172        printf("applying branch cut, sum is %g, bounds %g %g\n",sum,lb,ub);
5173      // possibly extend whichGenerator
5174      resizeWhichGenerator(numberViolated, numberViolated+1);
5175      // set whichgenerator (also serves as marker to say don't delete0
5176      whichGenerator_[numberViolated++]=-2;
5177    }
5178    double * newSolution = new double [numberColumns] ;
5179    double heuristicValue = getCutoff() ;
5180    int found = -1; // no solution found
5181
5182    for (int i = 0;i<numberCutGenerators_+numberHeuristics_;i++) {
5183      int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
5184      int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
5185      if (i<numberCutGenerators_) {
5186        bool generate = generator_[i]->normal();
5187        // skip if not optimal and should be (maybe a cut generator has fixed variables)
5188        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
5189          generate=false;
5190        if (generate) {
5191          bool mustResolve = 
5192            generator_[i]->generateCuts(theseCuts,fullScan,node) ;
5193#ifdef CBC_DEBUG
5194          {
5195            int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
5196            int k ;
5197            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
5198              OsiRowCut thisCut = theseCuts.rowCut(k) ;
5199              /* check size of elements.
5200                 We can allow smaller but this helps debug generators as it
5201                 is unsafe to have small elements */
5202              int n=thisCut.row().getNumElements();
5203              const int * column = thisCut.row().getIndices();
5204              const double * element = thisCut.row().getElements();
5205              //assert (n);
5206              for (int i=0;i<n;i++) {
5207                double value = element[i];
5208                assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
5209              }
5210            }
5211          }
5212#endif
5213          if (mustResolve) {
5214            feasible = resolve() ;
5215            if ((specialOptions_&1)!=0) {
5216              debugger = solver_->getRowCutDebugger() ;
5217              if (debugger) 
5218                onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
5219              else
5220                onOptimalPath=false;
5221              if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
5222                assert(feasible) ;
5223            }
5224            if (!feasible)
5225              break ;
5226          }
5227        }
5228      } else {
5229        // see if heuristic will do anything
5230        double saveValue = heuristicValue ;
5231        int ifSol = 
5232          heuristic_[i-numberCutGenerators_]->solution(heuristicValue,
5233                                                       newSolution,
5234                                                       theseCuts) ;
5235        if (ifSol>0) {
5236          // better solution found
5237          found = i-numberCutGenerators_ ;
5238          incrementUsed(newSolution);
5239        } else if (ifSol<0) {
5240          heuristicValue = saveValue ;
5241        }
5242      }
5243      int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
5244      int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
5245
5246      if ((specialOptions_&1)!=0) {
5247        if (onOptimalPath) {
5248          int k ;
5249          for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
5250            OsiRowCut thisCut = theseCuts.rowCut(k) ;
5251            if(debugger->invalidCut(thisCut)) {
5252              solver_->writeMps("badCut");
5253            }
5254            assert(!debugger->invalidCut(thisCut)) ;
5255          }
5256        }
5257      }
5258
5259/*
5260  The cut generator/heuristic has done its thing, and maybe it generated some
5261  cuts and/or a new solution.  Do a bit of bookkeeping: load
5262  whichGenerator[i] with the index of the generator responsible for a cut,
5263  and place cuts flagged as global in the global cut pool for the model.
5264
5265  lastNumberCuts is the sum of cuts added in previous iterations; it's the
5266  offset to the proper starting position in whichGenerator.
5267*/
5268      int numberBefore =
5269            numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
5270      int numberAfter =
5271            numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
5272      // possibly extend whichGenerator
5273      resizeWhichGenerator(numberBefore, numberAfter);
5274      int j ;
5275      if (fullScan) {
5276        // counts
5277        countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
5278      }
5279      countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
5280     
5281      for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
5282        whichGenerator_[numberBefore++] = i ;
5283        const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
5284        if (thisCut->lb()>thisCut->ub())
5285          violated=-2; // sub-problem is infeasible
5286        if (thisCut->globallyValid()) {
5287          // add to global list
5288          globalCuts_.insert(*thisCut) ;
5289        }
5290      }
5291      for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
5292        whichGenerator_[numberBefore++] = i ;
5293        const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
5294        if (thisCut->globallyValid()) {
5295          // add to global list
5296          globalCuts_.insert(*thisCut) ;
5297        }
5298      }
5299    }
5300    // If at root node and first pass do heuristics without cuts
5301    if (!numberNodes_&&currentPassNumber_==1) {
5302      // Save number solutions
5303      int saveNumberSolutions = numberSolutions_;
5304      int saveNumberHeuristicSolutions = numberHeuristicSolutions_;
5305      for (int i = 0;i<numberHeuristics_;i++) {
5306        // see if heuristic will do anything
5307        double saveValue = heuristicValue ;
5308        int ifSol = heuristic_[i]->solution(heuristicValue,
5309                                            newSolution);
5310        if (ifSol>0) {
5311          // better solution found
5312          found = i ;
5313          incrementUsed(newSolution);
5314          // increment number of solutions so other heuristics can test
5315          numberSolutions_++;
5316          numberHeuristicSolutions_++;
5317        } else {
5318          heuristicValue = saveValue ;
5319        }
5320      }
5321      // Restore number solutions
5322      numberSolutions_ = saveNumberSolutions;
5323      numberHeuristicSolutions_ = saveNumberHeuristicSolutions;
5324    }
5325    // Add in any violated saved cuts
5326    if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
5327      int numberOld = theseCuts.sizeRowCuts();
5328      int numberCuts = slackCuts.sizeRowCuts() ;
5329      int i;
5330      // possibly extend whichGenerator
5331      resizeWhichGenerator(numberOld, numberOld+numberCuts);
5332      for ( i = 0;i<numberCuts;i++) {
5333        const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
5334        if (thisCut->violated(solution)>100.0*primalTolerance) {
5335          if (messageHandler()->logLevel()>2)
5336            printf("Old cut added - violation %g\n",
5337                   thisCut->violated(solution)) ;
5338          whichGenerator_[numberOld++]=-1;
5339          theseCuts.insert(*thisCut) ;
5340        }
5341      }
5342    }
5343/*
5344  End of the loop to exercise each generator/heuristic.
5345
5346  Did any of the heuristics turn up a new solution? Record it before we free
5347  the vector.
5348*/
5349    if (found >= 0) { 
5350      phase_=4;
5351      incrementUsed(newSolution);
5352      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
5353      lastHeuristic_ = heuristic_[found];
5354      CbcTreeLocal * tree
5355          = dynamic_cast<CbcTreeLocal *> (tree_);
5356      if (tree)
5357        tree->passInSolution(bestSolution_,heuristicValue);
5358    }
5359    delete [] newSolution ;
5360
5361#if 0
5362    // switch on to get all cuts printed
5363    theseCuts.printCuts() ;
5364#endif
5365    int numberColumnCuts = theseCuts.sizeColCuts() ;
5366    int numberRowCuts = theseCuts.sizeRowCuts() ;
5367    if (violated>=0)
5368      violated = numberRowCuts + numberColumnCuts ;
5369/*
5370  Apply column cuts (aka bound tightening). This may be partially redundant
5371  for column cuts returned by CglProbing, as generateCuts installs bounds
5372  from CglProbing when it determines it can fix a variable.
5373
5374  TODO: Looks like the use of violated has evolved. The value set above is
5375        completely ignored. All that's left is violated == -1 indicates some
5376        cut is violated, violated == -2 indicates infeasibility. Only
5377        infeasibility warrants exceptional action.
5378
5379  TODO: Strikes me that this code will fail to detect infeasibility, because
5380        the breaks escape the inner loops but the outer loop keeps going.
5381        Infeasibility in an early cut will be overwritten if a later cut is
5382        merely violated.
5383*/
5384    if (numberColumnCuts) {
5385
5386#ifdef CBC_DEBUG
5387      double * oldLower = new double [numberColumns] ;
5388      double * oldUpper = new double [numberColumns] ;
5389      memcpy(oldLower,lower,numberColumns*sizeof(double)) ;
5390      memcpy(oldUpper,upper,numberColumns*sizeof(double)) ;
5391#endif
5392
5393      double integerTolerance = getDblParam(CbcIntegerTolerance) ;
5394      for (int i = 0;i<numberColumnCuts;i++) {
5395        const OsiColCut * thisCut = theseCuts.colCutPtr(i) ;
5396        const CoinPackedVector & lbs = thisCut->lbs() ;
5397        const CoinPackedVector & ubs = thisCut->ubs() ;
5398        int j ;
5399        int n ;
5400        const int * which ;
5401        const double * values ;
5402        n = lbs.getNumElements() ;
5403        which = lbs.getIndices() ;
5404        values = lbs.getElements() ;
5405        for (j = 0;j<n;j++) {
5406          int iColumn = which[j] ;
5407          double value = solution[iColumn] ;
5408#if CBC_DEBUG>1
5409          printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
5410                 solution[iColumn],oldUpper[iColumn],values[j]) ;
5411#endif
5412          solver_->setColLower(iColumn,values[j]) ;
5413          if (value<values[j]-integerTolerance)
5414            violated = -1 ;
5415          if (values[j]>upper[iColumn]+integerTolerance) {
5416            // infeasible
5417            violated = -2 ;
5418            break ;
5419          }
5420        }
5421        n = ubs.getNumElements() ;
5422        which = ubs.getIndices() ;
5423        values = ubs.getElements() ;
5424        for (j = 0;j<n;j++) {
5425          int iColumn = which[j] ;
5426          double value = solution[iColumn] ;
5427#if CBC_DEBUG>1
5428          printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
5429                 solution[iColumn],oldUpper[iColumn],values[j]) ;
5430#endif
5431          solver_->setColUpper(iColumn,values[j]) ;
5432          if (value>values[j]+integerTolerance)
5433            violated = -1 ;
5434          if (values[j]<lower[iColumn]-integerTolerance) {
5435            // infeasible
5436            violated = -2 ;
5437            break ;
5438          }
5439        }
5440      }
5441#ifdef CBC_DEBUG
5442      delete [] oldLower ;
5443      delete [] oldUpper ;
5444#endif
5445    }
5446/*
5447  End installation of column cuts. The break here escapes the numberTries
5448  loop.
5449*/
5450    if (violated == -2||!feasible) {
5451      // infeasible
5452      feasible = false ;
5453      violated = -2;
5454      if (!numberNodes_) 
5455        messageHandler()->message(CBC_INFEAS,
5456                                  messages())
5457                                    << CoinMessageEol ;
5458      break ;
5459    }
5460/*
5461  Now apply the row (constraint) cuts. This is a bit more work because we need
5462  to obtain and augment the current basis.
5463
5464  TODO: Why do this work, if there are no row cuts? The current basis will do
5465        just fine.
5466*/
5467    int numberRowsNow = solver_->getNumRows() ;
5468    assert(numberRowsNow == numberRowsAtStart+lastNumberCuts) ;
5469    int numberToAdd = theseCuts.sizeRowCuts() ;
5470    numberNewCuts_ = lastNumberCuts+numberToAdd ;
5471/*
5472  Now actually add the row cuts and reoptimise.
5473
5474  Install the cuts in the solver using applyRowCuts and
5475  augment the basis with the corresponding slack. We also add each row cut to
5476  the set of row cuts (cuts.insert()) supplied as a parameter. The new basis
5477  must be set with setWarmStart().
5478
5479  TODO: It's not clear to me why we can't separate this into two sections.
5480        The first would add the row cuts, and be executed only if row cuts
5481        need to be installed. The second would call resolve() and would be
5482        executed if either row or column cuts have been installed.
5483
5484  TODO: Seems to me the original code could allocate addCuts with size 0, if
5485        numberRowCuts was 0 and numberColumnCuts was nonzero. That might
5486        explain the memory fault noted in the comment by AJK.  Unfortunately,
5487        just commenting out the delete[] results in massive memory leaks. Try
5488        a revision to separate the row cut case. Why do we need addCuts at
5489        all? A typing issue, apparently: OsiCut vs. OsiRowCut.
5490
5491  TODO: It looks to me as if numberToAdd and numberRowCuts are identical at
5492        this point. Confirm & get rid of one of them.
5493
5494  TODO: Any reason why the three loops can't be consolidated?
5495*/
5496    if (numberRowCuts > 0 || numberColumnCuts > 0)
5497    { if (numberToAdd > 0)
5498      { int i ;
5499        // Faster to add all at once
5500        const OsiRowCut ** addCuts = new const OsiRowCut * [numberToAdd] ;
5501        for (i = 0 ; i < numberToAdd ; i++)
5502        { addCuts[i] = &theseCuts.rowCut(i) ; }
5503        solver_->applyRowCuts(numberToAdd,addCuts) ;
5504        // AJK this caused a memory fault on Win32
5505        // may be okay now with ** form
5506        delete [] addCuts ;
5507        for (i = 0 ; i < numberToAdd ; i++)
5508        { cuts.insert(theseCuts.rowCut(i)) ; }
5509        CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5510        assert(basis != NULL); // make sure not volume
5511        //basis->resize(numberRowsAtStart+numberNewCuts_,numberColumns) ;
5512        for (i = 0 ; i < numberToAdd ; i++) 
5513        { basis->setArtifStatus(numberRowsNow+i,
5514                                 CoinWarmStartBasis::basic) ; }
5515        if (solver_->setWarmStart(basis) == false)
5516        { throw CoinError("Fail setWarmStart() after cut installation.",
5517                          "solveWithCuts","CbcModel") ; }
5518        delete basis;
5519      }
5520      feasible = resolve() ;
5521      if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
5522        numberTries=0; // exit
5523#     ifdef CBC_DEBUG
5524      printf("Obj value after cuts %g %d rows\n",solver_->getObjValue(),
5525              solver_->getNumRows()) ;
5526      if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
5527        assert(feasible) ;
5528#     endif
5529    }
5530/*
5531  No cuts. Cut short the cut generation (numberTries) loop.
5532*/
5533    else
5534    { numberTries = 0 ; }
5535/*
5536  If the problem is still feasible, first, call takeOffCuts() to remove cuts
5537  that are now slack. takeOffCuts() will call the solver to reoptimise if
5538  that's needed to restore a valid solution.
5539
5540  Next, see if we should quit due to diminishing returns:
5541    * we've tried three rounds of cut generation and we're getting
5542      insufficient improvement in the objective; or
5543    * we generated no cuts; or
5544    * the solver declared optimality with 0 iterations after we added the
5545      cuts generated in this round.
5546  If we decide to keep going, prep for the next iteration.
5547
5548  It just seems more safe to tell takeOffCuts() to call resolve(), even if
5549  we're not continuing cut generation. Otherwise code executed between here
5550  and final disposition of the node will need to be careful not to access the
5551  lp solution. It can happen that we lose feasibility in takeOffCuts ---
5552  numerical jitters when the cutoff bound is epsilon less than the current
5553  best, and we're evaluating an alternative optimum.
5554
5555  TODO: After successive rounds of code motion, there seems no need to
5556        distinguish between the two checks for aborting the cut generation
5557        loop. Confirm and clean up.
5558*/
5559    if (feasible)
5560    { int cutIterations = solver_->getIterationCount() ;
5561      if (numberOldActiveCuts_+numberNewCuts_) {
5562        OsiCuts * saveCuts = node ? NULL : &slackCuts;
5563        takeOffCuts(cuts,resolveAfterTakeOffCuts_,saveCuts) ;
5564        if (solver_->isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_)
5565          { feasible = false ;
5566#       ifdef CBC_DEBUG
5567          double z = solver_->getObjValue() ;
5568          double cut = getCutoff() ;
5569          printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n",
5570                 z-cut,z,cut) ;
5571#       endif
5572          }
5573      }
5574      if (feasible)
5575        { numberRowsAtStart = numberOldActiveCuts_+numberRowsAtContinuous_ ;
5576        lastNumberCuts = numberNewCuts_ ;
5577        if (direction*solver_->getObjValue() < lastObjective+minimumDrop &&
5578            currentPassNumber_ >= 3)
5579          { numberTries = 0 ; }
5580        if (numberRowCuts+numberColumnCuts == 0 || cutIterations == 0)
5581          { break ; }
5582        if (numberTries > 0)
5583          { reducedCostFix() ;
5584          lastObjective = direction*solver_->getObjValue() ;
5585          lower = solver_->getColLower() ;
5586          upper = solver_->getColUpper() ;
5587          solution = solver_->getColSolution() ; } } }
5588/*
5589  We've lost feasibility --- this node won't be referencing the cuts we've
5590  been collecting, so decrement the reference counts.
5591
5592  TODO: Presumably this is in preparation for backtracking. Seems like it
5593        should be the `else' off the previous `if'.
5594*/
5595    if (!feasible)
5596    { int i ;
5597      for (i = 0;i<currentNumberCuts_;i++) {
5598        // take off node
5599        if (addedCuts_[i]) {
5600          if (!addedCuts_[i]->decrement())
5601            delete addedCuts_[i] ;
5602          addedCuts_[i] = NULL ;
5603        }
5604      }
5605      numberTries = 0 ;
5606    }
5607  } while (numberTries>0) ;
5608  //check feasibility.
5609  //If solution seems to be integer feasible calling setBestSolution
5610  //will eventually add extra global cuts which we need to install at
5611  //the nodes
5612
5613
5614  if(feasible&&solverCharacteristics_->solutionAddsCuts())  //check integer feasibility
5615  {
5616    bool integerFeasible = true;
5617    const double * save = testSolution_;
5618    testSolution_ = solver_->getColSolution();
5619    for (int i=0;i<numberObjects_ && integerFeasible;i++)
5620    {
5621      int preferredWay;
5622      double infeasibility = object_[i]->infeasibility(preferredWay);
5623      if(infeasibility)
5624        integerFeasible = false;
5625    }
5626    testSolution_ = save;
5627    if(integerFeasible)//update
5628    {
5629      double objValue = solver_->getObjValue();
5630      //int numberGlobalBefore = globalCuts_.sizeRowCuts();
5631      // SOLUTION2 so won't up cutoff or print message
5632      setBestSolution(CBC_SOLUTION2, objValue, 
5633                      solver_->getColSolution(),0);
5634#if 0
5635      int numberGlobalAfter = globalCuts_.sizeRowCuts();
5636      int numberToAdd = numberGlobalAfter - numberGlobalBefore;
5637      if(numberToAdd > 0)
5638        //We have added some cuts say they are tight at that node
5639        //Basis and lp should already have been updated
5640      {
5641        feasible = (solver_->isProvenOptimal() &&
5642                    !solver_->isDualObjectiveLimitReached()) ;
5643        if(feasible)
5644        {
5645          int numberCuts = numberNewCuts_ = cuts.sizeRowCuts();
5646      // possibly extend whichGenerator
5647          resizeWhichGenerator(numberCuts, numberToAdd+numberCuts);
5648         
5649          for (int i = numberGlobalBefore ; i < numberGlobalAfter ; i++)
5650            {       
5651              whichGenerator_[numberNewCuts_++]=-1;
5652              cuts.insert(globalCuts_.rowCut(i)) ;
5653            }
5654          numberNewCuts_ = lastNumberCuts+numberToAdd;
5655          //now take off the cuts which are not tight anymore
5656          takeOffCuts(cuts,resolveAfterTakeOffCuts_, NULL) ;
5657          if (solver_->isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_)
5658            { feasible = false ;}
5659        }
5660        if(!feasible) //node will be fathomed
5661          {
5662          for (int i = 0;i<currentNumberCuts_;i++)
5663            {
5664            // take off node
5665            if (addedCuts_[i])
5666              {
5667                if (!addedCuts_[i]->decrement())
5668                  delete addedCuts_[i] ;
5669                addedCuts_[i] = NULL ;
5670              }
5671            }
5672        }
5673      }
5674#endif
5675    }
5676  }
5677  // Reduced cost fix at end
5678  if (solver_->isProvenOptimal())
5679    reducedCostFix();
5680  // If at root node do heuristics
5681  if (!numberNodes_) {
5682    // First see if any cuts are slack
5683    int numberRows = solver_->getNumRows();
5684    int numberAdded = numberRows-numberRowsAtContinuous_;
5685    if (numberAdded) {
5686      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5687      assert(basis != NULL);
5688      int * added = new int[numberAdded];
5689      int nDelete = 0;
5690      for (int j=numberRowsAtContinuous_;j<numberRows;j++) {
5691        if (basis->getArtifStatus(j)==CoinWarmStartBasis::basic) {
5692          //printf("%d slack!\n",j);
5693          added[nDelete++]=j;
5694        }
5695      }
5696      if (nDelete) {
5697        solver_->deleteRows(nDelete,added);
5698      }
5699      delete [] added;
5700      delete basis ;
5701    }
5702    // mark so heuristics can tell
5703    int savePass=currentPassNumber_;
5704    currentPassNumber_=999999;
5705    double * newSolution = new double [numberColumns] ;
5706    double heuristicValue = getCutoff() ;
5707    int found = -1; // no solution found
5708    for (int i = 0;i<numberHeuristics_;i++) {
5709      // see if heuristic will do anything
5710      double saveValue = heuristicValue ;
5711      int ifSol = heuristic_[i]->solution(heuristicValue,
5712                                          newSolution);
5713      if (ifSol>0) {
5714        // better solution found
5715        found = i ;
5716        incrementUsed(newSolution);
5717      } else {
5718        heuristicValue = saveValue ;
5719      }
5720    }
5721    currentPassNumber_=savePass;
5722    if (found >= 0) { 
5723      phase_=4;
5724      incrementUsed(newSolution);
5725      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
5726      lastHeuristic_ = heuristic_[found];
5727    }
5728    delete [] newSolution ;
5729  }
5730  // Up change due to cuts
5731  if (feasible)
5732    sumChangeObjective2_ += solver_->getObjValue()*solver_->getObjSense()
5733      - objectiveValue ;
5734  //if ((numberNodes_%100)==0)
5735  //printf("XXb sum obj changed by %g\n",sumChangeObjective2_);
5736/*
5737  End of cut generation loop.
5738
5739  Now, consider if we want to disable or adjust the frequency of use for any
5740  of the cut generators. If the client specified a positive number for
5741  howOften, it will never change. If the original value was negative, it'll
5742  be converted to 1000000+|howOften|, and this value will be adjusted each
5743  time fullScan is true. Actual cut generation is performed every
5744  howOften%1000000 nodes; the 1000000 offset is just a convenient way to
5745  specify that the frequency is adjustable.
5746
5747  During cut generation, we recorded the number of cuts produced by each
5748  generator for this node. For all cuts, whichGenerator records the generator
5749  that produced a cut.
5750
5751  TODO: All this should probably be hidden in a method of the CbcCutGenerator
5752  class.
5753*/
5754#ifdef NODE_LOG
5755  int fatherNum = (node == NULL) ? -1 : node->nodeNumber();
5756  double value =  (node == NULL) ? -1 : node->branchingObject()->value();
5757  string bigOne = (solver_->getIterationCount() > 30) ? "*******" : "";
5758  string way = (node == NULL) ? "" : (node->branchingObject()->way())==1 ? "Down" : "Up";
5759  std::cout<<"Node "<<numberNodes_<<", father "<<fatherNum<<", #iterations "<<solver_->getIterationCount()<<", sol value : "<<solver_->getObjValue()<<std::endl;
5760#endif
5761  if (fullScan&&numberCutGenerators_) {
5762    /* If cuts just at root node then it will probably be faster to
5763       update matrix and leave all in */
5764    int willBeCutsInTree=0;
5765    double thisObjective = solver_->getObjValue()*direction ;
5766    // get sizes
5767    int numberRowsAdded = solver_->getNumRows()-numberRowsAtStart;
5768    CoinBigIndex numberElementsAdded =  solver_->getNumElements()-numberElementsAtStart ;
5769    double densityOld = ((double) numberElementsAtStart)/((double) numberRowsAtStart);
5770    double densityNew = numberRowsAdded ? ((double) (numberElementsAdded))/((double) numberRowsAdded)
5771      : 0.0;
5772    if (!numberNodes_) {
5773      if (numberRowsAdded)
5774        handler_->message(CBC_CUTS_STATS,messages_)
5775          <<numberRowsAdded
5776          <<densityNew
5777          <<CoinMessageEol ;
5778      if (thisObjective-startObjective<1.0e-5&&numberElementsAdded>0.2*numberElementsAtStart)
5779        willBeCutsInTree=-1;
5780    }
5781    if ((numberRowsAdded>100+0.5*numberRowsAtStart
5782         ||numberElementsAdded>0.5*numberElementsAtStart)
5783        &&(densityNew>200.0&&numberRowsAdded>100&&densityNew>2.0*densityOld)) {
5784      // much bigger
5785      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
5786      //willBeCutsInTree=-1;
5787      //printf("Cuts will be taken off , %d rows added with density %g\n",
5788      //     numberRowsAdded,densityNew);
5789    }
5790    if (densityNew>100.0&&numberRowsAdded>2&&densityNew>2.0*densityOld) {
5791      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
5792      //willBeCutsInTree=-2;
5793      //printf("Density says no cuts ? , %d rows added with density %g\n",
5794      //     numberRowsAdded,densityNew);
5795    }
5796    // Root node or every so often - see what to turn off
5797    int i ;
5798    for (i = 0;i<numberCutGenerators_;i++) {
5799      int howOften = generator_[i]->howOften() ;
5800      if (howOften>-90) 
5801        willBeCutsInTree=0;
5802    }
5803    if (!numberNodes_)
5804      handler_->message(CBC_ROOT,messages_)
5805        <<numberNewCuts_
5806        <<startObjective<<thisObjective
5807        <<currentPassNumber_
5808        <<CoinMessageEol ;
5809    int * count = new int[numberCutGenerators_] ;
5810    memset(count,0,numberCutGenerators_*sizeof(int)) ;
5811    int numberActiveGenerators=0;
5812    for (i = 0;i<numberNewCuts_;i++) {
5813      int iGenerator = whichGenerator_[i];
5814      if (iGenerator>=0&&iGenerator<numberCutGenerators_)
5815        count[iGenerator]++ ;
5816    }
5817    double totalCuts = 0.0 ;
5818    //#define JUST_ACTIVE
5819    for (i = 0;i<numberCutGenerators_;i++) {
5820      if (countRowCuts[i]||countColumnCuts[i])
5821        numberActiveGenerators++;
5822#ifdef JUST_ACTIVE
5823      totalCuts += count[i] + 5.0*countColumnCuts[i] ;
5824#else
5825      totalCuts += countRowCuts[i] + 5.0*countColumnCuts[i] ;
5826#endif
5827    }
5828    double small = (0.5* totalCuts) /
5829      ((double) numberActiveGenerators) ;
5830    for (i = 0;i<numberCutGenerators_;i++) {
5831      int howOften = generator_[i]->howOften() ;
5832      if (willBeCutsInTree<0&&howOften==-98)
5833        howOften =-99;
5834      if (howOften==-98&&generator_[i]->switchOffIfLessThan()<0) {
5835        if (thisObjective-startObjective<0.005*fabs(startObjective)+1.0e-5)
5836          howOften=-99; // switch off
5837        if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5
5838            &&5*solver_->getNumRows()<solver_->getNumCols())
5839          howOften=-99; // switch off
5840      }
5841      if (howOften<-99)
5842        continue ;
5843      if (howOften<0||howOften >= 1000000) {
5844        if( !numberNodes_) {
5845          // If small number switch mostly off
5846#ifdef JUST_ACTIVE
5847          double thisCuts = count[i] + 5.0*countColumnCuts[i] ;
5848#else
5849          double thisCuts = countRowCuts[i] + 5.0*countColumnCuts[i] ;
5850#endif
5851          if (!thisCuts||howOften == -99) {
5852            if (howOften == -99||howOften == -98) 
5853              howOften = -100 ;
5854            else
5855              howOften = 1000000+SCANCUTS; // wait until next time
5856          } else if (thisCuts<small) {
5857            int k = (int) sqrt(small/thisCuts) ;
5858            if (howOften!=-98)
5859              howOften = k+1000000 ;
5860            else
5861              howOften=-100;
5862          } else {
5863            howOften = 1+1000000 ;
5864            if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
5865              generator_[i]->setWhatDepth(5);
5866          }
5867        }
5868        // If cuts useless switch off
5869        if (numberNodes_>=100000&&sumChangeObjective1_>2.0e2*(sumChangeObjective2_+1.0e-12)) {
5870          howOften = 1000000+SCANCUTS; // wait until next time
5871          //printf("switch off cut %d due to lack of use\n",i);
5872        }
5873      }
5874      if (howOften>=0&&generator_[i]->generator()->mayGenerateRowCutsInTree())
5875        willBeCutsInTree=1;
5876       
5877      generator_[i]->setHowOften(howOften) ;
5878      if (howOften>=1000000&&howOften<2000000&&0) {
5879        // Go to depth
5880        int bias=1;
5881        if (howOften==1+1000000