source: stable/2.1/Cbc/src/CbcModel.cpp @ 943

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

so ctrl c will give correct result

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