source: trunk/CbcModel.cpp @ 186

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

heuristics

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