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

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

null currentNode_

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