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

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

for OsiSOS and add more stats to cut generators

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 291.7 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      // save statistics
4800      for (i = 0;i<numberCutGenerators_;i++) {
4801        generator_[i]->setNumberCutsAtRoot(generator_[i]->numberCutsInTotal());
4802        generator_[i]->setNumberActiveCutsAtRoot(generator_[i]->numberCutsActive());
4803      }
4804      // decide on pseudo cost strategy
4805      int howOften = iProbing>=0 ? generator_[iProbing]->howOften() : 0;
4806      if ((howOften %1000000)!=1) 
4807        howOften = 0;
4808      //if (howOften) {
4809      //CglProbing * probing = dynamic_cast<CglProbing*>(generator_[iProbing]->generator());
4810      //}
4811      howOften = 0;
4812      if (howOften) {
4813        printf("** method 1\n");
4814        //CglProbing * probing = dynamic_cast<CglProbing*>(generator_[iProbing]->generator());
4815        generator_[iProbing]->setWhatDepth(1);
4816        // could set no row cuts
4817        //if (thisObjective-startObjective<0.001*fabs(startObjective)+1.0e-5)
4818        // probing->setRowCuts(0);
4819        for (int i=0;i<numberObjects_;i++) {
4820          CbcSimpleIntegerDynamicPseudoCost * obj =
4821            dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
4822          if (obj) 
4823            obj->setMethod(1);
4824        }
4825      }
4826      if (willBeCutsInTree==-2)
4827        willBeCutsInTree=0;
4828      if( willBeCutsInTree<=0) {
4829        // Take off cuts
4830        cuts = OsiCuts();
4831        numberNewCuts_=0;
4832        if (!willBeCutsInTree) {
4833          // update size of problem
4834          numberRowsAtContinuous_ = solver_->getNumRows() ;
4835        } else {
4836          // take off cuts
4837          int numberRows = solver_->getNumRows();
4838          int numberAdded = numberRows-numberRowsAtContinuous_;
4839          if (numberAdded) {
4840            int * added = new int[numberAdded];
4841            for (int i=0;i<numberAdded;i++)
4842              added[i]=i+numberRowsAtContinuous_;
4843            solver_->deleteRows(numberAdded,added);
4844            delete [] added;
4845            // resolve so optimal
4846            resolve(solver_);
4847          }
4848        }
4849#ifdef COIN_HAS_CLP
4850        OsiClpSolverInterface * clpSolver
4851          = dynamic_cast<OsiClpSolverInterface *> (solver_);
4852        if (clpSolver) {
4853          // Maybe solver might like to know only column bounds will change
4854          //int options = clpSolver->specialOptions();
4855          //clpSolver->setSpecialOptions(options|128);
4856          clpSolver->synchronizeModel();
4857        }
4858#endif
4859      } else {
4860#ifdef COIN_HAS_CLP
4861        OsiClpSolverInterface * clpSolver
4862          = dynamic_cast<OsiClpSolverInterface *> (solver_);
4863        if (clpSolver) {
4864        // make sure factorization can't carry over
4865          int options = clpSolver->specialOptions();
4866          clpSolver->setSpecialOptions(options&(~8));
4867        }
4868#endif
4869      }
4870    }
4871  } else if (numberCutGenerators_) {
4872    int i;
4873    // add to counts anyway
4874    for (i = 0;i<numberCutGenerators_;i++) 
4875      generator_[i]->incrementNumberCutsInTotal(countRowCuts[i]);
4876    // What if not feasible as cuts may have helped
4877    if (feasible) {
4878      for (i = 0;i<numberNewCuts_;i++) {
4879        int iGenerator = whichGenerator_[i];
4880        if (iGenerator>=0)
4881          generator_[iGenerator]->incrementNumberCutsActive();
4882      }
4883    }
4884  }
4885
4886  delete [] countRowCuts ;
4887  delete [] countColumnCuts ;
4888
4889
4890#ifdef CHECK_CUT_COUNTS
4891  if (feasible)
4892  { 
4893    CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4894    printf("solveWithCuts: Number of rows at end (only active cuts) %d\n",
4895           numberRowsAtContinuous_+numberNewCuts_+numberOldActiveCuts_) ;
4896    basis->print() ; delete basis;}
4897#endif
4898#ifdef CBC_DEBUG
4899  if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
4900    assert(feasible) ;
4901#endif
4902# ifdef COIN_HAS_CLP
4903  if (clpSolver) 
4904    clpSolver->setSpecialOptions(saveClpOptions);
4905# endif
4906
4907  return feasible ; }
4908
4909
4910/*
4911  Remove slack cuts. We obtain a basis and scan it. Cuts with basic slacks
4912  are purged. If any cuts are purged, resolve() is called to restore the
4913  solution held in the solver.  If resolve() pivots, there's the possibility
4914  that a slack may be pivoted in (trust me :-), so the process iterates.
4915  Setting allowResolve to false will suppress reoptimisation (but see note
4916  below).
4917
4918  At the level of the solver's constraint system, loose cuts are really
4919  deleted.  There's an implicit assumption that deleteRows will also update
4920  the active basis in the solver.
4921
4922  At the level of nodes and models, it's more complicated.
4923
4924  New cuts exist only in the collection of cuts passed as a parameter. They
4925  are deleted from the collection and that's the end of them.
4926
4927  Older cuts have made it into addedCuts_. Two separate actions are needed.
4928  The reference count for the CbcCountRowCut object is decremented. If this
4929  count falls to 0, the node which owns the cut is located, the reference to
4930  the cut is removed, and then the cut object is destroyed (courtesy of the
4931  CbcCountRowCut destructor). We also need to set the addedCuts_ entry to
4932  NULL. This is important so that when it comes time to generate basis edits
4933  we can tell this cut was dropped from the basis during processing of the
4934  node.
4935
4936  NOTE: In general, it's necessary to call resolve() after purging slack
4937        cuts.  Deleting constraints constitutes a change in the problem, and
4938        an OSI is not required to maintain a valid solution when the problem
4939        is changed. But ... it's really useful to maintain the active basis,
4940        and the OSI is supposed to do that. (Yes, it's splitting hairs.) In
4941        some places, it's possible to know that the solution will never be
4942        consulted after this call, only the basis.  (E.g., this routine is
4943        called as a last act before generating info to place the node in the
4944        live set.) For such use, set allowResolve to false.
4945 
4946  TODO: No real harm would be done if we just ignored the rare occasion when
4947        the call to resolve() pivoted a slack back into the basis. It's a
4948        minor inefficiency, at worst. But it does break assertions which
4949        check that there are no loose cuts in the basis. It might be better
4950        to remove the assertions.
4951*/
4952
4953void
4954CbcModel::takeOffCuts (OsiCuts &newCuts,
4955                       bool allowResolve, OsiCuts * saveCuts)
4956
4957{ // int resolveIterations = 0 ;
4958  int firstOldCut = numberRowsAtContinuous_ ;
4959  int totalNumberCuts = numberNewCuts_+numberOldActiveCuts_ ;
4960  int *solverCutIndices = new int[totalNumberCuts] ;
4961  int *newCutIndices = new int[numberNewCuts_] ;
4962  const CoinWarmStartBasis* ws ;
4963  CoinWarmStartBasis::Status status ;
4964  bool needPurge = true ;
4965/*
4966  The outer loop allows repetition of purge in the event that reoptimisation
4967  changes the basis. To start an iteration, clear the deletion counts and grab
4968  the current basis.
4969*/
4970  while (needPurge)
4971  { int numberNewToDelete = 0 ;
4972    int numberOldToDelete = 0 ;
4973    int i ;
4974    ws = dynamic_cast<const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
4975/*
4976  Scan the basis entries of the old cuts generated prior to this round of cut
4977  generation.  Loose cuts are `removed' by decrementing their reference count
4978  and setting the addedCuts_ entry to NULL. (If the reference count falls to
4979  0, they're really deleted.  See CbcModel and CbcCountRowCut doc'n for
4980  principles of cut handling.)
4981*/
4982    int oldCutIndex = 0 ;
4983    for (i = 0 ; i < numberOldActiveCuts_ ; i++)
4984    { status = ws->getArtifStatus(i+firstOldCut) ;
4985      while (!addedCuts_[oldCutIndex]) oldCutIndex++ ;
4986      assert(oldCutIndex < currentNumberCuts_) ;
4987      // always leave if from nextRowCut_
4988      if (status == CoinWarmStartBasis::basic&&
4989          addedCuts_[oldCutIndex]->effectiveness()!=COIN_DBL_MAX)
4990      { solverCutIndices[numberOldToDelete++] = i+firstOldCut ;
4991        if (saveCuts) {
4992          // send to cut pool
4993          OsiRowCut * slackCut = addedCuts_[oldCutIndex];
4994          if (slackCut->effectiveness()!=-1.234) {
4995            slackCut->setEffectiveness(-1.234);
4996            saveCuts->insert(*slackCut);
4997          }
4998        }
4999        if (addedCuts_[oldCutIndex]->decrement() == 0)
5000          delete addedCuts_[oldCutIndex] ;
5001        addedCuts_[oldCutIndex] = NULL ;
5002        oldCutIndex++ ; }
5003      else
5004      { oldCutIndex++ ; } }
5005/*
5006  Scan the basis entries of the new cuts generated with this round of cut
5007  generation.  At this point, newCuts is the only record of the new cuts, so
5008  when we delete loose cuts from newCuts, they're really gone. newCuts is a
5009  vector, so it's most efficient to compress it (eraseRowCut) from back to
5010  front.
5011*/
5012    int firstNewCut = firstOldCut+numberOldActiveCuts_ ;
5013    int k = 0 ;
5014    for (i = 0 ; i < numberNewCuts_ ; i++)
5015    { status = ws->getArtifStatus(i+firstNewCut) ;
5016      if (status == CoinWarmStartBasis::basic&&whichGenerator_[i]!=-2)
5017      { solverCutIndices[numberNewToDelete+numberOldToDelete] = i+firstNewCut ;
5018        newCutIndices[numberNewToDelete++] = i ; }
5019      else
5020      { // save which generator did it
5021        whichGenerator_[k++] = whichGenerator_[i] ; } }
5022    delete ws ;
5023    for (i = numberNewToDelete-1 ; i >= 0 ; i--)
5024    { int iCut = newCutIndices[i] ;
5025      if (saveCuts) {
5026        // send to cut pool
5027        OsiRowCut * slackCut = newCuts.rowCutPtr(iCut);
5028        if (slackCut->effectiveness()!=-1.234) {
5029          slackCut->setEffectiveness(-1.234);
5030          saveCuts->insert(*slackCut);
5031        }
5032      }
5033      newCuts.eraseRowCut(iCut) ; }
5034/*
5035  Did we delete anything? If so, delete the cuts from the constraint system
5036  held in the solver and reoptimise unless we're forbidden to do so. If the
5037  call to resolve() results in pivots, there's the possibility we again have
5038  basic slacks. Repeat the purging loop.
5039*/
5040    if (numberNewToDelete+numberOldToDelete > 0)
5041    { solver_->deleteRows(numberNewToDelete+numberOldToDelete,
5042                          solverCutIndices) ;
5043      numberNewCuts_ -= numberNewToDelete ;
5044      numberOldActiveCuts_ -= numberOldToDelete ;
5045#     ifdef CBC_DEBUG
5046      printf("takeOffCuts: purged %d+%d cuts\n", numberOldToDelete,
5047             numberNewToDelete );
5048#     endif
5049      if (allowResolve)
5050      { 
5051        phase_=3;
5052        // can do quick optimality check
5053        int easy=2;
5054        solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
5055        resolve(solver_) ;
5056        setPointers(solver_);
5057        solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
5058        if (solver_->getIterationCount() == 0)
5059        { needPurge = false ; }
5060#       ifdef CBC_DEBUG
5061        else
5062          { printf( "Repeating purging loop. %d iters.\n",
5063                    solver_->getIterationCount()); }
5064#       endif
5065      }
5066      else
5067      { needPurge = false ; } }
5068    else
5069    { needPurge = false ; } }
5070/*
5071  Clean up and return.
5072*/
5073  delete [] solverCutIndices ;
5074  delete [] newCutIndices ;
5075}
5076
5077
5078
5079int
5080CbcModel::resolve(CbcNodeInfo * parent, int whereFrom)
5081{
5082  // We may have deliberately added in violated cuts - check to avoid message
5083  int iRow;
5084  int numberRows = solver_->getNumRows();
5085  const double * rowLower = solver_->getRowLower();
5086  const double * rowUpper = solver_->getRowUpper();
5087  bool feasible=true;
5088  for (iRow= numberRowsAtContinuous_;iRow<numberRows;iRow++) {
5089    if (rowLower[iRow]>rowUpper[iRow]+1.0e-8)
5090      feasible=false;
5091  }
5092  // Can't happen if strong branching as would have been found before
5093  if (!numberStrong_&&numberObjects_>numberIntegers_) {
5094    int iColumn;
5095    int numberColumns = solver_->getNumCols();
5096    const double * columnLower = solver_->getColLower();
5097    const double * columnUpper = solver_->getColUpper();
5098    for (iColumn= 0;iColumn<numberColumns;iColumn++) {
5099      if (columnLower[iColumn]>columnUpper[iColumn]+1.0e-5)
5100        feasible=false;
5101    }
5102  }
5103#if 0
5104  {
5105    int iColumn;
5106    int numberColumns = solver_->getNumCols();
5107    const double * columnLower = solver_->getColLower();
5108    const double * columnUpper = solver_->getColUpper();
5109    const double * sol = solver_->getColSolution();
5110    bool feasible=true;
5111    double xx[37];
5112    memset(xx,0,sizeof(xx));
5113    xx[6]=xx[7]=xx[9]=xx[18]=xx[26]=xx[36]=1.0;
5114    for (iColumn= 0;iColumn<numberColumns;iColumn++) {
5115      if (solver_->isInteger(iColumn)) {
5116        //printf("yy %d at %g (bounds %g, %g)\n",iColumn,
5117        //     sol[iColumn],columnLower[iColumn],columnUpper[iColumn]);
5118        if (columnLower[iColumn]>xx[iColumn]||
5119            columnUpper[iColumn]<xx[iColumn])
5120          feasible=false;
5121        //solver_->setColLower(iColumn,CoinMin(xx[iColumn],columnUpper[iColumn]));
5122        //solver_->setColUpper(iColumn,CoinMax(xx[iColumn],columnLower[iColumn]));
5123      }
5124    }
5125    if (feasible)
5126      printf("On Feasible path\n");
5127  }
5128#endif
5129/*
5130  Reoptimize. Consider the possibility that we should fathom on bounds. But be
5131  careful --- where the objective takes on integral values, we may want to keep
5132  a solution where the objective is right on the cutoff.
5133*/
5134  if (feasible)
5135    {
5136      resolve(solver_) ;
5137      numberIterations_ += solver_->getIterationCount() ;
5138      feasible = (solver_->isProvenOptimal() &&
5139                  !solver_->isDualObjectiveLimitReached()) ;
5140    }
5141  if (0&&feasible) {
5142    const double * lb = solver_->getColLower();
5143    const double * ub = solver_->getColUpper();
5144    const double * x = solver_->getColSolution();
5145    const double * dj = solver_->getReducedCost();
5146    int numberColumns = solver_->getNumCols();
5147    for (int i=0;i<numberColumns;i++) {
5148      if (dj[i]>1.0e-4&&ub[i]-lb[i]>1.0e-4&&x[i]>lb[i]+1.0e-4)
5149        printf("error %d %g %g %g %g\n",i,dj[i],lb[i],x[i],ub[i]);
5150      if (dj[i]<-1.0e-4&&ub[i]-lb[i]>1.0e-4&&x[i]<ub[i]-1.0e-4)
5151        printf("error %d %g %g %g %g\n",i,dj[i],lb[i],x[i],ub[i]);
5152    }
5153  } 
5154  if (!feasible&& continuousObjective_ <-1.0e30) {
5155    // at root node - double double check
5156    bool saveTakeHint;
5157    OsiHintStrength saveStrength;
5158    solver_->getHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength);
5159    if (saveTakeHint||saveStrength==OsiHintIgnore) {
5160      solver_->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
5161      resolve(solver_);
5162      solver_->setHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength);
5163      numberIterations_ += solver_->getIterationCount() ;
5164      feasible = solver_->isProvenOptimal();
5165      //      solver_->writeMps("infeas");
5166    }
5167  }
5168  setPointers(solver_);
5169  int returnStatus = feasible ? 1 : 0;
5170  if (strategy_) {
5171    // user can play clever tricks here
5172    int status = strategy_->status(this,parent,whereFrom);
5173    if (status>=0) {
5174      if (status==0)
5175        returnStatus = 1;
5176      else if(status==1)
5177        returnStatus=-1;
5178      else
5179        returnStatus=0;
5180    }
5181  }
5182  return returnStatus ;
5183}
5184
5185
5186/* Set up objects.  Only do ones whose length is in range.
5187   If makeEquality true then a new model may be returned if
5188   modifications had to be made, otherwise "this" is returned.
5189
5190   Could use Probing at continuous to extend objects
5191*/
5192CbcModel * 
5193CbcModel::findCliques(bool makeEquality,
5194                      int atLeastThisMany, int lessThanThis, int defaultValue)
5195{
5196  // No objects are allowed to exist
5197  assert(numberObjects_==numberIntegers_||!numberObjects_);
5198  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
5199  int numberRows = solver_->getNumRows();
5200  int numberColumns = solver_->getNumCols();
5201
5202  // We may want to add columns
5203  int numberSlacks=0;
5204  int * rows = new int[numberRows];
5205  double * element =new double[numberRows];
5206
5207  int iRow;
5208
5209  findIntegers(true);
5210  numberObjects_=numberIntegers_;
5211
5212  int numberCliques=0;
5213  OsiObject ** object = new OsiObject * [numberRows];
5214  int * which = new int[numberIntegers_];
5215  char * type = new char[numberIntegers_];
5216  int * lookup = new int[numberColumns];
5217  int i;
5218  for (i=0;i<numberColumns;i++) 
5219    lookup[i]=-1;
5220  for (i=0;i<numberIntegers_;i++) 
5221    lookup[integerVariable_[i]]=i;
5222
5223  // Statistics
5224  int totalP1=0,totalM1=0;
5225  int numberBig=0,totalBig=0;
5226  int numberFixed=0;
5227
5228  // Row copy
5229  const double * elementByRow = matrixByRow.getElements();
5230  const int * column = matrixByRow.getIndices();
5231  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5232  const int * rowLength = matrixByRow.getVectorLengths();
5233
5234  // Column lengths for slacks
5235  const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
5236
5237  const double * lower = getColLower();
5238  const double * upper = getColUpper();
5239  const double * rowLower = getRowLower();
5240  const double * rowUpper = getRowUpper();
5241
5242  for (iRow=0;iRow<numberRows;iRow++) {
5243    int numberP1=0, numberM1=0;
5244    int j;
5245    double upperValue=rowUpper[iRow];
5246    double lowerValue=rowLower[iRow];
5247    bool good=true;
5248    int slack = -1;
5249    for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
5250      int iColumn = column[j];
5251      int iInteger=lookup[iColumn];
5252      if (upper[iColumn]-lower[iColumn]<1.0e-8) {
5253        // fixed
5254        upperValue -= lower[iColumn]*elementByRow[j];
5255        lowerValue -= lower[iColumn]*elementByRow[j];
5256        continue;
5257      } else if (upper[iColumn]!=1.0||lower[iColumn]!=0.0) {
5258        good = false;
5259        break;
5260      } else if (iInteger<0) {
5261        good = false;
5262        break;
5263      } else {
5264        if (columnLength[iColumn]==1)
5265          slack = iInteger;
5266      }
5267      if (fabs(elementByRow[j])!=1.0) {
5268        good=false;
5269        break;
5270      } else if (elementByRow[j]>0.0) {
5271        which[numberP1++]=iInteger;
5272      } else {
5273        numberM1++;
5274        which[numberIntegers_-numberM1]=iInteger;
5275      }
5276    }
5277    int iUpper = (int) floor(upperValue+1.0e-5);
5278    int iLower = (int) ceil(lowerValue-1.0e-5);
5279    int state=0;
5280    if (upperValue<1.0e6) {
5281      if (iUpper==1-numberM1)
5282        state=1;
5283      else if (iUpper==-numberM1)
5284        state=2;
5285      else if (iUpper<-numberM1)
5286        state=3;
5287    }
5288    if (!state&&lowerValue>-1.0e6) {
5289      if (-iLower==1-numberP1)
5290        state=-1;
5291      else if (-iLower==-numberP1)
5292        state=-2;
5293      else if (-iLower<-numberP1)
5294        state=-3;
5295    }
5296    if (good&&state) {
5297      if (abs(state)==3) {
5298        // infeasible
5299        numberObjects_=-1;
5300        break;
5301      } else if (abs(state)==2) {
5302        // we can fix all
5303        numberFixed += numberP1+numberM1;
5304        if (state>0) {
5305          // fix all +1 at 0, -1 at 1
5306          for (i=0;i<numberP1;i++)
5307            solver_->setColUpper(integerVariable_[which[i]],0.0);
5308          for (i=0;i<numberM1;i++)
5309            solver_->setColLower(integerVariable_[which[numberIntegers_-i-1]],
5310                                 1.0);
5311        } else {
5312          // fix all +1 at 1, -1 at 0
5313          for (i=0;i<numberP1;i++)
5314            solver_->setColLower(integerVariable_[which[i]],1.0);
5315          for (i=0;i<numberM1;i++)
5316            solver_->setColUpper(integerVariable_[which[numberIntegers_-i-1]],
5317                                 0.0);
5318        }
5319      } else {
5320        int length = numberP1+numberM1;
5321        if (length >= atLeastThisMany&&length<lessThanThis) {
5322          // create object
5323          bool addOne=false;
5324          int objectType;
5325          if (iLower==iUpper) {
5326            objectType=1;
5327          } else {
5328            if (makeEquality) {
5329              objectType=1;
5330              element[numberSlacks]=state;
5331              rows[numberSlacks++]=iRow;
5332              addOne=true;
5333            } else {
5334              objectType=0;
5335            }
5336          }
5337          if (state>0) {
5338            totalP1 += numberP1;
5339            totalM1 += numberM1;
5340            for (i=0;i<numberP1;i++)
5341              type[i]=1;
5342            for (i=0;i<numberM1;i++) {
5343              which[numberP1]=which[numberIntegers_-i-1];
5344              type[numberP1++]=0;
5345            }
5346          } else {
5347            totalP1 += numberM1;
5348            totalM1 += numberP1;
5349            for (i=0;i<numberP1;i++)
5350              type[i]=0;
5351            for (i=0;i<numberM1;i++) {
5352              which[numberP1]=which[numberIntegers_-i-1];
5353              type[numberP1++]=1;
5354            }
5355          }
5356          if (addOne) {
5357            // add in slack
5358            which[numberP1]=numberIntegers_+numberSlacks-1;
5359            slack = numberP1;
5360            type[numberP1++]=1;
5361          } else if (slack >= 0) {
5362            for (i=0;i<numberP1;i++) {
5363              if (which[i]==slack) {
5364                slack=i;
5365              }
5366            }
5367          }
5368          object[numberCliques] = new CbcClique(this,objectType,numberP1,
5369                                              which,type,
5370                                               1000000+numberCliques,slack);
5371          numberCliques++;
5372        } else if (numberP1+numberM1 >= lessThanThis) {
5373          // too big
5374          numberBig++;
5375          totalBig += numberP1+numberM1;
5376        }
5377      }
5378    }
5379  }
5380  delete [] which;
5381  delete [] type;
5382  delete [] lookup;
5383  if (numberCliques<0) {
5384    printf("*** Problem infeasible\n");
5385  } else {
5386    if (numberCliques)
5387      printf("%d cliques of average size %g found, %d P1, %d M1\n",
5388             numberCliques,
5389             ((double)(totalP1+totalM1))/((double) numberCliques),
5390             totalP1,totalM1);
5391    else
5392      printf("No cliques found\n");
5393    if (numberBig)
5394      printf("%d large cliques ( >= %d) found, total %d\n",
5395             numberBig,lessThanThis,totalBig);
5396    if (numberFixed)
5397      printf("%d variables fixed\n",numberFixed);
5398  }
5399  if (numberCliques>0&&numberSlacks&&makeEquality) {
5400    printf("adding %d integer slacks\n",numberSlacks);
5401    // add variables to make equality rows
5402    int * temp = new int[numberIntegers_+numberSlacks];
5403    memcpy(temp,integerVariable_,numberIntegers_*sizeof(int));
5404    // Get new model
5405    CbcModel * newModel = new CbcModel(*this);
5406    OsiSolverInterface * newSolver = newModel->solver();
5407    for (i=0;i<numberSlacks;i++) {
5408      temp[i+numberIntegers_]=i+numberColumns;
5409      int iRow = rows[i];
5410      double value = element[i];
5411      double lowerValue = 0.0;
5412      double upperValue = 1.0;
5413      double objValue  = 0.0;
5414      CoinPackedVector column(1,&iRow,&value);
5415      newSolver->addCol(column,lowerValue,upperValue,objValue);
5416      // set integer
5417      newSolver->setInteger(numberColumns+i);
5418      if (value >0)
5419        newSolver->setRowLower(iRow,rowUpper[iRow]);
5420      else
5421        newSolver->setRowUpper(iRow,rowLower[iRow]);
5422    }
5423    // replace list of integers
5424    for (i=0;i<newModel->numberObjects_;i++)
5425      delete newModel->object_[i];
5426    newModel->numberObjects_ = 0;
5427    delete [] newModel->object_;
5428    newModel->object_=NULL;
5429    newModel->findIntegers(true); //Set up all integer objects
5430    for (i=0;i<numberIntegers_;i++) {
5431      newModel->modifiableObject(i)->setPriority(object_[i]->priority());
5432    }
5433    if (originalColumns_) {
5434      // old model had originalColumns
5435      delete [] newModel->originalColumns_;
5436      newModel->originalColumns_ = new int[numberColumns+numberSlacks];
5437      memcpy(newModel->originalColumns_,originalColumns_,numberColumns*sizeof(int));
5438      // mark as not in previous model
5439      for (i=numberColumns;i<numberColumns+numberSlacks;i++)
5440        newModel->originalColumns_[i]=-1;
5441    }
5442    delete [] rows;
5443    delete [] element;
5444    newModel->addObjects(numberCliques,object);
5445    for (;i<numberCliques;i++) 
5446      delete object[i];
5447    delete [] object;
5448    newModel->synchronizeModel();
5449    return newModel;
5450  } else {
5451    if (numberCliques>0) {
5452      addObjects(numberCliques,object);
5453      for (;i<numberCliques;i++) 
5454        delete object[i];
5455      synchronizeModel();
5456    }
5457    delete [] object;
5458    delete [] rows;
5459    delete [] element;
5460    return this;
5461  }
5462}
5463// Fill in useful estimates
5464void 
5465CbcModel::pseudoShadow(double * down, double * up)
5466{
5467  // Column copy of matrix
5468  const double * element = solver_->getMatrixByCol()->getElements();
5469  const int * row = solver_->getMatrixByCol()->getIndices();
5470  const CoinBigIndex * columnStart = solver_->getMatrixByCol()->getVectorStarts();
5471  const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
5472  const double *objective = solver_->getObjCoefficients() ;
5473  int numberColumns = solver_->getNumCols() ;
5474  double direction = solver_->getObjSense();
5475  int iColumn;
5476  const double * dual = cbcRowPrice_;
5477  down = new double[numberColumns];
5478  up = new double[numberColumns];
5479  double upSum=1.0e-20;
5480  double downSum = 1.0e-20;
5481  int numberIntegers=0;
5482  for (iColumn=0;iColumn<numberColumns;iColumn++) {
5483    CoinBigIndex start = columnStart[iColumn];
5484    CoinBigIndex end = start + columnLength[iColumn];
5485    double upValue = 0.0;
5486    double downValue = 0.0;
5487    double value = direction*objective[iColumn];
5488    if (value) {
5489      if (value>0.0)
5490        upValue += value;
5491      else
5492        downValue -= value;
5493    }
5494    for (CoinBigIndex j=start;j<end;j++) {
5495      int iRow = row[j];
5496      value = -dual[iRow];
5497      if (value) {
5498        value *= element[j];
5499        if (value>0.0)
5500          upValue += value;
5501        else
5502          downValue -= value;
5503      }
5504    }
5505    // use dj if bigger
5506    double dj = cbcReducedCost_[iColumn];
5507    upValue = CoinMax(upValue,dj);
5508    downValue = CoinMax(downValue,-dj);
5509    up[iColumn]=upValue;
5510    down[iColumn]=downValue;
5511    if (solver_->isInteger(iColumn)) {
5512      if (!numberNodes_)
5513        printf("%d - dj %g up %g down %g cost %g\n",
5514               iColumn,dj,upValue,downValue,objective[iColumn]);
5515      upSum += upValue;
5516      downSum += downValue;
5517      numberIntegers++;
5518    }
5519  }
5520  if (numberIntegers) {
5521    double smallDown = 0.01*(downSum/((double) numberIntegers));
5522    double smallUp = 0.01*(upSum/((double) numberIntegers));
5523    for (int i=0;i<numberObjects_;i++) {
5524      CbcSimpleIntegerDynamicPseudoCost * obj1 =
5525        dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[i]) ;
5526      if (obj1) {
5527        iColumn = obj1->columnNumber();
5528        double upPseudoCost = obj1->upDynamicPseudoCost();
5529        double saveUp = upPseudoCost;
5530        upPseudoCost = CoinMax(upPseudoCost,smallUp);
5531        upPseudoCost = CoinMax(upPseudoCost,up[iColumn]);
5532        upPseudoCost = CoinMax(upPseudoCost,0.1*down[iColumn]);
5533        obj1->setUpDynamicPseudoCost(upPseudoCost);
5534        if (upPseudoCost>saveUp&&!numberNodes_)
5535          printf("For %d up went from %g to %g\n",
5536                 iColumn,saveUp,upPseudoCost);
5537        double downPseudoCost = obj1->downDynamicPseudoCost();
5538        double saveDown = downPseudoCost;
5539        downPseudoCost = CoinMax(downPseudoCost,smallDown);
5540        downPseudoCost = CoinMax(downPseudoCost,down[iColumn]);
5541        downPseudoCost = CoinMax(downPseudoCost,0.1*down[iColumn]);
5542        obj1->setDownDynamicPseudoCost(downPseudoCost);
5543        if (downPseudoCost>saveDown&&!numberNodes_)
5544          printf("For %d down went from %g to %g\n",
5545                 iColumn,saveDown,downPseudoCost);
5546      }
5547    }
5548  }
5549  delete [] down;
5550  delete [] up;
5551}
5552
5553/*
5554  Set branching priorities.
5555
5556  Setting integer priorities looks pretty robust; the call to findIntegers
5557  makes sure that SimpleInteger objects are in place. Setting priorities for
5558  other objects is entirely dependent on their existence, and the routine may
5559  quietly fail in several directions.
5560*/
5561
5562void 
5563CbcModel::passInPriorities (const int * priorities,
5564                            bool ifObject)
5565{
5566  findIntegers(false);
5567  int i;
5568  if (priorities) {
5569    int i0=0;
5570    int i1=numberObjects_-1;
5571    if (ifObject) {
5572      for (i=numberIntegers_;i<numberObjects_;i++) {
5573        object_[i]->setPriority(priorities[i-numberIntegers_]);
5574      }
5575      i0=numberIntegers_;
5576    } else {
5577      for (i=0;i<numberIntegers_;i++) {
5578        object_[i]->setPriority(priorities[i]);
5579      }
5580      i1=numberIntegers_-1;
5581    }
5582    messageHandler()->message(CBC_PRIORITY,
5583                              messages())
5584                                << i0<<i1<<numberObjects_ << CoinMessageEol ;
5585  }
5586}
5587
5588// Delete all object information
5589void 
5590CbcModel::deleteObjects()
5591{
5592  int i;
5593  for (i=0;i<numberObjects_;i++)
5594    delete object_[i];
5595  delete [] object_;
5596  object_ = NULL;
5597  numberObjects_=0;
5598  findIntegers(true);
5599}
5600
5601/*!
5602  Ensure all attached objects (OsiObjects, heuristics, and cut
5603  generators) point to this model.
5604*/
5605void CbcModel::synchronizeModel()
5606{
5607  int i;
5608  for (i=0;i<numberHeuristics_;i++) 
5609    heuristic_[i]->setModel(this);
5610  for (i=0;i<numberObjects_;i++) {
5611    CbcObject * obj =
5612      dynamic_cast <CbcObject *>(object_[i]) ;
5613    if (obj)
5614      obj->setModel(this);
5615  }
5616  for (i=0;i<numberCutGenerators_;i++)
5617    generator_[i]->refreshModel(this);
5618}
5619
5620// Fill in integers and create objects
5621
5622/**
5623  The routine first does a scan to count the number of integer variables.
5624  It then creates an array, integerVariable_, to store the indices of the
5625  integer variables, and an array of `objects', one for each variable.
5626
5627  The scan is repeated, this time recording the index of each integer
5628  variable in integerVariable_, and creating an CbcSimpleInteger object that
5629  contains information about the integer variable. Initially, this is just
5630  the index and upper & lower bounds.
5631
5632  \todo
5633  Note the assumption in cbc that the first numberIntegers_ objects are
5634  CbcSimpleInteger. In particular, the code which handles the startAgain
5635  case assumes that if the object_ array exists it can simply replace the first
5636  numberInteger_ objects. This is arguably unsafe.
5637
5638  I am going to re-order if necessary
5639*/
5640
5641void 
5642CbcModel::findIntegers(bool startAgain,int type)
5643{
5644  assert(solver_);
5645/*
5646  No need to do this if we have previous information, unless forced to start
5647  over.
5648*/
5649  if (numberIntegers_&&!startAgain&&object_)
5650    return;
5651/*
5652  Clear out the old integer variable list, then count the number of integer
5653  variables.
5654*/
5655  delete [] integerVariable_;
5656  integerVariable_ = NULL;
5657  numberIntegers_=0;
5658  int numberColumns = getNumCols();
5659  int iColumn;
5660  for (iColumn=0;iColumn<numberColumns;iColumn++) {
5661    if (isInteger(iColumn)) 
5662      numberIntegers_++;
5663  }
5664  // Find out how many old non-integer objects there are
5665  int nObjects=0;
5666  OsiObject ** oldObject = object_;
5667  int iObject;
5668  for (iObject = 0;iObject<numberObjects_;iObject++) {
5669    CbcSimpleInteger * obj =
5670      dynamic_cast <CbcSimpleInteger *>(oldObject[iObject]) ;
5671    if (obj) 
5672      delete oldObject[iObject];
5673    else
5674      oldObject[nObjects++]=oldObject[iObject];
5675  }
5676  // See if there any SOS
5677#ifdef COIN_HAS_CLP
5678  if (!nObjects) {
5679    OsiClpSolverInterface * clpSolver
5680      = dynamic_cast<OsiClpSolverInterface *> (solver_);
5681    if (clpSolver&&clpSolver->numberSOS()) {
5682      // deal with sos
5683      const CoinSet * setInfo = clpSolver->setInfo();
5684      int numberSOS = clpSolver->numberSOS();
5685      nObjects=0;
5686      delete [] oldObject;
5687      oldObject = new OsiObject * [numberSOS];
5688      for (int i=0;i<numberSOS;i++) {
5689        int type = setInfo[i].setType();
5690        int n=setInfo[i].numberEntries();
5691        const int * which = setInfo[i].which();
5692        const double * weights = setInfo[i].weights();
5693        oldObject[nObjects++] = new CbcSOS(this,n,which,weights,i,type);
5694      }
5695    }
5696  }
5697#endif
5698   
5699/*
5700  Found any? Allocate an array to hold the indices of the integer variables.
5701  Make a large enough array for all objects
5702*/
5703  delete [] integerVariable_;
5704  object_ = new OsiObject * [numberIntegers_+nObjects];
5705  numberObjects_=numberIntegers_+nObjects;
5706  integerVariable_ = new int [numberIntegers_];
5707/*
5708  Walk the variables again, filling in the indices and creating objects for
5709  the integer variables. Initially, the objects hold the index and upper &
5710  lower bounds.
5711*/
5712  numberIntegers_=0;
5713  for (iColumn=0;iColumn<numberColumns;iColumn++) {
5714    if(isInteger(iColumn)) {
5715      if (!type)
5716        object_[numberIntegers_] =
5717          new CbcSimpleInteger(this,iColumn);
5718      else if (type==1)
5719        object_[numberIntegers_] =
5720          new CbcSimpleIntegerPseudoCost(this,numberIntegers_,iColumn,0.3);
5721      integerVariable_[numberIntegers_++]=iColumn;
5722    }
5723  }
5724  // Now append other objects
5725  memcpy(object_+numberIntegers_,oldObject,nObjects*sizeof(OsiObject *));
5726  // Delete old array (just array)
5727  delete [] oldObject;
5728 
5729  if (!numberObjects_)
5730      handler_->message(CBC_NOINT,messages_) << CoinMessageEol ;
5731}
5732/* If numberBeforeTrust >0 then we are going to use CbcBranchDynamic.
5733   Scan and convert CbcSimpleInteger objects
5734*/
5735void 
5736CbcModel::convertToDynamic()
5737{
5738  int iObject;
5739  const double * cost = solver_->getObjCoefficients();
5740  bool allDynamic=true;
5741  for (iObject = 0;iObject<numberObjects_;iObject++) {
5742    CbcSimpleInteger * obj1 =
5743      dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ;
5744    CbcSimpleIntegerPseudoCost * obj1a =
5745      dynamic_cast <CbcSimpleIntegerPseudoCost *>(object_[iObject]) ;
5746    CbcSimpleIntegerDynamicPseudoCost * obj2 =
5747      dynamic_cast <CbcSimpleIntegerDynamicPseudoCost *>(object_[iObject]) ;
5748    if (obj1&&!obj2) {
5749      // replace
5750      int iColumn = obj1->columnNumber();
5751      int priority = obj1->priority();
5752      int preferredWay = obj1->preferredWay();
5753      delete object_[iObject];
5754      double costValue = CoinMax(1.0e-5,fabs(cost[iColumn]));
5755      // treat as if will cost what it says up
5756      double upCost=costValue;
5757      // and balance at breakeven of 0.3
5758      double downCost=(0.7*upCost)/0.3;
5759      if (obj1a) {
5760        upCost=obj1a->upPseudoCost();
5761        downCost=obj1a->downPseudoCost();
5762      }
5763      CbcSimpleIntegerDynamicPseudoCost * newObject =
5764        new CbcSimpleIntegerDynamicPseudoCost(this,iColumn,1.0e0*downCost,1.0e0*upCost);
5765      newObject->setNumberBeforeTrust(numberBeforeTrust_);
5766      newObject->setPriority(priority);
5767      newObject->setPreferredWay(preferredWay);
5768      object_[iObject] = newObject;
5769    } else if (!obj2) {
5770      allDynamic=false;
5771    }
5772  }
5773  if (branchingMethod_) {
5774    if ((branchingMethod_->whichMethod()&1)==0&&!branchingMethod_->chooseMethod()) {
5775      // Need a method which can do better
5776      delete branchingMethod_;
5777      branchingMethod_=NULL;
5778    }
5779  }
5780  if (!branchingMethod_&&allDynamic) {
5781    // create one
5782    branchingMethod_ = new CbcBranchDynamicDecision();
5783  }
5784}
5785
5786/* Add in any object information (objects are cloned - owner can delete
5787   originals */
5788void 
5789CbcModel::addObjects(int numberObjects, CbcObject ** objects)
5790{
5791 // If integers but not enough objects fudge
5792  if (numberIntegers_>numberObjects_)
5793    findIntegers(true);
5794  /* But if incoming objects inherit from simple integer we just want
5795     to replace */
5796  int numberColumns = solver_->getNumCols();
5797  /** mark is -1 if not integer, >=0 if using existing simple integer and
5798      >=numberColumns if using new integer */
5799  int * mark = new int[numberColumns];
5800  int i;
5801  for (i=0;i<numberColumns;i++)
5802    mark[i]=-1;
5803  int newNumberObjects = numberObjects;
5804  int newIntegers=0;
5805  for (i=0;i<numberObjects;i++) { 
5806    CbcSimpleInteger * obj =
5807      dynamic_cast <CbcSimpleInteger *>(objects[i]) ;
5808    if (obj) {
5809      int iColumn = obj->columnNumber();
5810      mark[iColumn]=i+numberColumns;
5811      newIntegers++;
5812    }
5813  }
5814  // and existing
5815  for (i=0;i<numberObjects_;i++) { 
5816    CbcSimpleInteger * obj =
5817      dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
5818    if (obj) {
5819      int iColumn = obj->columnNumber();
5820      if (mark[iColumn]<0) {
5821        newIntegers++;
5822        newNumberObjects++;
5823        mark[iColumn]=i;
5824      }
5825    }
5826  } 
5827  delete [] integerVariable_;
5828  integerVariable_=NULL;
5829  if (newIntegers!=numberIntegers_) 
5830    printf("changing number of integers from %d to %d\n",
5831           numberIntegers_,newIntegers);
5832  numberIntegers_ = newIntegers;
5833  integerVariable_ = new int [numberIntegers_];
5834  OsiObject ** temp  = new OsiObject * [newNumberObjects];
5835  // Put integers first
5836  newIntegers=0;
5837  numberIntegers_=0;
5838  for (i=0;i<numberColumns;i++) {
5839    int which = mark[i];
5840    if (which>=0) {
5841      if (!isInteger(i)) {
5842        newIntegers++;
5843        solver_->setInteger(i);
5844      }
5845      if (which<numberColumns) {
5846        temp[numberIntegers_]=object_[which];
5847        object_[which]=NULL;
5848      } else {
5849        temp[numberIntegers_]=objects[which-numberColumns]->clone();
5850      }
5851      integerVariable_[numberIntegers_++]=i;
5852    }
5853  }
5854  if (newIntegers)
5855    printf("%d variables were declared integer\n",newIntegers);
5856  int n=numberIntegers_;
5857  // Now rest of old
5858  for (i=0;i<numberObjects_;i++) { 
5859    if (object_[i]) {
5860      CbcSimpleInteger * obj =
5861        dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
5862      if (obj) {
5863        delete object_[i];
5864      } else {
5865        temp[n++]=object_[i];
5866      }
5867    }
5868  }
5869  // and rest of new
5870  for (i=0;i<numberObjects;i++) { 
5871    CbcSimpleInteger * obj =
5872      dynamic_cast <CbcSimpleInteger *>(objects[i]) ;
5873    if (!obj) {
5874      temp[n]=objects[i]->clone();
5875      CbcObject * obj =
5876        dynamic_cast <CbcObject *>(temp[n]) ;
5877      if (obj)
5878        obj->setModel(this);
5879      n++;
5880    }
5881  }
5882  delete [] mark;
5883  delete [] object_;
5884  object_ = temp;
5885  assert (n==newNumberObjects);
5886  numberObjects_ = newNumberObjects;
5887}
5888/* Add in any object information (objects are cloned - owner can delete
5889   originals */
5890void 
5891CbcModel::addObjects(int numberObjects, OsiObject ** objects)
5892{
5893  // If integers but not enough objects fudge
5894  if (numberIntegers_>numberObjects_)
5895    findIntegers(true);
5896  /* But if incoming objects inherit from simple integer we just want
5897     to replace */
5898  int numberColumns = solver_->getNumCols();
5899  /** mark is -1 if not integer, >=0 if using existing simple integer and
5900      >=numberColumns if using new integer */
5901  int * mark = new int[numberColumns];
5902  int i;
5903  for (i=0;i<numberColumns;i++)
5904    mark[i]=-1;
5905  int newNumberObjects = numberObjects;
5906  int newIntegers=0;
5907  for (i=0;i<numberObjects;i++) { 
5908    CbcSimpleInteger * obj =
5909      <