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

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

try and speed up feasibility pump

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