source: trunk/CbcModel.cpp @ 188

Last change on this file since 188 was 188, checked in by forrest, 14 years ago

local tree

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