source: trunk/CbcModel.cpp @ 36

Last change on this file since 36 was 36, checked in by forrest, 17 years ago

start of CbcFathom? class

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