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

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

for infeasible problems

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