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

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

for Pierre

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