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

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

fixes for bonmin and other stuff

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