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

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

isContinuousUnbounded

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