source: trunk/CbcModel.cpp @ 201

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

mostly NDEBUG

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