source: trunk/CbcModel.cpp @ 67

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

modify message

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 169.0 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        handler_->printing(!numberNodes_&&generator_[i]->timing())
2916          <<generator_[i]->timeInCutGenerator();
2917        handler_->message()
2918          <<newFrequency
2919          <<CoinMessageEol ;
2920      }
2921    } 
2922    delete [] count ;
2923    if( !numberNodes_) {
2924      if( !willBeCutsInTree) {
2925        // Take off cuts
2926        cuts = OsiCuts();
2927        numberNewCuts=0;
2928        // update size of problem
2929        numberRowsAtContinuous_ = solver_->getNumRows() ;
2930      } else {
2931#ifdef COIN_USE_CLP
2932        OsiClpSolverInterface * clpSolver
2933          = dynamic_cast<OsiClpSolverInterface *> (solver_);
2934        if (clpSolver) {
2935        // make sure factorization can't carry over
2936          int options = clpSolver->specialOptions();
2937          if (options>=0)
2938            clpSolver->setSpecialOptions(clpSolver->specialOptions()&(~8));
2939        }
2940#endif
2941      }
2942    }
2943  } else if (numberCutGenerators_) {
2944    int i;
2945    // add to counts anyway
2946    for (i = 0;i<numberCutGenerators_;i++) 
2947      generator_[i]->incrementNumberCutsInTotal(countRowCuts[i]);
2948    // What if not feasible as cuts may have helped
2949    if (feasible) {
2950      for (i = 0;i<numberNewCuts;i++) {
2951        int iGenerator = whichGenerator[i];
2952        if (iGenerator>=0)
2953          generator_[iGenerator]->incrementNumberCutsActive();
2954      }
2955    }
2956  }
2957
2958  delete [] countRowCuts ;
2959  delete [] countColumnCuts ;
2960
2961
2962#ifdef CHECK_CUT_COUNTS
2963  if (feasible)
2964  { delete basis_ ;
2965    basis_ = dynamic_cast<CoinWarmStartBasis*>(solver_->getWarmStart()) ;
2966    printf("solveWithCuts: Number of rows at end (only active cuts) %d\n",
2967           numberRowsAtContinuous_+numberNewCuts+numberOldActiveCuts) ;
2968    basis_->print() ; }
2969#endif
2970#ifdef CBC_DEBUG
2971  if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
2972    assert(feasible) ;
2973#endif
2974
2975  return feasible ; }
2976
2977
2978/*
2979  Remove slack cuts. We obtain a basis and scan it. Cuts with basic slacks
2980  are purged. If any cuts are purged, resolve() is called to restore the
2981  solution held in the solver.  If resolve() pivots, there's the possibility
2982  that a slack may be pivoted in (trust me :-), so the process iterates.
2983  Setting allowResolve to false will suppress reoptimisation (but see note
2984  below).
2985
2986  At the level of the solver's constraint system, loose cuts are really
2987  deleted.  There's an implicit assumption that deleteRows will also update
2988  the active basis in the solver.
2989
2990  At the level of nodes and models, it's more complicated.
2991
2992  New cuts exist only in the collection of cuts passed as a parameter. They
2993  are deleted from the collection and that's the end of them.
2994
2995  Older cuts have made it into addedCuts_. Two separate actions are needed.
2996  The reference count for the CbcCountRowCut object is decremented. If this
2997  count falls to 0, the node which owns the cut is located, the reference to
2998  the cut is removed, and then the cut object is destroyed (courtesy of the
2999  CbcCountRowCut destructor). We also need to set the addedCuts_ entry to
3000  NULL. This is important so that when it comes time to generate basis edits
3001  we can tell this cut was dropped from the basis during processing of the
3002  node.
3003
3004  NOTE: In general, it's necessary to call resolve() after purging slack
3005        cuts.  Deleting constraints constitutes a change in the problem, and
3006        an OSI is not required to maintain a valid solution when the problem
3007        is changed. But ... it's really useful to maintain the active basis,
3008        and the OSI is supposed to do that. (Yes, it's splitting hairs.) In
3009        some places, it's possible to know that the solution will never be
3010        consulted after this call, only the basis.  (E.g., this routine is
3011        called as a last act before generating info to place the node in the
3012        live set.) For such use, set allowResolve to false.
3013 
3014  TODO: No real harm would be done if we just ignored the rare occasion when
3015        the call to resolve() pivoted a slack back into the basis. It's a
3016        minor inefficiency, at worst. But it does break assertions which
3017        check that there are no loose cuts in the basis. It might be better
3018        to remove the assertions.
3019*/
3020
3021void
3022CbcModel::takeOffCuts (OsiCuts &newCuts, int *whichGenerator,
3023                       int &numberOldActiveCuts, int &numberNewCuts,
3024                       bool allowResolve)
3025
3026{ // int resolveIterations = 0 ;
3027  int firstOldCut = numberRowsAtContinuous_ ;
3028  int totalNumberCuts = numberNewCuts+numberOldActiveCuts ;
3029  int *solverCutIndices = new int[totalNumberCuts] ;
3030  int *newCutIndices = new int[numberNewCuts] ;
3031  const CoinWarmStartBasis* ws ;
3032  CoinWarmStartBasis::Status status ;
3033  bool needPurge = true ;
3034/*
3035  The outer loop allows repetition of purge in the event that reoptimisation
3036  changes the basis. To start an iteration, clear the deletion counts and grab
3037  the current basis.
3038*/
3039  while (needPurge)
3040  { int numberNewToDelete = 0 ;
3041    int numberOldToDelete = 0 ;
3042    int i ;
3043    ws = dynamic_cast<const CoinWarmStartBasis*>(solver_->getWarmStart()) ;
3044/*
3045  Scan the basis entries of the old cuts generated prior to this round of cut
3046  generation.  Loose cuts are `removed' by decrementing their reference count
3047  and setting the addedCuts_ entry to NULL. (If the reference count falls to
3048  0, they're really deleted.  See CbcModel and CbcCountRowCut doc'n for
3049  principles of cut handling.)
3050*/
3051    int oldCutIndex = 0 ;
3052    for (i = 0 ; i < numberOldActiveCuts ; i++)
3053    { status = ws->getArtifStatus(i+firstOldCut) ;
3054      while (!addedCuts_[oldCutIndex]) oldCutIndex++ ;
3055      assert(oldCutIndex < currentNumberCuts_) ;
3056      if (status == CoinWarmStartBasis::basic)
3057      { solverCutIndices[numberOldToDelete++] = i+firstOldCut ;
3058        if (addedCuts_[oldCutIndex]->decrement() == 0)
3059          delete addedCuts_[oldCutIndex] ;
3060        addedCuts_[oldCutIndex] = NULL ;
3061        oldCutIndex++ ; }
3062      else
3063      { oldCutIndex++ ; } }
3064/*
3065  Scan the basis entries of the new cuts generated with this round of cut
3066  generation.  At this point, newCuts is the only record of the new cuts, so
3067  when we delete loose cuts from newCuts, they're really gone. newCuts is a
3068  vector, so it's most efficient to compress it (eraseRowCut) from back to
3069  front.
3070*/
3071    int firstNewCut = firstOldCut+numberOldActiveCuts ;
3072    int k = 0 ;
3073    for (i = 0 ; i < numberNewCuts ; i++)
3074    { status = ws->getArtifStatus(i+firstNewCut) ;
3075      if (status == CoinWarmStartBasis::basic&&whichGenerator[i]!=-2)
3076      { solverCutIndices[numberNewToDelete+numberOldToDelete] = i+firstNewCut ;
3077        newCutIndices[numberNewToDelete++] = i ; }
3078      else
3079      { // save which generator did it
3080        whichGenerator[k++] = whichGenerator[i] ; } }
3081    delete ws ;
3082    for (i = numberNewToDelete-1 ; i >= 0 ; i--)
3083    { int iCut = newCutIndices[i] ;
3084      newCuts.eraseRowCut(iCut) ; }
3085/*
3086  Did we delete anything? If so, delete the cuts from the constraint system
3087  held in the solver and reoptimise unless we're forbidden to do so. If the
3088  call to resolve() results in pivots, there's the possibility we again have
3089  basic slacks. Repeat the purging loop.
3090*/
3091    if (numberNewToDelete+numberOldToDelete > 0)
3092    { solver_->deleteRows(numberNewToDelete+numberOldToDelete,
3093                          solverCutIndices) ;
3094      numberNewCuts -= numberNewToDelete ;
3095      numberOldActiveCuts -= numberOldToDelete ;
3096#     ifdef CBC_DEBUG
3097      std::cout << "takeOffCuts: purged " << numberOldToDelete << "+"
3098                << numberNewToDelete << " cuts." << std::endl ;
3099#     endif
3100      if (allowResolve)
3101      { 
3102        phase_=3;
3103        solver_->resolve() ;
3104        if (solver_->getIterationCount() == 0)
3105        { needPurge = false ; }
3106#       ifdef CBC_DEBUG
3107        else
3108        { std::cout << "Repeating purging loop. "
3109                    << solver_->getIterationCount() << " iters."
3110                    << std::endl ; }
3111#       endif
3112      }
3113      else
3114      { needPurge = false ; } }
3115    else
3116    { needPurge = false ; } }
3117/*
3118  Clean up and return.
3119*/
3120  delete [] solverCutIndices ;
3121  delete [] newCutIndices ;
3122}
3123
3124
3125
3126bool
3127CbcModel::resolve()
3128{
3129  // We may have deliberately added in violated cuts - check to avoid message
3130  int iRow;
3131  int numberRows = solver_->getNumRows();
3132  const double * rowLower = solver_->getRowLower();
3133  const double * rowUpper = solver_->getRowUpper();
3134  bool feasible=true;
3135  for (iRow= numberRowsAtContinuous_;iRow<numberRows;iRow++) {
3136    if (rowLower[iRow]>rowUpper[iRow]+1.0e-8)
3137      feasible=false;
3138  }
3139/*
3140  Reoptimize. Consider the possibility that we should fathom on bounds. But be
3141  careful --- where the objective takes on integral values, we may want to keep
3142  a solution where the objective is right on the cutoff.
3143*/
3144  if (feasible)
3145  {
3146    solver_->resolve() ;
3147    numberIterations_ += getIterationCount() ;
3148    feasible = (solver_->isProvenOptimal() &&
3149                !solver_->isDualObjectiveLimitReached()) ; }
3150  if (!feasible&& continuousObjective_ <-1.0e30) {
3151    // at root node - double double check
3152    bool saveTakeHint;
3153    OsiHintStrength saveStrength;
3154    solver_->getHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength);
3155    if (saveTakeHint||saveStrength==OsiHintIgnore) {
3156      solver_->setHintParam(OsiDoDualInResolve,false,OsiHintDo) ;
3157      solver_->resolve();
3158      solver_->setHintParam(OsiDoDualInResolve,saveTakeHint,saveStrength);
3159      numberIterations_ += getIterationCount() ;
3160      feasible = solver_->isProvenOptimal();
3161    }
3162  }
3163  return feasible ; }
3164
3165
3166/* Set up objects.  Only do ones whose length is in range.
3167   If makeEquality true then a new model may be returned if
3168   modifications had to be made, otherwise "this" is returned.
3169
3170   Could use Probing at continuous to extend objects
3171*/
3172CbcModel * 
3173CbcModel::findCliques(bool makeEquality,
3174                      int atLeastThisMany, int lessThanThis, int defaultValue)
3175{
3176  // No objects are allowed to exist
3177  assert(numberObjects_==numberIntegers_||!numberObjects_);
3178  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3179  int numberRows = solver_->getNumRows();
3180  int numberColumns = solver_->getNumCols();
3181
3182  // We may want to add columns
3183  int numberSlacks=0;
3184  int * rows = new int[numberRows];
3185  double * element =new double[numberRows];
3186
3187  int iRow;
3188
3189  findIntegers(true);
3190  numberObjects_=numberIntegers_;
3191
3192  int numberCliques=0;
3193  CbcObject ** object = new CbcObject * [numberRows];
3194  int * which = new int[numberIntegers_];
3195  char * type = new char[numberIntegers_];
3196  int * lookup = new int[numberColumns];
3197  int i;
3198  for (i=0;i<numberColumns;i++) 
3199    lookup[i]=-1;
3200  for (i=0;i<numberIntegers_;i++) 
3201    lookup[integerVariable_[i]]=i;
3202
3203  // Statistics
3204  int totalP1=0,totalM1=0;
3205  int numberBig=0,totalBig=0;
3206  int numberFixed=0;
3207
3208  // Row copy
3209  const double * elementByRow = matrixByRow.getElements();
3210  const int * column = matrixByRow.getIndices();
3211  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3212  const int * rowLength = matrixByRow.getVectorLengths();
3213
3214  // Column lengths for slacks
3215  const int * columnLength = solver_->getMatrixByCol()->getVectorLengths();
3216
3217  const double * lower = getColLower();
3218  const double * upper = getColUpper();
3219  const double * rowLower = getRowLower();
3220  const double * rowUpper = getRowUpper();
3221
3222  for (iRow=0;iRow<numberRows;iRow++) {
3223    int numberP1=0, numberM1=0;
3224    int j;
3225    double upperValue=rowUpper[iRow];
3226    double lowerValue=rowLower[iRow];
3227    bool good=true;
3228    int slack = -1;
3229    for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3230      int iColumn = column[j];
3231      int iInteger=lookup[iColumn];
3232      if (upper[iColumn]-lower[iColumn]<1.0e-8) {
3233        // fixed
3234        upperValue -= lower[iColumn]*elementByRow[j];
3235        lowerValue -= lower[iColumn]*elementByRow[j];
3236        continue;
3237      } else if (upper[iColumn]!=1.0||lower[iColumn]!=0.0) {
3238        good = false;
3239        break;
3240      } else if (iInteger<0) {
3241        good = false;
3242        break;
3243      } else {
3244        if (columnLength[iColumn]==1)
3245          slack = iInteger;
3246      }
3247      if (fabs(elementByRow[j])!=1.0) {
3248        good=false;
3249        break;
3250      } else if (elementByRow[j]>0.0) {
3251        which[numberP1++]=iInteger;
3252      } else {
3253        numberM1++;
3254        which[numberIntegers_-numberM1]=iInteger;
3255      }
3256    }
3257    int iUpper = (int) floor(upperValue+1.0e-5);
3258    int iLower = (int) ceil(lowerValue-1.0e-5);
3259    int state=0;
3260    if (upperValue<1.0e6) {
3261      if (iUpper==1-numberM1)
3262        state=1;
3263      else if (iUpper==-numberM1)
3264        state=2;
3265      else if (iUpper<-numberM1)
3266        state=3;
3267    }
3268    if (!state&&lowerValue>-1.0e6) {
3269      if (-iLower==1-numberP1)
3270        state=-1;
3271      else if (-iLower==-numberP1)
3272        state=-2;
3273      else if (-iLower<-numberP1)
3274        state=-3;
3275    }
3276    if (good&&state) {
3277      if (abs(state)==3) {
3278        // infeasible
3279        numberObjects_=-1;
3280        break;
3281      } else if (abs(state)==2) {
3282        // we can fix all
3283        numberFixed += numberP1+numberM1;
3284        if (state>0) {
3285          // fix all +1 at 0, -1 at 1
3286          for (i=0;i<numberP1;i++)
3287            solver_->setColUpper(integerVariable_[which[i]],0.0);
3288          for (i=0;i<numberM1;i++)
3289            solver_->setColLower(integerVariable_[which[numberIntegers_-i-1]],
3290                                 1.0);
3291        } else {
3292          // fix all +1 at 1, -1 at 0
3293          for (i=0;i<numberP1;i++)
3294            solver_->setColLower(integerVariable_[which[i]],1.0);
3295          for (i=0;i<numberM1;i++)
3296            solver_->setColUpper(integerVariable_[which[numberIntegers_-i-1]],
3297                                 0.0);
3298        }
3299      } else {
3300        int length = numberP1+numberM1;
3301        if (length >= atLeastThisMany&&length<lessThanThis) {
3302          // create object
3303          bool addOne=false;
3304          int objectType;
3305          if (iLower==iUpper) {
3306            objectType=1;
3307          } else {
3308            if (makeEquality) {
3309              objectType=1;
3310              element[numberSlacks]=state;
3311              rows[numberSlacks++]=iRow;
3312              addOne=true;
3313            } else {
3314              objectType=0;
3315            }
3316          }
3317          if (state>0) {
3318            totalP1 += numberP1;
3319            totalM1 += numberM1;
3320            for (i=0;i<numberP1;i++)
3321              type[i]=1;
3322            for (i=0;i<numberM1;i++) {
3323              which[numberP1]=which[numberIntegers_-i-1];
3324              type[numberP1++]=0;
3325            }
3326          } else {
3327            totalP1 += numberM1;
3328            totalM1 += numberP1;
3329            for (i=0;i<numberP1;i++)
3330              type[i]=0;
3331            for (i=0;i<numberM1;i++) {
3332              which[numberP1]=which[numberIntegers_-i-1];
3333              type[numberP1++]=1;
3334            }
3335          }
3336          if (addOne) {
3337            // add in slack
3338            which[numberP1]=numberIntegers_+numberSlacks-1;
3339            slack = numberP1;
3340            type[numberP1++]=1;
3341          } else if (slack >= 0) {
3342            for (i=0;i<numberP1;i++) {
3343              if (which[i]==slack) {
3344                slack=i;
3345              }
3346            }
3347          }
3348          object[numberCliques++] = new CbcClique(this,objectType,numberP1,
3349                                              which,type,
3350                                               1000000+numberCliques,slack);
3351        } else if (numberP1+numberM1 >= lessThanThis) {
3352          // too big
3353          numberBig++;
3354          totalBig += numberP1+numberM1;
3355        }
3356      }
3357    }
3358  }
3359  delete [] which;
3360  delete [] type;
3361  delete [] lookup;
3362  if (numberCliques<0) {
3363    printf("*** Problem infeasible\n");
3364  } else {
3365    if (numberCliques)
3366      printf("%d cliques of average size %g found, %d P1, %d M1\n",
3367             numberCliques,
3368             ((double)(totalP1+totalM1))/((double) numberCliques),
3369             totalP1,totalM1);
3370    else
3371      printf("No cliques found\n");
3372    if (numberBig)
3373      printf("%d large cliques ( >= %d) found, total %d\n",
3374             numberBig,lessThanThis,totalBig);
3375    if (numberFixed)
3376      printf("%d variables fixed\n",numberFixed);
3377  }
3378  if (numberCliques>0&&numberSlacks&&makeEquality) {
3379    printf("adding %d integer slacks\n",numberSlacks);
3380    // add variables to make equality rows
3381    int * temp = new int[numberIntegers_+numberSlacks];
3382    memcpy(temp,integerVariable_,numberIntegers_*sizeof(int));
3383    // Get new model
3384    CbcModel * newModel = new CbcModel(*this);
3385    OsiSolverInterface * newSolver = newModel->solver();
3386    for (i=0;i<numberSlacks;i++) {
3387      temp[i+numberIntegers_]=i+numberColumns;
3388      int iRow = rows[i];
3389      double value = element[i];
3390      double lowerValue = 0.0;
3391      double upperValue = 1.0;
3392      double objValue  = 0.0;
3393      CoinPackedVector column(1,&iRow,&value);
3394      newSolver->addCol(column,lowerValue,upperValue,objValue);
3395      // set integer
3396      newSolver->setInteger(numberColumns+i);
3397      if (value >0)
3398        newSolver->setRowLower(iRow,rowUpper[iRow]);
3399      else
3400        newSolver->setRowUpper(iRow,rowLower[iRow]);
3401    }
3402    // replace list of integers
3403    for (i=0;i<newModel->numberObjects_;i++)
3404      delete newModel->object_[i];
3405    newModel->numberObjects_ = 0;
3406    delete [] newModel->object_;
3407    newModel->object_=NULL;
3408    newModel->findIntegers(true); //Set up all integer objects
3409    if (priority_) {
3410      // old model had priorities
3411      delete [] newModel->priority_;
3412      newModel->priority_ = new int[newModel->numberIntegers_+numberCliques];
3413      memcpy(newModel->priority_,priority_,numberIntegers_*sizeof(int));
3414      for (i=numberIntegers_;i<newModel->numberIntegers_+numberCliques;i++)
3415        newModel->priority_[i]=defaultValue;
3416    }
3417    if (originalColumns_) {
3418      // old model had originalColumns
3419      delete [] newModel->originalColumns_;
3420      newModel->originalColumns_ = new int[numberColumns+numberSlacks];
3421      memcpy(newModel->originalColumns_,originalColumns_,numberColumns*sizeof(int));
3422      // mark as not in previous model
3423      for (i=numberColumns;i<numberColumns+numberSlacks;i++)
3424        newModel->originalColumns_[i]=-1;
3425    }
3426    delete [] rows;
3427    delete [] element;
3428    newModel->addObjects(numberCliques,object);
3429    for (;i<numberCliques;i++) 
3430      delete object[i];
3431    delete [] object;
3432    newModel->synchronizeModel();
3433    return newModel;
3434  } else {
3435    if (numberCliques>0) {
3436      addObjects(numberCliques,object);
3437      for (;i<numberCliques;i++) 
3438        delete object[i];
3439      synchronizeModel();
3440    }
3441    delete [] object;
3442    if (priority_) {
3443      // model had priorities
3444      int * temp = new int[numberIntegers_+numberCliques];
3445      memcpy(temp,priority_,numberIntegers_*sizeof(int));
3446      delete [] priority_;
3447      priority_=temp;
3448      for (i=numberIntegers_;i<numberIntegers_+numberCliques;i++)
3449        priority_[i]=defaultValue;
3450    }
3451    delete [] rows;
3452    delete [] element;
3453    return this;
3454  }
3455}
3456
3457/*
3458  Set branching priorities.
3459
3460  Setting integer priorities looks pretty robust; the call to findIntegers
3461  makes sure that SimpleInteger objects are in place. Setting priorities for
3462  other objects is entirely dependent on their existence, and the routine may
3463  quietly fail in several directions.
3464*/
3465
3466void 
3467CbcModel::passInPriorities (const int * priorities,
3468                            bool ifObject, int defaultValue)
3469{
3470  findIntegers(false);
3471  int i;
3472  if (!priority_) {
3473    priority_ = new int[numberObjects_];
3474    for (i=0;i<numberObjects_;i++)
3475      priority_[i]=defaultValue;
3476  }
3477  if (priorities) {
3478    if (ifObject)
3479      memcpy(priority_+numberIntegers_,priorities,
3480             (numberObjects_-numberIntegers_)*sizeof(int));
3481    else
3482      memcpy(priority_,priorities,numberIntegers_*sizeof(int));
3483  }
3484}
3485
3486// Delete all object information
3487void 
3488CbcModel::deleteObjects()
3489{
3490  delete [] priority_;
3491  priority_=NULL;
3492  int i;
3493  for (i=0;i<numberObjects_;i++)
3494    delete object_[i];
3495  delete [] object_;
3496  object_ = NULL;
3497  numberObjects_=0;
3498  findIntegers(true);
3499}
3500
3501/*!
3502  Ensure all attached objects (CbcObjects, heuristics, and cut
3503  generators) point to this model.
3504*/
3505void CbcModel::synchronizeModel()
3506{
3507  int i;
3508  for (i=0;i<numberHeuristics_;i++) 
3509    heuristic_[i]->setModel(this);
3510  for (i=0;i<numberObjects_;i++)
3511    object_[i]->setModel(this);
3512  for (i=0;i<numberCutGenerators_;i++)
3513    generator_[i]->refreshModel(this);
3514}
3515
3516// Fill in integers and create objects
3517
3518/**
3519  The routine first does a scan to count the number of integer variables.
3520  It then creates an array, integerVariable_, to store the indices of the
3521  integer variables, and an array of `objects', one for each variable.
3522
3523  The scan is repeated, this time recording the index of each integer
3524  variable in integerVariable_, and creating an CbcSimpleInteger object that
3525  contains information about the integer variable. Initially, this is just
3526  the index and upper & lower bounds.
3527
3528  \todo
3529  Note the assumption in cbc that the first numberIntegers_ objects are
3530  CbcSimpleInteger. In particular, the code which handles the startAgain
3531  case assumes that if the object_ array exists it can simply replace the first
3532  numberInteger_ objects. This is arguably unsafe.
3533
3534  I am going to re-order if necessary
3535*/
3536
3537void 
3538CbcModel::findIntegers(bool startAgain)
3539{
3540  assert(solver_);
3541/*
3542  No need to do this if we have previous information, unless forced to start
3543  over.
3544*/
3545  if (numberIntegers_&&!startAgain&&object_)
3546    return;
3547/*
3548  Clear out the old integer variable list, then count the number of integer
3549  variables.
3550*/
3551  delete [] integerVariable_;
3552  numberIntegers_=0;
3553  int numberColumns = getNumCols();
3554  int iColumn;
3555  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3556    if (isInteger(iColumn)) 
3557      numberIntegers_++;
3558  }
3559  // Find out how many old non-integer objects there are
3560  int nObjects=0;
3561  CbcObject ** oldObject = object_;
3562  int iObject;
3563  for (iObject = 0;iObject<numberObjects_;iObject++) {
3564    CbcSimpleInteger * obj =
3565      dynamic_cast <CbcSimpleInteger *>(oldObject[iObject]) ;
3566    if (obj) 
3567      delete oldObject[iObject];
3568    else
3569      oldObject[nObjects++]=oldObject[iObject];
3570  }
3571   
3572/*
3573  Found any? Allocate an array to hold the indices of the integer variables.
3574  Make a large enough array for all objects
3575*/
3576  object_ = new CbcObject * [numberIntegers_+nObjects];
3577  numberObjects_=numberIntegers_+nObjects;;
3578  integerVariable_ = new int [numberIntegers_];
3579/*
3580  Walk the variables again, filling in the indices and creating objects for
3581  the integer variables. Initially, the objects hold the index and upper &
3582  lower bounds.
3583*/
3584  numberIntegers_=0;
3585  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3586    if(isInteger(iColumn)) {
3587      object_[numberIntegers_] =
3588        new CbcSimpleInteger(this,numberIntegers_,iColumn);
3589      integerVariable_[numberIntegers_++]=iColumn;
3590    }
3591  }
3592  // Now append other objects
3593  memcpy(object_+numberIntegers_,oldObject,nObjects*sizeof(CbcObject *));
3594  // Delete old array (just array)
3595  delete [] oldObject;
3596 
3597  if (!numberObjects_)
3598      handler_->message(CBC_NOINT,messages_) << CoinMessageEol ;
3599}
3600
3601/* Add in any object information (objects are cloned - owner can delete
3602   originals */
3603void 
3604CbcModel::addObjects(int numberObjects, CbcObject ** objects)
3605{
3606  int newNumberObjects= numberObjects+numberObjects_;
3607  CbcObject ** temp  = new CbcObject * [newNumberObjects];
3608  int i;
3609  for (i=0;i<numberObjects_;i++) 
3610    temp[i]=object_[i];
3611  for (;i<newNumberObjects;i++) { 
3612    temp[i]=objects[i-numberObjects_]->clone();
3613    temp[i]->setModel(this);
3614  }
3615  delete [] object_;
3616  object_ = temp;
3617  numberObjects_ = newNumberObjects;
3618}
3619
3620/**
3621  This routine sets the objective cutoff value used for fathoming and
3622  determining monotonic variables.
3623
3624  If the fathoming discipline is strict, a small tolerance is added to the
3625  new cutoff. This avoids problems due to roundoff when the target value
3626  is exact. The common example would be an IP with only integer variables in
3627  the objective. If the target is set to the exact value z of the optimum,
3628  it's possible to end up fathoming an ancestor of the solution because the
3629  solver returns z+epsilon.
3630
3631  Determining if strict fathoming is needed is best done by analysis.
3632  In cbc, that's analyseObjective. The default is false.
3633
3634  In cbc we always minimize so add epsilon
3635*/
3636
3637void CbcModel::setCutoff (double value)
3638
3639{
3640#if 0
3641  double tol = 0 ;
3642  int fathomStrict = getIntParam(CbcFathomDiscipline) ;
3643  if (fathomStrict == 1)
3644  { solver_->getDblParam(OsiDualTolerance,tol) ;
3645  tol = tol*(1+fabs(value)) ;
3646 
3647  value += tol ; }
3648#endif
3649  // Solvers know about direction
3650  double direction = solver_->getObjSense();
3651  solver_->setDblParam(OsiDualObjectiveLimit,value*direction); }
3652
3653
3654
3655/*
3656  Call this to really test if a valid solution can be feasible. The cutoff is
3657  passed in as a parameter so that we don't need to worry here after swapping
3658  solvers.  The solution is assumed to be numberColumns in size.  If
3659  fixVariables is true then the bounds of the continuous solver are updated.
3660  The routine returns the objective value determined by reoptimizing from
3661  scratch. If the solution is rejected, this will be worse than the cutoff.
3662
3663  TODO: There's an issue with getting the correct cutoff value: We update the
3664        cutoff in the regular solver, but not in continuousSolver_. But our only
3665        use for continuousSolver_ is verifying candidate solutions. Would it
3666        make sense to update the cutoff? Then we wouldn't need to step around
3667        isDualObjectiveLimitReached().
3668*/
3669double 
3670CbcModel::checkSolution (double cutoff, const double *solution,
3671                         bool fixVariables)
3672
3673{ int numberColumns = solver_->getNumCols();
3674
3675  /*
3676    Grab the continuous solver (the pristine copy of the problem, made before
3677    starting to work on the root node). Save the bounds on the variables.
3678    Install the solution passed as a parameter, and copy it to the model's
3679    currentSolution_.
3680   
3681    TODO: This is a belt-and-suspenders approach. Once the code has settled
3682          a bit, we can cast a critical eye here.
3683  */
3684  OsiSolverInterface * saveSolver = solver_;
3685  if (continuousSolver_)
3686    solver_ = continuousSolver_;
3687  // move solution to continuous copy
3688  solver_->setColSolution(solution);
3689  // Put current solution in safe place
3690  memcpy(currentSolution_,solver_->getColSolution(),
3691         numberColumns*sizeof(double));
3692  //solver_->messageHandler()->setLogLevel(4);
3693
3694  // save original bounds
3695  double * saveUpper = new double[numberColumns];
3696  double * saveLower = new double[numberColumns];
3697  memcpy(saveUpper,getColUpper(),numberColumns*sizeof(double));
3698  memcpy(saveLower,getColLower(),numberColumns*sizeof(double));
3699
3700  /*
3701    Run through the objects and use feasibleRegion() to set variable bounds
3702    so as to fix the variables specified in the objects at their value in this
3703    solution. Since the object list contains (at least) one object for every
3704    integer variable, this has the effect of fixing all integer variables.
3705  */
3706  int i;
3707  for (i=0;i<numberObjects_;i++)
3708    object_[i]->feasibleRegion();
3709  /*
3710    Remove any existing warm start information to be sure there is no
3711    residual influence on initialSolve().
3712  */
3713  CoinWarmStartBasis *slack =
3714      dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
3715  solver_->setWarmStart(slack);
3716  delete slack ;
3717  // Give a hint not to do scaling
3718  //bool saveTakeHint;
3719  //OsiHintStrength saveStrength;
3720  //bool gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
3721  //assert (gotHint);
3722  //solver_->setHintParam(OsiDoScale,false,OsiHintTry);
3723  solver_->initialSolve();
3724  //solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
3725  if (!solver_->isProvenOptimal())
3726    { std::cout << "checkSolution infeas! Retrying wihout scaling."
3727              << std::endl ;
3728    bool saveTakeHint;
3729    OsiHintStrength saveStrength;
3730    bool savePrintHint;
3731    solver_->writeMps("infeas");
3732    bool gotHint = (solver_->getHintParam(OsiDoReducePrint,savePrintHint,saveStrength));
3733    gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
3734    solver_->setHintParam(OsiDoScale,false,OsiHintTry);
3735    solver_->setHintParam(OsiDoReducePrint,false,OsiHintTry) ;
3736    solver_->initialSolve();
3737    solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
3738    solver_->setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ;
3739    }
3740  assert(solver_->isProvenOptimal());
3741  double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
3742
3743  /*
3744    Check that the solution still beats the objective cutoff.
3745
3746    If it passes, make a copy of the primal variable values and do some
3747    cleanup and checks:
3748    + Values of all variables are are within original bounds and values of
3749      all integer variables are within tolerance of integral.
3750    + There are no constraint violations.
3751    There really should be no need for the check against original bounds.
3752    Perhaps an opportunity for a sanity check?
3753  */
3754  if (solver_->isProvenOptimal() && objectiveValue <= cutoff)
3755  { solution = solver_->getColSolution() ;
3756    memcpy(currentSolution_,solution,numberColumns*sizeof(double)) ;
3757
3758    const double * rowLower = solver_->getRowLower() ;
3759    const double * rowUpper = solver_->getRowUpper() ;
3760    int numberRows = solver_->getNumRows() ;
3761    double *rowActivity = new double[numberRows] ;
3762    memset(rowActivity,0,numberRows*sizeof(double)) ;
3763
3764    double integerTolerance = getIntegerTolerance() ;
3765
3766    int iColumn;
3767    for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
3768    { double value = currentSolution_[iColumn] ;
3769      value = CoinMax(value, saveLower[iColumn]) ;
3770      value = CoinMin(value, saveUpper[iColumn]) ;
3771      if (solver_->isInteger(iColumn)) 
3772        assert(fabs(value-currentSolution_[iColumn]) <= integerTolerance) ;
3773      currentSolution_[iColumn] = value ; }
3774   
3775    solver_->getMatrixByCol()->times(currentSolution_,rowActivity) ;
3776    double primalTolerance ;
3777    solver_->getDblParam(OsiPrimalTolerance,primalTolerance) ;
3778    double largestInfeasibility =0.0;
3779    for (i=0 ; i < numberRows ; i++) {
3780      largestInfeasibility = CoinMax(largestInfeasibility,
3781                                 rowLower[i]-rowActivity[i]);
3782      largestInfeasibility = CoinMax(largestInfeasibility,
3783                                 rowActivity[i]-rowUpper[i]);
3784    }
3785    if (largestInfeasibility>100.0*primalTolerance)
3786      handler_->message(CBC_NOTFEAS3, messages_)
3787        << largestInfeasibility << CoinMessageEol ;
3788
3789    delete [] rowActivity ; }
3790  else
3791  { objectiveValue=1.0e50 ; }
3792  /*
3793    Regardless of what we think of the solution, we may need to restore the
3794    original bounds of the continuous solver. Unfortunately, const'ness
3795    prevents us from simply reversing the memcpy used to make these snapshots.
3796  */
3797  if (!fixVariables)
3798  { for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++)
3799    { solver_->setColLower(iColumn,saveLower[iColumn]) ;
3800      solver_->setColUpper(iColumn,saveUpper[iColumn]) ; } }
3801  delete [] saveLower;
3802  delete [] saveUpper;
3803   
3804  /*
3805    Restore the usual solver.
3806  */
3807  solver_=saveSolver;
3808  return objectiveValue;
3809}
3810
3811/*
3812  Call this routine from anywhere when a solution is found. The solution
3813  vector is assumed to contain one value for each structural variable.
3814
3815  The first action is to run checkSolution() to confirm the objective and
3816  feasibility. If this check causes the solution to be rejected, we're done.
3817  If fixVariables = true, the variable bounds held by the continuous solver
3818  will be left fixed to the values in the solution; otherwise they are
3819  restored to the original values.
3820
3821  If the solution is accepted, install it as the best solution.
3822
3823  The routine also contains a hook to run any cut generators that are flagged
3824  to run when a new solution is discovered. There's a potential hazard because
3825  the cut generators see the continuous solver >after< possible restoration of
3826  original bounds (which may well invalidate the solution).
3827*/
3828
3829void
3830CbcModel::setBestSolution (CBC_Message how,
3831                           double & objectiveValue, const double *solution,
3832                           bool fixVariables)
3833
3834{ double cutoff = getCutoff() ;
3835
3836/*
3837  Double check the solution to catch pretenders.
3838*/
3839  objectiveValue = checkSolution(cutoff,solution,fixVariables);
3840  if (objectiveValue > cutoff)
3841  { if (objectiveValue>1.0e30)
3842      handler_->message(CBC_NOTFEAS1, messages_) << CoinMessageEol ;
3843    else
3844      handler_->message(CBC_NOTFEAS2, messages_)
3845        << objectiveValue << cutoff << CoinMessageEol ; }
3846/*
3847  We have a winner. Install it as the new incumbent.
3848  Bump the objective cutoff value and solution counts. Give the user the
3849  good news.
3850*/
3851  else
3852  { bestObjective_ = objectiveValue;
3853    int numberColumns = solver_->getNumCols();
3854    if (!bestSolution_)
3855      bestSolution_ = new double[numberColumns];
3856    memcpy(bestSolution_,currentSolution_,numberColumns*sizeof(double));
3857
3858    cutoff = bestObjective_-dblParam_[CbcCutoffIncrement];
3859    // This is not correct - that way cutoff can go up if maximization
3860    //double direction = solver_->getObjSense();
3861    //setCutoff(cutoff*direction);
3862    setCutoff(cutoff);
3863
3864    if (how==CBC_ROUNDING)
3865      numberHeuristicSolutions_++;
3866    numberSolutions_++;
3867
3868    handler_->message(how,messages_)
3869      <<bestObjective_<<numberIterations_
3870      <<numberNodes_
3871      <<CoinMessageEol;
3872/*
3873  Now step through the cut generators and see if any of them are flagged to
3874  run when a new solution is discovered. Only global cuts are useful. (The
3875  solution being evaluated may not correspond to the current location in the
3876  search tree --- discovered by heuristic, for example.)
3877*/
3878    OsiCuts theseCuts;
3879    int i;
3880    int lastNumberCuts=0;
3881    for (i=0;i<numberCutGenerators_;i++) {
3882      if (generator_[i]->atSolution()) {
3883        generator_[i]->generateCuts(theseCuts,true,NULL);
3884        int numberCuts = theseCuts.sizeRowCuts();
3885        for (int j=lastNumberCuts;j<numberCuts;j++) {
3886          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j);
3887          if (thisCut->globallyValid()) {
3888            if ((specialOptions_&1)!=0) {
3889              /* As these are global cuts -
3890                 a) Always get debugger object
3891                 b) Not fatal error to cutoff optimal (if we have just got optimal)
3892              */
3893              const OsiRowCutDebugger *debugger = solver_->getRowCutDebuggerAlways() ;
3894              if (debugger) {
3895                if(debugger->invalidCut(*thisCut))
3896                  printf("ZZZZ Global cut - cuts off optimal solution!\n");
3897              }
3898            }
3899            // add to global list
3900            globalCuts_.insert(*thisCut);
3901            generator_[i]->incrementNumberCutsInTotal();
3902          }
3903        }
3904      }
3905    }
3906    int numberCuts = theseCuts.sizeColCuts();
3907    for (i=0;i<numberCuts;i++) {
3908      const OsiColCut * thisCut = theseCuts.colCutPtr(i);
3909      if (thisCut->globallyValid()) {
3910        // add to global list
3911        globalCuts_.insert(*thisCut);
3912      }
3913    }
3914  }
3915
3916  return ; }
3917
3918
3919/* Test the current solution for feasibility.
3920
3921   Calculate the number of standard integer infeasibilities, then scan the
3922   remaining objects to see if any of them report infeasibilities.
3923
3924   Currently (2003.08) the only object besides SimpleInteger is Clique, hence
3925   the comments about `odd ones' infeasibilities.
3926*/
3927bool 
3928CbcModel::feasibleSolution(int & numberIntegerInfeasibilities,
3929                        int & numberObjectInfeasibilities) const
3930{
3931  int numberUnsatisfied=0;
3932  double sumUnsatisfied=0.0;
3933  int preferredWay;
3934  double integerTolerance = getIntegerTolerance();
3935  int j;
3936  // Put current solution in safe place
3937  memcpy(currentSolution_,solver_->getColSolution(),
3938         solver_->getNumCols()*sizeof(double));
3939  for (j=0;j<numberIntegers_;j++) {
3940    const CbcObject * object = object_[j];
3941    double infeasibility = object->infeasibility(preferredWay);
3942    if (infeasibility>integerTolerance) {
3943      numberUnsatisfied++;
3944      sumUnsatisfied += infeasibility;
3945    }
3946  }
3947  numberIntegerInfeasibilities = numberUnsatisfied;
3948  for (;j<numberObjects_;j++) {
3949    const CbcObject * object = object_[j];
3950    double infeasibility = object->infeasibility(preferredWay);
3951    if (infeasibility>integerTolerance) {
3952      numberUnsatisfied++;
3953      sumUnsatisfied += infeasibility;
3954    }
3955  }
3956  numberObjectInfeasibilities = numberUnsatisfied-numberIntegerInfeasibilities;
3957  return (!numberUnsatisfied);
3958}
3959
3960/* For all vubs see if we can tighten bounds by solving Lp's
3961   type - 0 just vubs
3962   1 all (could be very slow)
3963   -1 just vubs where variable away from bound
3964   Returns false if not feasible
3965*/
3966bool 
3967CbcModel::tightenVubs(int type, bool allowMultipleBinary, double useCutoff)
3968{
3969
3970  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
3971  int numberRows = solver_->getNumRows();
3972  int numberColumns = solver_->getNumCols();
3973
3974  int iRow,iColumn;
3975
3976  // Row copy
3977  //const double * elementByRow = matrixByRow.getElements();
3978  const int * column = matrixByRow.getIndices();
3979  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
3980  const int * rowLength = matrixByRow.getVectorLengths();
3981
3982  const double * colUpper = solver_->getColUpper();
3983  const double * colLower = solver_->getColLower();
3984  //const double * rowUpper = solver_->getRowUpper();
3985  //const double * rowLower = solver_->getRowLower();
3986
3987  const double * objective = solver_->getObjCoefficients();
3988  //double direction = solver_->getObjSense();
3989  const double * colsol = solver_->getColSolution();
3990
3991  int numberVub=0;
3992  int * continuous = new int[numberColumns];
3993  if (type >= 0) {
3994    double * sort = new double[numberColumns];
3995    for (iRow=0;iRow<numberRows;iRow++) {
3996      int j;
3997      int numberBinary=0;
3998      int numberUnsatisfiedBinary=0;
3999      int numberContinuous=0;
4000      int iCont=-1;
4001      double weight=1.0e30;
4002      for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
4003        int iColumn = column[j];
4004        if (colUpper[iColumn]-colLower[iColumn]>1.0e-8) {
4005          if (solver_->isFreeBinary(iColumn)) {
4006            numberBinary++;
4007            /* For sort I make naive assumption:
4008               x - a * delta <=0 or
4009               -x + a * delta >= 0
4010            */
4011            if (colsol[iColumn]>colLower[iColumn]+1.0e-6&&
4012                colsol[iColumn]<colUpper[iColumn]-1.0e-6) {
4013              numberUnsatisfiedBinary++;
4014              weight = CoinMin(weight,fabs(objective[iColumn]));
4015            }
4016          } else {
4017            numberContinuous++;
4018            iCont=iColumn;
4019          }
4020        }
4021      }
4022      if (numberContinuous==1&&numberBinary) {
4023        if (numberBinary==1||allowMultipleBinary) {
4024          // treat as vub
4025          if (!numberUnsatisfiedBinary)
4026            weight=-1.0; // at end
4027          sort[numberVub]=-weight;
4028          continuous[numberVub++] = iCont;
4029        }
4030      }
4031    }
4032    if (type>0) {
4033      // take so many
4034      CoinSort_2(sort,sort+numberVub,continuous);
4035      numberVub = CoinMin(numberVub,type);
4036    }
4037    delete [] sort;
4038  } else {
4039    for (iColumn=0;iColumn<numberColumns;iColumn++) 
4040      continuous[iColumn]=iColumn;
4041    numberVub=numberColumns;
4042  }
4043  bool feasible = tightenVubs(numberVub,continuous,useCutoff);
4044  delete [] continuous;
4045
4046  return feasible;
4047}
4048// This version is just handed a list of variables
4049bool 
4050CbcModel::tightenVubs(int numberSolves, const int * which,
4051                      double useCutoff)
4052{
4053
4054  int numberColumns = solver_->getNumCols();
4055
4056  int iColumn;
4057 
4058  OsiSolverInterface * solver = solver_;
4059  double saveCutoff = getCutoff() ;
4060 
4061  double * objective = new double[numberColumns];
4062  memcpy(objective,solver_->getObjCoefficients(),numberColumns*sizeof(double));
4063  double direction = solver_->getObjSense();
4064 
4065  // add in objective if there is a cutoff
4066  if (useCutoff<1.0e30) {
4067    // get new version of model
4068    solver = solver_->clone();
4069    CoinPackedVector newRow;
4070    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4071      solver->setObjCoeff(iColumn,0.0); // zero out in new model
4072      if (objective[iColumn]) 
4073        newRow.insert(iColumn,direction * objective[iColumn]);
4074     
4075    }
4076    solver->addRow(newRow,-COIN_DBL_MAX,useCutoff);
4077    // signal no objective
4078    delete [] objective;
4079    objective=NULL;
4080  }
4081  setCutoff(COIN_DBL_MAX);
4082
4083
4084  bool * vub = new bool [numberColumns];
4085  int iVub;
4086
4087  // mark vub columns
4088  for (iColumn=0;iColumn<numberColumns;iColumn++) 
4089    vub[iColumn]=false;
4090  for (iVub=0;iVub<numberSolves;iVub++) 
4091    vub[which[iVub]]=true;
4092  OsiCuts cuts;
4093  // First tighten bounds anyway if CglProbing there
4094  CglProbing* generator = NULL;
4095  int iGen;
4096  for (iGen=0;iGen<numberCutGenerators_;iGen++) {
4097    generator = dynamic_cast<CglProbing*>(generator_[iGen]->generator());
4098    if (generator)
4099      break;
4100  }
4101  int numberFixed=0;
4102  int numberTightened=0;
4103  int numberFixedByProbing=0;
4104  int numberTightenedByProbing=0;
4105  int printFrequency = (numberSolves+19)/20; // up to 20 messages
4106  int save[4];
4107  if (generator) {
4108    // set to cheaper and then restore at end
4109    save[0]=generator->getMaxPass();
4110    save[1]=generator->getMaxProbe();
4111    save[2]=generator->getMaxLook();
4112    save[3]=generator->rowCuts();
4113    generator->setMaxPass(1);
4114    generator->setMaxProbe(10);
4115    generator->setMaxLook(50);
4116    generator->setRowCuts(0);
4117   
4118    // Probing - return tight column bounds
4119    generator->generateCutsAndModify(*solver,cuts);
4120    const double * tightLower = generator->tightLower();
4121    const double * lower = solver->getColLower();
4122    const double * tightUpper = generator->tightUpper();
4123    const double * upper = solver->getColUpper();
4124    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4125      double newUpper = tightUpper[iColumn];
4126      double newLower = tightLower[iColumn];
4127      if (newUpper<upper[iColumn]-1.0e-8*(fabs(upper[iColumn])+1)||
4128          newLower>lower[iColumn]+1.0e-8*(fabs(lower[iColumn])+1)) {
4129        if (newUpper<newLower) {
4130          fprintf(stderr,"Problem is infeasible\n");
4131          return false;
4132        }
4133        if (newUpper==newLower) {
4134          numberFixed++;
4135          numberFixedByProbing++;
4136          solver->setColLower(iColumn,newLower);
4137          solver->setColUpper(iColumn,newUpper);
4138          printf("Column %d, new bounds %g %g\n",iColumn,
4139                 newLower,newUpper);
4140        } else if (vub[iColumn]) {
4141          numberTightened++;
4142          numberTightenedByProbing++;
4143          if (!solver->isInteger(iColumn)) {
4144            // relax
4145            newLower=CoinMax(lower[iColumn],
4146                                    newLower
4147                                    -1.0e-5*(fabs(lower[iColumn])+1));
4148            newUpper=CoinMin(upper[iColumn],
4149                                    newUpper
4150                                    +1.0e-5*(fabs(upper[iColumn])+1));
4151          }
4152          solver->setColLower(iColumn,newLower);
4153          solver->setColUpper(iColumn,newUpper);
4154        }
4155      }
4156    }
4157  }
4158  CoinWarmStart * ws = solver->getWarmStart();
4159  double * solution = new double [numberColumns];
4160  memcpy(solution,solver->getColSolution(),numberColumns*sizeof(double));
4161  for (iColumn=0;iColumn<numberColumns;iColumn++) 
4162    solver->setObjCoeff(iColumn,0.0);
4163  //solver->messageHandler()->setLogLevel(2);
4164  for (iVub=0;iVub<numberSolves;iVub++) {
4165    iColumn = which[iVub];
4166    int iTry;
4167    for (iTry=0;iTry<2;iTry++) {
4168      double saveUpper = solver->getColUpper()[iColumn];
4169      double saveLower = solver->getColLower()[iColumn];
4170      double value;
4171      if (iTry==1) {
4172        // try all way up
4173        solver->setObjCoeff(iColumn,-1.0);
4174      } else {
4175        // try all way down
4176        solver->setObjCoeff(iColumn,1.0);
4177      }
4178      solver->initialSolve();
4179      value = solver->getColSolution()[iColumn];
4180      bool change=false;
4181      if (iTry==1) {
4182        if (value<saveUpper-1.0e-4) {
4183          if (solver->isInteger(iColumn)) {
4184            value = floor(value+0.00001);
4185          } else {
4186            // relax a bit
4187            value=CoinMin(saveUpper,value+1.0e-5*(fabs(saveUpper)+1));
4188          }
4189          if (value-saveLower<1.0e-7) 
4190            value = saveLower; // make sure exactly same
4191          solver->setColUpper(iColumn,value);
4192          saveUpper=value;
4193          change=true;
4194        }
4195      } else {
4196        if (value>saveLower+1.0e-4) {
4197          if (solver->isInteger(iColumn)) {
4198            value = ceil(value-0.00001);
4199          } else {
4200            // relax a bit
4201            value=CoinMax(saveLower,value-1.0e-5*(fabs(saveLower)+1));
4202          }
4203          if (saveUpper-value<1.0e-7) 
4204            value = saveUpper; // make sure exactly same
4205          solver->setColLower(iColumn,value);
4206          saveLower=value;
4207          change=true;
4208        }
4209      }
4210      solver->setObjCoeff(iColumn,0.0);
4211      if (change) {
4212        if (saveUpper==saveLower) 
4213          numberFixed++;
4214        else
4215          numberTightened++;
4216        int saveFixed=numberFixed;
4217       
4218        int jColumn;
4219        if (generator) {
4220          // Probing - return tight column bounds
4221          cuts = OsiCuts();
4222          generator->generateCutsAndModify(*solver,cuts);
4223          const double * tightLower = generator->tightLower();
4224          const double * lower = solver->getColLower();
4225          const double * tightUpper = generator->tightUpper();
4226          const double * upper = solver->getColUpper();
4227          for (jColumn=0;jColumn<numberColumns;jColumn++) {
4228            double newUpper = tightUpper[jColumn];
4229            double newLower = tightLower[jColumn];
4230            if (newUpper<upper[jColumn]-1.0e-8*(fabs(upper[jColumn])+1)||
4231                newLower>lower[jColumn]+1.0e-8*(fabs(lower[jColumn])+1)) {
4232              if (newUpper<newLower) {
4233                fprintf(stderr,"Problem is infeasible\n");
4234                return false;
4235              }
4236              if (newUpper==newLower) {
4237                numberFixed++;
4238                numberFixedByProbing++;
4239                solver->setColLower(jColumn,newLower);
4240                solver->setColUpper(jColumn,newUpper);
4241              } else if (vub[jColumn]) {
4242                numberTightened++;
4243                numberTightenedByProbing++;
4244                if (!solver->isInteger(jColumn)) {
4245                  // relax
4246                  newLower=CoinMax(lower[jColumn],
4247                               newLower
4248                               -1.0e-5*(fabs(lower[jColumn])+1));
4249                  newUpper=CoinMin(upper[jColumn],
4250                               newUpper
4251                               +1.0e-5*(fabs(upper[jColumn])+1));
4252                }
4253                solver->setColLower(jColumn,newLower);
4254                solver->setColUpper(jColumn,newUpper);
4255              }
4256            }
4257          }
4258        }
4259        if (numberFixed>saveFixed) {
4260          // original solution may not be feasible
4261          // go back to true costs to solve if exists
4262          if (objective) {
4263            for (jColumn=0;jColumn<numberColumns;jColumn++) 
4264              solver->setObjCoeff(jColumn,objective[jColumn]);
4265          }
4266          solver->setColSolution(solution);
4267          solver->setWarmStart(ws);
4268          solver->resolve();
4269          if (!solver->isProvenOptimal()) {
4270            fprintf(stderr,"Problem is infeasible\n");
4271            return false;
4272          }
4273          delete ws;
4274          ws = solver->getWarmStart();
4275          memcpy(solution,solver->getColSolution(),
4276                 numberColumns*sizeof(double));
4277          for (jColumn=0;jColumn<numberColumns;jColumn++) 
4278            solver->setObjCoeff(jColumn,0.0);
4279        }
4280      }
4281      solver->setColSolution(solution);
4282      solver->setWarmStart(ws);
4283    }
4284    if (iVub%printFrequency==0) 
4285      handler_->message(CBC_VUB_PASS,messages_)
4286        <<iVub+1<<numberFixed<<numberTightened
4287        <<CoinMessageEol;
4288  }
4289  handler_->message(CBC_VUB_END,messages_)
4290    <<numberFixed<<numberTightened
4291    <<CoinMessageEol;
4292  delete ws;
4293  delete [] solution;
4294  // go back to true costs to solve if exists
4295  if (objective) {
4296    for (iColumn=0;iColumn<numberColumns;iColumn++) 
4297      solver_->setObjCoeff(iColumn,objective[iColumn]);
4298    delete [] objective;
4299  }
4300  delete [] vub;
4301  if (generator) {
4302    /*printf("Probing fixed %d and tightened %d\n",
4303           numberFixedByProbing,
4304           numberTightenedByProbing);*/
4305    if (generator_[iGen]->howOften()==-1&&
4306        (numberFixedByProbing+numberTightenedByProbing)*5>
4307        (numberFixed+numberTightened))
4308      generator_[iGen]->setHowOften(1000000+1);
4309    generator->setMaxPass(save[0]);
4310    generator->setMaxProbe(save[1]);
4311    generator->setMaxLook(save[2]);
4312    generator->setRowCuts(save[3]);
4313  }
4314
4315  if (solver!=solver_) {
4316    // move bounds across
4317    const double * lower = solver->getColLower();
4318    const double * upper = solver->getColUpper();
4319    const double * lowerOrig = solver_->getColLower();
4320    const double * upperOrig = solver_->getColUpper();
4321    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4322      solver_->setColLower(iColumn,CoinMax(lower[iColumn],lowerOrig[iColumn]));
4323      solver_->setColUpper(iColumn,CoinMin(upper[iColumn],upperOrig[iColumn]));
4324    }
4325    delete solver;
4326  }
4327  setCutoff(saveCutoff);
4328  return true;
4329}
4330/*
4331   Do Integer Presolve. Returns new model.
4332   I have to work out cleanest way of getting solution to
4333   original problem at end.  So this is very preliminary.
4334*/
4335CbcModel * 
4336CbcModel::integerPresolve(bool weak)
4337{
4338  status_ = 0;
4339  // solve LP
4340  //solver_->writeMps("bad");
4341  bool feasible = resolve();
4342
4343  CbcModel * newModel = NULL;
4344  if (feasible) {
4345
4346    // get a new model
4347    newModel = new CbcModel(*this);
4348    newModel->messageHandler()->setLogLevel(messageHandler()->logLevel());
4349
4350    feasible = newModel->integerPresolveThisModel(solver_,weak);
4351  }
4352  if (!feasible) {
4353    handler_->message(CBC_INFEAS,messages_)
4354    <<CoinMessageEol;
4355    status_ = 0;
4356    delete newModel;
4357    return NULL;
4358  } else {
4359    newModel->synchronizeModel(); // make sure everything that needs solver has it
4360    return newModel;
4361  }
4362}
4363/*
4364   Do Integer Presolve - destroying current model
4365*/
4366bool 
4367CbcModel::integerPresolveThisModel(OsiSolverInterface * originalSolver,
4368                                   bool weak)
4369{
4370  status_ = 0;
4371  // solve LP
4372  bool feasible = resolve();
4373
4374  bestObjective_=1.0e50;
4375  numberSolutions_=0;
4376  numberHeuristicSolutions_=0;
4377  double cutoff = getCutoff() ;
4378  double direction = solver_->getObjSense();
4379  if (cutoff < 1.0e20&&direction<0.0)
4380    messageHandler()->message(CBC_CUTOFF_WARNING1,
4381                                    messages())
4382                                      << cutoff << -cutoff << CoinMessageEol ;
4383  if (cutoff > bestObjective_)
4384    cutoff = bestObjective_ ;
4385  setCutoff(cutoff) ;
4386  int iColumn;
4387  int numberColumns = getNumCols();
4388  int originalNumberColumns = numberColumns;
4389  currentPassNumber_=0;
4390  synchronizeModel(); // make sure everything that needs solver has it
4391  // just point to solver_
4392  delete continuousSolver_;
4393  continuousSolver_ = solver_;
4394  // get a copy of original so we can fix bounds
4395  OsiSolverInterface * cleanModel = originalSolver->clone();
4396#ifdef CBC_DEBUG
4397  std::string problemName;
4398  cleanModel->getStrParam(OsiProbName,problemName);
4399  printf("Problem name - %s\n",problemName.c_str());
4400  cleanModel->activateRowCutDebugger(problemName.c_str());
4401  const OsiRowCutDebugger * debugger = cleanModel->getRowCutDebugger();
4402#endif
4403
4404  // array which points from original columns to presolved
4405  int * original = new int[numberColumns];
4406  // arrays giving bounds - only ones found by probing
4407  // rest will be found by presolve
4408  double * originalLower = new double[numberColumns];
4409  double * originalUpper = new double[numberColumns];
4410  {
4411    const double * lower = getColLower();
4412    const double * upper = getColUpper();
4413    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4414      original[iColumn]=iColumn;
4415      originalLower[iColumn] = lower[iColumn];
4416      originalUpper[iColumn] = upper[iColumn];
4417    }
4418  }
4419  findIntegers(true);
4420  // save original integers
4421  int * originalPriority = NULL;
4422  int * originalIntegers = new int[numberIntegers_];
4423  int originalNumberIntegers = numberIntegers_;
4424  memcpy(originalIntegers,integerVariable_,numberIntegers_*sizeof(int));
4425
4426  if (priority_) {
4427    originalPriority = new int[numberIntegers_];
4428    memcpy(originalPriority,priority_,numberIntegers_*sizeof(int));
4429    delete [] priority_;
4430    priority_=NULL;
4431  }
4432  int todo=20;
4433  if (weak)
4434    todo=1;
4435  while (currentPassNumber_<todo) {
4436   
4437    currentPassNumber_++;
4438    numberSolutions_=0;
4439    // this will be set false to break out of loop with presolved problem
4440    bool doIntegerPresolve=(currentPassNumber_!=20);
4441   
4442    // Current number of free integer variables
4443    // Get increment in solutions
4444    {
4445      const double * objective = cleanModel->getObjCoefficients();
4446      const double * lower = cleanModel->getColLower();
4447      const double * upper = cleanModel->getColUpper();
4448      double maximumCost=0.0;
4449      bool possibleMultiple=true;
4450      int numberChanged=0;
4451      for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
4452        if (originalUpper[iColumn]>originalLower[iColumn]) {
4453          if( cleanModel->isInteger(iColumn)) {
4454            maximumCost = CoinMax(maximumCost,fabs(objective[iColumn]));
4455          } else if (objective[iColumn]) {
4456            possibleMultiple=false;
4457          }
4458        }
4459        if (originalUpper[iColumn]<upper[iColumn]) {
4460#ifdef CBC_DEBUG
4461          printf("Changing upper bound on %d from %g to %g\n",
4462                 iColumn,upper[iColumn],originalUpper[iColumn]);
4463#endif
4464          cleanModel->setColUpper(iColumn,originalUpper[iColumn]);
4465          numberChanged++;
4466        }
4467        if (originalLower[iColumn]>lower[iColumn]) {
4468#ifdef CBC_DEBUG
4469          printf("Changing lower bound on %d from %g to %g\n",
4470                 iColumn,lower[iColumn],originalLower[iColumn]);
4471#endif
4472          cleanModel->setColLower(iColumn,originalLower[iColumn]);
4473          numberChanged++;
4474        }
4475      }
4476      // if first pass - always try
4477      if (currentPassNumber_==1)
4478        numberChanged += 1;
4479      if (possibleMultiple&&maximumCost) {
4480        int increment=0; 
4481        double multiplier = 2520.0;
4482        while (10.0*multiplier*maximumCost<1.0e8)
4483          multiplier *= 10.0;
4484        for (int j =0;j<originalNumberIntegers;j++) {
4485          iColumn = originalIntegers[j];
4486          if (originalUpper[iColumn]>originalLower[iColumn]) {
4487            if(objective[iColumn]) {
4488              double value = fabs(objective[iColumn])*multiplier;
4489              int nearest = (int) floor(value+0.5);
4490              if (fabs(value-floor(value+0.5))>1.0e-8) {
4491                increment=0;
4492                break; // no good
4493              } else if (!increment) {
4494                // first
4495                increment=nearest;
4496              } else {
4497                increment = gcd(increment,nearest);
4498              }
4499            }
4500          }
4501        }
4502        if (increment) {
4503          double value = increment;
4504          value /= multiplier;
4505          if (value*0.999>dblParam_[CbcCutoffIncrement]) {
4506            messageHandler()->message(CBC_INTEGERINCREMENT,messages())
4507              <<value
4508              <<CoinMessageEol;
4509            dblParam_[CbcCutoffIncrement]=value*0.999;
4510          }
4511        }
4512      }
4513      if (!numberChanged) {
4514        doIntegerPresolve=false; // not doing any better
4515      }
4516    }
4517#ifdef CBC_DEBUG
4518    if (debugger) 
4519      assert(debugger->onOptimalPath(*cleanModel));
4520#endif
4521#ifdef COIN_USE_CLP
4522    // do presolve - for now just clp but easy to get osi interface
4523    OsiClpSolverInterface * clpSolver
4524      = dynamic_cast<OsiClpSolverInterface *> (cleanModel);
4525    if (clpSolver) {
4526      ClpSimplex * clp = clpSolver->getModelPtr();
4527      clp->messageHandler()->setLogLevel(cleanModel->messageHandler()->logLevel());
4528      ClpPresolve pinfo;
4529      //printf("integerPresolve - temp switch off doubletons\n");
4530      //pinfo.setPresolveActions(4);
4531      ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e-8);
4532      if (!model2) {
4533        // presolve found to be infeasible
4534        feasible=false;
4535      } else {
4536        // update original array
4537        const int * originalColumns = pinfo.originalColumns();
4538        // just slot in new solver
4539        OsiClpSolverInterface * temp = new OsiClpSolverInterface(model2,true);
4540        numberColumns = temp->getNumCols();
4541        for (iColumn=0;iColumn<originalNumberColumns;iColumn++)
4542          original[iColumn]=-1;
4543        for (iColumn=0;iColumn<numberColumns;iColumn++)
4544          original[originalColumns[iColumn]]=iColumn;
4545        // copy parameters
4546        temp->copyParameters(*solver_);
4547        // and specialized ones
4548        temp->setSpecialOptions(clpSolver->specialOptions());
4549        delete solver_;
4550        solver_ = temp;
4551        setCutoff(cutoff);
4552        deleteObjects();
4553        if (!numberObjects_) {
4554          // Nothing left
4555          doIntegerPresolve=false;
4556          weak=true;
4557          break;
4558        }
4559        synchronizeModel(); // make sure everything that needs solver has it
4560        // just point to solver_
4561        continuousSolver_ = solver_;
4562        feasible=resolve();
4563        if (!feasible||!doIntegerPresolve||weak) break;
4564        // see if we can get solution by heuristics
4565        int found=-1;
4566        int iHeuristic;
4567        double * newSolution = new double [numberColumns];
4568        double heuristicValue=getCutoff();
4569        for (iHeuristic=0;iHeuristic<numberHeuristics_;iHeuristic++) {
4570          double saveValue=heuristicValue;
4571          int ifSol = heuristic_[iHeuristic]->solution(heuristicValue,
4572                                                       newSolution);
4573          if (ifSol>0) {
4574            // better solution found
4575            found=iHeuristic;
4576          } else if (ifSol<0) {
4577            heuristicValue = saveValue;
4578          }
4579        }
4580        if (found >= 0) {
4581          // We probably already have a current solution, but just in case ...
4582          int numberColumns = getNumCols() ;
4583          if (!currentSolution_)
4584            currentSolution_ = new double[numberColumns] ;
4585          // better solution save
4586          setBestSolution(CBC_ROUNDING,heuristicValue,
4587                          newSolution);
4588          // update cutoff
4589          cutoff = getCutoff();
4590        }
4591        delete [] newSolution;
4592        // Space for type of cuts
4593        int maximumWhich=1000;
4594        int * whichGenerator = new int[maximumWhich];
4595        // save number of rows
4596        numberRowsAtContinuous_ = getNumRows();
4597        maximumNumberCuts_=0;
4598        currentNumberCuts_=0;
4599        delete [] addedCuts_;
4600        addedCuts_ = NULL;
4601       
4602        // maximum depth for tree walkback
4603        maximumDepth_=10;
4604        delete [] walkback_;
4605        walkback_ = new CbcNodeInfo * [maximumDepth_];
4606       
4607        OsiCuts cuts;
4608        int numberOldActiveCuts=0;
4609        int numberNewCuts = 0;
4610        feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
4611                                 NULL,numberOldActiveCuts,numberNewCuts,
4612                                 maximumWhich, whichGenerator);
4613        currentNumberCuts_=numberNewCuts;
4614        delete [] whichGenerator;
4615        delete [] walkback_;
4616        walkback_ = NULL;
4617        delete [] addedCuts_;
4618        addedCuts_=NULL;
4619        if (feasible) {
4620          // fix anything in original which integer presolve fixed
4621          // for now just integers
4622          const double * lower = solver_->getColLower();
4623          const double * upper = solver_->getColUpper();
4624          int i;
4625          for (i=0;i<originalNumberIntegers;i++) {
4626            iColumn = originalIntegers[i];
4627            int jColumn = original[iColumn];
4628            if (jColumn >= 0) {
4629              if (upper[jColumn]<originalUpper[iColumn]) 
4630                originalUpper[iColumn]  = upper[jColumn];
4631              if (lower[jColumn]>originalLower[iColumn]) 
4632                originalLower[iColumn]  = lower[jColumn];
4633            }
4634          }
4635        }
4636      }
4637    }
4638#endif
4639    if (!feasible||!doIntegerPresolve) {
4640      break;
4641    }
4642  }
4643  //solver_->writeMps("xx");
4644  delete cleanModel;
4645  // create new priorities and save list of original columns
4646  if (originalPriority) {
4647    priority_ = new int[numberIntegers_];
4648    int i;
4649    int number=0;
4650    // integer variables are in same order if they exist at all
4651    for (i=0;i<originalNumberIntegers;i++) {
4652      iColumn = originalIntegers[i];
4653      int newColumn=original[iColumn];
4654      if (newColumn >= 0) 
4655        priority_[number++]=originalPriority[i];
4656    }
4657    assert (number==numberIntegers_);
4658    delete [] originalPriority;
4659  }
4660  delete [] originalIntegers;
4661  numberColumns = getNumCols();
4662  delete [] originalColumns_;
4663  originalColumns_ = new int[numberColumns];
4664  numberColumns=0;
4665  for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
4666    int jColumn = original[iColumn];
4667    if (jColumn >= 0) 
4668      originalColumns_[numberColumns++]=iColumn;
4669  }
4670  delete [] original;
4671  delete [] originalLower;
4672  delete [] originalUpper;
4673 
4674  deleteObjects();
4675  synchronizeModel(); // make sure everything that needs solver has it
4676  continuousSolver_=NULL;
4677  currentNumberCuts_=0;
4678  return feasible;
4679}
4680// Put back information into original model - after integerpresolve
4681void 
4682CbcModel::originalModel(CbcModel * presolvedModel,bool weak)
4683{
4684  solver_->copyParameters(*(presolvedModel->solver_));
4685  bestObjective_ = presolvedModel->bestObjective_;
4686  delete [] bestSolution_;
4687  findIntegers(true);
4688  if (presolvedModel->bestSolution_) {
4689    int numberColumns = getNumCols();
4690    int numberOtherColumns = presolvedModel->getNumCols();
4691    //bestSolution_ = new double[numberColumns];
4692    // set up map
4693    int * back = new int[numberColumns];
4694    int i;
4695    for (i=0;i<numberColumns;i++)
4696      back[i]=-1;
4697    for (i=0;i<numberOtherColumns;i++)
4698      back[presolvedModel->originalColumns_[i]]=i;
4699    int iColumn;
4700    // set ones in presolved model to values
4701    double * otherSolution = presolvedModel->bestSolution_;
4702    //const double * lower = getColLower();
4703    for (i=0;i<numberIntegers_;i++) {
4704      iColumn = integerVariable_[i];
4705      int jColumn = back[iColumn];
4706      //bestSolution_[iColumn]=lower[iColumn];
4707      if (jColumn >= 0) {
4708        double value=floor(otherSolution[jColumn]+0.5);
4709        solver_->setColLower(iColumn,value);
4710        solver_->setColUpper(iColumn,value);
4711        //bestSolution_[iColumn]=value;
4712      }
4713    }
4714    delete [] back;
4715#if 0
4716    // ** looks as if presolve needs more intelligence
4717    // do presolve - for now just clp but easy to get osi interface
4718    OsiClpSolverInterface * clpSolver
4719      = dynamic_cast<OsiClpSolverInterface *> (solver_);
4720    assert (clpSolver);
4721    ClpSimplex * clp = clpSolver->getModelPtr();
4722    Presolve pinfo;
4723    ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e-8);
4724    model2->primal(1);
4725    pinfo.postsolve(true);
4726    const double * solution = solver_->getColSolution();
4727    for (i=0;i<numberIntegers_;i++) {
4728      iColumn = integerVariable_[i];
4729      double value=floor(solution[iColumn]+0.5);
4730      solver_->setColLower(iColumn,value);
4731      solver_->setColUpper(iColumn,value);
4732    }
4733#else
4734    if (!weak) {
4735      // for now give up
4736      int save = numberCutGenerators_;
4737      numberCutGenerators_=0;
4738      bestObjective_=1.0e100;
4739      branchAndBound();
4740      numberCutGenerators_=save;
4741    }
4742#endif
4743    if (bestSolution_) {
4744      // solve problem
4745      resolve();
4746      // should be feasible
4747      int numberIntegerInfeasibilities;
4748      int numberObjectInfeasibilities;
4749      if (!currentSolution_)
4750        currentSolution_ = new double[numberColumns] ;
4751      assert(feasibleSolution(numberIntegerInfeasibilities,
4752                              numberObjectInfeasibilities));
4753    }
4754  } else {
4755    bestSolution_=NULL;
4756  }
4757  numberSolutions_=presolvedModel->numberSolutions_;
4758  numberHeuristicSolutions_=presolvedModel->numberHeuristicSolutions_;
4759  numberNodes_ = presolvedModel->numberNodes_;
4760  numberIterations_ = presolvedModel->numberIterations_;
4761  status_ = presolvedModel->status_;
4762  synchronizeModel();
4763} 
4764// Pass in Message handler (not deleted at end)
4765void 
4766CbcModel::passInMessageHandler(CoinMessageHandler * handler)
4767{
4768  if (defaultHandler_) {
4769    delete handler_;
4770    handler_ = NULL;
4771  }
4772  defaultHandler_=false;
4773  handler_=handler;
4774}
4775void 
4776CbcModel::passInTreeHandler(CbcTree & tree)
4777{
4778  delete tree_;
4779  tree_ = tree.clone();
4780}
4781// Make sure region there
4782void 
4783CbcModel::reserveCurrentSolution()
4784{
4785  int numberColumns = getNumCols() ;
4786  if (!currentSolution_)
4787    currentSolution_ = new double[numberColumns] ;
4788}
4789/* For passing in an CbcModel to do a sub Tree (with derived tree handlers).
4790   Passed in model must exist for duration of branch and bound
4791*/
4792void 
4793CbcModel::passInSubTreeModel(CbcModel & model)
4794{
4795  subTreeModel_=&model;
4796}
4797// For retrieving a copy of subtree model with given OsiSolver or NULL
4798CbcModel * 
4799CbcModel::subTreeModel(OsiSolverInterface * solver) const
4800{
4801  const CbcModel * subModel=subTreeModel_;
4802  if (!subModel)
4803    subModel=this;
4804  // Get new copy
4805  CbcModel * newModel = new CbcModel(*subModel);
4806  if (solver)
4807    newModel->assignSolver(solver);
4808  return newModel;
4809}
4810//#############################################################################
4811// Set/Get Application Data
4812// This is a pointer that the application can store into and retrieve
4813// from the solverInterface.
4814// This field is the application to optionally define and use.
4815//#############################################################################
4816
4817void CbcModel::setApplicationData(void * appData)
4818{
4819  appData_ = appData;
4820}
4821//-----------------------------------------------------------------------------
4822void * CbcModel::getApplicationData() const
4823{
4824  return appData_;
4825}
4826/*  create a submodel from partially fixed problem
4827
4828The method creates a new clean model with given bounds.
4829*/
4830CbcModel * 
4831CbcModel::cleanModel(const double * lower, const double * upper)
4832{
4833  OsiSolverInterface * solver = continuousSolver_->clone();
4834
4835  int numberIntegers = numberIntegers_;
4836  const int * integerVariable = integerVariable_;
4837 
4838  int i;
4839  for (i=0;i<numberIntegers;i++) {
4840    int iColumn=integerVariable[i];
4841    const CbcObject * object = object_[i];
4842    const CbcSimpleInteger * integerObject = 
4843      dynamic_cast<const  CbcSimpleInteger *> (object);
4844    assert(integerObject);
4845    // get original bounds
4846    double originalLower = integerObject->originalLowerBound();
4847    double originalUpper = integerObject->originalUpperBound();
4848    solver->setColLower(iColumn,CoinMax(lower[iColumn],originalLower));
4849    solver->setColUpper(iColumn,CoinMin(upper[iColumn],originalUpper));
4850  }
4851  CbcModel * model = new CbcModel(*solver);
4852  // off some messages
4853  if (handler_->logLevel()<=1) {
4854    model->messagesPointer()->setDetailMessage(3,9);
4855    model->messagesPointer()->setDetailMessage(3,6);
4856    model->messagesPointer()->setDetailMessage(3,4);
4857    model->messagesPointer()->setDetailMessage(3,1);
4858    model->messagesPointer()->setDetailMessage(3,13);
4859    model->messagesPointer()->setDetailMessage(3,14);
4860    model->messagesPointer()->setDetailMessage(3,3007);
4861  }
4862  // Cuts
4863  for ( i = 0;i<numberCutGenerators_;i++) {
4864    int howOften = generator_[i]->howOftenInSub();
4865    if (howOften>-100) {
4866      CbcCutGenerator * generator = virginGenerator_[i];
4867      CglCutGenerator * cglGenerator = generator->generator();
4868      model->addCutGenerator(cglGenerator,howOften,
4869                              generator->cutGeneratorName(),
4870                              generator->normal(),
4871                              generator->atSolution(),
4872                              generator->whenInfeasible(),
4873                              -100, generator->whatDepthInSub(),-1);
4874    }
4875  }
4876  double cutoff = getCutoff();
4877  model->setCutoff(cutoff);
4878  return model;
4879}
4880/* Invoke the branch & cut algorithm on partially fixed problem
4881   
4882   The method uses a subModel created by cleanModel. The search
4883   ends when the tree is exhausted or maximum nodes is reached.
4884
4885   If better solution found then it is saved.
4886   
4887   Returns 0 if search completed and solution, 1 if not completed and solution,
4888   2 if completed and no solution, 3 if not completed and no solution.
4889   
4890   Normally okay to do subModel immediately followed by subBranchandBound
4891   (== other form of subBranchAndBound)
4892   but may need to get at model for advanced features.
4893   
4894   Deletes model
4895   
4896*/
4897 
4898int 
4899CbcModel::subBranchAndBound(CbcModel * model,
4900                            CbcModel * presolvedModel,
4901                            int maximumNodes)
4902{
4903  int i;
4904  double cutoff=model->getCutoff();
4905  CbcModel * model2;
4906  if (presolvedModel) 
4907    model2=presolvedModel;
4908  else
4909    model2=model;
4910  // Do complete search
4911 
4912  for (i=0;i<numberHeuristics_;i++) {
4913    model2->addHeuristic(heuristic_[i]);
4914    model2->heuristic(i)->resetModel(model2);
4915  }
4916  // Definition of node choice
4917  model2->setNodeComparison(nodeCompare_->clone());
4918  //model2->solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
4919  model2->messageHandler()->setLogLevel(CoinMax(0,handler_->logLevel()-1));
4920  //model2->solver()->messageHandler()->setLogLevel(2);
4921  model2->setMaximumCutPassesAtRoot(maximumCutPassesAtRoot_);
4922  model2->setPrintFrequency(50);
4923  model2->setIntParam(CbcModel::CbcMaxNumNode,maximumNodes);
4924  model2->branchAndBound();
4925  delete model2->nodeComparison();
4926  if (model2->getMinimizationObjValue()>cutoff) {
4927    // no good
4928    if (model!=model2)
4929      delete model2;
4930    delete model;
4931    return 2;
4932  }
4933  if (model!=model2) {
4934    // get back solution
4935    model->originalModel(model2,false);
4936    delete model2;
4937  }
4938  int status;
4939  if (model->getMinimizationObjValue()<cutoff&&model->bestSolution()) {
4940    double objValue = model->getObjValue();
4941    const double * solution = model->bestSolution();
4942    setBestSolution(CBC_TREE_SOL,objValue,solution);
4943    status = 0;
4944  } else {
4945    status=2;
4946  }
4947  if (model->status())
4948    status ++ ; // not finished search
4949  delete model;
4950  return status;
4951}
4952/* Invoke the branch & cut algorithm on partially fixed problem
4953   
4954The method creates a new model with given bounds, presolves it
4955then proceeds to explore the branch & cut search tree. The search
4956ends when the tree is exhausted or maximum nodes is reached.
4957Returns 0 if search completed and solution, 1 if not completed and solution,
49582 if completed and no solution, 3 if not completed and no solution.
4959*/
4960int 
4961CbcModel::subBranchAndBound(const double * lower, const double * upper,
4962                            int maximumNodes)
4963{
4964  OsiSolverInterface * solver = continuousSolver_->clone();
4965
4966  int numberIntegers = numberIntegers_;
4967  const int * integerVariable = integerVariable_;
4968 
4969  int i;
4970  for (i=0;i<numberIntegers;i++) {
4971    int iColumn=integerVariable[i];
4972    const CbcObject * object = object_[i];
4973    const CbcSimpleInteger * integerObject = 
4974      dynamic_cast<const  CbcSimpleInteger *> (object);
4975    assert(integerObject);
4976    // get original bounds
4977    double originalLower = integerObject->originalLowerBound();
4978    double originalUpper = integerObject->originalUpperBound();
4979    solver->setColLower(iColumn,CoinMax(lower[iColumn],originalLower));
4980    solver->setColUpper(iColumn,CoinMin(upper[iColumn],originalUpper));
4981  }
4982  CbcModel model(*solver);
4983  // off some messages
4984  if (handler_->logLevel()<=1) {
4985    model.messagesPointer()->setDetailMessage(3,9);
4986    model.messagesPointer()->setDetailMessage(3,6);
4987    model.messagesPointer()->setDetailMessage(3,4);
4988    model.messagesPointer()->setDetailMessage(3,1);
4989    model.messagesPointer()->setDetailMessage(3,3007);
4990  }
4991  double cutoff = getCutoff();
4992  model.setCutoff(cutoff);
4993  // integer presolve
4994  CbcModel * model2 = model.integerPresolve(false);
4995  if (!model2||!model2->getNumRows()) {
4996    delete model2;
4997    delete solver;
4998    return 2;
4999  }
5000  if (handler_->logLevel()>1)
5001    printf("Reduced model has %d rows and %d columns\n",
5002           model2->getNumRows(),model2->getNumCols());
5003  // Do complete search
5004 
5005  // Cuts
5006  for ( i = 0;i<numberCutGenerators_;i++) {
5007    int howOften = generator_[i]->howOftenInSub();
5008    if (howOften>-100) {
5009      CbcCutGenerator * generator = virginGenerator_[i];
5010      CglCutGenerator * cglGenerator = generator->generator();
5011      model2->addCutGenerator(cglGenerator,howOften,
5012                              generator->cutGeneratorName(),
5013                              generator->normal(),
5014                              generator->atSolution(),
5015                              generator->whenInfeasible(),
5016                              -100, generator->whatDepthInSub(),-1);
5017    }
5018  }
5019  for (i=0;i<numberHeuristics_;i++) {
5020    model2->addHeuristic(heuristic_[i]);
5021    model2->heuristic(i)->resetModel(model2);
5022  }
5023  // Definition of node choice
5024  model2->setNodeComparison(nodeCompare_->clone());
5025  //model2->solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
5026  model2->messageHandler()->setLogLevel(CoinMax(0,handler_->logLevel()-1));
5027  //model2->solver()->messageHandler()->setLogLevel(2);
5028  model2->setMaximumCutPassesAtRoot(maximumCutPassesAtRoot_);
5029  model2->setPrintFrequency(50);
5030  model2->setIntParam(CbcModel::CbcMaxNumNode,maximumNodes);
5031  model2->branchAndBound();
5032  delete model2->nodeComparison();
5033  if (model2->getMinimizationObjValue()>cutoff) {
5034    // no good
5035    delete model2;
5036    delete solver;
5037    return 2;
5038  }
5039  // get back solution
5040  model.originalModel(model2,false);
5041  delete model2;
5042  int status;
5043  if (model.getMinimizationObjValue()<cutoff&&model.bestSolution()) {
5044    double objValue = model.getObjValue();
5045    const double * solution = model.bestSolution();
5046    setBestSolution(CBC_TREE_SOL,objValue,solution);
5047    status = 0;
5048  } else {
5049    status=2;
5050  }
5051  if (model.status())
5052    status ++ ; // not finished search
5053  delete solver;
5054  return status;
5055}
5056// Set a pointer to a row cut which will be added instead of normal branching.
5057void 
5058CbcModel::setNextRowCut(const OsiRowCut & cut)
5059{ 
5060  nextRowCut_=new OsiRowCut(cut);
5061}
5062/* Process root node and return a strengthened model
5063   
5064The method assumes that initialSolve() has been called to solve the
5065LP relaxation. It processes the root node and then returns a pointer
5066to the strengthened model (or NULL if infeasible)
5067*/
5068OsiSolverInterface * 
5069CbcModel::strengthenedModel()
5070{
5071/*
5072  Assume we're done, and see if we're proven wrong.
5073*/
5074/*
5075  Scan the variables, noting the integer variables. Create an
5076  CbcSimpleInteger object for each integer variable.
5077*/
5078  findIntegers(false) ;
5079/*
5080  Ensure that objects on the lists of CbcObjects, heuristics, and cut
5081  generators attached to this model all refer to this model.
5082*/
5083  synchronizeModel() ;
5084
5085  // Set so we can tell we are in initial phase in resolve
5086  continuousObjective_ = -COIN_DBL_MAX ;
5087/*
5088  Solve the relaxation.
5089
5090  Apparently there are circumstances where this will be non-trivial --- i.e.,
5091  we've done something since initialSolve that's trashed the solution to the
5092  continuous relaxation.
5093*/
5094  bool feasible = resolve() ;
5095/*
5096  If the linear relaxation of the root is infeasible, bail out now. Otherwise,
5097  continue with processing the root node.
5098*/
5099  if (!feasible)
5100  { handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
5101    return NULL; }
5102  // Save objective (just so user can access it)
5103  originalContinuousObjective_ = solver_->getObjValue();
5104
5105/*
5106  Begin setup to process a feasible root node.
5107*/
5108  bestObjective_ = CoinMin(bestObjective_,1.0e50) ;
5109  numberSolutions_ = 0 ;
5110  numberHeuristicSolutions_ = 0 ;
5111  // Everything is minimization
5112  double cutoff=getCutoff() ;
5113  double direction = solver_->getObjSense() ;
5114  if (cutoff < 1.0e20&&direction<0.0)
5115    messageHandler()->message(CBC_CUTOFF_WARNING1,
5116                                    messages())
5117                                      << cutoff << -cutoff << CoinMessageEol ;
5118  if (cutoff > bestObjective_)
5119    cutoff = bestObjective_ ;
5120  setCutoff(cutoff) ;
5121/*
5122  We probably already have a current solution, but just in case ...
5123*/
5124  int numberColumns = getNumCols() ;
5125  if (!currentSolution_)
5126    currentSolution_ = new double[numberColumns] ;
5127/*
5128  Create a copy of the solver, thus capturing the original (root node)
5129  constraint system (aka the continuous system).
5130*/
5131  continuousSolver_ = solver_->clone() ;
5132  numberRowsAtContinuous_ = getNumRows() ;
5133/*
5134  Check the objective to see if we can deduce a nontrivial increment. If
5135  it's better than the current value for CbcCutoffIncrement, it'll be
5136  installed.
5137*/
5138  analyzeObjective() ;
5139/*
5140  Set up for cut generation. addedCuts_ holds the cuts which are relevant for
5141  the active subproblem. whichGenerator will be used to record the generator
5142  that produced a given cut.
5143*/
5144  int maximumWhich = 1000 ;
5145  int * whichGenerator = new int[maximumWhich] ;
5146  maximumNumberCuts_ = 0 ;
5147  currentNumberCuts_ = 0 ;
5148  delete [] addedCuts_ ;
5149  addedCuts_ = NULL ;
5150  /* 
5151  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
5152  lifting. It will iterate a generate/reoptimise loop (including reduced cost
5153  fixing) until no cuts are generated, the change in objective falls off,  or
5154  the limit on the number of rounds of cut generation is exceeded.
5155
5156  At the end of all this, any cuts will be recorded in cuts and also
5157  installed in the solver's constraint system. We'll have reoptimised, and
5158  removed any slack cuts (numberOldActiveCuts and numberNewCuts have been
5159  adjusted accordingly).
5160
5161  Tell cut generators they can be a bit more aggressive at root node
5162
5163*/
5164  int iCutGenerator;
5165  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
5166    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
5167    generator->setAggressiveness(generator->getAggressiveness()+100);
5168  }
5169  OsiCuts cuts ;
5170  int numberOldActiveCuts = 0 ;
5171  int numberNewCuts = 0 ;
5172  { int iObject ;
5173    int preferredWay ;
5174    int numberUnsatisfied = 0 ;
5175    double integerTolerance = getIntegerTolerance() ;
5176
5177    memcpy(currentSolution_,solver_->getColSolution(),
5178           numberColumns*sizeof(double)) ;
5179
5180    for (iObject = 0 ; iObject < numberObjects_ ; iObject++)
5181    { double infeasibility =
5182          object_[iObject]->infeasibility(preferredWay) ;
5183      if (infeasibility > integerTolerance) numberUnsatisfied++ ; }
5184    if (numberUnsatisfied)
5185    { feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
5186                               NULL,numberOldActiveCuts,numberNewCuts,
5187                               maximumWhich, whichGenerator) ; } }
5188/*
5189  We've taken the continuous relaxation as far as we can.
5190*/
5191
5192  OsiSolverInterface * newSolver=NULL;
5193  if (feasible) {
5194    // make copy of current solver
5195    newSolver = solver_->clone();
5196  }
5197/*
5198  Clean up dangling objects. continuousSolver_ may already be toast.
5199*/
5200  delete [] whichGenerator ;
5201  delete [] walkback_ ;
5202  walkback_ = NULL ;
5203  delete [] addedCuts_ ;
5204  addedCuts_ = NULL ;
5205  if (continuousSolver_)
5206  { delete continuousSolver_ ;
5207    continuousSolver_ = NULL ; }
5208/*
5209  Destroy global cuts by replacing with an empty OsiCuts object.
5210*/
5211  globalCuts_= OsiCuts() ;
5212 
5213  return newSolver; 
5214}
5215// Just update objectiveValue
5216void CbcModel::setBestObjectiveValue( double objectiveValue)
5217{
5218  bestObjective_=objectiveValue;
5219}
5220double 
5221CbcModel::getBestPossibleObjValue() const
5222{ 
5223  return CoinMin(bestPossibleObjective_,bestObjective_) * solver_->getObjSense() ;
5224}
Note: See TracBrowser for help on using the repository browser.