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

Last change on this file since 354 was 354, checked in by ladanyi, 14 years ago

finished switch to subversion

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