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

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

to make easier to use as function

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