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

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

for nonlinear

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