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

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

unbounded changes from devel

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