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

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

for bonmin

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