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

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

compiler error on MS

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