source: branches/devel/Cbc/src/CbcModel.cpp @ 473

Last change on this file since 473 was 473, checked in by forrest, 13 years ago

for osibranching

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