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

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

trying to improve performance and fix message problem in CbcSolver?

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