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

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

trying to move ampl stuff away

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