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

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

try and fix hang on max nodes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 357.5 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    //bool stopped = (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
2229    //        numberSolutions_ < intParam_[CbcMaxNumSol] &&
2230    //        totalTime < dblParam_[CbcMaximumSeconds] &&
2231    //        !stoppedOnGap_&&!eventHappened_));
2232    for (i=0;i<numberThreads_;i++) {
2233      while (threadInfo[i].returnCode==0) {
2234        pthread_cond_signal(threadInfo[i].condition2); // unlock
2235        pthread_mutex_lock(&condition_mutex);
2236        struct timespec absTime;
2237        clock_gettime(CLOCK_REALTIME,&absTime);
2238        absTime.tv_nsec += 1000000; // millisecond
2239        if (absTime.tv_nsec>=1000000000) {
2240          absTime.tv_nsec -= 1000000000;
2241          absTime.tv_sec++;
2242        }
2243        pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
2244        clock_gettime(CLOCK_REALTIME,&absTime);
2245        pthread_mutex_unlock(&condition_mutex);
2246      }
2247      pthread_cond_signal(threadInfo[i].condition2); // unlock
2248      pthread_mutex_lock(&condition_mutex); // not sure necessary but have had one hang on interrupt
2249      threadModel[i]->numberThreads_=0; // say exit
2250      threadInfo[i].returnCode=0;
2251      pthread_mutex_unlock(&condition_mutex);
2252      pthread_cond_signal(threadInfo[i].condition2); // unlock
2253      //if (!stopped)
2254        pthread_join(threadId[i],NULL);
2255        //else
2256        //pthread_kill(threadId[i]); // kill rather than try and synchronize
2257      threadModel[i]->moveToModel(this,2);
2258      pthread_mutex_destroy (threadInfo[i].mutex2);
2259      pthread_cond_destroy (threadInfo[i].condition2);
2260      assert (threadInfo[i].numberTimesLocked==threadInfo[i].numberTimesUnlocked);
2261      handler_->message(CBC_THREAD_STATS,messages_)
2262        <<"Thread";
2263      handler_->printing(true)
2264        <<i<<threadCount[i]<<threadInfo[i].timeWaitingToStart;
2265      handler_->printing(goodTimer)<<threadInfo[i].timeInThread;
2266      handler_->printing(false)<<0.0;
2267      handler_->printing(true)<<threadInfo[i].numberTimesLocked
2268        <<threadInfo[i].timeLocked<<threadInfo[i].timeWaitingToLock
2269        <<CoinMessageEol;
2270    }
2271    assert (threadInfo[numberThreads_].numberTimesLocked==threadInfo[numberThreads_].numberTimesUnlocked);
2272    handler_->message(CBC_THREAD_STATS,messages_)
2273      <<"Main thread";
2274    handler_->printing(false)<<0<<0<<0.0;
2275    handler_->printing(false)<<0.0;
2276    handler_->printing(true)<<timeWaiting;
2277    handler_->printing(true)<<threadInfo[numberThreads_].numberTimesLocked
2278      <<threadInfo[numberThreads_].timeLocked<<threadInfo[numberThreads_].timeWaitingToLock
2279      <<CoinMessageEol;
2280    pthread_mutex_destroy (&mutex);
2281    pthread_cond_destroy (&condition_main);
2282    pthread_mutex_destroy (&condition_mutex);
2283    // delete models (here in case some point to others)
2284    for (i=0;i<numberThreads_;i++) {
2285      delete threadModel[i];
2286    }
2287    delete [] mutex2;
2288    delete [] condition2;
2289    delete [] threadId;
2290    delete [] threadInfo;
2291    delete [] threadModel;
2292    delete [] threadCount;
2293    mutex_=NULL;
2294    // adjust time to allow for children on some systems
2295    dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
2296  }
2297#endif
2298/*
2299  End of the non-abort actions. The next block of code is executed if we've
2300  aborted because we hit one of the limits. Clean up by deleting the live set
2301  and break out of the node processing loop. Note that on an abort, node may
2302  have been pushed back onto the tree for further processing, in which case
2303  it'll be deleted in cleanTree. We need to check.
2304*/
2305    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
2306        numberSolutions_ < intParam_[CbcMaxNumSol] &&
2307        totalTime < dblParam_[CbcMaximumSeconds] &&
2308        !stoppedOnGap_&&!eventHappened_)) {
2309      if (tree_->size())
2310        tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
2311      delete nextRowCut_;
2312      if (stoppedOnGap_)
2313        { messageHandler()->message(CBC_GAP,messages())
2314          << bestObjective_-bestPossibleObjective_
2315          << dblParam_[CbcAllowableGap]
2316          << dblParam_[CbcAllowableFractionGap]*100.0
2317          << CoinMessageEol ;
2318        secondaryStatus_ = 2;
2319        status_ = 0 ; }
2320        else
2321          if (isNodeLimitReached())
2322            { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
2323            secondaryStatus_ = 3;
2324            status_ = 1 ; }
2325          else
2326        if (totalTime >= dblParam_[CbcMaximumSeconds])
2327          { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
2328          secondaryStatus_ = 4;
2329          status_ = 1 ; }
2330        else
2331          if (eventHappened_)
2332            { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
2333            secondaryStatus_ = 5;
2334            status_ = 5 ; }
2335          else
2336            { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
2337            secondaryStatus_ = 6;
2338            status_ = 1 ; }
2339    }
2340/*
2341  That's it, we've exhausted the search tree, or broken out of the loop because
2342  we hit some limit on evaluation.
2343
2344  We may have got an intelligent tree so give it one more chance
2345*/
2346  // Tell solver we are not in Branch and Cut
2347  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
2348  tree_->endSearch();
2349  //  If we did any sub trees - did we give up on any?
2350  if ( numberStoppedSubTrees_)
2351    status_=1;
2352  if (!status_) {
2353    // Set best possible unless stopped on gap
2354    if(secondaryStatus_ != 2)
2355      bestPossibleObjective_=bestObjective_;
2356    handler_->message(CBC_END_GOOD,messages_)
2357      << bestObjective_ << numberIterations_ << numberNodes_<<getCurrentSeconds()
2358      << CoinMessageEol ;
2359  } else {
2360    handler_->message(CBC_END,messages_)
2361      << bestObjective_ <<bestPossibleObjective_
2362      << numberIterations_ << numberNodes_<<getCurrentSeconds()
2363      << CoinMessageEol ;
2364  }
2365  if (numberStrongIterations_)
2366    handler_->message(CBC_STRONG_STATS,messages_)
2367      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
2368      << strongInfo_[1] << CoinMessageEol ;
2369  handler_->message(CBC_OTHER_STATS,messages_)
2370      << maximumDepthActual_
2371      << numberDJFixed_ << CoinMessageEol ;
2372  if (statistics_) {
2373    // report in some way
2374    int * lookup = new int[numberObjects_];
2375    int i;
2376    for (i=0;i<numberObjects_;i++) 
2377      lookup[i]=-1;
2378    bool goodIds=true;
2379    for (i=0;i<numberObjects_;i++) {
2380      int iColumn = object_[i]->columnNumber();
2381      if(iColumn>=0&&iColumn<numberColumns) {
2382        if (lookup[i]==-1) {
2383          lookup[i]=iColumn;
2384        } else {
2385          goodIds=false;
2386          break;
2387        }
2388      } else {
2389        goodIds=false;
2390        break;
2391      }
2392    }
2393    if (!goodIds) {
2394      delete [] lookup;
2395      lookup=NULL;
2396    }
2397    if (doStatistics==3) {
2398      printf("  node parent depth column   value                    obj      inf\n");
2399      for ( i=0;i<numberNodes2_;i++) {
2400        statistics_[i]->print(lookup);
2401      }
2402    }
2403    if (doStatistics>1) {
2404      // Find last solution
2405      int k;
2406      for (k=numberNodes2_-1;k>=0;k--) {
2407        if (statistics_[k]->endingObjective()!=COIN_DBL_MAX&&
2408            !statistics_[k]->endingInfeasibility())
2409          break;
2410      }
2411      if (k>=0) {
2412        int depth=statistics_[k]->depth();
2413        int * which = new int[depth+1];
2414        for (i=depth;i>=0;i--) {
2415          which[i]=k;
2416          k=statistics_[k]->parentNode();
2417        }
2418        printf("  node parent depth column   value                    obj      inf\n");
2419        for (i=0;i<=depth;i++) {
2420          statistics_[which[i]]->print(lookup);
2421        }
2422        delete [] which;
2423      }
2424    }
2425    // now summary
2426    int maxDepth=0;
2427    double averageSolutionDepth=0.0;
2428    int numberSolutions=0;
2429    double averageCutoffDepth=0.0;
2430    double averageSolvedDepth=0.0;
2431    int numberCutoff=0;
2432    int numberDown=0;
2433    int numberFirstDown=0;
2434    double averageInfDown=0.0;
2435    double averageObjDown=0.0;
2436    int numberCutoffDown=0;
2437    int numberUp=0;
2438    int numberFirstUp=0;
2439    double averageInfUp=0.0;
2440    double averageObjUp=0.0;
2441    int numberCutoffUp=0;
2442    double averageNumberIterations1=0.0;
2443    double averageValue=0.0;
2444    for ( i=0;i<numberNodes2_;i++) {
2445      int depth =  statistics_[i]->depth(); 
2446      int way =  statistics_[i]->way(); 
2447      double value = statistics_[i]->value(); 
2448      double startingObjective =  statistics_[i]->startingObjective(); 
2449      int startingInfeasibility = statistics_[i]->startingInfeasibility(); 
2450      double endingObjective = statistics_[i]->endingObjective(); 
2451      int endingInfeasibility = statistics_[i]->endingInfeasibility(); 
2452      maxDepth = CoinMax(depth,maxDepth);
2453      // Only for completed
2454      averageNumberIterations1 += statistics_[i]->numberIterations();
2455      averageValue += value;
2456      if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) {
2457        numberSolutions++;
2458        averageSolutionDepth += depth;
2459      }
2460      if (endingObjective==COIN_DBL_MAX) {
2461        numberCutoff++;
2462        averageCutoffDepth += depth;
2463        if (way<0) {
2464          numberDown++;
2465          numberCutoffDown++;
2466          if (way==-1)
2467            numberFirstDown++;
2468        } else {
2469          numberUp++;
2470          numberCutoffUp++;
2471          if (way==1)
2472            numberFirstUp++;
2473        }
2474      } else {
2475        averageSolvedDepth += depth;
2476        if (way<0) {
2477          numberDown++;
2478          averageInfDown += startingInfeasibility-endingInfeasibility;
2479          averageObjDown += endingObjective-startingObjective;
2480          if (way==-1)
2481            numberFirstDown++;
2482        } else {
2483          numberUp++;
2484          averageInfUp += startingInfeasibility-endingInfeasibility;
2485          averageObjUp += endingObjective-startingObjective;
2486          if (way==1)
2487            numberFirstUp++;
2488        }
2489      }
2490    }
2491    // Now print
2492    if (numberSolutions)
2493      averageSolutionDepth /= (double) numberSolutions;
2494    int numberSolved = numberNodes2_-numberCutoff;
2495    double averageNumberIterations2=numberIterations_-averageNumberIterations1
2496      -numberIterationsAtContinuous;
2497    if(numberCutoff) {
2498      averageCutoffDepth /= (double) numberCutoff;
2499      averageNumberIterations2 /= (double) numberCutoff;
2500    }
2501    if (numberNodes2_) 
2502      averageValue /= (double) numberNodes2_;
2503    if (numberSolved) {
2504      averageNumberIterations1 /= (double) numberSolved;
2505      averageSolvedDepth /= (double) numberSolved;
2506    }
2507    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
2508           numberSolutions,averageSolutionDepth);
2509    printf("average value of variable being branched on was %g\n",
2510           averageValue);
2511    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
2512           numberCutoff,averageCutoffDepth,averageNumberIterations2);
2513    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
2514           numberSolved,averageSolvedDepth,averageNumberIterations1);
2515    if (numberDown) {
2516      averageInfDown /= (double) numberDown;
2517      averageObjDown /= (double) numberDown;
2518    }
2519    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
2520           numberDown,numberFirstDown,numberDown-numberFirstDown,numberCutoffDown,
2521           averageInfDown,averageObjDown);
2522    if (numberUp) {
2523      averageInfUp /= (double) numberUp;
2524      averageObjUp /= (double) numberUp;
2525    }
2526    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
2527           numberUp,numberFirstUp,numberUp-numberFirstUp,numberCutoffUp,
2528           averageInfUp,averageObjUp);
2529    for ( i=0;i<numberNodes2_;i++) 
2530      delete statistics_[i];
2531    delete [] statistics_;
2532    statistics_=NULL;
2533    maximumStatistics_=0;
2534    delete [] lookup;
2535  }
2536/*
2537  If we think we have a solution, restore and confirm it with a call to
2538  setBestSolution().  We need to reset the cutoff value so as not to fathom
2539  the solution on bounds.  Note that calling setBestSolution( ..., true)
2540  leaves the continuousSolver_ bounds vectors fixed at the solution value.
2541
2542  Running resolve() here is a failsafe --- setBestSolution has already
2543  reoptimised using the continuousSolver_. If for some reason we fail to
2544  prove optimality, run the problem again after instructing the solver to
2545  tell us more.
2546
2547  If all looks good, replace solver_ with continuousSolver_, so that the
2548  outside world will be able to obtain information about the solution using
2549  public methods.
2550*/
2551  if (bestSolution_&&(solverCharacteristics_->solverType()<2||solverCharacteristics_->solverType()==4)) 
2552  { setCutoff(1.0e50) ; // As best solution should be worse than cutoff
2553    phase_=5;
2554    double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
2555    if ((specialOptions_&4)!=0)
2556      bestObjective_ += 100.0*increment+1.0e-3; // only set if we are going to solve
2557    setBestSolution(CBC_END_SOLUTION,bestObjective_,bestSolution_,true) ;
2558    continuousSolver_->resolve() ;
2559    if (!continuousSolver_->isProvenOptimal())
2560    { continuousSolver_->messageHandler()->setLogLevel(2) ;
2561      continuousSolver_->initialSolve() ; }
2562    delete solver_ ;
2563    // above deletes solverCharacteristics_
2564    solverCharacteristics_ = NULL;
2565    solver_ = continuousSolver_ ;
2566    setPointers(solver_);
2567    continuousSolver_ = NULL ; }
2568/*
2569  Clean up dangling objects. continuousSolver_ may already be toast.
2570*/
2571  delete lastws ;
2572  delete [] whichGenerator_ ;
2573  whichGenerator_=NULL;
2574  delete [] lowerBefore ;
2575  delete [] upperBefore ;
2576  delete [] walkback_ ;
2577  walkback_ = NULL ;
2578  delete [] addedCuts_ ;
2579  addedCuts_ = NULL ;
2580  //delete persistentInfo;
2581  // Get rid of characteristics
2582  solverCharacteristics_=NULL;
2583  if (continuousSolver_)
2584  { delete continuousSolver_ ;
2585    continuousSolver_ = NULL ; }
2586/*
2587  Destroy global cuts by replacing with an empty OsiCuts object.
2588*/
2589  globalCuts_= OsiCuts() ;
2590  if (!bestSolution_) {
2591    // make sure lp solver is infeasible
2592    int numberColumns = solver_->getNumCols();
2593    const double * columnLower = solver_->getColLower();
2594    int iColumn;
2595    for (iColumn=0;iColumn<numberColumns;iColumn++) {
2596      if (solver_->isInteger(iColumn)) 
2597        solver_->setColUpper(iColumn,columnLower[iColumn]);
2598    }
2599    solver_->initialSolve();
2600  }
2601  if (strategy_&&strategy_->preProcessState()>0) {
2602    // undo preprocessing
2603    CglPreProcess * process = strategy_->process();
2604    assert (process);
2605    int n = originalSolver->getNumCols();
2606    if (bestSolution_) {
2607      delete [] bestSolution_;
2608      bestSolution_ = new double [n];
2609      process->postProcess(*solver_);
2610    }
2611    strategy_->deletePreProcess();
2612    // Solution now back in originalSolver
2613    delete solver_;
2614    solver_=originalSolver;
2615    if (bestSolution_)
2616      memcpy(bestSolution_,solver_->getColSolution(),n*sizeof(double));
2617    // put back original objects if there were any
2618    if (originalObject) {
2619      int iColumn;
2620      assert (ownObjects_);
2621      for (iColumn=0;iColumn<numberObjects_;iColumn++) 
2622        delete object_[iColumn];
2623      delete [] object_;
2624      numberObjects_ = numberOriginalObjects;
2625      object_=originalObject;
2626      delete [] integerVariable_;
2627      numberIntegers_=0;
2628      for (iColumn=0;iColumn<n;iColumn++) {
2629        if (solver_->isInteger(iColumn))
2630          numberIntegers_++;
2631      }
2632      integerVariable_ = new int[numberIntegers_];
2633      numberIntegers_=0;
2634      for (iColumn=0;iColumn<n;iColumn++) {
2635        if (solver_->isInteger(iColumn))
2636            integerVariable_[numberIntegers_++]=iColumn;
2637      }
2638    }
2639  }
2640#ifdef CLP_QUICK_OPTIONS
2641  {
2642    OsiClpSolverInterface * clpSolver
2643      = dynamic_cast<OsiClpSolverInterface *> (solver_);
2644    if (clpSolver) {
2645      // Try and re-use regions   
2646      ClpSimplex * simplex = clpSolver->getModelPtr();
2647      simplex->setPersistenceFlag(0);
2648      clpSolver->deleteScaleFactors();
2649      clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~131072));
2650      simplex->setSpecialOptions(simplex->specialOptions()&(~131072));
2651      simplex->setAlphaAccuracy(-1.0);
2652      //clpSolver->setSpecialOptions((clpSolver->specialOptions()&~128)|65536);
2653    }
2654  }
2655#endif
2656  return ;
2657 }
2658
2659
2660// Solve the initial LP relaxation
2661void 
2662CbcModel::initialSolve() 
2663{
2664  assert (solver_);
2665  assert (!solverCharacteristics_);
2666  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
2667  if (solverCharacteristics) {
2668    solverCharacteristics_ = solverCharacteristics;
2669  } else {
2670    // replace in solver
2671    OsiBabSolver defaultC;
2672    solver_->setAuxiliaryInfo(&defaultC);
2673    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
2674  }
2675  solverCharacteristics_->setSolver(solver_);
2676  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
2677  solver_->initialSolve();
2678  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
2679  // But set up so Jon Lee will be happy
2680  status_=-1;
2681  secondaryStatus_ = -1;
2682  originalContinuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
2683  delete [] continuousSolution_;
2684  continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
2685                                             solver_->getNumCols());
2686  setPointers(solver_);
2687  solverCharacteristics_ = NULL;
2688}
2689
2690/*! \brief Get an empty basis object
2691
2692  Return an empty CoinWarmStartBasis object with the requested capacity,
2693  appropriate for the current solver. The object is cloned from the object
2694  cached as emptyWarmStart_. If there is no cached object, the routine
2695  queries the solver for a warm start object, empties it, and caches the
2696  result.
2697*/
2698
2699CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
2700
2701{ CoinWarmStartBasis *emptyBasis ;
2702/*
2703  Acquire an empty basis object, if we don't yet have one.
2704*/
2705  if (emptyWarmStart_ == 0)
2706  { if (solver_ == 0)
2707    { throw CoinError("Cannot construct basis without solver!",
2708                      "getEmptyBasis","CbcModel") ; }
2709    emptyBasis =
2710        dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
2711    if (emptyBasis == 0)
2712    { throw CoinError(
2713        "Solver does not appear to use a basis-oriented warm start.",
2714        "getEmptyBasis","CbcModel") ; }
2715    emptyBasis->setSize(0,0) ;
2716    emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ; }
2717/*
2718  Clone the empty basis object, resize it as requested, and return.
2719*/
2720  emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
2721  assert(emptyBasis) ;
2722  if (ns != 0 || na != 0) emptyBasis->setSize(ns,na) ;
2723
2724  return (emptyBasis) ; }
2725   
2726
2727/** Default Constructor
2728
2729  Creates an empty model without an associated solver.
2730*/
2731CbcModel::CbcModel() 
2732
2733:
2734  solver_(NULL),
2735  ourSolver_(true),
2736  continuousSolver_(NULL),
2737  referenceSolver_(NULL),
2738  defaultHandler_(true),
2739  emptyWarmStart_(NULL),
2740  bestObjective_(COIN_DBL_MAX),
2741  bestPossibleObjective_(COIN_DBL_MAX),
2742  sumChangeObjective1_(0.0),
2743  sumChangeObjective2_(0.0),
2744  bestSolution_(NULL),
2745  currentSolution_(NULL),
2746  testSolution_(NULL),
2747  minimumDrop_(1.0e-4),
2748  numberSolutions_(0),
2749  stateOfSearch_(0),
2750  hotstartSolution_(NULL),
2751  hotstartPriorities_(NULL),
2752  numberHeuristicSolutions_(0),
2753  numberNodes_(0),
2754  numberNodes2_(0),
2755  numberIterations_(0),
2756  status_(-1),
2757  secondaryStatus_(-1),
2758  numberIntegers_(0),
2759  numberRowsAtContinuous_(0),
2760  maximumNumberCuts_(0),
2761  phase_(0),
2762  currentNumberCuts_(0),
2763  maximumDepth_(0),
2764  walkback_(NULL),
2765  addedCuts_(NULL),
2766  nextRowCut_(NULL),
2767  currentNode_(NULL),
2768  integerVariable_(NULL),
2769  integerInfo_(NULL),
2770  continuousSolution_(NULL),
2771  usedInSolution_(NULL),
2772  specialOptions_(0),
2773  subTreeModel_(NULL),
2774  numberStoppedSubTrees_(0),
2775  mutex_(NULL),
2776  presolve_(0),
2777  numberStrong_(5),
2778  numberBeforeTrust_(0),
2779  numberPenalties_(20),
2780  penaltyScaleFactor_(3.0),
2781  numberAnalyzeIterations_(0),
2782  analyzeResults_(NULL),
2783  numberInfeasibleNodes_(0),
2784  problemType_(0),
2785  printFrequency_(0),
2786  numberCutGenerators_(0),
2787  generator_(NULL),
2788  virginGenerator_(NULL),
2789  numberHeuristics_(0),
2790  heuristic_(NULL),
2791  lastHeuristic_(NULL),
2792  eventHandler_(0),
2793  numberObjects_(0),
2794  object_(NULL),
2795  ownObjects_(true),
2796  originalColumns_(NULL),
2797  howOftenGlobalScan_(1),
2798  numberGlobalViolations_(0),
2799  continuousObjective_(COIN_DBL_MAX),
2800  originalContinuousObjective_(COIN_DBL_MAX),
2801  continuousInfeasibilities_(INT_MAX),
2802  maximumCutPassesAtRoot_(20),
2803  maximumCutPasses_(10),
2804  preferredWay_(0),
2805  currentPassNumber_(0),
2806  maximumWhich_(1000),
2807  whichGenerator_(NULL),
2808  maximumStatistics_(0),
2809  statistics_(NULL),
2810  maximumDepthActual_(0),
2811  numberDJFixed_(0.0),
2812  probingInfo_(NULL),
2813  numberFixedAtRoot_(0),
2814  numberFixedNow_(0),
2815  stoppedOnGap_(false),
2816  eventHappened_(false),
2817  numberLongStrong_(0),
2818  numberOldActiveCuts_(0),
2819  numberNewCuts_(0),
2820  sizeMiniTree_(0),
2821  searchStrategy_(-1),
2822  numberStrongIterations_(0),
2823  resolveAfterTakeOffCuts_(true),
2824#if NEW_UPDATE_OBJECT>1
2825  numberUpdateItems_(0),
2826  maximumNumberUpdateItems_(0),
2827  updateItems_(NULL),
2828#endif
2829  numberThreads_(0),
2830  threadMode_(0)
2831{
2832  memset(intParam_,0,sizeof(intParam_));
2833  intParam_[CbcMaxNumNode] = 2147483647;
2834  intParam_[CbcMaxNumSol] = 9999999;
2835  intParam_[CbcFathomDiscipline] = 0;
2836
2837  dblParam_[CbcIntegerTolerance] = 1e-6;
2838  dblParam_[CbcInfeasibilityWeight] = 0.0;
2839  dblParam_[CbcCutoffIncrement] = 1e-5;
2840  dblParam_[CbcAllowableGap] = 1.0e-10;
2841  dblParam_[CbcAllowableFractionGap] = 0.0;
2842  dblParam_[CbcMaximumSeconds] = 1.0e100;
2843  dblParam_[CbcCurrentCutoff] = 1.0e100;
2844  dblParam_[CbcOptimizationDirection] = 1.0;
2845  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
2846  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
2847  dblParam_[CbcStartSeconds] = 0.0;
2848  strongInfo_[0]=0;
2849  strongInfo_[1]=0;
2850  strongInfo_[2]=0;
2851  solverCharacteristics_ = NULL;
2852  nodeCompare_=new CbcCompareDefault();;
2853  problemFeasibility_=new CbcFeasibilityBase();
2854  tree_= new CbcTree();
2855  branchingMethod_=NULL;
2856  cutModifier_=NULL;
2857  strategy_=NULL;
2858  parentModel_=NULL;
2859  cbcColLower_ = NULL;
2860  cbcColUpper_ = NULL;
2861  cbcRowLower_ = NULL;
2862  cbcRowUpper_ = NULL;
2863  cbcColSolution_ = NULL;
2864  cbcRowPrice_ = NULL;
2865  cbcReducedCost_ = NULL;
2866  cbcRowActivity_ = NULL;
2867  appData_=NULL;
2868  handler_ = new CoinMessageHandler();
2869  handler_->setLogLevel(2);
2870  messages_ = CbcMessage();
2871  eventHandler_ = new CbcEventHandler() ;
2872}
2873
2874/** Constructor from solver.
2875
2876  Creates a model complete with a clone of the solver passed as a parameter.
2877*/
2878
2879CbcModel::CbcModel(const OsiSolverInterface &rhs)
2880:
2881  continuousSolver_(NULL),
2882  referenceSolver_(NULL),
2883  defaultHandler_(true),
2884  emptyWarmStart_(NULL),
2885  bestObjective_(COIN_DBL_MAX),
2886  bestPossibleObjective_(COIN_DBL_MAX),
2887  sumChangeObjective1_(0.0),
2888  sumChangeObjective2_(0.0),
2889  minimumDrop_(1.0e-4),
2890  numberSolutions_(0),
2891  stateOfSearch_(0),
2892  hotstartSolution_(NULL),
2893  hotstartPriorities_(NULL),
2894  numberHeuristicSolutions_(0),
2895  numberNodes_(0),
2896  numberNodes2_(0),
2897  numberIterations_(0),
2898  status_(-1),
2899  secondaryStatus_(-1),
2900  numberRowsAtContinuous_(0),
2901  maximumNumberCuts_(0),
2902  phase_(0),
2903  currentNumberCuts_(0),
2904  maximumDepth_(0),
2905  walkback_(NULL),
2906  addedCuts_(NULL),
2907  nextRowCut_(NULL),
2908  currentNode_(NULL),
2909  integerInfo_(NULL),
2910  specialOptions_(0),
2911  subTreeModel_(NULL),
2912  numberStoppedSubTrees_(0),
2913  mutex_(NULL),
2914  presolve_(0),
2915  numberStrong_(5),
2916  numberBeforeTrust_(0),
2917  numberPenalties_(20),
2918  penaltyScaleFactor_(3.0),
2919  numberAnalyzeIterations_(0),
2920  analyzeResults_(NULL),
2921  numberInfeasibleNodes_(0),
2922  problemType_(0),
2923  printFrequency_(0),
2924  numberCutGenerators_(0),
2925  generator_(NULL),
2926  virginGenerator_(NULL),
2927  numberHeuristics_(0),
2928  heuristic_(NULL),
2929  lastHeuristic_(NULL),
2930  eventHandler_(0),
2931  numberObjects_(0),
2932  object_(NULL),
2933  ownObjects_(true),
2934  originalColumns_(NULL),
2935  howOftenGlobalScan_(1),
2936  numberGlobalViolations_(0),
2937  continuousObjective_(COIN_DBL_MAX),
2938  originalContinuousObjective_(COIN_DBL_MAX),
2939  continuousInfeasibilities_(INT_MAX),
2940  maximumCutPassesAtRoot_(20),
2941  maximumCutPasses_(10),
2942  preferredWay_(0),
2943  currentPassNumber_(0),
2944  maximumWhich_(1000),
2945  whichGenerator_(NULL),
2946  maximumStatistics_(0),
2947  statistics_(NULL),
2948  maximumDepthActual_(0),
2949  numberDJFixed_(0.0),
2950  probingInfo_(NULL),
2951  numberFixedAtRoot_(0),
2952  numberFixedNow_(0),
2953  stoppedOnGap_(false),
2954  eventHappened_(false),
2955  numberLongStrong_(0),
2956  numberOldActiveCuts_(0),
2957  numberNewCuts_(0),
2958  sizeMiniTree_(0),
2959  searchStrategy_(-1),
2960  numberStrongIterations_(0),
2961  resolveAfterTakeOffCuts_(true),
2962#if NEW_UPDATE_OBJECT>1
2963  numberUpdateItems_(0),
2964  maximumNumberUpdateItems_(0),
2965  updateItems_(NULL),
2966#endif
2967  numberThreads_(0),
2968  threadMode_(0)
2969{
2970  memset(intParam_,0,sizeof(intParam_));
2971  intParam_[CbcMaxNumNode] = 2147483647;
2972  intParam_[CbcMaxNumSol] = 9999999;
2973  intParam_[CbcFathomDiscipline] = 0;
2974
2975  dblParam_[CbcIntegerTolerance] = 1e-6;
2976  dblParam_[CbcInfeasibilityWeight] = 0.0;
2977  dblParam_[CbcCutoffIncrement] = 1e-5;
2978  dblParam_[CbcAllowableGap] = 1.0e-10;
2979  dblParam_[CbcAllowableFractionGap] = 0.0;
2980  dblParam_[CbcMaximumSeconds] = 1.0e100;
2981  dblParam_[CbcCurrentCutoff] = 1.0e100;
2982  dblParam_[CbcOptimizationDirection] = 1.0;
2983  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
2984  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
2985  dblParam_[CbcStartSeconds] = 0.0;
2986  strongInfo_[0]=0;
2987  strongInfo_[1]=0;
2988  strongInfo_[2]=0;
2989  solverCharacteristics_ = NULL;
2990
2991  nodeCompare_=new CbcCompareDefault();;
2992  problemFeasibility_=new CbcFeasibilityBase();
2993  tree_= new CbcTree();
2994  branchingMethod_=NULL;
2995  cutModifier_=NULL;
2996  strategy_=NULL;
2997  parentModel_=NULL;
2998  appData_=NULL;
2999  handler_ = new CoinMessageHandler();
3000  handler_->setLogLevel(2);
3001  messages_ = CbcMessage();
3002  eventHandler_ = new CbcEventHandler() ;
3003  solver_ = rhs.clone();
3004  referenceSolver_ = solver_->clone();
3005  ourSolver_ = true ;
3006  cbcColLower_ = NULL;
3007  cbcColUpper_ = NULL;
3008  cbcRowLower_ = NULL;
3009  cbcRowUpper_ = NULL;
3010  cbcColSolution_ = NULL;
3011  cbcRowPrice_ = NULL;
3012  cbcReducedCost_ = NULL;
3013  cbcRowActivity_ = NULL;
3014
3015  // Initialize solution and integer variable vectors
3016  bestSolution_ = NULL; // to say no solution found
3017  numberIntegers_=0;
3018  int numberColumns = solver_->getNumCols();
3019  int iColumn;
3020  if (numberColumns) {
3021    // Space for current solution
3022    currentSolution_ = new double[numberColumns];
3023    continuousSolution_ = new double[numberColumns];
3024    usedInSolution_ = new int[numberColumns];
3025    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3026      if( solver_->isInteger(iColumn)) 
3027        numberIntegers_++;
3028    }
3029  } else {
3030    // empty model
3031    currentSolution_=NULL;
3032    continuousSolution_=NULL;
3033    usedInSolution_=NULL;
3034  }
3035  testSolution_=currentSolution_;
3036  if (numberIntegers_) {
3037    integerVariable_ = new int [numberIntegers_];
3038    numberIntegers_=0;
3039    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3040      if( solver_->isInteger(iColumn)) 
3041        integerVariable_[numberIntegers_++]=iColumn;
3042    }
3043  } else {
3044    integerVariable_ = NULL;
3045  }
3046}
3047
3048/*
3049  Assign a solver to the model (model assumes ownership)
3050
3051  The integer variable vector is initialized if it's not already present.
3052  If deleteSolver then current solver deleted (if model owned)
3053
3054  Assuming ownership matches usage in OsiSolverInterface
3055  (cf. assignProblem, loadProblem).
3056
3057  TODO: What to do about solver parameters? A simple copy likely won't do it,
3058        because the SI must push the settings into the underlying solver. In
3059        the context of switching solvers in cbc, this means that command line
3060        settings will get lost. Stash the command line somewhere and reread it
3061        here, maybe?
3062 
3063  TODO: More generally, how much state should be transferred from the old
3064        solver to the new solver? Best perhaps to see how usage develops.
3065        What's done here mimics the CbcModel(OsiSolverInterface) constructor.
3066*/
3067void
3068CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver)
3069
3070{
3071  // resize best solution if exists
3072  if (bestSolution_) {
3073    int nOld = solver_->getNumCols();
3074    int nNew = solver->getNumCols();
3075    if (nNew>nOld) {
3076      double * temp = new double[nNew];
3077      memcpy(temp,bestSolution_,nOld*sizeof(double));
3078      memset(temp+nOld,0,(nNew-nOld)*sizeof(double));
3079      delete [] bestSolution_;
3080      bestSolution_=temp;
3081    }
3082  }
3083  // Keep the current message level for solver (if solver exists)
3084  if (solver_)
3085    solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
3086
3087  if (ourSolver_&&deleteSolver) delete solver_ ;
3088  solver_ = solver;
3089  solver = NULL ;
3090  ourSolver_ = true ;
3091/*
3092  Basis information is solver-specific.
3093*/
3094  if (emptyWarmStart_)
3095  { delete emptyWarmStart_  ;
3096    emptyWarmStart_ = 0 ; }
3097/*
3098  Initialize integer variable vector.
3099*/
3100  numberIntegers_=0;
3101  int numberColumns = solver_->getNumCols();
3102  int iColumn;
3103  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3104    if( solver_->isInteger(iColumn)) 
3105      numberIntegers_++;
3106  }
3107  if (numberIntegers_) {
3108    delete [] integerVariable_;
3109    integerVariable_ = new int [numberIntegers_];
3110    numberIntegers_=0;
3111    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3112      if( solver_->isInteger(iColumn)) 
3113        integerVariable_[numberIntegers_++]=iColumn;
3114    }
3115  } else {
3116    integerVariable_ = NULL;
3117  }
3118
3119  return ;
3120}
3121
3122// Copy constructor.
3123
3124CbcModel::CbcModel(const CbcModel & rhs, bool noTree)
3125:
3126  continuousSolver_(NULL),
3127  referenceSolver_(NULL),
3128  defaultHandler_(rhs.defaultHandler_),
3129  emptyWarmStart_(NULL),
3130  bestObjective_(rhs.bestObjective_),
3131  bestPossibleObjective_(rhs.bestPossibleObjective_),
3132  sumChangeObjective1_(rhs.sumChangeObjective1_),
3133  sumChangeObjective2_(rhs.sumChangeObjective2_),
3134  minimumDrop_(rhs.minimumDrop_),
3135  numberSolutions_(rhs.numberSolutions_),
3136  stateOfSearch_(rhs.stateOfSearch_),
3137  numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
3138  numberNodes_(rhs.numberNodes_),
3139  numberNodes2_(rhs.numberNodes2_),
3140  numberIterations_(rhs.numberIterations_),
3141  status_(rhs.status_),
3142  secondaryStatus_(rhs.secondaryStatus_),
3143  specialOptions_(rhs.specialOptions_),
3144  subTreeModel_(rhs.subTreeModel_),
3145  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
3146  mutex_(NULL),
3147  presolve_(rhs.presolve_),
3148  numberStrong_(rhs.numberStrong_),
3149  numberBeforeTrust_(rhs.numberBeforeTrust_),
3150  numberPenalties_(rhs.numberPenalties_),
3151  penaltyScaleFactor_(rhs.penaltyScaleFactor_),
3152  numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
3153  analyzeResults_(NULL),
3154  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
3155  problemType_(rhs.problemType_),
3156  printFrequency_(rhs.printFrequency_),
3157  howOftenGlobalScan_(rhs.howOftenGlobalScan_),
3158  numberGlobalViolations_(rhs.numberGlobalViolations_),
3159  continuousObjective_(rhs.continuousObjective_),
3160  originalContinuousObjective_(rhs.originalContinuousObjective_),
3161  continuousInfeasibilities_(rhs.continuousInfeasibilities_),
3162  maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
3163  maximumCutPasses_( rhs.maximumCutPasses_),
3164  preferredWay_(rhs.preferredWay_),
3165  currentPassNumber_(rhs.currentPassNumber_),
3166  maximumWhich_(rhs.maximumWhich_),
3167  whichGenerator_(NULL),
3168  maximumStatistics_(0),
3169  statistics_(NULL),
3170  maximumDepthActual_(0),
3171  numberDJFixed_(0.0),
3172  probingInfo_(NULL),
3173  numberFixedAtRoot_(rhs.numberFixedAtRoot_),
3174  numberFixedNow_(rhs.numberFixedNow_),
3175  stoppedOnGap_(rhs.stoppedOnGap_),
3176  eventHappened_(rhs.eventHappened_),
3177  numberLongStrong_(rhs.numberLongStrong_),
3178  numberOldActiveCuts_(rhs.numberOldActiveCuts_),
3179  numberNewCuts_(rhs.numberNewCuts_),
3180  sizeMiniTree_(rhs.sizeMiniTree_),
3181  searchStrategy_(rhs.searchStrategy_),
3182  numberStrongIterations_(rhs.numberStrongIterations_),
3183  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
3184#if NEW_UPDATE_OBJECT>1
3185  numberUpdateItems_(rhs.numberUpdateItems_),
3186  maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_),
3187  updateItems_(NULL),
3188#endif
3189  numberThreads_(rhs.numberThreads_),
3190  threadMode_(rhs.threadMode_)
3191{
3192  memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
3193  memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
3194  strongInfo_[0]=rhs.strongInfo_[0];
3195  strongInfo_[1]=rhs.strongInfo_[1];
3196  strongInfo_[2]=rhs.strongInfo_[2];
3197  solverCharacteristics_ = NULL;
3198  if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
3199  if (defaultHandler_) {
3200    handler_ = new CoinMessageHandler();
3201    handler_->setLogLevel(2);
3202  } else {
3203    handler_ = rhs.handler_;
3204  }
3205  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3206  numberCutGenerators_ = rhs.numberCutGenerators_;
3207  if (numberCutGenerators_) {
3208    generator_ = new CbcCutGenerator * [numberCutGenerators_];
3209    virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
3210    int i;
3211    for (i=0;i<numberCutGenerators_;i++) {
3212      generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
3213      virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
3214    }
3215  } else {
3216    generator_=NULL;
3217    virginGenerator_=NULL;
3218  }
3219  if (!noTree)
3220    globalCuts_ = rhs.globalCuts_;
3221  numberHeuristics_ = rhs.numberHeuristics_;
3222  if (numberHeuristics_) {
3223    heuristic_ = new CbcHeuristic * [numberHeuristics_];
3224    int i;
3225    for (i=0;i<numberHeuristics_;i++) {
3226      heuristic_[i]=rhs.heuristic_[i]->clone();
3227    }
3228  } else {
3229    heuristic_=NULL;
3230  }
3231  lastHeuristic_ = NULL;
3232  if (rhs.eventHandler_)
3233  { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; }
3234  else
3235  { eventHandler_ = NULL ; }
3236  ownObjects_ = rhs.ownObjects_;
3237  if (ownObjects_) {
3238    numberObjects_=rhs.numberObjects_;
3239    if (numberObjects_) {
3240      object_ = new OsiObject * [numberObjects_];
3241      int i;
3242      for (i=0;i<numberObjects_;i++) {
3243        object_[i]=(rhs.object_[i])->clone();
3244        CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
3245        // Could be OsiObjects
3246        if (obj)
3247          obj->setModel(this);
3248      }
3249    } else {
3250      object_=NULL;
3251    }
3252  } else {
3253    // assume will be redone
3254    numberObjects_=0;
3255    object_=NULL;
3256  }
3257  if (rhs.referenceSolver_)
3258    referenceSolver_ = rhs.referenceSolver_->clone();
3259  else
3260    referenceSolver_=NULL;
3261  if (!noTree||!rhs.continuousSolver_)
3262    solver_ = rhs.solver_->clone();
3263  else
3264    solver_ = rhs.continuousSolver_->clone();
3265  if (rhs.originalColumns_) {
3266    int numberColumns = solver_->getNumCols();
3267    originalColumns_= new int [numberColumns];
3268    memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
3269  } else {
3270    originalColumns_=NULL;
3271  }
3272#if NEW_UPDATE_OBJECT>1
3273  if (maximumNumberUpdateItems_) {
3274    updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
3275    for (int i=0;i<maximumNumberUpdateItems_;i++)
3276      updateItems_[i] = rhs.updateItems_[i];
3277  }
3278#endif
3279  if (maximumWhich_&&rhs.whichGenerator_)
3280    whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
3281  nodeCompare_=rhs.nodeCompare_->clone();
3282  problemFeasibility_=rhs.problemFeasibility_->clone();
3283  tree_= rhs.tree_->clone();
3284  if (rhs.branchingMethod_)
3285    branchingMethod_=rhs.branchingMethod_->clone();
3286  else
3287    branchingMethod_=NULL;
3288  if (rhs.cutModifier_)
3289    cutModifier_=rhs.cutModifier_->clone();
3290  else
3291    cutModifier_=NULL;
3292  cbcColLower_ = NULL;
3293  cbcColUpper_ = NULL;
3294  cbcRowLower_ = NULL;
3295  cbcRowUpper_ = NULL;
3296  cbcColSolution_ = NULL;
3297  cbcRowPrice_ = NULL;
3298  cbcReducedCost_ = NULL;
3299  cbcRowActivity_ = NULL;
3300  if (rhs.strategy_)
3301    strategy_=rhs.strategy_->clone();
3302  else
3303    strategy_=NULL;
3304  parentModel_=rhs.parentModel_;
3305  appData_=rhs.appData_;
3306  messages_ = rhs.messages_;
3307  ourSolver_ = true ;
3308  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3309  numberIntegers_=rhs.numberIntegers_;
3310  if (numberIntegers_) {
3311    integerVariable_ = new int [numberIntegers_];
3312    memcpy(integerVariable_,rhs.integerVariable_,numberIntegers_*sizeof(int));
3313    integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
3314  } else {
3315    integerVariable_ = NULL;
3316    integerInfo_=NULL;
3317  }
3318  if (rhs.hotstartSolution_) {
3319    int numberColumns = solver_->getNumCols();
3320    hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
3321    hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
3322  } else {
3323    hotstartSolution_ = NULL;
3324    hotstartPriorities_ =NULL;
3325  }
3326  if (rhs.bestSolution_&&!noTree) {
3327    int numberColumns = solver_->getNumCols();
3328    bestSolution_ = new double[numberColumns];
3329    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
3330  } else {
3331    bestSolution_=NULL;
3332  }
3333  if (!noTree) {
3334    int numberColumns = solver_->getNumCols();
3335    currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns);
3336    continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns);
3337    usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns);
3338  } else {
3339    currentSolution_=NULL;
3340    continuousSolution_=NULL;
3341    usedInSolution_=NULL;
3342  }
3343  testSolution_=currentSolution_;
3344  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
3345  maximumNumberCuts_=rhs.maximumNumberCuts_;
3346  phase_ = rhs.phase_;
3347  currentNumberCuts_=rhs.currentNumberCuts_;
3348  maximumDepth_= rhs.maximumDepth_;
3349  if (noTree) {
3350    bestObjective_ = COIN_DBL_MAX;
3351    numberSolutions_ =0;
3352    stateOfSearch_= 0;
3353    numberHeuristicSolutions_=0;
3354    numberNodes_=0;
3355    numberNodes2_=0;
3356    numberIterations_=0;
3357    status_=0;
3358    subTreeModel_=NULL;
3359    numberStoppedSubTrees_=0;
3360    continuousObjective_=COIN_DBL_MAX;
3361    originalContinuousObjective_=COIN_DBL_MAX;
3362    continuousInfeasibilities_=INT_MAX;
3363    maximumNumberCuts_=0;
3364    tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_);
3365    bestPossibleObjective_ = COIN_DBL_MAX;
3366  }
3367  // These are only used as temporary arrays so need not be filled
3368  if (maximumNumberCuts_) {
3369    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3370  } else {
3371    addedCuts_ = NULL;
3372  }
3373  nextRowCut_ = NULL;
3374  currentNode_ = NULL;
3375  if (maximumDepth_)
3376    walkback_ = new CbcNodeInfo * [maximumDepth_];
3377  else
3378    walkback_ = NULL;
3379  synchronizeModel();
3380}
3381 
3382// Assignment operator
3383CbcModel & 
3384CbcModel::operator=(const CbcModel& rhs)
3385{
3386  if (this!=&rhs) {
3387    gutsOfDestructor();
3388    if (defaultHandler_)
3389    { delete handler_;
3390      handler_ = NULL; }
3391    defaultHandler_ = rhs.defaultHandler_;
3392    if (defaultHandler_)
3393    { handler_ = new CoinMessageHandler();
3394      handler_->setLogLevel(2); }
3395    else
3396    { handler_ = rhs.handler_; }
3397    messages_ = rhs.messages_;
3398    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3399
3400    delete solver_;
3401    if (rhs.solver_)
3402    { solver_ = rhs.solver_->clone() ; }
3403    else
3404    { solver_ = 0 ; }
3405    ourSolver_ = true ;
3406    delete continuousSolver_ ;
3407    if (rhs.continuousSolver_)
3408    { solver_ = rhs.continuousSolver_->clone() ; }
3409    else
3410    { continuousSolver_ = 0 ; }
3411
3412    if (rhs.referenceSolver_)
3413    { solver_ = rhs.referenceSolver_->clone() ; }
3414    else
3415    { referenceSolver_ = NULL ; }
3416
3417    delete emptyWarmStart_ ;
3418    if (rhs.emptyWarmStart_)
3419    { emptyWarmStart_ = rhs.emptyWarmStart_->clone() ; }
3420    else
3421    { emptyWarmStart_ = 0 ; }
3422
3423    bestObjective_ = rhs.bestObjective_;
3424    bestPossibleObjective_=rhs.bestPossibleObjective_;
3425    sumChangeObjective1_=rhs.sumChangeObjective1_;
3426    sumChangeObjective2_=rhs.sumChangeObjective2_;
3427    delete [] bestSolution_;
3428    if (rhs.bestSolution_) {
3429      int numberColumns = rhs.getNumCols();
3430      bestSolution_ = new double[numberColumns];
3431      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
3432    } else {
3433      bestSolution_=NULL;
3434    }
3435    int numberColumns = solver_->getNumCols();
3436    currentSolution_ = CoinCopyOfArray(rhs.currentSolution_,numberColumns);
3437    continuousSolution_ = CoinCopyOfArray(rhs.continuousSolution_,numberColumns);
3438    usedInSolution_ = CoinCopyOfArray(rhs.usedInSolution_,numberColumns);
3439    testSolution_=currentSolution_;
3440    minimumDrop_ = rhs.minimumDrop_;
3441    numberSolutions_=rhs.numberSolutions_;
3442    stateOfSearch_= rhs.stateOfSearch_;
3443    numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
3444    numberNodes_ = rhs.numberNodes_;
3445    numberNodes2_ = rhs.numberNodes2_;
3446    numberIterations_ = rhs.numberIterations_;
3447    status_ = rhs.status_;
3448    secondaryStatus_ = rhs.secondaryStatus_;
3449    specialOptions_ = rhs.specialOptions_;
3450    subTreeModel_ = rhs.subTreeModel_;
3451    numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
3452    mutex_ = NULL;
3453    presolve_ = rhs.presolve_;
3454    numberStrong_ = rhs.numberStrong_;
3455    numberBeforeTrust_ = rhs.numberBeforeTrust_;
3456    numberPenalties_ = rhs.numberPenalties_;
3457    penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
3458    numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
3459    delete [] analyzeResults_;
3460    analyzeResults_ = NULL;
3461    numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
3462    problemType_ = rhs.problemType_;
3463    printFrequency_ = rhs.printFrequency_;
3464    howOftenGlobalScan_=rhs.howOftenGlobalScan_;
3465    numberGlobalViolations_=rhs.numberGlobalViolations_;
3466    continuousObjective_=rhs.continuousObjective_;
3467    originalContinuousObjective_ = rhs.originalContinuousObjective_;
3468    continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
3469    maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
3470    maximumCutPasses_ = rhs.maximumCutPasses_;
3471    preferredWay_ = rhs.preferredWay_;
3472    currentPassNumber_ = rhs.currentPassNumber_;
3473    memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
3474    memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
3475    globalCuts_ = rhs.globalCuts_;
3476    int i;
3477    for (i=0;i<numberCutGenerators_;i++) {
3478      delete generator_[i];
3479      delete virginGenerator_[i];
3480    }
3481    delete [] generator_;
3482    delete [] virginGenerator_;
3483    delete [] heuristic_;
3484    maximumWhich_ = rhs.maximumWhich_;
3485    delete [] whichGenerator_;
3486    whichGenerator_ = NULL;
3487    if (maximumWhich_&&rhs.whichGenerator_)
3488      whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
3489    for (i=0;i<maximumStatistics_;i++)
3490      delete statistics_[i];
3491    delete [] statistics_;
3492    maximumStatistics_ = 0;
3493    statistics_ = NULL;
3494    delete probingInfo_;
3495    probingInfo_=NULL;
3496    numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
3497    numberFixedNow_ = rhs.numberFixedNow_;
3498    stoppedOnGap_ = rhs.stoppedOnGap_;
3499    eventHappened_ = rhs.eventHappened_;
3500    numberLongStrong_ = rhs.numberLongStrong_;
3501    numberOldActiveCuts_ = rhs.numberOldActiveCuts_;
3502    numberNewCuts_ = rhs.numberNewCuts_;
3503    resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_;
3504#if NEW_UPDATE_OBJECT>1
3505    numberUpdateItems_ = rhs.numberUpdateItems_;
3506    maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_;
3507    delete [] updateItems_;
3508    if (maximumNumberUpdateItems_) {
3509      updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
3510      for (i=0;i<maximumNumberUpdateItems_;i++)
3511        updateItems_[i] = rhs.updateItems_[i];
3512    } else {
3513      updateItems_ = NULL;
3514    }
3515#endif
3516    numberThreads_ = rhs.numberThreads_;
3517    threadMode_ = rhs.threadMode_;
3518    sizeMiniTree_ = rhs.sizeMiniTree_;
3519    searchStrategy_ = rhs.searchStrategy_;
3520    numberStrongIterations_ = rhs.numberStrongIterations_;
3521    strongInfo_[0]=rhs.strongInfo_[0];
3522    strongInfo_[1]=rhs.strongInfo_[1];
3523    strongInfo_[2]=rhs.strongInfo_[2];
3524    solverCharacteristics_ = NULL;
3525    lastHeuristic_ = NULL;
3526    numberCutGenerators_ = rhs.numberCutGenerators_;
3527    if (numberCutGenerators_) {
3528      generator_ = new CbcCutGenerator * [numberCutGenerators_];
3529      virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
3530      int i;
3531      for (i=0;i<numberCutGenerators_;i++) {
3532        generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
3533        virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
3534      }
3535    } else {
3536      generator_=NULL;
3537      virginGenerator_=NULL;
3538    }
3539    numberHeuristics_ = rhs.numberHeuristics_;
3540    if (numberHeuristics_) {
3541      heuristic_ = new CbcHeuristic * [numberHeuristics_];
3542      memcpy(heuristic_,rhs.heuristic_,
3543             numberHeuristics_*sizeof(CbcHeuristic *));
3544    } else {
3545      heuristic_=NULL;
3546    }
3547    lastHeuristic_ = NULL;
3548    if (eventHandler_)
3549      delete eventHandler_ ;
3550    if (rhs.eventHandler_)
3551    { eventHandler_ = new CbcEventHandler(*rhs.eventHandler_) ; }
3552    else
3553    { eventHandler_ = NULL ; }
3554    if (ownObjects_) {
3555      for (i=0;i<numberObjects_;i++)
3556        delete object_[i];
3557      delete [] object_;
3558      numberObjects_=rhs.numberObjects_;
3559      if (numberObjects_) {
3560        object_ = new OsiObject * [numberObjects_];
3561        int i;
3562        for (i=0;i<numberObjects_;i++) 
3563          object_[i]=(rhs.object_[i])->clone();
3564      } else {
3565        object_=NULL;
3566    }
3567    } else {
3568      // assume will be redone
3569      numberObjects_=0;
3570      object_=NULL;
3571    }
3572    delete [] originalColumns_;
3573    if (rhs.originalColumns_) {
3574      int numberColumns = rhs.getNumCols();
3575      originalColumns_= new int [numberColumns];
3576      memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
3577    } else {
3578      originalColumns_=NULL;
3579    }
3580    nodeCompare_=rhs.nodeCompare_->clone();
3581    problemFeasibility_=rhs.problemFeasibility_->clone();
3582    delete tree_;
3583    tree_= rhs.tree_->clone();
3584    if (rhs.branchingMethod_)
3585      branchingMethod_=rhs.branchingMethod_->clone();
3586    else
3587      branchingMethod_=NULL;
3588    if (rhs.cutModifier_)
3589      cutModifier_=rhs.cutModifier_->clone();
3590    else
3591      cutModifier_=NULL;
3592    delete strategy_;
3593    if (rhs.strategy_)
3594      strategy_=rhs.strategy_->clone();
3595    else
3596      strategy_=NULL;
3597    parentModel_=rhs.parentModel_;
3598    appData_=rhs.appData_;
3599
3600    delete [] integerVariable_;
3601    numberIntegers_=rhs.numberIntegers_;
3602    if (numberIntegers_) {
3603      integerVariable_ = new int [numberIntegers_];
3604      memcpy(integerVariable_,rhs.integerVariable_,
3605             numberIntegers_*sizeof(int));
3606      integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
3607    } else {
3608      integerVariable_ = NULL;
3609      integerInfo_=NULL;
3610    }
3611    if (rhs.hotstartSolution_) {
3612      int numberColumns = solver_->getNumCols();
3613      hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
3614      hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
3615    } else {
3616      hotstartSolution_ = NULL;
3617      hotstartPriorities_ =NULL;
3618    }
3619    numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
3620    maximumNumberCuts_=rhs.maximumNumberCuts_;
3621    phase_ = rhs.phase_;
3622    currentNumberCuts_=rhs.currentNumberCuts_;
3623    maximumDepth_= rhs.maximumDepth_;
3624    delete [] addedCuts_;
3625    delete [] walkback_;
3626    // These are only used as temporary arrays so need not be filled
3627    if (maximumNumberCuts_) {
3628      addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
3629    } else {
3630      addedCuts_ = NULL;
3631    }
3632    nextRowCut_ = NULL;
3633    currentNode_ = NULL;
3634    if (maximumDepth_)
3635      walkback_ = new CbcNodeInfo * [maximumDepth_];
3636    else
3637      walkback_ = NULL;
3638    synchronizeModel();
3639    cbcColLower_ = NULL;
3640    cbcColUpper_ = NULL;
3641    cbcRowLower_ = NULL;
3642    cbcRowUpper_ = NULL;
3643    cbcColSolution_ = NULL;
3644    cbcRowPrice_ = NULL;
3645    cbcReducedCost_ = NULL;
3646    cbcRowActivity_ = NULL;
3647  }
3648  return *this;
3649}
3650 
3651// Destructor
3652CbcModel::~CbcModel ()
3653{
3654  if (defaultHandler_) {
3655    delete handler_;
3656    handler_ = NULL;
3657  }
3658  delete tree_;
3659  tree_=NULL;
3660  if (ourSolver_) delete solver_;
3661  gutsOfDestructor();
3662  delete eventHandler_ ;
3663  eventHandler_ = NULL ;
3664}
3665// Clears out as much as possible (except solver)
3666void 
3667CbcModel::gutsOfDestructor()
3668{
3669  delete referenceSolver_;
3670  referenceSolver_=NULL;
3671  int i;
3672  for (i=0;i<numberCutGenerators_;i++) {
3673    delete generator_[i];
3674    delete virginGenerator_[i];
3675  }
3676  delete [] generator_;
3677  delete [] virginGenerator_;
3678  generator_=NULL;
3679  virginGenerator_=NULL;
3680  for (i=0;i<numberHeuristics_;i++)
3681    delete heuristic_[i];
3682  delete [] heuristic_;
3683  heuristic_=NULL;
3684  delete nodeCompare_;
3685  nodeCompare_=NULL;
3686  delete problemFeasibility_;
3687  problemFeasibility_=NULL;
3688  delete [] originalColumns_;
3689  originalColumns_=NULL;
3690  delete strategy_;
3691#if NEW_UPDATE_OBJECT>1
3692  delete [] updateItems_;
3693  updateItems_=NULL;
3694  numberUpdateItems_=0;
3695  maximumNumberUpdateItems_=0;
3696#endif
3697  gutsOfDestructor2();
3698}
3699// Clears out enough to reset CbcModel
3700void 
3701CbcModel::gutsOfDestructor2()
3702{
3703  delete [] integerInfo_;
3704  integerInfo_=NULL;
3705  delete [] integerVariable_;
3706  integerVariable_=NULL;
3707  int i;
3708  if (ownObjects_) {
3709    for (i=0;i<numberObjects_;i++)
3710      delete object_[i];
3711    delete [] object_;
3712  }
3713  ownObjects_=true;
3714  object_=NULL;
3715  numberIntegers_=0;
3716  numberObjects_=0;
3717  delete emptyWarmStart_ ;
3718  emptyWarmStart_ =NULL;
3719  delete continuousSolver_;
3720  continuousSolver_=NULL;
3721  delete [] bestSolution_;
3722  bestSolution_=NULL;
3723  delete [] currentSolution_;
3724  currentSolution_=NULL;
3725  delete [] continuousSolution_;
3726  continuousSolution_=NULL;
3727  solverCharacteristics_=NULL;
3728  delete [] usedInSolution_;
3729  usedInSolution_ = NULL;
3730  testSolution_=NULL;
3731  lastHeuristic_ = NULL;
3732  delete [] addedCuts_;
3733  addedCuts_=NULL;
3734  nextRowCut_ = NULL;
3735  currentNode_ = NULL;
3736  delete [] walkback_;
3737  walkback_=NULL;
3738  delete [] whichGenerator_;
3739  whichGenerator_ = NULL;
3740  for (i=0;i<maximumStatistics_;i++)
3741    delete statistics_[i];
3742  delete [] statistics_;
3743  statistics_=NULL;
3744  maximumDepthActual_ = 0;
3745  numberDJFixed_ =0.0;
3746  delete probingInfo_;
3747  probingInfo_ = NULL;
3748  maximumStatistics_=0;
3749  delete [] analyzeResults_;
3750  analyzeResults_=NULL;
3751  // Below here is whatever consensus is
3752  ourSolver_=true;
3753  bestObjective_=COIN_DBL_MAX;
3754  bestPossibleObjective_=COIN_DBL_MAX;
3755  sumChangeObjective1_=0.0;
3756  sumChangeObjective2_=0.0;
3757  numberSolutions_=0;
3758  stateOfSearch_=0;
3759  delete [] hotstartSolution_;
3760  hotstartSolution_=NULL;
3761  delete [] hotstartPriorities_;
3762  hotstartPriorities_=NULL;
3763  numberHeuristicSolutions_=0;
3764  numberNodes_=0;
3765  numberNodes2_=0;
3766  numberIterations_=0;
3767  status_=-1;
3768  secondaryStatus_=-1;
3769  maximumNumberCuts_=0;
3770  phase_=0;
3771  currentNumberCuts_=0;
3772  maximumDepth_=0;
3773  nextRowCut_=NULL;
3774  currentNode_=NULL;
3775  delete branchingMethod_;
3776  branchingMethod_=NULL;
3777  delete cutModifier_;
3778  cutModifier_=NULL;
3779  // clear out tree
3780  if (tree_&&tree_->size())
3781    tree_->cleanTree(this, -1.0e100,bestPossibleObjective_) ;
3782  subTreeModel_=NULL;
3783  numberStoppedSubTrees_=0;
3784  numberInfeasibleNodes_=0;
3785  numberGlobalViolations_=0;
3786  continuousObjective_=0.0;
3787  originalContinuousObjective_=0.0;
3788  continuousInfeasibilities_=0;
3789  numberFixedAtRoot_=0;
3790  numberFixedNow_=0;
3791  stoppedOnGap_=false;
3792  eventHappened_=false;
3793  numberLongStrong_=0;
3794  numberOldActiveCuts_=0;
3795  numberNewCuts_=0;
3796  searchStrategy_=-1;
3797  numberStrongIterations_=0;
3798  // Parameters which need to be reset
3799  dblParam_[CbcCutoffIncrement] = 1e-5;
3800  dblParam_[CbcCurrentCutoff] = 1.0e100;
3801  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
3802  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
3803}
3804// Save a copy of the current solver so can be reset to
3805void 
3806CbcModel::saveReferenceSolver()
3807{
3808  delete referenceSolver_;
3809  referenceSolver_= solver_->clone();
3810}
3811
3812// Uses a copy of reference solver to be current solver
3813void 
3814CbcModel::resetToReferenceSolver()
3815{
3816  delete solver_;
3817  solver_ = referenceSolver_->clone();
3818  // clear many things
3819  gutsOfDestructor2();
3820  // Reset cutoff
3821  // Solvers know about direction
3822  double direction = solver_->getObjSense();
3823  double value;
3824  solver_->getDblParam(OsiDualObjectiveLimit,value); 
3825  setCutoff(value*direction);
3826}
3827
3828// Are there a numerical difficulties?
3829bool 
3830CbcModel::isAbandoned() const
3831{
3832  return status_ == 2;
3833}
3834// Is optimality proven?
3835bool 
3836CbcModel::isProvenOptimal() const
3837{
3838  if (!status_ && bestObjective_<1.0e30)
3839    return true;
3840  else
3841    return false;
3842}
3843// Is  infeasiblity proven (or none better than cutoff)?
3844bool 
3845CbcModel::isProvenInfeasible() const
3846{
3847  if (!status_ && bestObjective_>=1.0e30)
3848    return true;
3849  else
3850    return false;
3851}
3852// Was continuous solution unbounded
3853bool 
3854CbcModel::isContinuousUnbounded() const
3855{
3856  if (!status_ && secondaryStatus_==7)
3857    return true;
3858  else
3859    return false;
3860}
3861// Was continuous solution unbounded
3862bool 
3863CbcModel::isProvenDualInfeasible() const
3864{
3865  if (!status_ && secondaryStatus_==7)
3866    return true;
3867  else
3868    return false;
3869}
3870// Node limit reached?
3871bool 
3872CbcModel::isNodeLimitReached() const
3873{
3874  return numberNodes_ >= intParam_[CbcMaxNumNode];
3875}
3876// Time limit reached?
3877bool 
3878CbcModel::isSecondsLimitReached() const
3879{
3880  if (status_==1&&secondaryStatus_==4)
3881    return true;
3882  else
3883    return false;
3884}
3885// Solution limit reached?
3886bool 
3887CbcModel::isSolutionLimitReached() const
3888{
3889  return numberSolutions_ >= intParam_[CbcMaxNumSol];
3890}
3891// Set language
3892void 
3893CbcModel::newLanguage(CoinMessages::Language language)
3894{
3895  messages_ = CbcMessage(language);
3896}
3897void 
3898CbcModel::setNumberStrong(int number)
3899{
3900  if (number<0)
3901    numberStrong_=0;
3902   else
3903    numberStrong_=number;
3904}
3905void 
3906CbcModel::setNumberBeforeTrust(int number)
3907{
3908  if (number<-1) {
3909    numberBeforeTrust_=0;
3910  } else {
3911    numberBeforeTrust_=number;
3912    //numberStrong_ = CoinMax(numberStrong_,1);
3913  }
3914}
3915void 
3916CbcModel::setNumberPenalties(int number)
3917{
3918  if (number<=0) {
3919    numberPenalties_=0;
3920  } else {
3921    numberPenalties_=number;
3922  }
3923}
3924void 
3925CbcModel::setPenaltyScaleFactor(double value)
3926{
3927  if (value<=0) {
3928    penaltyScaleFactor_=3.0;
3929  } else {
3930    penaltyScaleFactor_=value;
3931  }
3932}
3933void 
3934CbcModel::setHowOftenGlobalScan(int number)
3935{
3936  if (number<-1)
3937    howOftenGlobalScan_=0;
3938   else
3939    howOftenGlobalScan_=number;
3940}
3941
3942// Add one generator
3943void 
3944CbcModel::addCutGenerator(CglCutGenerator * generator,
3945                          int howOften, const char * name,
3946                          bool normal, bool atSolution,
3947                          bool whenInfeasible,int howOftenInSub,
3948                          int whatDepth, int whatDepthInSub)
3949{
3950  CbcCutGenerator ** temp = generator_;
3951  generator_ = new CbcCutGenerator * [numberCutGenerators_+1];
3952  memcpy(generator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
3953  delete[] temp ;
3954  generator_[numberCutGenerators_]= 
3955    new CbcCutGenerator(this,generator, howOften, name,
3956                        normal,atSolution,whenInfeasible,howOftenInSub,
3957                        whatDepth, whatDepthInSub);
3958  // and before any cahnges
3959  temp = virginGenerator_;
3960  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_+1];
3961  memcpy(virginGenerator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
3962  delete[] temp ;
3963  virginGenerator_[numberCutGenerators_++]= 
3964    new CbcCutGenerator(this,generator, howOften, name,
3965                        normal,atSolution,whenInfeasible,howOftenInSub,
3966                        whatDepth, whatDepthInSub);
3967                                                         
3968}
3969// Add one heuristic
3970void 
3971CbcModel::addHeuristic(CbcHeuristic * generator, const char *name)
3972{
3973  CbcHeuristic ** temp = heuristic_;
3974  heuristic_ = new CbcHeuristic * [numberHeuristics_+1];
3975  memcpy(heuristic_,temp,numberHeuristics_*sizeof(CbcHeuristic *));
3976  delete [] temp;
3977  heuristic_[numberHeuristics_]=generator->clone();
3978  if (name)
3979  { heuristic_[numberHeuristics_]->setHeuristicName(name) ; }
3980  numberHeuristics_++ ;
3981}
3982
3983/*
3984  The last subproblem handled by the solver is not necessarily related to the
3985  one being recreated, so the first action is to remove all cuts from the
3986  constraint system.  Next, traverse the tree from node to the root to
3987  determine the basis size required for this subproblem and create an empty
3988  basis with the right capacity.  Finally, traverse the tree from root to
3989  node, adjusting bounds in the constraint system, adjusting the basis, and
3990  collecting the cuts that must be added to the constraint system.
3991  applyToModel does the heavy lifting.
3992
3993  addCuts1 is used in contexts where all that's desired is the list of cuts:
3994  the node is already fathomed, and we're collecting cuts so that we can
3995  adjust reference counts as we prune nodes. Arguably the two functions
3996  should be separated. The culprit is applyToModel, which performs cut
3997  collection and model adjustment.
3998
3999  Certainly in the contexts where all we need is a list of cuts, there's no
4000  point in passing in a valid basis --- an empty basis will do just fine.
4001*/
4002void CbcModel::addCuts1 (CbcNode * node, CoinWarmStartBasis *&lastws)
4003{ int i;
4004  int nNode=0;
4005  int numberColumns = getNumCols();
4006  CbcNodeInfo * nodeInfo = node->nodeInfo();
4007
4008/*
4009  Remove all cuts from the constraint system.
4010  (original comment includes ``see note below for later efficiency'', but
4011  the reference isn't clear to me).
4012*/
4013  int currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_;
4014  int *which = new int[currentNumberCuts];
4015  for (i = 0 ; i < currentNumberCuts ; i++)
4016    which[i] = i+numberRowsAtContinuous_;
4017  solver_->deleteRows(currentNumberCuts,which);
4018  delete [] which;
4019/*
4020  Accumulate the path from node to the root in walkback_, and accumulate a
4021  cut count in currentNumberCuts.
4022
4023  original comment: when working then just unwind until where new node joins
4024  old node (for cuts?)
4025*/
4026  currentNumberCuts = 0;
4027  while (nodeInfo) {
4028    //printf("nNode = %d, nodeInfo = %x\n",nNode,nodeInfo);
4029    walkback_[nNode++]=nodeInfo;
4030    currentNumberCuts += nodeInfo->numberCuts() ;
4031    nodeInfo = nodeInfo->parent() ;
4032    if (nNode==maximumDepth_) {
4033      maximumDepth_ *= 2;
4034      CbcNodeInfo ** temp = new CbcNodeInfo * [maximumDepth_];
4035      for (i=0;i<nNode;i++) 
4036        temp[i] = walkback_[i];
4037      delete [] walkback_;
4038      walkback_ = temp;
4039    }
4040  }
4041/*
4042  Create an empty basis with sufficient capacity for the constraint system
4043  we'll construct: original system plus cuts. Make sure we have capacity to
4044  record those cuts in addedCuts_.
4045
4046  The method of adjusting the basis at a FullNodeInfo object (the root, for
4047  example) is to use a copy constructor to duplicate the basis held in the
4048  nodeInfo, then resize it and return the new basis object. Guaranteed,
4049  lastws will point to a different basis when it returns. We pass in a basis
4050  because we need the parameter to return the allocated basis, and it's an
4051  easy way to pass in the size. But we take a hit for memory allocation.
4052*/
4053  currentNumberCuts_=currentNumberCuts;
4054  if (currentNumberCuts > maximumNumberCuts_) {
4055    maximumNumberCuts_ = currentNumberCuts;
4056    delete [] addedCuts_;
4057    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4058  }
4059  lastws->setSize(numberColumns,numberRowsAtContinuous_+currentNumberCuts);
4060/*
4061  This last bit of code traverses the path collected in walkback_ from the
4062  root back to node. At the end of the loop,
4063   * lastws will be an appropriate basis for node;
4064   * variable bounds in the constraint system will be set to be correct for
4065     node; and
4066   * addedCuts_ will be set to a list of cuts that need to be added to the
4067     constraint system at node.
4068  applyToModel does all the heavy lifting.
4069*/
4070  currentNumberCuts=0;
4071  while (nNode) {
4072    --nNode;
4073    walkback_[nNode]->applyToModel(this,lastws,addedCuts_,currentNumberCuts);
4074  }
4075}
4076
4077/*
4078  adjustCuts might be a better name: If the node is feasible, we sift through
4079  the cuts collected by addCuts1, add the ones that are tight and omit the
4080  ones that are loose. If the node is infeasible, we just adjust the
4081  reference counts to reflect that we're about to prune this node and its
4082  descendants.
4083*/
4084int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws,bool canFix)
4085{
4086/*
4087  addCuts1 performs step 1 of restoring the subproblem at this node; see the
4088  comments there.
4089*/
4090  addCuts1(node,lastws);
4091  int i;
4092  int numberColumns = getNumCols();
4093  CbcNodeInfo * nodeInfo = node->nodeInfo();
4094  double cutoff = getCutoff() ;
4095  int currentNumberCuts=currentNumberCuts_;
4096  if (canFix) {
4097    bool feasible=true;
4098    const double *lower = solver_->getColLower() ;
4099    const double *upper = solver_->getColUpper() ;
4100    double * newLower = analyzeResults_;
4101    double * objLower = newLower+numberIntegers_;
4102    double * newUpper = objLower+numberIntegers_;
4103    double * objUpper = newUpper+numberIntegers_;
4104    int n=0;
4105    for (i=0;i<numberIntegers_;i++) {
4106      int iColumn = integerVariable_[i];
4107      bool changed=false;
4108      double lo = 0.0;
4109      double up = 0.0;
4110      if (objLower[i]>cutoff) {
4111        lo = lower[iColumn];
4112        up = upper[iColumn];
4113        if (lo<newLower[i]) {
4114          lo = newLower[i];
4115          solver_->setColLower(iColumn,lo);
4116          changed=true;
4117          n++;
4118        }
4119        if (objUpper[i]>cutoff) {
4120          if (up>newUpper[i]) {
4121            up = newUpper[i];
4122            solver_->setColUpper(iColumn,up);
4123            changed=true;
4124            n++;
4125          }
4126        }
4127      } else if (objUpper[i]>cutoff) {
4128        lo = lower[iColumn];
4129        up = upper[iColumn];
4130        if (up>newUpper[i]) {
4131          up = newUpper[i];
4132          solver_->setColUpper(iColumn,up);
4133          changed=true;
4134          n++;
4135        }
4136      }
4137      if (changed&&lo>up) {
4138        feasible=false;
4139        break;
4140      }
4141    }
4142    if (!feasible) {
4143      printf("analysis says node infeas\n");
4144      cutoff=-COIN_DBL_MAX;
4145    }
4146  }
4147/*
4148  If the node can't be fathomed by bound, reinstall tight cuts in the
4149  constraint system. Even if there are no cuts, we'll want to set the
4150  reconstructed basis in the solver.
4151*/
4152  if (node->objectiveValue() < cutoff||numberThreads_)
4153  { 
4154#   ifdef CBC_CHECK_BASIS
4155    printf("addCuts: expanded basis; rows %d+%d\n",
4156           numberRowsAtContinuous_,currentNumberCuts);
4157    lastws->print();
4158#   endif
4159/*
4160  Adjust the basis and constraint system so that we retain only active cuts.
4161  There are three steps:
4162    1) Scan the basis. Sort the cuts into effective cuts to be kept and
4163       loose cuts to be dropped.
4164    2) Drop the loose cuts and resize the basis to fit.
4165    3) Install the tight cuts in the constraint system (applyRowCuts) and
4166       and install the basis (setWarmStart).
4167  Use of compressRows conveys we're compressing the basis and not just
4168  tweaking the artificialStatus_ array.
4169*/
4170    if (currentNumberCuts > 0) {
4171      int numberToAdd = 0;
4172      const OsiRowCut **addCuts;
4173      int numberToDrop = 0 ;
4174      int *cutsToDrop ;
4175      addCuts = new const OsiRowCut* [currentNumberCuts];
4176      cutsToDrop = new int[currentNumberCuts] ;
4177      assert (currentNumberCuts+numberRowsAtContinuous_<=lastws->getNumArtificial());
4178      for (i=0;i<currentNumberCuts;i++) {
4179        CoinWarmStartBasis::Status status = 
4180          lastws->getArtifStatus(i+numberRowsAtContinuous_);
4181        if (addedCuts_[i] &&
4182            (status != CoinWarmStartBasis::basic ||
4183             addedCuts_[i]->effectiveness()==COIN_DBL_MAX)) {
4184#         ifdef CHECK_CUT_COUNTS
4185          printf("Using cut %d %x as row %d\n",i,addedCuts_[i],
4186                 numberRowsAtContinuous_+numberToAdd);
4187#         endif
4188          addCuts[numberToAdd++] = new OsiRowCut(*addedCuts_[i]);
4189        } else {
4190#         ifdef CHECK_CUT_COUNTS
4191          printf("Dropping cut %d %x\n",i,addedCuts_[i]);
4192#         endif
4193          addedCuts_[i]=NULL;
4194          cutsToDrop[numberToDrop++] = numberRowsAtContinuous_+i ;
4195        }
4196      }
4197      int numberRowsNow=numberRowsAtContinuous_+numberToAdd;
4198      lastws->compressRows(numberToDrop,cutsToDrop) ;
4199      lastws->resize(numberRowsNow,numberColumns);
4200      solver_->applyRowCuts(numberToAdd,addCuts);
4201#     ifdef CBC_CHECK_BASIS
4202      printf("addCuts: stripped basis; rows %d + %d\n",
4203             numberRowsAtContinuous_,numberToAdd);
4204      lastws->print();
4205#     endif
4206      for (i=0;i<numberToAdd;i++)
4207        delete addCuts[i];
4208      delete [] addCuts;
4209      delete [] cutsToDrop ;
4210    }
4211/*
4212  Set the basis in the solver.
4213*/
4214    solver_->setWarmStart(lastws);
4215#if 0
4216    if ((numberNodes_%printFrequency_)==0) {
4217      printf("Objective %g, depth %d, unsatisfied %d\n",
4218             node->objectiveValue(),
4219             node->depth(),node->numberUnsatisfied());
4220    }
4221#endif
4222/*
4223  Clean up and we're out of here.
4224*/
4225    numberNodes_++;
4226    return 0;
4227  } 
4228/*
4229  This node has been fathomed by bound as we try to revive it out of the live
4230  set. Adjust the cut reference counts to reflect that we no longer need to
4231  explore the remaining branch arms, hence they will no longer reference any
4232  cuts. Cuts whose reference count falls to zero are deleted. 
4233*/
4234  else
4235  { int i;
4236    if (currentNumberCuts) {
4237      lockThread();
4238      int numberLeft = nodeInfo->numberBranchesLeft();
4239      for (i = 0 ; i < currentNumberCuts ; i++)
4240        { if (addedCuts_[i])
4241          { if (!addedCuts_[i]->decrement(numberLeft))
4242            { delete addedCuts_[i];
4243            addedCuts_[i] = NULL; } } }
4244      unlockThread();
4245    }
4246    return 1 ; }
4247}
4248
4249
4250/*
4251  Perform reduced cost fixing on integer variables.
4252
4253  The variables in question are already nonbasic at bound. We're just nailing
4254  down the current situation.
4255*/
4256
4257int CbcModel::reducedCostFix ()
4258
4259{
4260  if(!solverCharacteristics_->reducedCostsAccurate())
4261    return 0; //NLP
4262  double cutoff = getCutoff() ;
4263  double direction = solver_->getObjSense() ;
4264  double gap = cutoff - solver_->getObjValue()*direction ;
4265  double tolerance;
4266  solver_->getDblParam(OsiDualTolerance,tolerance) ;
4267  if (gap<=0.0)
4268    return 0;
4269  gap += 100.0*tolerance;
4270  double integerTolerance = getDblParam(CbcIntegerTolerance) ;
4271
4272  const double *lower = solver_->getColLower() ;
4273  const double *upper = solver_->getColUpper() ;
4274  const double *solution = solver_->getColSolution() ;
4275  const double *reducedCost = solver_->getReducedCost() ;
4276
4277  int numberFixed = 0 ;
4278
4279# ifdef COIN_HAS_CLP
4280  OsiClpSolverInterface * clpSolver
4281    = dynamic_cast<OsiClpSolverInterface *> (solver_);
4282  ClpSimplex * clpSimplex=NULL;
4283  if (clpSolver) 
4284    clpSimplex = clpSolver->getModelPtr();
4285# endif
4286  for (int i = 0 ; i < numberIntegers_ ; i++)
4287  { int iColumn = integerVariable_[i] ;
4288    double djValue = direction*reducedCost[iColumn] ;
4289    if (upper[iColumn]-lower[iColumn] > integerTolerance)
4290    { if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap)
4291      { solver_->setColUpper(iColumn,lower[iColumn]) ;
4292#ifdef COIN_HAS_CLP
4293      // may just have been fixed before
4294      if (clpSimplex)
4295        assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
4296               clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4297#endif
4298      numberFixed++ ; }
4299      else
4300      if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap)
4301      { solver_->setColLower(iColumn,upper[iColumn]) ;
4302#ifdef COIN_HAS_CLP
4303      // may just have been fixed before
4304      if (clpSimplex)
4305        assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
4306               clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
4307#endif
4308      numberFixed++ ; } } }
4309  numberDJFixed_ += numberFixed;
4310 
4311  return numberFixed; }
4312
4313// Collect coding to replace whichGenerator
4314void
4315CbcModel::resizeWhichGenerator(int numberNow, int numberAfter)
4316{
4317  if (numberAfter > maximumWhich_) {
4318    maximumWhich_ = CoinMax(maximumWhich_*2+100,numberAfter) ;
4319    int * temp = new int[2*maximumWhich_] ;
4320    memcpy(temp,whichGenerator_,numberNow*sizeof(int)) ;
4321    delete [] whichGenerator_ ;
4322    whichGenerator_ = temp ;
4323    memset(whichGenerator_+numberNow,0,(maximumWhich_-numberNow)*sizeof(int));
4324  }
4325}
4326
4327/** Solve the model using cuts
4328
4329  This version takes off redundant cuts from node.
4330  Returns true if feasible.
4331
4332  \todo
4333  Why do I need to resolve the problem? What has been done between the last
4334  relaxation and calling solveWithCuts?
4335
4336  If numberTries == 0 then user did not want any cuts.
4337*/
4338
4339bool 
4340CbcModel::solveWithCuts (OsiCuts &cuts, int numberTries, CbcNode *node)
4341/*
4342  Parameters:
4343    numberTries: (i) the maximum number of iterations for this round of cut
4344                     generation; if negative then we don't mind if drop is tiny.
4345   
4346    cuts:       (o) all cuts generated in this round of cut generation
4347
4348    node: (i)     So we can update dynamic pseudo costs
4349*/
4350                       
4351
4352{
4353# ifdef COIN_HAS_CLP
4354  OsiClpSolverInterface * clpSolver
4355    = dynamic_cast<OsiClpSolverInterface *> (solver_);
4356  int saveClpOptions=0;
4357  if (clpSolver) 
4358    saveClpOptions = clpSolver->specialOptions();
4359# endif
4360#ifdef CBC_THREAD
4361  CbcModel ** threadModel = NULL;
4362  pthread_t * threadId = NULL;
4363  pthread_cond_t condition_main;
4364  pthread_mutex_t condition_mutex;
4365  pthread_mutex_t * mutex2 = NULL;
4366  pthread_cond_t * condition2 = NULL;
4367  threadStruct * threadInfo = NULL;
4368  void * saveMutex = NULL;
4369  if (numberThreads_&&(threadMode_&2)!=0&&!numberNodes_) {
4370    threadId = new pthread_t [numberThreads_];
4371    pthread_cond_init(&condition_main,NULL);
4372    pthread_mutex_init(&condition_mutex,NULL);
4373    threadModel = new CbcModel * [numberThreads_];
4374    threadInfo = new threadStruct [numberThreads_+1];
4375    mutex2 = new pthread_mutex_t [numberThreads_];
4376    condition2 = new pthread_cond_t [numberThreads_];
4377    saveMutex = mutex_;
4378    for (int i=0;i<numberThreads_;i++) {
4379      pthread_mutex_init(mutex2+i,NULL);
4380      pthread_cond_init(condition2+i,NULL);
4381      threadId[i]=0;
4382      threadModel[i]=new CbcModel;
4383      threadModel[i]->generator_ = new CbcCutGenerator * [1];
4384      delete threadModel[i]->solver_;
4385      threadModel[i]->solver_=NULL;
4386      threadModel[i]->numberThreads_=numberThreads_;
4387      mutex_ = (void *) (threadInfo+i);
4388      threadInfo[i].thisModel=(CbcModel *) threadModel[i];
4389      threadInfo[i].baseModel=this;
4390      threadInfo[i].threadIdOfBase=pthread_self();
4391      threadInfo[i].mutex2=mutex2+i;
4392      threadInfo[i].condition2=condition2+i;
4393      threadInfo[i].returnCode=-1;
4394      pthread_create(threadId+i,NULL,doCutsThread,threadInfo+i);
4395    }
4396    // Do a partial one for base model
4397    threadInfo[numberThreads_].baseModel=this;
4398    mutex_ = (void *) (threadInfo+numberThreads_);
4399    threadInfo[numberThreads_].condition2=&condition_main;
4400    threadInfo[numberThreads_].mutex2=&condition_mutex;
4401  }
4402#endif
4403  bool feasible = true ;
4404  int lastNumberCuts = 0 ;
4405  double lastObjective = -1.0e100 ;
4406  int violated = 0 ;
4407  int numberRowsAtStart = solver_->getNumRows() ;
4408  //printf("solver had %d rows\n",numberRowsAtStart);
4409  int numberColumns = solver_->getNumCols() ;
4410  CoinBigIndex numberElementsAtStart = solver_->getNumElements();
4411
4412  numberOldActiveCuts_ = numberRowsAtStart-numberRowsAtContinuous_ ;
4413  numberNewCuts_ = 0 ;
4414
4415  bool onOptimalPath = false ;
4416  const OsiRowCutDebugger *debugger = NULL;
4417  if ((specialOptions_&1)!=0) {
4418    /*
4419      See OsiRowCutDebugger for details. In a nutshell, make sure that current
4420      variable values do not conflict with a known optimal solution. (Obviously
4421      this can be fooled when there are multiple solutions.)
4422    */
4423    debugger = solver_->getRowCutDebugger() ;
4424    if (debugger) 
4425      onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
4426  }
4427  OsiCuts slackCuts;
4428/*
4429  Resolve the problem. If we've lost feasibility, might as well bail out right
4430  after the debug stuff. The resolve will also refresh cached copies of the
4431  solver solution (cbcColLower_, ...) held by CbcModel.
4432*/
4433  double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
4434  if (node)
4435    objectiveValue= node->objectiveValue();
4436  int returnCode = resolve(node ? node->nodeInfo() : NULL,1);
4437#ifdef COIN_DEVELOP
4438  //if (!solver_->getIterationCount()&&solver_->isProvenOptimal())
4439  //printf("zero iterations on first solve of branch\n");
4440#endif
4441  if (node&&node->nodeInfo()&&!node->nodeInfo()->numberBranchesLeft())
4442    node->nodeInfo()->allBranchesGone(); // can clean up
4443  feasible = returnCode  != 0 ;
4444  if (returnCode<0)
4445    numberTries=0;
4446  if (problemFeasibility_->feasible(this,0)<0) {
4447    feasible=false; // pretend infeasible
4448  }
4449 
4450#if NEW_UPDATE_OBJECT==0
4451  // Update branching information if wanted
4452  if(node &&branchingMethod_)
4453    branchingMethod_->updateInformation(solver_,node);
4454#elif NEW_UPDATE_OBJECT<2
4455  // Update branching information if wanted
4456  if(node &&branchingMethod_) {
4457    OsiBranchingObject * bobj = node->modifiableBranchingObject();
4458    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
4459    if (cbcobj) {
4460      CbcObject * object = cbcobj->object();
4461      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
4462      object->updateInformation(update);
4463    } else {
4464      branchingMethod_->updateInformation(solver_,node);
4465    }
4466  }
4467#else
4468  // Update branching information if wanted
4469  if(node &&branchingMethod_) {
4470    OsiBranchingObject * bobj = node->modifiableBranchingObject();
4471    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
4472    if (cbcobj) {
4473      CbcObject * object = cbcobj->object();
4474      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
4475      // have to compute object number as not saved
4476      CbcSimpleInteger * simpleObject =
4477          dynamic_cast <CbcSimpleInteger *>(object) ;
4478      int iObject;
4479      int iColumn = simpleObject->columnNumber();
4480      for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
4481        simpleObject =
4482          dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ;
4483        if (simpleObject->columnNumber()==iColumn)
4484          break;
4485      }
4486      assert (iObject<numberObjects_);
4487      update.objectNumber_ = iObject;
4488      addUpdateInformation(update);
4489    } else {
4490      branchingMethod_->updateInformation(solver_,node);
4491    }
4492  }
4493#endif
4494
4495#ifdef CBC_DEBUG
4496  if (feasible)
4497  { printf("Obj value %g (%s) %d rows\n",solver_->getObjValue(),
4498           (solver_->isProvenOptimal())?"proven":"unproven",
4499           solver_->getNumRows()) ; }
4500 
4501  else
4502  { printf("Infeasible %d rows\n",solver_->getNumRows()) ; }
4503#endif
4504  if ((specialOptions_&1)!=0) {
4505/*
4506  If the RowCutDebugger said we were compatible with the optimal solution,
4507  and now we're suddenly infeasible, we might be confused. Then again, we
4508  may have fathomed by bound, heading for a rediscovery of an optimal solution.
4509*/
4510    if (onOptimalPath && !solver_->isDualObjectiveLimitReached()) {
4511      if (!feasible) {
4512        solver_->writeMps("infeas");
4513        CoinWarmStartBasis *slack =
4514          dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
4515        solver_->setWarmStart(slack);
4516        delete slack ;
4517        solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
4518        solver_->initialSolve();
4519      }
4520      assert(feasible) ;
4521    }
4522  }
4523
4524  if (!feasible) {
4525    numberInfeasibleNodes_++;
4526# ifdef COIN_HAS_CLP
4527    if (clpSolver) 
4528    clpSolver->setSpecialOptions(saveClpOptions);
4529# endif
4530    return (false) ;
4531  }
4532  sumChangeObjective1_ += solver_->getObjValue()*solver_->getObjSense()
4533    - objectiveValue ;
4534  if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] )
4535    numberTries=0; // exit
4536  //if ((numberNodes_%100)==0)
4537  //printf("XXa sum obj changed by %g\n",sumChangeObjective1_);
4538  objectiveValue = solver_->getObjValue()*solver_->getObjSense();
4539  // Return at once if numberTries zero
4540  if (!numberTries) {
4541    cuts=OsiCuts();
4542    numberNewCuts_=0;
4543# ifdef COIN_HAS_CLP
4544    if (clpSolver) 
4545    clpSolver->setSpecialOptions(saveClpOptions);
4546# endif
4547    return true;
4548  }
4549/*
4550  Do reduced cost fixing.
4551*/
4552  reducedCostFix() ;
4553/*
4554  Set up for at most numberTries rounds of cut generation. If numberTries is
4555  negative, we'll ignore the minimumDrop_ cutoff and keep generating cuts for
4556  the specified number of rounds.
4557*/
4558  double minimumDrop = minimumDrop_ ;
4559  if (numberTries<0)
4560  { numberTries = -numberTries ;
4561    minimumDrop = -1.0 ; }
4562/*
4563  Is it time to scan the cuts in order to remove redundant cuts? If so, set
4564  up to do it.
4565*/
4566# define SCANCUTS 100 
4567  int *countColumnCuts = NULL ;
4568  // Always accumulate row cut counts
4569  int * countRowCuts =new int[numberCutGenerators_+numberHeuristics_] ;
4570  memset(countRowCuts,0,
4571         (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ;
4572  bool fullScan = false ;
4573  if ((numberNodes_%SCANCUTS) == 0)
4574  { fullScan = true ;
4575    countColumnCuts = new int[numberCutGenerators_+numberHeuristics_] ;
4576    memset(countColumnCuts,0,
4577           (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ; }
4578
4579  double direction = solver_->getObjSense() ;
4580  double startObjective = solver_->getObjValue()*direction ;
4581
4582  currentPassNumber_ = 0 ;
4583  double primalTolerance = 1.0e-7 ;
4584  // We may need to keep going on
4585  bool keepGoing=false;
4586/*
4587  Begin cut generation loop. Cuts generated during each iteration are
4588  collected in theseCuts. The loop can be divided into four phases:
4589   1) Prep: Fix variables using reduced cost. In the first iteration only,
4590      consider scanning globalCuts_ and activating any applicable cuts.
4591   2) Cut Generation: Call each generator and heuristic registered in the
4592      generator_ and heuristic_ arrays. Newly generated global cuts are
4593      copied to globalCuts_ at this time.
4594   3) Cut Installation and Reoptimisation: Install column and row cuts in
4595      the solver. Copy row cuts to cuts (parameter). Reoptimise.
4596   4) Cut Purging: takeOffCuts() removes inactive cuts from the solver, and
4597      does the necessary bookkeeping in the model.
4598*/
4599  do
4600  { currentPassNumber_++ ;
4601    numberTries-- ;
4602    if (numberTries<0&&keepGoing) {
4603      // switch off all normal ones
4604      for (int i = 0;i<numberCutGenerators_;i++) {
4605        if (!generator_[i]->mustCallAgain())
4606          generator_[i]->setSwitchedOff(true);
4607      }
4608    }
4609    keepGoing=false;
4610    OsiCuts theseCuts ;
4611/*
4612  Scan previously generated global column and row cuts to see if any are
4613  useful.
4614*/
4615    int numberViolated=0;
4616    if (currentPassNumber_ == 1 && howOftenGlobalScan_ > 0 &&
4617        (numberNodes_%howOftenGlobalScan_) == 0)
4618    { int numberCuts = globalCuts_.sizeColCuts() ;
4619      int i;
4620      // possibly extend whichGenerator
4621      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
4622      for ( i = 0 ; i < numberCuts ; i++)
4623      { OsiColCut *thisCut = globalCuts_.colCutPtr(i) ;
4624        if (thisCut->violated(cbcColSolution_)>primalTolerance) {
4625          printf("Global cut added - violation %g\n",
4626                 thisCut->violated(cbcColSolution_)) ;
4627          whichGenerator_[numberViolated++]=-1;
4628#ifndef GLOBAL_CUTS_JUST_POINTERS
4629          theseCuts.insert(*thisCut) ;
4630#else
4631          theseCuts.insert(thisCut) ;
4632#endif
4633        }
4634      }
4635      numberCuts = globalCuts_.sizeRowCuts() ;
4636      // possibly extend whichGenerator
4637      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
4638      for ( i = 0;i<numberCuts;i++) {
4639        OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ;
4640        if (thisCut->violated(cbcColSolution_)>primalTolerance) {
4641          //printf("Global cut added - violation %g\n",
4642          // thisCut->violated(cbcColSolution_)) ;
4643          whichGenerator_[numberViolated++]=-1;
4644#ifndef GLOBAL_CUTS_JUST_POINTERS
4645          theseCuts.insert(*thisCut) ;
4646#else
4647          theseCuts.insert(thisCut) ;
4648#endif
4649        }
4650      }
4651      numberGlobalViolations_+=numberViolated;
4652    }
4653/*
4654  Generate new cuts (global and/or local) and/or apply heuristics.  If
4655  CglProbing is used, then it should be first as it can fix continuous
4656  variables.
4657
4658  At present, CglProbing is the only case where generateCuts will return
4659  true. generateCuts actually modifies variable bounds in the solver when
4660  CglProbing indicates that it can fix a variable. Reoptimisation is required
4661  to take full advantage.
4662
4663  The need to resolve here should only happen after a heuristic solution.
4664  (Note default OSI implementation of optimalBasisIsAvailable always returns
4665  false.)
4666*/
4667    if (solverCharacteristics_->warmStart()&&
4668        !solver_->optimalBasisIsAvailable()) {
4669      //printf("XXXXYY no opt basis\n");
4670      resolve(node ? node->nodeInfo() : NULL,3);
4671    }
4672    if (nextRowCut_) {
4673      // branch was a cut - add it
4674      theseCuts.insert(*nextRowCut_);
4675      if (handler_->logLevel()>1)
4676        nextRowCut_->print();
4677      const OsiRowCut * cut=nextRowCut_;
4678      double lb = cut->lb();
4679      double ub = cut->ub();
4680      int n=cut->row().getNumElements();
4681      const int * column = cut->row().getIndices();
4682      const double * element = cut->row().getElements();
4683      double sum=0.0;
4684      for (int i=0;i<n;i++) {
4685        int iColumn = column[i];
4686        double value = element[i];
4687        //if (cbcColSolution_[iColumn]>1.0e-7)
4688        //printf("value of %d is %g\n",iColumn,cbcColSolution_[iColumn]);
4689        sum += value * cbcColSolution_[iColumn];
4690      }
4691      delete nextRowCut_;
4692      nextRowCut_=NULL;
4693      if (handler_->logLevel()>1)
4694        printf("applying branch cut, sum is %g, bounds %g %g\n",sum,lb,ub);
4695      // possibly extend whichGenerator
4696      resizeWhichGenerator(numberViolated, numberViolated+1);
4697      // set whichgenerator (also serves as marker to say don't delete0
4698      whichGenerator_[numberViolated++]=-2;
4699    }
4700
4701    // reset probing info
4702    if (probingInfo_)
4703      probingInfo_->initializeFixing();
4704    int i;
4705    if ((threadMode_&2)==0||numberNodes_) {
4706      for (i = 0;i<numberCutGenerators_;i++) {
4707        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
4708        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
4709        bool generate = generator_[i]->normal();
4710        // skip if not optimal and should be (maybe a cut generator has fixed variables)
4711        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
4712          generate=false;
4713        if (generator_[i]->switchedOff())
4714          generate=false;;
4715        if (generate) {
4716          bool mustResolve = 
4717            generator_[i]->generateCuts(theseCuts,fullScan,solver_,node) ;
4718          int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
4719          if(numberRowCutsBefore < numberRowCutsAfter &&
4720             generator_[i]->mustCallAgain())
4721            keepGoing=true; // say must go round
4722          // Check last cut to see if infeasible
4723          if(numberRowCutsBefore < numberRowCutsAfter) {
4724            const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
4725            if (thisCut->lb()>thisCut->ub()) {
4726              feasible = false; // sub-problem is infeasible
4727              break;
4728            }
4729          }
4730#ifdef CBC_DEBUG
4731          {
4732            int k ;
4733            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
4734              OsiRowCut thisCut = theseCuts.rowCut(k) ;
4735              /* check size of elements.
4736                 We can allow smaller but this helps debug generators as it
4737                 is unsafe to have small elements */
4738              int n=thisCut.row().getNumElements();
4739              const int * column = thisCut.row().getIndices();
4740              const double * element = thisCut.row().getElements();
4741              //assert (n);
4742              for (int i=0;i<n;i++) {
4743                double value = element[i];
4744                assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
4745              }
4746            }
4747          }
4748#endif
4749          if (mustResolve) {
4750            int returncode = resolve(node ? node->nodeInfo() : NULL,2);
4751            feasible = returnCode  != 0 ;
4752            if (returncode<0)
4753              numberTries=0;
4754            if ((specialOptions_&1)!=0) {
4755              debugger = solver_->getRowCutDebugger() ;
4756              if (debugger) 
4757                onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
4758              else
4759                onOptimalPath=false;
4760              if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
4761                assert(feasible) ;
4762            }
4763            if (!feasible)
4764              break ;
4765          }
4766        }
4767        int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
4768        int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
4769       
4770        if ((specialOptions_&1)!=0) {
4771          if (onOptimalPath) {
4772            int k ;
4773            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
4774              OsiRowCut thisCut = theseCuts.rowCut(k) ;
4775              if(debugger->invalidCut(thisCut)) {
4776                solver_->writeMps("badCut");
4777#ifdef NDEBUG
4778                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
4779                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
4780                abort();
4781#endif
4782              }
4783              assert(!debugger->invalidCut(thisCut)) ;
4784            }
4785          }
4786        }
4787/*
4788  The cut generator has done its thing, and maybe it generated some
4789  cuts.  Do a bit of bookkeeping: load
4790  whichGenerator[i] with the index of the generator responsible for a cut,
4791  and place cuts flagged as global in the global cut pool for the model.
4792
4793  lastNumberCuts is the sum of cuts added in previous iterations; it's the
4794  offset to the proper starting position in whichGenerator.
4795*/
4796        int numberBefore =
4797          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
4798        int numberAfter =
4799          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
4800        // possibly extend whichGenerator
4801        resizeWhichGenerator(numberBefore, numberAfter);
4802        int j ;
4803        if (fullScan) {
4804          // counts
4805          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
4806        }
4807        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
4808       
4809        bool dodgyCuts=false;
4810        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
4811          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
4812          if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) {
4813            dodgyCuts=true;
4814            break;
4815          }
4816          whichGenerator_[numberBefore++] = i ;
4817          if (thisCut->lb()>thisCut->ub())
4818            violated=-2; // sub-problem is infeasible
4819          if (thisCut->globallyValid()) {
4820            // add to global list
4821            OsiRowCut newCut(*thisCut);
4822            newCut.setGloballyValid(2);
4823            newCut.mutableRow().setTestForDuplicateIndex(false);
4824            globalCuts_.insert(newCut) ;
4825          }
4826        }
4827        if (dodgyCuts) {
4828          for (int k=numberRowCutsAfter-1;k>=j;k--) {
4829            const OsiRowCut * thisCut = theseCuts.rowCutPtr(k) ;
4830            if (thisCut->lb()>thisCut->ub())
4831              violated=-2; // sub-problem is infeasible
4832            if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) 
4833              theseCuts.eraseRowCut(k);
4834          }
4835          numberRowCutsAfter = theseCuts.sizeRowCuts() ;
4836          for (;j<numberRowCutsAfter;j++) {
4837            const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
4838            whichGenerator_[numberBefore++] = i ;
4839            if (thisCut->globallyValid()) {
4840              // add to global list
4841              OsiRowCut newCut(*thisCut);
4842              newCut.setGloballyValid(2);
4843              newCut.mutableRow().setTestForDuplicateIndex(false);
4844              globalCuts_.insert(newCut) ;
4845            }
4846          }
4847        }
4848        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
4849          whichGenerator_[numberBefore++] = i ;
4850          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
4851          if (thisCut->globallyValid()) {
4852            // add to global list
4853            OsiColCut newCut(*thisCut);
4854            newCut.setGloballyValid(2);
4855            globalCuts_.insert(newCut) ;
4856          }
4857        }
4858      }
4859      // Add in any violated saved cuts
4860      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
4861        int numberOld = theseCuts.sizeRowCuts();
4862        int numberCuts = slackCuts.sizeRowCuts() ;
4863        int i;
4864        // possibly extend whichGenerator
4865        resizeWhichGenerator(numberOld, numberOld+numberCuts);
4866        for ( i = 0;i<numberCuts;i++) {
4867          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
4868          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
4869            if (messageHandler()->logLevel()>2)
4870              printf("Old cut added - violation %g\n",
4871                     thisCut->violated(cbcColSolution_)) ;
4872            whichGenerator_[numberOld++]=-1;
4873            theseCuts.insert(*thisCut) ;
4874          }
4875        }
4876      }
4877    } else {
4878      // do cuts independently
4879      OsiCuts * eachCuts = new OsiCuts [numberCutGenerators_];;
4880#ifdef CBC_THREAD
4881      if (!threadModel) {
4882#endif
4883        // generate cuts
4884        for (i = 0;i<numberCutGenerators_;i++) {
4885          bool generate = generator_[i]->normal();
4886          // skip if not optimal and should be (maybe a cut generator has fixed variables)
4887          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
4888            generate=false;
4889          if (generator_[i]->switchedOff())
4890            generate=false;;
4891          if (generate) 
4892            generator_[i]->generateCuts(eachCuts[i],fullScan,solver_,node) ;
4893        }
4894#ifdef CBC_THREAD
4895      } else {
4896        for (i=0;i<numberThreads_;i++) {
4897          // set solver here after cloning
4898          threadModel[i]->solver_=solver_->clone();
4899          threadModel[i]->numberNodes_ = (fullScan) ? 1 : 0;
4900        }
4901        // generate cuts
4902        for (i = 0;i<numberCutGenerators_;i++) {
4903          bool generate = generator_[i]->normal();
4904          // skip if not optimal and should be (maybe a cut generator has fixed variables)
4905          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
4906            generate=false;
4907          if (generator_[i]->switchedOff())
4908            generate=false;;
4909          if (generate) { 
4910            bool finished=false;
4911            int iThread=-1;
4912            // see if any available
4913            for (iThread=0;iThread<numberThreads_;iThread++) {
4914              if (threadInfo[iThread].returnCode) {
4915                finished=true;
4916                break;
4917              } else if (threadInfo[iThread].returnCode==0) {
4918                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4919              }
4920            }
4921            while (!finished) {
4922              pthread_mutex_lock(&condition_mutex);
4923              struct timespec absTime;
4924              clock_gettime(CLOCK_REALTIME,&absTime);
4925              absTime.tv_nsec += 1000000; // millisecond
4926              if (absTime.tv_nsec>=1000000000) {
4927                absTime.tv_nsec -= 1000000000;
4928                absTime.tv_sec++;
4929              }
4930              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
4931              pthread_mutex_unlock(&condition_mutex);
4932              for (iThread=0;iThread<numberThreads_;iThread++) {
4933                if (threadInfo[iThread].returnCode>0) {
4934                  finished=true;
4935                  break;
4936                } else if (threadInfo[iThread].returnCode==0) {
4937                  pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4938                }
4939              }
4940            }
4941            assert (iThread<numberThreads_);
4942            assert (threadInfo[iThread].returnCode);
4943            threadModel[iThread]->generator_[0]=generator_[i];
4944            threadModel[iThread]->object_ = (OsiObject **) (eachCuts+i);
4945            // allow to start
4946            threadInfo[iThread].returnCode=0;
4947            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4948          }
4949        }
4950        // wait
4951        for (int iThread=0;iThread<numberThreads_;iThread++) {
4952          if (threadInfo[iThread].returnCode==0) {
4953            bool finished=false;
4954            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4955            while (!finished) {
4956              pthread_mutex_lock(&condition_mutex);
4957              struct timespec absTime;
4958              clock_gettime(CLOCK_REALTIME,&absTime);
4959              absTime.tv_nsec += 1000000; // millisecond
4960              if (absTime.tv_nsec>=1000000000) {
4961                absTime.tv_nsec -= 1000000000;
4962                absTime.tv_sec++;
4963              }
4964              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
4965              pthread_mutex_unlock(&condition_mutex);
4966              if (threadInfo[iThread].returnCode>0) {
4967                finished=true;
4968                break;
4969              } else if (threadInfo[iThread].returnCode==0) {
4970                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
4971              }
4972            }
4973          }
4974          assert (threadInfo[iThread].returnCode);
4975          // say available
4976          threadInfo[iThread].returnCode=-1;
4977          delete threadModel[iThread]->solver_;
4978          threadModel[iThread]->solver_=NULL;
4979        }
4980      }
4981#endif
4982      // Now put together
4983      for (i = 0;i<numberCutGenerators_;i++) {
4984        // add column cuts
4985        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
4986        int numberColumnCuts = eachCuts[i].sizeColCuts();
4987        int numberColumnCutsAfter = numberColumnCutsBefore
4988          + numberColumnCuts;
4989        int j;
4990        for (j=0;j<numberColumnCuts;j++) {
4991          theseCuts.insert(eachCuts[i].colCut(j));
4992        }
4993        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
4994        int numberRowCuts = eachCuts[i].sizeRowCuts();
4995        int numberRowCutsAfter = numberRowCutsBefore
4996          + numberRowCuts;
4997        if (numberRowCuts) {
4998          for (j=0;j<numberRowCuts;j++) {
4999            const OsiRowCut * thisCut = eachCuts[i].rowCutPtr(j) ;
5000            if (thisCut->lb()<=1.0e10&&thisCut->ub()>=-1.0e10) 
5001              theseCuts.insert(eachCuts[i].rowCut(j));
5002          }
5003          if (generator_[i]->mustCallAgain())
5004            keepGoing=true; // say must go round
5005          // Check last cut to see if infeasible
5006          const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
5007          if (thisCut->lb()>thisCut->ub()) {
5008            feasible = false; // sub-problem is infeasible
5009            break;
5010          }
5011        }
5012#ifdef CBC_DEBUG
5013        {
5014          int k ;
5015          for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
5016            OsiRowCut thisCut = theseCuts.rowCut(k) ;
5017            /* check size of elements.
5018               We can allow smaller but this helps debug generators as it
5019               is unsafe to have small elements */
5020            int n=thisCut.row().getNumElements();
5021            const int * column = thisCut.row().getIndices();
5022            const double * element = thisCut.row().getElements();
5023            //assert (n);
5024            for (int i=0;i<n;i++) {
5025              double value = element[i];
5026              assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
5027            }
5028          }
5029        }
5030#endif
5031        if ((specialOptions_&1)!=0) {
5032          if (onOptimalPath) {
5033            int k ;
5034            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
5035              OsiRowCut thisCut = theseCuts.rowCut(k) ;
5036              if(debugger->invalidCut(thisCut)) {
5037                solver_->writeMps("badCut");
5038#ifdef NDEBUG
5039                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
5040                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
5041                abort();
5042#endif
5043              }
5044              assert(!debugger->invalidCut(thisCut)) ;
5045            }
5046          }
5047        }
5048/*
5049  The cut generator has done its thing, and maybe it generated some
5050  cuts.  Do a bit of bookkeeping: load
5051  whichGenerator[i] with the index of the generator responsible for a cut,
5052  and place cuts flagged as global in the global cut pool for the model.
5053
5054  lastNumberCuts is the sum of cuts added in previous iterations; it's the
5055  offset to the proper starting position in whichGenerator.
5056*/
5057        int numberBefore =
5058          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
5059        int numberAfter =
5060          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
5061        // possibly extend whichGenerator
5062        resizeWhichGenerator(numberBefore, numberAfter);
5063        if (fullScan) {
5064          // counts
5065          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
5066        }
5067        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
5068       
5069        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
5070          whichGenerator_[numberBefore++] = i ;
5071          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
5072          if (thisCut->lb()>thisCut->ub())
5073            violated=-2; // sub-problem is infeasible
5074          if (thisCut->globallyValid()) {
5075            // add to global list
5076            OsiRowCut newCut(*thisCut);
5077            newCut.setGloballyValid(2);
5078            newCut.mutableRow().setTestForDuplicateIndex(false);
5079            globalCuts_.insert(newCut) ;
5080          }
5081        }
5082        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
5083          whichGenerator_[numberBefore++] = i ;
5084          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
5085          if (thisCut->globallyValid()) {
5086            // add to global list
5087            OsiColCut newCut(*thisCut);
5088            newCut.setGloballyValid(2);
5089            globalCuts_.insert(newCut) ;
5090          }
5091        }
5092      }
5093      // Add in any violated saved cuts
5094      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
5095        int numberOld = theseCuts.sizeRowCuts();
5096        int numberCuts = slackCuts.sizeRowCuts() ;
5097        int i;
5098        // possibly extend whichGenerator
5099        resizeWhichGenerator(numberOld, numberOld+numberCuts);
5100        for ( i = 0;i<numberCuts;i++) {
5101          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
5102          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
5103            if (messageHandler()->logLevel()>2)
5104              printf("Old cut added - violation %g\n",
5105                     thisCut->violated(cbcColSolution_)) ;
5106            whichGenerator_[numberOld++]=-1;
5107            theseCuts.insert(*thisCut) ;
5108          }
5109        }
5110      }
5111      delete [] eachCuts;
5112    }
5113    //if (!feasible)
5114    //break;
5115/*
5116  End of the loop to exercise each generator - try heuristics
5117  - unless at root node and first pass
5118*/
5119    if (numberNodes_||currentPassNumber_!=1) {
5120      double * newSolution = new double [numberColumns] ;
5121      double heuristicValue = getCutoff() ;
5122      int found = -1; // no solution found
5123      for (i = 0;i<numberHeuristics_;i++) {
5124        // see if heuristic will do anything
5125        double saveValue = heuristicValue ;
5126        int ifSol = 
5127          heuristic_[i]->solution(heuristicValue,
5128                                  newSolution,
5129                                  theseCuts) ;
5130        if (ifSol>0) {
5131          // better solution found
5132          found = i ;
5133          incrementUsed(newSolution);
5134        } else if (ifSol<0) {
5135          heuristicValue = saveValue ;
5136        }
5137      }
5138/*
5139  Did any of the heuristics turn up a new solution? Record it before we free
5140  the vector.
5141*/
5142      if (found >= 0) { 
5143        phase_=4;
5144        incrementUsed(newSolution);
5145        lastHeuristic_ = heuristic_[found];
5146        setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
5147        CbcTreeLocal * tree
5148          = dynamic_cast<CbcTreeLocal *> (tree_);
5149        if (tree)
5150          tree->passInSolution(bestSolution_,heuristicValue);
5151      }
5152      delete [] newSolution ;
5153    }
5154
5155#if 0
5156    // switch on to get all cuts printed
5157    theseCuts.printCuts() ;
5158#endif
5159    int numberColumnCuts = theseCuts.sizeColCuts() ;
5160    int numberRowCuts = theseCuts.sizeRowCuts() ;
5161    if (violated>=0)
5162      violated = numberRowCuts + numberColumnCuts ;
5163/*
5164  Apply column cuts (aka bound tightening). This may be partially redundant
5165  for column cuts returned by CglProbing, as generateCuts installs bounds
5166  from CglProbing when it determines it can fix a variable.
5167
5168  TODO: Looks like the use of violated has evolved. The value set above is
5169        completely ignored. All that's left is violated == -1 indicates some
5170        cut is violated, violated == -2 indicates infeasibility. Only
5171        infeasibility warrants exceptional action.
5172
5173  TODO: Strikes me that this code will fail to detect infeasibility, because
5174        the breaks escape the inner loops but the outer loop keeps going.
5175        Infeasibility in an early cut will be overwritten if a later cut is
5176        merely violated.
5177*/
5178    if (numberColumnCuts) {
5179
5180#ifdef CBC_DEBUG
5181      double * oldLower = new double [numberColumns] ;
5182      double * oldUpper = new double [numberColumns] ;
5183      memcpy(oldLower,cbcColLower_,numberColumns*sizeof(double)) ;
5184      memcpy(oldUpper,cbcColUpper_,numberColumns*sizeof(double)) ;
5185#endif
5186
5187      double integerTolerance = getDblParam(CbcIntegerTolerance) ;
5188      for (int i = 0;i<numberColumnCuts;i++) {
5189        const OsiColCut * thisCut = theseCuts.colCutPtr(i) ;
5190        const CoinPackedVector & lbs = thisCut->lbs() ;
5191        const CoinPackedVector & ubs = thisCut->ubs() ;
5192        int j ;
5193        int n ;
5194        const int * which ;
5195        const double * values ;
5196        n = lbs.getNumElements() ;
5197        which = lbs.getIndices() ;
5198        values = lbs.getElements() ;
5199        for (j = 0;j<n;j++) {
5200          int iColumn = which[j] ;
5201          double value = cbcColSolution_[iColumn] ;
5202#if CBC_DEBUG>1
5203          printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
5204                 cbcColSolution_[iColumn],oldUpper[iColumn],values[j]) ;
5205#endif
5206          solver_->setColLower(iColumn,values[j]) ;
5207          if (value<values[j]-integerTolerance)
5208            violated = -1 ;
5209          if (values[j]>cbcColUpper_[iColumn]+integerTolerance) {
5210            // infeasible
5211            violated = -2 ;
5212            break ;
5213          }
5214        }
5215        n = ubs.getNumElements() ;
5216        which = ubs.getIndices() ;
5217        values = ubs.getElements() ;
5218        for (j = 0;j<n;j++) {
5219          int iColumn = which[j] ;
5220          double value = cbcColSolution_[iColumn] ;
5221#if CBC_DEBUG>1
5222          printf("%d %g %g %g %g\n",iColumn,oldLower[iColumn],
5223                 cbcColSolution_[iColumn],oldUpper[iColumn],values[j]) ;
5224#endif
5225          solver_->setColUpper(iColumn,values[j]) ;
5226          if (value>values[j]+integerTolerance)
5227            violated = -1 ;
5228          if (values[j]<cbcColLower_[iColumn]-integerTolerance) {
5229            // infeasible
5230            violated = -2 ;
5231            break ;
5232          }
5233        }
5234      }
5235#ifdef CBC_DEBUG
5236      delete [] oldLower ;
5237      delete [] oldUpper ;
5238#endif
5239    }
5240/*
5241  End installation of column cuts. The break here escapes the numberTries
5242  loop.
5243*/
5244    if (violated == -2||!feasible) {
5245      // infeasible
5246      feasible = false ;
5247      violated = -2;
5248      if (!numberNodes_) 
5249        messageHandler()->message(CBC_INFEAS,
5250                                  messages())
5251                                    << CoinMessageEol ;
5252      break ;
5253    }
5254/*
5255  Now apply the row (constraint) cuts. This is a bit more work because we need
5256  to obtain and augment the current basis.
5257
5258  TODO: Why do this work, if there are no row cuts? The current basis will do
5259        just fine.
5260*/
5261    int numberRowsNow = solver_->getNumRows() ;
5262    assert(numberRowsNow == numberRowsAtStart+lastNumberCuts) ;
5263    int numberToAdd = theseCuts.sizeRowCuts() ;
5264    numberNewCuts_ = lastNumberCuts+numberToAdd ;
5265/*
5266  Now actually add the row cuts and reoptimise.
5267
5268  Install the cuts in the solver using applyRowCuts and
5269  augment the basis with the corresponding slack. We also add each row cut to
5270  the set of row cuts (cuts.insert()) supplied as a parameter. The new basis
5271  must be set with setWarmStart().
5272
5273  TODO: Seems to me the original code could allocate addCuts with size 0, if
5274        numberRowCuts was 0 and numberColumnCuts was nonzero. That might
5275        explain the memory fault noted in the comment by AJK.  Unfortunately,
5276        just commenting out the delete[] results in massive memory leaks. Try
5277        a revision to separate the row cut case. Why do we need addCuts at
5278        all? A typing issue, apparently: OsiCut vs. OsiRowCut.
5279
5280  TODO: It looks to me as if numberToAdd and numberRowCuts are identical at
5281        this point. Confirm & get rid of one of them.
5282
5283  TODO: Any reason why the three loops can't be consolidated?
5284*/
5285    if (numberRowCuts > 0 || numberColumnCuts > 0)
5286    { if (numberToAdd > 0)
5287      { int i ;
5288        // Faster to add all at once
5289        const OsiRowCut ** addCuts = new const OsiRowCut * [numberToAdd] ;
5290        for (i = 0 ; i < numberToAdd ; i++)
5291        { addCuts[i] = &theseCuts.rowCut(i) ; }
5292        solver_->applyRowCuts(numberToAdd,addCuts) ;
5293        // AJK this caused a memory fault on Win32
5294        // may be okay now with ** form
5295        delete [] addCuts ;
5296        for (i = 0 ; i < numberToAdd ; i++)
5297        { cuts.insert(theseCuts.rowCut(i)) ; }
5298        CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5299        assert(basis != NULL); // make sure not volume
5300/* dylp bug
5301
5302  Consistent size used by OsiDylp as sanity check. Implicit resize seen
5303  as an error. Hence this call to resize is necessary.
5304*/
5305        basis->resize(numberRowsAtStart+numberNewCuts_,numberColumns) ;
5306        for (i = 0 ; i < numberToAdd ; i++) 
5307        { basis->setArtifStatus(numberRowsNow+i,
5308                                 CoinWarmStartBasis::basic) ; }
5309        if (solver_->setWarmStart(basis) == false)
5310        { throw CoinError("Fail setWarmStart() after cut installation.",
5311                          "solveWithCuts","CbcModel") ; }
5312        delete basis;
5313      }
5314      feasible = resolve(node ? node->nodeInfo() : NULL,2) ;
5315      if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] )
5316        numberTries=0; // exit
5317#     ifdef CBC_DEBUG
5318      printf("Obj value after cuts %g %d rows\n",solver_->getObjValue(),
5319              solver_->getNumRows()) ;
5320      if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
5321        assert(feasible) ;
5322#     endif
5323    }
5324/*
5325  No cuts. Cut short the cut generation (numberTries) loop.
5326*/
5327    else
5328    { numberTries = 0 ; }
5329/*
5330  If the problem is still feasible, first, call takeOffCuts() to remove cuts
5331  that are now slack. takeOffCuts() will call the solver to reoptimise if
5332  that's needed to restore a valid solution.
5333
5334  Next, see if we should quit due to diminishing returns:
5335    * we've tried three rounds of cut generation and we're getting
5336      insufficient improvement in the objective; or
5337    * we generated no cuts; or
5338    * the solver declared optimality with 0 iterations after we added the
5339      cuts generated in this round.
5340  If we decide to keep going, prep for the next iteration.
5341
5342  It just seems more safe to tell takeOffCuts() to call resolve(), even if
5343  we're not continuing cut generation. Otherwise code executed between here
5344  and final disposition of the node will need to be careful not to access the
5345  lp solution. It can happen that we lose feasibility in takeOffCuts ---
5346  numerical jitters when the cutoff bound is epsilon less than the current
5347  best, and we're evaluating an alternative optimum.
5348
5349  TODO: After successive rounds of code motion, there seems no need to
5350        distinguish between the two checks for aborting the cut generation
5351        loop. Confirm and clean up.
5352*/
5353    if (feasible)
5354    { int cutIterations = solver_->getIterationCount() ;
5355      if (numberOldActiveCuts_+numberNewCuts_) {
5356        OsiCuts * saveCuts = node ? NULL : &slackCuts;
5357        takeOffCuts(cuts,resolveAfterTakeOffCuts_,saveCuts) ;
5358        if (solver_->isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_)
5359          { feasible = false ;
5360#       ifdef CBC_DEBUG
5361          double z = solver_->getObjValue() ;
5362          double cut = getCutoff() ;
5363          printf("Lost feasibility by %g in takeOffCuts; z = %g, cutoff = %g\n",
5364                 z-cut,z,cut) ;
5365#       endif
5366          }
5367      }
5368      if (feasible) {
5369        numberRowsAtStart = numberOldActiveCuts_+numberRowsAtContinuous_ ;
5370        lastNumberCuts = numberNewCuts_ ;
5371        if (direction*solver_->getObjValue() < lastObjective+minimumDrop &&
5372            currentPassNumber_ >= 3 && !keepGoing)
5373          { numberTries = 0 ; }
5374        if (numberRowCuts+numberColumnCuts == 0 || cutIterations == 0)
5375          { break ; }
5376        if (numberTries > 0) {
5377          reducedCostFix() ;
5378          lastObjective = direction*solver_->getObjValue() ;
5379        }
5380      }
5381    }
5382/*
5383  We've lost feasibility --- this node won't be referencing the cuts we've
5384  been collecting, so decrement the reference counts.
5385*/
5386    if (!feasible)
5387    { int i ;
5388      if (currentNumberCuts_) {
5389        lockThread();
5390        for (i = 0;i<currentNumberCuts_;i++) {
5391          // take off node
5392          if (addedCuts_[i]) {
5393            if (!addedCuts_[i]->decrement())
5394              delete addedCuts_[i] ;
5395            addedCuts_[i] = NULL ;
5396          }
5397        }
5398      }
5399      numberTries = 0 ;
5400    }
5401  } while (numberTries>0||keepGoing) ;
5402  {
5403    // switch on
5404    for (int i = 0;i<numberCutGenerators_;i++) 
5405      generator_[i]->setSwitchedOff(false);
5406  }
5407  //check feasibility.
5408  //If solution seems to be integer feasible calling setBestSolution
5409  //will eventually add extra global cuts which we need to install at
5410  //the nodes
5411
5412
5413  if(feasible&&solverCharacteristics_->solutionAddsCuts())  //check integer feasibility
5414  {
5415    bool integerFeasible = true;
5416    const double * save = testSolution_;
5417    testSolution_ = solver_->getColSolution();
5418    // point to useful information
5419    OsiBranchingInformation usefulInfo=usefulInformation();
5420    for (int i=0;i<numberObjects_ && integerFeasible;i++)
5421    {
5422      int preferredWay;
5423      double infeasibility = object_[i]->infeasibility(&usefulInfo,preferredWay);
5424      if(infeasibility)
5425        integerFeasible = false;
5426    }
5427    testSolution_ = save;
5428    if(integerFeasible)//update
5429    {
5430      double objValue = solver_->getObjValue();
5431      int numberGlobalBefore = globalCuts_.sizeRowCuts();
5432      // SOLUTION2 so won't up cutoff or print message
5433      setBestSolution(CBC_SOLUTION2, objValue, 
5434                      solver_->getColSolution(),0);
5435      int numberGlobalAfter = globalCuts_.sizeRowCuts();
5436      int numberToAdd = numberGlobalAfter - numberGlobalBefore;
5437      if(numberToAdd > 0)
5438        //We have added some cuts say they are tight at that node
5439        //Basis and lp should already have been updated
5440      {
5441        feasible = (solver_->isProvenOptimal() &&
5442                    !solver_->isDualObjectiveLimitReached()) ;
5443        if(feasible)
5444        {
5445          int numberCuts = numberNewCuts_ = cuts.sizeRowCuts();
5446      // possibly extend whichGenerator
5447          resizeWhichGenerator(numberCuts, numberToAdd+numberCuts);
5448         
5449          for (int i = numberGlobalBefore ; i < numberGlobalAfter ; i++)
5450            {       
5451              whichGenerator_[numberNewCuts_++]=-1;
5452#ifndef GLOBAL_CUTS_JUST_POINTERS
5453              cuts.insert(globalCuts_.rowCut(i)) ; 
5454#else
5455              OsiRowCut * rowCutPointer = globalCuts_.rowCutPtr(i);
5456              cuts.insert(rowCutPointer) ; 
5457#endif
5458            }
5459          numberNewCuts_ = lastNumberCuts+numberToAdd;
5460          //now take off the cuts which are not tight anymore
5461          takeOffCuts(cuts,resolveAfterTakeOffCuts_, NULL) ;
5462          if (solver_->isDualObjectiveLimitReached()&&resolveAfterTakeOffCuts_)
5463            { feasible = false ;}
5464        }
5465        if(!feasible) //node will be fathomed
5466          { 
5467          for (int i = 0;i<currentNumberCuts_;i++) 
5468            {
5469            // take off node
5470            if (addedCuts_[i])
5471              {
5472                if (!addedCuts_[i]->decrement())
5473                  delete addedCuts_[i] ;
5474                addedCuts_[i] = NULL ;
5475              }
5476            }
5477        }
5478      }
5479    }
5480  }
5481/*
5482  Reduced cost fix at end. Must also check feasible, in case we've popped out
5483  because a generator indicated we're infeasible.
5484*/
5485  if (feasible && solver_->isProvenOptimal())
5486    reducedCostFix();
5487  // If at root node do heuristics
5488  if (!numberNodes_) {
5489    // First see if any cuts are slack
5490    int numberRows = solver_->getNumRows();
5491    int numberAdded = numberRows-numberRowsAtContinuous_;
5492    if (numberAdded) {
5493      CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
5494      assert(basis != NULL);
5495      int * added = new int[numberAdded];
5496      int nDelete = 0;
5497      for (int j=numberRowsAtContinuous_;j<numberRows;j++) {
5498        if (basis->getArtifStatus(j)==CoinWarmStartBasis::basic) {
5499          //printf("%d slack!\n",j);
5500          added[nDelete++]=j;
5501        }
5502      }
5503      if (nDelete) {
5504        solver_->deleteRows(nDelete,added);
5505      }
5506      delete [] added;
5507      delete basis ;
5508    }
5509    // mark so heuristics can tell
5510    int savePass=currentPassNumber_;
5511    currentPassNumber_=999999;
5512    double * newSolution = new double [numberColumns] ;
5513    double heuristicValue = getCutoff() ;
5514    int found = -1; // no solution found
5515    for (int i = 0;i<numberHeuristics_;i++) {
5516      // see if heuristic will do anything
5517      double saveValue = heuristicValue ;
5518      int ifSol = heuristic_[i]->solution(heuristicValue,
5519                                          newSolution);
5520      if (ifSol>0) {
5521        // better solution found
5522        found = i ;
5523        incrementUsed(newSolution);
5524      } else {
5525        heuristicValue = saveValue ;
5526      }
5527    }
5528    currentPassNumber_=savePass;
5529    if (found >= 0) { 
5530      phase_=4;
5531      incrementUsed(newSolution);
5532      lastHeuristic_ = heuristic_[found];
5533      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
5534    }
5535    delete [] newSolution ;
5536  }
5537  // Up change due to cuts
5538  if (feasible)
5539    sumChangeObjective2_ += solver_->getObjValue()*solver_->getObjSense()
5540      - objectiveValue ;
5541  //if ((numberNodes_%100)==0)
5542  //printf("XXb sum obj changed by %g\n",sumChangeObjective2_);
5543/*
5544  End of cut generation loop.
5545
5546  Now, consider if we want to disable or adjust the frequency of use for any
5547  of the cut generators. If the client specified a positive number for
5548  howOften, it will never change. If the original value was negative, it'll
5549  be converted to 1000000+|howOften|, and this value will be adjusted each
5550  time fullScan is true. Actual cut generation is performed every
5551  howOften%1000000 nodes; the 1000000 offset is just a convenient way to
5552  specify that the frequency is adjustable.
5553
5554  During cut generation, we recorded the number of cuts produced by each
5555  generator for this node. For all cuts, whichGenerator records the generator
5556  that produced a cut.
5557
5558  TODO: All this should probably be hidden in a method of the CbcCutGenerator
5559  class.
5560
5561  TODO: Can the loop that scans over whichGenerator to accumulate per generator
5562        counts be replaced by values in countRowCuts and countColumnCuts?
5563*/
5564#ifdef NODE_LOG
5565  int fatherNum = (node == NULL) ? -1 : node->nodeNumber();
5566  double value =  (node == NULL) ? -1 : node->branchingObject()->value();
5567  string bigOne = (solver_->getIterationCount() > 30) ? "*******" : "";
5568  string way = (node == NULL) ? "" : (node->branchingObject()->way())==1 ? "Down" : "Up";
5569  std::cout<<"Node "<<numberNodes_<<", father "<<fatherNum<<", #iterations "<<solver_->getIterationCount()<<", sol value : "<<solver_->getObjValue()<<std::endl;
5570#endif
5571  if (fullScan&&numberCutGenerators_) {
5572    /* If cuts just at root node then it will probably be faster to
5573       update matrix and leave all in */
5574    int willBeCutsInTree=0;
5575    double thisObjective = solver_->getObjValue()*direction ;
5576    // get sizes
5577    int numberRowsAdded = solver_->getNumRows()-numberRowsAtStart;
5578    CoinBigIndex numberElementsAdded =  solver_->getNumElements()-numberElementsAtStart ;
5579    double densityOld = ((double) numberElementsAtStart)/((double) numberRowsAtStart);
5580    double densityNew = numberRowsAdded ? ((double) (numberElementsAdded))/((double) numberRowsAdded)
5581      : 0.0;
5582    if (!numberNodes_) {
5583      if (numberRowsAdded)
5584        handler_->message(CBC_CUTS_STATS,messages_)
5585          <<numberRowsAdded
5586          <<densityNew
5587          <<CoinMessageEol ;
5588      if (thisObjective-startObjective<1.0e-5&&numberElementsAdded>0.2*numberElementsAtStart)
5589        willBeCutsInTree=-1;
5590    }
5591    if ((numberRowsAdded>100+0.5*numberRowsAtStart
5592         ||numberElementsAdded>0.5*numberElementsAtStart)
5593        &&(densityNew>200.0&&numberRowsAdded>100&&densityNew>2.0*densityOld)) {
5594      // much bigger
5595      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
5596      //willBeCutsInTree=-1;
5597      //printf("Cuts will be taken off , %d rows added with density %g\n",
5598      //     numberRowsAdded,densityNew);
5599    }
5600    if (densityNew>100.0&&numberRowsAdded>2&&densityNew>2.0*densityOld) {
5601      //if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5)
5602      //willBeCutsInTree=-2;
5603      //printf("Density says no cuts ? , %d rows added with density %g\n",
5604      //     numberRowsAdded,densityNew);
5605    }
5606    // Root node or every so often - see what to turn off
5607    int i ;
5608    for (i = 0;i<numberCutGenerators_;i++) {
5609      int howOften = generator_[i]->howOften() ;
5610      if (howOften>-90) 
5611        willBeCutsInTree=0;
5612    }
5613    if (!numberNodes_)
5614      handler_->message(CBC_ROOT,messages_)
5615        <<numberNewCuts_
5616        <<startObjective<<thisObjective
5617        <<currentPassNumber_
5618        <<CoinMessageEol ;
5619    int * count = new int[numberCutGenerators_] ;
5620    memset(count,0,numberCutGenerators_*sizeof(int)) ;
5621    int numberActiveGenerators=0;
5622    for (i = 0;i<numberNewCuts_;i++) {
5623      int iGenerator = whichGenerator_[i];
5624      if (iGenerator>=0&&iGenerator<numberCutGenerators_)
5625        count[iGenerator]++ ;
5626    }
5627    double totalCuts = 0.0 ;
5628    //#define JUST_ACTIVE
5629    for (i = 0;i<numberCutGenerators_;i++) {
5630      if (countRowCuts[i]||countColumnCuts[i])
5631        numberActiveGenerators++;
5632#ifdef JUST_ACTIVE
5633      totalCuts += count[i] + 5.0*countColumnCuts[i] ;
5634#else
5635      totalCuts += countRowCuts[i] + 5.0*countColumnCuts[i] ;
5636#endif
5637    }
5638    int iProbing=-1;
5639    double small = (0.5* totalCuts) /
5640      ((double) numberActiveGenerators) ;
5641    for (i = 0;i<numberCutGenerators_;i++) {
5642      int howOften = generator_[i]->howOften() ;
5643      /*  Probing can be set to just do column cuts in treee.
5644          But if doing good then leave as on */
5645      bool probingWasOnBut=false;
5646      CglProbing * probing = dynamic_cast<CglProbing*>(generator_[i]->generator());
5647      if (probing&&!numberNodes_) {
5648        iProbing = i;
5649        if (probing->rowCuts()==-3) {
5650          probingWasOnBut=true;
5651          howOften = -98;
5652          probing->setRowCuts(3);
5653        }
5654      }
5655      if (willBeCutsInTree<0&&howOften==-98)
5656        howOften =-99;
5657      if (howOften==-98&&generator_[i]->switchOffIfLessThan()<0) {
5658        if (thisObjective-startObjective<0.005*fabs(startObjective)+1.0e-5)
5659          howOften=-99; // switch off
5660        if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5
5661            &&5*solver_->getNumRows()<solver_->getNumCols())
5662          howOften=-99; // switch off
5663      }
5664      if (howOften<-99)
5665        continue ;
5666      if (howOften<0||howOften >= 1000000) {
5667        if( !numberNodes_) {
5668          // If small number switch mostly off
5669#ifdef JUST_ACTIVE
5670          double thisCuts = count[i] + 5.0*countColumnCuts[i] ;
5671#else
5672          double thisCuts = countRowCuts[i] + 5.0*countColumnCuts[i] ;
5673#endif
5674          if (!thisCuts||howOften == -99) {
5675            if (howOften == -99||howOften == -98) 
5676              howOften = -100 ;
5677            else
5678              howOften = 1000000+SCANCUTS; // wait until next time
5679          } else if (thisCuts<small) {
5680            if (howOften!=1&&!probingWasOnBut) {
5681              if (generator_[i]->whatDepth()<0||howOften!=-1) {
5682                int k = (int) sqrt(small/thisCuts) ;
5683                if (howOften!=-98)
5684                  howOften = k+1000000 ;
5685                else
5686                  howOften=-100;
5687              } else {
5688                howOften=1;
5689              }
5690            } else {
5691              howOften=1;
5692              // allow cuts
5693              probingWasOnBut=false;
5694            }
5695          } else {
5696            if (thisObjective-startObjective<0.1*fabs(startObjective)+1.0e-5&&generator_[i]->whatDepth()<0)
5697              generator_[i]->setWhatDepth(5);
5698            howOften = 1+1000000 ;
5699          }
5700        }
5701        // If cuts useless switch off
5702        if (numberNodes_>=100000&&sumChangeObjective1_>2.0e2*(sumChangeObjective2_+1.0e-12)) {
5703          howOften = 1000000+SCANCUTS; // wait until next time
5704          //printf("switch off cut %d due to lack of use\n",i);
5705        }
5706      }
5707      if (!numberNodes_) {
5708        if (probingWasOnBut&&howOften==-100) {
5709          probing->setRowCuts(-3);
5710          howOften=1;
5711        }
5712        if (howOften==1)
5713          generator_[i]->setWhatDepth(1);
5714       
5715        if (howOften>=0&&generator_[i]->generator()->mayGenerateRowCutsInTree())
5716          willBeCutsInTree=1;
5717      }
5718      generator_[i]->setHowOften(howOften) ;
5719      if (howOften>=1000000&&howOften<2000000&&0) {
5720        // Go to depth
5721        int bias=1;
5722        if (howOften==1+1000000)
5723          generator_[i]->setWhatDepth(bias+1);
5724        else if (howOften<=10+1000000)
5725          generator_[i]->setWhatDepth(bias+2);
5726        else
5727          generator_[i]->setWhatDepth(bias+1000);
5728      }
5729      int newFrequency = generator_[i]->howOften()%1000000 ;
5730      // increment cut counts
5731      generator_[i]->incrementNumberCutsInTotal(countRowCuts[i]);
5732      generator_[i]->incrementNumberCutsActive(count[i]);
5733      if (handler_->logLevel()>1||!numberNodes_) {
5734        handler_->message(CBC_GENERATOR,messages_)
5735          <<i
5736          <<generator_[i]->cutGeneratorName()
5737          <<countRowCuts[i]<<count[i]
5738          <<countColumnCuts[i];
5739        handler_->printing(!numberNodes_&&generator_[i]->timing())
5740          <<generator_[i]->timeInCutGenerator();
5741        handler_->message()
5742          <<newFrequency
5743          <<CoinMessageEol ;
5744      }
5745    } 
5746    delete [] count ;
5747    if( !numberNodes_) {
5748      // save statistics
5749      for (i = 0;i<numberCutGenerators_;i++) {
5750        generator_[i]->setNumberCutsAtRoot(generator_[i]->numberCutsInTotal());
5751        generator_[i]->setNumberActiveCutsAtRoot(generator_[i]->numberCutsActive());
5752      }
5753      // decide on pseudo cost strategy
5754      int howOften = iProbing>=0 ? generator_[iProbing]->howOften() : 0;
5755      if ((howOften %1000000)!=1) 
5756        howOften = 0;
5757      //if (howOften) {
5758      //CglProbing * probing = dynamic_cast<CglProbing*>(generator_[iProbing]->generator());
5759      //}
5760      howOften = 0;
5761      if (howOften) {
5762        printf("** method 1\n");
5763        //CglProbing * probing = dynamic_cast<CglProbing*>(generator_[iProbing]->generator());
5764        generator_[iProbing]->setWhatDepth(1);
5765        // could set no row cuts
5766        //if (thisObjective-startObjective<0.001*fabs(startObjective)+1.0e-5)
5767        // probing->setRowCuts(0);
5768        for (int i=0;i<numberObjects_;i++) {
5769          CbcSimpleIntegerDynamicPseudoCost * obj =
5770            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
5771          if (obj) 
5772            obj->setMethod(1);
5773        }
5774      }
5775      if (willBeCutsInTree==-2)
5776        willBeCutsInTree=0;
5777      if( willBeCutsInTree<=0) {
5778        // Take off cuts
5779        cuts = OsiCuts();
5780        numberNewCuts_=0;
5781        if (!willBeCutsInTree) {
5782          // update size of problem
5783          numberRowsAtContinuous_ = solver_->getNumRows() ;
5784        } else {
5785          // take off cuts
5786          int numberRows = solver_->getNumRows();
5787          int numberAdded = numberRows-numberRowsAtContinuous_;
5788          if (numberAdded) {
5789            int * added = new int[numberAdded];
5790            for (int i=0;i<numberAdded;i++)
5791              added[i]=i+numberRowsAtContinuous_;
5792            solver_->deleteRows(numberAdded,added);
5793            delete [] added;
5794            // resolve so optimal
5795            resolve(solver_);
5796          }
5797        }
5798#ifdef COIN_HAS_CLP
5799        OsiClpSolverInterface * clpSolver
5800          = dynamic_cast<OsiClpSolverInterface *> (solver_);
5801        if (clpSolver) {
5802          // Maybe solver might like to know only column bounds will change
5803          //int options = clpSolver->specialOptions();
5804          //clpSolver->setSpecialOptions(options|128);
5805          clpSolver->synchronizeModel();
5806        }
5807#endif
5808      } else {
5809#ifdef COIN_HAS_CLP
5810        OsiClpSolverInterface * clpSolver
5811          = dynamic_cast