source: trunk/Cbc/src/CbcModel.cpp @ 434

Last change on this file since 434 was 434, checked in by forrest, 15 years ago

ampl fix

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