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

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

timing

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