source: branches/devel/Cbc/src/CbcModel.cpp @ 648

Last change on this file since 648 was 648, checked in by forrest, 12 years ago

for parallel

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