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

Last change on this file since 325 was 325, checked in by andreasw, 14 years ago

changed Config.h behavior

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