source: trunk/Cbc/src/CbcModel.cpp @ 310

Last change on this file since 310 was 310, checked in by andreasw, 13 years ago

first commit for autotools conversion to be able to move more files

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