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

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

for cut generators

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