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

Last change on this file since 669 was 669, checked in by forrest, 12 years ago

still trying to get round compiler error

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