source: branches/devel/Cbc/src/CbcModel.cpp @ 417

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

for allCuts

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