source: trunk/CbcModel.cpp @ 165

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

CbcFeasibility?

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