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

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

allow dense solver

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