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

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

faster parallel?

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