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

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

OA bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 422.3 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,true);
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      // make sure handler will be deleted
2997      threadModel[i]->defaultHandler_=true;
2998      delete threadModel[i];
2999    }
3000    delete [] mutex2;
3001    delete [] condition2;
3002    delete [] threadId;
3003    delete [] threadInfo;
3004    delete [] threadModel;
3005    delete [] threadCount;
3006    mutex_=NULL;
3007    // adjust time to allow for children on some systems
3008    dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
3009  }
3010#endif
3011/*
3012  End of the non-abort actions. The next block of code is executed if we've
3013  aborted because we hit one of the limits. Clean up by deleting the live set
3014  and break out of the node processing loop. Note that on an abort, node may
3015  have been pushed back onto the tree for further processing, in which case
3016  it'll be deleted in cleanTree. We need to check.
3017*/
3018    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
3019        numberSolutions_ < intParam_[CbcMaxNumSol] &&
3020        totalTime < dblParam_[CbcMaximumSeconds] &&
3021        !stoppedOnGap_&&!eventHappened_)) {
3022      if (tree_->size())
3023        tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
3024      delete nextRowCut_;
3025      if (stoppedOnGap_)
3026        { messageHandler()->message(CBC_GAP,messages())
3027          << bestObjective_-bestPossibleObjective_
3028          << dblParam_[CbcAllowableGap]
3029          << dblParam_[CbcAllowableFractionGap]*100.0
3030          << CoinMessageEol ;
3031        secondaryStatus_ = 2;
3032        status_ = 0 ; }
3033        else
3034          if (isNodeLimitReached())
3035            { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
3036            secondaryStatus_ = 3;
3037            status_ = 1 ; }
3038          else
3039        if (totalTime >= dblParam_[CbcMaximumSeconds])
3040          { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
3041          secondaryStatus_ = 4;
3042          status_ = 1 ; }
3043        else
3044          if (eventHappened_)
3045            { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
3046            secondaryStatus_ = 5;
3047            status_ = 5 ; }
3048          else
3049            { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
3050            secondaryStatus_ = 6;
3051            status_ = 1 ; }
3052    }
3053/*
3054  That's it, we've exhausted the search tree, or broken out of the loop because
3055  we hit some limit on evaluation.
3056
3057  We may have got an intelligent tree so give it one more chance
3058*/
3059  // Tell solver we are not in Branch and Cut
3060  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
3061  tree_->endSearch();
3062  //  If we did any sub trees - did we give up on any?
3063  if ( numberStoppedSubTrees_)
3064    status_=1;
3065  if (!status_) {
3066    // Set best possible unless stopped on gap
3067    if(secondaryStatus_ != 2)
3068      bestPossibleObjective_=bestObjective_;
3069    handler_->message(CBC_END_GOOD,messages_)
3070      << bestObjective_ << numberIterations_ << numberNodes_<<getCurrentSeconds()
3071      << CoinMessageEol ;
3072  } else {
3073    handler_->message(CBC_END,messages_)
3074      << bestObjective_ <<bestPossibleObjective_
3075      << numberIterations_ << numberNodes_<<getCurrentSeconds()
3076      << CoinMessageEol ;
3077  }
3078  if (numberStrongIterations_)
3079    handler_->message(CBC_STRONG_STATS,messages_)
3080      << strongInfo_[0] << numberStrongIterations_ << strongInfo_[2]
3081      << strongInfo_[1] << CoinMessageEol ;
3082  handler_->message(CBC_OTHER_STATS,messages_)
3083      << maximumDepthActual_
3084      << numberDJFixed_ << CoinMessageEol ;
3085  if (doStatistics==100) {
3086    for (int i=0;i<numberObjects_;i++) {
3087      CbcSimpleIntegerDynamicPseudoCost * obj =
3088        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
3089      if (obj)
3090        obj->print();
3091    }
3092  }
3093  if (statistics_) {
3094    // report in some way
3095    int * lookup = new int[numberObjects_];
3096    int i;
3097    for (i=0;i<numberObjects_;i++) 
3098      lookup[i]=-1;
3099    bool goodIds=true;
3100    for (i=0;i<numberObjects_;i++) {
3101      int iColumn = object_[i]->columnNumber();
3102      if(iColumn>=0&&iColumn<numberColumns) {
3103        if (lookup[i]==-1) {
3104          lookup[i]=iColumn;
3105        } else {
3106          goodIds=false;
3107          break;
3108        }
3109      } else {
3110        goodIds=false;
3111        break;
3112      }
3113    }
3114    if (!goodIds) {
3115      delete [] lookup;
3116      lookup=NULL;
3117    }
3118    if (doStatistics==3) {
3119      printf("  node parent depth column   value                    obj      inf\n");
3120      for ( i=0;i<numberNodes2_;i++) {
3121        statistics_[i]->print(lookup);
3122      }
3123    }
3124    if (doStatistics>1) {
3125      // Find last solution
3126      int k;
3127      for (k=numberNodes2_-1;k>=0;k--) {
3128        if (statistics_[k]->endingObjective()!=COIN_DBL_MAX&&
3129            !statistics_[k]->endingInfeasibility())
3130          break;
3131      }
3132      if (k>=0) {
3133        int depth=statistics_[k]->depth();
3134        int * which = new int[depth+1];
3135        for (i=depth;i>=0;i--) {
3136          which[i]=k;
3137          k=statistics_[k]->parentNode();
3138        }
3139        printf("  node parent depth column   value                    obj      inf\n");
3140        for (i=0;i<=depth;i++) {
3141          statistics_[which[i]]->print(lookup);
3142        }
3143        delete [] which;
3144      }
3145    }
3146    // now summary
3147    int maxDepth=0;
3148    double averageSolutionDepth=0.0;
3149    int numberSolutions=0;
3150    double averageCutoffDepth=0.0;
3151    double averageSolvedDepth=0.0;
3152    int numberCutoff=0;
3153    int numberDown=0;
3154    int numberFirstDown=0;
3155    double averageInfDown=0.0;
3156    double averageObjDown=0.0;
3157    int numberCutoffDown=0;
3158    int numberUp=0;
3159    int numberFirstUp=0;
3160    double averageInfUp=0.0;
3161    double averageObjUp=0.0;
3162    int numberCutoffUp=0;
3163    double averageNumberIterations1=0.0;
3164    double averageValue=0.0;
3165    for ( i=0;i<numberNodes2_;i++) {
3166      int depth =  statistics_[i]->depth(); 
3167      int way =  statistics_[i]->way(); 
3168      double value = statistics_[i]->value(); 
3169      double startingObjective =  statistics_[i]->startingObjective(); 
3170      int startingInfeasibility = statistics_[i]->startingInfeasibility(); 
3171      double endingObjective = statistics_[i]->endingObjective(); 
3172      int endingInfeasibility = statistics_[i]->endingInfeasibility(); 
3173      maxDepth = CoinMax(depth,maxDepth);
3174      // Only for completed
3175      averageNumberIterations1 += statistics_[i]->numberIterations();
3176      averageValue += value;
3177      if (endingObjective!=COIN_DBL_MAX&&!endingInfeasibility) {
3178        numberSolutions++;
3179        averageSolutionDepth += depth;
3180      }
3181      if (endingObjective==COIN_DBL_MAX) {
3182        numberCutoff++;
3183        averageCutoffDepth += depth;
3184        if (way<0) {
3185          numberDown++;
3186          numberCutoffDown++;
3187          if (way==-1)
3188            numberFirstDown++;
3189        } else {
3190          numberUp++;
3191          numberCutoffUp++;
3192          if (way==1)
3193            numberFirstUp++;
3194        }
3195      } else {
3196        averageSolvedDepth += depth;
3197        if (way<0) {
3198          numberDown++;
3199          averageInfDown += startingInfeasibility-endingInfeasibility;
3200          averageObjDown += endingObjective-startingObjective;
3201          if (way==-1)
3202            numberFirstDown++;
3203        } else {
3204          numberUp++;
3205          averageInfUp += startingInfeasibility-endingInfeasibility;
3206          averageObjUp += endingObjective-startingObjective;
3207          if (way==1)
3208            numberFirstUp++;
3209        }
3210      }
3211    }
3212    // Now print
3213    if (numberSolutions)
3214      averageSolutionDepth /= (double) numberSolutions;
3215    int numberSolved = numberNodes2_-numberCutoff;
3216    double averageNumberIterations2=numberIterations_-averageNumberIterations1
3217      -numberIterationsAtContinuous;
3218    if(numberCutoff) {
3219      averageCutoffDepth /= (double) numberCutoff;
3220      averageNumberIterations2 /= (double) numberCutoff;
3221    }
3222    if (numberNodes2_) 
3223      averageValue /= (double) numberNodes2_;
3224    if (numberSolved) {
3225      averageNumberIterations1 /= (double) numberSolved;
3226      averageSolvedDepth /= (double) numberSolved;
3227    }
3228    printf("%d solution(s) were found (by branching) at an average depth of %g\n",
3229           numberSolutions,averageSolutionDepth);
3230    printf("average value of variable being branched on was %g\n",
3231           averageValue);
3232    printf("%d nodes were cutoff at an average depth of %g with iteration count of %g\n",
3233           numberCutoff,averageCutoffDepth,averageNumberIterations2);
3234    printf("%d nodes were solved at an average depth of %g with iteration count of %g\n",
3235           numberSolved,averageSolvedDepth,averageNumberIterations1);
3236    if (numberDown) {
3237      averageInfDown /= (double) numberDown;
3238      averageObjDown /= (double) numberDown;
3239    }
3240    printf("Down %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
3241           numberDown,numberFirstDown,numberDown-numberFirstDown,numberCutoffDown,
3242           averageInfDown,averageObjDown);
3243    if (numberUp) {
3244      averageInfUp /= (double) numberUp;
3245      averageObjUp /= (double) numberUp;
3246    }
3247    printf("Up %d nodes (%d first, %d second) - %d cutoff, rest decrease numinf %g increase obj %g\n",
3248           numberUp,numberFirstUp,numberUp-numberFirstUp,numberCutoffUp,
3249           averageInfUp,averageObjUp);
3250    for ( i=0;i<numberNodes2_;i++) 
3251      delete statistics_[i];
3252    delete [] statistics_;
3253    statistics_=NULL;
3254    maximumStatistics_=0;
3255    delete [] lookup;
3256  }
3257/*
3258  If we think we have a solution, restore and confirm it with a call to
3259  setBestSolution().  We need to reset the cutoff value so as not to fathom
3260  the solution on bounds.  Note that calling setBestSolution( ..., true)
3261  leaves the continuousSolver_ bounds vectors fixed at the solution value.
3262
3263  Running resolve() here is a failsafe --- setBestSolution has already
3264  reoptimised using the continuousSolver_. If for some reason we fail to
3265  prove optimality, run the problem again after instructing the solver to
3266  tell us more.
3267
3268  If all looks good, replace solver_ with continuousSolver_, so that the
3269  outside world will be able to obtain information about the solution using
3270  public methods.
3271*/
3272  if (bestSolution_&&(solverCharacteristics_->solverType()<2||solverCharacteristics_->solverType()==4)) 
3273  { setCutoff(1.0e50) ; // As best solution should be worse than cutoff
3274    phase_=5;
3275    double increment = getDblParam(CbcModel::CbcCutoffIncrement) ;
3276    if ((specialOptions_&4)==0)
3277      bestObjective_ += 100.0*increment+1.0e-3; // only set if we are going to solve
3278    setBestSolution(CBC_END_SOLUTION,bestObjective_,bestSolution_,true) ;
3279    continuousSolver_->resolve() ;
3280    if (!continuousSolver_->isProvenOptimal())
3281    { continuousSolver_->messageHandler()->setLogLevel(2) ;
3282      continuousSolver_->initialSolve() ; }
3283    delete solver_ ;
3284    // above deletes solverCharacteristics_
3285    solverCharacteristics_ = NULL;
3286    solver_ = continuousSolver_ ;
3287    setPointers(solver_);
3288    continuousSolver_ = NULL ; }
3289/*
3290  Clean up dangling objects. continuousSolver_ may already be toast.
3291*/
3292  delete lastws ;
3293  if (saveObjects) {
3294    for (int i=0;i<numberObjects_;i++)
3295      delete saveObjects[i];
3296    delete [] saveObjects;
3297  }
3298  delete [] whichGenerator_ ;
3299  whichGenerator_=NULL;
3300  delete [] lowerBefore ;
3301  delete [] upperBefore ;
3302  delete [] walkback_ ;
3303  walkback_ = NULL ;
3304  delete [] addedCuts_ ;
3305  addedCuts_ = NULL ;
3306  //delete persistentInfo;
3307  // Get rid of characteristics
3308  solverCharacteristics_=NULL;
3309  if (continuousSolver_)
3310  { delete continuousSolver_ ;
3311    continuousSolver_ = NULL ; }
3312/*
3313  Destroy global cuts by replacing with an empty OsiCuts object.
3314*/
3315  globalCuts_= OsiCuts() ;
3316  if (!bestSolution_) {
3317    // make sure lp solver is infeasible
3318    int numberColumns = solver_->getNumCols();
3319    const double * columnLower = solver_->getColLower();
3320    int iColumn;
3321    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3322      if (solver_->isInteger(iColumn)) 
3323        solver_->setColUpper(iColumn,columnLower[iColumn]);
3324    }
3325    solver_->initialSolve();
3326  }
3327  if (strategy_&&strategy_->preProcessState()>0) {
3328    // undo preprocessing
3329    CglPreProcess * process = strategy_->process();
3330    assert (process);
3331    int n = originalSolver->getNumCols();
3332    if (bestSolution_) {
3333      delete [] bestSolution_;
3334      bestSolution_ = new double [n];
3335      process->postProcess(*solver_);
3336    }
3337    strategy_->deletePreProcess();
3338    // Solution now back in originalSolver
3339    delete solver_;
3340    solver_=originalSolver;
3341    if (bestSolution_)
3342      memcpy(bestSolution_,solver_->getColSolution(),n*sizeof(double));
3343    // put back original objects if there were any
3344    if (originalObject) {
3345      int iColumn;
3346      assert (ownObjects_);
3347      for (iColumn=0;iColumn<numberObjects_;iColumn++) 
3348        delete object_[iColumn];
3349      delete [] object_;
3350      numberObjects_ = numberOriginalObjects;
3351      object_=originalObject;
3352      delete [] integerVariable_;
3353      numberIntegers_=0;
3354      for (iColumn=0;iColumn<n;iColumn++) {
3355        if (solver_->isInteger(iColumn))
3356          numberIntegers_++;
3357      }
3358      integerVariable_ = new int[numberIntegers_];
3359      numberIntegers_=0;
3360      for (iColumn=0;iColumn<n;iColumn++) {
3361        if (solver_->isInteger(iColumn))
3362            integerVariable_[numberIntegers_++]=iColumn;
3363      }
3364    }
3365  }
3366#ifdef CLP_QUICK_OPTIONS
3367  {
3368    OsiClpSolverInterface * clpSolver
3369      = dynamic_cast<OsiClpSolverInterface *> (solver_);
3370    if (clpSolver) {
3371      // Try and re-use regions   
3372      ClpSimplex * simplex = clpSolver->getModelPtr();
3373      simplex->setPersistenceFlag(0);
3374      clpSolver->deleteScaleFactors();
3375      clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~131072));
3376      simplex->setSpecialOptions(simplex->specialOptions()&(~131072));
3377      simplex->setAlphaAccuracy(-1.0);
3378      //clpSolver->setSpecialOptions((clpSolver->specialOptions()&~128)|65536);
3379    }
3380  }
3381#endif
3382  return ;
3383 }
3384
3385
3386// Solve the initial LP relaxation
3387void 
3388CbcModel::initialSolve() 
3389{
3390  assert (solver_);
3391  assert (!solverCharacteristics_);
3392  OsiBabSolver * solverCharacteristics = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
3393  if (solverCharacteristics) {
3394    solverCharacteristics_ = solverCharacteristics;
3395  } else {
3396    // replace in solver
3397    OsiBabSolver defaultC;
3398    solver_->setAuxiliaryInfo(&defaultC);
3399    solverCharacteristics_ = dynamic_cast<OsiBabSolver *> (solver_->getAuxiliaryInfo());
3400  }
3401  solverCharacteristics_->setSolver(solver_);
3402  solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
3403  solver_->initialSolve();
3404  solver_->setHintParam(OsiDoInBranchAndCut,false,OsiHintDo,NULL) ;
3405  // But set up so Jon Lee will be happy
3406  status_=-1;
3407  secondaryStatus_ = -1;
3408  originalContinuousObjective_ = solver_->getObjValue()*solver_->getObjSense();
3409  delete [] continuousSolution_;
3410  continuousSolution_ = CoinCopyOfArray(solver_->getColSolution(),
3411                                             solver_->getNumCols());
3412  setPointers(solver_);
3413  solverCharacteristics_ = NULL;
3414}
3415
3416/*! \brief Get an empty basis object
3417
3418  Return an empty CoinWarmStartBasis object with the requested capacity,
3419  appropriate for the current solver. The object is cloned from the object
3420  cached as emptyWarmStart_. If there is no cached object, the routine
3421  queries the solver for a warm start object, empties it, and caches the
3422  result.
3423*/
3424
3425CoinWarmStartBasis *CbcModel::getEmptyBasis (int ns, int na) const
3426
3427{ CoinWarmStartBasis *emptyBasis ;
3428/*
3429  Acquire an empty basis object, if we don't yet have one.
3430*/
3431  if (emptyWarmStart_ == 0)
3432  { if (solver_ == 0)
3433    { throw CoinError("Cannot construct basis without solver!",
3434                      "getEmptyBasis","CbcModel") ; }
3435    emptyBasis =
3436        dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
3437    if (emptyBasis == 0)
3438    { throw CoinError(
3439        "Solver does not appear to use a basis-oriented warm start.",
3440        "getEmptyBasis","CbcModel") ; }
3441    emptyBasis->setSize(0,0) ;
3442    emptyWarmStart_ = dynamic_cast<CoinWarmStart *>(emptyBasis) ; }
3443/*
3444  Clone the empty basis object, resize it as requested, and return.
3445*/
3446  emptyBasis = dynamic_cast<CoinWarmStartBasis *>(emptyWarmStart_->clone()) ;
3447  assert(emptyBasis) ;
3448  if (ns != 0 || na != 0) emptyBasis->setSize(ns,na) ;
3449
3450  return (emptyBasis) ; }
3451   
3452
3453/** Default Constructor
3454
3455  Creates an empty model without an associated solver.
3456*/
3457CbcModel::CbcModel() 
3458
3459:
3460  solver_(NULL),
3461  ownership_(0x80000000),
3462  continuousSolver_(NULL),
3463  referenceSolver_(NULL),
3464  defaultHandler_(true),
3465  emptyWarmStart_(NULL),
3466  bestObjective_(COIN_DBL_MAX),
3467  bestPossibleObjective_(COIN_DBL_MAX),
3468  sumChangeObjective1_(0.0),
3469  sumChangeObjective2_(0.0),
3470  bestSolution_(NULL),
3471  currentSolution_(NULL),
3472  testSolution_(NULL),
3473  minimumDrop_(1.0e-4),
3474  numberSolutions_(0),
3475  stateOfSearch_(0),
3476  hotstartSolution_(NULL),
3477  hotstartPriorities_(NULL),
3478  numberHeuristicSolutions_(0),
3479  numberNodes_(0),
3480  numberNodes2_(0),
3481  numberIterations_(0),
3482  status_(-1),
3483  secondaryStatus_(-1),
3484  numberIntegers_(0),
3485  numberRowsAtContinuous_(0),
3486  maximumNumberCuts_(0),
3487  phase_(0),
3488  currentNumberCuts_(0),
3489  maximumDepth_(0),
3490  walkback_(NULL),
3491  addedCuts_(NULL),
3492  nextRowCut_(NULL),
3493  currentNode_(NULL),
3494  integerVariable_(NULL),
3495  integerInfo_(NULL),
3496  continuousSolution_(NULL),
3497  usedInSolution_(NULL),
3498  specialOptions_(0),
3499  subTreeModel_(NULL),
3500  numberStoppedSubTrees_(0),
3501  mutex_(NULL),
3502  presolve_(0),
3503  numberStrong_(5),
3504  numberBeforeTrust_(10),
3505  numberPenalties_(20),
3506  stopNumberIterations_(-1),
3507  penaltyScaleFactor_(3.0),
3508  numberAnalyzeIterations_(0),
3509  analyzeResults_(NULL),
3510  numberInfeasibleNodes_(0),
3511  problemType_(0),
3512  printFrequency_(0),
3513  numberCutGenerators_(0),
3514  generator_(NULL),
3515  virginGenerator_(NULL),
3516  numberHeuristics_(0),
3517  heuristic_(NULL),
3518  lastHeuristic_(NULL),
3519  eventHandler_(0),
3520  numberObjects_(0),
3521  object_(NULL),
3522  ownObjects_(true),
3523  originalColumns_(NULL),
3524  howOftenGlobalScan_(1),
3525  numberGlobalViolations_(0),
3526  continuousObjective_(COIN_DBL_MAX),
3527  originalContinuousObjective_(COIN_DBL_MAX),
3528  continuousInfeasibilities_(COIN_INT_MAX),
3529  maximumCutPassesAtRoot_(20),
3530  maximumCutPasses_(10),
3531  preferredWay_(0),
3532  currentPassNumber_(0),
3533  maximumWhich_(1000),
3534  maximumRows_(0),
3535  whichGenerator_(NULL),
3536  maximumStatistics_(0),
3537  statistics_(NULL),
3538  maximumDepthActual_(0),
3539  numberDJFixed_(0.0),
3540  probingInfo_(NULL),
3541  numberFixedAtRoot_(0),
3542  numberFixedNow_(0),
3543  stoppedOnGap_(false),
3544  eventHappened_(false),
3545  numberLongStrong_(0),
3546  numberOldActiveCuts_(0),
3547  numberNewCuts_(0),
3548  sizeMiniTree_(0),
3549  searchStrategy_(-1),
3550  numberStrongIterations_(0),
3551  resolveAfterTakeOffCuts_(true),
3552#if NEW_UPDATE_OBJECT>1
3553  numberUpdateItems_(0),
3554  maximumNumberUpdateItems_(0),
3555  updateItems_(NULL),
3556#endif
3557  numberThreads_(0),
3558  threadMode_(0)
3559{
3560  memset(intParam_,0,sizeof(intParam_));
3561  intParam_[CbcMaxNumNode] = 2147483647;
3562  intParam_[CbcMaxNumSol] = 9999999;
3563  intParam_[CbcFathomDiscipline] = 0;
3564
3565  dblParam_[CbcIntegerTolerance] = 1e-6;
3566  dblParam_[CbcInfeasibilityWeight] = 0.0;
3567  dblParam_[CbcCutoffIncrement] = 1e-5;
3568  dblParam_[CbcAllowableGap] = 1.0e-10;
3569  dblParam_[CbcAllowableFractionGap] = 0.0;
3570  dblParam_[CbcMaximumSeconds] = 1.0e100;
3571  dblParam_[CbcCurrentCutoff] = 1.0e100;
3572  dblParam_[CbcOptimizationDirection] = 1.0;
3573  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
3574  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
3575  dblParam_[CbcStartSeconds] = 0.0;
3576  strongInfo_[0]=0;
3577  strongInfo_[1]=0;
3578  strongInfo_[2]=0;
3579  solverCharacteristics_ = NULL;
3580  nodeCompare_=new CbcCompareDefault();;
3581  problemFeasibility_=new CbcFeasibilityBase();
3582  tree_= new CbcTree();
3583  branchingMethod_=NULL;
3584  cutModifier_=NULL;
3585  strategy_=NULL;
3586  parentModel_=NULL;
3587  cbcColLower_ = NULL;
3588  cbcColUpper_ = NULL;
3589  cbcRowLower_ = NULL;
3590  cbcRowUpper_ = NULL;
3591  cbcColSolution_ = NULL;
3592  cbcRowPrice_ = NULL;
3593  cbcReducedCost_ = NULL;
3594  cbcRowActivity_ = NULL;
3595  appData_=NULL;
3596  handler_ = new CoinMessageHandler();
3597  handler_->setLogLevel(2);
3598  messages_ = CbcMessage();
3599  eventHandler_ = new CbcEventHandler() ;
3600}
3601
3602/** Constructor from solver.
3603
3604  Creates a model complete with a clone of the solver passed as a parameter.
3605*/
3606
3607CbcModel::CbcModel(const OsiSolverInterface &rhs)
3608:
3609  continuousSolver_(NULL),
3610  referenceSolver_(NULL),
3611  defaultHandler_(true),
3612  emptyWarmStart_(NULL),
3613  bestObjective_(COIN_DBL_MAX),
3614  bestPossibleObjective_(COIN_DBL_MAX),
3615  sumChangeObjective1_(0.0),
3616  sumChangeObjective2_(0.0),
3617  minimumDrop_(1.0e-4),
3618  numberSolutions_(0),
3619  stateOfSearch_(0),
3620  hotstartSolution_(NULL),
3621  hotstartPriorities_(NULL),
3622  numberHeuristicSolutions_(0),
3623  numberNodes_(0),
3624  numberNodes2_(0),
3625  numberIterations_(0),
3626  status_(-1),
3627  secondaryStatus_(-1),
3628  numberRowsAtContinuous_(0),
3629  maximumNumberCuts_(0),
3630  phase_(0),
3631  currentNumberCuts_(0),
3632  maximumDepth_(0),
3633  walkback_(NULL),
3634  addedCuts_(NULL),
3635  nextRowCut_(NULL),
3636  currentNode_(NULL),
3637  integerInfo_(NULL),
3638  specialOptions_(0),
3639  subTreeModel_(NULL),
3640  numberStoppedSubTrees_(0),
3641  mutex_(NULL),
3642  presolve_(0),
3643  numberStrong_(5),
3644  numberBeforeTrust_(10),
3645  numberPenalties_(20),
3646  stopNumberIterations_(-1),
3647  penaltyScaleFactor_(3.0),
3648  numberAnalyzeIterations_(0),
3649  analyzeResults_(NULL),
3650  numberInfeasibleNodes_(0),
3651  problemType_(0),
3652  printFrequency_(0),
3653  numberCutGenerators_(0),
3654  generator_(NULL),
3655  virginGenerator_(NULL),
3656  numberHeuristics_(0),
3657  heuristic_(NULL),
3658  lastHeuristic_(NULL),
3659  eventHandler_(0),
3660  numberObjects_(0),
3661  object_(NULL),
3662  ownObjects_(true),
3663  originalColumns_(NULL),
3664  howOftenGlobalScan_(1),
3665  numberGlobalViolations_(0),
3666  continuousObjective_(COIN_DBL_MAX),
3667  originalContinuousObjective_(COIN_DBL_MAX),
3668  continuousInfeasibilities_(COIN_INT_MAX),
3669  maximumCutPassesAtRoot_(20),
3670  maximumCutPasses_(10),
3671  preferredWay_(0),
3672  currentPassNumber_(0),
3673  maximumWhich_(1000),
3674  maximumRows_(0),
3675  whichGenerator_(NULL),
3676  maximumStatistics_(0),
3677  statistics_(NULL),
3678  maximumDepthActual_(0),
3679  numberDJFixed_(0.0),
3680  probingInfo_(NULL),
3681  numberFixedAtRoot_(0),
3682  numberFixedNow_(0),
3683  stoppedOnGap_(false),
3684  eventHappened_(false),
3685  numberLongStrong_(0),
3686  numberOldActiveCuts_(0),
3687  numberNewCuts_(0),
3688  sizeMiniTree_(0),
3689  searchStrategy_(-1),
3690  numberStrongIterations_(0),
3691  resolveAfterTakeOffCuts_(true),
3692#if NEW_UPDATE_OBJECT>1
3693  numberUpdateItems_(0),
3694  maximumNumberUpdateItems_(0),
3695  updateItems_(NULL),
3696#endif
3697  numberThreads_(0),
3698  threadMode_(0)
3699{
3700  memset(intParam_,0,sizeof(intParam_));
3701  intParam_[CbcMaxNumNode] = 2147483647;
3702  intParam_[CbcMaxNumSol] = 9999999;
3703  intParam_[CbcFathomDiscipline] = 0;
3704
3705  dblParam_[CbcIntegerTolerance] = 1e-6;
3706  dblParam_[CbcInfeasibilityWeight] = 0.0;
3707  dblParam_[CbcCutoffIncrement] = 1e-5;
3708  dblParam_[CbcAllowableGap] = 1.0e-10;
3709  dblParam_[CbcAllowableFractionGap] = 0.0;
3710  dblParam_[CbcMaximumSeconds] = 1.0e100;
3711  dblParam_[CbcCurrentCutoff] = 1.0e100;
3712  dblParam_[CbcOptimizationDirection] = 1.0;
3713  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
3714  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
3715  dblParam_[CbcStartSeconds] = 0.0;
3716  strongInfo_[0]=0;
3717  strongInfo_[1]=0;
3718  strongInfo_[2]=0;
3719  solverCharacteristics_ = NULL;
3720
3721  nodeCompare_=new CbcCompareDefault();;
3722  problemFeasibility_=new CbcFeasibilityBase();
3723  tree_= new CbcTree();
3724  branchingMethod_=NULL;
3725  cutModifier_=NULL;
3726  strategy_=NULL;
3727  parentModel_=NULL;
3728  appData_=NULL;
3729  handler_ = new CoinMessageHandler();
3730  handler_->setLogLevel(2);
3731  messages_ = CbcMessage();
3732  eventHandler_ = new CbcEventHandler() ;
3733  solver_ = rhs.clone();
3734  referenceSolver_ = solver_->clone();
3735  ownership_ = 0x80000000;
3736  cbcColLower_ = NULL;
3737  cbcColUpper_ = NULL;
3738  cbcRowLower_ = NULL;
3739  cbcRowUpper_ = NULL;
3740  cbcColSolution_ = NULL;
3741  cbcRowPrice_ = NULL;
3742  cbcReducedCost_ = NULL;
3743  cbcRowActivity_ = NULL;
3744
3745  // Initialize solution and integer variable vectors
3746  bestSolution_ = NULL; // to say no solution found
3747  numberIntegers_=0;
3748  int numberColumns = solver_->getNumCols();
3749  int iColumn;
3750  if (numberColumns) {
3751    // Space for current solution
3752    currentSolution_ = new double[numberColumns];
3753    continuousSolution_ = new double[numberColumns];
3754    usedInSolution_ = new int[numberColumns];
3755    CoinZeroN(usedInSolution_,numberColumns);
3756    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3757      if( solver_->isInteger(iColumn)) 
3758        numberIntegers_++;
3759    }
3760  } else {
3761    // empty model
3762    currentSolution_=NULL;
3763    continuousSolution_=NULL;
3764    usedInSolution_=NULL;
3765  }
3766  testSolution_=currentSolution_;
3767  if (numberIntegers_) {
3768    integerVariable_ = new int [numberIntegers_];
3769    numberIntegers_=0;
3770    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3771      if( solver_->isInteger(iColumn)) 
3772        integerVariable_[numberIntegers_++]=iColumn;
3773    }
3774  } else {
3775    integerVariable_ = NULL;
3776  }
3777}
3778
3779/*
3780  Assign a solver to the model (model assumes ownership)
3781
3782  The integer variable vector is initialized if it's not already present.
3783  If deleteSolver then current solver deleted (if model owned)
3784
3785  Assuming ownership matches usage in OsiSolverInterface
3786  (cf. assignProblem, loadProblem).
3787
3788  TODO: What to do about solver parameters? A simple copy likely won't do it,
3789        because the SI must push the settings into the underlying solver. In
3790        the context of switching solvers in cbc, this means that command line
3791        settings will get lost. Stash the command line somewhere and reread it
3792        here, maybe?
3793 
3794  TODO: More generally, how much state should be transferred from the old
3795        solver to the new solver? Best perhaps to see how usage develops.
3796        What's done here mimics the CbcModel(OsiSolverInterface) constructor.
3797*/
3798void
3799CbcModel::assignSolver(OsiSolverInterface *&solver, bool deleteSolver)
3800
3801{
3802  // resize best solution if exists
3803  if (bestSolution_&&solver&&solver_) {
3804    int nOld = solver_->getNumCols();
3805    int nNew = solver->getNumCols();
3806    if (nNew>nOld) {
3807      double * temp = new double[nNew];
3808      memcpy(temp,bestSolution_,nOld*sizeof(double));
3809      memset(temp+nOld,0,(nNew-nOld)*sizeof(double));
3810      delete [] bestSolution_;
3811      bestSolution_=temp;
3812    }
3813  }
3814  // Keep the current message level for solver (if solver exists)
3815  if (solver_)
3816    solver->messageHandler()->setLogLevel(solver_->messageHandler()->logLevel()) ;
3817
3818  if (modelOwnsSolver()&&deleteSolver) delete solver_ ;
3819  solver_ = solver;
3820  solver = NULL ;
3821  setModelOwnsSolver(true) ;
3822/*
3823  Basis information is solver-specific.
3824*/
3825  if (emptyWarmStart_)
3826  { delete emptyWarmStart_  ;
3827    emptyWarmStart_ = 0 ; }
3828  bestSolutionBasis_ = CoinWarmStartBasis();
3829/*
3830  Initialize integer variable vector.
3831*/
3832  numberIntegers_=0;
3833  int numberColumns = solver_->getNumCols();
3834  int iColumn;
3835  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3836    if( solver_->isInteger(iColumn)) 
3837      numberIntegers_++;
3838  }
3839  delete [] integerVariable_;
3840  if (numberIntegers_) {
3841    integerVariable_ = new int [numberIntegers_];
3842    numberIntegers_=0;
3843    for (iColumn=0;iColumn<numberColumns;iColumn++) {
3844      if( solver_->isInteger(iColumn)) 
3845        integerVariable_[numberIntegers_++]=iColumn;
3846    }
3847  } else {
3848    integerVariable_ = NULL;
3849  }
3850
3851  return ;
3852}
3853
3854// Copy constructor.
3855
3856CbcModel::CbcModel(const CbcModel & rhs, bool cloneHandler)
3857:
3858  continuousSolver_(NULL),
3859  referenceSolver_(NULL),
3860  defaultHandler_(rhs.defaultHandler_),
3861  emptyWarmStart_(NULL),
3862  bestObjective_(rhs.bestObjective_),
3863  bestPossibleObjective_(rhs.bestPossibleObjective_),
3864  sumChangeObjective1_(rhs.sumChangeObjective1_),
3865  sumChangeObjective2_(rhs.sumChangeObjective2_),
3866  minimumDrop_(rhs.minimumDrop_),
3867  numberSolutions_(rhs.numberSolutions_),
3868  stateOfSearch_(rhs.stateOfSearch_),
3869  numberHeuristicSolutions_(rhs.numberHeuristicSolutions_),
3870  numberNodes_(rhs.numberNodes_),
3871  numberNodes2_(rhs.numberNodes2_),
3872  numberIterations_(rhs.numberIterations_),
3873  status_(rhs.status_),
3874  secondaryStatus_(rhs.secondaryStatus_),
3875  specialOptions_(rhs.specialOptions_),
3876  subTreeModel_(rhs.subTreeModel_),
3877  numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
3878  mutex_(NULL),
3879  presolve_(rhs.presolve_),
3880  numberStrong_(rhs.numberStrong_),
3881  numberBeforeTrust_(rhs.numberBeforeTrust_),
3882  numberPenalties_(rhs.numberPenalties_),
3883  stopNumberIterations_(rhs.stopNumberIterations_),
3884  penaltyScaleFactor_(rhs.penaltyScaleFactor_),
3885  numberAnalyzeIterations_(rhs.numberAnalyzeIterations_),
3886  analyzeResults_(NULL),
3887  numberInfeasibleNodes_(rhs.numberInfeasibleNodes_),
3888  problemType_(rhs.problemType_),
3889  printFrequency_(rhs.printFrequency_),
3890  howOftenGlobalScan_(rhs.howOftenGlobalScan_),
3891  numberGlobalViolations_(rhs.numberGlobalViolations_),
3892  continuousObjective_(rhs.continuousObjective_),
3893  originalContinuousObjective_(rhs.originalContinuousObjective_),
3894  continuousInfeasibilities_(rhs.continuousInfeasibilities_),
3895  maximumCutPassesAtRoot_(rhs.maximumCutPassesAtRoot_),
3896  maximumCutPasses_( rhs.maximumCutPasses_),
3897  preferredWay_(rhs.preferredWay_),
3898  currentPassNumber_(rhs.currentPassNumber_),
3899  maximumWhich_(rhs.maximumWhich_),
3900  maximumRows_(0),
3901  whichGenerator_(NULL),
3902  maximumStatistics_(0),
3903  statistics_(NULL),
3904  maximumDepthActual_(0),
3905  numberDJFixed_(0.0),
3906  probingInfo_(NULL),
3907  numberFixedAtRoot_(rhs.numberFixedAtRoot_),
3908  numberFixedNow_(rhs.numberFixedNow_),
3909  stoppedOnGap_(rhs.stoppedOnGap_),
3910  eventHappened_(rhs.eventHappened_),
3911  numberLongStrong_(rhs.numberLongStrong_),
3912  numberOldActiveCuts_(rhs.numberOldActiveCuts_),
3913  numberNewCuts_(rhs.numberNewCuts_),
3914  sizeMiniTree_(rhs.sizeMiniTree_),
3915  searchStrategy_(rhs.searchStrategy_),
3916  numberStrongIterations_(rhs.numberStrongIterations_),
3917  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
3918#if NEW_UPDATE_OBJECT>1
3919  numberUpdateItems_(rhs.numberUpdateItems_),
3920  maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_),
3921  updateItems_(NULL),
3922#endif
3923  numberThreads_(rhs.numberThreads_),
3924  threadMode_(rhs.threadMode_)
3925{
3926  memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
3927  memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
3928  strongInfo_[0]=rhs.strongInfo_[0];
3929  strongInfo_[1]=rhs.strongInfo_[1];
3930  strongInfo_[2]=rhs.strongInfo_[2];
3931  solverCharacteristics_ = NULL;
3932  if (rhs.emptyWarmStart_) emptyWarmStart_ = rhs.emptyWarmStart_->clone() ;
3933  bool noTree=false;
3934  if (defaultHandler_||cloneHandler) {
3935    handler_ = new CoinMessageHandler();
3936    handler_->setLogLevel(2);
3937  } else {
3938    handler_ = rhs.handler_;
3939  }
3940  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
3941  numberCutGenerators_ = rhs.numberCutGenerators_;
3942  if (numberCutGenerators_) {
3943    generator_ = new CbcCutGenerator * [numberCutGenerators_];
3944    virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
3945    int i;
3946    for (i=0;i<numberCutGenerators_;i++) {
3947      generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
3948      virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
3949    }
3950  } else {
3951    generator_=NULL;
3952    virginGenerator_=NULL;
3953  }
3954  if (!noTree)
3955    globalCuts_ = rhs.globalCuts_;
3956  numberHeuristics_ = rhs.numberHeuristics_;
3957  if (numberHeuristics_) {
3958    heuristic_ = new CbcHeuristic * [numberHeuristics_];
3959    int i;
3960    for (i=0;i<numberHeuristics_;i++) {
3961      heuristic_[i]=rhs.heuristic_[i]->clone();
3962    }
3963  } else {
3964    heuristic_=NULL;
3965  }
3966  lastHeuristic_ = NULL;
3967  if (rhs.eventHandler_)
3968    { eventHandler_ = rhs.eventHandler_->clone() ; }
3969  else
3970  { eventHandler_ = NULL ; }
3971  ownObjects_ = rhs.ownObjects_;
3972  if (ownObjects_) {
3973    numberObjects_=rhs.numberObjects_;
3974    if (numberObjects_) {
3975      object_ = new OsiObject * [numberObjects_];
3976      int i;
3977      for (i=0;i<numberObjects_;i++) {
3978        object_[i]=(rhs.object_[i])->clone();
3979        CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
3980        // Could be OsiObjects
3981        if (obj)
3982          obj->setModel(this);
3983      }
3984    } else {
3985      object_=NULL;
3986    }
3987  } else {
3988    // assume will be redone
3989    numberObjects_=0;
3990    object_=NULL;
3991  }
3992  if (rhs.referenceSolver_)
3993    referenceSolver_ = rhs.referenceSolver_->clone();
3994  else
3995    referenceSolver_=NULL;
3996  if (!noTree||!rhs.continuousSolver_)
3997    solver_ = rhs.solver_->clone();
3998  else
3999    solver_ = rhs.continuousSolver_->clone();
4000  if (rhs.originalColumns_) {
4001    int numberColumns = solver_->getNumCols();
4002    originalColumns_= new int [numberColumns];
4003    memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
4004  } else {
4005    originalColumns_=NULL;
4006  }
4007#if NEW_UPDATE_OBJECT>1
4008  if (maximumNumberUpdateItems_) {
4009    updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
4010    for (int i=0;i<maximumNumberUpdateItems_;i++)
4011      updateItems_[i] = rhs.updateItems_[i];
4012  }
4013#endif
4014  if (maximumWhich_&&rhs.whichGenerator_)
4015    whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
4016  nodeCompare_=rhs.nodeCompare_->clone();
4017  problemFeasibility_=rhs.problemFeasibility_->clone();
4018  tree_= rhs.tree_->clone();
4019  if (rhs.branchingMethod_)
4020    branchingMethod_=rhs.branchingMethod_->clone();
4021  else
4022    branchingMethod_=NULL;
4023  if (rhs.cutModifier_)
4024    cutModifier_=rhs.cutModifier_->clone();
4025  else
4026    cutModifier_=NULL;
4027  cbcColLower_ = NULL;
4028  cbcColUpper_ = NULL;
4029  cbcRowLower_ = NULL;
4030  cbcRowUpper_ = NULL;
4031  cbcColSolution_ = NULL;
4032  cbcRowPrice_ = NULL;
4033  cbcReducedCost_ = NULL;
4034  cbcRowActivity_ = NULL;
4035  if (rhs.strategy_)
4036    strategy_=rhs.strategy_->clone();
4037  else
4038    strategy_=NULL;
4039  parentModel_=rhs.parentModel_;
4040  appData_=rhs.appData_;
4041  messages_ = rhs.messages_;
4042  ownership_ = 0x80000000;
4043  messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4044  numberIntegers_=rhs.numberIntegers_;
4045  if (numberIntegers_) {
4046    integerVariable_ = new int [numberIntegers_];
4047    memcpy(integerVariable_,rhs.integerVariable_,numberIntegers_*sizeof(int));
4048    integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,solver_->getNumCols());
4049  } else {
4050    integerVariable_ = NULL;
4051    integerInfo_=NULL;
4052  }
4053  if (rhs.hotstartSolution_) {
4054    int numberColumns = solver_->getNumCols();
4055    hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
4056    hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
4057  } else {
4058    hotstartSolution_ = NULL;
4059    hotstartPriorities_ =NULL;
4060  }
4061  if (rhs.bestSolution_&&!noTree) {
4062    int numberColumns = solver_->getNumCols();
4063    bestSolution_ = new double[numberColumns];
4064    memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
4065  } else {
4066    bestSolution_=NULL;
4067  }
4068  if (!noTree) {
4069    int numberColumns = solver_->getNumCols();
4070    // Space for current solution
4071    currentSolution_ = new double[numberColumns];
4072    continuousSolution_ = new double[numberColumns];
4073    usedInSolution_ = new int[numberColumns];
4074    CoinZeroN(usedInSolution_,numberColumns);
4075  } else {
4076    currentSolution_=NULL;
4077    continuousSolution_=NULL;
4078    usedInSolution_=NULL;
4079  }
4080  testSolution_=currentSolution_;
4081  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
4082  maximumNumberCuts_=rhs.maximumNumberCuts_;
4083  phase_ = rhs.phase_;
4084  currentNumberCuts_=rhs.currentNumberCuts_;
4085  maximumDepth_= rhs.maximumDepth_;
4086  if (noTree) {
4087    bestObjective_ = COIN_DBL_MAX;
4088    numberSolutions_ =0;
4089    stateOfSearch_= 0;
4090    numberHeuristicSolutions_=0;
4091    numberNodes_=0;
4092    numberNodes2_=0;
4093    numberIterations_=0;
4094    status_=0;
4095    subTreeModel_=NULL;
4096    numberStoppedSubTrees_=0;
4097    continuousObjective_=COIN_DBL_MAX;
4098    originalContinuousObjective_=COIN_DBL_MAX;
4099    continuousInfeasibilities_=COIN_INT_MAX;
4100    maximumNumberCuts_=0;
4101    tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_);
4102    bestPossibleObjective_ = COIN_DBL_MAX;
4103  }
4104  // These are only used as temporary arrays so need not be filled
4105  if (maximumNumberCuts_) {
4106    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4107  } else {
4108    addedCuts_ = NULL;
4109  }
4110  bestSolutionBasis_ = rhs.bestSolutionBasis_;
4111  nextRowCut_ = NULL;
4112  currentNode_ = NULL;
4113  if (maximumDepth_)
4114    walkback_ = new CbcNodeInfo * [maximumDepth_];
4115  else
4116    walkback_ = NULL;
4117  synchronizeModel();
4118  if (cloneHandler&&!defaultHandler_) {
4119    delete handler_;
4120    CoinMessageHandler * handler = rhs.handler_->clone();
4121    passInMessageHandler(handler);
4122  }
4123}
4124 
4125// Assignment operator
4126CbcModel & 
4127CbcModel::operator=(const CbcModel& rhs)
4128{
4129  if (this!=&rhs) {
4130    if (modelOwnsSolver()) {
4131      delete solver_;
4132      solver_=NULL;
4133    }
4134    gutsOfDestructor();
4135    if (defaultHandler_)
4136    { delete handler_;
4137      handler_ = NULL; }
4138    defaultHandler_ = rhs.defaultHandler_;
4139    if (defaultHandler_)
4140    { handler_ = new CoinMessageHandler();
4141      handler_->setLogLevel(2); }
4142    else
4143    { handler_ = rhs.handler_; }
4144    messages_ = rhs.messages_;
4145    messageHandler()->setLogLevel(rhs.messageHandler()->logLevel());
4146    if (rhs.solver_)
4147    { solver_ = rhs.solver_->clone() ; }
4148    else
4149    { solver_ = 0 ; }
4150    ownership_ = 0x80000000;
4151    delete continuousSolver_ ;
4152    if (rhs.continuousSolver_)
4153    { continuousSolver_ = rhs.continuousSolver_->clone() ; }
4154    else
4155    { continuousSolver_ = 0 ; }
4156    delete referenceSolver_;
4157    if (rhs.referenceSolver_)
4158    { referenceSolver_ = rhs.referenceSolver_->clone() ; }
4159    else
4160    { referenceSolver_ = NULL ; }
4161
4162    delete emptyWarmStart_ ;
4163    if (rhs.emptyWarmStart_)
4164    { emptyWarmStart_ = rhs.emptyWarmStart_->clone() ; }
4165    else
4166    { emptyWarmStart_ = 0 ; }
4167
4168    bestObjective_ = rhs.bestObjective_;
4169    bestPossibleObjective_=rhs.bestPossibleObjective_;
4170    sumChangeObjective1_=rhs.sumChangeObjective1_;
4171    sumChangeObjective2_=rhs.sumChangeObjective2_;
4172    delete [] bestSolution_;
4173    if (rhs.bestSolution_) {
4174      int numberColumns = rhs.getNumCols();
4175      bestSolution_ = new double[numberColumns];
4176      memcpy(bestSolution_,rhs.bestSolution_,numberColumns*sizeof(double));
4177    } else {
4178      bestSolution_=NULL;
4179    }
4180    int numberColumns = rhs.getNumCols();
4181    if (numberColumns) {
4182      // Space for current solution
4183      currentSolution_ = new double[numberColumns];
4184      continuousSolution_ = new double[numberColumns];
4185      usedInSolution_ = new int[numberColumns];
4186      CoinZeroN(usedInSolution_,numberColumns);
4187    } else {
4188      currentSolution_=NULL;
4189      continuousSolution_=NULL;
4190      usedInSolution_=NULL;
4191    }
4192    testSolution_=currentSolution_;
4193    minimumDrop_ = rhs.minimumDrop_;
4194    numberSolutions_=rhs.numberSolutions_;
4195    stateOfSearch_= rhs.stateOfSearch_;
4196    numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
4197    numberNodes_ = rhs.numberNodes_;
4198    numberNodes2_ = rhs.numberNodes2_;
4199    numberIterations_ = rhs.numberIterations_;
4200    status_ = rhs.status_;
4201    secondaryStatus_ = rhs.secondaryStatus_;
4202    specialOptions_ = rhs.specialOptions_;
4203    subTreeModel_ = rhs.subTreeModel_;
4204    numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
4205    mutex_ = NULL;
4206    presolve_ = rhs.presolve_;
4207    numberStrong_ = rhs.numberStrong_;
4208    numberBeforeTrust_ = rhs.numberBeforeTrust_;
4209    numberPenalties_ = rhs.numberPenalties_;
4210    stopNumberIterations_ = rhs.stopNumberIterations_;
4211    penaltyScaleFactor_ = rhs.penaltyScaleFactor_;
4212    numberAnalyzeIterations_ = rhs.numberAnalyzeIterations_;
4213    delete [] analyzeResults_;
4214    analyzeResults_ = NULL;
4215    numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
4216    problemType_ = rhs.problemType_;
4217    printFrequency_ = rhs.printFrequency_;
4218    howOftenGlobalScan_=rhs.howOftenGlobalScan_;
4219    numberGlobalViolations_=rhs.numberGlobalViolations_;
4220    continuousObjective_=rhs.continuousObjective_;
4221    originalContinuousObjective_ = rhs.originalContinuousObjective_;
4222    continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
4223    maximumCutPassesAtRoot_ = rhs.maximumCutPassesAtRoot_;
4224    maximumCutPasses_ = rhs.maximumCutPasses_;
4225    preferredWay_ = rhs.preferredWay_;
4226    currentPassNumber_ = rhs.currentPassNumber_;
4227    memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
4228    memcpy(dblParam_,rhs.dblParam_,sizeof(dblParam_));
4229    globalCuts_ = rhs.globalCuts_;
4230    int i;
4231    for (i=0;i<numberCutGenerators_;i++) {
4232      delete generator_[i];
4233      delete virginGenerator_[i];
4234    }
4235    delete [] generator_;
4236    delete [] virginGenerator_;
4237    delete [] heuristic_;
4238    maximumWhich_ = rhs.maximumWhich_;
4239    delete [] whichGenerator_;
4240    whichGenerator_ = NULL;
4241    if (maximumWhich_&&rhs.whichGenerator_)
4242      whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
4243    maximumRows_=0;
4244    workingBasis_ = CoinWarmStartBasis();
4245    for (i=0;i<maximumStatistics_;i++)
4246      delete statistics_[i];
4247    delete [] statistics_;
4248    maximumStatistics_ = 0;
4249    statistics_ = NULL;
4250    delete probingInfo_;
4251    probingInfo_=NULL;
4252    numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
4253    numberFixedNow_ = rhs.numberFixedNow_;
4254    stoppedOnGap_ = rhs.stoppedOnGap_;
4255    eventHappened_ = rhs.eventHappened_;
4256    numberLongStrong_ = rhs.numberLongStrong_;
4257    numberOldActiveCuts_ = rhs.numberOldActiveCuts_;
4258    numberNewCuts_ = rhs.numberNewCuts_;
4259    resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_;
4260#if NEW_UPDATE_OBJECT>1
4261    numberUpdateItems_ = rhs.numberUpdateItems_;
4262    maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_;
4263    delete [] updateItems_;
4264    if (maximumNumberUpdateItems_) {
4265      updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
4266      for (i=0;i<maximumNumberUpdateItems_;i++)
4267        updateItems_[i] = rhs.updateItems_[i];
4268    } else {
4269      updateItems_ = NULL;
4270    }
4271#endif
4272    numberThreads_ = rhs.numberThreads_;
4273    threadMode_ = rhs.threadMode_;
4274    sizeMiniTree_ = rhs.sizeMiniTree_;
4275    searchStrategy_ = rhs.searchStrategy_;
4276    numberStrongIterations_ = rhs.numberStrongIterations_;
4277    strongInfo_[0]=rhs.strongInfo_[0];
4278    strongInfo_[1]=rhs.strongInfo_[1];
4279    strongInfo_[2]=rhs.strongInfo_[2];
4280    solverCharacteristics_ = NULL;
4281    lastHeuristic_ = NULL;
4282    numberCutGenerators_ = rhs.numberCutGenerators_;
4283    if (numberCutGenerators_) {
4284      generator_ = new CbcCutGenerator * [numberCutGenerators_];
4285      virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_];
4286      int i;
4287      for (i=0;i<numberCutGenerators_;i++) {
4288        generator_[i]=new CbcCutGenerator(*rhs.generator_[i]);
4289        virginGenerator_[i]=new CbcCutGenerator(*rhs.virginGenerator_[i]);
4290      }
4291    } else {
4292      generator_=NULL;
4293      virginGenerator_=NULL;
4294    }
4295    numberHeuristics_ = rhs.numberHeuristics_;
4296    if (numberHeuristics_) {
4297      heuristic_ = new CbcHeuristic * [numberHeuristics_];
4298      memcpy(heuristic_,rhs.heuristic_,
4299             numberHeuristics_*sizeof(CbcHeuristic *));
4300    } else {
4301      heuristic_=NULL;
4302    }
4303    lastHeuristic_ = NULL;
4304    if (eventHandler_)
4305      delete eventHandler_ ;
4306    if (rhs.eventHandler_)
4307      { eventHandler_ = rhs.eventHandler_->clone() ; }
4308    else
4309    { eventHandler_ = NULL ; }
4310    if (ownObjects_) {
4311      for (i=0;i<numberObjects_;i++)
4312        delete object_[i];
4313      delete [] object_;
4314      numberObjects_=rhs.numberObjects_;
4315      if (numberObjects_) {
4316        object_ = new OsiObject * [numberObjects_];
4317        int i;
4318        for (i=0;i<numberObjects_;i++) 
4319          object_[i]=(rhs.object_[i])->clone();
4320      } else {
4321        object_=NULL;
4322    }
4323    } else {
4324      // assume will be redone
4325      numberObjects_=0;
4326      object_=NULL;
4327    }
4328    delete [] originalColumns_;
4329    if (rhs.originalColumns_) {
4330      int numberColumns = rhs.getNumCols();
4331      originalColumns_= new int [numberColumns];
4332      memcpy(originalColumns_,rhs.originalColumns_,numberColumns*sizeof(int));
4333    } else {
4334      originalColumns_=NULL;
4335    }
4336    nodeCompare_=rhs.nodeCompare_->clone();
4337    problemFeasibility_=rhs.problemFeasibility_->clone();
4338    delete tree_;
4339    tree_= rhs.tree_->clone();
4340    if (rhs.branchingMethod_)
4341      branchingMethod_=rhs.branchingMethod_->clone();
4342    else
4343      branchingMethod_=NULL;
4344    if (rhs.cutModifier_)
4345      cutModifier_=rhs.cutModifier_->clone();
4346    else
4347      cutModifier_=NULL;
4348    delete strategy_;
4349    if (rhs.strategy_)
4350      strategy_=rhs.strategy_->clone();
4351    else
4352      strategy_=NULL;
4353    parentModel_=rhs.parentModel_;
4354    appData_=rhs.appData_;
4355
4356    delete [] integerVariable_;
4357    numberIntegers_=rhs.numberIntegers_;
4358    if (numberIntegers_) {
4359      integerVariable_ = new int [numberIntegers_];
4360      memcpy(integerVariable_,rhs.integerVariable_,
4361             numberIntegers_*sizeof(int));
4362      integerInfo_ = CoinCopyOfArray(rhs.integerInfo_,rhs.getNumCols());
4363    } else {
4364      integerVariable_ = NULL;
4365      integerInfo_=NULL;
4366    }
4367    if (rhs.hotstartSolution_) {
4368      int numberColumns = solver_->getNumCols();
4369      hotstartSolution_ = CoinCopyOfArray(rhs.hotstartSolution_,numberColumns);
4370      hotstartPriorities_ = CoinCopyOfArray(rhs.hotstartPriorities_,numberColumns);
4371    } else {
4372      hotstartSolution_ = NULL;
4373      hotstartPriorities_ =NULL;
4374    }
4375    numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
4376    maximumNumberCuts_=rhs.maximumNumberCuts_;
4377    phase_ = rhs.phase_;
4378    currentNumberCuts_=rhs.currentNumberCuts_;
4379    maximumDepth_= rhs.maximumDepth_;
4380    delete [] addedCuts_;
4381    delete [] walkback_;
4382    // These are only used as temporary arrays so need not be filled
4383    if (maximumNumberCuts_) {
4384      addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4385    } else {
4386      addedCuts_ = NULL;
4387    }
4388    bestSolutionBasis_ = rhs.bestSolutionBasis_;
4389    nextRowCut_ = NULL;
4390    currentNode_ = NULL;
4391    if (maximumDepth_)
4392      walkback_ = new CbcNodeInfo * [maximumDepth_];
4393    else
4394      walkback_ = NULL;
4395    synchronizeModel();
4396    cbcColLower_ = NULL;
4397    cbcColUpper_ = NULL;
4398    cbcRowLower_ = NULL;
4399    cbcRowUpper_ = NULL;
4400    cbcColSolution_ = NULL;
4401    cbcRowPrice_ = NULL;
4402    cbcReducedCost_ = NULL;
4403    cbcRowActivity_ = NULL;
4404  }
4405  return *this;
4406}
4407// Destructor
4408CbcModel::~CbcModel ()
4409{
4410  if (defaultHandler_) {
4411    delete handler_;
4412    handler_ = NULL;
4413  }
4414  delete tree_;
4415  tree_=NULL;
4416  if (modelOwnsSolver()) {
4417    delete solver_;
4418    solver_ = NULL;
4419  }
4420  gutsOfDestructor();
4421  delete eventHandler_ ;
4422  eventHandler_ = NULL ;
4423}
4424// Clears out as much as possible (except solver)
4425void 
4426CbcModel::gutsOfDestructor()
4427{
4428  delete referenceSolver_;
4429  referenceSolver_=NULL;
4430  int i;
4431  for (i=0;i<numberCutGenerators_;i++) {
4432    delete generator_[i];
4433    delete virginGenerator_[i];
4434  }
4435  delete [] generator_;
4436  delete [] virginGenerator_;
4437  generator_=NULL;
4438  virginGenerator_=NULL;
4439  for (i=0;i<numberHeuristics_;i++)
4440    delete heuristic_[i];
4441  delete [] heuristic_;
4442  heuristic_=NULL;
4443  delete nodeCompare_;
4444  nodeCompare_=NULL;
4445  delete problemFeasibility_;
4446  problemFeasibility_=NULL;
4447  delete [] originalColumns_;
4448  originalColumns_=NULL;
4449  delete strategy_;
4450#if NEW_UPDATE_OBJECT>1
4451  delete [] updateItems_;
4452  updateItems_=NULL;
4453  numberUpdateItems_=0;
4454  maximumNumberUpdateItems_=0;
4455#endif
4456  gutsOfDestructor2();
4457}
4458// Clears out enough to reset CbcModel
4459void 
4460CbcModel::gutsOfDestructor2()
4461{
4462  delete [] integerInfo_;
4463  integerInfo_=NULL;
4464  delete [] integerVariable_;
4465  integerVariable_=NULL;
4466  int i;
4467  if (ownObjects_) {
4468    for (i=0;i<numberObjects_;i++)
4469      delete object_[i];
4470    delete [] object_;
4471  }
4472  ownObjects_=true;
4473  object_=NULL;
4474  numberIntegers_=0;
4475  numberObjects_=0;
4476  // Below here is whatever consensus is
4477  ownership_ = 0x80000000;
4478  delete branchingMethod_;
4479  branchingMethod_=NULL;
4480  delete cutModifier_;
4481  cutModifier_=NULL;
4482  resetModel();
4483}
4484// Clears out enough to reset CbcModel
4485void 
4486CbcModel::resetModel()
4487{
4488  delete emptyWarmStart_ ;
4489  emptyWarmStart_ =NULL;
4490  delete continuousSolver_;
4491  continuousSolver_=NULL;
4492  delete [] bestSolution_;
4493  bestSolution_=NULL;
4494  delete [] currentSolution_;
4495  currentSolution_=NULL;
4496  delete [] continuousSolution_;
4497  continuousSolution_=NULL;
4498  solverCharacteristics_=NULL;
4499  delete [] usedInSolution_;
4500  usedInSolution_ = NULL;
4501  testSolution_=NULL;
4502  lastHeuristic_ = NULL;
4503  delete [] addedCuts_;
4504  addedCuts_=NULL;
4505  nextRowCut_ = NULL;
4506  currentNode_ = NULL;
4507  delete [] walkback_;
4508  walkback_=NULL;
4509  delete [] whichGenerator_;
4510  whichGenerator_ = NULL;
4511  for (int i=0;i<maximumStatistics_;i++)
4512    delete statistics_[i];
4513  delete [] statistics_;
4514  statistics_=NULL;
4515  maximumDepthActual_ = 0;
4516  numberDJFixed_ =0.0;
4517  delete probingInfo_;
4518  probingInfo_ = NULL;
4519  maximumStatistics_=0;
4520  delete [] analyzeResults_;
4521  analyzeResults_=NULL;
4522  bestObjective_=COIN_DBL_MAX;
4523  bestPossibleObjective_=COIN_DBL_MAX;
4524  sumChangeObjective1_=0.0;
4525  sumChangeObjective2_=0.0;
4526  numberSolutions_=0;
4527  stateOfSearch_=0;
4528  delete [] hotstartSolution_;
4529  hotstartSolution_=NULL;
4530  delete [] hotstartPriorities_;
4531  hotstartPriorities_=NULL;
4532  numberHeuristicSolutions_=0;
4533  numberNodes_=0;
4534  numberNodes2_=0;
4535  numberIterations_=0;
4536  status_=-1;
4537  secondaryStatus_=-1;
4538  maximumNumberCuts_=0;
4539  phase_=0;
4540  currentNumberCuts_=0;
4541  maximumDepth_=0;
4542  nextRowCut_=NULL;
4543  currentNode_=NULL;
4544  // clear out tree
4545  if (tree_&&tree_->size())
4546    tree_->cleanTree(this, -1.0e100,bestPossibleObjective_) ;
4547  subTreeModel_=NULL;
4548  numberStoppedSubTrees_=0;
4549  numberInfeasibleNodes_=0;
4550  numberGlobalViolations_=0;
4551  continuousObjective_=0.0;
4552  originalContinuousObjective_=0.0;
4553  continuousInfeasibilities_=0;
4554  numberFixedAtRoot_=0;
4555  numberFixedNow_=0;
4556  stoppedOnGap_=false;
4557  eventHappened_=false;
4558  numberLongStrong_=0;
4559  numberOldActiveCuts_=0;
4560  numberNewCuts_=0;
4561  searchStrategy_=-1;
4562  numberStrongIterations_=0;
4563  // Parameters which need to be reset
4564  setCutoff(COIN_DBL_MAX);
4565  dblParam_[CbcCutoffIncrement] = 1e-5;
4566  dblParam_[CbcCurrentCutoff] = 1.0e100;
4567  dblParam_[CbcCurrentObjectiveValue] = 1.0e100;
4568  dblParam_[CbcCurrentMinimizationObjectiveValue] = 1.0e100;
4569}
4570// Move status, nodes etc etc across
4571void 
4572CbcModel::moveInfo(const CbcModel & rhs)
4573{
4574  bestObjective_ = rhs.bestObjective_;
4575  bestPossibleObjective_=rhs.bestPossibleObjective_;
4576  numberSolutions_=rhs.numberSolutions_;
4577  numberHeuristicSolutions_=rhs.numberHeuristicSolutions_;
4578  numberNodes_ = rhs.numberNodes_;
4579  // In case of interrupt
4580  intParam_[0] = rhs.intParam_[0];
4581  numberNodes2_ = rhs.numberNodes2_;
4582  numberIterations_ = rhs.numberIterations_;
4583  status_ = rhs.status_;
4584  secondaryStatus_ = rhs.secondaryStatus_;
4585  numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
4586  numberInfeasibleNodes_ = rhs.numberInfeasibleNodes_;
4587  continuousObjective_=rhs.continuousObjective_;
4588  originalContinuousObjective_ = rhs.originalContinuousObjective_;
4589  continuousInfeasibilities_ = rhs.continuousInfeasibilities_;
4590  numberFixedAtRoot_ = rhs.numberFixedAtRoot_;
4591  numberFixedNow_ = rhs.numberFixedNow_;
4592  stoppedOnGap_ = rhs.stoppedOnGap_;
4593  eventHappened_ = rhs.eventHappened_;
4594  numberLongStrong_ = rhs.numberLongStrong_;
4595  numberStrongIterations_ = rhs.numberStrongIterations_;
4596  strongInfo_[0]=rhs.strongInfo_[0];
4597  strongInfo_[1]=rhs.strongInfo_[1];
4598  strongInfo_[2]=rhs.strongInfo_[2];
4599  numberRowsAtContinuous_ = rhs.numberRowsAtContinuous_;
4600  maximumDepth_= rhs.maximumDepth_;
4601}
4602// Save a copy of the current solver so can be reset to
4603void 
4604CbcModel::saveReferenceSolver()
4605{
4606  delete referenceSolver_;
4607  referenceSolver_= solver_->clone();
4608}
4609
4610// Uses a copy of reference solver to be current solver
4611void 
4612CbcModel::resetToReferenceSolver()
4613{
4614  delete solver_;
4615  solver_ = referenceSolver_->clone();
4616  // clear many things
4617  gutsOfDestructor2();
4618  // Reset cutoff
4619  // Solvers know about direction
4620  double direction = solver_->getObjSense();
4621  double value;
4622  solver_->getDblParam(OsiDualObjectiveLimit,value); 
4623  setCutoff(value*direction);
4624}
4625
4626// Are there a numerical difficulties?
4627bool 
4628CbcModel::isAbandoned() const
4629{
4630  return status_ == 2;
4631}
4632// Is optimality proven?
4633bool 
4634CbcModel::isProvenOptimal() const
4635{
4636  if (!status_ && bestObjective_<1.0e30)
4637    return true;
4638  else
4639    return false;
4640}
4641// Is  infeasiblity proven (or none better than cutoff)?
4642bool 
4643CbcModel::isProvenInfeasible() const
4644{
4645  if (!status_ && bestObjective_>=1.0e30)
4646    return true;
4647  else
4648    return false;
4649}
4650// Was continuous solution unbounded
4651bool 
4652CbcModel::isContinuousUnbounded() const
4653{
4654  if (!status_ && secondaryStatus_==7)
4655    return true;
4656  else
4657    return false;
4658}
4659// Was continuous solution unbounded
4660bool 
4661CbcModel::isProvenDualInfeasible() const
4662{
4663  if (!status_ && secondaryStatus_==7)
4664    return true;
4665  else
4666    return false;
4667}
4668// Node limit reached?
4669bool 
4670CbcModel::isNodeLimitReached() const
4671{
4672  return numberNodes_ >= intParam_[CbcMaxNumNode];
4673}
4674// Time limit reached?
4675bool 
4676CbcModel::isSecondsLimitReached() const
4677{
4678  if (status_==1&&secondaryStatus_==4)
4679    return true;
4680  else
4681    return false;
4682}
4683// Solution limit reached?
4684bool 
4685CbcModel::isSolutionLimitReached() const
4686{
4687  return numberSolutions_ >= intParam_[CbcMaxNumSol];
4688}
4689// Set language
4690void 
4691CbcModel::newLanguage(CoinMessages::Language language)
4692{
4693  messages_ = CbcMessage(language);
4694}
4695void 
4696CbcModel::setNumberStrong(int number)
4697{
4698  if (number<0)
4699    numberStrong_=0;
4700   else
4701    numberStrong_=number;
4702}
4703void 
4704CbcModel::setNumberBeforeTrust(int number)
4705{
4706  if (number<-3) {
4707    numberBeforeTrust_=0;
4708  } else {
4709    numberBeforeTrust_=number;
4710    //numberStrong_ = CoinMax(numberStrong_,1);
4711  }
4712}
4713void 
4714CbcModel::setNumberPenalties(int number)
4715{
4716  if (number<=0) {
4717    numberPenalties_=0;
4718  } else {
4719    numberPenalties_=number;
4720  }
4721}
4722void 
4723CbcModel::setPenaltyScaleFactor(double value)
4724{
4725  if (value<=0) {
4726    penaltyScaleFactor_=3.0;
4727  } else {
4728    penaltyScaleFactor_=value;
4729  }
4730}
4731void 
4732CbcModel::setHowOftenGlobalScan(int number)
4733{
4734  if (number<-1)
4735    howOftenGlobalScan_=0;
4736   else
4737    howOftenGlobalScan_=number;
4738}
4739
4740// Add one generator
4741void 
4742CbcModel::addCutGenerator(CglCutGenerator * generator,
4743                          int howOften, const char * name,
4744                          bool normal, bool atSolution,
4745                          bool whenInfeasible,int howOftenInSub,
4746                          int whatDepth, int whatDepthInSub)
4747{
4748  CbcCutGenerator ** temp = generator_;
4749  generator_ = new CbcCutGenerator * [numberCutGenerators_+1];
4750  memcpy(generator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
4751  delete[] temp ;
4752  generator_[numberCutGenerators_]= 
4753    new CbcCutGenerator(this,generator, howOften, name,
4754                        normal,atSolution,whenInfeasible,howOftenInSub,
4755                        whatDepth, whatDepthInSub);
4756  // and before any cahnges
4757  temp = virginGenerator_;
4758  virginGenerator_ = new CbcCutGenerator * [numberCutGenerators_+1];
4759  memcpy(virginGenerator_,temp,numberCutGenerators_*sizeof(CbcCutGenerator *));
4760  delete[] temp ;
4761  virginGenerator_[numberCutGenerators_++]= 
4762    new CbcCutGenerator(this,generator, howOften, name,
4763                        normal,atSolution,whenInfeasible,howOftenInSub,
4764                        whatDepth, whatDepthInSub);
4765                                                         
4766}
4767// Add one heuristic
4768void 
4769CbcModel::addHeuristic(CbcHeuristic * generator, const char *name)
4770{
4771  CbcHeuristic ** temp = heuristic_;
4772  heuristic_ = new CbcHeuristic * [numberHeuristics_+1];
4773  memcpy(heuristic_,temp,numberHeuristics_*sizeof(CbcHeuristic *));
4774  delete [] temp;
4775  heuristic_[numberHeuristics_]=generator->clone();
4776  if (name)
4777  { heuristic_[numberHeuristics_]->setHeuristicName(name) ; }
4778  heuristic_[numberHeuristics_]->setSeed(987654321+numberHeuristics_);
4779  numberHeuristics_++ ;
4780}
4781
4782/*
4783  The last subproblem handled by the solver is not necessarily related to the
4784  one being recreated, so the first action is to remove all cuts from the
4785  constraint system.  Next, traverse the tree from node to the root to
4786  determine the basis size required for this subproblem and create an empty
4787  basis with the right capacity.  Finally, traverse the tree from root to
4788  node, adjusting bounds in the constraint system, adjusting the basis, and
4789  collecting the cuts that must be added to the constraint system.
4790  applyToModel does the heavy lifting.
4791
4792  addCuts1 is used in contexts where all that's desired is the list of cuts:
4793  the node is already fathomed, and we're collecting cuts so that we can
4794  adjust reference counts as we prune nodes. Arguably the two functions
4795  should be separated. The culprit is applyToModel, which performs cut
4796  collection and model adjustment.
4797
4798  Certainly in the contexts where all we need is a list of cuts, there's no
4799  point in passing in a valid basis --- an empty basis will do just fine.
4800*/
4801void CbcModel::addCuts1 (CbcNode * node, CoinWarmStartBasis *&lastws)
4802{ int i;
4803  int nNode=0;
4804  int numberColumns = getNumCols();
4805  CbcNodeInfo * nodeInfo = node->nodeInfo();
4806
4807/*
4808  Remove all cuts from the constraint system.
4809  (original comment includes ``see note below for later efficiency'', but
4810  the reference isn't clear to me).
4811*/
4812  if (numberThreads_<=0) {
4813    solver_->restoreBaseModel(numberRowsAtContinuous_);
4814  } else {
4815    // *** Fix later
4816    int numberCuts = solver_->getNumRows()-numberRowsAtContinuous_;
4817    int *which = new int[numberCuts];
4818    for (i = 0 ; i < numberCuts ; i++)
4819      which[i] = i+numberRowsAtContinuous_;
4820    solver_->deleteRows(numberCuts,which);
4821    delete [] which;
4822  }
4823/*
4824  Accumulate the path from node to the root in walkback_, and accumulate a
4825  cut count in currentNumberCuts.
4826
4827  original comment: when working then just unwind until where new node joins
4828  old node (for cuts?)
4829*/
4830  int currentNumberCuts = 0;
4831  while (nodeInfo) {
4832    //printf("nNode = %d, nodeInfo = %x\n",nNode,nodeInfo);
4833    walkback_[nNode++]=nodeInfo;
4834    currentNumberCuts += nodeInfo->numberCuts() ;
4835    nodeInfo = nodeInfo->parent() ;
4836    if (nNode==maximumDepth_) {
4837      maximumDepth_ *= 2;
4838      CbcNodeInfo ** temp = new CbcNodeInfo * [maximumDepth_];
4839      for (i=0;i<nNode;i++) 
4840        temp[i] = walkback_[i];
4841      delete [] walkback_;
4842      walkback_ = temp;
4843    }
4844  }
4845/*
4846  Create an empty basis with sufficient capacity for the constraint system
4847  we'll construct: original system plus cuts. Make sure we have capacity to
4848  record those cuts in addedCuts_.
4849
4850  The method of adjusting the basis at a FullNodeInfo object (the root, for
4851  example) is to use a copy constructor to duplicate the basis held in the
4852  nodeInfo, then resize it and return the new basis object. Guaranteed,
4853  lastws will point to a different basis when it returns. We pass in a basis
4854  because we need the parameter to return the allocated basis, and it's an
4855  easy way to pass in the size. But we take a hit for memory allocation.
4856*/
4857  currentNumberCuts_=currentNumberCuts;
4858  if (currentNumberCuts > maximumNumberCuts_) {
4859    maximumNumberCuts_ = currentNumberCuts;
4860    delete [] addedCuts_;
4861    addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
4862  }
4863  lastws->setSize(numberColumns,numberRowsAtContinuous_+currentNumberCuts);
4864/*
4865  This last bit of code traverses the path collected in walkback_ from the
4866  root back to node. At the end of the loop,
4867   * lastws will be an appropriate basis for node;
4868   * variable bounds in the constraint system will be set to be correct for
4869     node; and
4870   * addedCuts_ will be set to a list of cuts that need to be added to the
4871     constraint system at node.
4872  applyToModel does all the heavy lifting.
4873*/
4874  currentNumberCuts=0;
4875  //#define CBC_PRINT2
4876#ifdef CBC_PRINT2
4877  printf("Starting bounds at node %d\n",numberNodes_);
4878#endif
4879  while (nNode) {
4880    --nNode;
4881    walkback_[nNode]->applyToModel(this,lastws,addedCuts_,currentNumberCuts);
4882  }
4883  if (0) {
4884    int numberDebugValues=18;
4885    double * debugValues = new double[numberDebugValues];
4886    CoinZeroN(debugValues,numberDebugValues);
4887    debugValues[1]=6.0;
4888    debugValues[3]=60.0;
4889    debugValues[4]=6.0;
4890    debugValues[6]=60.0;
4891    debugValues[7]=16.0;
4892    debugValues[9]=70.0;
4893    debugValues[10]=7.0;
4894    debugValues[12]=70.0;
4895    debugValues[13]=12.0;
4896    debugValues[15]=75.0;
4897    int nBad=0;
4898    for (int j=0;j<numberColumns;j++) {
4899      if (integerInfo_[j]) {
4900        if(solver_->getColLower()[j]>debugValues[j]||
4901           solver_->getColUpper()[j]<debugValues[j]) {
4902          printf("** (%g) ** ",debugValues[j]);
4903          nBad++;
4904        }
4905        printf("%d bounds %g %g\n",j,solver_->getColLower()[j],solver_->getColUpper()[j]);
4906      }
4907    }
4908    if (nBad)
4909      printf("%d BAD\n",nBad);
4910    else
4911      printf("OKAY\n");
4912    delete [] debugValues;
4913  }
4914}
4915
4916/*
4917  adjustCuts might be a better name: If the node is feasible, we sift through
4918  the cuts collected by addCuts1, add the ones that are tight and omit the
4919  ones that are loose. If the node is infeasible, we just adjust the
4920  reference counts to reflect that we're about to prune this node and its
4921  descendants.
4922*/
4923int CbcModel::addCuts (CbcNode *node, CoinWarmStartBasis *&lastws,bool canFix)
4924{
4925/*
4926  addCuts1 performs step 1 of restoring the subproblem at this node; see the
4927  comments there.
4928*/
4929  addCuts1(node,lastws);
4930  int i;
4931  int numberColumns = getNumCols();
4932  if (solver_->getNumRows()>maximumRows_) {
4933    maximumRows_ = solver_->getNumRows();
4934    workingBasis_.resize(maximumRows_,numberColumns);
4935  }
4936  CbcNodeInfo * nodeInfo = node->nodeInfo();
4937  double cutoff = getCutoff() ;
4938  int currentNumberCuts=currentNumberCuts_;
4939  if (canFix) {
4940    bool feasible=true;
4941    const double *lower = solver_->getColLower() ;
4942    const double *upper = solver_->getColUpper() ;
4943    double * newLower = analyzeResults_;
4944    double * objLower = newLower+numberIntegers_;
4945    double * newUpper = objLower+numberIntegers_;
4946    double * objUpper = newUpper+numberIntegers_;
4947    int n=0;
4948    for (i=0;i<numberIntegers_;i++) {
4949      int iColumn = integerVariable_[i];
4950      bool changed=false;
4951      double lo = 0.0;
4952      double up = 0.0;
4953      if (objLower[i]>cutoff) {
4954        lo = lower[iColumn];
4955        up = upper[iColumn];
4956        if (lo<newLower[i]) {
4957          lo = newLower[i];
4958          solver_->setColLower(iColumn,lo);
4959          changed=true;
4960          n++;
4961        }
4962        if (objUpper[i]>cutoff) {
4963          if (up>newUpper[i]) {
4964            up = newUpper[i];
4965            solver_->setColUpper(iColumn,up);
4966            changed=true;
4967            n++;
4968          }
4969        }
4970      } else if (objUpper[i]>cutoff) {
4971        lo = lower[iColumn];
4972        up = upper[iColumn];
4973        if (up>newUpper[i]) {
4974          up = newUpper[i];
4975          solver_->setColUpper(iColumn,up);
4976          changed=true;
4977          n++;
4978        }
4979      }
4980      if (changed&&lo>up) {
4981        feasible=false;
4982        break;
4983      }
4984    }
4985    if (!feasible) {
4986      printf("analysis says node infeas\n");
4987      cutoff=-COIN_DBL_MAX;
4988    }
4989  }
4990/*
4991  If the node can't be fathomed by bound, reinstall tight cuts in the
4992  constraint system. Even if there are no cuts, we'll want to set the
4993  reconstructed basis in the solver.
4994*/
4995  if (node->objectiveValue() < cutoff||numberThreads_)
4996  { 
4997#   ifdef CBC_CHECK_BASIS
4998    printf("addCuts: expanded basis; rows %d+%d\n",
4999           numberRowsAtContinuous_,currentNumberCuts);
5000    lastws->print();
5001#   endif
5002/*
5003  Adjust the basis and constraint system so that we retain only active cuts.
5004  There are three steps:
5005    1) Scan the basis. Sort the cuts into effective cuts to be kept and
5006       loose cuts to be dropped.
5007    2) Drop the loose cuts and resize the basis to fit.
5008    3) Install the tight cuts in the constraint system (applyRowCuts) and
5009       and install the basis (setWarmStart).
5010  Use of compressRows conveys we're compressing the basis and not just
5011  tweaking the artificialStatus_ array.
5012*/
5013    if (currentNumberCuts > 0) {
5014      int numberToAdd = 0;
5015      const OsiRowCut **addCuts;
5016      int numberToDrop = 0 ;
5017      int *cutsToDrop ;
5018      addCuts = new const OsiRowCut* [currentNumberCuts];
5019      cutsToDrop = new int[currentNumberCuts] ;
5020      assert (currentNumberCuts+numberRowsAtContinuous_<=lastws->getNumArtificial());
5021      for (i=0;i<currentNumberCuts;i++) {
5022        CoinWarmStartBasis::Status status = 
5023          lastws->getArtifStatus(i+numberRowsAtContinuous_);
5024        if (addedCuts_[i] &&
5025            (status != CoinWarmStartBasis::basic ||
5026             addedCuts_[i]->effectiveness()==COIN_DBL_MAX)) {
5027#         ifdef CHECK_CUT_COUNTS
5028          printf("Using cut %d %x as row %d\n",i,addedCuts_[i],
5029                 numberRowsAtContinuous_+numberToAdd);
5030#         endif
5031          addCuts[numberToAdd++] = addedCuts_[i];
5032        } else {
5033#         ifdef CHECK_CUT_COUNTS
5034          printf("Dropping cut %d %x\n",i,addedCuts_[i]);
5035#         endif
5036          addedCuts_[i]=NULL;
5037          cutsToDrop[numberToDrop++] = numberRowsAtContinuous_+i ;
5038        }
5039      }
5040      int numberRowsNow=numberRowsAtContinuous_+numberToAdd;
5041      lastws->compressRows(numberToDrop,cutsToDrop) ;
5042      lastws->resize(numberRowsNow,numberColumns);
5043      solver_->applyRowCuts(numberToAdd,addCuts);
5044#     ifdef CBC_CHECK_BASIS
5045      printf("addCuts: stripped basis; rows %d + %d\n",
5046             numberRowsAtContinuous_,numberToAdd);
5047      lastws->print();
5048#     endif
5049      //for (i=0;i<numberToAdd;i++)
5050      //delete addCuts[i];
5051      delete [] addCuts;
5052      delete [] cutsToDrop ;
5053    }
5054/*
5055  Set the basis in the solver.
5056*/
5057    solver_->setWarmStart(lastws);
5058#if 0
5059    if ((numberNodes_%printFrequency_)==0) {
5060      printf("Objective %g, depth %d, unsatisfied %d\n",
5061             node->objectiveValue(),
5062             node->depth(),node->numberUnsatisfied());
5063    }
5064#endif
5065/*
5066  Clean up and we're out of here.
5067*/
5068    numberNodes_++;
5069    return 0;
5070  } 
5071/*
5072  This node has been fathomed by bound as we try to revive it out of the live
5073  set. Adjust the cut reference counts to reflect that we no longer need to
5074  explore the remaining branch arms, hence they will no longer reference any
5075  cuts. Cuts whose reference count falls to zero are deleted. 
5076*/
5077  else
5078  { int i;
5079    if (currentNumberCuts) {
5080#ifndef CBC_DETERMINISTIC_THREAD
5081      lockThread();
5082#endif
5083      int numberLeft = nodeInfo->numberBranchesLeft();
5084      for (i = 0 ; i < currentNumberCuts ; i++)
5085        { if (addedCuts_[i])
5086          { if (!addedCuts_[i]->decrement(numberLeft))
5087            { delete addedCuts_[i];
5088            addedCuts_[i] = NULL; } } }
5089#ifndef CBC_DETERMINISTIC_THREAD
5090      unlockThread();
5091#endif
5092    }
5093    return 1 ; }
5094}
5095
5096
5097/*
5098  Perform reduced cost fixing on integer variables.
5099
5100  The variables in question are already nonbasic at bound. We're just nailing
5101  down the current situation.
5102*/
5103
5104int CbcModel::reducedCostFix ()
5105
5106{
5107  if(!solverCharacteristics_->reducedCostsAccurate())
5108    return 0; //NLP
5109  double cutoff = getCutoff() ;
5110  double direction = solver_->getObjSense() ;
5111  double gap = cutoff - solver_->getObjValue()*direction ;
5112  double tolerance;
5113  solver_->getDblParam(OsiDualTolerance,tolerance) ;
5114  if (gap<=0.0)
5115    return 0;
5116  gap += 100.0*tolerance;
5117  double integerTolerance = getDblParam(CbcIntegerTolerance) ;
5118
5119  const double *lower = solver_->getColLower() ;
5120  const double *upper = solver_->getColUpper() ;
5121  const double *solution = solver_->getColSolution() ;
5122  const double *reducedCost = solver_->getReducedCost() ;
5123
5124  int numberFixed = 0 ;
5125
5126# ifdef COIN_HAS_CLP
5127  OsiClpSolverInterface * clpSolver
5128    = dynamic_cast<OsiClpSolverInterface *> (solver_);
5129  ClpSimplex * clpSimplex=NULL;
5130  if (clpSolver) 
5131    clpSimplex = clpSolver->getModelPtr();
5132# endif
5133  for (int i = 0 ; i < numberIntegers_ ; i++)
5134  { int iColumn = integerVariable_[i] ;
5135    double djValue = direction*reducedCost[iColumn] ;
5136    if (upper[iColumn]-lower[iColumn] > integerTolerance)
5137    { if (solution[iColumn] < lower[iColumn]+integerTolerance && djValue > gap)
5138      { solver_->setColUpper(iColumn,lower[iColumn]) ;
5139#ifdef COIN_HAS_CLP
5140      // may just have been fixed before
5141      if (clpSimplex)
5142        assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atLowerBound||
5143               clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
5144#endif
5145      numberFixed++ ; }
5146      else
5147      if (solution[iColumn] > upper[iColumn]-integerTolerance && -djValue > gap)
5148      { solver_->setColLower(iColumn,upper[iColumn]) ;
5149#ifdef COIN_HAS_CLP
5150      // may just have been fixed before
5151      if (clpSimplex)
5152        assert(clpSimplex->getColumnStatus(iColumn)==ClpSimplex::atUpperBound||
5153               clpSimplex->getColumnStatus(iColumn)==ClpSimplex::isFixed);
5154#endif
5155      numberFixed++ ; } } }
5156  numberDJFixed_ += numberFixed;
5157  return numberFixed; }
5158
5159// Collect coding to replace whichGenerator
5160void
5161CbcModel::resizeWhichGenerator(int numberNow, int numberAfter)
5162{
5163  if (numberAfter > maximumWhich_) {
5164    maximumWhich_ = CoinMax(maximumWhich_*2+100,numberAfter) ;
5165    int * temp = new int[2*maximumWhich_] ;
5166    memcpy(temp,whichGenerator_,numberNow*sizeof(int)) ;
5167    delete [] whichGenerator_ ;
5168    whichGenerator_ = temp ;
5169    memset(whichGenerator_+numberNow,0,(maximumWhich_-numberNow)*sizeof(int));
5170  }
5171}
5172
5173/** Solve the model using cuts
5174
5175  This version takes off redundant cuts from node.
5176  Returns true if feasible.
5177
5178  \todo
5179  Why do I need to resolve the problem? What has been done between the last
5180  relaxation and calling solveWithCuts?
5181
5182  If numberTries == 0 then user did not want any cuts.
5183*/
5184
5185bool 
5186CbcModel::solveWithCuts (OsiCuts &cuts, int numberTries, CbcNode *node)
5187/*
5188  Parameters:
5189    numberTries: (i) the maximum number of iterations for this round of cut
5190                     generation; if negative then we don't mind if drop is tiny.
5191   
5192    cuts:       (o) all cuts generated in this round of cut generation
5193
5194    node: (i)     So we can update dynamic pseudo costs
5195*/
5196
5197
5198{
5199# ifdef COIN_HAS_CLP
5200  OsiClpSolverInterface * clpSolver
5201    = dynamic_cast<OsiClpSolverInterface *> (solver_);
5202  int saveClpOptions=0;
5203  if (clpSolver) 
5204    saveClpOptions = clpSolver->specialOptions();
5205# endif
5206  //solver_->writeMps("saved");
5207#ifdef CBC_THREAD
5208  CbcModel ** threadModel = NULL;
5209  Coin_pthread_t * threadId = NULL;
5210  pthread_cond_t condition_main;
5211  pthread_mutex_t condition_mutex;
5212  pthread_mutex_t * mutex2 = NULL;
5213  pthread_cond_t * condition2 = NULL;
5214  threadStruct * threadInfo = NULL;
5215  void * saveMutex = NULL;
5216  if (numberThreads_&&(threadMode_&2)!=0&&!numberNodes_) {
5217    threadId = new Coin_pthread_t [numberThreads_];
5218    pthread_cond_init(&condition_main,NULL);
5219    pthread_mutex_init(&condition_mutex,NULL);
5220    threadModel = new CbcModel * [numberThreads_];
5221    threadInfo = new threadStruct [numberThreads_+1];
5222    mutex2 = new pthread_mutex_t [numberThreads_];
5223    condition2 = new pthread_cond_t [numberThreads_];
5224    saveMutex = mutex_;
5225    for (int i=0;i<numberThreads_;i++) {
5226      pthread_mutex_init(mutex2+i,NULL);
5227      pthread_cond_init(condition2+i,NULL);
5228      threadId[i].status =0;
5229      threadModel[i]=new CbcModel;
5230      threadModel[i]->generator_ = new CbcCutGenerator * [1];
5231      delete threadModel[i]->solver_;
5232      threadModel[i]->solver_=NULL;
5233      threadModel[i]->numberThreads_=numberThreads_;
5234      mutex_ = (void *) (threadInfo+i);
5235      threadInfo[i].thisModel=(CbcModel *) threadModel[i];
5236      threadInfo[i].baseModel=this;
5237          threadInfo[i].threadIdOfBase.thr=pthread_self();
5238      threadInfo[i].mutex2=mutex2+i;
5239      threadInfo[i].condition2=condition2+i;
5240      threadInfo[i].returnCode=-1;
5241          pthread_create(&threadId[i].thr,NULL,doCutsThread,threadInfo+i);
5242          threadId[i].status = 1;
5243
5244    }
5245    // Do a partial one for base model
5246    threadInfo[numberThreads_].baseModel=this;
5247    mutex_ = (void *) (threadInfo+numberThreads_);
5248    threadInfo[numberThreads_].condition2=&condition_main;
5249    threadInfo[numberThreads_].mutex2=&condition_mutex;
5250  }
5251#endif
5252  bool feasible = true ;
5253  int lastNumberCuts = 0 ;
5254  double lastObjective = -1.0e100 ;
5255  int violated = 0 ;
5256  int numberRowsAtStart = solver_->getNumRows() ;
5257  //printf("solver had %d rows\n",numberRowsAtStart);
5258  int numberColumns = solver_->getNumCols() ;
5259  CoinBigIndex numberElementsAtStart = solver_->getNumElements();
5260
5261  numberOldActiveCuts_ = numberRowsAtStart-numberRowsAtContinuous_ ;
5262  numberNewCuts_ = 0 ;
5263
5264  bool onOptimalPath = false ;
5265  const OsiRowCutDebugger *debugger = NULL;
5266  if ((specialOptions_&1)!=0) {
5267    /*
5268      See OsiRowCutDebugger for details. In a nutshell, make sure that current
5269      variable values do not conflict with a known optimal solution. (Obviously
5270      this can be fooled when there are multiple solutions.)
5271    */
5272    debugger = solver_->getRowCutDebugger() ;
5273    if (debugger) 
5274      onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
5275  }
5276  OsiCuts slackCuts;
5277/*
5278  Resolve the problem. If we've lost feasibility, might as well bail out right
5279  after the debug stuff. The resolve will also refresh cached copies of the
5280  solver solution (cbcColLower_, ...) held by CbcModel.
5281*/
5282  double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
5283  if (node)
5284    objectiveValue= node->objectiveValue();
5285  int returnCode = resolve(node ? node->nodeInfo() : NULL,1);
5286#ifdef COIN_DEVELOP
5287  //if (!solver_->getIterationCount()&&solver_->isProvenOptimal())
5288  //printf("zero iterations on first solve of branch\n");
5289#endif
5290  if (node&&node->nodeInfo()&&!node->nodeInfo()->numberBranchesLeft())
5291    node->nodeInfo()->allBranchesGone(); // can clean up
5292  feasible = returnCode  != 0 ;
5293  if (returnCode<0)
5294    numberTries=0;
5295  if (problemFeasibility_->feasible(this,0)<0) {
5296    feasible=false; // pretend infeasible
5297  }
5298 
5299#if NEW_UPDATE_OBJECT==0
5300  // Update branching information if wanted
5301  if(node &&branchingMethod_)
5302    branchingMethod_->updateInformation(solver_,node);
5303#elif NEW_UPDATE_OBJECT<2
5304  // Update branching information if wanted
5305  if(node &&branchingMethod_) {
5306    OsiBranchingObject * bobj = node->modifiableBranchingObject();
5307    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
5308    if (cbcobj) {
5309      CbcObject * object = cbcobj->object();
5310      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
5311      object->updateInformation(update);
5312    } else {
5313      branchingMethod_->updateInformation(solver_,node);
5314    }
5315  }
5316#else
5317  // Update branching information if wanted
5318  if(node &&branchingMethod_) {
5319    OsiBranchingObject * bobj = node->modifiableBranchingObject();
5320    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
5321    if (cbcobj&&cbcobj->object()) {
5322      CbcObject * object = cbcobj->object();
5323      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
5324      // have to compute object number as not saved
5325      CbcSimpleInteger * simpleObject =
5326          dynamic_cast <CbcSimpleInteger *>(object) ;
5327      int iObject;
5328      int iColumn = simpleObject->columnNumber();
5329      for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
5330        simpleObject =
5331          dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ;
5332        if (simpleObject->columnNumber()==iColumn)
5333          break;
5334      }
5335      assert (iObject<numberObjects_);
5336      update.objectNumber_ = iObject;
5337      addUpdateInformation(update);
5338    } else {
5339      OsiIntegerBranchingObject * obj = dynamic_cast<OsiIntegerBranchingObject *> (bobj);
5340      if (obj) {
5341        const OsiObject * object = obj->originalObject();
5342        // have to compute object number as not saved
5343        int iObject;
5344        int iColumn = object->columnNumber();
5345        for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
5346          if (object_[iObject]->columnNumber()==iColumn)
5347            break;
5348        }
5349        assert (iObject<numberObjects_);
5350        int branch = obj->firstBranch();
5351        if (obj->branchIndex()==2)
5352          branch = 1-branch;
5353        assert (branch==0||branch==1);
5354        double originalValue=node->objectiveValue();
5355        double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
5356        double changeInObjective = CoinMax(0.0,objectiveValue-originalValue);
5357        int iStatus = (feasible) ? 0 : 0;
5358        double value = obj->value();
5359        double movement;
5360        if (branch)
5361          movement = ceil(value)-value;
5362        else
5363          movement = value -floor(value);
5364#if 0
5365        // OUT as much too complicated - we are not at a natural hotstart place
5366        OsiBranchingInformation usefulInfo=usefulInformation();
5367        // hotInfo is meant for BEFORE a branch so we need to fool
5368        // was much simpler with alternate method
5369        double save[3];
5370        save[0]=usefulInfo.lower_[iColumn];
5371        save[1]=usefulInfo.solution_[iColumn];
5372        save[2]=usefulInfo.upper_[iColumn];
5373        usefulInfo.lower_[iColumn]=floor(value);
5374        usefulInfo.solution_[iColumn]=value;
5375        usefulInfo.upper_[iColumn]=ceil(value);
5376        OsiHotInfo hotInfo(solver_,&usefulInfo,&object,0);
5377        usefulInfo.lower_[iColumn]=save[0];
5378        usefulInfo.solution_[iColumn]=save[1];
5379        usefulInfo.upper_[iColumn]=save[2];
5380        if (branch) {
5381          hotInfo.setUpStatus(iStatus);
5382          hotInfo.setUpChange(changeInObjective);
5383          //object->setUpEstimate(movement);
5384        } else {
5385          hotInfo.setDownStatus(iStatus);
5386          hotInfo.setDownChange(changeInObjective);
5387          //object->setDownEstimate(movement);
5388        }
5389        branchingMethod_->chooseMethod()->updateInformation(&usefulInfo,branch,&hotInfo);
5390#else
5391        branchingMethod_->chooseMethod()->updateInformation(iObject,branch,changeInObjective,
5392                                                            movement,iStatus);
5393#endif
5394      }
5395    }
5396  }
5397#endif
5398
5399#ifdef CBC_DEBUG
5400  if (feasible)
5401  { printf("Obj value %g (%s) %d rows\n",solver_->getObjValue(),
5402           (solver_->isProvenOptimal())?"proven":"unproven",
5403           solver_->getNumRows()) ; }
5404 
5405  else
5406  { printf("Infeasible %d rows\n",solver_->getNumRows()) ; }
5407#endif
5408  if ((specialOptions_&1)!=0) {
5409/*
5410  If the RowCutDebugger said we were compatible with the optimal solution,
5411  and now we're suddenly infeasible, we might be confused. Then again, we
5412  may have fathomed by bound, heading for a rediscovery of an optimal solution.
5413*/
5414    if (onOptimalPath && !solver_->isDualObjectiveLimitReached()) {
5415      if (!feasible) {
5416        solver_->writeMps("infeas");
5417        CoinWarmStartBasis *slack =
5418          dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
5419        solver_->setWarmStart(slack);
5420        delete slack ;
5421        solver_->setHintParam(OsiDoReducePrint,false,OsiHintDo,0) ;
5422        solver_->initialSolve();
5423      }
5424      assert(feasible) ;
5425    }
5426  }
5427
5428  if (!feasible) {
5429    numberInfeasibleNodes_++;
5430# ifdef COIN_HAS_CLP
5431    if (clpSolver) 
5432    clpSolver->setSpecialOptions(saveClpOptions);
5433# endif
5434    return (false) ;
5435  }
5436  sumChangeObjective1_ += solver_->getObjValue()*solver_->getObjSense()
5437    - objectiveValue ;
5438  if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] )
5439    numberTries=0; // exit
5440  //if ((numberNodes_%100)==0)
5441  //printf("XXa sum obj changed by %g\n",sumChangeObjective1_);
5442  objectiveValue = solver_->getObjValue()*solver_->getObjSense();
5443  // Return at once if numberTries zero
5444  if (!numberTries) {
5445    cuts=OsiCuts();
5446    numberNewCuts_=0;
5447# ifdef COIN_HAS_CLP
5448    if (clpSolver) 
5449    clpSolver->setSpecialOptions(saveClpOptions);
5450# endif
5451    return true;
5452  }
5453/*
5454  Do reduced cost fixing.
5455*/
5456  reducedCostFix() ;
5457/*
5458  Set up for at most numberTries rounds of cut generation. If numberTries is
5459  negative, we'll ignore the minimumDrop_ cutoff and keep generating cuts for
5460  the specified number of rounds.
5461*/
5462  double minimumDrop = minimumDrop_ ;
5463  if (numberTries<0)
5464  { numberTries = -numberTries ;
5465    minimumDrop = -1.0 ; }
5466/*
5467  Is it time to scan the cuts in order to remove redundant cuts? If so, set
5468  up to do it.
5469*/
5470# define SCANCUTS 100 
5471  int *countColumnCuts = NULL ;
5472  // Always accumulate row cut counts
5473  int * countRowCuts =new int[numberCutGenerators_+numberHeuristics_] ;
5474  memset(countRowCuts,0,
5475         (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ;
5476  bool fullScan = false ;
5477  if ((numberNodes_%SCANCUTS) == 0)
5478  { fullScan = true ;
5479    countColumnCuts = new int[numberCutGenerators_+numberHeuristics_] ;
5480    memset(countColumnCuts,0,
5481           (numberCutGenerators_+numberHeuristics_)*sizeof(int)) ; }
5482
5483  double direction = solver_->getObjSense() ;
5484  double startObjective = solver_->getObjValue()*direction ;
5485
5486  currentPassNumber_ = 0 ;
5487  double primalTolerance = 1.0e-7 ;
5488  // We may need to keep going on
5489  bool keepGoing=false;
5490/*
5491  Begin cut generation loop. Cuts generated during each iteration are
5492  collected in theseCuts. The loop can be divided into four phases:
5493   1) Prep: Fix variables using reduced cost. In the first iteration only,
5494      consider scanning globalCuts_ and activating any applicable cuts.
5495   2) Cut Generation: Call each generator and heuristic registered in the
5496      generator_ and heuristic_ arrays. Newly generated global cuts are
5497      copied to globalCuts_ at this time.
5498   3) Cut Installation and Reoptimisation: Install column and row cuts in
5499      the solver. Copy row cuts to cuts (parameter). Reoptimise.
5500   4) Cut Purging: takeOffCuts() removes inactive cuts from the solver, and
5501      does the necessary bookkeeping in the model.
5502*/
5503  do
5504  { currentPassNumber_++ ;
5505    numberTries-- ;
5506    if (numberTries<0&&keepGoing) {
5507      // switch off all normal ones
5508      for (int i = 0;i<numberCutGenerators_;i++) {
5509        if (!generator_[i]->mustCallAgain())
5510          generator_[i]->setSwitchedOff(true);
5511      }
5512    }
5513    keepGoing=false;
5514    OsiCuts theseCuts ;
5515/*
5516  Scan previously generated global column and row cuts to see if any are
5517  useful.
5518*/
5519    int numberViolated=0;
5520    if (currentPassNumber_ == 1 && howOftenGlobalScan_ > 0 &&
5521        (numberNodes_%howOftenGlobalScan_) == 0)
5522    { int numberCuts = globalCuts_.sizeColCuts() ;
5523      int i;
5524      // possibly extend whichGenerator
5525      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
5526      for ( i = 0 ; i < numberCuts ; i++)
5527      { OsiColCut *thisCut = globalCuts_.colCutPtr(i) ;
5528        if (thisCut->violated(cbcColSolution_)>primalTolerance) {
5529          printf("Global cut added - violation %g\n",
5530                 thisCut->violated(cbcColSolution_)) ;
5531          whichGenerator_[numberViolated++]=-1;
5532#ifndef GLOBAL_CUTS_JUST_POINTERS
5533          theseCuts.insert(*thisCut) ;
5534#else
5535          theseCuts.insert(thisCut) ;
5536#endif
5537        }
5538      }
5539      numberCuts = globalCuts_.sizeRowCuts() ;
5540      // possibly extend whichGenerator
5541      resizeWhichGenerator(numberViolated, numberViolated+numberCuts);
5542      for ( i = 0;i<numberCuts;i++) {
5543        OsiRowCut * thisCut = globalCuts_.rowCutPtr(i) ;
5544        if (thisCut->violated(cbcColSolution_)>primalTolerance) {
5545          //printf("Global cut added - violation %g\n",
5546          // thisCut->violated(cbcColSolution_)) ;
5547          whichGenerator_[numberViolated++]=-1;
5548#ifndef GLOBAL_CUTS_JUST_POINTERS
5549          theseCuts.insert(*thisCut) ;
5550#else
5551          theseCuts.insert(thisCut) ;
5552#endif
5553        }
5554      }
5555      numberGlobalViolations_+=numberViolated;
5556    }
5557/*
5558  Generate new cuts (global and/or local) and/or apply heuristics.  If
5559  CglProbing is used, then it should be first as it can fix continuous
5560  variables.
5561
5562  At present, CglProbing is the only case where generateCuts will return
5563  true. generateCuts actually modifies variable bounds in the solver when
5564  CglProbing indicates that it can fix a variable. Reoptimisation is required
5565  to take full advantage.
5566
5567  The need to resolve here should only happen after a heuristic solution.
5568  (Note default OSI implementation of optimalBasisIsAvailable always returns
5569  false.)
5570*/
5571    if (solverCharacteristics_->warmStart()&&
5572        !solver_->optimalBasisIsAvailable()) {
5573      //printf("XXXXYY no opt basis\n");
5574      resolve(node ? node->nodeInfo() : NULL,3);
5575    }
5576    if (nextRowCut_) {
5577      // branch was a cut - add it
5578      theseCuts.insert(*nextRowCut_);
5579      if (handler_->logLevel()>1)
5580        nextRowCut_->print();
5581      const OsiRowCut * cut=nextRowCut_;
5582      double lb = cut->lb();
5583      double ub = cut->ub();
5584      int n=cut->row().getNumElements();
5585      const int * column = cut->row().getIndices();
5586      const double * element = cut->row().getElements();
5587      double sum=0.0;
5588      for (int i=0;i<n;i++) {
5589        int iColumn = column[i];
5590        double value = element[i];
5591        //if (cbcColSolution_[iColumn]>1.0e-7)
5592        //printf("value of %d is %g\n",iColumn,cbcColSolution_[iColumn]);
5593        sum += value * cbcColSolution_[iColumn];
5594      }
5595      delete nextRowCut_;
5596      nextRowCut_=NULL;
5597      if (handler_->logLevel()>1)
5598        printf("applying branch cut, sum is %g, bounds %g %g\n",sum,lb,ub);
5599      // possibly extend whichGenerator
5600      resizeWhichGenerator(numberViolated, numberViolated+1);
5601      // set whichgenerator (also serves as marker to say don't delete0
5602      whichGenerator_[numberViolated++]=-2;
5603    }
5604
5605    // reset probing info
5606    //if (probingInfo_)
5607    //probingInfo_->initializeFixing();
5608    int i;
5609    if ((threadMode_&2)==0||numberNodes_) {
5610      for (i = 0;i<numberCutGenerators_;i++) {
5611        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
5612        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
5613        int numberRowCutsAfter = numberRowCutsBefore;
5614        int numberColumnCutsAfter = numberColumnCutsBefore;
5615        bool generate = generator_[i]->normal();
5616        // skip if not optimal and should be (maybe a cut generator has fixed variables)
5617        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
5618          generate=false;
5619        if (generator_[i]->switchedOff())
5620          generate=false;;
5621        if (generate) {
5622          bool mustResolve = 
5623            generator_[i]->generateCuts(theseCuts,fullScan,solver_,node) ;
5624          numberRowCutsAfter = theseCuts.sizeRowCuts() ;
5625          if(numberRowCutsBefore < numberRowCutsAfter &&
5626             generator_[i]->mustCallAgain())
5627            keepGoing=true; // say must go round
5628          // Check last cut to see if infeasible
5629          if(numberRowCutsBefore < numberRowCutsAfter) {
5630            const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
5631            if (thisCut->lb()>thisCut->ub()) {
5632              feasible = false; // sub-problem is infeasible
5633              break;
5634            }
5635          }
5636#ifdef CBC_DEBUG
5637          {
5638            int k ;
5639            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
5640              OsiRowCut thisCut = theseCuts.rowCut(k) ;
5641              /* check size of elements.
5642                 We can allow smaller but this helps debug generators as it
5643                 is unsafe to have small elements */
5644              int n=thisCut.row().getNumElements();
5645              const int * column = thisCut.row().getIndices();
5646              const double * element = thisCut.row().getElements();
5647              //assert (n);
5648              for (int i=0;i<n;i++) {
5649                double value = element[i];
5650                assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
5651              }
5652            }
5653          }
5654#endif
5655          if (mustResolve) {
5656            int returncode = resolve(node ? node->nodeInfo() : NULL,2);
5657            feasible = returnCode  != 0 ;
5658            if (returncode<0)
5659              numberTries=0;
5660            if ((specialOptions_&1)!=0) {
5661              debugger = solver_->getRowCutDebugger() ;
5662              if (debugger) 
5663                onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
5664              else
5665                onOptimalPath=false;
5666              if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
5667                assert(feasible) ;
5668            }
5669            if (!feasible)
5670              break ;
5671          }
5672        }
5673        numberRowCutsAfter = theseCuts.sizeRowCuts() ;
5674        numberColumnCutsAfter = theseCuts.sizeColCuts() ;
5675
5676        if ((specialOptions_&1)!=0) {
5677          if (onOptimalPath) {
5678            int k ;
5679            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
5680              OsiRowCut thisCut = theseCuts.rowCut(k) ;
5681              if(debugger->invalidCut(thisCut)) {
5682                solver_->writeMps("badCut");
5683#ifdef NDEBUG
5684                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
5685                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
5686                const double *lower = getColLower() ;
5687                const double *upper = getColUpper() ;
5688                int numberColumns = solver_->getNumCols();
5689                for (int i=0;i<numberColumns;i++)
5690                  printf("%d bounds %g,%g\n",i,lower[i],upper[i]);
5691                abort();
5692#endif
5693              }
5694              assert(!debugger->invalidCut(thisCut)) ;
5695            }
5696          }
5697        }
5698/*
5699  The cut generator has done its thing, and maybe it generated some
5700  cuts.  Do a bit of bookkeeping: load
5701  whichGenerator[i] with the index of the generator responsible for a cut,
5702  and place cuts flagged as global in the global cut pool for the model.
5703
5704  lastNumberCuts is the sum of cuts added in previous iterations; it's the
5705  offset to the proper starting position in whichGenerator.
5706*/
5707        int numberBefore =
5708          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
5709        int numberAfter =
5710          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
5711        // possibly extend whichGenerator
5712        resizeWhichGenerator(numberBefore, numberAfter);
5713        int j ;
5714        if (fullScan) {
5715          // counts
5716          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
5717        }
5718        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
5719
5720        bool dodgyCuts=false;
5721        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
5722          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
5723          if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) {
5724            dodgyCuts=true;
5725            break;
5726          }
5727          whichGenerator_[numberBefore++] = i ;
5728          if (thisCut->lb()>thisCut->ub())
5729            violated=-2; // sub-problem is infeasible
5730          if (thisCut->globallyValid()) {
5731            // add to global list
5732            OsiRowCut newCut(*thisCut);
5733            newCut.setGloballyValid(true);
5734            newCut.mutableRow().setTestForDuplicateIndex(false);
5735            globalCuts_.insert(newCut) ;
5736          }
5737        }
5738        if (dodgyCuts) {
5739          for (int k=numberRowCutsAfter-1;k>=j;k--) {
5740            const OsiRowCut * thisCut = theseCuts.rowCutPtr(k) ;
5741            if (thisCut->lb()>thisCut->ub())
5742              violated=-2; // sub-problem is infeasible
5743            if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) 
5744              theseCuts.eraseRowCut(k);
5745          }
5746          numberRowCutsAfter = theseCuts.sizeRowCuts() ;
5747          for (;j<numberRowCutsAfter;j++) {
5748            const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
5749            whichGenerator_[numberBefore++] = i ;
5750            if (thisCut->globallyValid()) {
5751              // add to global list
5752              OsiRowCut newCut(*thisCut);
5753              newCut.setGloballyValid(true);
5754              newCut.mutableRow().setTestForDuplicateIndex(false);
5755              globalCuts_.insert(newCut) ;
5756            }
5757          }
5758        }
5759        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
5760          whichGenerator_[numberBefore++] = i ;
5761          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
5762          if (thisCut->globallyValid()) {
5763            // add to global list
5764            OsiColCut newCut(*thisCut);
5765            newCut.setGloballyValid(true);
5766            globalCuts_.insert(newCut) ;
5767          }
5768        }
5769      }
5770      // Add in any violated saved cuts
5771      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
5772        int numberOld = theseCuts.sizeRowCuts()+lastNumberCuts;
5773        int numberCuts = slackCuts.sizeRowCuts() ;
5774        int i;
5775        // possibly extend whichGenerator
5776        resizeWhichGenerator(numberOld, numberOld+numberCuts);
5777        for ( i = 0;i<numberCuts;i++) {
5778          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
5779          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
5780            if (messageHandler()->logLevel()>2)
5781              printf("Old cut added - violation %g\n",
5782                     thisCut->violated(cbcColSolution_)) ;
5783            whichGenerator_[numberOld++]=-1;
5784            theseCuts.insert(*thisCut) ;
5785          }
5786        }
5787      }
5788    } else {
5789      // do cuts independently
5790      OsiCuts * eachCuts = new OsiCuts [numberCutGenerators_];;
5791#ifdef CBC_THREAD
5792      if (!threadModel) {
5793#endif
5794        // generate cuts
5795        for (i = 0;i<numberCutGenerators_;i++) {
5796          bool generate = generator_[i]->normal();
5797          // skip if not optimal and should be (maybe a cut generator has fixed variables)
5798          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
5799            generate=false;
5800          if (generator_[i]->switchedOff())
5801            generate=false;;
5802          if (generate) 
5803            generator_[i]->generateCuts(eachCuts[i],fullScan,solver_,node) ;
5804        }
5805#ifdef CBC_THREAD
5806      } else {
5807        for (i=0;i<numberThreads_;i++) {
5808          // set solver here after cloning
5809          threadModel[i]->solver_=solver_->clone();
5810          threadModel[i]->numberNodes_ = (fullScan) ? 1 : 0;
5811        }
5812        // generate cuts
5813        for (i = 0;i<numberCutGenerators_;i++) {
5814          bool generate = generator_[i]->normal();
5815          // skip if not optimal and should be (maybe a cut generator has fixed variables)
5816          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
5817            generate=false;
5818          if (generator_[i]->switchedOff())
5819            generate=false;;
5820          if (generate) { 
5821            bool finished=false;
5822            int iThread=-1;
5823            // see if any available
5824            for (iThread=0;iThread<numberThreads_;iThread++) {
5825              if (threadInfo[iThread].returnCode) {
5826                finished=true;
5827                break;
5828              } else if (threadInfo[iThread].returnCode==0) {