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

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

on way to CbcSolver? class

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