source: trunk/CbcModel.cpp @ 66

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

few additional methods

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