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

Last change on this file since 414 was 414, checked in by forrest, 14 years ago

mostly try and simplify CbcOr?...

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