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

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

add ampl stuff after heuristic

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