source: trunk/CbcModel.cpp @ 48

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

skip messages etc

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