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

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

bit of message missing

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