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

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

fix nullinfo bug

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