source: trunk/CbcModel.cpp @ 154

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

back to old node count

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