source: trunk/Cbc/src/CbcModel.cpp @ 758

Last change on this file since 758 was 758, checked in by forrest, 14 years ago

fpump - make dj fixing optional

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