source: trunk/CbcModel.cpp @ 129

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

major changes

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