source: trunk/CbcModel.cpp @ 84

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

new message

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