source: stable/1.1/Cbc/src/CbcModel.cpp @ 584

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

fix seg fault

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