source: trunk/CbcModel.cpp @ 108

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

get around bad coding in CbcMain?

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