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

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

towards common use with other solvers

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