source: trunk/CbcModel.cpp @ 182

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

more parameters

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