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

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

fix bug with 12 cpus

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