source: trunk/CbcModel.cpp @ 45

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

correct minor bug

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