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

Last change on this file since 837 was 837, checked in by forrest, 13 years ago

memory leak

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