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

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

for nonlinear

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