source: trunk/CbcModel.cpp @ 295

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

for ampl

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