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

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

event handling

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