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

Last change on this file since 1015 was 1015, checked in by forrest, 11 years ago

add ifdefs for future exploration

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