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

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

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