source: trunk/CbcModel.cpp @ 82

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

correct for pseudocosts

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 171.1 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    if (ifObject)
3477      memcpy(priority_+numberIntegers_,priorities,
3478             (numberObjects_-numberIntegers_)*sizeof(int));
3479    else
3480      memcpy(priority_,priorities,numberIntegers_*sizeof(int));
3481  }
3482}
3483
3484// Delete all object information
3485void 
3486CbcModel::deleteObjects()
3487{
3488  delete [] priority_;
3489  priority_=NULL;
3490  int i;
3491  for (i=0;i<numberObjects_;i++)
3492    delete object_[i];
3493  delete [] object_;
3494  object_ = NULL;
3495  numberObjects_=0;
3496  findIntegers(true);
3497}
3498
3499/*!
3500  Ensure all attached objects (CbcObjects, heuristics, and cut
3501  generators) point to this model.
3502*/
3503void CbcModel::synchronizeModel()
3504{
3505  int i;
3506  for (i=0;i<numberHeuristics_;i++) 
3507    heuristic_[i]->setModel(this);
3508  for (i=0;i<numberObjects_;i++)
3509    object_[i]->setModel(this);
3510  for (i=0;i<numberCutGenerators_;i++)
3511    generator_[i]->refreshModel(this);
3512}
3513
3514// Fill in integers and create objects
3515
3516/**
3517  The routine first does a scan to count the number of integer variables.
3518  It then creates an array, integerVariable_, to store the indices of the
3519  integer variables, and an array of `objects', one for each variable.
3520
3521  The scan is repeated, this time recording the index of each integer
3522  variable in integerVariable_, and creating an CbcSimpleInteger object that
3523  contains information about the integer variable. Initially, this is just
3524  the index and upper & lower bounds.
3525
3526  \todo
3527  Note the assumption in cbc that the first numberIntegers_ objects are
3528  CbcSimpleInteger. In particular, the code which handles the startAgain
3529  case assumes that if the object_ array exists it can simply replace the first
3530  numberInteger_ objects. This is arguably unsafe.
3531
3532  I am going to re-order if necessary
3533*/
3534
3535void 
3536CbcModel::findIntegers(bool startAgain)
3537{
3538  assert(solver_);
3539/*
3540  No need to do this if we have previous information, unless forced to start
3541  over.
3542*/
3543  if (numberIntegers_&&!startAgain&&object_)
3544    return;
3545/*
3546  Clear out the old integer variable list, then count the number of integer
3547  variables.
3548*/
3549  delete [] integerVariable_;
3550  numberIntegers_=0;
3551  int numberColumns = getNumCols();
3552  int iColumn;
3553  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3554    if (isInteger(iColumn)) 
3555      numberIntegers_++;
3556  }
3557  // Find out how many old non-integer objects there are
3558  int nObjects=0;
3559  CbcObject ** oldObject = object_;
3560  int iObject;
3561  for (iObject = 0;iObject<numberObjects_;iObject++) {
3562    CbcSimpleInteger * obj =
3563      dynamic_cast <CbcSimpleInteger *>(oldObject[iObject]) ;
3564    if (obj) 
3565      delete oldObject[iObject];
3566    else
3567      oldObject[nObjects++]=oldObject[iObject];
3568  }
3569   
3570/*
3571  Found any? Allocate an array to hold the indices of the integer variables.
3572  Make a large enough array for all objects
3573*/
3574  object_ = new CbcObject * [numberIntegers_+nObjects];
3575  numberObjects_=numberIntegers_+nObjects;;
3576  integerVariable_ = new int [numberIntegers_];
3577/*
3578  Walk the variables again, filling in the indices and creating objects for
3579  the integer variables. Initially, the objects hold the index and upper &
3580  lower bounds.
3581*/
3582  numberIntegers_=0;
3583  for (iColumn=0;iColumn<numberColumns;iColumn++) {
3584    if(isInteger(iColumn)) {
3585      object_[numberIntegers_] =
3586        new CbcSimpleInteger(this,numberIntegers_,iColumn);
3587      integerVariable_[numberIntegers_++]=iColumn;
3588    }
3589  }
3590  // Now append other objects
3591  memcpy(object_+numberIntegers_,oldObject,nObjects*sizeof(CbcObject *));
3592  // Delete old array (just array)
3593  delete [] oldObject;
3594 
3595  if (!numberObjects_)
3596      handler_->message(CBC_NOINT,messages_) << CoinMessageEol ;
3597}
3598
3599/* Add in any object information (objects are cloned - owner can delete
3600   originals */
3601void 
3602CbcModel::addObjects(int numberObjects, CbcObject ** objects)
3603{
3604  // If integers but not enough objects fudge
3605  if (numberIntegers_>numberObjects_)
3606    findIntegers(true);
3607  /* But if incoming objects inherit from simple integer we just want
3608     to replace */
3609  int numberColumns = solver_->getNumCols();
3610  /** mark is -1 if not integer, >=0 if using existing simple integer and
3611      >=numberColumns if using new integer */
3612  int * mark = new int[numberColumns];
3613  int i;
3614  for (i=0;i<numberColumns;i++)
3615    mark[i]=-1;
3616  int newNumberObjects = numberObjects;
3617  int newIntegers=0;
3618  for (i=0;i<numberObjects;i++) { 
3619    CbcSimpleInteger * obj =
3620      dynamic_cast <CbcSimpleInteger *>(objects[i]) ;
3621    if (obj) {
3622      int iColumn = obj->columnNumber();
3623      mark[iColumn]=i+numberColumns;
3624      newIntegers++;
3625    }
3626  }
3627  // and existing
3628  for (i=0;i<numberObjects_;i++) { 
3629    CbcSimpleInteger * obj =
3630      dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
3631    if (obj) {
3632      int iColumn = obj->columnNumber();
3633      if (mark[iColumn]<0) {
3634        newIntegers++;
3635        newNumberObjects++;
3636        mark[iColumn]=i;
3637      }
3638    }
3639  } 
3640  delete [] integerVariable_;
3641  integerVariable_=NULL;
3642  if (newIntegers!=numberIntegers_) 
3643    printf("changing number of integers from %d to %d\n",
3644           numberIntegers_,newIntegers);
3645  numberIntegers_ = newIntegers;
3646  integerVariable_ = new int [numberIntegers_];
3647  CbcObject ** temp  = new CbcObject * [newNumberObjects];
3648  // Put integers first
3649  newIntegers=0;
3650  numberIntegers_=0;
3651  for (i=0;i<numberColumns;i++) {
3652    int which = mark[i];
3653    if (which>=0) {
3654      if (!isInteger(i)) {
3655        newIntegers++;
3656        solver_->setInteger(i);
3657      }
3658      if (which<numberColumns) {
3659        temp[numberIntegers_]=object_[which];
3660        object_[which]=NULL;
3661      } else {
3662        temp[numberIntegers_]=objects[which-numberColumns]->clone();
3663        temp[numberIntegers_]->setModel(this);
3664      }
3665      integerVariable_[numberIntegers_++]=i;
3666    }
3667  }
3668  if (newIntegers)
3669    printf("%d variables were declared integer\n",newIntegers);
3670  int n=numberIntegers_;
3671  // Now rest of old
3672  for (i=0;i<numberObjects_;i++) { 
3673    if (object_[i]) {
3674      CbcSimpleInteger * obj =
3675        dynamic_cast <CbcSimpleInteger *>(object_[i]) ;
3676      if (obj) {
3677        delete object_[i];
3678      } else {
3679        temp[n++]=object_[i];
3680      }
3681    }
3682  }
3683  // and rest of new
3684  for (i=0;i<numberObjects;i++) { 
3685    CbcSimpleInteger * obj =
3686      dynamic_cast <CbcSimpleInteger *>(objects[i]) ;
3687    if (!obj) {
3688      temp[n]=objects[i]->clone();
3689      temp[n++]->setModel(this);
3690    }
3691  }
3692  delete [] mark;
3693  delete [] object_;
3694  object_ = temp;
3695  assert (n==newNumberObjects);
3696  numberObjects_ = newNumberObjects;
3697}
3698
3699/**
3700  This routine sets the objective cutoff value used for fathoming and
3701  determining monotonic variables.
3702
3703  If the fathoming discipline is strict, a small tolerance is added to the
3704  new cutoff. This avoids problems due to roundoff when the target value
3705  is exact. The common example would be an IP with only integer variables in
3706  the objective. If the target is set to the exact value z of the optimum,
3707  it's possible to end up fathoming an ancestor of the solution because the
3708  solver returns z+epsilon.
3709
3710  Determining if strict fathoming is needed is best done by analysis.
3711  In cbc, that's analyseObjective. The default is false.
3712
3713  In cbc we always minimize so add epsilon
3714*/
3715
3716void CbcModel::setCutoff (double value)
3717
3718{
3719#if 0
3720  double tol = 0 ;
3721  int fathomStrict = getIntParam(CbcFathomDiscipline) ;
3722  if (fathomStrict == 1)
3723  { solver_->getDblParam(OsiDualTolerance,tol) ;
3724  tol = tol*(1+fabs(value)) ;
3725 
3726  value += tol ; }
3727#endif
3728  // Solvers know about direction
3729  double direction = solver_->getObjSense();
3730  solver_->setDblParam(OsiDualObjectiveLimit,value*direction); }
3731
3732
3733
3734/*
3735  Call this to really test if a valid solution can be feasible. The cutoff is
3736  passed in as a parameter so that we don't need to worry here after swapping
3737  solvers.  The solution is assumed to be numberColumns in size.  If
3738  fixVariables is true then the bounds of the continuous solver are updated.
3739  The routine returns the objective value determined by reoptimizing from
3740  scratch. If the solution is rejected, this will be worse than the cutoff.
3741
3742  TODO: There's an issue with getting the correct cutoff value: We update the
3743        cutoff in the regular solver, but not in continuousSolver_. But our only
3744        use for continuousSolver_ is verifying candidate solutions. Would it
3745        make sense to update the cutoff? Then we wouldn't need to step around
3746        isDualObjectiveLimitReached().
3747*/
3748double 
3749CbcModel::checkSolution (double cutoff, const double *solution,
3750                         bool fixVariables)
3751
3752{ int numberColumns = solver_->getNumCols();
3753
3754  /*
3755    Grab the continuous solver (the pristine copy of the problem, made before
3756    starting to work on the root node). Save the bounds on the variables.
3757    Install the solution passed as a parameter, and copy it to the model's
3758    currentSolution_.
3759   
3760    TODO: This is a belt-and-suspenders approach. Once the code has settled
3761          a bit, we can cast a critical eye here.
3762  */
3763  OsiSolverInterface * saveSolver = solver_;
3764  if (continuousSolver_)
3765    solver_ = continuousSolver_;
3766  // move solution to continuous copy
3767  solver_->setColSolution(solution);
3768  // Put current solution in safe place
3769  memcpy(currentSolution_,solver_->getColSolution(),
3770         numberColumns*sizeof(double));
3771  //solver_->messageHandler()->setLogLevel(4);
3772
3773  // save original bounds
3774  double * saveUpper = new double[numberColumns];
3775  double * saveLower = new double[numberColumns];
3776  memcpy(saveUpper,getColUpper(),numberColumns*sizeof(double));
3777  memcpy(saveLower,getColLower(),numberColumns*sizeof(double));
3778
3779  /*
3780    Run through the objects and use feasibleRegion() to set variable bounds
3781    so as to fix the variables specified in the objects at their value in this
3782    solution. Since the object list contains (at least) one object for every
3783    integer variable, this has the effect of fixing all integer variables.
3784  */
3785  int i;
3786  for (i=0;i<numberObjects_;i++)
3787    object_[i]->feasibleRegion();
3788  /*
3789    Remove any existing warm start information to be sure there is no
3790    residual influence on initialSolve().
3791  */
3792  CoinWarmStartBasis *slack =
3793      dynamic_cast<CoinWarmStartBasis *>(solver_->getEmptyWarmStart()) ;
3794  solver_->setWarmStart(slack);
3795  delete slack ;
3796  // Give a hint not to do scaling
3797  //bool saveTakeHint;
3798  //OsiHintStrength saveStrength;
3799  //bool gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
3800  //assert (gotHint);
3801  //solver_->setHintParam(OsiDoScale,false,OsiHintTry);
3802  solver_->initialSolve();
3803  //solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
3804  if (!solver_->isProvenOptimal())
3805    { std::cout << "checkSolution infeas! Retrying wihout scaling."
3806              << std::endl ;
3807    bool saveTakeHint;
3808    OsiHintStrength saveStrength;
3809    bool savePrintHint;
3810    solver_->writeMps("infeas");
3811    bool gotHint = (solver_->getHintParam(OsiDoReducePrint,savePrintHint,saveStrength));
3812    gotHint = (solver_->getHintParam(OsiDoScale,saveTakeHint,saveStrength));
3813    solver_->setHintParam(OsiDoScale,false,OsiHintTry);
3814    solver_->setHintParam(OsiDoReducePrint,false,OsiHintTry) ;
3815    solver_->initialSolve();
3816    solver_->setHintParam(OsiDoScale,saveTakeHint,saveStrength);
3817    solver_->setHintParam(OsiDoReducePrint,savePrintHint,OsiHintTry) ;
3818    }
3819  assert(solver_->isProvenOptimal());
3820  double objectiveValue = solver_->getObjValue()*solver_->getObjSense();
3821
3822  /*
3823    Check that the solution still beats the objective cutoff.
3824
3825    If it passes, make a copy of the primal variable values and do some
3826    cleanup and checks:
3827    + Values of all variables are are within original bounds and values of
3828      all integer variables are within tolerance of integral.
3829    + There are no constraint violations.
3830    There really should be no need for the check against original bounds.
3831    Perhaps an opportunity for a sanity check?
3832  */
3833  if (solver_->isProvenOptimal() && objectiveValue <= cutoff)
3834  { solution = solver_->getColSolution() ;
3835    memcpy(currentSolution_,solution,numberColumns*sizeof(double)) ;
3836
3837    const double * rowLower = solver_->getRowLower() ;
3838    const double * rowUpper = solver_->getRowUpper() ;
3839    int numberRows = solver_->getNumRows() ;
3840    double *rowActivity = new double[numberRows] ;
3841    memset(rowActivity,0,numberRows*sizeof(double)) ;
3842
3843    double integerTolerance = getIntegerTolerance() ;
3844
3845    int iColumn;
3846    for (iColumn = 0 ; iColumn < numberColumns ; iColumn++)
3847    { double value = currentSolution_[iColumn] ;
3848      value = CoinMax(value, saveLower[iColumn]) ;
3849      value = CoinMin(value, saveUpper[iColumn]) ;
3850      if (solver_->isInteger(iColumn)) 
3851        assert(fabs(value-currentSolution_[iColumn]) <= integerTolerance) ;
3852      currentSolution_[iColumn] = value ; }
3853   
3854    solver_->getMatrixByCol()->times(currentSolution_,rowActivity) ;
3855    double primalTolerance ;
3856    solver_->getDblParam(OsiPrimalTolerance,primalTolerance) ;
3857    double largestInfeasibility =0.0;
3858    for (i=0 ; i < numberRows ; i++) {
3859      largestInfeasibility = CoinMax(largestInfeasibility,
3860                                 rowLower[i]-rowActivity[i]);
3861      largestInfeasibility = CoinMax(largestInfeasibility,
3862                                 rowActivity[i]-rowUpper[i]);
3863    }
3864    if (largestInfeasibility>100.0*primalTolerance)
3865      handler_->message(CBC_NOTFEAS3, messages_)
3866        << largestInfeasibility << CoinMessageEol ;
3867
3868    delete [] rowActivity ; }
3869  else
3870  { objectiveValue=1.0e50 ; }
3871  /*
3872    Regardless of what we think of the solution, we may need to restore the
3873    original bounds of the continuous solver. Unfortunately, const'ness
3874    prevents us from simply reversing the memcpy used to make these snapshots.
3875  */
3876  if (!fixVariables)
3877  { for (int iColumn = 0 ; iColumn < numberColumns ; iColumn++)
3878    { solver_->setColLower(iColumn,saveLower[iColumn]) ;
3879      solver_->setColUpper(iColumn,saveUpper[iColumn]) ; } }
3880  delete [] saveLower;
3881  delete [] saveUpper;
3882   
3883  /*
3884    Restore the usual solver.
3885  */
3886  solver_=saveSolver;
3887  return objectiveValue;
3888}
3889
3890/*
3891  Call this routine from anywhere when a solution is found. The solution
3892  vector is assumed to contain one value for each structural variable.
3893
3894  The first action is to run checkSolution() to confirm the objective and
3895  feasibility. If this check causes the solution to be rejected, we're done.
3896  If fixVariables = true, the variable bounds held by the continuous solver
3897  will be left fixed to the values in the solution; otherwise they are
3898  restored to the original values.
3899
3900  If the solution is accepted, install it as the best solution.
3901
3902  The routine also contains a hook to run any cut generators that are flagged
3903  to run when a new solution is discovered. There's a potential hazard because
3904  the cut generators see the continuous solver >after< possible restoration of
3905  original bounds (which may well invalidate the solution).
3906*/
3907
3908void
3909CbcModel::setBestSolution (CBC_Message how,
3910                           double & objectiveValue, const double *solution,
3911                           bool fixVariables)
3912
3913{ double cutoff = getCutoff() ;
3914
3915/*
3916  Double check the solution to catch pretenders.
3917*/
3918  objectiveValue = checkSolution(cutoff,solution,fixVariables);
3919  if (objectiveValue > cutoff)
3920  { if (objectiveValue>1.0e30)
3921      handler_->message(CBC_NOTFEAS1, messages_) << CoinMessageEol ;
3922    else
3923      handler_->message(CBC_NOTFEAS2, messages_)
3924        << objectiveValue << cutoff << CoinMessageEol ; }
3925/*
3926  We have a winner. Install it as the new incumbent.
3927  Bump the objective cutoff value and solution counts. Give the user the
3928  good news.
3929*/
3930  else
3931  { bestObjective_ = objectiveValue;
3932    int numberColumns = solver_->getNumCols();
3933    if (!bestSolution_)
3934      bestSolution_ = new double[numberColumns];
3935    memcpy(bestSolution_,currentSolution_,numberColumns*sizeof(double));
3936
3937    cutoff = bestObjective_-dblParam_[CbcCutoffIncrement];
3938    // This is not correct - that way cutoff can go up if maximization
3939    //double direction = solver_->getObjSense();
3940    //setCutoff(cutoff*direction);
3941    setCutoff(cutoff);
3942
3943    if (how==CBC_ROUNDING)
3944      numberHeuristicSolutions_++;
3945    numberSolutions_++;
3946
3947    handler_->message(how,messages_)
3948      <<bestObjective_<<numberIterations_
3949      <<numberNodes_
3950      <<CoinMessageEol;
3951/*
3952  Now step through the cut generators and see if any of them are flagged to
3953  run when a new solution is discovered. Only global cuts are useful. (The
3954  solution being evaluated may not correspond to the current location in the
3955  search tree --- discovered by heuristic, for example.)
3956*/
3957    OsiCuts theseCuts;
3958    int i;
3959    int lastNumberCuts=0;
3960    for (i=0;i<numberCutGenerators_;i++) {
3961      if (generator_[i]->atSolution()) {
3962        generator_[i]->generateCuts(theseCuts,true,NULL);
3963        int numberCuts = theseCuts.sizeRowCuts();
3964        for (int j=lastNumberCuts;j<numberCuts;j++) {
3965          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j);
3966          if (thisCut->globallyValid()) {
3967            if ((specialOptions_&1)!=0) {
3968              /* As these are global cuts -
3969                 a) Always get debugger object
3970                 b) Not fatal error to cutoff optimal (if we have just got optimal)
3971              */
3972              const OsiRowCutDebugger *debugger = solver_->getRowCutDebuggerAlways() ;
3973              if (debugger) {
3974                if(debugger->invalidCut(*thisCut))
3975                  printf("ZZZZ Global cut - cuts off optimal solution!\n");
3976              }
3977            }
3978            // add to global list
3979            globalCuts_.insert(*thisCut);
3980            generator_[i]->incrementNumberCutsInTotal();
3981          }
3982        }
3983      }
3984    }
3985    int numberCuts = theseCuts.sizeColCuts();
3986    for (i=0;i<numberCuts;i++) {
3987      const OsiColCut * thisCut = theseCuts.colCutPtr(i);
3988      if (thisCut->globallyValid()) {
3989        // add to global list
3990        globalCuts_.insert(*thisCut);
3991      }
3992    }
3993  }
3994
3995  return ; }
3996
3997
3998/* Test the current solution for feasibility.
3999
4000   Calculate the number of standard integer infeasibilities, then scan the
4001   remaining objects to see if any of them report infeasibilities.
4002
4003   Currently (2003.08) the only object besides SimpleInteger is Clique, hence
4004   the comments about `odd ones' infeasibilities.
4005*/
4006bool 
4007CbcModel::feasibleSolution(int & numberIntegerInfeasibilities,
4008                        int & numberObjectInfeasibilities) const
4009{
4010  int numberUnsatisfied=0;
4011  double sumUnsatisfied=0.0;
4012  int preferredWay;
4013  int j;
4014  // Put current solution in safe place
4015  memcpy(currentSolution_,solver_->getColSolution(),
4016         solver_->getNumCols()*sizeof(double));
4017  for (j=0;j<numberIntegers_;j++) {
4018    const CbcObject * object = object_[j];
4019    double infeasibility = object->infeasibility(preferredWay);
4020    if (infeasibility) {
4021      assert (infeasibility>0);
4022      numberUnsatisfied++;
4023      sumUnsatisfied += infeasibility;
4024    }
4025  }
4026  numberIntegerInfeasibilities = numberUnsatisfied;
4027  for (;j<numberObjects_;j++) {
4028    const CbcObject * object = object_[j];
4029    double infeasibility = object->infeasibility(preferredWay);
4030    if (infeasibility) {
4031      assert (infeasibility>0);
4032      numberUnsatisfied++;
4033      sumUnsatisfied += infeasibility;
4034    }
4035  }
4036  numberObjectInfeasibilities = numberUnsatisfied-numberIntegerInfeasibilities;
4037  return (!numberUnsatisfied);
4038}
4039
4040/* For all vubs see if we can tighten bounds by solving Lp's
4041   type - 0 just vubs
4042   1 all (could be very slow)
4043   -1 just vubs where variable away from bound
4044   Returns false if not feasible
4045*/
4046bool 
4047CbcModel::tightenVubs(int type, bool allowMultipleBinary, double useCutoff)
4048{
4049
4050  CoinPackedMatrix matrixByRow(*solver_->getMatrixByRow());
4051  int numberRows = solver_->getNumRows();
4052  int numberColumns = solver_->getNumCols();
4053
4054  int iRow,iColumn;
4055
4056  // Row copy
4057  //const double * elementByRow = matrixByRow.getElements();
4058  const int * column = matrixByRow.getIndices();
4059  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
4060  const int * rowLength = matrixByRow.getVectorLengths();
4061
4062  const double * colUpper = solver_->getColUpper();
4063  const double * colLower = solver_->getColLower();
4064  //const double * rowUpper = solver_->getRowUpper();
4065  //const double * rowLower = solver_->getRowLower();
4066
4067  const double * objective = solver_->getObjCoefficients();
4068  //double direction = solver_->getObjSense();
4069  const double * colsol = solver_->getColSolution();
4070
4071  int numberVub=0;
4072  int * continuous = new int[numberColumns];
4073  if (type >= 0) {
4074    double * sort = new double[numberColumns];
4075    for (iRow=0;iRow<numberRows;iRow++) {
4076      int j;
4077      int numberBinary=0;
4078      int numberUnsatisfiedBinary=0;
4079      int numberContinuous=0;
4080      int iCont=-1;
4081      double weight=1.0e30;
4082      for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
4083        int iColumn = column[j];
4084        if (colUpper[iColumn]-colLower[iColumn]>1.0e-8) {
4085          if (solver_->isFreeBinary(iColumn)) {
4086            numberBinary++;
4087            /* For sort I make naive assumption:
4088               x - a * delta <=0 or
4089               -x + a * delta >= 0
4090            */
4091            if (colsol[iColumn]>colLower[iColumn]+1.0e-6&&
4092                colsol[iColumn]<colUpper[iColumn]-1.0e-6) {
4093              numberUnsatisfiedBinary++;
4094              weight = CoinMin(weight,fabs(objective[iColumn]));
4095            }
4096          } else {
4097            numberContinuous++;
4098            iCont=iColumn;
4099          }
4100        }
4101      }
4102      if (numberContinuous==1&&numberBinary) {
4103        if (numberBinary==1||allowMultipleBinary) {
4104          // treat as vub
4105          if (!numberUnsatisfiedBinary)
4106            weight=-1.0; // at end
4107          sort[numberVub]=-weight;
4108          continuous[numberVub++] = iCont;
4109        }
4110      }
4111    }
4112    if (type>0) {
4113      // take so many
4114      CoinSort_2(sort,sort+numberVub,continuous);
4115      numberVub = CoinMin(numberVub,type);
4116    }
4117    delete [] sort;
4118  } else {
4119    for (iColumn=0;iColumn<numberColumns;iColumn++) 
4120      continuous[iColumn]=iColumn;
4121    numberVub=numberColumns;
4122  }
4123  bool feasible = tightenVubs(numberVub,continuous,useCutoff);
4124  delete [] continuous;
4125
4126  return feasible;
4127}
4128// This version is just handed a list of variables
4129bool 
4130CbcModel::tightenVubs(int numberSolves, const int * which,
4131                      double useCutoff)
4132{
4133
4134  int numberColumns = solver_->getNumCols();
4135
4136  int iColumn;
4137 
4138  OsiSolverInterface * solver = solver_;
4139  double saveCutoff = getCutoff() ;
4140 
4141  double * objective = new double[numberColumns];
4142  memcpy(objective,solver_->getObjCoefficients(),numberColumns*sizeof(double));
4143  double direction = solver_->getObjSense();
4144 
4145  // add in objective if there is a cutoff
4146  if (useCutoff<1.0e30) {
4147    // get new version of model
4148    solver = solver_->clone();
4149    CoinPackedVector newRow;
4150    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4151      solver->setObjCoeff(iColumn,0.0); // zero out in new model
4152      if (objective[iColumn]) 
4153        newRow.insert(iColumn,direction * objective[iColumn]);
4154     
4155    }
4156    solver->addRow(newRow,-COIN_DBL_MAX,useCutoff);
4157    // signal no objective
4158    delete [] objective;
4159    objective=NULL;
4160  }
4161  setCutoff(COIN_DBL_MAX);
4162
4163
4164  bool * vub = new bool [numberColumns];
4165  int iVub;
4166
4167  // mark vub columns
4168  for (iColumn=0;iColumn<numberColumns;iColumn++) 
4169    vub[iColumn]=false;
4170  for (iVub=0;iVub<numberSolves;iVub++) 
4171    vub[which[iVub]]=true;
4172  OsiCuts cuts;
4173  // First tighten bounds anyway if CglProbing there
4174  CglProbing* generator = NULL;
4175  int iGen;
4176  for (iGen=0;iGen<numberCutGenerators_;iGen++) {
4177    generator = dynamic_cast<CglProbing*>(generator_[iGen]->generator());
4178    if (generator)
4179      break;
4180  }
4181  int numberFixed=0;
4182  int numberTightened=0;
4183  int numberFixedByProbing=0;
4184  int numberTightenedByProbing=0;
4185  int printFrequency = (numberSolves+19)/20; // up to 20 messages
4186  int save[4];
4187  if (generator) {
4188    // set to cheaper and then restore at end
4189    save[0]=generator->getMaxPass();
4190    save[1]=generator->getMaxProbe();
4191    save[2]=generator->getMaxLook();
4192    save[3]=generator->rowCuts();
4193    generator->setMaxPass(1);
4194    generator->setMaxProbe(10);
4195    generator->setMaxLook(50);
4196    generator->setRowCuts(0);
4197   
4198    // Probing - return tight column bounds
4199    generator->generateCutsAndModify(*solver,cuts);
4200    const double * tightLower = generator->tightLower();
4201    const double * lower = solver->getColLower();
4202    const double * tightUpper = generator->tightUpper();
4203    const double * upper = solver->getColUpper();
4204    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4205      double newUpper = tightUpper[iColumn];
4206      double newLower = tightLower[iColumn];
4207      if (newUpper<upper[iColumn]-1.0e-8*(fabs(upper[iColumn])+1)||
4208          newLower>lower[iColumn]+1.0e-8*(fabs(lower[iColumn])+1)) {
4209        if (newUpper<newLower) {
4210          fprintf(stderr,"Problem is infeasible\n");
4211          return false;
4212        }
4213        if (newUpper==newLower) {
4214          numberFixed++;
4215          numberFixedByProbing++;
4216          solver->setColLower(iColumn,newLower);
4217          solver->setColUpper(iColumn,newUpper);
4218          printf("Column %d, new bounds %g %g\n",iColumn,
4219                 newLower,newUpper);
4220        } else if (vub[iColumn]) {
4221          numberTightened++;
4222          numberTightenedByProbing++;
4223          if (!solver->isInteger(iColumn)) {
4224            // relax
4225            newLower=CoinMax(lower[iColumn],
4226                                    newLower
4227                                    -1.0e-5*(fabs(lower[iColumn])+1));
4228            newUpper=CoinMin(upper[iColumn],
4229                                    newUpper
4230                                    +1.0e-5*(fabs(upper[iColumn])+1));
4231          }
4232          solver->setColLower(iColumn,newLower);
4233          solver->setColUpper(iColumn,newUpper);
4234        }
4235      }
4236    }
4237  }
4238  CoinWarmStart * ws = solver->getWarmStart();
4239  double * solution = new double [numberColumns];
4240  memcpy(solution,solver->getColSolution(),numberColumns*sizeof(double));
4241  for (iColumn=0;iColumn<numberColumns;iColumn++) 
4242    solver->setObjCoeff(iColumn,0.0);
4243  //solver->messageHandler()->setLogLevel(2);
4244  for (iVub=0;iVub<numberSolves;iVub++) {
4245    iColumn = which[iVub];
4246    int iTry;
4247    for (iTry=0;iTry<2;iTry++) {
4248      double saveUpper = solver->getColUpper()[iColumn];
4249      double saveLower = solver->getColLower()[iColumn];
4250      double value;
4251      if (iTry==1) {
4252        // try all way up
4253        solver->setObjCoeff(iColumn,-1.0);
4254      } else {
4255        // try all way down
4256        solver->setObjCoeff(iColumn,1.0);
4257      }
4258      solver->initialSolve();
4259      value = solver->getColSolution()[iColumn];
4260      bool change=false;
4261      if (iTry==1) {
4262        if (value<saveUpper-1.0e-4) {
4263          if (solver->isInteger(iColumn)) {
4264            value = floor(value+0.00001);
4265          } else {
4266            // relax a bit
4267            value=CoinMin(saveUpper,value+1.0e-5*(fabs(saveUpper)+1));
4268          }
4269          if (value-saveLower<1.0e-7) 
4270            value = saveLower; // make sure exactly same
4271          solver->setColUpper(iColumn,value);
4272          saveUpper=value;
4273          change=true;
4274        }
4275      } else {
4276        if (value>saveLower+1.0e-4) {
4277          if (solver->isInteger(iColumn)) {
4278            value = ceil(value-0.00001);
4279          } else {
4280            // relax a bit
4281            value=CoinMax(saveLower,value-1.0e-5*(fabs(saveLower)+1));
4282          }
4283          if (saveUpper-value<1.0e-7) 
4284            value = saveUpper; // make sure exactly same
4285          solver->setColLower(iColumn,value);
4286          saveLower=value;
4287          change=true;
4288        }
4289      }
4290      solver->setObjCoeff(iColumn,0.0);
4291      if (change) {
4292        if (saveUpper==saveLower) 
4293          numberFixed++;
4294        else
4295          numberTightened++;
4296        int saveFixed=numberFixed;
4297       
4298        int jColumn;
4299        if (generator) {
4300          // Probing - return tight column bounds
4301          cuts = OsiCuts();
4302          generator->generateCutsAndModify(*solver,cuts);
4303          const double * tightLower = generator->tightLower();
4304          const double * lower = solver->getColLower();
4305          const double * tightUpper = generator->tightUpper();
4306          const double * upper = solver->getColUpper();
4307          for (jColumn=0;jColumn<numberColumns;jColumn++) {
4308            double newUpper = tightUpper[jColumn];
4309            double newLower = tightLower[jColumn];
4310            if (newUpper<upper[jColumn]-1.0e-8*(fabs(upper[jColumn])+1)||
4311                newLower>lower[jColumn]+1.0e-8*(fabs(lower[jColumn])+1)) {
4312              if (newUpper<newLower) {
4313                fprintf(stderr,"Problem is infeasible\n");
4314                return false;
4315              }
4316              if (newUpper==newLower) {
4317                numberFixed++;
4318                numberFixedByProbing++;
4319                solver->setColLower(jColumn,newLower);
4320                solver->setColUpper(jColumn,newUpper);
4321              } else if (vub[jColumn]) {
4322                numberTightened++;
4323                numberTightenedByProbing++;
4324                if (!solver->isInteger(jColumn)) {
4325                  // relax
4326                  newLower=CoinMax(lower[jColumn],
4327                               newLower
4328                               -1.0e-5*(fabs(lower[jColumn])+1));
4329                  newUpper=CoinMin(upper[jColumn],
4330                               newUpper
4331                               +1.0e-5*(fabs(upper[jColumn])+1));
4332                }
4333                solver->setColLower(jColumn,newLower);
4334                solver->setColUpper(jColumn,newUpper);
4335              }
4336            }
4337          }
4338        }
4339        if (numberFixed>saveFixed) {
4340          // original solution may not be feasible
4341          // go back to true costs to solve if exists
4342          if (objective) {
4343            for (jColumn=0;jColumn<numberColumns;jColumn++) 
4344              solver->setObjCoeff(jColumn,objective[jColumn]);
4345          }
4346          solver->setColSolution(solution);
4347          solver->setWarmStart(ws);
4348          solver->resolve();
4349          if (!solver->isProvenOptimal()) {
4350            fprintf(stderr,"Problem is infeasible\n");
4351            return false;
4352          }
4353          delete ws;
4354          ws = solver->getWarmStart();
4355          memcpy(solution,solver->getColSolution(),
4356                 numberColumns*sizeof(double));
4357          for (jColumn=0;jColumn<numberColumns;jColumn++) 
4358            solver->setObjCoeff(jColumn,0.0);
4359        }
4360      }
4361      solver->setColSolution(solution);
4362      solver->setWarmStart(ws);
4363    }
4364    if (iVub%printFrequency==0) 
4365      handler_->message(CBC_VUB_PASS,messages_)
4366        <<iVub+1<<numberFixed<<numberTightened
4367        <<CoinMessageEol;
4368  }
4369  handler_->message(CBC_VUB_END,messages_)
4370    <<numberFixed<<numberTightened
4371    <<CoinMessageEol;
4372  delete ws;
4373  delete [] solution;
4374  // go back to true costs to solve if exists
4375  if (objective) {
4376    for (iColumn=0;iColumn<numberColumns;iColumn++) 
4377      solver_->setObjCoeff(iColumn,objective[iColumn]);
4378    delete [] objective;
4379  }
4380  delete [] vub;
4381  if (generator) {
4382    /*printf("Probing fixed %d and tightened %d\n",
4383           numberFixedByProbing,
4384           numberTightenedByProbing);*/
4385    if (generator_[iGen]->howOften()==-1&&
4386        (numberFixedByProbing+numberTightenedByProbing)*5>
4387        (numberFixed+numberTightened))
4388      generator_[iGen]->setHowOften(1000000+1);
4389    generator->setMaxPass(save[0]);
4390    generator->setMaxProbe(save[1]);
4391    generator->setMaxLook(save[2]);
4392    generator->setRowCuts(save[3]);
4393  }
4394
4395  if (solver!=solver_) {
4396    // move bounds across
4397    const double * lower = solver->getColLower();
4398    const double * upper = solver->getColUpper();
4399    const double * lowerOrig = solver_->getColLower();
4400    const double * upperOrig = solver_->getColUpper();
4401    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4402      solver_->setColLower(iColumn,CoinMax(lower[iColumn],lowerOrig[iColumn]));
4403      solver_->setColUpper(iColumn,CoinMin(upper[iColumn],upperOrig[iColumn]));
4404    }
4405    delete solver;
4406  }
4407  setCutoff(saveCutoff);
4408  return true;
4409}
4410/*
4411   Do Integer Presolve. Returns new model.
4412   I have to work out cleanest way of getting solution to
4413   original problem at end.  So this is very preliminary.
4414*/
4415CbcModel * 
4416CbcModel::integerPresolve(bool weak)
4417{
4418  status_ = 0;
4419  // solve LP
4420  //solver_->writeMps("bad");
4421  bool feasible = resolve();
4422
4423  CbcModel * newModel = NULL;
4424  if (feasible) {
4425
4426    // get a new model
4427    newModel = new CbcModel(*this);
4428    newModel->messageHandler()->setLogLevel(messageHandler()->logLevel());
4429
4430    feasible = newModel->integerPresolveThisModel(solver_,weak);
4431  }
4432  if (!feasible) {
4433    handler_->message(CBC_INFEAS,messages_)
4434    <<CoinMessageEol;
4435    status_ = 0;
4436    delete newModel;
4437    return NULL;
4438  } else {
4439    newModel->synchronizeModel(); // make sure everything that needs solver has it
4440    return newModel;
4441  }
4442}
4443/*
4444   Do Integer Presolve - destroying current model
4445*/
4446bool 
4447CbcModel::integerPresolveThisModel(OsiSolverInterface * originalSolver,
4448                                   bool weak)
4449{
4450  status_ = 0;
4451  // solve LP
4452  bool feasible = resolve();
4453
4454  bestObjective_=1.0e50;
4455  numberSolutions_=0;
4456  numberHeuristicSolutions_=0;
4457  double cutoff = getCutoff() ;
4458  double direction = solver_->getObjSense();
4459  if (cutoff < 1.0e20&&direction<0.0)
4460    messageHandler()->message(CBC_CUTOFF_WARNING1,
4461                                    messages())
4462                                      << cutoff << -cutoff << CoinMessageEol ;
4463  if (cutoff > bestObjective_)
4464    cutoff = bestObjective_ ;
4465  setCutoff(cutoff) ;
4466  int iColumn;
4467  int numberColumns = getNumCols();
4468  int originalNumberColumns = numberColumns;
4469  currentPassNumber_=0;
4470  synchronizeModel(); // make sure everything that needs solver has it
4471  // just point to solver_
4472  delete continuousSolver_;
4473  continuousSolver_ = solver_;
4474  // get a copy of original so we can fix bounds
4475  OsiSolverInterface * cleanModel = originalSolver->clone();
4476#ifdef CBC_DEBUG
4477  std::string problemName;
4478  cleanModel->getStrParam(OsiProbName,problemName);
4479  printf("Problem name - %s\n",problemName.c_str());
4480  cleanModel->activateRowCutDebugger(problemName.c_str());
4481  const OsiRowCutDebugger * debugger = cleanModel->getRowCutDebugger();
4482#endif
4483
4484  // array which points from original columns to presolved
4485  int * original = new int[numberColumns];
4486  // arrays giving bounds - only ones found by probing
4487  // rest will be found by presolve
4488  double * originalLower = new double[numberColumns];
4489  double * originalUpper = new double[numberColumns];
4490  {
4491    const double * lower = getColLower();
4492    const double * upper = getColUpper();
4493    for (iColumn=0;iColumn<numberColumns;iColumn++) {
4494      original[iColumn]=iColumn;
4495      originalLower[iColumn] = lower[iColumn];
4496      originalUpper[iColumn] = upper[iColumn];
4497    }
4498  }
4499  findIntegers(true);
4500  // save original integers
4501  int * originalPriority = NULL;
4502  int * originalIntegers = new int[numberIntegers_];
4503  int originalNumberIntegers = numberIntegers_;
4504  memcpy(originalIntegers,integerVariable_,numberIntegers_*sizeof(int));
4505
4506  if (priority_) {
4507    originalPriority = new int[numberIntegers_];
4508    memcpy(originalPriority,priority_,numberIntegers_*sizeof(int));
4509    delete [] priority_;
4510    priority_=NULL;
4511  }
4512  int todo=20;
4513  if (weak)
4514    todo=1;
4515  while (currentPassNumber_<todo) {
4516   
4517    currentPassNumber_++;
4518    numberSolutions_=0;
4519    // this will be set false to break out of loop with presolved problem
4520    bool doIntegerPresolve=(currentPassNumber_!=20);
4521   
4522    // Current number of free integer variables
4523    // Get increment in solutions
4524    {
4525      const double * objective = cleanModel->getObjCoefficients();
4526      const double * lower = cleanModel->getColLower();
4527      const double * upper = cleanModel->getColUpper();
4528      double maximumCost=0.0;
4529      bool possibleMultiple=true;
4530      int numberChanged=0;
4531      for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
4532        if (originalUpper[iColumn]>originalLower[iColumn]) {
4533          if( cleanModel->isInteger(iColumn)) {
4534            maximumCost = CoinMax(maximumCost,fabs(objective[iColumn]));
4535          } else if (objective[iColumn]) {
4536            possibleMultiple=false;
4537          }
4538        }
4539        if (originalUpper[iColumn]<upper[iColumn]) {
4540#ifdef CBC_DEBUG
4541          printf("Changing upper bound on %d from %g to %g\n",
4542                 iColumn,upper[iColumn],originalUpper[iColumn]);
4543#endif
4544          cleanModel->setColUpper(iColumn,originalUpper[iColumn]);
4545          numberChanged++;
4546        }
4547        if (originalLower[iColumn]>lower[iColumn]) {
4548#ifdef CBC_DEBUG
4549          printf("Changing lower bound on %d from %g to %g\n",
4550                 iColumn,lower[iColumn],originalLower[iColumn]);
4551#endif
4552          cleanModel->setColLower(iColumn,originalLower[iColumn]);
4553          numberChanged++;
4554        }
4555      }
4556      // if first pass - always try
4557      if (currentPassNumber_==1)
4558        numberChanged += 1;
4559      if (possibleMultiple&&maximumCost) {
4560        int increment=0; 
4561        double multiplier = 2520.0;
4562        while (10.0*multiplier*maximumCost<1.0e8)
4563          multiplier *= 10.0;
4564        for (int j =0;j<originalNumberIntegers;j++) {
4565          iColumn = originalIntegers[j];
4566          if (originalUpper[iColumn]>originalLower[iColumn]) {
4567            if(objective[iColumn]) {
4568              double value = fabs(objective[iColumn])*multiplier;
4569              int nearest = (int) floor(value+0.5);
4570              if (fabs(value-floor(value+0.5))>1.0e-8) {
4571                increment=0;
4572                break; // no good
4573              } else if (!increment) {
4574                // first
4575                increment=nearest;
4576              } else {
4577                increment = gcd(increment,nearest);
4578              }
4579            }
4580          }
4581        }
4582        if (increment) {
4583          double value = increment;
4584          value /= multiplier;
4585          if (value*0.999>dblParam_[CbcCutoffIncrement]) {
4586            messageHandler()->message(CBC_INTEGERINCREMENT,messages())
4587              <<value
4588              <<CoinMessageEol;
4589            dblParam_[CbcCutoffIncrement]=value*0.999;
4590          }
4591        }
4592      }
4593      if (!numberChanged) {
4594        doIntegerPresolve=false; // not doing any better
4595      }
4596    }
4597#ifdef CBC_DEBUG
4598    if (debugger) 
4599      assert(debugger->onOptimalPath(*cleanModel));
4600#endif
4601#ifdef COIN_USE_CLP
4602    // do presolve - for now just clp but easy to get osi interface
4603    OsiClpSolverInterface * clpSolver
4604      = dynamic_cast<OsiClpSolverInterface *> (cleanModel);
4605    if (clpSolver) {
4606      ClpSimplex * clp = clpSolver->getModelPtr();
4607      clp->messageHandler()->setLogLevel(cleanModel->messageHandler()->logLevel());
4608      ClpPresolve pinfo;
4609      //printf("integerPresolve - temp switch off doubletons\n");
4610      //pinfo.setPresolveActions(4);
4611      ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e-8);
4612      if (!model2) {
4613        // presolve found to be infeasible
4614        feasible=false;
4615      } else {
4616        // update original array
4617        const int * originalColumns = pinfo.originalColumns();
4618        // just slot in new solver
4619        OsiClpSolverInterface * temp = new OsiClpSolverInterface(model2,true);
4620        numberColumns = temp->getNumCols();
4621        for (iColumn=0;iColumn<originalNumberColumns;iColumn++)
4622          original[iColumn]=-1;
4623        for (iColumn=0;iColumn<numberColumns;iColumn++)
4624          original[originalColumns[iColumn]]=iColumn;
4625        // copy parameters
4626        temp->copyParameters(*solver_);
4627        // and specialized ones
4628        temp->setSpecialOptions(clpSolver->specialOptions());
4629        delete solver_;
4630        solver_ = temp;
4631        setCutoff(cutoff);
4632        deleteObjects();
4633        if (!numberObjects_) {
4634          // Nothing left
4635          doIntegerPresolve=false;
4636          weak=true;
4637          break;
4638        }
4639        synchronizeModel(); // make sure everything that needs solver has it
4640        // just point to solver_
4641        continuousSolver_ = solver_;
4642        feasible=resolve();
4643        if (!feasible||!doIntegerPresolve||weak) break;
4644        // see if we can get solution by heuristics
4645        int found=-1;
4646        int iHeuristic;
4647        double * newSolution = new double [numberColumns];
4648        double heuristicValue=getCutoff();
4649        for (iHeuristic=0;iHeuristic<numberHeuristics_;iHeuristic++) {
4650          double saveValue=heuristicValue;
4651          int ifSol = heuristic_[iHeuristic]->solution(heuristicValue,
4652                                                       newSolution);
4653          if (ifSol>0) {
4654            // better solution found
4655            found=iHeuristic;
4656          } else if (ifSol<0) {
4657            heuristicValue = saveValue;
4658          }
4659        }
4660        if (found >= 0) {
4661          // We probably already have a current solution, but just in case ...
4662          int numberColumns = getNumCols() ;
4663          if (!currentSolution_)
4664            currentSolution_ = new double[numberColumns] ;
4665          // better solution save
4666          setBestSolution(CBC_ROUNDING,heuristicValue,
4667                          newSolution);
4668          // update cutoff
4669          cutoff = getCutoff();
4670        }
4671        delete [] newSolution;
4672        // Space for type of cuts
4673        int maximumWhich=1000;
4674        int * whichGenerator = new int[maximumWhich];
4675        // save number of rows
4676        numberRowsAtContinuous_ = getNumRows();
4677        maximumNumberCuts_=0;
4678        currentNumberCuts_=0;
4679        delete [] addedCuts_;
4680        addedCuts_ = NULL;
4681       
4682        // maximum depth for tree walkback
4683        maximumDepth_=10;
4684        delete [] walkback_;
4685        walkback_ = new CbcNodeInfo * [maximumDepth_];
4686       
4687        OsiCuts cuts;
4688        int numberOldActiveCuts=0;
4689        int numberNewCuts = 0;
4690        feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
4691                                 NULL,numberOldActiveCuts,numberNewCuts,
4692                                 maximumWhich, whichGenerator);
4693        currentNumberCuts_=numberNewCuts;
4694        delete [] whichGenerator;
4695        delete [] walkback_;
4696        walkback_ = NULL;
4697        delete [] addedCuts_;
4698        addedCuts_=NULL;
4699        if (feasible) {
4700          // fix anything in original which integer presolve fixed
4701          // for now just integers
4702          const double * lower = solver_->getColLower();
4703          const double * upper = solver_->getColUpper();
4704          int i;
4705          for (i=0;i<originalNumberIntegers;i++) {
4706            iColumn = originalIntegers[i];
4707            int jColumn = original[iColumn];
4708            if (jColumn >= 0) {
4709              if (upper[jColumn]<originalUpper[iColumn]) 
4710                originalUpper[iColumn]  = upper[jColumn];
4711              if (lower[jColumn]>originalLower[iColumn]) 
4712                originalLower[iColumn]  = lower[jColumn];
4713            }
4714          }
4715        }
4716      }
4717    }
4718#endif
4719    if (!feasible||!doIntegerPresolve) {
4720      break;
4721    }
4722  }
4723  //solver_->writeMps("xx");
4724  delete cleanModel;
4725  // create new priorities and save list of original columns
4726  if (originalPriority) {
4727    priority_ = new int[numberIntegers_];
4728    int i;
4729    int number=0;
4730    // integer variables are in same order if they exist at all
4731    for (i=0;i<originalNumberIntegers;i++) {
4732      iColumn = originalIntegers[i];
4733      int newColumn=original[iColumn];
4734      if (newColumn >= 0) 
4735        priority_[number++]=originalPriority[i];
4736    }
4737    assert (number==numberIntegers_);
4738    delete [] originalPriority;
4739  }
4740  delete [] originalIntegers;
4741  numberColumns = getNumCols();
4742  delete [] originalColumns_;
4743  originalColumns_ = new int[numberColumns];
4744  numberColumns=0;
4745  for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
4746    int jColumn = original[iColumn];
4747    if (jColumn >= 0) 
4748      originalColumns_[numberColumns++]=iColumn;
4749  }
4750  delete [] original;
4751  delete [] originalLower;
4752  delete [] originalUpper;
4753 
4754  deleteObjects();
4755  synchronizeModel(); // make sure everything that needs solver has it
4756  continuousSolver_=NULL;
4757  currentNumberCuts_=0;
4758  return feasible;
4759}
4760// Put back information into original model - after integerpresolve
4761void 
4762CbcModel::originalModel(CbcModel * presolvedModel,bool weak)
4763{
4764  solver_->copyParameters(*(presolvedModel->solver_));
4765  bestObjective_ = presolvedModel->bestObjective_;
4766  delete [] bestSolution_;
4767  findIntegers(true);
4768  if (presolvedModel->bestSolution_) {
4769    int numberColumns = getNumCols();
4770    int numberOtherColumns = presolvedModel->getNumCols();
4771    //bestSolution_ = new double[numberColumns];
4772    // set up map
4773    int * back = new int[numberColumns];
4774    int i;
4775    for (i=0;i<numberColumns;i++)
4776      back[i]=-1;
4777    for (i=0;i<numberOtherColumns;i++)
4778      back[presolvedModel->originalColumns_[i]]=i;
4779    int iColumn;
4780    // set ones in presolved model to values
4781    double * otherSolution = presolvedModel->bestSolution_;
4782    //const double * lower = getColLower();
4783    for (i=0;i<numberIntegers_;i++) {
4784      iColumn = integerVariable_[i];
4785      int jColumn = back[iColumn];
4786      //bestSolution_[iColumn]=lower[iColumn];
4787      if (jColumn >= 0) {
4788        double value=floor(otherSolution[jColumn]+0.5);
4789        solver_->setColLower(iColumn,value);
4790        solver_->setColUpper(iColumn,value);
4791        //bestSolution_[iColumn]=value;
4792      }
4793    }
4794    delete [] back;
4795#if 0
4796    // ** looks as if presolve needs more intelligence
4797    // do presolve - for now just clp but easy to get osi interface
4798    OsiClpSolverInterface * clpSolver
4799      = dynamic_cast<OsiClpSolverInterface *> (solver_);
4800    assert (clpSolver);
4801    ClpSimplex * clp = clpSolver->getModelPtr();
4802    Presolve pinfo;
4803    ClpSimplex * model2 = pinfo.presolvedModel(*clp,1.0e-8);
4804    model2->primal(1);
4805    pinfo.postsolve(true);
4806    const double * solution = solver_->getColSolution();
4807    for (i=0;i<numberIntegers_;i++) {
4808      iColumn = integerVariable_[i];
4809      double value=floor(solution[iColumn]+0.5);
4810      solver_->setColLower(iColumn,value);
4811      solver_->setColUpper(iColumn,value);
4812    }
4813#else
4814    if (!weak) {
4815      // for now give up
4816      int save = numberCutGenerators_;
4817      numberCutGenerators_=0;
4818      bestObjective_=1.0e100;
4819      branchAndBound();
4820      numberCutGenerators_=save;
4821    }
4822#endif
4823    if (bestSolution_) {
4824      // solve problem
4825      resolve();
4826      // should be feasible
4827      int numberIntegerInfeasibilities;
4828      int numberObjectInfeasibilities;
4829      if (!currentSolution_)
4830        currentSolution_ = new double[numberColumns] ;
4831      assert(feasibleSolution(numberIntegerInfeasibilities,
4832                              numberObjectInfeasibilities));
4833    }
4834  } else {
4835    bestSolution_=NULL;
4836  }
4837  numberSolutions_=presolvedModel->numberSolutions_;
4838  numberHeuristicSolutions_=presolvedModel->numberHeuristicSolutions_;
4839  numberNodes_ = presolvedModel->numberNodes_;
4840  numberIterations_ = presolvedModel->numberIterations_;
4841  status_ = presolvedModel->status_;
4842  synchronizeModel();
4843} 
4844// Pass in Message handler (not deleted at end)
4845void 
4846CbcModel::passInMessageHandler(CoinMessageHandler * handler)
4847{
4848  if (defaultHandler_) {
4849    delete handler_;
4850    handler_ = NULL;
4851  }
4852  defaultHandler_=false;
4853  handler_=handler;
4854}
4855void 
4856CbcModel::passInTreeHandler(CbcTree & tree)
4857{
4858  delete tree_;
4859  tree_ = tree.clone();
4860}
4861// Make sure region there
4862void 
4863CbcModel::reserveCurrentSolution()
4864{
4865  int numberColumns = getNumCols() ;
4866  if (!currentSolution_)
4867    currentSolution_ = new double[numberColumns] ;
4868}
4869/* For passing in an CbcModel to do a sub Tree (with derived tree handlers).
4870   Passed in model must exist for duration of branch and bound
4871*/
4872void 
4873CbcModel::passInSubTreeModel(CbcModel & model)
4874{
4875  subTreeModel_=&model;
4876}
4877// For retrieving a copy of subtree model with given OsiSolver or NULL
4878CbcModel * 
4879CbcModel::subTreeModel(OsiSolverInterface * solver) const
4880{
4881  const CbcModel * subModel=subTreeModel_;
4882  if (!subModel)
4883    subModel=this;
4884  // Get new copy
4885  CbcModel * newModel = new CbcModel(*subModel);
4886  if (solver)
4887    newModel->assignSolver(solver);
4888  return newModel;
4889}
4890//#############################################################################
4891// Set/Get Application Data
4892// This is a pointer that the application can store into and retrieve
4893// from the solverInterface.
4894// This field is the application to optionally define and use.
4895//#############################################################################
4896
4897void CbcModel::setApplicationData(void * appData)
4898{
4899  appData_ = appData;
4900}
4901//-----------------------------------------------------------------------------
4902void * CbcModel::getApplicationData() const
4903{
4904  return appData_;
4905}
4906/*  create a submodel from partially fixed problem
4907
4908The method creates a new clean model with given bounds.
4909*/
4910CbcModel * 
4911CbcModel::cleanModel(const double * lower, const double * upper)
4912{
4913  OsiSolverInterface * solver = continuousSolver_->clone();
4914
4915  int numberIntegers = numberIntegers_;
4916  const int * integerVariable = integerVariable_;
4917 
4918  int i;
4919  for (i=0;i<numberIntegers;i++) {
4920    int iColumn=integerVariable[i];
4921    const CbcObject * object = object_[i];
4922    const CbcSimpleInteger * integerObject = 
4923      dynamic_cast<const  CbcSimpleInteger *> (object);
4924    assert(integerObject);
4925    // get original bounds
4926    double originalLower = integerObject->originalLowerBound();
4927    double originalUpper = integerObject->originalUpperBound();
4928    solver->setColLower(iColumn,CoinMax(lower[iColumn],originalLower));
4929    solver->setColUpper(iColumn,CoinMin(upper[iColumn],originalUpper));
4930  }
4931  CbcModel * model = new CbcModel(*solver);
4932  // off some messages
4933  if (handler_->logLevel()<=1) {
4934    model->messagesPointer()->setDetailMessage(3,9);
4935    model->messagesPointer()->setDetailMessage(3,6);
4936    model->messagesPointer()->setDetailMessage(3,4);
4937    model->messagesPointer()->setDetailMessage(3,1);
4938    model->messagesPointer()->setDetailMessage(3,13);
4939    model->messagesPointer()->setDetailMessage(3,14);
4940    model->messagesPointer()->setDetailMessage(3,3007);
4941  }
4942  // Cuts
4943  for ( i = 0;i<numberCutGenerators_;i++) {
4944    int howOften = generator_[i]->howOftenInSub();
4945    if (howOften>-100) {
4946      CbcCutGenerator * generator = virginGenerator_[i];
4947      CglCutGenerator * cglGenerator = generator->generator();
4948      model->addCutGenerator(cglGenerator,howOften,
4949                              generator->cutGeneratorName(),
4950                              generator->normal(),
4951                              generator->atSolution(),
4952                              generator->whenInfeasible(),
4953                              -100, generator->whatDepthInSub(),-1);
4954    }
4955  }
4956  double cutoff = getCutoff();
4957  model->setCutoff(cutoff);
4958  return model;
4959}
4960/* Invoke the branch & cut algorithm on partially fixed problem
4961   
4962   The method uses a subModel created by cleanModel. The search
4963   ends when the tree is exhausted or maximum nodes is reached.
4964
4965   If better solution found then it is saved.
4966   
4967   Returns 0 if search completed and solution, 1 if not completed and solution,
4968   2 if completed and no solution, 3 if not completed and no solution.
4969   
4970   Normally okay to do subModel immediately followed by subBranchandBound
4971   (== other form of subBranchAndBound)
4972   but may need to get at model for advanced features.
4973   
4974   Deletes model
4975   
4976*/
4977 
4978int 
4979CbcModel::subBranchAndBound(CbcModel * model,
4980                            CbcModel * presolvedModel,
4981                            int maximumNodes)
4982{
4983  int i;
4984  double cutoff=model->getCutoff();
4985  CbcModel * model2;
4986  if (presolvedModel) 
4987    model2=presolvedModel;
4988  else
4989    model2=model;
4990  // Do complete search
4991 
4992  for (i=0;i<numberHeuristics_;i++) {
4993    model2->addHeuristic(heuristic_[i]);
4994    model2->heuristic(i)->resetModel(model2);
4995  }
4996  // Definition of node choice
4997  model2->setNodeComparison(nodeCompare_->clone());
4998  //model2->solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
4999  model2->messageHandler()->setLogLevel(CoinMax(0,handler_->logLevel()-1));
5000  //model2->solver()->messageHandler()->setLogLevel(2);
5001  model2->setMaximumCutPassesAtRoot(maximumCutPassesAtRoot_);
5002  model2->setPrintFrequency(50);
5003  model2->setIntParam(CbcModel::CbcMaxNumNode,maximumNodes);
5004  model2->branchAndBound();
5005  delete model2->nodeComparison();
5006  if (model2->getMinimizationObjValue()>cutoff) {
5007    // no good
5008    if (model!=model2)
5009      delete model2;
5010    delete model;
5011    return 2;
5012  }
5013  if (model!=model2) {
5014    // get back solution
5015    model->originalModel(model2,false);
5016    delete model2;
5017  }
5018  int status;
5019  if (model->getMinimizationObjValue()<cutoff&&model->bestSolution()) {
5020    double objValue = model->getObjValue();
5021    const double * solution = model->bestSolution();
5022    setBestSolution(CBC_TREE_SOL,objValue,solution);
5023    status = 0;
5024  } else {
5025    status=2;
5026  }
5027  if (model->status())
5028    status ++ ; // not finished search
5029  delete model;
5030  return status;
5031}
5032/* Invoke the branch & cut algorithm on partially fixed problem
5033   
5034The method creates a new model with given bounds, presolves it
5035then proceeds to explore the branch & cut search tree. The search
5036ends when the tree is exhausted or maximum nodes is reached.
5037Returns 0 if search completed and solution, 1 if not completed and solution,
50382 if completed and no solution, 3 if not completed and no solution.
5039*/
5040int 
5041CbcModel::subBranchAndBound(const double * lower, const double * upper,
5042                            int maximumNodes)
5043{
5044  OsiSolverInterface * solver = continuousSolver_->clone();
5045
5046  int numberIntegers = numberIntegers_;
5047  const int * integerVariable = integerVariable_;
5048 
5049  int i;
5050  for (i=0;i<numberIntegers;i++) {
5051    int iColumn=integerVariable[i];
5052    const CbcObject * object = object_[i];
5053    const CbcSimpleInteger * integerObject = 
5054      dynamic_cast<const  CbcSimpleInteger *> (object);
5055    assert(integerObject);
5056    // get original bounds
5057    double originalLower = integerObject->originalLowerBound();
5058    double originalUpper = integerObject->originalUpperBound();
5059    solver->setColLower(iColumn,CoinMax(lower[iColumn],originalLower));
5060    solver->setColUpper(iColumn,CoinMin(upper[iColumn],originalUpper));
5061  }
5062  CbcModel model(*solver);
5063  // off some messages
5064  if (handler_->logLevel()<=1) {
5065    model.messagesPointer()->setDetailMessage(3,9);
5066    model.messagesPointer()->setDetailMessage(3,6);
5067    model.messagesPointer()->setDetailMessage(3,4);
5068    model.messagesPointer()->setDetailMessage(3,1);
5069    model.messagesPointer()->setDetailMessage(3,3007);
5070  }
5071  double cutoff = getCutoff();
5072  model.setCutoff(cutoff);
5073  // integer presolve
5074  CbcModel * model2 = model.integerPresolve(false);
5075  if (!model2||!model2->getNumRows()) {
5076    delete model2;
5077    delete solver;
5078    return 2;
5079  }
5080  if (handler_->logLevel()>1)
5081    printf("Reduced model has %d rows and %d columns\n",
5082           model2->getNumRows(),model2->getNumCols());
5083  // Do complete search
5084 
5085  // Cuts
5086  for ( i = 0;i<numberCutGenerators_;i++) {
5087    int howOften = generator_[i]->howOftenInSub();
5088    if (howOften>-100) {
5089      CbcCutGenerator * generator = virginGenerator_[i];
5090      CglCutGenerator * cglGenerator = generator->generator();
5091      model2->addCutGenerator(cglGenerator,howOften,
5092                              generator->cutGeneratorName(),
5093                              generator->normal(),
5094                              generator->atSolution(),
5095                              generator->whenInfeasible(),
5096                              -100, generator->whatDepthInSub(),-1);
5097    }
5098  }
5099  for (i=0;i<numberHeuristics_;i++) {
5100    model2->addHeuristic(heuristic_[i]);
5101    model2->heuristic(i)->resetModel(model2);
5102  }
5103  // Definition of node choice
5104  model2->setNodeComparison(nodeCompare_->clone());
5105  //model2->solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
5106  model2->messageHandler()->setLogLevel(CoinMax(0,handler_->logLevel()-1));
5107  //model2->solver()->messageHandler()->setLogLevel(2);
5108  model2->setMaximumCutPassesAtRoot(maximumCutPassesAtRoot_);
5109  model2->setPrintFrequency(50);
5110  model2->setIntParam(CbcModel::CbcMaxNumNode,maximumNodes);
5111  model2->branchAndBound();
5112  delete model2->nodeComparison();
5113  if (model2->getMinimizationObjValue()>cutoff) {
5114    // no good
5115    delete model2;
5116    delete solver;
5117    return 2;
5118  }
5119  // get back solution
5120  model.originalModel(model2,false);
5121  delete model2;
5122  int status;
5123  if (model.getMinimizationObjValue()<cutoff&&model.bestSolution()) {
5124    double objValue = model.getObjValue();
5125    const double * solution = model.bestSolution();
5126    setBestSolution(CBC_TREE_SOL,objValue,solution);
5127    status = 0;
5128  } else {
5129    status=2;
5130  }
5131  if (model.status())
5132    status ++ ; // not finished search
5133  delete solver;
5134  return status;
5135}
5136// Set a pointer to a row cut which will be added instead of normal branching.
5137void 
5138CbcModel::setNextRowCut(const OsiRowCut & cut)
5139{ 
5140  nextRowCut_=new OsiRowCut(cut);
5141}
5142/* Process root node and return a strengthened model
5143   
5144The method assumes that initialSolve() has been called to solve the
5145LP relaxation. It processes the root node and then returns a pointer
5146to the strengthened model (or NULL if infeasible)
5147*/
5148OsiSolverInterface * 
5149CbcModel::strengthenedModel()
5150{
5151/*
5152  Assume we're done, and see if we're proven wrong.
5153*/
5154/*
5155  Scan the variables, noting the integer variables. Create an
5156  CbcSimpleInteger object for each integer variable.
5157*/
5158  findIntegers(false) ;
5159/*
5160  Ensure that objects on the lists of CbcObjects, heuristics, and cut
5161  generators attached to this model all refer to this model.
5162*/
5163  synchronizeModel() ;
5164
5165  // Set so we can tell we are in initial phase in resolve
5166  continuousObjective_ = -COIN_DBL_MAX ;
5167/*
5168  Solve the relaxation.
5169
5170  Apparently there are circumstances where this will be non-trivial --- i.e.,
5171  we've done something since initialSolve that's trashed the solution to the
5172  continuous relaxation.
5173*/
5174  bool feasible = resolve() ;
5175/*
5176  If the linear relaxation of the root is infeasible, bail out now. Otherwise,
5177  continue with processing the root node.
5178*/
5179  if (!feasible)
5180  { handler_->message(CBC_INFEAS,messages_)<< CoinMessageEol ;
5181    return NULL; }
5182  // Save objective (just so user can access it)
5183  originalContinuousObjective_ = solver_->getObjValue();
5184
5185/*
5186  Begin setup to process a feasible root node.
5187*/
5188  bestObjective_ = CoinMin(bestObjective_,1.0e50) ;
5189  numberSolutions_ = 0 ;
5190  numberHeuristicSolutions_ = 0 ;
5191  // Everything is minimization
5192  double cutoff=getCutoff() ;
5193  double direction = solver_->getObjSense() ;
5194  if (cutoff < 1.0e20&&direction<0.0)
5195    messageHandler()->message(CBC_CUTOFF_WARNING1,
5196                                    messages())
5197                                      << cutoff << -cutoff << CoinMessageEol ;
5198  if (cutoff > bestObjective_)
5199    cutoff = bestObjective_ ;
5200  setCutoff(cutoff) ;
5201/*
5202  We probably already have a current solution, but just in case ...
5203*/
5204  int numberColumns = getNumCols() ;
5205  if (!currentSolution_)
5206    currentSolution_ = new double[numberColumns] ;
5207/*
5208  Create a copy of the solver, thus capturing the original (root node)
5209  constraint system (aka the continuous system).
5210*/
5211  continuousSolver_ = solver_->clone() ;
5212  numberRowsAtContinuous_ = getNumRows() ;
5213/*
5214  Check the objective to see if we can deduce a nontrivial increment. If
5215  it's better than the current value for CbcCutoffIncrement, it'll be
5216  installed.
5217*/
5218  analyzeObjective() ;
5219/*
5220  Set up for cut generation. addedCuts_ holds the cuts which are relevant for
5221  the active subproblem. whichGenerator will be used to record the generator
5222  that produced a given cut.
5223*/
5224  int maximumWhich = 1000 ;
5225  int * whichGenerator = new int[maximumWhich] ;
5226  maximumNumberCuts_ = 0 ;
5227  currentNumberCuts_ = 0 ;
5228  delete [] addedCuts_ ;
5229  addedCuts_ = NULL ;
5230  /* 
5231  Generate cuts at the root node and reoptimise. solveWithCuts does the heavy
5232  lifting. It will iterate a generate/reoptimise loop (including reduced cost
5233  fixing) until no cuts are generated, the change in objective falls off,  or
5234  the limit on the number of rounds of cut generation is exceeded.
5235
5236  At the end of all this, any cuts will be recorded in cuts and also
5237  installed in the solver's constraint system. We'll have reoptimised, and
5238  removed any slack cuts (numberOldActiveCuts and numberNewCuts have been
5239  adjusted accordingly).
5240
5241  Tell cut generators they can be a bit more aggressive at root node
5242
5243*/
5244  int iCutGenerator;
5245  for (iCutGenerator = 0;iCutGenerator<numberCutGenerators_;iCutGenerator++) {
5246    CglCutGenerator * generator = generator_[iCutGenerator]->generator();
5247    generator->setAggressiveness(generator->getAggressiveness()+100);
5248  }
5249  OsiCuts cuts ;
5250  int numberOldActiveCuts = 0 ;
5251  int numberNewCuts = 0 ;
5252  { int iObject ;
5253    int preferredWay ;
5254    int numberUnsatisfied = 0 ;
5255    memcpy(currentSolution_,solver_->getColSolution(),
5256           numberColumns*sizeof(double)) ;
5257
5258    for (iObject = 0 ; iObject < numberObjects_ ; iObject++)
5259    { double infeasibility =
5260          object_[iObject]->infeasibility(preferredWay) ;
5261      if (infeasibility) numberUnsatisfied++ ; }
5262    if (numberUnsatisfied)
5263    { feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,
5264                               NULL,numberOldActiveCuts,numberNewCuts,
5265                               maximumWhich, whichGenerator) ; } }
5266/*
5267  We've taken the continuous relaxation as far as we can.
5268*/
5269
5270  OsiSolverInterface * newSolver=NULL;
5271  if (feasible) {
5272    // make copy of current solver
5273    newSolver = solver_->clone();
5274  }
5275/*
5276  Clean up dangling objects. continuousSolver_ may already be toast.
5277*/
5278  delete [] whichGenerator ;
5279  delete [] walkback_ ;
5280  walkback_ = NULL ;
5281  delete [] addedCuts_ ;
5282  addedCuts_ = NULL ;
5283  if (continuousSolver_)
5284  { delete continuousSolver_ ;
5285    continuousSolver_ = NULL ; }
5286/*
5287  Destroy global cuts by replacing with an empty OsiCuts object.
5288*/
5289  globalCuts_= OsiCuts() ;
5290 
5291  return newSolver; 
5292}
5293// Just update objectiveValue
5294void CbcModel::setBestObjectiveValue( double objectiveValue)
5295{
5296  bestObjective_=objectiveValue;
5297}
5298double 
5299CbcModel::getBestPossibleObjValue() const
5300{ 
5301  return CoinMin(bestPossibleObjective_,bestObjective_) * solver_->getObjSense() ;
5302}
Note: See TracBrowser for help on using the repository browser.