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

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

many changes

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