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

Last change on this file since 633 was 633, checked in by forrest, 12 years ago

for cgl infeas cut

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