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

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

for threads and hotstart

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