source: trunk/CbcModel.cpp @ 216

Last change on this file since 216 was 216, checked in by forrest, 15 years ago

add branching stuff

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