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

Last change on this file since 616 was 616, checked in by forrest, 12 years ago

exit a bit earlier on easy models

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