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

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

for Kish yet again

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