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

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

possible thread loop?

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