source: trunk/CbcModel.cpp @ 277

Last change on this file since 277 was 277, checked in by lou, 14 years ago

Revise cbc build process for better control of solvers included in build
(COIN_USE_XXX -> CBC_USE_XXX). Add CbcEventHandler? for independence from
clp.

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