source: trunk/CbcModel.cpp @ 171

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

redow maxWhich

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