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

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

some changes e.g. DINS added and a bit of printing

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 461.5 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/*
1525  Used to generate bound edits for CbcPartialNodeInfo.
1526*/
1527  double * lowerBefore = new double [numberColumns] ;
1528  double * upperBefore = new double [numberColumns] ;
1529/*
1530 
1531  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
1532  lifting. It will iterate a generate/reoptimise loop (including reduced cost
1533  fixing) until no cuts are generated, the change in objective falls off,  or
1534  the limit on the number of rounds of cut generation is exceeded.
1535
1536  At the end of all this, any cuts will be recorded in cuts and also
1537  installed in the solver's constraint system. We'll have reoptimised, and
1538  removed any slack cuts (numberOldActiveCuts_ and numberNewCuts_ have been
1539  adjusted accordingly).
1540
1541  Tell cut generators they can be a bit more aggressive at root node
1542
1543  TODO: Why don't we make a copy of the solution after solveWithCuts?
1544  TODO: If numberUnsatisfied == 0, don't we have a solution?
1545*/
1546  phase_=1;
1547  int iCutGenerator;
1548  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
1549    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
1550    generator->setAggressiveness(generator->getAggressiveness()+100);
1551  }
1552  OsiCuts cuts ;
1553  int anyAction = -1 ;
1554  numberOldActiveCuts_ = 0 ;
1555  numberNewCuts_ = 0 ;
1556  // Array to mark solution
1557  delete [] usedInSolution_;
1558  usedInSolution_ = new int[numberColumns];
1559  CoinZeroN(usedInSolution_,numberColumns);
1560/*
1561  For printing totals and for CbcNode (numberNodes_)
1562*/
1563  numberIterations_ = 0 ;
1564  numberNodes_ = 0 ;
1565  numberNodes2_ = 0 ;
1566  maximumStatistics_=0;
1567  maximumDepthActual_=0;
1568  numberDJFixed_=0.0;
1569  // Do heuristics
1570  doHeuristicsAtRoot();
1571  if ( intParam_[CbcMaxNumNode] < 0)
1572    eventHappened_=true; // stop as fast as possible
1573  statistics_ = NULL;
1574  // Do on switch
1575  if (doStatistics>0&&doStatistics<100) {
1576    maximumStatistics_=10000;
1577    statistics_ = new CbcStatistics * [maximumStatistics_];
1578    memset(statistics_,0,maximumStatistics_*sizeof(CbcStatistics *));
1579  }
1580
1581  { int iObject ;
1582    int preferredWay ;
1583    int numberUnsatisfied = 0 ;
1584    delete [] currentSolution_;
1585    currentSolution_ = new double [numberColumns];
1586    testSolution_ = currentSolution_;
1587    memcpy(currentSolution_,solver_->getColSolution(),
1588           numberColumns*sizeof(double)) ;
1589    // point to useful information
1590    OsiBranchingInformation usefulInfo=usefulInformation();
1591
1592    for (iObject = 0 ; iObject < numberObjects_ ; iObject++)
1593    { double infeasibility =
1594          object_[iObject]->infeasibility(&usefulInfo,preferredWay) ;
1595      if (infeasibility ) numberUnsatisfied++ ; }
1596    // replace solverType
1597    if(solverCharacteristics_->tryCuts())  {
1598
1599      if (numberUnsatisfied)   {
1600        // User event
1601        if (!eventHappened_)
1602          feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
1603                                   NULL);
1604        else
1605          feasible=false;
1606      } else if (solverCharacteristics_->solutionAddsCuts()||
1607                 solverCharacteristics_->alwaysTryCutsAtRootNode()) {
1608        // may generate cuts and turn the solution
1609        //to an infeasible one
1610        feasible = solveWithCuts(cuts, 1,
1611                                 NULL);
1612      }
1613    }
1614    // check extra info on feasibility
1615    if (!solverCharacteristics_->mipFeasible())
1616      feasible = false;
1617  }
1618  // make cut generators less aggressive
1619  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
1620    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
1621    generator->setAggressiveness(generator->getAggressiveness()-100);
1622  }
1623  currentNumberCuts_ = numberNewCuts_ ;
1624  // See if can stop on gap
1625  stoppedOnGap_ = false ;
1626  bestPossibleObjective_ = solver_->getObjValue()*solver_->getObjSense();
1627  double testGap = CoinMax(dblParam_[CbcAllowableGap],
1628                           CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_))
1629                           *dblParam_[CbcAllowableFractionGap]);
1630  if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
1631    if (bestPossibleObjective_<getCutoff())
1632      stoppedOnGap_ = true ;
1633    feasible = false;
1634  }
1635  // User event
1636  if (eventHappened_)
1637    feasible=false;
1638  if(fastNodeDepth_>=0&&!parentModel_) {
1639    // add in a general depth object
1640    CbcObject * obj = 
1641      new CbcGeneralDepth(this,fastNodeDepth_);
1642    addObjects(1,&obj);
1643    // mark as done
1644    fastNodeDepth_ += 1000000;
1645    delete obj;
1646    // fake number of objects
1647    numberObjects_--;
1648  }
1649#ifdef COIN_HAS_CLP
1650#ifdef NO_CRUNCH
1651 {
1652   OsiClpSolverInterface * clpSolver
1653     = dynamic_cast<OsiClpSolverInterface *> (solver_);
1654   if (clpSolver&&!parentModel_) {
1655     ClpSimplex * clpSimplex = clpSolver->getModelPtr();
1656     clpSimplex->setSpecialOptions(clpSimplex->specialOptions()|131072);
1657     //clpSimplex->startPermanentArrays();
1658     clpSimplex->setPersistenceFlag(2);
1659   }
1660 }
1661#endif
1662#endif
1663 // Save copy of solver
1664 OsiSolverInterface * saveSolver = NULL;
1665 if (!parentModel_&&(specialOptions_&(512+2048))!=0)
1666   saveSolver = solver_->clone();
1667 if (!parentModel_&&(specialOptions_&2048)!=0) {
1668   // See if worth trying reduction
1669   bool tryNewSearch=solverCharacteristics_->reducedCostsAccurate();
1670   int numberColumns = getNumCols();
1671   if (tryNewSearch) {
1672     double cutoff = getCutoff() ;
1673     saveSolver->resolve();
1674     double direction = saveSolver->getObjSense() ;
1675     double gap = cutoff - saveSolver->getObjValue()*direction ;
1676     double tolerance;
1677     saveSolver->getDblParam(OsiDualTolerance,tolerance) ;
1678     if (gap<=0.0)
1679       gap = tolerance; 
1680     gap += 100.0*tolerance;
1681     double integerTolerance = getDblParam(CbcIntegerTolerance) ;
1682     
1683     const double *lower = saveSolver->getColLower() ;
1684     const double *upper = saveSolver->getColUpper() ;
1685     const double *solution = saveSolver->getColSolution() ;
1686     const double *reducedCost = saveSolver->getReducedCost() ;
1687     
1688     int numberFixed = 0 ;
1689     int numberFixed2=0;
1690     for (int i = 0 ; i < numberIntegers_ ; i++) {
1691       int iColumn = integerVariable_[i] ;
1692       double djValue = direction*reducedCost[iColumn] ;
1693       if (upper[iColumn]-lower[iColumn] > integerTolerance) {
1694         if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap) {
1695           saveSolver->setColUpper(iColumn,lower[iColumn]) ;
1696           numberFixed++ ;
1697         } else if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap) {
1698           saveSolver->setColLower(iColumn,upper[iColumn]) ;
1699           numberFixed++ ;
1700         }
1701       } else {
1702         numberFixed2++;
1703       }
1704     }
1705#ifdef COIN_DEVELOP
1706     if ((specialOptions_&1)!=0) {
1707       const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
1708       if (debugger) { 
1709         printf("Contains optimal\n") ;
1710         saveSolver->writeMps("reduced");
1711       } else {
1712         abort();
1713       }
1714     }
1715     printf("Restart could fix %d integers (%d already fixed)\n",
1716            numberFixed+numberFixed2,numberFixed2);
1717#endif
1718     numberFixed += numberFixed2;
1719     if (numberFixed*20<numberColumns)
1720       tryNewSearch=false;
1721   }
1722   if (tryNewSearch) {
1723     // back to solver without cuts?
1724#if 0
1725     OsiSolverInterface * solver2 = continuousSolver_->clone();
1726#else
1727     OsiSolverInterface * solver2 = saveSolver->clone();
1728#endif
1729     const double *lower = saveSolver->getColLower() ;
1730     const double *upper = saveSolver->getColUpper() ;
1731     for (int i = 0 ; i < numberIntegers_ ; i++) {
1732       int iColumn = integerVariable_[i] ;
1733       solver2->setColLower(iColumn,lower[iColumn]);
1734       solver2->setColUpper(iColumn,upper[iColumn]);
1735     }
1736     // swap
1737     delete saveSolver;
1738     saveSolver=solver2;
1739     double * newSolution = new double[numberColumns];
1740     double objectiveValue=cutoff;
1741     CbcSerendipity heuristic(*this);
1742     if (bestSolution_)
1743       heuristic.setInputSolution(bestSolution_,bestObjective_);
1744     heuristic.setFractionSmall(0.9);
1745     heuristic.setFeasibilityPumpOptions(1008013);
1746     // Use numberNodes to say how many are original rows
1747     heuristic.setNumberNodes(continuousSolver_->getNumRows());
1748#ifdef COIN_DEVELOP
1749     if (continuousSolver_->getNumRows()<
1750         solver_->getNumRows())
1751       printf("%d rows added ZZZZZ\n",
1752              solver_->getNumRows()-continuousSolver_->getNumRows());
1753#endif
1754     int returnCode= heuristic.smallBranchAndBound(saveSolver,
1755                                                   -1,newSolution,
1756                                                   objectiveValue,
1757                                                   cutoff,"Reduce");
1758     if (returnCode==-1) {
1759#ifdef COIN_DEVELOP
1760       printf("Restart - not small enough to do search after fixing\n");
1761#endif
1762       delete [] newSolution;
1763     } else {
1764       assert (returnCode>=0);
1765       if ((returnCode&1)!=0) {
1766         // increment number of solutions so other heuristics can test
1767         numberSolutions_++;
1768         numberHeuristicSolutions_++;
1769         lastHeuristic_ = NULL;
1770         setBestSolution(CBC_ROUNDING,objectiveValue,newSolution) ;
1771       }
1772       delete [] newSolution;
1773       feasible=false; // stop search
1774     }
1775   } 
1776 }
1777/*
1778  We've taken the continuous relaxation as far as we can. Time to branch.
1779  The first order of business is to actually create a node. chooseBranch
1780  currently uses strong branching to evaluate branch object candidates,
1781  unless forced back to simple branching. If chooseBranch concludes that a
1782  branching candidate is monotone (anyAction == -1) or infeasible (anyAction
1783  == -2) when forced to integer values, it returns here immediately.
1784
1785  Monotone variables trigger a call to resolve(). If the problem remains
1786  feasible, try again to choose a branching variable. At the end of the loop,
1787  resolved == true indicates that some variables were fixed.
1788
1789  Loss of feasibility will result in the deletion of newNode.
1790*/
1791
1792  bool resolved = false ;
1793  CbcNode *newNode = NULL ;
1794  numberFixedAtRoot_=0;
1795  numberFixedNow_=0;
1796  int numberIterationsAtContinuous = numberIterations_;
1797  //solverCharacteristics_->setSolver(solver_);
1798  if (feasible) {
1799    if (probingInfo_) {
1800      int number01 = probingInfo_->numberIntegers();
1801      //const fixEntry * entry = probingInfo_->fixEntries();
1802      const int * toZero = probingInfo_->toZero();
1803      //const int * toOne = probingInfo_->toOne();
1804      //const int * integerVariable = probingInfo_->integerVariable();
1805      if (toZero[number01]) {
1806        for (int i = 0;i<numberCutGenerators_;i++) {
1807          CglFakeClique * clique = dynamic_cast<CglFakeClique*>(generator_[i]->generator());
1808          if (clique) {
1809            OsiSolverInterface * fakeSolver = probingInfo_->analyze(*solver_,1);
1810            if (fakeSolver) {
1811              printf("Probing fake solver has %d rows\n",fakeSolver->getNumRows());
1812              //if (fakeSolver)
1813              //fakeSolver->writeMps("bad");
1814              if (generator_[i]->numberCutsInTotal())
1815                generator_[i]->setHowOften(1);
1816            }
1817            clique->assignSolver(fakeSolver);
1818            //stored->setProbingInfo(probingInfo_);
1819            break;
1820          }
1821        }
1822      }
1823      delete probingInfo_;
1824      probingInfo_=NULL;
1825    }
1826    newNode = new CbcNode ;
1827    // Set objective value (not so obvious if NLP etc)
1828    setObjectiveValue(newNode,NULL);
1829    anyAction = -1 ;
1830    // To make depth available we may need a fake node
1831    CbcNode fakeNode;
1832    if (!currentNode_) {
1833      // Not true if sub trees assert (!numberNodes_);
1834      currentNode_=&fakeNode;
1835    }
1836    phase_=3;
1837    // only allow 1000 passes
1838    int numberPassesLeft=1000;
1839    // This is first crude step
1840    if (numberAnalyzeIterations_) {
1841      delete [] analyzeResults_;
1842      analyzeResults_ = new double [4*numberIntegers_];
1843      numberFixedAtRoot_=newNode->analyze(this,analyzeResults_);
1844      if (numberFixedAtRoot_>0) {
1845        printf("%d fixed by analysis\n",numberFixedAtRoot_);
1846        setPointers(solver_);
1847        numberFixedNow_ = numberFixedAtRoot_;
1848      } else if (numberFixedAtRoot_<0) {
1849        printf("analysis found to be infeasible\n");
1850        anyAction=-2;
1851        delete newNode ;
1852        newNode = NULL ;
1853        feasible = false ;
1854      }
1855    }
1856    OsiSolverBranch * branches = NULL;
1857    anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved,
1858                             NULL,NULL,NULL,branches);
1859    if (anyAction == -2||newNode->objectiveValue() >= cutoff) {
1860      if (anyAction != -2) {
1861        // zap parent nodeInfo
1862#ifdef COIN_DEVELOP
1863        printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
1864#endif
1865        if (newNode->nodeInfo())
1866          newNode->nodeInfo()->nullParent();
1867      }
1868      delete newNode ;
1869      newNode = NULL ;
1870      feasible = false ;
1871    }
1872  }
1873/*
1874  At this point, the root subproblem is infeasible or fathomed by bound
1875  (newNode == NULL), or we're live with an objective value that satisfies the
1876  current objective cutoff.
1877*/
1878  assert (!newNode || newNode->objectiveValue() <= cutoff) ;
1879  // Save address of root node as we don't want to delete it
1880  // initialize for print out
1881  int lastDepth=0;
1882  int lastUnsatisfied=0;
1883  if (newNode)
1884    lastUnsatisfied=newNode->numberUnsatisfied();
1885/*
1886  The common case is that the lp relaxation is feasible but doesn't satisfy
1887  integrality (i.e., newNode->branchingObject(), indicating we've been able to
1888  select a branching variable). Remove any cuts that have gone slack due to
1889  forcing monotone variables. Then tack on an CbcFullNodeInfo object and full
1890  basis (via createInfo()) and stash the new cuts in the nodeInfo (via
1891  addCuts()). If, by some miracle, we have an integral solution at the root
1892  (newNode->branchingObject() is NULL), takeOffCuts() will ensure that the solver holds
1893  a valid solution for use by setBestSolution().
1894*/
1895  CoinWarmStartBasis *lastws = NULL ;
1896  if (feasible && newNode->branchingObject())
1897  { if (resolved)
1898    { takeOffCuts(cuts,false,NULL) ;
1899#     ifdef CHECK_CUT_COUNTS
1900      { printf("Number of rows after chooseBranch fix (root)"
1901               "(active only) %d\n",
1902                numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ;
1903        const CoinWarmStartBasis* debugws =
1904          dynamic_cast <const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
1905        debugws->print() ;
1906        delete debugws ; }
1907#     endif
1908    }
1909  //newNode->createInfo(this,NULL,NULL,NULL,NULL,0,0) ;
1910    //newNode->nodeInfo()->addCuts(cuts,
1911    //                   newNode->numberBranches(),whichGenerator_) ;
1912    if (lastws) delete lastws ;
1913    lastws = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
1914  }
1915/*
1916  Continuous data to be used later
1917*/
1918  continuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
1919  continuousInfeasibilities_ = 0 ;
1920  if (newNode)
1921  { continuousObjective_ = newNode->objectiveValue() ;
1922    delete [] continuousSolution_;
1923    continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
1924                                             numberColumns);
1925    continuousInfeasibilities_ = newNode->numberUnsatisfied() ; }
1926/*
1927  Bound may have changed so reset in objects
1928*/
1929  { int i ;
1930    for (i = 0;i < numberObjects_;i++)
1931      object_[i]->resetBounds(solver_) ; }
1932/*
1933  Feasible? Then we should have either a live node prepped for future
1934  expansion (indicated by variable() >= 0), or (miracle of miracles) an
1935  integral solution at the root node.
1936
1937  initializeInfo sets the reference counts in the nodeInfo object.  Since
1938  this node is still live, push it onto the heap that holds the live set.
1939*/
1940  double bestValue = 0.0 ;
1941  if (newNode) {
1942    bestValue = newNode->objectiveValue();
1943    if (newNode->branchingObject()) {
1944      newNode->initializeInfo() ;
1945      tree_->push(newNode) ;
1946      if (statistics_) {
1947        if (numberNodes2_==maximumStatistics_) {
1948          maximumStatistics_ = 2*maximumStatistics_;
1949          CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
1950          memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
1951          memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
1952          delete [] statistics_;
1953          statistics_=temp;
1954        }
1955        assert (!statistics_[numberNodes2_]);
1956        statistics_[numberNodes2_]=new CbcStatistics(newNode,this);
1957      }
1958      numberNodes2_++;
1959#     ifdef CHECK_NODE
1960      printf("Node %x on tree\n",newNode) ;
1961#     endif
1962    } else {
1963      // continuous is integer
1964      double objectiveValue = newNode->objectiveValue();
1965      setBestSolution(CBC_SOLUTION,objectiveValue,
1966                      solver_->getColSolution()) ;
1967      delete newNode ;
1968      newNode = NULL ;
1969    }
1970  }
1971
1972  if (printFrequency_ <= 0) {
1973    printFrequency_ = 1000 ;
1974    if (getNumCols() > 2000)
1975      printFrequency_ = 100 ;
1976  }
1977  /*
1978    It is possible that strong branching fixes one variable and then the code goes round
1979    again and again.  This can take too long.  So we need to warn user - just once.
1980  */
1981  numberLongStrong_=0;
1982  double totalTime = 0.0;
1983  CbcNode * createdNode=NULL;
1984#ifdef CBC_THREAD
1985  CbcModel ** threadModel = NULL;
1986  Coin_pthread_t * threadId = NULL;
1987  int * threadCount = NULL;
1988  pthread_mutex_t mutex;
1989  pthread_cond_t condition_main;
1990  pthread_mutex_t condition_mutex;
1991  pthread_mutex_t * mutex2 = NULL;
1992  pthread_cond_t * condition2 = NULL;
1993  threadStruct * threadInfo = NULL;
1994#ifdef CBC_NORMAL_THREAD
1995  bool locked=false;
1996#endif
1997  int threadStats[6];
1998#ifdef CBC_DETERMINISTIC_THREAD
1999  int defaultParallelIterations=500;
2000  int defaultParallelNodes=10;
2001#endif
2002  memset(threadStats,0,sizeof(threadStats));
2003  double timeWaiting=0.0;
2004  // For now just one model
2005  if (numberThreads_) {
2006    nodeCompare_->sayThreaded(); // need to use addresses
2007    threadId = new Coin_pthread_t [numberThreads_];
2008    threadCount = new int [numberThreads_];
2009    CoinZeroN(threadCount,numberThreads_);
2010    pthread_mutex_init(&mutex,NULL);
2011    pthread_cond_init(&condition_main,NULL);
2012    pthread_mutex_init(&condition_mutex,NULL);
2013    threadModel = new CbcModel * [numberThreads_+1];
2014    threadInfo = new threadStruct [numberThreads_+1];
2015    mutex2 = new pthread_mutex_t [numberThreads_];
2016    condition2 = new pthread_cond_t [numberThreads_];
2017#ifdef CBC_DETERMINISTIC_THREAD
2018    // May need for deterministic
2019    saveObjects=new OsiObject * [numberObjects_];
2020    for (int i=0;i<numberObjects_;i++) {
2021      saveObjects[i] = object_[i]->clone();
2022    }
2023#endif
2024    // we don't want a strategy object
2025    CbcStrategy * saveStrategy = strategy_;
2026    strategy_ = NULL;
2027    for (int i=0;i<numberThreads_;i++) {
2028      pthread_mutex_init(mutex2+i,NULL);
2029      pthread_cond_init(condition2+i,NULL);
2030      threadId[i].status=0;
2031      threadInfo[i].baseModel=this;
2032      threadModel[i]=new CbcModel(*this,true);
2033#ifdef COIN_HAS_CLP
2034      // Solver may need to know about model
2035      CbcModel * thisModel = threadModel[i];
2036      CbcOsiSolver * solver =
2037              dynamic_cast<CbcOsiSolver *>(thisModel->solver()) ;
2038      if (solver)
2039        solver->setCbcModel(thisModel);
2040#endif
2041      mutex_ = (void *) (threadInfo+i);
2042      threadModel[i]->moveToModel(this,-1);
2043      threadInfo[i].thisModel=threadModel[i];
2044      threadInfo[i].node=NULL;
2045      threadInfo[i].createdNode=NULL;
2046      threadInfo[i].threadIdOfBase.thr=pthread_self();
2047      threadInfo[i].mutex=&mutex;
2048      threadInfo[i].mutex2=mutex2+i;
2049      threadInfo[i].condition2=condition2+i;
2050      threadInfo[i].returnCode=-1;
2051      threadInfo[i].timeLocked=0.0;
2052      threadInfo[i].timeWaitingToLock=0.0;
2053      threadInfo[i].timeWaitingToStart=0.0;
2054      threadInfo[i].timeInThread=0.0;
2055      threadInfo[i].numberTimesLocked=0;
2056      threadInfo[i].numberTimesUnlocked=0;
2057      threadInfo[i].numberTimesWaitingToStart=0;
2058      threadInfo[i].dantzigState=0; // 0 unset, -1 waiting to be set, 1 set
2059      threadInfo[i].locked=false;
2060#if CBC_THREAD_DEBUG
2061      threadInfo[i].threadNumber=i+2;
2062#endif
2063#ifdef CBC_DETERMINISTIC_THREAD
2064      threadInfo[i].delNode = NULL;
2065      threadInfo[i].maxDeleteNode=0;
2066      threadInfo[i].nDeleteNode=0;
2067      threadInfo[i].nodesThisTime=0;
2068      threadInfo[i].iterationsThisTime=0;
2069#endif
2070      pthread_create(&(threadId[i].thr),NULL,doNodesThread,threadInfo+i);
2071      threadId[i].status = 1;
2072    }
2073    strategy_ = saveStrategy;
2074    // Do a partial one for base model
2075    threadInfo[numberThreads_].baseModel=this;
2076    threadModel[numberThreads_]=this;
2077    mutex_ = (void *) (threadInfo+numberThreads_);
2078    threadInfo[numberThreads_].node=NULL;
2079    threadInfo[numberThreads_].mutex=&mutex;
2080    threadInfo[numberThreads_].condition2=&condition_main;
2081    threadInfo[numberThreads_].mutex2=&condition_mutex;
2082    threadInfo[numberThreads_].timeLocked=0.0;
2083    threadInfo[numberThreads_].timeWaitingToLock=0.0;
2084    threadInfo[numberThreads_].numberTimesLocked=0;
2085    threadInfo[numberThreads_].numberTimesUnlocked=0;
2086    threadInfo[numberThreads_].locked=false;
2087#if CBC_THREAD_DEBUG
2088    threadInfo[numberThreads_].threadNumber=1;
2089#endif
2090  }
2091#endif
2092#ifdef COIN_HAS_CLP
2093  {
2094    OsiClpSolverInterface * clpSolver
2095      = dynamic_cast<OsiClpSolverInterface *> (solver_);
2096    if (clpSolver&&!parentModel_) {
2097      clpSolver->computeLargestAway();
2098    }
2099  }
2100#endif
2101/*
2102  At last, the actual branch-and-cut search loop, which will iterate until
2103  the live set is empty or we hit some limit (integrality gap, time, node
2104  count, etc.). The overall flow is to rebuild a subproblem, reoptimise using
2105  solveWithCuts(), choose a branching pattern with chooseBranch(), and finally
2106  add the node to the live set.
2107
2108  The first action is to winnow the live set to remove nodes which are worse
2109  than the current objective cutoff.
2110*/
2111  if (solver_->getRowCutDebuggerAlways()) {
2112    OsiRowCutDebugger * debuggerX = const_cast<OsiRowCutDebugger *> (solver_->getRowCutDebuggerAlways());
2113    const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
2114    if (!debugger) {
2115      // infeasible!!
2116      printf("before search\n");
2117      const double * lower = solver_->getColLower();
2118      const double * upper = solver_->getColUpper();
2119      const double * solution = debuggerX->optimalSolution();
2120      int numberColumns = solver_->getNumCols();
2121      for (int i=0;i<numberColumns;i++) {
2122        if (solver_->isInteger(i)) {
2123          if (solution[i]<lower[i]-1.0e-6||solution[i]>upper[i]+1.0e-6)
2124            printf("**** ");
2125          printf("%d %g <= %g <= %g\n",
2126                 i,lower[i],solution[i],upper[i]);
2127        }
2128      }
2129      //abort();
2130    }
2131  }
2132#ifdef CBC_DETERMINISTIC_THREAD
2133#define MAX_DEL_NODE 1
2134  CbcNode * delNode[MAX_DEL_NODE+1];
2135  int nDeleteNode=0;
2136  bool goneParallel=false;
2137#endif
2138  // For Printing etc when parallel
2139  int lastEvery1000=0;
2140  int lastPrintEvery=0;
2141  while (true) {
2142#ifdef CBC_NORMAL_THREAD
2143    if (!locked) {
2144      lockThread();
2145      locked=true;
2146    }
2147#endif
2148#ifdef COIN_HAS_CLP
2149    // Possible change of pivot method
2150    if(!savePivotMethod&&!parentModel_) {
2151      OsiClpSolverInterface * clpSolver
2152        = dynamic_cast<OsiClpSolverInterface *> (solver_);
2153      if (clpSolver&&numberNodes_>=100&&numberNodes_<200) {
2154        if (numberIterations_<numberNodes_*20) {
2155          ClpSimplex * simplex = clpSolver->getModelPtr();
2156          ClpDualRowPivot * pivotMethod=simplex->dualRowPivot();
2157          ClpDualRowDantzig * pivot =
2158            dynamic_cast< ClpDualRowDantzig*>(pivotMethod);
2159          if (!pivot) {
2160            savePivotMethod = pivotMethod->clone(true);
2161            ClpDualRowDantzig dantzig;
2162            simplex->setDualRowPivotAlgorithm(dantzig);
2163#ifdef COIN_DEVELOP
2164            printf("%d node, %d iterations ->Dantzig\n",numberNodes_,
2165                   numberIterations_);
2166#endif
2167#ifdef CBC_THREAD
2168            for (int i=0;i<numberThreads_;i++) {
2169              threadInfo[i].dantzigState=-1;
2170            }
2171#endif
2172          }
2173        }
2174      }
2175    }
2176#endif
2177    if (tree_->empty()) {
2178#ifdef CBC_NORMAL_THREAD
2179      if (numberThreads_) {
2180#ifdef COIN_DEVELOP
2181        printf("empty\n");
2182#endif
2183        // may still be outstanding nodes
2184        int iThread;
2185        for (iThread=0;iThread<numberThreads_;iThread++) {
2186          if (threadId[iThread].status) {
2187            if (threadInfo[iThread].returnCode==0) 
2188              break;
2189          }
2190        }
2191        if (iThread<numberThreads_) {
2192#ifdef COIN_DEVELOP
2193          printf("waiting for thread %d code 0\n",iThread);
2194#endif
2195#ifndef CBC_DETERMINISTIC_THREAD
2196          unlockThread();
2197#endif
2198          locked = false;
2199          pthread_cond_signal(threadInfo[iThread].condition2); // unlock in case
2200          while (true) {
2201            pthread_mutex_lock(&condition_mutex);
2202            struct timespec absTime;
2203            clock_gettime(CLOCK_REALTIME,&absTime);
2204            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
2205            absTime.tv_nsec += 1000000; // millisecond
2206            if (absTime.tv_nsec>=1000000000) {
2207              absTime.tv_nsec -= 1000000000;
2208              absTime.tv_sec++;
2209            }
2210            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
2211            clock_gettime(CLOCK_REALTIME,&absTime);
2212            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
2213            timeWaiting += time2-time;
2214            pthread_mutex_unlock(&condition_mutex);
2215            if (threadInfo[iThread].returnCode!=0) 
2216              break;
2217            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
2218          }
2219          threadModel[iThread]->moveToModel(this,1);
2220          assert (threadInfo[iThread].returnCode==1);
2221          if (threadInfo[iThread].dantzigState==-1) {
2222            // 0 unset, -1 waiting to be set, 1 set
2223            threadInfo[iThread].dantzigState=1;
2224            CbcModel * model = threadInfo[iThread].thisModel;
2225            OsiClpSolverInterface * clpSolver2
2226              = dynamic_cast<OsiClpSolverInterface *> (model->solver());
2227            assert (clpSolver2);
2228            ClpSimplex * simplex2 = clpSolver2->getModelPtr();
2229            ClpDualRowDantzig dantzig;
2230            simplex2->setDualRowPivotAlgorithm(dantzig);
2231          }
2232          // say available
2233          threadInfo[iThread].returnCode=-1;
2234          threadStats[4]++;
2235#ifdef COIN_DEVELOP
2236          printf("thread %d code now -1\n",iThread);
2237#endif
2238          continue;
2239        } else {
2240#ifdef COIN_DEVELOP
2241          printf("no threads at code 0 \n");
2242#endif
2243          // now check if any have just finished
2244          for (iThread=0;iThread<numberThreads_;iThread++) {
2245            if (threadId[iThread].status) {
2246              if (threadInfo[iThread].returnCode==1) 
2247                break;
2248            }
2249          }
2250          if (iThread<numberThreads_) {
2251#ifndef CBC_DETERMINISTIC_THREAD
2252            unlockThread();
2253#endif
2254            locked = false;
2255            threadModel[iThread]->moveToModel(this,1);
2256            assert (threadInfo[iThread].returnCode==1);
2257            // say available
2258            threadInfo[iThread].returnCode=-1;
2259            threadStats[4]++;
2260#ifdef COIN_DEVELOP
2261            printf("thread %d code now -1\n",iThread);
2262#endif
2263            continue;
2264          }
2265        }
2266        if (!tree_->empty()) {
2267#ifdef COIN_DEVELOP
2268          printf("tree not empty!!!!!!\n");
2269#endif
2270          continue;
2271        }
2272        for (iThread=0;iThread<numberThreads_;iThread++) {
2273          if (threadId[iThread].status) {
2274            if (threadInfo[iThread].returnCode!=-1) { 
2275              printf("bad end of tree\n");
2276              abort();
2277            }
2278          }
2279        }
2280#ifdef COIN_DEVELOP
2281        printf("finished ************\n");
2282#endif
2283      }
2284#ifndef CBC_DETERMINISTIC_THREAD
2285      unlockThread();
2286#endif
2287      locked=false; // not needed as break
2288#endif
2289      break;
2290    }
2291#ifdef CBC_NORMAL_THREAD
2292    unlockThread();
2293    locked = false;
2294#endif
2295    // If done 100 nodes see if worth trying reduction
2296    if (numberNodes_==100&&saveSolver) {
2297      bool tryNewSearch=solverCharacteristics_->reducedCostsAccurate();
2298      int numberColumns = getNumCols();
2299      if (tryNewSearch) {
2300        double cutoff = getCutoff() ;
2301        saveSolver->resolve();
2302        double direction = saveSolver->getObjSense() ;
2303        double gap = cutoff - saveSolver->getObjValue()*direction ;
2304        double tolerance;
2305        saveSolver->getDblParam(OsiDualTolerance,tolerance) ;
2306        if (gap<=0.0)
2307          gap = tolerance; 
2308        gap += 100.0*tolerance;
2309        double integerTolerance = getDblParam(CbcIntegerTolerance) ;
2310       
2311        const double *lower = saveSolver->getColLower() ;
2312        const double *upper = saveSolver->getColUpper() ;
2313        const double *solution = saveSolver->getColSolution() ;
2314        const double *reducedCost = saveSolver->getReducedCost() ;
2315       
2316        int numberFixed = 0 ;
2317        int numberFixed2=0;
2318#ifdef COIN_DEVELOP
2319        printf("gap %g\n",gap);
2320#endif
2321        for (int i = 0 ; i < numberIntegers_ ; i++) {
2322          int iColumn = integerVariable_[i] ;
2323          double djValue = direction*reducedCost[iColumn] ;
2324          if (upper[iColumn]-lower[iColumn] > integerTolerance) {
2325            if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap) {
2326              //printf("%d to lb on dj of %g - bounds %g %g\n",
2327              //     iColumn,djValue,lower[iColumn],upper[iColumn]);
2328              saveSolver->setColUpper(iColumn,lower[iColumn]) ;
2329              numberFixed++ ;
2330            } else if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap) {
2331              //printf("%d to ub on dj of %g - bounds %g %g\n",
2332              //     iColumn,djValue,lower[iColumn],upper[iColumn]);
2333              saveSolver->setColLower(iColumn,upper[iColumn]) ;
2334              numberFixed++ ;
2335            }
2336          } else {
2337            //printf("%d has dj of %g - already fixed to %g\n",
2338            //     iColumn,djValue,lower[iColumn]);
2339            numberFixed2++;
2340          }
2341        }
2342#ifdef COIN_DEVELOP
2343        if ((specialOptions_&1)!=0) {
2344          const OsiRowCutDebugger *debugger = saveSolver->getRowCutDebugger() ;
2345          if (debugger) { 
2346            printf("Contains optimal\n") ;
2347            saveSolver->writeMps("reduced");
2348          } else {
2349            abort();
2350          }
2351        }
2352        printf("Restart could fix %d integers (%d already fixed)\n",
2353               numberFixed+numberFixed2,numberFixed2);
2354#endif
2355        numberFixed += numberFixed2;
2356        if (numberFixed*10<numberColumns)
2357          tryNewSearch=false;
2358      }
2359      if (tryNewSearch) {
2360        // back to solver without cuts?
2361#if 0
2362        OsiSolverInterface * solver2 = continuousSolver_->clone();
2363#else
2364        OsiSolverInterface * solver2 = saveSolver->clone();
2365#endif
2366        const double *lower = saveSolver->getColLower() ;
2367        const double *upper = saveSolver->getColUpper() ;
2368        for (int i = 0 ; i < numberIntegers_ ; i++) {
2369          int iColumn = integerVariable_[i] ;
2370          solver2->setColLower(iColumn,lower[iColumn]);
2371          solver2->setColUpper(iColumn,upper[iColumn]);
2372        }
2373        // swap
2374        delete saveSolver;
2375        saveSolver=solver2;
2376        double * newSolution = new double[numberColumns];
2377        double objectiveValue=cutoff;
2378        CbcSerendipity heuristic(*this);
2379        if (bestSolution_)
2380          heuristic.setInputSolution(bestSolution_,bestObjective_);
2381        heuristic.setFractionSmall(0.6);
2382        heuristic.setFeasibilityPumpOptions(1008013);
2383        // Use numberNodes to say how many are original rows
2384        heuristic.setNumberNodes(continuousSolver_->getNumRows());
2385#ifdef COIN_DEVELOP
2386        if (continuousSolver_->getNumRows()<
2387            solver_->getNumRows())
2388          printf("%d rows added ZZZZZ\n",
2389                 solver_->getNumRows()-continuousSolver_->getNumRows());
2390#endif
2391        int returnCode= heuristic.smallBranchAndBound(saveSolver,
2392                                                      -1,newSolution,
2393                                                      objectiveValue,
2394                                                      cutoff,"Reduce");
2395        if (returnCode==-1) {
2396#ifdef COIN_DEVELOP
2397          printf("Restart - not small enough to do search after fixing\n");
2398#endif
2399          delete [] newSolution;
2400        } else {
2401          assert (returnCode>=0);
2402          if ((returnCode&1)!=0) {
2403            // increment number of solutions so other heuristics can test
2404            numberSolutions_++;
2405            numberHeuristicSolutions_++;
2406            lastHeuristic_ = NULL;
2407            setBestSolution(CBC_ROUNDING,objectiveValue,newSolution) ;
2408          }
2409          delete [] newSolution;
2410          break;
2411        }
2412      } 
2413      delete saveSolver;
2414      saveSolver=NULL;
2415    }
2416/*
2417  Check for abort on limits: node count, solution count, time, integrality gap.
2418*/
2419    totalTime = getCurrentSeconds() ;
2420    double maxSeconds = getMaximumSeconds();
2421    if (parentModel_)
2422      maxSeconds=CoinMin(maxSeconds,parentModel_->getMaximumSeconds());
2423    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
2424          numberSolutions_ < intParam_[CbcMaxNumSol] &&
2425          totalTime < maxSeconds &&
2426          !stoppedOnGap_&&!eventHappened_&&(maximumNumberIterations_<0||
2427                                            numberIterations_<maximumNumberIterations_))) {
2428      // out of loop
2429      break;
2430    }
2431#ifdef BONMIN
2432    assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
2433#endif
2434    if (cutoff > getCutoff()) {
2435      double newCutoff = getCutoff();
2436      if (analyzeResults_) {
2437        // see if we could fix any (more)
2438        int n=0;
2439        double * newLower = analyzeResults_;
2440        double * objLower = newLower+numberIntegers_;
2441        double * newUpper = objLower+numberIntegers_;
2442        double * objUpper = newUpper+numberIntegers_;
2443        for (int i=0;i<numberIntegers_;i++) {
2444          if (objLower[i]>newCutoff) {
2445            n++;
2446            if (objUpper[i]>newCutoff) {
2447              newCutoff = -COIN_DBL_MAX;
2448              break;
2449            }
2450          } else if (objUpper[i]>newCutoff) {
2451            n++;
2452          }
2453        }
2454        if (newCutoff==-COIN_DBL_MAX) {
2455          printf("Root analysis says finished\n");
2456        } else if (n>numberFixedNow_) {
2457          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
2458          numberFixedNow_=n;
2459        }
2460      }
2461      if (eventHandler) {
2462        if (!eventHandler->event(CbcEventHandler::solution)) {
2463          eventHappened_=true; // exit
2464        }
2465      }
2466#ifndef CBC_DETERMINISTIC_THREAD
2467      lockThread();
2468#endif
2469      // Do from deepest
2470      tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ;
2471      nodeCompare_->newSolution(this) ;
2472      nodeCompare_->newSolution(this,continuousObjective_,
2473                                continuousInfeasibilities_) ;
2474      tree_->setComparison(*nodeCompare_) ;
2475      if (tree_->empty()) {
2476#ifndef CBC_DETERMINISTIC_THREAD
2477        unlockThread();
2478#endif
2479        // For threads we need to check further
2480        //break; // finished
2481        continue;
2482      }
2483#ifndef CBC_DETERMINISTIC_THREAD
2484      unlockThread();
2485#endif
2486    }
2487    cutoff = getCutoff() ;
2488/*
2489    Periodic activities: Opportunities to
2490    + tweak the nodeCompare criteria,
2491    + check if we've closed the integrality gap enough to quit,
2492    + print a summary line to let the user know we're working
2493*/
2494    if (numberNodes_>=lastEvery1000) {
2495#ifndef CBC_DETERMINISTIC_THREAD
2496      lockThread();
2497#endif
2498#ifdef COIN_HAS_CLP
2499      // Possible change of pivot method
2500      if(!savePivotMethod&&!parentModel_) {
2501        OsiClpSolverInterface * clpSolver
2502          = dynamic_cast<OsiClpSolverInterface *> (solver_);
2503        if (clpSolver&&numberNodes_>=1000&&numberNodes_<2000) {
2504          if (numberIterations_<numberNodes_*20) {
2505            ClpSimplex * simplex = clpSolver->getModelPtr();
2506            ClpDualRowPivot * pivotMethod=simplex->dualRowPivot();
2507            ClpDualRowDantzig * pivot =
2508              dynamic_cast< ClpDualRowDantzig*>(pivotMethod);
2509            if (!pivot) {
2510              savePivotMethod = pivotMethod->clone(true);
2511              ClpDualRowDantzig dantzig;
2512              simplex->setDualRowPivotAlgorithm(dantzig);
2513#ifdef COIN_DEVELOP
2514              printf("%d node, %d iterations ->Dantzig\n",numberNodes_,
2515                     numberIterations_);
2516#endif
2517#ifdef CBC_THREAD
2518              for (int i=0;i<numberThreads_;i++) {
2519                threadInfo[i].dantzigState=-1;
2520              }
2521#endif
2522            }
2523          }
2524        }
2525      }
2526#endif
2527      lastEvery1000 = numberNodes_ + 1000;
2528      bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ;
2529#ifdef CHECK_CUT_SIZE
2530      verifyCutSize (tree_, *this);
2531#endif
2532      // redo tree if wanted
2533      if (redoTree)
2534        tree_->setComparison(*nodeCompare_) ;
2535#ifndef CBC_DETERMINISTIC_THREAD
2536      unlockThread();
2537#endif
2538    }
2539    if (saveCompare&&!hotstartSolution_) {
2540      // hotstart switched off
2541      delete nodeCompare_; // off depth first
2542      nodeCompare_=saveCompare;
2543      saveCompare=NULL;
2544      // redo tree
2545#ifndef CBC_DETERMINISTIC_THREAD
2546      lockThread();
2547#endif
2548      tree_->setComparison(*nodeCompare_) ;
2549#ifndef CBC_DETERMINISTIC_THREAD
2550      unlockThread();
2551#endif
2552    }
2553    if (numberNodes_>=lastPrintEvery) {
2554      lastPrintEvery = numberNodes_ + printFrequency_;
2555#ifdef CBC_INSTRUMENT
2556      if (0) {
2557        printf("==Start instrument\n");
2558        for (int iObject=0;iObject<numberObjects_;iObject++) {
2559          CbcSimpleIntegerDynamicPseudoCost * obj =
2560            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[iObject]) ;
2561          if (obj)
2562            obj->print();
2563        }
2564        printf("==End instrument\n");
2565      }
2566#endif
2567#ifndef CBC_DETERMINISTIC_THREAD
2568      lockThread();
2569#endif
2570      int nNodes = tree_->size() ;
2571
2572      //MODIF PIERRE
2573      bestPossibleObjective_ = tree_->getBestPossibleObjective();
2574#ifndef CBC_DETERMINISTIC_THREAD
2575      unlockThread();
2576#endif
2577      if (!intParam_[CbcPrinting]) {
2578        messageHandler()->message(CBC_STATUS,messages())
2579          << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
2580          <<getCurrentSeconds()
2581          << CoinMessageEol ;
2582      } else if (intParam_[CbcPrinting]==1) {
2583        messageHandler()->message(CBC_STATUS2,messages())
2584          << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
2585          <<lastDepth<<lastUnsatisfied<<numberIterations_
2586          <<getCurrentSeconds()
2587          << CoinMessageEol ;
2588      } else if (!numberExtraIterations_) {
2589        messageHandler()->message(CBC_STATUS2,messages())
2590          << numberNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
2591          <<lastDepth<<lastUnsatisfied<<numberIterations_
2592          <<getCurrentSeconds()
2593          << CoinMessageEol ;
2594      } else {
2595        messageHandler()->message(CBC_STATUS3,messages())
2596          << numberNodes_<<numberExtraNodes_<< nNodes<< bestObjective_<< bestPossibleObjective_
2597          <<lastDepth<<lastUnsatisfied<<numberIterations_<<numberExtraIterations_
2598          <<getCurrentSeconds()
2599          << CoinMessageEol ;
2600      }
2601      if (!eventHandler->event(CbcEventHandler::treeStatus)) {
2602        eventHappened_=true; // exit
2603      }
2604    }
2605    // See if can stop on gap
2606    double testGap = CoinMax(dblParam_[CbcAllowableGap],
2607                             CoinMax(fabs(bestObjective_),fabs(bestPossibleObjective_))
2608                             *dblParam_[CbcAllowableFractionGap]);
2609    if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
2610      stoppedOnGap_ = true ;
2611    }
2612   
2613#ifdef CHECK_NODE_FULL
2614    verifyTreeNodes(tree_,*this) ;
2615#   endif
2616#   ifdef CHECK_CUT_COUNTS
2617    verifyCutCounts(tree_,*this) ;
2618#   endif
2619/*
2620  Now we come to the meat of the loop. To create the active subproblem, we'll
2621  pop the most promising node in the live set, rebuild the subproblem it
2622  represents, and then execute the current arm of the branch to create the
2623  active subproblem.
2624*/
2625#ifdef BACK_TO_OLD_WAY //old way without threads #ifndef CBC_THREAD
2626    CbcNode *node = tree_->bestNode(cutoff) ;
2627    // Possible one on tree worse than cutoff
2628    if (!node||node->objectiveValue()>cutoff)
2629      continue;
2630    int currentNumberCuts = 0 ;
2631    currentNode_=node; // so can be accessed elsewhere
2632#ifdef CBC_DEBUG
2633    printf("%d unsat, way %d, obj %g est %g\n",
2634           node->numberUnsatisfied(),node->way(),node->objectiveValue(),
2635           node->guessedObjectiveValue());
2636#endif
2637#if NEW_UPDATE_OBJECT==0
2638    // Save clone in branching decision
2639    if(branchingMethod_)
2640      branchingMethod_->saveBranchingObject(node->modifiableBranchingObject());
2641#endif
2642    // Say not on optimal path
2643    bool onOptimalPath=false;
2644#   ifdef CHECK_NODE
2645    printf("Node %x popped from tree - %d left, %d count\n",node,
2646           node->nodeInfo()->numberBranchesLeft(),
2647           node->nodeInfo()->numberPointingToThis()) ;
2648    printf("\tdepth = %d, z =  %g, unsat = %d, var = %d.\n",
2649           node->depth(),node->objectiveValue(),
2650           node->numberUnsatisfied(),
2651           node->columnNumber()) ;
2652#   endif
2653    lastDepth=node->depth();
2654    lastUnsatisfied=node->numberUnsatisfied();
2655
2656/*
2657  Rebuild the subproblem for this node:  Call addCuts() to adjust the model
2658  to recreate the subproblem for this node (set proper variable bounds, add
2659  cuts, create a basis).  This may result in the problem being fathomed by
2660  bound or infeasibility. Returns 1 if node is fathomed.
2661  Execute the current arm of the branch: If the problem survives, save the
2662  resulting variable bounds and call branch() to modify variable bounds
2663  according to the current arm of the branching object. If we're processing
2664  the final arm of the branching object, flag the node for removal from the
2665  live set.
2666*/
2667    CbcNodeInfo * nodeInfo = node->nodeInfo() ;
2668    newNode = NULL ;
2669    int branchesLeft=0;
2670    if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_))
2671    { int i ;
2672      const double * lower = getColLower() ;
2673      const double * upper = getColUpper() ;
2674      for (i = 0 ; i < numberColumns ; i++)
2675      { lowerBefore[i]= lower[i] ;
2676        upperBefore[i]= upper[i] ; }
2677      if ((solverCharacteristics_->extraCharacteristics()&2)!=0) {
2678        solverCharacteristics_->setBeforeLower(lowerBefore);
2679        solverCharacteristics_->setBeforeUpper(upperBefore);
2680      }
2681      if (messageHandler()->logLevel()>2)
2682        node->modifiableBranchingObject()->print();
2683      if (!useOsiBranching) 
2684        branchesLeft = node->branch(NULL); // old way
2685      else
2686        branchesLeft = node->branch(solver_); // new way
2687      if (branchesLeft) {
2688        // set nodenumber correctly
2689        node->nodeInfo()->setNodeNumber(numberNodes2_);
2690        tree_->push(node) ;
2691        if (statistics_) {
2692          if (numberNodes2_==maximumStatistics_) {
2693            maximumStatistics_ = 2*maximumStatistics_;
2694            CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
2695            memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
2696            memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
2697            delete [] statistics_;
2698            statistics_=temp;
2699          }
2700          assert (!statistics_[numberNodes2_]);
2701          statistics_[numberNodes2_]=new CbcStatistics(node,this);
2702        }
2703        numberNodes2_++;
2704        //nodeOnTree=true; // back on tree
2705        //deleteNode = false ;
2706#       ifdef CHECK_NODE
2707        printf("Node %x pushed back on tree - %d left, %d count\n",node,
2708               nodeInfo->numberBranchesLeft(),
2709               nodeInfo->numberPointingToThis()) ;
2710#       endif
2711      } else {
2712        //deleteNode = true ;
2713        if (!nodeInfo->numberBranchesLeft())
2714          nodeInfo->allBranchesGone(); // can clean up
2715      }
2716      if ((specialOptions_&1)!=0) {
2717        /*
2718          This doesn't work as intended --- getRowCutDebugger will return null
2719          unless the current feasible solution region includes the optimal solution
2720          that RowCutDebugger knows. There's no way to tell inactive from off the
2721          optimal path.
2722        */
2723        const OsiRowCutDebugger *debugger = solver_->getRowCutDebugger() ;
2724        if (debugger) {
2725          onOptimalPath=true;
2726          printf("On optimal path\n") ;
2727        }
2728      }
2729     
2730/*
2731  Reoptimize, possibly generating cuts and/or using heuristics to find
2732  solutions.  Cut reference counts are unaffected unless we lose feasibility,
2733  in which case solveWithCuts() will make the adjustment.
2734*/
2735      phase_=2;
2736      cuts = OsiCuts() ;
2737      currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_ ;
2738      int saveNumber = numberIterations_;
2739      if(solverCharacteristics_->solutionAddsCuts()) {
2740        int returnCode=resolve(node ? node->nodeInfo() : NULL,1);
2741        feasible = returnCode != 0;
2742        if (feasible) {
2743          int iObject ;
2744          int preferredWay ;
2745          int numberUnsatisfied = 0 ;
2746          memcpy(currentSolution_,solver_->getColSolution(),
2747                 numberColumns*sizeof(double)) ;
2748          // point to useful information
2749          OsiBranchingInformation usefulInfo=usefulInformation();
2750         
2751          for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
2752            double infeasibility =
2753              object_[iObject]->infeasibility(&usefulInfo,preferredWay) ;
2754            if (infeasibility ) numberUnsatisfied++ ;
2755          }
2756          if (returnCode>0) {
2757            if (numberUnsatisfied)   {
2758              feasible = solveWithCuts(cuts,maximumCutPasses_,node);
2759            } else {
2760              // may generate cuts and turn the solution
2761              //to an infeasible one
2762              feasible = solveWithCuts(cuts, 1,
2763                                       node);
2764#if 0
2765              currentNumberCuts_ = cuts.sizeRowCuts();
2766              if (currentNumberCuts_ >= maximumNumberCuts_) {
2767                maximumNumberCuts_ = currentNumberCuts;
2768                delete [] addedCuts_;
2769                addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
2770              }
2771#endif
2772            }
2773          }
2774          // check extra info on feasibility
2775          if (!solverCharacteristics_->mipFeasible()) {
2776            feasible = false;
2777            solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
2778          }
2779        }
2780      } else {
2781        // normal
2782        //int zzzzzz=0;
2783        //if (zzzzzz)
2784        //solver_->writeMps("before");
2785        feasible = solveWithCuts(cuts,maximumCutPasses_,node);
2786      }
2787      if ((specialOptions_&1)!=0&&onOptimalPath) {
2788        if (!solver_->getRowCutDebugger()) {
2789          if (solver_->getRowCutDebuggerAlways()->optimalValue()<
2790              getCutoff()-1.0e-5) {
2791            // dj fix did something???
2792            solver_->writeMpsNative("infeas2.mps",NULL,NULL,2);
2793            solver_->getRowCutDebuggerAlways()->printOptimalSolution(*solver_);
2794            assert (solver_->getRowCutDebugger()) ;
2795          }
2796        }
2797      }
2798      if (statistics_) {
2799        assert (numberNodes2_);
2800        assert (statistics_[numberNodes2_-1]);
2801        assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
2802        statistics_[numberNodes2_-1]->endOfBranch(numberIterations_-saveNumber,
2803                                               feasible ? solver_->getObjValue()
2804                                               : COIN_DBL_MAX);
2805      }
2806/*
2807  Are we still feasible? If so, create a node and do the work to attach a
2808  branching object, reoptimising as needed if chooseBranch() identifies
2809  monotone objects.
2810
2811  Finally, attach a partial nodeInfo object and store away any cuts that we
2812  created back in solveWithCuts. addCuts() will initialise the reference
2813  counts for these new cuts.
2814
2815  This next test can be problematic if we've discovered an
2816  alternate equivalent answer and subsequently fathom the solution
2817  known to the row cut debugger due to bounds.
2818*/
2819        if (onOptimalPath) {
2820          bool objLim = solver_->isDualObjectiveLimitReached() ;
2821          if (!feasible && !objLim) {
2822            printf("infeas2\n");
2823            solver_->writeMpsNative("infeas.mps",NULL,NULL,2);
2824            solver_->getRowCutDebuggerAlways()->printOptimalSolution(*solver_);
2825            CoinWarmStartBasis *slack =
2826              dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
2827            solver_->setWarmStart(slack);
2828            delete slack ;
2829            solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
2830            solver_->initialSolve();
2831            assert (!solver_->isProvenOptimal());
2832          }
2833          assert (feasible || objLim);
2834        }
2835        bool checkingNode=false;
2836        if (feasible) {
2837          newNode = new CbcNode ;//Regular node of the tree
2838          // Set objective value (not so obvious if NLP etc)
2839          setObjectiveValue(newNode,node);
2840          anyAction =-1 ;
2841          resolved = false ;
2842          if (newNode->objectiveValue() >= getCutoff()) 
2843            anyAction=-2;
2844          // only allow at most a few passes
2845          int numberPassesLeft=5;
2846          checkingNode=true;
2847        OsiSolverBranch * branches=NULL;
2848        // point to useful information
2849        anyAction = chooseBranch(newNode, numberPassesLeft,node, cuts,resolved,
2850                                 lastws, lowerBefore, upperBefore, branches);
2851/*
2852  If we end up infeasible, we can delete the new node immediately. Since this
2853  node won't be needing the cuts we collected, decrement the reference counts.
2854  If we are feasible, then we'll be placing this node into the live set, so
2855  increment the reference count in the current (parent) nodeInfo.
2856*/
2857        if (anyAction == -2)
2858          { delete newNode ;
2859          newNode = NULL ;
2860          // say strong doing well
2861          if (checkingNode)
2862            setSpecialOptions(specialOptions_|8);
2863          for (i = 0 ; i < currentNumberCuts_ ; i++)
2864            { if (addedCuts_[i])
2865              { if (!addedCuts_[i]->decrement(1))
2866                delete addedCuts_[i] ; } } }
2867        else
2868          { nodeInfo->increment() ;
2869          if ((numberNodes_%20)==0) {
2870            // say strong not doing as well
2871            setSpecialOptions(specialOptions_&~8);
2872          }
2873        }
2874        }
2875/*
2876  At this point, there are three possibilities:
2877    * newNode is live and will require further branching to resolve
2878      (variable() >= 0). Increment the cut reference counts by
2879      numberBranches() to allow for use by children of this node, and
2880      decrement by 1 because we've executed one arm of the branch of our
2881      parent (consuming one reference). Before we push newNode onto the
2882      search tree, try for a heuristic solution.
2883    * We have a solution, in which case newNode is non-null but we have no
2884      branching variable. Decrement the cut counts and save the solution.
2885    * The node was found to be infeasible, in which case it's already been
2886      deleted, and newNode is null.
2887*/
2888        if (!eventHandler->event(CbcEventHandler::node)) {
2889          eventHappened_=true; // exit
2890        }
2891        assert (!newNode || newNode->objectiveValue() <= getCutoff()) ;
2892        if (statistics_) {
2893          assert (numberNodes2_);
2894          assert (statistics_[numberNodes2_-1]);
2895          assert (statistics_[numberNodes2_-1]->node()==numberNodes2_-1);
2896          if (newNode)
2897            statistics_[numberNodes2_-1]->updateInfeasibility(newNode->numberUnsatisfied());
2898          else
2899            statistics_[numberNodes2_-1]->sayInfeasible();
2900        }
2901        if (newNode) {
2902          if (newNode->branchingObject() == NULL&&solverCharacteristics_->solverType()==4) {
2903            // need to check if any cuts would do anything
2904            OsiCuts theseCuts;
2905            // reset probing info
2906            //if (probingInfo_)
2907            //probingInfo_->initializeFixing();
2908            for (int i = 0;i<numberCutGenerators_;i++) {
2909              bool generate = generator_[i]->normal();
2910              // skip if not optimal and should be (maybe a cut generator has fixed variables)
2911              if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
2912                generate=false;
2913              if (!generator_[i]->mustCallAgain())
2914                generate=false; // only special cuts
2915              if (generate) {
2916                generator_[i]->generateCuts(theseCuts,1,solver_,NULL) ;
2917                int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
2918                if (numberRowCutsAfter) {
2919                  // need dummy branch
2920                  newNode->setBranchingObject(new CbcDummyBranchingObject(this));
2921                  newNode->nodeInfo()->initializeInfo(1);
2922                  break;
2923                }
2924              }
2925            }
2926          }
2927          if (newNode->branchingObject())
2928          { handler_->message(CBC_BRANCH,messages_)
2929               << numberNodes_<< newNode->objectiveValue()
2930               << newNode->numberUnsatisfied()<< newNode->depth()
2931               << CoinMessageEol ;
2932            // Increment cut counts (taking off current)
2933            int numberLeft = newNode->numberBranches() ;
2934            for (i = 0;i < currentNumberCuts_;i++)
2935            { if (addedCuts_[i])
2936              {
2937#               ifdef CHECK_CUT_COUNTS
2938                printf("Count on cut %x increased by %d\n",addedCuts_[i],
2939                        numberLeft-1) ;
2940#               endif
2941                addedCuts_[i]->increment(numberLeft-1) ; } }
2942
2943            double estValue = newNode->guessedObjectiveValue() ;
2944            int found = -1 ;
2945            // no - overhead on small problems solver_->resolve() ;     // double check current optimal
2946            // assert (!solver_->getIterationCount());
2947            double * newSolution = new double [numberColumns] ;
2948            double heurValue = getCutoff() ;
2949            int iHeur ;
2950            for (iHeur = 0 ; iHeur < numberHeuristics_ ; iHeur++) {
2951              double saveValue = heurValue ;
2952              int ifSol = heuristic_[iHeur]->solution(heurValue,newSolution) ;
2953              if (ifSol > 0) {
2954                // new solution found
2955                heuristic_[iHeur]->incrementNumberSolutionsFound();
2956                found = iHeur ;
2957                incrementUsed(newSolution);
2958                lastHeuristic_ = heuristic_[found];
2959                setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
2960              } else if (ifSol < 0) {
2961                // just returning an estimate
2962                estValue = CoinMin(heurValue,estValue) ;
2963                heurValue = saveValue ;
2964              }
2965            }
2966            delete [] newSolution ;
2967            newNode->setGuessedObjectiveValue(estValue) ;
2968            tree_->push(newNode) ;
2969            if (statistics_) {
2970              if (numberNodes2_==maximumStatistics_) {
2971                maximumStatistics_ = 2*maximumStatistics_;
2972                CbcStatistics ** temp = new CbcStatistics * [maximumStatistics_];
2973                memset(temp,0,maximumStatistics_*sizeof(CbcStatistics *));
2974                memcpy(temp,statistics_,numberNodes2_*sizeof(CbcStatistics *));
2975                delete [] statistics_;
2976                statistics_=temp;
2977              }
2978              assert (!statistics_[numberNodes2_]);
2979              statistics_[numberNodes2_]=new CbcStatistics(newNode,this);
2980            }
2981            numberNodes2_++;
2982#           ifdef CHECK_NODE
2983            printf("Node %x pushed on tree c\n",newNode) ;
2984#           endif
2985          }
2986          else
2987          { 
2988            if(solverCharacteristics_ && //we may be in a non standard bab
2989               solverCharacteristics_->solutionAddsCuts()// we are in some kind of OA based bab.
2990               )
2991              {
2992                std::cerr<<"You should never get here"<<std::endl;
2993                throw CoinError("Nodes should not be fathomed on integer infeasibility in this setting",
2994                                "branchAndBound","CbcModel") ;
2995              }
2996            for (i = 0 ; i < currentNumberCuts_ ; i++)
2997            { if (addedCuts_[i])
2998              { if (!addedCuts_[i]->decrement(1))
2999                  delete addedCuts_[i] ; } }
3000          double objectiveValue = newNode->objectiveValue();
3001            setBestSolution(CBC_SOLUTION,objectiveValue,
3002                            solver_->getColSolution()) ;
3003            lastHeuristic_ = NULL;
3004            incrementUsed(solver_->getColSolution());
3005            //assert(nodeInfo->numberPointingToThis() <= 2) ;
3006            // avoid accidental pruning, if newNode was final branch arm
3007            nodeInfo->increment();
3008            delete newNode ;
3009            nodeInfo->decrement() ; } }
3010/*
3011  This node has been completely expanded and can be removed from the live
3012  set.
3013*/
3014      if (branchesLeft)
3015      { 
3016      }
3017      else
3018      { 
3019        if (!nodeInfo->numberBranchesLeft())
3020          nodeInfo->allBranchesGone(); // can clean up
3021        delete node ; }
3022    } else {
3023      // add cuts found to be infeasible (on bound)!
3024      abort();
3025      delete node;
3026    }
3027/*
3028  Delete cuts to get back to the original system.
3029
3030  I'm thinking this is redundant --- the call to addCuts that conditions entry
3031  to this code block also performs this action.
3032*/
3033      int numberToDelete = getNumRows()-numberRowsAtContinuous_ ;
3034      if (numberToDelete)
3035      { int * delRows = new int[numberToDelete] ;
3036        int i ;
3037        for (i = 0 ; i < numberToDelete ; i++)
3038        { delRows[i] = i+numberRowsAtContinuous_ ; }
3039        solver_->deleteRows(numberToDelete,delRows) ;
3040        delete [] delRows ; }
3041#else // end of not CBC_THREAD
3042#ifndef CBC_DETERMINISTIC_THREAD
3043      CbcNode *node = tree_->bestNode(cutoff) ;
3044      // Possible one on tree worse than cutoff
3045      if (!node||node->objectiveValue()>cutoff) 
3046        continue;
3047    if (!numberThreads_) {
3048#else
3049      if (!numberThreads_||(tree_->size()<5*numberThreads_&&!goneParallel)) {
3050      CbcNode *node = tree_->bestNode(cutoff) ;
3051      // Possible one on tree worse than cutoff
3052      if (!node||node->objectiveValue()>cutoff)
3053        continue;
3054#endif
3055      doOneNode(this,node,createdNode);
3056#ifdef CBC_DETERMINISTIC_THREAD
3057        assert (createdNode);
3058        if (!createdNode->active()) {
3059          //if (createdNode->nodeInfo()) {
3060          //createdNode->nodeInfo()->throwAway();
3061          //}
3062          delete createdNode;
3063          createdNode=NULL;
3064        } else {
3065          // Say one more pointing to this
3066          node->nodeInfo()->increment() ;
3067          tree_->push(createdNode) ;
3068        }
3069        //if (node) {
3070        //assert (node->active());
3071        if (node->active()) {
3072          assert (node->nodeInfo());
3073          if (node->nodeInfo()->numberBranchesLeft()) {
3074            tree_->push(node) ;
3075          } else {
3076            node->setActive(false);
3077          }
3078        } else {
3079          if (node->nodeInfo()) {
3080            if (!node->nodeInfo()->numberBranchesLeft())
3081              node->nodeInfo()->allBranchesGone(); // can clean up
3082            // So will delete underlying stuff
3083            node->setActive(true);
3084          }
3085          delNode[nDeleteNode++]=node;
3086          node=NULL;
3087        } 
3088        if (nDeleteNode>=MAX_DEL_NODE) {
3089          for (int i=0;i<nDeleteNode;i++) {
3090            //printf("trying to del %d %x\n",i,delNode[i]);
3091            delete delNode[i];
3092            //printf("done to del %d %x\n",i,delNode[i]);
3093          }
3094          nDeleteNode=0;
3095        }
3096#endif
3097      } else {
3098#ifdef CBC_NORMAL_THREAD
3099        threadStats[0]++;
3100        //need to think
3101        int iThread;
3102        // Start one off if any available
3103        for (iThread=0;iThread<numberThreads_;iThread++) {
3104          if (threadInfo[iThread].returnCode==-1) 
3105            break;
3106        }
3107        if (iThread<numberThreads_) {
3108          threadInfo[iThread].node=node;
3109          assert (threadInfo[iThread].returnCode==-1);
3110          // say in use
3111          threadInfo[iThread].returnCode=0;
3112          threadModel[iThread]->moveToModel(this,0);
3113          pthread_cond_signal(threadInfo[iThread].condition2); // unlock
3114          threadCount[iThread]++;
3115        }
3116#ifndef CBC_DETERMINISTIC_THREAD
3117        lockThread();
3118#endif
3119        locked=true;
3120        // see if any finished
3121        for (iThread=0;iThread<numberThreads_;iThread++) {
3122          if (threadInfo[iThread].returnCode>0) 
3123            break;
3124        }
3125#ifndef CBC_DETERMINISTIC_THREAD
3126        unlockThread();
3127#endif
3128        locked=false;
3129        if (iThread<numberThreads_) {
3130          threadModel[iThread]->moveToModel(this,1);
3131          assert (threadInfo[iThread].returnCode==1);
3132          // say available
3133          threadInfo[iThread].returnCode=-1;
3134          // carry on
3135          threadStats[3]++;
3136        } else {
3137          // Start one off if any available
3138          for (iThread=0;iThread<numberThreads_;iThread++) {
3139            if (threadInfo[iThread].returnCode==-1) 
3140              break;
3141          }
3142          if (iThread<numberThreads_) {
3143#ifndef CBC_DETERMINISTIC_THREAD
3144            lockThread();
3145#endif
3146            locked=true;
3147            // If any on tree get
3148            if (!tree_->empty()) {
3149              //node = tree_->bestNode(cutoff) ;
3150              //assert (node);
3151              threadStats[1]++;
3152              continue; // ** get another node
3153            }
3154#ifndef CBC_DETERMINISTIC_THREAD
3155            unlockThread();
3156#endif
3157            locked=false;
3158          }
3159          // wait (for debug could sleep and use test)
3160          bool finished=false;
3161          while (!finished) {
3162            pthread_mutex_lock(&condition_mutex);
3163            struct timespec absTime;
3164            clock_gettime(CLOCK_REALTIME,&absTime);
3165            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
3166            absTime.tv_nsec += 1000000; // millisecond
3167            if (absTime.tv_nsec>=1000000000) {
3168              absTime.tv_nsec -= 1000000000;
3169              absTime.tv_sec++;
3170            }
3171            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
3172            clock_gettime(CLOCK_REALTIME,&absTime);
3173            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
3174            timeWaiting += time2-time;
3175            pthread_mutex_unlock(&condition_mutex);
3176            for (iThread=0;iThread<numberThreads_;iThread++) {
3177              if (threadInfo[iThread].returnCode>0) {
3178                finished=true;
3179                break;
3180              } else if (threadInfo[iThread].returnCode==0) {
3181                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
3182              }
3183            }
3184          }
3185          assert (iThread<numberThreads_);
3186          threadModel[iThread]->moveToModel(this,1);
3187          node = threadInfo[iThread].node;
3188          threadInfo[iThread].node=NULL;
3189          assert (threadInfo[iThread].returnCode==1);
3190          // say available
3191          threadInfo[iThread].returnCode=-1;
3192          // carry on
3193          threadStats[2]++;
3194        }
3195#else
3196        // Deterministic parallel
3197#ifndef CBC_DETERMINISTIC_THREAD
3198        abort();
3199#endif
3200#ifdef CBC_THREAD
3201        int saveTreeSize = tree_->size();
3202        goneParallel=true;
3203        int nAffected=splitModel(numberThreads_,threadModel,defaultParallelNodes);
3204        int saveTreeSize2 = tree_->size();
3205        int iThread;
3206        // do all until finished
3207        for (iThread=0;iThread<numberThreads_;iThread++) {
3208          // obviously tune
3209          threadInfo[iThread].nDeleteNode=defaultParallelIterations;
3210        }
3211        // Save current state
3212        int iObject;
3213        for (iObject=0;iObject<numberObjects_;iObject++) {
3214          saveObjects[iObject]->updateBefore(object_[iObject]);
3215        }
3216        for (iThread=0;iThread<numberThreads_;iThread++) {
3217          threadInfo[iThread].returnCode=0;
3218          pthread_cond_signal(threadInfo[iThread].condition2); // unlock
3219#if 0
3220          //wait!!
3221          bool finished=false;
3222          while (!finished) {
3223            pthread_mutex_lock(&condition_mutex);
3224            struct timespec absTime;
3225            clock_gettime(CLOCK_REALTIME,&absTime);
3226            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
3227            absTime.tv_nsec += 1000000; // millisecond
3228            if (absTime.tv_nsec>=1000000000) {
3229              absTime.tv_nsec -= 1000000000;
3230              absTime.tv_sec++;
3231            }
3232            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
3233            clock_gettime(CLOCK_REALTIME,&absTime);
3234            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
3235            timeWaiting += time2-time;
3236            pthread_mutex_unlock(&condition_mutex);
3237            finished=true;
3238            if (threadInfo[iThread].returnCode<=0) {
3239              finished=false;
3240            }
3241          }
3242#endif
3243        }
3244        // wait
3245        bool finished=false;
3246        while (!finished) {
3247          pthread_mutex_lock(&condition_mutex);
3248          struct timespec absTime;
3249          clock_gettime(CLOCK_REALTIME,&absTime);
3250          double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
3251          absTime.tv_nsec += 1000000; // millisecond
3252          if (absTime.tv_nsec>=1000000000) {
3253            absTime.tv_nsec -= 1000000000;
3254            absTime.tv_sec++;
3255          }
3256          pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
3257          clock_gettime(CLOCK_REALTIME,&absTime);
3258          double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
3259          timeWaiting += time2-time;
3260          pthread_mutex_unlock(&condition_mutex);
3261          finished=true;
3262          for (iThread=0;iThread<numberThreads_;iThread++) {
3263            if (threadInfo[iThread].returnCode<=0) {
3264              finished=false;
3265            }
3266          }
3267        }
3268        // Unmark marked
3269        for (int i=0;i<nAffected;i++) {
3270          walkback_[i]->unmark();
3271        }
3272        assert (saveTreeSize2 == tree_->size());
3273        if (0) { 
3274          // put back cut counts
3275          for (int i=0;i<nAffected;i++) {
3276            walkback_[i]->decrementCuts(1000000);
3277          }
3278        }
3279#ifndef NDEBUG
3280        for (iObject=0;iObject<numberObjects_;iObject++) {
3281          CbcSimpleIntegerDynamicPseudoCost * obj =
3282            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[iObject]) ;
3283          CbcSimpleIntegerDynamicPseudoCost * obj2 =
3284            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(saveObjects[iObject]) ;
3285          assert (obj->same(obj2));
3286        }
3287#endif
3288        int iModel;
3289        double scaleFactor=1.0;
3290        for (iModel=0;iModel<numberThreads_;iModel++) {
3291          //printf("model %d tree size %d\n",iModel,threadModel[iModel]->tree_->size());
3292          if (saveTreeSize>4*numberThreads_*defaultParallelNodes) {
3293            if (!threadModel[iModel]->tree_->size()) {
3294              scaleFactor *= 1.05;
3295            }
3296          }
3297          threadModel[iModel]->moveToModel(this,11);
3298          // Update base model
3299          OsiObject ** threadObject = threadModel[iModel]->object_;
3300          for (iObject=0;iObject<numberObjects_;iObject++) {
3301            object_[iObject]->updateAfter(threadObject[iObject],saveObjects[iObject]);
3302          }
3303        }
3304        if (scaleFactor!=1.0) {
3305          int newNumber = (int) (defaultParallelNodes * scaleFactor+0.5001);
3306          if (newNumber*2<defaultParallelIterations) {
3307            printf("Changing tree size from %d to %d\n",
3308                   defaultParallelNodes,newNumber);
3309            defaultParallelNodes = newNumber;
3310          }
3311        }
3312        printf("Tree sizes %d %d %d - affected %d\n",saveTreeSize,saveTreeSize2,tree_->size(),nAffected);
3313        // later remember random may not be thread neutral
3314#endif
3315#endif
3316      }
3317      //lastDepth=node->depth();
3318      //lastUnsatisfied=node->numberUnsatisfied();
3319#endif // end of CBC_THREAD
3320  }
3321#ifdef CBC_DETERMINISTIC_THREAD
3322  if (nDeleteNode) {
3323    for (int i=0;i<nDeleteNode;i++) {
3324      delete delNode[i];
3325    }
3326    nDeleteNode=0;
3327  }
3328#endif
3329#ifdef CBC_THREAD
3330  if (numberThreads_) {
3331    //printf("stats ");
3332    //for (unsigned int j=0;j<sizeof(threadStats)/sizeof(int);j++)
3333    //printf("%d ",threadStats[j]);
3334    //printf("\n");
3335    int i;
3336    // Seems to be bug in CoinCpu on Linux - does threads as well despite documentation
3337    double time=0.0;
3338    for (i=0;i<numberThreads_;i++) 
3339      time += threadInfo[i].timeInThread;
3340    bool goodTimer = time<(getCurrentSeconds());
3341    //bool stopped = (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
3342    //        numberSolutions_ < intParam_[CbcMaxNumSol] &&
3343    //        totalTime < dblParam_[CbcMaximumSeconds] &&
3344    //        !stoppedOnGap_&&!eventHappened_));
3345    for (i=0;i<numberThreads_;i++) {
3346      while (threadInfo[i].returnCode==0) {
3347        pthread_cond_signal(threadInfo[i].condition2); // unlock
3348        pthread_mutex_lock(&condition_mutex);
3349        struct timespec absTime;
3350        clock_gettime(CLOCK_REALTIME,&absTime);
3351        absTime.tv_nsec += 1000000; // millisecond
3352        if (absTime.tv_nsec>=1000000000) {
3353          absTime.tv_nsec -= 1000000000;
3354          absTime.tv_sec++;
3355        }
3356        pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
3357        clock_gettime(CLOCK_REALTIME,&absTime);
3358        pthread_mutex_unlock(&condition_mutex);
3359      }
3360      pthread_cond_signal(threadInfo[i].condition2); // unlock
3361      pthread_mutex_lock(&condition_mutex); // not sure necessary but have had one hang on interrupt
3362      threadModel[i]->numberThreads_=0; // say exit
3363#ifdef CBC_DETERMINISTIC_THREAD
3364      delete [] threadInfo[i].delNode;
3365#endif
3366      threadInfo[i].returnCode=0;
3367      pthread_mutex_unlock(&condition_mutex);
3368      pthread_cond_signal(threadInfo[i].condition2); // unlock
3369      //if (!stopped)
3370      //pthread_join(threadId[i],NULL);
3371      int returnCode;
3372      returnCode=pthread_join(threadId[i].thr,NULL);
3373      threadId[i].status = 0;
3374      assert (!returnCode);
3375        //else
3376        //pthread_kill(threadId[i]); // kill rather than try and synchronize
3377      threadModel[i]->moveToModel(this,2);
3378      pthread_mutex_destroy (threadInfo[i].mutex2);
3379      pthread_cond_destroy (threadInfo[i].condition2);
3380      assert (threadInfo[i].numberTimesLocked==threadInfo[i].numberTimesUnlocked);
3381      handler_->message(CBC_THREAD_STATS,messages_)
3382        <<"Thread";
3383      handler_->printing(true)
3384        <<i<<threadCount[i]<<threadInfo[i].timeWaitingToStart;
3385      handler_->printing(goodTimer)<<threadInfo[i].timeInThread;
3386      handler_->printing(false)<<0.0;
3387      handler_->printing(true)<<threadInfo[i].numberTimesLocked
3388        <<threadInfo[i].timeLocked<<threadInfo[i].timeWaitingToLock
3389        <<CoinMessageEol;
3390    }
3391    assert (threadInfo[numberThreads_].numberTimesLocked==threadInfo[numberThreads_].numberTimesUnlocked);
3392    handler_->message(CBC_THREAD_STATS,messages_)
3393      <<"Main thread";
3394    handler_->printing(false)<<0<<0<<0.0;
3395    handler_->printing(false)<<0.0;
3396    handler_->printing(true)<<timeWaiting;
3397    handler_->printing(true)<<threadInfo[numberThreads_].numberTimesLocked
3398      <<threadInfo[numberThreads_].timeLocked<<threadInfo[numberThreads_].timeWaitingToLock
3399      <<CoinMessageEol;
3400    pthread_mutex_destroy (&mutex);
3401    pthread_cond_destroy (&condition_main);
3402    pthread_mutex_destroy (&condition_mutex);
3403    // delete models (here in case some point to others)
3404    for (i=0;i<numberThreads_;i++) {
3405      // make sure handler will be deleted
3406      threadModel[i]->defaultHandler_=true;
3407      delete threadModel[i];
3408    }
3409    delete [] mutex2;
3410    delete [] condition2;
3411    delete [] threadId;
3412    delete [] threadInfo;
3413    delete [] threadModel;
3414    delete [] threadCount;
3415    mutex_=NULL;
3416    // adjust time to allow for children on some systems
3417    dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
3418  }
3419#endif
3420/*
3421  End of the non-abort actions. The next block of code is executed if we've
3422  aborted because we hit one of the limits. Clean up by deleting the live set
3423  and break out of the node processing loop. Note that on an abort, node may
3424  have been pushed back onto the tree for further processing, in which case
3425  it'll be deleted in cleanTree. We need to check.
3426*/
3427    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
3428        numberSolutions_ < intParam_[CbcMaxNumSol] &&
3429        totalTime < dblParam_[CbcMaximumSeconds] &&
3430        !stoppedOnGap_&&!eventHappened_)) {
3431      if (tree_->size()) {
3432        double dummyBest;
3433        tree_->cleanTree(this,-COIN_DBL_MAX,dummyBest) ;
3434      }
3435      delete nextRowCut_;
3436      if (stoppedOnGap_)
3437        { messageHandler()->message(CBC_GAP,messages())
3438          << bestObjective_-bestPossibleObjective_
3439          << dblParam_[CbcAllowableGap]
3440          << dblParam_[CbcAllowableFractionGap]*100.0
3441          << CoinMessageEol ;
3442        secondaryStatus_ = 2;
3443        status_ = 0 ; }
3444        else
3445          if (isNodeLimitReached())
3446            { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
3447            secondaryStatus_ = 3;
3448            status_ = 1 ; }
3449          else
3450        if (totalTime >= dblParam_[CbcMaximumSeconds])
3451          { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
3452          secondaryStatus_ = 4;
3453          status_ = 1 ; }
3454        else
3455          if (eventHappened_)
3456            { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
3457            secondaryStatus_ = 5;
3458            status_ = 5 ; }
3459          else
3460            { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
3461            secondaryStatus_ = 6;
3462            status_ = 1 ; }
3463    }
3464/*
3465  That's it, we've exhausted the search tree, or broken out of the loop because
3466  we hit some limit on evaluation.
3467
3468  We may have got an intelligent tree so give it one more chance
3469*/
3470  // Tell solver we are not in Branch and Cut
3471  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
3472  tree_->endSearch();
3473  //  If we did any sub trees - did we give up on any?
3474  if ( numberStoppedSubTrees_)
3475    status_=1;
3476  if (!status_) {
3477    // Set best possible unless stopped on gap
3478    if(secondaryStatus_ != 2)
3479      bestPossibleObjective_=bestObjective_;
3480    handler_->message(CBC_END_GOOD,messages_)
3481      << bestObjective_ << numberIterations_ << numberNodes_<<getCurrentSeconds()
3482      << CoinMessageEol ;
3483  } else {
3484    handler_->message(CBC_END,messages_)
3485      << bestObjective_ <<bestPossibleObjective_
3486      << numberIterations_ << numberNodes_<<getCurrentSeconds()
3487      << CoinMessageEol ;
3488  }
3489  if (numberStrongIterations_)
3490    handler_->message(CBC_STRONG_STATS,messages_)
3491      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
3492      << strongInfo_[1] << CoinMessageEol ;
3493  if (!numberExtraNodes_) 
3494    handler_->message(CBC_OTHER_STATS,messages_)
3495      << maximumDepthActual_
3496      << numberDJFixed_ << CoinMessageEol ;
3497  else
3498    handler_->message(CBC_OTHER_STATS2,messages_)
3499      << maximumDepthActual_
3500      << numberDJFixed_ << numberExtraNodes_<<numberExtraIterations_
3501      <<CoinMessageEol ;
3502  if (doStatistics==100) {
3503    for (int i=0;i<numberObjects_;i++) {
3504      CbcSimpleIntegerDynamicPseudoCost * obj =
3505        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
3506      if (obj)
3507        obj->print();
3508    }
3509  }
3510  if (statistics_) {
3511    // report in some way
3512    int * lookup = new int[numberObjects_];
3513    int i;
3514    for (i=0;i<numberObjects_;i++) 
3515      lookup[i]=-1;
3516    bool goodIds=true;
3517    for (i=0;i<numberObjects_;i++) {
3518      int iColumn = object_[i]->columnNumber();
3519      if(iColumn>=0&&iColumn<numberColumns) {
3520        if (lookup[i]==-1) {
3521          lookup[i]=iColumn;
3522        } else {
3523          goodIds=false;
3524          break;
3525        }
3526      } else {
3527        goodIds=false;
3528        break;
3529      }
3530    }
3531    if (!goodIds) {
3532      delete [] lookup;
3533      lookup=NULL;
3534    }
3535    if (doStatistics==3) {
3536      printf("  node parent depth column   value                    obj      inf\n");
3537      for ( i=0;i<numberNodes2_;i++) {
3538        statistics_[i]->print(lookup);
3539      }
3540    }
3541    if (doStatistics>1) {
3542      // Find last solution
3543      int k;
3544      for (k=numberNodes2_-1;k>=0;k--) {
3545        if (statistics_[k]->endingObjective()!=COIN_DBL_MAX&&
3546            !statistics_[k]->endingInfeasibility())
3547          break;
3548      }
3549      if (k>=0) {
3550        int depth=statistics_[k]->depth();
3551        int * which = new int[depth+1];
3552        for (i=depth;i>=0;i--) {
3553          which[i]=k;
3554          k=statistics_[k]->parentNode();
3555        }
3556        printf("  node parent depth column   value                    obj      inf\n");
3557        for (i=0;i<=depth;i++) {
3558          statistics_[which[i]]->print(lookup);
3559        }
3560        delete [] which;
3561      }
3562    }
3563    // now summary
3564    int maxDepth=0;
3565    double averageSolutionDepth=0.0;
3566    int numberSolutions=0;
3567    double averageCutoffDepth=0.0;
3568    double averageSolvedDepth=0.0;
3569    int numberCutoff=0;
3570    int numberDown=0;
3571    int numberFirstDown=0;
3572    double averageInfDown=0.0;
3573    double averageObjDown=0.0;
3574    int numberCutoffDown=0;
3575    int numberUp=0;
3576    int numberFirstUp=0;
3577    double averageInfUp=0.0;
3578    double averageObjUp=0.0;
3579    int numberCutoffUp=0;
3580    double averageNumberIterations1=0.0;
3581    double averageValue=0.0;
3582    for ( i=0;i<numberNodes2_;i++) {
3583      int depth =  statistics_[i]->depth(); 
3584      int way =  statistics_[i]->way(); 
3585      double value = statistics_[i]->value(); 
3586      double startingObjective =  statistics_[i]->startingObjective(); 
3587      int startingInfeasibility = statistics_[i]->startingInfeasibility(); 
3588      double endingObjective = statistics_[i]->endingObjective(); 
3589      int endingInfeasibility = statistics_[i]->endingInfeasibility(); 
3590      maxDepth = CoinMax(depth,maxDepth);
3591      // Only for completed
3592      averageNumberIterations1 += statistics_[i]->numberIterations();
3593      averageValue += value;
3594      if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) {
3595        numberSolutions++;
3596        averageSolutionDepth += depth;
3597      }
3598      if (endingObjective==COIN_DBL_MAX) {
3599        numberCutoff++;
3600        averageCutoffDepth += depth;
3601        if (way<0) {
3602          numberDown++;
3603          numberCutoffDown++;
3604          if (way==-1)
3605            numberFirstDown++;
3606        } else {
3607          numberUp++;
3608          numberCutoffUp++;
3609          if (way==1)
3610            numberFirstUp++;
3611        }
3612      } else {
3613        averageSolvedDepth += depth;
3614        if (way<0) {
3615          numberDown++;
3616          averageInfDown += startingInfeasibility-endingInfeasibility;
3617          averageObjDown += endingObjective-startingObjective;
3618          if (way==-1)
3619            numberFirstDown++;
3620        } else {
3621          numberUp++;
3622          averageInfUp += startingInfeasibility-endingInfeasibility;
3623          averageObjUp += endingObjective-startingObjective;
3624          if (way==1)
3625            numberFirstUp++;
3626        }
3627      }
3628    }
3629    // Now print
3630    if (numberSolutions)
3631      averageSolutionDepth /= (double) numberSolutions;
3632    int numberSolved = numberNodes2_-numberCutoff;
3633    double averageNumberIterations2=numberIterations_-averageNumberIterations1
3634      -numberIterationsAtContinuous;
3635    if(numberCutoff) {
3636      averageCutoffDepth /= (double) numberCutoff;
3637      averageNumberIterations2 /= (double) numberCutoff;
3638    }
3639    if (numberNodes2_) 
3640      averageValue /= (double) numberNodes2_;
3641    if (numberSolved) {
3642      averageNumberIterations1 /= (double) numberSolved;
3643      averageSolvedDepth /= (double) numberSolved;
3644    }
3645    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
3646           numberSolutions,averageSolutionDepth);
3647    printf("average value of variable being branched on was %g\n",
3648           averageValue);
3649    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
3650           numberCutoff,averageCutoffDepth,averageNumberIterations2);
3651    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
3652           numberSolved,averageSolvedDepth,averageNumberIterations1);
3653    if (numberDown) {
3654      averageInfDown /= (double) numberDown;
3655      averageObjDown /= (double) numberDown;
3656    }
3657    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
3658           numberDown,numberFirstDown,numberDown-numberFirstDown,numberCutoffDown,
3659           averageInfDown,averageObjDown);
3660    if (numberUp) {
3661      averageInfUp /= (double) numberUp;
3662      averageObjUp /= (double) numberUp;
3663    }
3664    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
3665           numberUp,numberFirstUp,numberUp-numberFirstUp,numberCutoffUp,
3666           averageInfUp,averageObjUp);
3667    for ( i=0;i<numberNodes2_;i++) 
3668      delete statistics_[i];
3669    delete [] statistics_;
3670    statistics_=NULL;
3671    maximumStatistics_=0;
3672    delete [] lookup;
3673  }
3674/*
3675  If we think we have a solution, restore and confirm it with a call to
3676  setBestSolution().  We need to reset the cutoff value so as not to fathom
3677  the solution on bounds.  Note that calling setBestSolution( ..., true)
3678  leaves the continuousSolver_ bounds vectors fixed at the solution value.
3679
3680  Running resolve() here is a failsafe --- setBestSolution has already
3681  reoptimised using the continuousSolver_. If for some reason we fail to
3682  prove optimality, run the problem again after instructing the solver to
3683  tell us more.
3684
3685  If all looks good, replace solver_ with continuousSolver_, so that the
3686  outside world will be able to obtain information about the solution using
3687  public methods.
3688*/
3689  if (bestSolution_&&(solverCharacteristics_->solverType()<2||solverCharacteristics_->solverType()==4)) 
3690  { setCutoff(1.0e50) ; // As best solution should be worse than cutoff
3691    phase_=5;
3692    double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
3693    if ((specialOptions_&4)==0)
3694      bestObjective_ += 100.0*increment+1.0e-3; // only set if we are going to solve
3695    setBestSolution(CBC_END_SOLUTION,bestObjective_,bestSolution_,1) ;
3696    continuousSolver_->resolve() ;
3697    if (!continuousSolver_->isProvenOptimal())
3698    { continuousSolver_->messageHandler()->setLogLevel(2) ;
3699      continuousSolver_->initialSolve() ; }
3700    delete solver_ ;
3701    // above deletes solverCharacteristics_
3702    solverCharacteristics_ = NULL;
3703    solver_ = continuousSolver_ ;
3704    setPointers(solver_);
3705    continuousSolver_ = NULL ; }
3706/*
3707  Clean up dangling objects. continuousSolver_ may already be toast.
3708*/
3709  delete lastws ;
3710  if (saveObjects) {
3711    for (int i=0;i<numberObjects_;i++)
3712      delete saveObjects[i];
3713    delete [] saveObjects;
3714  }
3715  numberStrong_ = saveNumberStrong;
3716  numberBeforeTrust_ = saveNumberBeforeTrust;
3717  delete [] whichGenerator_ ;
3718  whichGenerator_=NULL;
3719  delete [] lowerBefore ;
3720  delete [] upperBefore ;
3721  delete [] walkback_ ;
3722  walkback_ = NULL ;
3723  delete [] addedCuts_ ;
3724  addedCuts_ = NULL ;
3725  //delete persistentInfo;
3726  // Get rid of characteristics
3727  solverCharacteristics_=NULL;
3728  if (continuousSolver_)
3729  { delete continuousSolver_ ;
3730    continuousSolver_ = NULL ; }
3731/*
3732  Destroy global cuts by replacing with an empty OsiCuts object.
3733*/
3734  globalCuts_= OsiCuts() ;
3735  if (!bestSolution_) {
3736    // make sure lp solver is infeasible
3737    int numberColumns = solver_->getNumCols();
3738    const double * columnLower = solver_->getColLower();
3739    int iColumn;
3740    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3741      if (solver_->isInteger(iColumn)) 
3742        solver_->setColUpper(iColumn,columnLower[iColumn]);
3743    }
3744    solver_->initialSolve();
3745  }
3746#ifdef COIN_HAS_CLP
3747  {
3748    OsiClpSolverInterface * clpSolver
3749      = dynamic_cast<OsiClpSolverInterface *> (solver_);
3750    if (clpSolver) {
3751      // Possible restore of pivot method
3752      if(savePivotMethod) {
3753        // model may have changed
3754        savePivotMethod->setModel(NULL);
3755        clpSolver->getModelPtr()->setDualRowPivotAlgorithm(*savePivotMethod);
3756        delete savePivotMethod;
3757      }
3758      clpSolver->setLargestAway(-1.0);
3759    }
3760  }
3761#endif
3762  if(fastNodeDepth_>=1000000&&!parentModel_) {
3763    // delete object off end
3764    delete object_[numberObjects_];
3765    fastNodeDepth_ -= 1000000;
3766  }
3767  delete saveSolver;
3768#if COIN_DEVELOP>1
3769  void print_fac_stats();
3770  //if (!parentModel_)
3771  //  print_fac_stats();
3772#endif
3773  if (strategy_&&strategy_->preProcessState()>0) {
3774    // undo preprocessing
3775    CglPreProcess * process = strategy_->process();
3776    assert (process);
3777    int n = originalSolver->getNumCols();
3778    if (bestSolution_) {
3779      delete [] bestSolution_;
3780      bestSolution_ = new double [n];
3781      process->postProcess(*solver_);
3782    }
3783    strategy_->deletePreProcess();
3784    // Solution now back in originalSolver
3785    delete solver_;
3786    solver_=originalSolver;
3787    if (bestSolution_)
3788      memcpy(bestSolution_,solver_->getColSolution(),n*sizeof(double));
3789    // put back original objects if there were any
3790    if (originalObject) {
3791      int iColumn;
3792      assert (ownObjects_);
3793      for (iColumn=0;iColumn<numberObjects_;iColumn++) 
3794        delete object_[iColumn];
3795      delete [] object_;
3796      numberObjects_ = numberOriginalObjects;
3797      object_=originalObject;
3798      delete [] integerVariable_;
3799      numberIntegers_=0;
3800      for (iColumn=0;iColumn<n;iColumn++) {
3801        if (solver_->isInteger(iColumn))
3802          numberIntegers_++;
3803      }
3804      integerVariable_ = new int[numberIntegers_];
3805      numberIntegers_=0;
3806      for (iColumn=0;iColumn<n;iColumn++) {
3807        if (solver_->isInteger(iColumn))
3808            integerVariable_[numberIntegers_++]=iColumn;
3809      }
3810    }
3811  }
3812#ifdef CLP_QUICK_OPTIONS
3813  {
3814    OsiClpSolverInterface * clpSolver
3815      = dynamic_cast<OsiClpSolverInterface *> (solver_);
3816    if (clpSolver) {
3817      // Try and re-use regions   
3818      ClpSimplex * simplex = clpSolver->getModelPtr();
3819      simplex->setPersistenceFlag(0);
3820      clpSolver->deleteScaleFactors();
3821      clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~131072));
3822      simplex->setSpecialOptions(simplex->specialOptions()&(~131072));
3823      simplex->setAlphaAccuracy(-1.0);
3824      //clpSolver->setSpecialOptions((clpSolver->specialOptions()&~128)|65536);
3825    }
3826  }
3827#endif
3828  return ;
3829 }
3830
3831
3832// Solve the initial LP relaxation
3833void 
3834CbcModel::initialSolve() 
3835{
3836  assert (solver_);
3837  assert (!solverCharacteristics_);
3838  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
3839  if (solverCharacteristics) {
3840    solverCharacteristics_ = solverCharacteristics;
3841  } else {
3842    // replace in solver
3843    OsiBabSolver defaultC;
3844    solver_->setAuxiliaryInfo(&defaultC);
3845    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
3846  }
3847  solverCharacteristics_->setSolver(solver_);
3848  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3849  solver_->initialSolve();
3850  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
3851  // But set up so Jon Lee will be happy
3852  status_=-1;
3853  secondaryStatus_ = -1;
3854  originalContinuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
3855  delete [] continuousSolution_;
3856  continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
3857                                             solver_->getNumCols());
3858  setPointers(solver_);
3859  solverCharacteristics_ = NULL;
3860}
3861
3862/*! \brief Get an empty basis object
3863
3864  Return an empty CoinWarmStartBasis object with the requested capacity,
3865  appropriate for the current solver. The object is cloned from the object
3866  cached as emptyWarmStart_. If there is no cached object, the routine
3867  queries the solver for a warm start object, empties it, and caches the
3868  result.
3869*/
3870
3871CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
3872
3873{ CoinWarmStartBasis *emptyBasis ;
3874/*
3875  Acquire an empty basis object, if we don't yet have one.
3876*/
3877  if (emptyWarmStart_ == 0)
3878  { if (solver_ == 0)
3879    { throw CoinError("Cannot construct basis without solver!",
3880                      "getEmptyBasis","CbcModel") ; }
3881    emptyBasis =
3882        dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
3883    if (emptyBasis == 0)
3884    { throw CoinError(
3885        "Solver does not appear to use a basis-oriented warm start.",
3886        "getEmptyBasis","CbcModel") ; }
3887    emptyBasis->setSize(0,0) ;
3888    emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ; }
3889/*
3890  Clone the empty basis object, resize it as requested, and return.
3891*/
3892  emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
3893  assert(emptyBasis) ;
3894  if (ns != 0 || na != 0) emptyBasis->setSize(ns,na) ;
3895
3896  return (emptyBasis) ; }
3897   
3898
3899/** Default Constructor
3900
3901  Creates an empty model without an associated solver.
3902*/
3903CbcModel::CbcModel() 
3904
3905:
3906  solver_(NULL),
3907  ownership_(0x80000000),
3908  continuousSolver_(NULL),
3909  referenceSolver_(NULL),
3910  defaultHandler_(true),
3911  emptyWarmStart_(NULL),
3912  bestObjective_(COIN_DBL_MAX),
3913  bestPossibleObjective_(COIN_DBL_MAX),
3914  sumChangeObjective1_(0.0),
3915  sumChangeObjective2_(0.0),
3916  bestSolution_(NULL),
3917  currentSolution_(NULL),
3918  testSolution_(NULL),
3919  minimumDrop_(1.0e-4),
3920  numberSolutions_(0),
3921  stateOfSearch_(0),
3922  hotstartSolution_(NULL),
3923  hotstartPriorities_(NULL),
3924  numberHeuristicSolutions_(0),
3925  numberNodes_(0),
3926  numberNodes2_(0),
3927  numberIterations_(0),
3928  status_(-1),
3929  secondaryStatus_(-1),
3930  numberIntegers_(0),
3931  numberRowsAtContinuous_(0),
3932  maximumNumberCuts_(0),
3933  phase_(0),
3934  currentNumberCuts_(0),
3935  maximumDepth_(0),
3936  walkback_(NULL),
3937  addedCuts_(NULL),
3938  nextRowCut_(NULL),
3939  currentNode_(NULL),
3940  integerVariable_(NULL),
3941  integerInfo_(NULL),
3942  continuousSolution_(NULL),
3943  usedInSolution_(NULL),
3944  specialOptions_(0),
3945  subTreeModel_(NULL),
3946  numberStoppedSubTrees_(0),
3947  mutex_(NULL),
3948  presolve_(0),
3949  numberStrong_(5),
3950  numberBeforeTrust_(10),
3951  numberPenalties_(20),
3952  stopNumberIterations_(-1),
3953  penaltyScaleFactor_(3.0),
3954  numberAnalyzeIterations_(0),
3955  analyzeResults_(NULL),
3956  numberInfeasibleNodes_(0),
3957  problemType_(0),
3958  printFrequency_(0),
3959  numberCutGenerators_(0),
3960  generator_(NULL),
3961  virginGenerator_(NULL),
3962  numberHeuristics_(0),
3963  heuristic_(NULL),
3964  lastHeuristic_(NULL),
3965# ifdef COIN_HAS_CLP
3966  fastNodeDepth_(-1),
3967#endif
3968  eventHandler_(NULL),
3969  numberObjects_(0),
3970  object_(NULL),
3971  ownObjects_(true),
3972  originalColumns_(NULL),
3973  howOftenGlobalScan_(1),
3974  numberGlobalViolations_(0),
3975  numberExtraIterations_(0),
3976  numberExtraNodes_(0),
3977  continuousObjective_(COIN_DBL_MAX),
3978  originalContinuousObjective_(COIN_DBL_MAX),
3979  continuousInfeasibilities_(COIN_INT_MAX),
3980  maximumCutPassesAtRoot_(20),
3981  maximumCutPasses_(10),
3982  preferredWay_(0),
3983  currentPassNumber_(0),
3984  maximumWhich_(1000),
3985  maximumRows_(0),
3986  currentDepth_(0),
3987  whichGenerator_(NULL),
3988  maximumStatistics_(0),
3989  statistics_(NULL),
3990  maximumDepthActual_(0),
3991  numberDJFixed_(0.0),
3992  probingInfo_(NULL),
3993  numberFixedAtRoot_(0),
3994  numberFixedNow_(0),
3995  stoppedOnGap_(false),
3996  eventHappened_(false),
3997  numberLongStrong_(0),
3998  numberOldActiveCuts_(0),
3999  numberNewCuts_(0),
4000  sizeMiniTree_(0),
4001  searchStrategy_(-1),
4002  numberStrongIterations_(0),
4003  resolveAfterTakeOffCuts_(true),
4004  maximumNumberIterations_(-1),
4005#if NEW_UPDATE_OBJECT>1
4006  numberUpdateItems_(0),
4007  maximumNumberUpdateItems_(0),
4008  updateItems_(NULL),
4009#endif
4010  numberThreads_(0),
4011  threadMode_(0)
4012{
4013  memset(intParam_,0,sizeof(intParam_));
4014  intParam_[CbcMaxNumNode] = 2147483647;
4015  intParam_[CbcMaxNumSol] = 9999999;
4016  intParam_[CbcFathomDiscipline] = 0;
4017
4018  dblParam_[CbcIntegerTolerance] = 1e-6;
4019  dblParam_[CbcInfeasibilityWeight] = 0.0;
4020  dblParam_[CbcCutoffIncrement] = 1e-5;
4021  dblParam_[CbcAllowableGap] = 1.0e-10;
4022  dblParam_[CbcAllowableFractionGap] = 0.0;
4023  dblParam_[CbcMaximumSeconds] = 1.0e100;
4024  dblParam_[CbcCurrentCutoff] = 1.0e100;
4025  dblParam_[CbcOptimizationDirection] = 1.0;
4026  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4027  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4028  dblParam_[CbcStartSeconds] = 0.0;
4029  strongInfo_[0]=0;
4030  strongInfo_[1]=0;
4031  strongInfo_[2]=0;
4032  solverCharacteristics_ = NULL;
4033  nodeCompare_=new CbcCompareDefault();;
4034  problemFeasibility_=new CbcFeasibilityBase();
4035  tree_= new CbcTree();
4036  branchingMethod_=NULL;
4037  cutModifier_=NULL;
4038  strategy_=NULL;
4039  parentModel_=NULL;
4040  cbcColLower_ = NULL;
4041  cbcColUpper_ = NULL;
4042  cbcRowLower_ = NULL;
4043  cbcRowUpper_ = NULL;
4044  cbcColSolution_ = NULL;
4045  cbcRowPrice_ = NULL;
4046  cbcReducedCost_ = NULL;
4047  cbcRowActivity_ = NULL;
4048  appData_=NULL;
4049  handler_ = new CoinMessageHandler();
4050  handler_->setLogLevel(2);
4051  messages_ = CbcMessage();
4052  eventHandler_ = new CbcEventHandler() ;
4053}
4054
4055/** Constructor from solver.
4056
4057  Creates a model complete with a clone of the solver passed as a parameter.
4058*/
4059
4060CbcModel::CbcModel(const OsiSolverInterface &rhs)
4061:
4062  continuousSolver_(NULL),
4063  referenceSolver_(NULL),
4064  defaultHandler_(true),
4065  emptyWarmStart_(NULL),
4066  bestObjective_(COIN_DBL_MAX),
4067  bestPossibleObjective_(COIN_DBL_MAX),
4068  sumChangeObjective1_(0.0),
4069  sumChangeObjective2_(0.0),
4070  minimumDrop_(1.0e-4),
4071  numberSolutions_(0),
4072  stateOfSearch_(0),
4073  hotstartSolution_(NULL),
4074  hotstartPriorities_(NULL),
4075  numberHeuristicSolutions_(0),
4076  numberNodes_(0),
4077  numberNodes2_(0),
4078  numberIterations_(0),
4079  status_(-1),
4080  secondaryStatus_(-1),
4081  numberRowsAtContinuous_(0),
4082  maximumNumberCuts_(0),
4083  phase_(0),
4084  currentNumberCuts_(0),
4085  maximumDepth_(0),
4086  walkback_(NULL),
4087  addedCuts_(NULL),
4088  nextRowCut_(NULL),
4089  currentNode_(NULL),
4090  integerInfo_(NULL),
4091  specialOptions_(0),
4092  subTreeModel_(NULL),
4093  numberStoppedSubTrees_(0),
4094  mutex_(NULL),
4095  presolve_(0),
4096  numberStrong_(5),
4097  numberBeforeTrust_(10),
4098  numberPenalties_(20),
4099  stopNumberIterations_(-1),
4100  penaltyScaleFactor_(3.0),
4101  numberAnalyzeIterations_(0),
4102  analyzeResults_(NULL),
4103  numberInfeasibleNodes_(0),
4104  problemType_(0),
4105  printFrequency_(0),
4106  numberCutGenerators_(0),
4107  generator_(NULL),
4108  virginGenerator_(NULL),
4109  numberHeuristics_(0),
4110  heuristic_(NULL),
4111  lastHeuristic_(NULL),
4112# ifdef COIN_HAS_CLP
4113  fastNodeDepth_(-1),
4114#endif
4115  eventHandler_(NULL),
4116  numberObjects_(0),
4117  object_(NULL),
4118  ownObjects_(true),
4119  originalColumns_(NULL),
4120  howOftenGlobalScan_(1),
4121  numberGlobalViolations_(0),
4122  numberExtraIterations_(0),
4123  numberExtraNodes_(0),
4124  continuousObjective_(COIN_DBL_MAX),
4125  originalContinuousObjective_(COIN_DBL_MAX),
4126  continuousInfeasibilities_(COIN_INT_MAX),
4127  maximumCutPassesAtRoot_(20),
4128  maximumCutPasses_(10),
4129  preferredWay_(0),
4130  currentPassNumber_(0),
4131  maximumWhich_(1000),
4132  maximumRows_(0),
4133  currentDepth_(0),
4134  whichGenerator_(NULL),
4135  maximumStatistics_(0),
4136  statistics_(NULL),
4137  maximumDepthActual_(0),
4138  numberDJFixed_(0.0),
4139  probingInfo_(NULL),
4140  numberFixedAtRoot_(0),
4141  numberFixedNow_(0),
4142  stoppedOnGap_(false),
4143  eventHappened_(false),
4144  numberLongStrong_(0),
4145  numberOldActiveCuts_(0),
4146  numberNewCuts_(0),
4147  sizeMiniTree_(0),
4148  searchStrategy_(-1),
4149  numberStrongIterations_(0),
4150  resolveAfterTakeOffCuts_(true),
4151  maximumNumberIterations_(-1),
4152#if NEW_UPDATE_OBJECT>1
4153  numberUpdateItems_(0),
4154  maximumNumberUpdateItems_(0),
4155  updateItems_(NULL),
4156#endif
4157  numberThreads_(0),
4158  threadMode_(0)
4159{
4160  memset(intParam_,0,sizeof(intParam_));
4161  intParam_[CbcMaxNumNode] = 2147483647;
4162  intParam_[CbcMaxNumSol] = 9999999;
4163  intParam_[CbcFathomDiscipline] = 0;
4164
4165  dblParam_[CbcIntegerTolerance] = 1e-6;
4166  dblParam_[CbcInfeasibilityWeight] = 0.0;
4167  dblParam_[CbcCutoffIncrement] = 1e-5;
4168  dblParam_[CbcAllowableGap] = 1.0e-10;
4169  dblParam_[CbcAllowableFractionGap] = 0.0;
4170  dblParam_[CbcMaximumSeconds] = 1.0e100;
4171  dblParam_[CbcCurrentCutoff] = 1.0e100;
4172  dblParam_[CbcOptimizationDirection] = 1.0;
4173  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4174  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4175  dblParam_[CbcStartSeconds] = 0.0;
4176  strongInfo_[0]=0;
4177  strongInfo_[1]=0;
4178  strongInfo_[2]=0;
4179  solverCharacteristics_ = NULL;
4180
4181  nodeCompare_=new CbcCompareDefault();;
4182  problemFeasibility_=new CbcFeasibilityBase();
4183  tree_= new CbcTree();
4184  branchingMethod_=NULL;
4185  cutModifier_=NULL;
4186  strategy_=NULL;
4187  parentModel_=NULL;
4188  appData_=NULL;
4189  handler_ = new CoinMessageHandler();
4190  handler_->setLogLevel(2);
4191  messages_ = CbcMessage();
4192  eventHandler_ = new CbcEventHandler() ;
4193  solver_ = rhs.clone();
4194  referenceSolver_ = solver_->clone();
4195  ownership_ = 0x80000000;
4196  cbcColLower_ = NULL;
4197  cbcColUpper_ = NULL;
4198  cbcRowLower_ = NULL;
4199  cbcRowUpper_ = NULL;
4200  cbcColSolution_ = NULL;
4201  cbcRowPrice_ = NULL;
4202  cbcReducedCost_ = NULL;
4203  cbcRowActivity_ = NULL;
4204
4205  // Initialize solution and integer variable vectors
4206  bestSolution_ = NULL; // to say no solution found
4207  numberIntegers_=0;
4208  int numberColumns = solver_->getNumCols();
4209  int iColumn;
4210  if (numberColumns) {
4211    // Space for current solution
4212    currentSolution_ = new double[numberColumns];
4213    continuousSolution_ = new double[numberColumns];
4214    usedInSolution_ = new int[numberColumns];
4215    CoinZeroN(usedInSolution_,numberColumns);
4216    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4217      if( solver_->isInteger(iColumn)) 
4218        numberIntegers_++;
4219    }
4220  } else {
4221    // empty model
4222    currentSolution_=NULL;
4223    continuousSolution_=NULL;
4224    usedInSolution_=NULL;
4225  }
4226  testSolution_=currentSolution_;
4227  if (numberIntegers_) {
4228    integerVariable_ = new int [numberIntegers_];
4229    numberIntegers_=0;
4230    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4231      if( solver_->isInteger(iColumn)) 
4232        integerVariable_[numberIntegers_++]=iColumn;
4233    }
4234  } else {
4235    integerVariable_ = NULL;
4236  }
4237}
4238
4239/*
4240  Assign a solver to the model (model assumes ownership)
4241
4242  The integer variable vector is initialized if it's not already present.
4243  If deleteSolver then current solver deleted (if model owned)
4244
4245  Assuming ownership matches usage in OsiSolverInterface
4246  (cf. assignProblem, loadProblem).
4247
4248  TODO: What to do about solver parameters? A simple copy likely won't do it,
4249        because the SI must push the settings into the underlying solver. In
4250        the context of switching solvers in cbc, this means that command line
4251        settings will get lost. Stash the command line somewhere and reread it
4252        here, maybe?
4253 
4254  TODO: More generally, how much state should be transferred from the old
4255        solver to the new solver? Best perhaps to see how usage develops.
4256        What's done here mimics the CbcModel(OsiSolverInterface) constructor.
4257*/
4258void
4259CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver)
4260
4261{
4262  // resize best solution if exists
4263  if (bestSolution_&&solver&&solver_) {
4264    int nOld = solver_->getNumCols();
4265    int nNew = solver->getNumCols();
4266    if (nNew>nOld) {
4267      double * temp = new double[nNew];
4268      memcpy(temp,bestSolution_,nOld*sizeof(double));
4269      memset(temp+nOld,0,(nNew-nOld)*sizeof(double));
4270      delete [] bestSolution_;
4271      bestSolution_=temp;
4272    }
4273  }
4274  // Keep the current message level for solver (if solver exists)
4275  if (solver_)
4276    solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
4277
4278  if (modelOwnsSolver()&&deleteSolver) delete solver_ ;
4279  solver_ = solver;
4280  solver = NULL ;
4281  setModelOwnsSolver(true) ;
4282/*
4283  Basis information is solver-specific.
4284*/
4285  if (emptyWarmStart_)
4286  { delete emptyWarmStart_  ;
4287    emptyWarmStart_ = 0 ; }
4288  bestSolutionBasis_ = CoinWarmStartBasis();
4289/*
4290  Initialize integer variable vector.
4291*/
4292  numberIntegers_=0;
4293  int numberColumns = solver_->getNumCols();
4294  int iColumn;
4295  for (iColumn=0;iColumn<numberColumns;iColumn++) {
4296    if( solver_->isInteger(iColumn)) 
4297      numberIntegers_++;
4298  }
4299  delete [] integerVariable_;
4300  if (numberIntegers_) {
4301    integerVariable_ = new int [numberIntegers_];
4302    numberIntegers_=0;
4303    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4304      if( solver_->isInteger(iColumn)) 
4305        integerVariable_[numberIntegers_++]=iColumn;
4306    }
4307  } else {
4308    integerVariable_ = NULL;
4309  }
4310
4311  return ;
4312}
4313
4314// Copy constructor.
4315
4316CbcModel::CbcModel(const CbcModel & rhs, bool cloneHandler)
4317:
4318  continuousSolver_(NULL),
4319  referenceSolver_(NULL),
4320  defaultHandler_(rhs.defaultHandler_),
4321  emptyWarmStart_(NULL),
4322  bestObjective_(rhs.bestObjective_),
4323  bestPossibleObjective_(rhs.bestPossibleObjective_),
4324  sumChangeObjective1_(rhs.sumChangeObjective1_),
4325  sumChangeObjective2_(rhs.sumChangeObjective2_),
4326  minimumDrop_(rhs.minimumDrop_),
4327  numberSolutions_(rhs.numberSolutions_),
4328  stateOfSearch_(rhs.stateOfSearch_),
4329  numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
4330  numberNodes_(rhs.numberNodes_),
4331  numberNodes2_(rhs.numberNodes2_),
4332  numberIterations_(rhs.numberIterations_),
4333  status_(rhs.status_),
4334  secondaryStatus_(rhs.secondaryStatus_),
4335  specialOptions_(rhs.specialOptions_),
4336  subTreeModel_(rhs.subTreeModel_),
4337  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
4338  mutex_(NULL),
4339  presolve_(rhs.presolve_),
4340  numberStrong_(rhs.numberStrong_),
4341  numberBeforeTrust_(rhs.numberBeforeTrust_),
4342  numberPenalties_(rhs.numberPenalties_),
4343  stopNumberIterations_(rhs.stopNumberIterations_),
4344  penaltyScaleFactor_(rhs.penaltyScaleFactor_),
4345  numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
4346  analyzeResults_(NULL),
4347  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
4348  problemType_(rhs.problemType_),
4349  printFrequency_(rhs.printFrequency_),
4350# ifdef COIN_HAS_CLP
4351  fastNodeDepth_(rhs.fastNodeDepth_),
4352#endif
4353  howOftenGlobalScan_(rhs.howOftenGlobalScan_),
4354  numberGlobalViolations_(rhs.numberGlobalViolations_),
4355  numberExtraIterations_(rhs.numberExtraIterations_),
4356  numberExtraNodes_(rhs.numberExtraNodes_),
4357  continuousObjective_(rhs.continuousObjective_),
4358  originalContinuousObjective_(rhs.originalContinuousObjective_),
4359  continuousInfeasibilities_(rhs.continuousInfeasibilities_),
4360  maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
4361  maximumCutPasses_( rhs.maximumCutPasses_),
4362  preferredWay_(rhs.preferredWay_),
4363  currentPassNumber_(rhs.currentPassNumber_),
4364  maximumWhich_(rhs.maximumWhich_),
4365  maximumRows_(0),
4366  currentDepth_(0),
4367  whichGenerator_(NULL),
4368  maximumStatistics_(0),
4369  statistics_(NULL),
4370  maximumDepthActual_(0),
4371  numberDJFixed_(0.0),
4372  probingInfo_(NULL),
4373  numberFixedAtRoot_(rhs.numberFixedAtRoot_),
4374  numberFixedNow_(rhs.numberFixedNow_),
4375  stoppedOnGap_(rhs.stoppedOnGap_),
4376  eventHappened_(rhs.eventHappened_),
4377  numberLongStrong_(rhs.numberLongStrong_),
4378  numberOldActiveCuts_(rhs.numberOldActiveCuts_),
4379  numberNewCuts_(rhs.numberNewCuts_),
4380  sizeMiniTree_(rhs.sizeMiniTree_),
4381  searchStrategy_(rhs.searchStrategy_),
4382  numberStrongIterations_(rhs.numberStrongIterations_),
4383  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
4384  maximumNumberIterations_(rhs.maximumNumberIterations_),
4385#if NEW_UPDATE_OBJECT>1
4386  numberUpdateItems_(rhs.numberUpdateItems_),
4387  maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_),
4388  updateItems_(NULL),
4389#endif
4390  numberThreads_(rhs.numberThreads_),
4391  threadMode_(rhs.threadMode_)
4392{
4393  memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
4394  memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
4395  strongInfo_[0]=rhs.strongInfo_[0];
4396  strongInfo_[1]=rhs.strongInfo_[1];
4397  strongInfo_[2]=rhs.strongInfo_[2];
4398  solverCharacteristics_ = NULL;
4399  if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
4400  if (defaultHandler_||cloneHandler) {
4401    handler_ = new CoinMessageHandler();
4402    handler_->setLogLevel(2);
4403  } else {
4404    handler_ = rhs.handler_;
4405  }
4406  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4407  numberCutGenerators_ = rhs.numberCutGenerators_;
4408  if (numberCutGenerators_) {
4409    generator_ = new CbcCutGenerator * [numberCutGenerators_];
4410    virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
4411    int i;
4412    for (i=0;i<numberCutGenerators_;i++) {
4413      generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
4414      virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
4415    }
4416  } else {
4417    generator_=NULL;
4418    virginGenerator_=NULL;
4419  }
4420  globalCuts_ = rhs.globalCuts_;
4421  numberHeuristics_ = rhs.numberHeuristics_;
4422  if (numberHeuristics_) {
4423    heuristic_ = new CbcHeuristic * [numberHeuristics_];
4424    int i;
4425    for (i=0;i<numberHeuristics_;i++) {
4426      heuristic_[i]=rhs.heuristic_[i]->clone();
4427    }
4428  } else {
4429    heuristic_=NULL;
4430  }
4431  lastHeuristic_ = NULL;
4432  if (rhs.eventHandler_)
4433    { eventHandler_ = rhs.eventHandler_->clone() ; }
4434  else
4435  { eventHandler_ = NULL ; }
4436  ownObjects_ = rhs.ownObjects_;
4437  if (ownObjects_) {
4438    numberObjects_=rhs.numberObjects_;
4439    if (numberObjects_) {
4440      object_ = new OsiObject * [numberObjects_];
4441      int i;
4442      for (i=0;i<numberObjects_;i++) {
4443        object_[i]=(rhs.object_[i])->clone();
4444        CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
4445        // Could be OsiObjects
4446        if (obj)
4447          obj->setModel(this);
4448      }
4449    } else {
4450      object_=NULL;
4451    }
4452  } else {
4453    // assume will be redone
4454    numberObjects_=0;
4455    object_=NULL;
4456  }
4457  if (rhs.referenceSolver_)
4458    referenceSolver_ = rhs.referenceSolver_->clone();
4459  else
4460    referenceSolver_=NULL;
4461  solver_ = rhs.solver_->clone();
4462  if (rhs.originalColumns_) {
4463    int numberColumns = solver_->getNumCols();
4464    originalColumns_= new int [numberColumns];
4465    memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
4466  } else {
4467    originalColumns_=NULL;
4468  }
4469#if NEW_UPDATE_OBJECT>1
4470  if (maximumNumberUpdateItems_) {
4471    updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
4472    for (int i=0;i<maximumNumberUpdateItems_;i++)
4473      updateItems_[i] = rhs.updateItems_[i];
4474  }
4475#endif
4476  if (maximumWhich_&&rhs.whichGenerator_)
4477    whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
4478  nodeCompare_=rhs.nodeCompare_->clone();
4479  problemFeasibility_=rhs.problemFeasibility_->clone();
4480  tree_= rhs.tree_->clone();
4481  if (rhs.branchingMethod_)
4482    branchingMethod_=rhs.branchingMethod_->clone();
4483  else
4484    branchingMethod_=NULL;
4485  if (rhs.cutModifier_)
4486    cutModifier_=rhs.cutModifier_->clone();
4487  else
4488    cutModifier_=NULL;
4489  cbcColLower_ = NULL;
4490  cbcColUpper_ = NULL;
4491  cbcRowLower_ = NULL;
4492  cbcRowUpper_ = NULL;
4493  cbcColSolution_ = NULL;
4494  cbcRowPrice_ = NULL;
4495  cbcReducedCost_ = NULL;
4496  cbcRowActivity_ = NULL;
4497  if (rhs.strategy_)
4498    strategy_=rhs.strategy_->clone();
4499  else
4500    strategy_=NULL;
4501  parentModel_=rhs.parentModel_;
4502  appData_=rhs.appData_;
4503  messages_ = rhs.messages_;
4504  ownership_ = 0x80000000;
4505  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4506  numberIntegers_=rhs.numberIntegers_;
4507  randomNumberGenerator_ = rhs.randomNumberGenerator_;
4508  if (numberIntegers_) {
4509    integerVariable_ = new int [numberIntegers_];
4510    memcpy(integerVariable_,rhs.integerVariable_,numberIntegers_*sizeof(int));
4511    integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
4512  } else {
4513    integerVariable_ = NULL;
4514    integerInfo_=NULL;
4515  }
4516  if (rhs.hotstartSolution_) {
4517    int numberColumns = solver_->getNumCols();
4518    hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
4519    hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
4520  } else {
4521    hotstartSolution_ = NULL;
4522    hotstartPriorities_ =NULL;
4523  }
4524  if (rhs.bestSolution_) {
4525    int numberColumns = solver_->getNumCols();
4526    bestSolution_ = new double[numberColumns];
4527    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
4528  } else {
4529    bestSolution_=NULL;
4530  }
4531  int numberColumns = solver_->getNumCols();
4532  // Space for current solution
4533  currentSolution_ = new double[numberColumns];
4534  continuousSolution_ = new double[numberColumns];
4535  usedInSolution_ = new int[numberColumns];
4536  CoinZeroN(usedInSolution_,numberColumns);
4537  testSolution_=currentSolution_;
4538  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
4539  maximumNumberCuts_=rhs.maximumNumberCuts_;
4540  phase_ = rhs.phase_;
4541  currentNumberCuts_=rhs.currentNumberCuts_;
4542  maximumDepth_= rhs.maximumDepth_;
4543  // These are only used as temporary arrays so need not be filled
4544  if (maximumNumberCuts_) {
4545    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4546  } else {
4547    addedCuts_ = NULL;
4548  }
4549  bestSolutionBasis_ = rhs.bestSolutionBasis_;
4550  nextRowCut_ = NULL;
4551  currentNode_ = NULL;
4552  if (maximumDepth_)
4553    walkback_ = new CbcNodeInfo * [maximumDepth_];
4554  else
4555    walkback_ = NULL;
4556  synchronizeModel();
4557  if (cloneHandler) {
4558    delete handler_;
4559    CoinMessageHandler * handler = rhs.handler_->clone();
4560    passInMessageHandler(handler);
4561  }
4562}
4563 
4564// Assignment operator
4565CbcModel & 
4566CbcModel::operator=(const CbcModel& rhs)
4567{
4568  if (this!=&rhs) {
4569    if (modelOwnsSolver()) {
4570      delete solver_;
4571      solver_=NULL;
4572    }
4573    gutsOfDestructor();
4574    if (defaultHandler_)
4575    { delete handler_;
4576      handler_ = NULL; }
4577    defaultHandler_ = rhs.defaultHandler_;
4578    if (defaultHandler_)
4579    { handler_ = new CoinMessageHandler();
4580      handler_->setLogLevel(2); }
4581    else
4582    { handler_ = rhs.handler_; }
4583    messages_ = rhs.messages_;
4584    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4585    if (rhs.solver_)
4586    { solver_ = rhs.solver_->clone() ; }
4587    else
4588    { solver_ = 0 ; }
4589    ownership_ = 0x80000000;
4590    delete continuousSolver_ ;
4591    if (rhs.continuousSolver_)
4592    { continuousSolver_ = rhs.continuousSolver_->clone() ; }
4593    else
4594    { continuousSolver_ = 0 ; }
4595    delete referenceSolver_;
4596    if (rhs.referenceSolver_)
4597    { referenceSolver_ = rhs.referenceSolver_->clone() ; }
4598    else
4599    { referenceSolver_ = NULL ; }
4600
4601    delete emptyWarmStart_ ;
4602    if (rhs.emptyWarmStart_)
4603    { emptyWarmStart_ = rhs.emptyWarmStart_->clone() ; }
4604    else
4605    { emptyWarmStart_ = 0 ; }
4606
4607    bestObjective_ = rhs.bestObjective_;
4608    bestPossibleObjective_=rhs.bestPossibleObjective_;
4609    sumChangeObjective1_=rhs.sumChangeObjective1_;
4610    sumChangeObjective2_=rhs.sumChangeObjective2_;
4611    delete [] bestSolution_;
4612    if (rhs.bestSolution_) {
4613      int numberColumns = rhs.getNumCols();
4614      bestSolution_ = new double[numberColumns];
4615      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
4616    } else {
4617      bestSolution_=NULL;
4618    }
4619    int numberColumns = rhs.getNumCols();
4620    if (numberColumns) {
4621      // Space for current solution
4622      currentSolution_ = new double[numberColumns];
4623      continuousSolution_ = new double[numberColumns];
4624      usedInSolution_ = new int[numberColumns];
4625      CoinZeroN(usedInSolution_,numberColumns);
4626    } else {
4627      currentSolution_=NULL;
4628      continuousSolution_=NULL;
4629      usedInSolution_=NULL;
4630    }
4631    testSolution_=currentSolution_;
4632    minimumDrop_ = rhs.minimumDrop_;
4633    numberSolutions_=rhs.numberSolutions_;
4634    stateOfSearch_= rhs.stateOfSearch_;
4635    numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
4636    numberNodes_ = rhs.numberNodes_;
4637    numberNodes2_ = rhs.numberNodes2_;
4638    numberIterations_ = rhs.numberIterations_;
4639    status_ = rhs.status_;
4640    secondaryStatus_ = rhs.secondaryStatus_;
4641    specialOptions_ = rhs.specialOptions_;
4642    subTreeModel_ = rhs.subTreeModel_;
4643    numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
4644    mutex_ = NULL;
4645    presolve_ = rhs.presolve_;
4646    numberStrong_ = rhs.numberStrong_;
4647    numberBeforeTrust_ = rhs.numberBeforeTrust_;
4648    numberPenalties_ = rhs.numberPenalties_;
4649    stopNumberIterations_ = rhs.stopNumberIterations_;
4650    penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
4651    numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
4652    delete [] analyzeResults_;
4653    analyzeResults_ = NULL;
4654    numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
4655    problemType_ = rhs.problemType_;
4656    printFrequency_ = rhs.printFrequency_;
4657    howOftenGlobalScan_=rhs.howOftenGlobalScan_;
4658    numberGlobalViolations_=rhs.numberGlobalViolations_;
4659    numberExtraIterations_ = rhs.numberExtraIterations_;
4660    numberExtraNodes_ = rhs.numberExtraNodes_;
4661    continuousObjective_=rhs.continuousObjective_;
4662    originalContinuousObjective_ = rhs.originalContinuousObjective_;
4663    continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
4664    maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
4665    maximumCutPasses_ = rhs.maximumCutPasses_;
4666    preferredWay_ = rhs.preferredWay_;
4667    currentPassNumber_ = rhs.currentPassNumber_;
4668    memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
4669    memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
4670    globalCuts_ = rhs.globalCuts_;
4671    int i;
4672    for (i=0;i<numberCutGenerators_;i++) {
4673      delete generator_[i];
4674      delete virginGenerator_[i];
4675    }
4676    delete [] generator_;
4677    delete [] virginGenerator_;
4678    delete [] heuristic_;
4679    maximumWhich_ = rhs.maximumWhich_;
4680    delete [] whichGenerator_;
4681    whichGenerator_ = NULL;
4682    if (maximumWhich_&&rhs.whichGenerator_)
4683      whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
4684    maximumRows_=0;
4685    currentDepth_ = 0;
4686    randomNumberGenerator_ = rhs.randomNumberGenerator_;
4687    workingBasis_ = CoinWarmStartBasis();
4688    for (i=0;i<maximumStatistics_;i++)
4689      delete statistics_[i];
4690    delete [] statistics_;
4691    maximumStatistics_ = 0;
4692    statistics_ = NULL;
4693    delete probingInfo_;
4694    probingInfo_=NULL;
4695    numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
4696    numberFixedNow_ = rhs.numberFixedNow_;
4697    stoppedOnGap_ = rhs.stoppedOnGap_;
4698    eventHappened_ = rhs.eventHappened_;
4699    numberLongStrong_ = rhs.numberLongStrong_;
4700    numberOldActiveCuts_ = rhs.numberOldActiveCuts_;
4701    numberNewCuts_ = rhs.numberNewCuts_;
4702    resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_;
4703    maximumNumberIterations_ = rhs.maximumNumberIterations_;
4704#if NEW_UPDATE_OBJECT>1
4705    numberUpdateItems_ = rhs.numberUpdateItems_;
4706    maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_;
4707    delete [] updateItems_;
4708    if (maximumNumberUpdateItems_) {
4709      updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
4710      for (i=0;i<maximumNumberUpdateItems_;i++)
4711        updateItems_[i] = rhs.updateItems_[i];
4712    } else {
4713      updateItems_ = NULL;
4714    }
4715#endif
4716    numberThreads_ = rhs.numberThreads_;
4717    threadMode_ = rhs.threadMode_;
4718    sizeMiniTree_ = rhs.sizeMiniTree_;
4719    searchStrategy_ = rhs.searchStrategy_;
4720    numberStrongIterations_ = rhs.numberStrongIterations_;
4721    strongInfo_[0]=rhs.strongInfo_[0];
4722    strongInfo_[1]=rhs.strongInfo_[1];
4723    strongInfo_[2]=rhs.strongInfo_[2];
4724    solverCharacteristics_ = NULL;
4725    lastHeuristic_ = NULL;
4726    numberCutGenerators_ = rhs.numberCutGenerators_;
4727    if (numberCutGenerators_) {
4728      generator_ = new CbcCutGenerator * [numberCutGenerators_];
4729      virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
4730      int i;
4731      for (i=0;i<numberCutGenerators_;i++) {
4732        generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
4733        virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
4734      }
4735    } else {
4736      generator_=NULL;
4737      virginGenerator_=NULL;
4738    }
4739    numberHeuristics_ = rhs.numberHeuristics_;
4740    if (numberHeuristics_) {
4741      heuristic_ = new CbcHeuristic * [numberHeuristics_];
4742      memcpy(heuristic_,rhs.heuristic_,
4743             numberHeuristics_*sizeof(CbcHeuristic *));
4744    } else {
4745      heuristic_=NULL;
4746    }
4747    lastHeuristic_ = NULL;
4748    if (eventHandler_)
4749      delete eventHandler_ ;
4750    if (rhs.eventHandler_)
4751      { eventHandler_ = rhs.eventHandler_->clone() ; }
4752    else
4753    { eventHandler_ = NULL ; }
4754# ifdef COIN_HAS_CLP
4755    fastNodeDepth_ = rhs.fastNodeDepth_;
4756#endif
4757    if (ownObjects_) {
4758      for (i=0;i<numberObjects_;i++)
4759        delete object_[i];
4760      delete [] object_;
4761      numberObjects_=rhs.numberObjects_;
4762      if (numberObjects_) {
4763        object_ = new OsiObject * [numberObjects_];
4764        int i;
4765        for (i=0;i<numberObjects_;i++) 
4766          object_[i]=(rhs.object_[i])->clone();
4767      } else {
4768        object_=NULL;
4769    }
4770    } else {
4771      // assume will be redone
4772      numberObjects_=0;
4773      object_=NULL;
4774    }
4775    delete [] originalColumns_;
4776    if (rhs.originalColumns_) {
4777      int numberColumns = rhs.getNumCols();
4778      originalColumns_= new int [numberColumns];
4779      memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
4780    } else {
4781      originalColumns_=NULL;
4782    }
4783    nodeCompare_=rhs.nodeCompare_->clone();
4784    problemFeasibility_=rhs.problemFeasibility_->clone();
4785    delete tree_;
4786    tree_= rhs.tree_->clone();
4787    if (rhs.branchingMethod_)
4788      branchingMethod_=rhs.branchingMethod_->clone();
4789    else
4790      branchingMethod_=NULL;
4791    if (rhs.cutModifier_)
4792      cutModifier_=rhs.cutModifier_->clone();
4793    else
4794      cutModifier_=NULL;
4795    delete strategy_;
4796    if (rhs.strategy_)
4797      strategy_=rhs.strategy_->clone();
4798    else
4799      strategy_=NULL;
4800    parentModel_=rhs.parentModel_;
4801    appData_=rhs.appData_;
4802
4803    delete [] integerVariable_;
4804    numberIntegers_=rhs.numberIntegers_;
4805    if (numberIntegers_) {
4806      integerVariable_ = new int [numberIntegers_];
4807      memcpy(integerVariable_,rhs.integerVariable_,
4808             numberIntegers_*sizeof(int));
4809      integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,rhs.getNumCols());
4810    } else {
4811      integerVariable_ = NULL;
4812      integerInfo_=NULL;
4813    }
4814    if (rhs.hotstartSolution_) {
4815      int numberColumns = solver_->getNumCols();
4816      hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
4817      hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
4818    } else {
4819      hotstartSolution_ = NULL;
4820      hotstartPriorities_ =NULL;
4821    }
4822    numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
4823    maximumNumberCuts_=rhs.maximumNumberCuts_;
4824    phase_ = rhs.phase_;
4825    currentNumberCuts_=rhs.currentNumberCuts_;
4826    maximumDepth_= rhs.maximumDepth_;
4827    delete [] addedCuts_;
4828    delete [] walkback_;
4829    // These are only used as temporary arrays so need not be filled
4830    if (maximumNumberCuts_) {
4831      addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4832    } else {
4833      addedCuts_ = NULL;
4834    }
4835    bestSolutionBasis_ = rhs.bestSolutionBasis_;
4836    nextRowCut_ = NULL;
4837    currentNode_ = NULL;
4838    if (maximumDepth_)
4839      walkback_ = new CbcNodeInfo * [maximumDepth_];
4840    else
4841      walkback_ = NULL;
4842    synchronizeModel();
4843    cbcColLower_ = NULL;
4844    cbcColUpper_ = NULL;
4845    cbcRowLower_ = NULL;
4846    cbcRowUpper_ = NULL;
4847    cbcColSolution_ = NULL;
4848    cbcRowPrice_ = NULL;
4849    cbcReducedCost_ = NULL;
4850    cbcRowActivity_ = NULL;
4851  }
4852  return *this;
4853}
4854// Destructor
4855CbcModel::~CbcModel ()
4856{
4857  if (defaultHandler_) {
4858    delete handler_;
4859    handler_ = NULL;
4860  }
4861  delete tree_;
4862  tree_=NULL;
4863  if (modelOwnsSolver()) {
4864    delete solver_;
4865    solver_ = NULL;
4866  }
4867  gutsOfDestructor();
4868  delete eventHandler_ ;
4869  eventHandler_ = NULL ;
4870}
4871// Clears out as much as possible (except solver)
4872void 
4873CbcModel::gutsOfDestructor()
4874{
4875  delete referenceSolver_;
4876  referenceSolver_=NULL;
4877  int i;
4878  for (i=0;i<numberCutGenerators_;i++) {
4879    delete generator_[i];
4880    delete virginGenerator_[i];
4881  }
4882  delete [] generator_;
4883  delete [] virginGenerator_;
4884  generator_=NULL;
4885  virginGenerator_=NULL;
4886  for (i=0;i<numberHeuristics_;i++)
4887    delete heuristic_[i];
4888  delete [] heuristic_;
4889  heuristic_=NULL;
4890  delete nodeCompare_;
4891  nodeCompare_=NULL;
4892  delete problemFeasibility_;
4893  problemFeasibility_=NULL;
4894  delete [] originalColumns_;
4895  originalColumns_=NULL;
4896  delete strategy_;
4897#if NEW_UPDATE_OBJECT>1
4898  delete [] updateItems_;
4899  updateItems_=NULL;
4900  numberUpdateItems_=0;
4901  maximumNumberUpdateItems_=0;
4902#endif
4903  gutsOfDestructor2();
4904}
4905// Clears out enough to reset CbcModel
4906void 
4907CbcModel::gutsOfDestructor2()
4908{
4909  delete [] integerInfo_;
4910  integerInfo_=NULL;
4911  delete [] integerVariable_;
4912  integerVariable_=NULL;
4913  int i;
4914  if (ownObjects_) {
4915    for (i=0;i<numberObjects_;i++)
4916      delete object_[i];
4917    delete [] object_;
4918  }
4919  ownObjects_=true;
4920  object_=NULL;
4921  numberIntegers_=0;
4922  numberObjects_=0;
4923  // Below here is whatever consensus is
4924  ownership_ = 0x80000000;
4925  delete branchingMethod_;
4926  branchingMethod_=NULL;
4927  delete cutModifier_;
4928  cutModifier_=NULL;
4929  resetModel();
4930}
4931// Clears out enough to reset CbcModel
4932void 
4933CbcModel::resetModel()
4934{
4935  delete emptyWarmStart_ ;
4936  emptyWarmStart_ =NULL;
4937  delete continuousSolver_;
4938  continuousSolver_=NULL;
4939  delete [] bestSolution_;
4940  bestSolution_=NULL;
4941  delete [] currentSolution_;
4942  currentSolution_=NULL;
4943  delete [] continuousSolution_;
4944  continuousSolution_=NULL;
4945  solverCharacteristics_=NULL;
4946  delete [] usedInSolution_;
4947  usedInSolution_ = NULL;
4948  testSolution_=NULL;
4949  lastHeuristic_ = NULL;
4950  delete [] addedCuts_;
4951  addedCuts_=NULL;
4952  nextRowCut_ = NULL;
4953  currentNode_ = NULL;
4954  delete [] walkback_;
4955  walkback_=NULL;
4956  delete [] whichGenerator_;
4957  whichGenerator_ = NULL;
4958  for (int i=0;i<maximumStatistics_;i++)
4959    delete statistics_[i];
4960  delete [] statistics_;
4961  statistics_=NULL;
4962  maximumDepthActual_ = 0;
4963  numberDJFixed_ =0.0;
4964  delete probingInfo_;
4965  probingInfo_ = NULL;
4966  maximumStatistics_=0;
4967  delete [] analyzeResults_;
4968  analyzeResults_=NULL;
4969  bestObjective_=COIN_DBL_MAX;
4970  bestPossibleObjective_=COIN_DBL_MAX;
4971  sumChangeObjective1_=0.0;
4972  sumChangeObjective2_=0.0;
4973  numberSolutions_=0;
4974  stateOfSearch_=0;
4975  delete [] hotstartSolution_;
4976  hotstartSolution_=NULL;
4977  delete [] hotstartPriorities_;
4978  hotstartPriorities_=NULL;
4979  numberHeuristicSolutions_=0;
4980  numberNodes_=0;
4981  numberNodes2_=0;
4982  numberIterations_=0;
4983  status_=-1;
4984  secondaryStatus_=-1;
4985  maximumNumberCuts_=0;
4986  phase_=0;
4987  currentNumberCuts_=0;
4988  maximumDepth_=0;
4989  nextRowCut_=NULL;
4990  currentNode_=NULL;
4991  // clear out tree
4992  if (tree_&&tree_->size())
4993    tree_->cleanTree(this, -1.0e100,bestPossibleObjective_) ;
4994  subTreeModel_=NULL;
4995  numberStoppedSubTrees_=0;
4996  numberInfeasibleNodes_=0;
4997  numberGlobalViolations_=0;
4998  numberExtraIterations_ = 0;
4999  numberExtraNodes_ = 0;
5000  continuousObjective_=0.0;
5001  originalContinuousObjective_=0.0;
5002  continuousInfeasibilities_=0;
5003  numberFixedAtRoot_=0;
5004  numberFixedNow_=0;
5005  stoppedOnGap_=false;
5006  eventHappened_=false;
5007  numberLongStrong_=0;
5008  numberOldActiveCuts_=0;
5009  numberNewCuts_=0;
5010  searchStrategy_=-1;
5011  numberStrongIterations_=0;
5012  // Parameters which need to be reset
5013  setCutoff(COIN_DBL_MAX);
5014  dblParam_[CbcCutoffIncrement] = 1e-5;
5015  dblParam_[CbcCurrentCutoff] = 1.0e100;
5016  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
5017  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
5018}
5019/* Most of copy constructor
5020      mode - 0 copy but don't delete before
5021             1 copy and delete before
5022             2 copy and delete before (but use virgin generators)
5023*/
5024void 
5025CbcModel::gutsOfCopy(const CbcModel & rhs,int mode)
5026{
5027  minimumDrop_ = rhs.minimumDrop_;
5028  specialOptions_ = rhs.specialOptions_;
5029  numberStrong_ = rhs.numberStrong_;
5030  numberBeforeTrust_ = rhs.numberBeforeTrust_;
5031  numberPenalties_ = rhs.numberPenalties_;
5032  printFrequency_ = rhs.printFrequency_;
5033# ifdef COIN_HAS_CLP
5034  fastNodeDepth_ = rhs.fastNodeDepth_;
5035#endif
5036  howOftenGlobalScan_ = rhs.howOftenGlobalScan_;
5037  maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
5038  maximumCutPasses_ =  rhs.maximumCutPasses_;
5039  preferredWay_ = rhs.preferredWay_;
5040  resolveAfterTakeOffCuts_ = rhs.resolveAfterTakeOffCuts_;
5041  maximumNumberIterations_ = rhs.maximumNumberIterations_;
5042  numberThreads_ = rhs.numberThreads_;
5043  threadMode_ = rhs.threadMode_;
5044  memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
5045  memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
5046  int i;
5047  if (mode) {
5048    for (i=0;i<numberCutGenerators_;i++) {
5049      delete generator_[i];
5050      delete virginGenerator_[i];
5051    }
5052    delete [] generator_;
5053    delete [] virginGenerator_;
5054    for (i=0;i<numberHeuristics_;i++) {
5055      delete heuristic_[i];
5056    }
5057    delete [] heuristic_;
5058    delete eventHandler_;
5059    delete branchingMethod_;
5060  }
5061  numberCutGenerators_ = rhs.numberCutGenerators_;
5062  if (numberCutGenerators_) {
5063    generator_ = new CbcCutGenerator * [numberCutGenerators_];
5064    virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
5065    int i;
5066    for (i=0;i<numberCutGenerators_;i++) {
5067      if (mode<2)
5068        generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
5069      else
5070        generator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
5071      virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
5072    }
5073  } else {
5074    generator_=NULL;
5075    virginGenerator_=NULL;
5076  }
5077  numberHeuristics_ = rhs.numberHeuristics_;
5078  if (numberHeuristics_) {
5079    heuristic_ = new CbcHeuristic * [numberHeuristics_];
5080    int i;
5081    for (i=0;i<numberHeuristics_;i++) {
5082      heuristic_[i]=rhs.heuristic_[i]->clone();
5083    }
5084  } else {
5085    heuristic_=NULL;
5086  }
5087  if (rhs.eventHandler_)
5088    eventHandler_ = rhs.eventHandler_->clone() ; 
5089  else
5090    eventHandler_ = NULL ; 
5091  if (rhs.branchingMethod_)
5092    branchingMethod_=rhs.branchingMethod_->clone();
5093  else
5094    branchingMethod_=NULL;
5095  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
5096  synchronizeModel();
5097}
5098// Move status, nodes etc etc across
5099void 
5100CbcModel::moveInfo(const CbcModel & rhs)
5101{
5102  bestObjective_ = rhs.bestObjective_;
5103  bestPossibleObjective_=rhs.bestPossibleObjective_;
5104  numberSolutions_=rhs.numberSolutions_;
5105  numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
5106  numberNodes_ = rhs.numberNodes_;
5107  numberNodes2_ = rhs.numberNodes2_;
5108  numberIterations_ = rhs.numberIterations_;
5109  status_ = rhs.status_;
5110  secondaryStatus_ = rhs.secondaryStatus_;
5111  numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
5112  numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
5113  continuousObjective_=rhs.continuousObjective_;
5114  originalContinuousObjective_ = rhs.originalContinuousObjective_;
5115  continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
5116  numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
5117  numberFixedNow_ = rhs.numberFixedNow_;
5118  stoppedOnGap_ = rhs.stoppedOnGap_;
5119  eventHappened_ = rhs.eventHappened_;
5120  numberLongStrong_ = rhs.numberLongStrong_;
5121  numberStrongIterations_ = rhs.numberStrongIterations_;
5122  strongInfo_[0]=rhs.strongInfo_[0];
5123  strongInfo_[1]=rhs.strongInfo_[1];
5124  strongInfo_[2]=rhs.strongInfo_[2];
5125  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
5126  maximumDepth_= rhs.maximumDepth_;
5127}
5128// Save a copy of the current solver so can be reset to
5129void 
5130CbcModel::saveReferenceSolver()
5131{
5132  delete referenceSolver_;
5133  referenceSolver_= solver_->clone();
5134}
5135
5136// Uses a copy of reference solver to be current solver
5137void 
5138CbcModel::resetToReferenceSolver()
5139{
5140  delete solver_;
5141  solver_ = referenceSolver_->clone();
5142  // clear many things
5143  gutsOfDestructor2();
5144  // Reset cutoff
5145  // Solvers know about direction
5146  double direction = solver_->getObjSense();
5147  double value;
5148  solver_->getDblParam(OsiDualObjectiveLimit,value); 
5149  setCutoff(value*direction);
5150}
5151
5152// Are there a numerical difficulties?
5153bool 
5154CbcModel::isAbandoned() const
5155{
5156  return status_ == 2;
5157}
5158// Is optimality proven?
5159bool 
5160CbcModel::isProvenOptimal() const
5161{
5162  if (!status_ && bestObjective_<1.0e30)
5163    return true;
5164  else
5165    return false;
5166}
5167// Is  infeasiblity proven (or none better than cutoff)?
5168bool 
5169CbcModel::isProvenInfeasible() const
5170{
5171  if (!status_ && bestObjective_>=1.0e30)
5172    return true;
5173  else
5174    return false;
5175}
5176// Was continuous solution unbounded
5177bool 
5178CbcModel::isContinuousUnbounded() const
5179{
5180  if (!status_ && secondaryStatus_==7)
5181    return true;
5182  else
5183    return false;
5184}
5185// Was continuous solution unbounded
5186bool 
5187CbcModel::isProvenDualInfeasible() const
5188{
5189  if (!status_ && secondaryStatus_==7)
5190    return true;
5191  else
5192    return false;
5193}
5194// Node limit reached?
5195bool 
5196CbcModel::isNodeLimitReached() const
5197{
5198  return numberNodes_ >= intParam_[CbcMaxNumNode];
5199}
5200// Time limit reached?
5201bool 
5202CbcModel::isSecondsLimitReached() const
5203{
5204  if (status_==1&&secondaryStatus_==4)
5205    return true;
5206  else
5207    return false;
5208}
5209// Solution limit reached?
5210bool 
5211CbcModel::isSolutionLimitReached() const
5212{
5213  return numberSolutions_ >= intParam_[CbcMaxNumSol];
5214}
5215// Set language
5216void 
5217CbcModel::newLanguage(CoinMessages::Language language)
5218{
5219  messages_ = CbcMessage(language);
5220}
5221void 
5222CbcModel::setNumberStrong(int number)
5223{
5224  if (number<0)
5225    numberStrong_=0;
5226   else
5227    numberStrong_=number;
5228}
5229void 
5230CbcModel::setNumberBeforeTrust(int number)
5231{
5232  if (number<-3) {
5233    numberBeforeTrust_=0;
5234  } else {
5235    numberBeforeTrust_=number;
5236    //numberStrong_ = CoinMax(numberStrong_,1);
5237  }
5238}
5239void 
5240CbcModel::setNumberPenalties(int number)
5241{
5242  if (number<=0) {
5243    numberPenalties_=0;
5244  } else {
5245    numberPenalties_=number;
5246  }
5247}
5248void 
5249CbcModel::setPenaltyScaleFactor(double value)
5250{
5251  if (value<=0) {
5252    penaltyScaleFactor_=3.0;
5253  } else {
5254    penaltyScaleFactor_=value;
5255  }
5256}
5257void 
5258CbcModel::setHowOftenGlobalScan(int number)
5259{
5260  if (number<-1)
5261    howOftenGlobalScan_=0;
5262   else
5263    howOftenGlobalScan_=number;
5264}
5265
5266// Add one generator
5267void 
5268CbcModel::addCutGenerator(CglCutGenerator * generator,
5269                          int howOften, const char * name,
5270                          bool normal, bool atSolution,
5271                          bool whenInfeasible,int howOftenInSub,
5272                          int whatDepth, int whatDepthInSub)
5273{
5274  CbcCutGenerator ** temp = generator_;
5275  generator_ = new CbcCutGenerator * [numberCutGenerators_+1];
5276  memcpy(generator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
5277  delete[] temp ;
5278  generator_[numberCutGenerators_]= 
5279    new CbcCutGenerator(this,generator, howOften, name,
5280                        normal,atSolution,whenInfeasible,howOftenInSub,
5281                        whatDepth, whatDepthInSub);
5282  // and before any changes
5283  temp = virginGenerator_;
5284  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_+1];
5285  memcpy(virginGenerator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
5286  delete[] temp ;
5287  virginGenerator_[numberCutGenerators_++]= 
5288    new CbcCutGenerator(this,generator, howOften, name,
5289                        normal,atSolution,whenInfeasible,howOftenInSub,
5290                        whatDepth, whatDepthInSub);
5291                                                         
5292}
5293// Add one heuristic
5294void 
5295CbcModel::addHeuristic(CbcHeuristic * generator, const char *name,
5296                       int before)
5297{
5298  CbcHeuristic ** temp = heuristic_;
5299  heuristic_ = new CbcHeuristic * [numberHeuristics_+1];
5300  memcpy(heuristic_,temp,numberHeuristics_*sizeof(CbcHeuristic *));
5301  delete [] temp;
5302  int where;
5303  if (before<0||before>=numberHeuristics_) {
5304    where=numberHeuristics_;
5305  } else {
5306    // move up
5307    for (int i=numberHeuristics_;i>before;i--) 
5308      heuristic_[i]=heuristic_[i-1];
5309    where=before;
5310  }
5311  heuristic_[where]=generator->clone();
5312  if (name)
5313    heuristic_[where]->setHeuristicName(name) ; 
5314  heuristic_[where]->setSeed(987654321+where);
5315  numberHeuristics_++ ;
5316}
5317
5318/*
5319  The last subproblem handled by the solver is not necessarily related to the
5320  one being recreated, so the first action is to remove all cuts from the
5321  constraint system.  Next, traverse the tree from node to the root to
5322  determine the basis size required for this subproblem and create an empty
5323  basis with the right capacity.  Finally, traverse the tree from root to
5324  node, adjusting bounds in the constraint system, adjusting the basis, and
5325  collecting the cuts that must be added to the constraint system.
5326  applyToModel does the heavy lifting.
5327
5328  addCuts1 is used in contexts where all that's desired is the list of cuts:
5329  the node is already fathomed, and we're collecting cuts so that we can
5330  adjust reference counts as we prune nodes. Arguably the two functions
5331  should be separated. The culprit is applyToModel, which performs cut
5332  collection and model adjustment.
5333
5334  Certainly in the contexts where all we need is a list of cuts, there's no
5335  point in passing in a valid basis --- an empty basis will do just fine.
5336*/
5337void CbcModel::addCuts1 (CbcNode * node, CoinWarmStartBasis *&lastws)
5338{ int i;
5339  int nNode=0;
5340  int numberColumns = getNumCols();
5341  CbcNodeInfo * nodeInfo = node->nodeInfo();
5342
5343/*
5344  Remove all cuts from the constraint system.
5345  (original comment includes ``see note below for later efficiency'', but
5346  the reference isn't clear to me).
5347*/
5348  int currentNumberCuts ;
5349  if (numberThreads_<=0) {
5350    solver_->restoreBaseModel(numberRowsAtContinuous_);
5351  } else {
5352    // *** Fix later
5353    currentNumberCuts = solver_->getNumRows()-numberRowsAtContinuous_;
5354    int *which = new int[currentNumberCuts];
5355    for (i = 0 ; i < currentNumberCuts ; i++)
5356      which[i] = i+numberRowsAtContinuous_;
5357    solver_->deleteRows(currentNumberCuts,which);
5358    delete [] which;
5359  }
5360/*
5361  Accumulate the path from node to the root in walkback_, and accumulate a
5362  cut count in currentNumberCuts.
5363
5364  original comment: when working then just unwind until where new node joins
5365  old node (for cuts?)
5366*/
5367  currentNumberCuts = 0;
5368  while (nodeInfo) {
5369    //printf("nNode = %d, nodeInfo = %x\n",nNode,nodeInfo);
5370    walkback_[nNode++]=nodeInfo;
5371    currentNumberCuts += nodeInfo->numberCuts() ;
5372    nodeInfo = nodeInfo->parent() ;
5373    if (nNode==maximumDepth_) {
5374      maximumDepth_ *= 2;
5375      CbcNodeInfo ** temp = new CbcNodeInfo * [maximumDepth_];
5376      for (i=0;i<nNode;i++) 
5377        temp[i] = walkback_[i];
5378      delete [] walkback_;
5379      walkback_ = temp;
5380    }
5381  }
5382/*
5383  Create an empty basis with sufficient capacity for the constraint system
5384  we'll construct: original system plus cuts. Make sure we have capacity to
5385  record those cuts in addedCuts_.
5386
5387  The method of adjusting the basis at a FullNodeInfo object (the root, for
5388  example) is to use a copy constructor to duplicate the basis held in the
5389  nodeInfo, then resize it and return the new basis object. Guaranteed,
5390  lastws will point to a different basis when it returns. We pass in a basis
5391  because we need the parameter to return the allocated basis, and it's an
5392  easy way to pass in the size. But we take a hit for memory allocation.
5393*/
5394  currentNumberCuts_=currentNumberCuts;
5395  if (currentNumberCuts > maximumNumberCuts_) {
5396    maximumNumberCuts_ = currentNumberCuts;
5397    delete [] addedCuts_;
5398    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
5399  }
5400  lastws->setSize(numberColumns,numberRowsAtContinuous_+currentNumberCuts);
5401/*
5402  This last bit of code traverses the path collected in walkback_ from the
5403  root back to node. At the end of the loop,
5404   * lastws will be an appropriate basis for node;
5405   * variable bounds in the constraint system will be set to be correct for
5406     node; and
5407   * addedCuts_ will be set to a list of cuts that need to be added to the
5408     constraint system at node.
5409  applyToModel does all the heavy lifting.
5410*/
5411  currentNumberCuts=0;
5412  //#define CBC_PRINT2
5413#ifdef CBC_PRINT2
5414  printf("Starting bounds at node %d\n",numberNodes_);
5415#endif
5416  currentDepth_=nNode;
5417  while (nNode) {
5418    --nNode;
5419    walkback_[nNode]->applyToModel(this,lastws,addedCuts_,currentNumberCuts);
5420  }
5421  if (!lastws->fullBasis()) {
5422#ifdef COIN_DEVELOP
5423    printf("******* bad basis\n");
5424#endif
5425    int numberRows = lastws->getNumArtificial();
5426    int i;
5427    for (i=0;i<numberRows;i++)
5428      lastws->setArtifStatus(i,CoinWarmStartBasis::basic);
5429    int numberColumns = lastws->getNumStructural();
5430    for (i=0;i<numberColumns;i++) {
5431      if (lastws->getStructStatus(i)==CoinWarmStartBasis::basic)
5432        lastws->setStructStatus(i,CoinWarmStartBasis::atLowerBound);
5433    }
5434#if 0
5435  } else {
5436    // OPTION - take off slack cuts
5437    // First see if any cuts are slack
5438    int numberAdded = currentNumberCuts;
5439    if (saveNode<2&&false) {
5440      printf("nNode %d cuts %d\n",saveNode,currentNumberCuts);
5441      for (int i=0;i<currentNumberCuts;i++)
5442        addedCuts_[i]->print();
5443    }
5444    if (numberAdded&&saveNode<5&&!parentModel_) {
5445#if 0
5446      currentNumberCuts=0;
5447      for (int j=numberRowsAtContinuous_;
5448           j<numberAdded+numberRowsAtContinuous_;j++) {
5449        CoinWarmStartBasis::Status status = lastws->getArtifStatus(j);
5450        if (status!=CoinWarmStartBasis::basic) {
5451          lastws->setArtifStatus(currentNumberCuts+numberRowsAtContinuous_,
5452                                 status);
5453          addedCuts_[currentNumberCuts++]=addedCuts_[j-numberRowsAtContinuous_];
5454        }
5455      }
5456      if (currentNumberCuts<numberAdded) {
5457        printf("deleting %d rows\n",numberAdded-currentNumberCuts);
5458        lastws->resize(currentNumberCuts+numberRowsAtContinuous_,
5459                       lastws->getNumStructural());
5460        currentNumberCuts_=currentNumberCuts;
5461      }
5462#else
5463      int nDelete=0;
5464      for (int j=numberRowsAtContinuous_;
5465           j<numberAdded+numberRowsAtContinuous_;j++) {
5466        CoinWarmStartBasis::Status status = lastws->getArtifStatus(j);
5467        if (status==CoinWarmStartBasis::basic)
5468          nDelete++;
5469      }
5470      if (nDelete)
5471        printf("depth %d can delete %d\n",saveNode-1,nDelete);
5472#endif
5473    }
5474#endif
5475  }
5476  if (0) {
5477    int numberDebugValues=18;
5478    double * debugValues = new double[numberDebugValues];
5479    CoinZeroN(debugValues,numberDebugValues);
5480    debugValues[1]=6.0;
5481    debugValues[3]=60.0;
5482    debugValues[4]=6.0;
5483    debugValues[6]=60.0;
5484    debugValues[7]=16.0;
5485    debugValues[9]=70.0;
5486    debugValues[10]=7.0;
5487    debugValues[12]=70.0;
5488    debugValues[13]=12.0;
5489    debugValues[15]=75.0;
5490    int nBad=0;
5491    for (int j=0;j<numberColumns;j++) {
5492      if (integerInfo_[j]) {
5493        if(solver_->getColLower()[j]>debugValues[j]||
5494           solver_->getColUpper()[j]<debugValues[j]) {
5495          printf("** (%g) ** ",debugValues[j]);
5496          nBad++;
5497        }
5498        printf("%d bounds %g %g\n",j,solver_->getColLower()[j],solver_->getColUpper()[j]);
5499      }
5500    }
5501    if (nBad)
5502      printf("%d BAD\n",nBad);
5503    else
5504      printf("OKAY\n");
5505    delete [] debugValues;
5506  }
5507}
5508
5509/*
5510  adjustCuts might be a better name: If the node is feasible, we sift through
5511  the cuts collected by addCuts1, add the ones that are tight and omit the
5512  ones that are loose. If the node is infeasible, we just adjust the
5513  reference counts to reflect that we're about to prune this node and its
5514  descendants.
5515*/
5516int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws,bool canFix)
5517{
5518/*
5519  addCuts1 performs step 1 of restoring the subproblem at this node; see the
5520  comments there.
5521*/
5522  addCuts1(node,lastws);
5523  int i;
5524  int numberColumns = getNumCols();
5525  if (solver_->getNumRows()>maximumRows_) {
5526    maximumRows_ = solver_->getNumRows();
5527    workingBasis_.resize(maximumRows_,numberColumns);
5528  }
5529  CbcNodeInfo * nodeInfo = node->nodeInfo();
5530  double cutoff = getCutoff() ;
5531  int currentNumberCuts=currentNumberCuts_;
5532  if (canFix) {
5533    bool feasible=true;
5534    const double *lower = solver_->getColLower() ;
5535    const double *upper = solver_->getColUpper() ;
5536    double * newLower = analyzeResults_;
5537    double * objLower = newLower+numberIntegers_;
5538    double * newUpper = objLower+numberIntegers_;
5539    double * objUpper = newUpper+numberIntegers_;
5540    int n=0;
5541    for (i=0;i<numberIntegers_;i++) {
5542      int iColumn = integerVariable_[i];
5543      bool changed=false;
5544      double lo = 0.0;
5545      double up = 0.0;
5546      if (objLower[i]>cutoff) {
5547        lo = lower[iColumn];
5548        up = upper[iColumn];
5549        if (lo<newLower[i]) {
5550          lo = newLower[i];
5551          solver_->setColLower(iColumn,lo);
5552          changed=true;
5553          n++;
5554        }
5555        if (objUpper[i]>cutoff) {
5556          if (up>newUpper[i]) {
5557            up = newUpper[i];
5558            solver_->setColUpper(iColumn,up);
5559            changed=true;
5560            n++;
5561          }
5562        }
5563      } else if (objUpper[i]>cutoff) {
5564        lo = lower[iColumn];
5565        up = upper[iColumn];
5566        if (up>newUpper[i]) {
5567          up = newUpper[i];
5568          solver_->setColUpper(iColumn,up);
5569          changed=true;
5570          n++;
5571        }
5572      }
5573      if (changed&&lo>up) {
5574        feasible=false;
5575        break;
5576      }
5577    }
5578    if (!feasible) {
5579      printf("analysis says node infeas\n");
5580      cutoff=-COIN_DBL_MAX;
5581    }
5582  }
5583/*
5584  If the node can't be fathomed by bound, reinstall tight cuts in the
5585  constraint system. Even if there are no cuts, we'll want to set the
5586  reconstructed basis in the solver.
5587*/
5588  if (node->objectiveValue() < cutoff||numberThreads_)
5589  { 
5590    //#   define CBC_CHECK_BASIS
5591#   ifdef CBC_CHECK_BASIS
5592    printf("addCuts: expanded basis; rows %d+%d\n",
5593           numberRowsAtContinuous_,currentNumberCuts);
5594    lastws->print();
5595#   endif
5596/*
5597  Adjust the basis and constraint system so that we retain only active cuts.
5598  There are three steps:
5599    1) Scan the basis. Sort the cuts into effective cuts to be kept and
5600       loose cuts to be dropped.
5601    2) Drop the loose cuts and resize the basis to fit.
5602    3) Install the tight cuts in the constraint system (applyRowCuts) and
5603       and install the basis (setWarmStart).
5604  Use of compressRows conveys we're compressing the basis and not just
5605  tweaking the artificialStatus_ array.
5606*/
5607    if (currentNumberCuts > 0) {
5608      int numberToAdd = 0;
5609      const OsiRowCut **addCuts;
5610      int numberToDrop = 0 ;
5611      int *cutsToDrop ;
5612      addCuts = new const OsiRowCut* [currentNumberCuts];
5613      cutsToDrop = new int[currentNumberCuts] ;
5614      assert (currentNumberCuts+numberRowsAtContinuous_<=lastws->getNumArtificial());
5615      for (i=0;i<currentNumberCuts;i++) {
5616        CoinWarmStartBasis::Status status = 
5617          lastws->getArtifStatus(i+numberRowsAtContinuous_);
5618        if (addedCuts_[i] &&
5619            (status != CoinWarmStartBasis::basic ||
5620             (addedCuts_[i]->effectiveness()>1.0e10&&
5621              !addedCuts_[i]->canDropCut(solver_,i+numberRowsAtContinuous_)))) {
5622#         ifdef CHECK_CUT_COUNTS
5623          printf("Using cut %d %x as row %d\n",i,addedCuts_[i],
5624                 numberRowsAtContinuous_+numberToAdd);
5625#         endif
5626          addCuts[numberToAdd++] = addedCuts_[i];
5627        } else {
5628#         ifdef CHECK_CUT_COUNTS
5629          printf("Dropping cut %d %x\n",i,addedCuts_[i]);
5630#         endif
5631          addedCuts_[i]=NULL;
5632          cutsToDrop[numberToDrop++] = numberRowsAtContinuous_+i ;
5633        }
5634      }
5635      int numberRowsNow=numberRowsAtContinuous_+numberToAdd;
5636      lastws->compressRows(numberToDrop,cutsToDrop) ;
5637      lastws->resize(numberRowsNow,numberColumns);
5638      solver_->applyRowCuts(numberToAdd,addCuts);
5639#     ifdef CBC_CHECK_BASIS
5640      printf("addCuts: stripped basis; rows %d + %d\n",
5641             numberRowsAtContinuous_,numberToAdd);
5642      lastws->print();
5643#     endif
5644      //for (i=0;i<numberToAdd;i++)
5645      //delete addCuts[i];
5646      delete [] addCuts;
5647      delete [] cutsToDrop ;
5648    }
5649/*
5650  Set the basis in the solver.
5651*/
5652    solver_->setWarmStart(lastws);
5653#if 0
5654    if ((numberNodes_%printFrequency_)==0) {
5655      printf("Objective %g, depth %d, unsatisfied %d\n",
5656             node->objectiveValue(),
5657             node->depth(),node->numberUnsatisfied());
5658    }
5659#endif
5660/*
5661  Clean up and we're out of here.
5662*/
5663    numberNodes_++;
5664    return 0;
5665  } 
5666/*
5667  This node has been fathomed by bound as we try to revive it out of the live
5668  set. Adjust the cut reference counts to reflect that we no longer need to
5669  explore the remaining branch arms, hence they will no longer reference any
5670  cuts. Cuts whose reference count falls to zero are deleted. 
5671*/
5672  else
5673  { int i;
5674    if (currentNumberCuts) {
5675#ifndef CBC_DETERMINISTIC_THREAD
5676      lockThread();
5677#endif
5678      int numberLeft = nodeInfo->numberBranchesLeft();
5679      for (i = 0 ; i < currentNumberCuts ; i++)
5680        { if (addedCuts_[i])
5681          { if (!addedCuts_[i]->decrement(numberLeft))
5682            { delete addedCuts_[i];
5683            addedCuts_[i] = NULL; } } }
5684#ifndef CBC_DETERMINISTIC_THREAD
5685      unlockThread();
5686#endif
5687    }
5688    return 1 ; }
5689}
5690
5691
5692/*
5693  Perform reduced cost fixing on integer variables.
5694
5695  The variables in question are already nonbasic at bound. We're just nailing
5696  down the current situation.
5697*/
5698#if 0
5699int CbcModel::reducedCostFix ()
5700
5701{
5702  if(!solverCharacteristics_->reducedCostsAccurate())
5703    return 0; //NLP
5704  double cutoff = getCutoff() ;
5705  double direction = solver_->getObjSense() ;
5706  double gap = cutoff - solver_->getObjValue()*direction ;
5707  double tolerance;
5708  solver_->getDblParam(OsiDualTolerance,tolerance) ;
5709  if (gap<=0.0)
5710    gap = tolerance; //return 0;
5711  gap += 100.0*tolerance;
5712  double integerTolerance = getDblParam(CbcIntegerTolerance) ;
5713
5714  const double *lower = solver_->getColLower() ;
5715  const double *upper = solver_->getColUpper() ;
5716  const double *solution = solver_->getColSolution() ;
5717  const double *reducedCost = solver_->getReducedCost() ;
5718
5719  int numberFixed = 0 ;
5720
5721# ifdef COIN_HAS_CLP
5722  OsiClpSolverInterface * clpSolver
5723    = dynamic_cast<OsiClpSolverInterface *> (solver_);
5724  ClpSimplex * clpSimplex=NULL;
5725  if (clpSolver)
5726    clpSimplex = clpSolver->getModelPtr();
5727# endif
5728  for (int i = 0 ; i < numberIntegers_ ; i++) {
5729    int iColumn = integerVariable_[i] ;
5730    double djValue = direction*reducedCost[iColumn] ;
5731    if (upper[iColumn]-lower[iColumn] > integerTolerance) {
5732      if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap) {
5733#ifdef COIN_HAS_CLP
5734        // may just have been fixed before
5735        if (clpSimplex) {
5736          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
5737#ifdef COIN_DEVELOP
5738            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
5739                   iColumn,clpSimplex->getColumnStatus(iColumn),
5740                   djValue,gap,lower[iColumn],upper[iColumn]);
5741#endif
5742          } else {         
5743            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
5744                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
5745          }
5746        }
5747#endif
5748        solver_->setColUpper(iColumn,lower[iColumn]) ;
5749        numberFixed++ ;
5750      } else if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap) {
5751#ifdef COIN_HAS_CLP
5752        // may just have been fixed before
5753        if (clpSimplex) {
5754          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
5755#ifdef COIN_DEVELOP
5756            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
5757                   iColumn,clpSimplex->getColumnStatus(iColumn),
5758                   djValue,gap,lower[iColumn],upper[iColumn]);
5759#endif
5760          } else {         
5761            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
5762                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
5763          }
5764        }
5765#endif
5766        solver_->setColLower(iColumn,upper[iColumn]) ;
5767        numberFixed++ ;
5768      }
5769    }
5770  }
5771  numberDJFixed_ += numberFixed;
5772  return numberFixed; }
5773#else
5774int CbcModel::reducedCostFix ()
5775
5776{
5777  if(!solverCharacteristics_->reducedCostsAccurate())
5778    return 0; //NLP
5779  double cutoff = getCutoff() ;
5780  double direction = solver_->getObjSense() ;
5781  double gap = cutoff - solver_->getObjValue()*direction ;
5782  double tolerance;
5783  solver_->getDblParam(OsiDualTolerance,tolerance) ;
5784  if (gap<=0.0)
5785    gap = tolerance; //return 0;
5786  gap += 100.0*tolerance;
5787  double integerTolerance = getDblParam(CbcIntegerTolerance) ;
5788
5789  const double *lower = solver_->getColLower() ;
5790  const double *upper = solver_->getColUpper() ;
5791  const double *solution = solver_->getColSolution() ;
5792  const double *reducedCost = solver_->getReducedCost() ;
5793
5794  int numberFixed = 0 ;
5795  int numberTightened = 0 ;
5796
5797# ifdef COIN_HAS_CLP
5798  OsiClpSolverInterface * clpSolver
5799    = dynamic_cast<OsiClpSolverInterface *> (solver_);
5800  ClpSimplex * clpSimplex=NULL;
5801  if (clpSolver)
5802    clpSimplex = clpSolver->getModelPtr();
5803# endif
5804  for (int i = 0 ; i < numberIntegers_ ; i++) {
5805    int iColumn = integerVariable_[i] ;
5806    double djValue = direction*reducedCost[iColumn] ;
5807    double boundGap=upper[iColumn]-lower[iColumn];
5808    if (boundGap > integerTolerance) {
5809      if (solution[iColumn] < lower[iColumn]+integerTolerance
5810          && djValue*boundGap > gap) {
5811#ifdef COIN_HAS_CLP
5812        // may just have been fixed before
5813        if (clpSimplex) {
5814          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
5815#ifdef COIN_DEVELOP
5816            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
5817                   iColumn,clpSimplex->getColumnStatus(iColumn),
5818                   djValue,gap,lower[iColumn],upper[iColumn]);
5819#endif
5820          } else {         
5821            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
5822                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
5823          }
5824        }
5825#endif
5826        double newBound=lower[iColumn];
5827        if (boundGap>1.99) {
5828          boundGap = gap/djValue+1.0e-4*boundGap;
5829          newBound=lower[iColumn]+floor(boundGap);
5830          numberTightened++;
5831          //if (newBound)
5832          //printf("tighter - gap %g dj %g newBound %g\n",
5833          //   gap,djValue,newBound);
5834        }
5835        solver_->setColUpper(iColumn,newBound) ;
5836        numberFixed++ ;
5837      } else if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > boundGap*gap) {
5838#ifdef COIN_HAS_CLP
5839        // may just have been fixed before
5840        if (clpSimplex) {
5841          if (clpSimplex->getColumnStatus(iColumn)==ClpSimplex::basic) {
5842#ifdef COIN_DEVELOP
5843            printf("DJfix %d has status of %d, dj of %g gap %g, bounds %g %g\n",
5844                   iColumn,clpSimplex->getColumnStatus(iColumn),
5845                   djValue,gap,lower[iColumn],upper[iColumn]);
5846#endif
5847          } else {         
5848            assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
5849                   clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
5850          }
5851        }
5852#endif
5853        double newBound=upper[iColumn];
5854        if (boundGap>1.99) {
5855          boundGap = -gap/djValue+1.0e-4*boundGap;
5856          newBound=upper[iColumn]-floor(boundGap);
5857          //if (newBound)
5858          //printf("tighter - gap %g dj %g newBound %g\n",
5859          //   gap,djValue,newBound);
5860          numberTightened++;
5861        }
5862        solver_->setColLower(iColumn,newBound) ;
5863        numberFixed++ ;
5864      }
5865    }
5866  }
5867  numberDJFixed_ += numberFixed-numberTightened;
5868#ifdef COIN_DEVELOP
5869  //if (numberTightened)
5870  //printf("%d tightened, %d others fixed\n",numberTightened,
5871  //   numberFixed-numberTightened);
5872#endif
5873  return numberFixed; }
5874#endif
5875// Collect coding to replace whichGenerator
5876void
5877CbcModel::resizeWhichGenerator(int numberNow, int numberAfter)
5878{
5879  if (numberAfter > maximumWhich_) {
5880    maximumWhich_ = CoinMax(maximumWhich_*2+100,numberAfter) ;
5881    int * temp = new int[2*maximumWhich_] ;
5882    memcpy(temp,whichGenerator_,numberNow*sizeof(int)) ;
5883    delete [] whichGenerator_ ;
5884    whichGenerator_ = temp ;
5885    memset(whichGenerator_+numberNow,0,(maximumWhich_-numberNow)*sizeof(int));
5886  }
5887}
5888
5889/** Solve the model using cuts
5890
5891  This version takes off redundant cuts from node.
5892  Returns true if feasible.
5893
5894  \todo
5895  Why do I need to resolve the problem? What has been done between the last
5896  relaxation and calling solveWithCuts?
5897
5898  If numberTries == 0 then user did not want any cuts.
5899*/
5900
5901bool 
5902CbcModel::solveWithCuts (OsiCuts &cuts, int numberTries, CbcNode *node)
5903/*
5904  Parameters:
5905    numberTries: (i) the maximum number of iterations for this round of cut
5906                     generation; if negative then we don't mind if drop is tiny.
5907   
5908    cuts:       (o) all cuts generated in this round of cut generation
5909
5910    node: (i)     So we can update dynamic pseudo costs
5911*/
5912                       
5913
5914{
5915# ifdef COIN_HAS_CLP
5916  O