source: trunk/CbcModel.cpp @ 273

Last change on this file since 273 was 273, checked in by forrest, 14 years ago

another stupidity

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