source: trunk/CbcModel.cpp @ 195

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

stuff

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