source: trunk/CbcModel.cpp @ 146

Last change on this file since 146 was 146, checked in by forrest, 16 years ago

for node stats

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