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

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

extend fpump heuristic and allow more probing

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