source: trunk/CbcModel.cpp @ 96

Last change on this file since 96 was 96, checked in by forrest, 15 years ago

mostly fix overflow in priority

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