source: branches/devel/Cbc/src/CbcModel.cpp @ 477

Last change on this file since 477 was 477, checked in by forrest, 13 years ago

for nonlinear and start moving to OsiTree?
afor n

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