source: trunk/CbcModel.cpp @ 299

Last change on this file since 299 was 299, checked in by lou, 14 years ago

addCuts: addCuts: rewrote loop that drops loose cuts; compressRows makes
clear we're not just tweaking status. solveWithCuts: use CbcModel?'s cached
bounds & solution; automatically refreshed on resolve().
FULL_DEBUG -> CBC_CHECK_BASIS

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