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

Last change on this file since 391 was 391, checked in by forrest, 14 years ago

for ampl and version

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