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

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

try and fix heap problem

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