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

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

handler should have had _

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