source: trunk/CbcModel.cpp @ 202

Last change on this file since 202 was 202, checked in by forrest, 16 years ago

try and improve

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