source: trunk/CbcModel.cpp @ 97

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

for OsiCbc?

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