source: trunk/CbcModel.cpp @ 162

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

issecondslimit

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