source: trunk/CbcModel.cpp @ 160

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

for Jon

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