source: trunk/Cbc/src/CbcHeuristic.cpp

Last change on this file was 2645, checked in by samuelbrito, 3 months ago

conflict based preprocessing and cut separators

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 110.5 KB
Line 
1/* $Id: CbcHeuristic.cpp 2645 2019-07-30 18:32:57Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#pragma warning(disable : 4786)
9#endif
10
11#include "CbcConfig.h"
12
13#include <cassert>
14#include <cstdlib>
15#include <cmath>
16#include <cfloat>
17
18//#define PRINT_DEBUG
19#ifdef COIN_HAS_CLP
20#include "OsiClpSolverInterface.hpp"
21#endif
22#include "CbcModel.hpp"
23#include "CbcMessage.hpp"
24#include "CbcHeuristic.hpp"
25#include "CbcHeuristicFPump.hpp"
26#include "CbcHeuristicRINS.hpp"
27#include "CbcEventHandler.hpp"
28#include "CbcStrategy.hpp"
29#include "CglPreProcess.hpp"
30#include "CglGomory.hpp"
31#include "CglProbing.hpp"
32#include "OsiAuxInfo.hpp"
33#include "OsiRowCutDebugger.hpp"
34#include "OsiPresolve.hpp"
35#include "CbcBranchActual.hpp"
36#include "CbcCutGenerator.hpp"
37#include "CoinMpsIO.hpp"
38//==============================================================================
39
40CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode &rhs)
41{
42  numObjects_ = rhs.numObjects_;
43  brObj_ = new CbcBranchingObject *[numObjects_];
44  for (int i = 0; i < numObjects_; ++i) {
45    brObj_[i] = rhs.brObj_[i]->clone();
46  }
47}
48
49void CbcHeuristicNodeList::gutsOfDelete()
50{
51  for (int i = (static_cast< int >(nodes_.size())) - 1; i >= 0; --i) {
52    delete nodes_[i];
53  }
54}
55
56void CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList &rhs)
57{
58  append(rhs);
59}
60
61CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList &rhs)
62{
63  gutsOfCopy(rhs);
64}
65
66CbcHeuristicNodeList &CbcHeuristicNodeList::operator=(const CbcHeuristicNodeList &rhs)
67{
68  if (this != &rhs) {
69    gutsOfDelete();
70    gutsOfCopy(rhs);
71  }
72  return *this;
73}
74
75CbcHeuristicNodeList::~CbcHeuristicNodeList()
76{
77  gutsOfDelete();
78}
79
80void CbcHeuristicNodeList::append(CbcHeuristicNode *&node)
81{
82  nodes_.push_back(node);
83  node = NULL;
84}
85
86void CbcHeuristicNodeList::append(const CbcHeuristicNodeList &nodes)
87{
88  nodes_.reserve(nodes_.size() + nodes.size());
89  for (int i = 0; i < nodes.size(); ++i) {
90    CbcHeuristicNode *node = new CbcHeuristicNode(*nodes.node(i));
91    append(node);
92  }
93}
94
95//==============================================================================
96#define DEFAULT_WHERE ((255 - 2 - 16) * (1 + 256))
97// Default Constructor
98CbcHeuristic::CbcHeuristic()
99  : model_(NULL)
100  , when_(2)
101  , numberNodes_(200)
102  , feasibilityPumpOptions_(-1)
103  , fractionSmall_(1.0)
104  , heuristicName_("Unknown")
105  , howOften_(1)
106  , decayFactor_(0.0)
107  , switches_(0)
108  , whereFrom_(DEFAULT_WHERE)
109  , shallowDepth_(1)
110  , howOftenShallow_(1)
111  , numInvocationsInShallow_(0)
112  , numInvocationsInDeep_(0)
113  , lastRunDeep_(0)
114  , numRuns_(0)
115  , minDistanceToRun_(1)
116  , runNodes_()
117  , numCouldRun_(0)
118  , numberSolutionsFound_(0)
119  , numberNodesDone_(0)
120  , inputSolution_(NULL)
121{
122  // As CbcHeuristic virtual need to modify .cpp if above change
123}
124
125// Constructor from model
126CbcHeuristic::CbcHeuristic(CbcModel &model)
127  : model_(&model)
128  , when_(2)
129  , numberNodes_(200)
130  , feasibilityPumpOptions_(-1)
131  , fractionSmall_(1.0)
132  , heuristicName_("Unknown")
133  , howOften_(1)
134  , decayFactor_(0.0)
135  , switches_(0)
136  , whereFrom_(DEFAULT_WHERE)
137  , shallowDepth_(1)
138  , howOftenShallow_(1)
139  , numInvocationsInShallow_(0)
140  , numInvocationsInDeep_(0)
141  , lastRunDeep_(0)
142  , numRuns_(0)
143  , minDistanceToRun_(1)
144  , runNodes_()
145  , numCouldRun_(0)
146  , numberSolutionsFound_(0)
147  , numberNodesDone_(0)
148  , inputSolution_(NULL)
149{
150}
151
152void CbcHeuristic::gutsOfCopy(const CbcHeuristic &rhs)
153{
154  model_ = rhs.model_;
155  when_ = rhs.when_;
156  numberNodes_ = rhs.numberNodes_;
157  feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
158  fractionSmall_ = rhs.fractionSmall_;
159  randomNumberGenerator_ = rhs.randomNumberGenerator_;
160  heuristicName_ = rhs.heuristicName_;
161  howOften_ = rhs.howOften_;
162  decayFactor_ = rhs.decayFactor_;
163  switches_ = rhs.switches_;
164  whereFrom_ = rhs.whereFrom_;
165  shallowDepth_ = rhs.shallowDepth_;
166  howOftenShallow_ = rhs.howOftenShallow_;
167  numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
168  numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
169  lastRunDeep_ = rhs.lastRunDeep_;
170  numRuns_ = rhs.numRuns_;
171  numCouldRun_ = rhs.numCouldRun_;
172  minDistanceToRun_ = rhs.minDistanceToRun_;
173  runNodes_ = rhs.runNodes_;
174  numberSolutionsFound_ = rhs.numberSolutionsFound_;
175  numberNodesDone_ = rhs.numberNodesDone_;
176  if (rhs.inputSolution_) {
177    int numberColumns = model_->getNumCols();
178    setInputSolution(rhs.inputSolution_, rhs.inputSolution_[numberColumns]);
179  }
180}
181// Copy constructor
182CbcHeuristic::CbcHeuristic(const CbcHeuristic &rhs)
183{
184  inputSolution_ = NULL;
185  gutsOfCopy(rhs);
186}
187
188// Assignment operator
189CbcHeuristic &
190CbcHeuristic::operator=(const CbcHeuristic &rhs)
191{
192  if (this != &rhs) {
193    gutsOfDelete();
194    gutsOfCopy(rhs);
195  }
196  return *this;
197}
198
199void CbcHeurDebugNodes(CbcModel *model_)
200{
201  CbcNode *node = model_->currentNode();
202  CbcNodeInfo *nodeInfo = node->nodeInfo();
203  std::cout << "===============================================================\n";
204  while (nodeInfo) {
205    const CbcNode *node = nodeInfo->owner();
206    printf("nodeinfo: node %i\n", nodeInfo->nodeNumber());
207    {
208      const CbcIntegerBranchingObject *brPrint = dynamic_cast< const CbcIntegerBranchingObject * >(nodeInfo->parentBranch());
209      if (!brPrint) {
210        printf("    parentBranch: NULL\n");
211      } else {
212        const double *downBounds = brPrint->downBounds();
213        const double *upBounds = brPrint->upBounds();
214        int variable = brPrint->variable();
215        int way = brPrint->way();
216        printf("   parentBranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
217          variable, static_cast< int >(downBounds[0]), static_cast< int >(downBounds[1]),
218          static_cast< int >(upBounds[0]), static_cast< int >(upBounds[1]), way);
219      }
220    }
221    if (!node) {
222      printf("    owner: NULL\n");
223    } else {
224      printf("    owner: node %i depth %i onTree %i active %i",
225        node->nodeNumber(), node->depth(), node->onTree(), node->active());
226      const OsiBranchingObject *osibr = nodeInfo->owner()->branchingObject();
227      const CbcBranchingObject *cbcbr = dynamic_cast< const CbcBranchingObject * >(osibr);
228      const CbcIntegerBranchingObject *brPrint = dynamic_cast< const CbcIntegerBranchingObject * >(cbcbr);
229      if (!brPrint) {
230        printf("        ownerBranch: NULL\n");
231      } else {
232        const double *downBounds = brPrint->downBounds();
233        const double *upBounds = brPrint->upBounds();
234        int variable = brPrint->variable();
235        int way = brPrint->way();
236        printf("        ownerbranch: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
237          variable, static_cast< int >(downBounds[0]), static_cast< int >(downBounds[1]),
238          static_cast< int >(upBounds[0]), static_cast< int >(upBounds[1]), way);
239      }
240    }
241    nodeInfo = nodeInfo->parent();
242  }
243}
244
245void CbcHeuristic::debugNodes()
246{
247  CbcHeurDebugNodes(model_);
248}
249
250void CbcHeuristic::printDistanceToNodes()
251{
252  const CbcNode *currentNode = model_->currentNode();
253  if (currentNode != NULL) {
254    CbcHeuristicNode *nodeDesc = new CbcHeuristicNode(*model_);
255    for (int i = runNodes_.size() - 1; i >= 0; --i) {
256      nodeDesc->distance(runNodes_.node(i));
257    }
258    runNodes_.append(nodeDesc);
259  }
260}
261
262bool CbcHeuristic::shouldHeurRun(int whereFrom)
263{
264  assert(whereFrom >= 0 && whereFrom < 16);
265  // take off 8 (code - likes new solution)
266  whereFrom &= 7;
267  if ((whereFrom_ & (1 << whereFrom)) == 0)
268    return false;
269    // No longer used for original purpose - so use for ever run at all JJF
270#ifndef JJF_ONE
271  // Don't run if hot start or no rows!
272  if (model_ && (model_->hotstartSolution() || !model_->getNumRows()))
273    return false;
274  else
275    return true;
276#else
277#ifdef JJF_ZERO
278  const CbcNode *currentNode = model_->currentNode();
279  if (currentNode == NULL) {
280    return false;
281  }
282
283  debugNodes();
284  //   return false;
285
286  const int depth = currentNode->depth();
287#else
288  int depth = model_->currentDepth();
289#endif
290
291  const int nodeCount = model_->getNodeCount(); // FIXME: check that this is
292  // correct in parallel
293
294  if (nodeCount == 0 || depth <= shallowDepth_) {
295    // what to do when we are in the shallow part of the tree
296    if (model_->getCurrentPassNumber() == 1) {
297      // first time in the node...
298      numInvocationsInShallow_ = 0;
299    }
300    ++numInvocationsInShallow_;
301    // Very large howOftenShallow_ will give the original test:
302    // (model_->getCurrentPassNumber() != 1)
303    //    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
304    if ((numInvocationsInShallow_ % howOftenShallow_) != 0) {
305      return false;
306    }
307    // LL: should we save these nodes in the list of nodes where the heur was
308    // LL: run?
309#ifndef JJF_ONE
310    if (currentNode != NULL) {
311      // Get where we are and create the appropriate CbcHeuristicNode object
312      CbcHeuristicNode *nodeDesc = new CbcHeuristicNode(*model_);
313      runNodes_.append(nodeDesc);
314    }
315#endif
316  } else {
317    // deeper in the tree
318    if (model_->getCurrentPassNumber() == 1) {
319      // first time in the node...
320      ++numInvocationsInDeep_;
321    }
322    if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
323      return false;
324    }
325    if (model_->getCurrentPassNumber() > 1) {
326      // Run the heuristic only when first entering the node.
327      // LL: I don't think this is right. It should run just before strong
328      // LL: branching, I believe.
329      return false;
330    }
331    // Get where we are and create the appropriate CbcHeuristicNode object
332    CbcHeuristicNode *nodeDesc = new CbcHeuristicNode(*model_);
333    //#ifdef PRINT_DEBUG
334#ifndef JJF_ONE
335    const double minDistanceToRun = 1.5 * log((double)depth) / log((double)2);
336#else
337    const double minDistanceToRun = minDistanceToRun_;
338#endif
339#ifdef PRINT_DEBUG
340    double minDistance = nodeDesc->minDistance(runNodes_);
341    std::cout << "minDistance = " << minDistance
342              << ", minDistanceToRun = " << minDistanceToRun << std::endl;
343#endif
344    if (nodeDesc->minDistanceIsSmall(runNodes_, minDistanceToRun)) {
345      delete nodeDesc;
346      return false;
347    }
348    runNodes_.append(nodeDesc);
349    lastRunDeep_ = numInvocationsInDeep_;
350    //    ++lastRunDeep_;
351  }
352  ++numRuns_;
353  return true;
354#endif
355}
356
357bool CbcHeuristic::shouldHeurRun_randomChoice()
358{
359  if (!when_)
360    return false;
361  int depth = model_->currentDepth();
362  // when_ -999 is special marker to force to run
363  if (depth != 0 && when_ != -999) {
364    const double numerator = depth * depth;
365    const double denominator = exp(depth * log(2.0));
366    double probability = numerator / denominator;
367    double randomNumber = randomNumberGenerator_.randomDouble();
368    int when = when_ % 100;
369    if (when > 2 && when < 8) {
370      /* JJF adjustments
371            3 only at root and if no solution
372            4 only at root and if this heuristic has not got solution
373            5 decay (but only if no solution)
374            6 if depth <3 or decay
375            7 run up to 2 times if solution found 4 otherwise
376            */
377      switch (when) {
378      case 3:
379      default:
380        if (model_->bestSolution())
381          probability = -1.0;
382        break;
383      case 4:
384        if (numberSolutionsFound_)
385          probability = -1.0;
386        break;
387      case 5:
388        assert(decayFactor_);
389        if (model_->bestSolution()) {
390          probability = -1.0;
391        } else if (numCouldRun_ > 1000) {
392          decayFactor_ *= 0.99;
393          probability *= decayFactor_;
394        }
395        break;
396      case 6:
397        if (depth >= 3) {
398          if ((numCouldRun_ % howOften_) == 0 && numberSolutionsFound_ * howOften_ < numCouldRun_) {
399            //#define COIN_DEVELOP
400#ifdef COIN_DEVELOP
401            int old = howOften_;
402#endif
403            howOften_ = CoinMin(CoinMax(static_cast< int >(howOften_ * 1.1), howOften_ + 1), 1000000);
404#ifdef COIN_DEVELOP
405            printf("Howoften changed from %d to %d for %s\n",
406              old, howOften_, heuristicName_.c_str());
407#endif
408          }
409          probability = 1.0 / howOften_;
410          if (model_->bestSolution())
411            probability *= 0.5;
412        } else {
413          probability = 1.1;
414        }
415        break;
416      case 7:
417        if ((model_->bestSolution() && numRuns_ >= 2) || numRuns_ >= 4)
418          probability = -1.0;
419        break;
420      }
421    }
422    if (randomNumber > probability)
423      return false;
424
425    if (model_->getCurrentPassNumber() > 1)
426      return false;
427#ifdef COIN_DEVELOP
428    printf("Running %s, random %g probability %g\n",
429      heuristicName_.c_str(), randomNumber, probability);
430#endif
431  } else {
432#ifdef COIN_DEVELOP
433    printf("Running %s, depth %d when %d\n",
434      heuristicName_.c_str(), depth, when_);
435#endif
436  }
437  ++numRuns_;
438  return true;
439}
440
441// Resets stuff if model changes
442void CbcHeuristic::resetModel(CbcModel *model)
443{
444  model_ = model;
445}
446// Set seed
447void CbcHeuristic::setSeed(int value)
448{
449  if (value == 0) {
450    double time = fabs(CoinGetTimeOfDay());
451    while (time >= COIN_INT_MAX)
452      time *= 0.5;
453    value = static_cast< int >(time);
454    char printArray[100];
455    sprintf(printArray, "using time of day seed was changed from %d to %d",
456      randomNumberGenerator_.getSeed(), value);
457    if (model_)
458      model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
459        << printArray
460        << CoinMessageEol;
461  }
462  randomNumberGenerator_.setSeed(value);
463}
464// Get seed
465int CbcHeuristic::getSeed() const
466{
467  return randomNumberGenerator_.getSeed();
468}
469
470// Create C++ lines to get to current state
471void CbcHeuristic::generateCpp(FILE *fp, const char *heuristic)
472{
473  // hard coded as CbcHeuristic virtual
474  if (when_ != 2)
475    fprintf(fp, "3  %s.setWhen(%d);\n", heuristic, when_);
476  else
477    fprintf(fp, "4  %s.setWhen(%d);\n", heuristic, when_);
478  if (numberNodes_ != 200)
479    fprintf(fp, "3  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
480  else
481    fprintf(fp, "4  %s.setNumberNodes(%d);\n", heuristic, numberNodes_);
482  if (feasibilityPumpOptions_ != -1)
483    fprintf(fp, "3  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
484  else
485    fprintf(fp, "4  %s.setFeasibilityPumpOptions(%d);\n", heuristic, feasibilityPumpOptions_);
486  if (fractionSmall_ != 1.0)
487    fprintf(fp, "3  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
488  else
489    fprintf(fp, "4  %s.setFractionSmall(%g);\n", heuristic, fractionSmall_);
490  if (heuristicName_ != "Unknown")
491    fprintf(fp, "3  %s.setHeuristicName(\"%s\");\n",
492      heuristic, heuristicName_.c_str());
493  else
494    fprintf(fp, "4  %s.setHeuristicName(\"%s\");\n",
495      heuristic, heuristicName_.c_str());
496  if (decayFactor_ != 0.0)
497    fprintf(fp, "3  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
498  else
499    fprintf(fp, "4  %s.setDecayFactor(%g);\n", heuristic, decayFactor_);
500  if (switches_ != 0)
501    fprintf(fp, "3  %s.setSwitches(%d);\n", heuristic, switches_);
502  else
503    fprintf(fp, "4  %s.setSwitches(%d);\n", heuristic, switches_);
504  if (whereFrom_ != DEFAULT_WHERE)
505    fprintf(fp, "3  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
506  else
507    fprintf(fp, "4  %s.setWhereFrom(%d);\n", heuristic, whereFrom_);
508  if (shallowDepth_ != 1)
509    fprintf(fp, "3  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
510  else
511    fprintf(fp, "4  %s.setShallowDepth(%d);\n", heuristic, shallowDepth_);
512  if (howOftenShallow_ != 1)
513    fprintf(fp, "3  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
514  else
515    fprintf(fp, "4  %s.setHowOftenShallow(%d);\n", heuristic, howOftenShallow_);
516  if (minDistanceToRun_ != 1)
517    fprintf(fp, "3  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
518  else
519    fprintf(fp, "4  %s.setMinDistanceToRun(%d);\n", heuristic, minDistanceToRun_);
520}
521// Destructor
522CbcHeuristic::~CbcHeuristic()
523{
524  delete[] inputSolution_;
525}
526
527// update model
528void CbcHeuristic::setModel(CbcModel *model)
529{
530  model_ = model;
531}
532/* Clone but ..
533   type 0 clone solver, 1 clone continuous solver
534   Add 2 to say without integer variables which are at low priority
535   Add 4 to say quite likely infeasible so give up easily.*/
536OsiSolverInterface *
537CbcHeuristic::cloneBut(int type)
538{
539  OsiSolverInterface *solver;
540  if ((type & 1) == 0 || !model_->continuousSolver())
541    solver = model_->solver()->clone();
542  else
543    solver = model_->continuousSolver()->clone();
544#ifdef COIN_HAS_CLP
545  OsiClpSolverInterface *clpSolver
546    = dynamic_cast< OsiClpSolverInterface * >(solver);
547#endif
548  if ((type & 2) != 0) {
549    int n = model_->numberObjects();
550    int priority = model_->continuousPriority();
551    if (priority < COIN_INT_MAX) {
552      for (int i = 0; i < n; i++) {
553        const OsiObject *obj = model_->object(i);
554        const CbcSimpleInteger *thisOne = dynamic_cast< const CbcSimpleInteger * >(obj);
555        if (thisOne) {
556          int iColumn = thisOne->columnNumber();
557          if (thisOne->priority() >= priority)
558            solver->setContinuous(iColumn);
559        }
560      }
561    }
562#ifdef COIN_HAS_CLP
563    if (clpSolver) {
564      for (int i = 0; i < n; i++) {
565        const OsiObject *obj = model_->object(i);
566        const CbcSimpleInteger *thisOne = dynamic_cast< const CbcSimpleInteger * >(obj);
567        if (thisOne) {
568          int iColumn = thisOne->columnNumber();
569          if (clpSolver->isOptionalInteger(iColumn))
570            clpSolver->setContinuous(iColumn);
571        }
572      }
573    }
574#endif
575  }
576#ifdef COIN_HAS_CLP
577  if ((type & 4) != 0 && clpSolver) {
578    int options = clpSolver->getModelPtr()->moreSpecialOptions();
579    clpSolver->getModelPtr()->setMoreSpecialOptions(options | 64);
580  }
581  if (clpSolver) {
582    // take out zero cost integers which will be integer anyway
583    const double *rowLower = clpSolver->getRowLower();
584    const double *rowUpper = clpSolver->getRowUpper();
585    const double *objective = clpSolver->getObjCoefficients();
586    int numberRows = clpSolver->getNumRows();
587    const CoinPackedMatrix *matrixByRow = clpSolver->getMatrixByRow();
588    const double *elementByRow = matrixByRow->getElements();
589    const int *column = matrixByRow->getIndices();
590    const CoinBigIndex *rowStart = matrixByRow->getVectorStarts();
591    const int *rowLength = matrixByRow->getVectorLengths();
592    const CoinPackedMatrix *matrixByColumn = clpSolver->getMatrixByCol();
593    //const double * element = matrixByColumn->getElements();
594    //const int *row = matrixByColumn->getIndices();
595    //const CoinBigIndex *columnStart = matrixByColumn->getVectorStarts();
596    const int *columnLength = matrixByColumn->getVectorLengths();
597    for (int i = 0; i < numberRows; i++) {
598      if (rowLower[i] != floor(rowLower[i]) ||
599          rowUpper[i] != floor(rowUpper[i]))
600        continue;
601      int jColumn = -1;
602      for (CoinBigIndex k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
603        int iColumn = column[k];
604        double value = elementByRow[k];
605        if (!clpSolver->isInteger(iColumn) || floor(value) != value) {
606          jColumn = -2;
607          break;
608        } else if (!objective[iColumn] && columnLength[iColumn] == 1) {
609          jColumn = iColumn;
610        }
611      }
612      if (jColumn>=0)
613        clpSolver->setContinuous(jColumn);
614    }
615  }
616#endif
617  return solver;
618}
619// Whether to exit at once on gap
620bool CbcHeuristic::exitNow(double bestObjective) const
621{
622  if ((switches_ & 2048) != 0) {
623    // exit may be forced - but unset for next time
624    switches_ &= ~2048;
625    if ((switches_ & 1024) != 0)
626      return true;
627  } else if ((switches_ & 1) == 0) {
628    return false;
629  }
630  // See if can stop on gap
631  OsiSolverInterface *solver = model_->solver();
632  double bestPossibleObjective = solver->getObjValue() * solver->getObjSense();
633  double absGap = CoinMax(model_->getAllowableGap(),
634    model_->getHeuristicGap());
635  double fracGap = CoinMax(model_->getAllowableFractionGap(),
636    model_->getHeuristicFractionGap());
637  double testGap = CoinMax(absGap, fracGap * CoinMax(fabs(bestObjective), fabs(bestPossibleObjective)));
638
639  if (bestObjective - bestPossibleObjective < testGap
640    && model_->getCutoffIncrement() >= 0.0) {
641    return true;
642  } else {
643    return false;
644  }
645}
646#ifdef HISTORY_STATISTICS
647extern bool getHistoryStatistics_;
648#endif
649static double sizeRatio(int numberRowsNow, int numberColumnsNow,
650  int numberRowsStart, int numberColumnsStart)
651{
652  double valueNow;
653  if (numberRowsNow * 10 > numberColumnsNow || numberColumnsNow < 200) {
654    valueNow = 2 * numberRowsNow + numberColumnsNow;
655  } else {
656    // long and thin - rows are more important
657    if (numberRowsNow * 40 > numberColumnsNow)
658      valueNow = 10 * numberRowsNow + numberColumnsNow;
659    else
660      valueNow = 200 * numberRowsNow + numberColumnsNow;
661  }
662  double valueStart;
663  if (numberRowsStart * 10 > numberColumnsStart || numberColumnsStart < 200) {
664    valueStart = 2 * numberRowsStart + numberColumnsStart;
665  } else {
666    // long and thin - rows are more important
667    if (numberRowsStart * 40 > numberColumnsStart)
668      valueStart = 10 * numberRowsStart + numberColumnsStart;
669    else
670      valueStart = 200 * numberRowsStart + numberColumnsStart;
671  }
672  //printf("sizeProblem Now %g, %d rows, %d columns\nsizeProblem Start %g, %d rows, %d columns\n",
673  // valueNow,numberRowsNow,numberColumnsNow,
674  // valueStart,numberRowsStart,numberColumnsStart);
675  if (10 * numberRowsNow < 8 * numberRowsStart || 10 * numberColumnsNow < 7 * numberColumnsStart)
676    return valueNow / valueStart;
677  else if (10 * numberRowsNow < 9 * numberRowsStart)
678    return 1.1 * (valueNow / valueStart);
679  else if (numberRowsNow < numberRowsStart)
680    return 1.5 * (valueNow / valueStart);
681  else
682    return 2.0 * (valueNow / valueStart);
683}
684
685//static int saveModel=0;
686// Do mini branch and bound (return 1 if solution)
687int CbcHeuristic::smallBranchAndBound(OsiSolverInterface *solver, int numberNodes,
688  double *newSolution, double &newSolutionValue,
689  double cutoff, std::string name) const
690{
691  CbcEventHandler *eventHandler = model_->getEventHandler();
692  // Use this fraction
693  double fractionSmall = fractionSmall_;
694  int maximumSolutions = model_->getMaximumSolutions();
695  int iterationMultiplier = 100;
696  if (eventHandler) {
697    typedef struct {
698      double fractionSmall;
699      double spareDouble[3];
700      OsiSolverInterface *solver;
701      void *sparePointer[2];
702      int numberNodes;
703      int maximumSolutions;
704      int iterationMultiplier;
705      int howOften;
706      int spareInt[3];
707    } SmallMod;
708    SmallMod temp;
709    temp.solver = solver;
710    temp.fractionSmall = fractionSmall;
711    temp.numberNodes = numberNodes;
712    temp.iterationMultiplier = iterationMultiplier;
713    temp.howOften = howOften_;
714    temp.maximumSolutions = maximumSolutions;
715    CbcEventHandler::CbcAction status = eventHandler->event(CbcEventHandler::smallBranchAndBound,
716      &temp);
717    if (status == CbcEventHandler::killSolution)
718      return -1;
719    if (status == CbcEventHandler::takeAction) {
720      fractionSmall = temp.fractionSmall;
721      numberNodes = temp.numberNodes;
722      iterationMultiplier = temp.iterationMultiplier;
723      howOften_ = temp.howOften;
724      maximumSolutions = temp.maximumSolutions;
725    }
726  }
727#if 0 
728  if (saveModel || model_->getMaximumSolutions()==100) {
729    printf("writing model\n");
730    solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
731  }
732#endif
733  // size before
734  int shiftRows = 0;
735  if (numberNodes < 0)
736    shiftRows = solver->getNumRows() - numberNodes_;
737  int numberRowsStart = solver->getNumRows() - shiftRows;
738  int numberColumnsStart = solver->getNumCols();
739#ifdef CLP_INVESTIGATE
740  printf("%s has %d rows, %d columns\n",
741    name.c_str(), solver->getNumRows(), solver->getNumCols());
742#endif
743  double before = 2 * numberRowsStart + numberColumnsStart;
744  if (before > 40000.0) {
745    // fairly large - be more conservative
746    double multiplier = 1.0 - 0.3 * CoinMin(100000.0, before - 40000.0) / 100000.0;
747    if (multiplier < 1.0) {
748      fractionSmall *= multiplier;
749#ifdef CLP_INVESTIGATE
750      printf("changing fractionSmall from %g to %g for %s\n",
751        fractionSmall_, fractionSmall, name.c_str());
752#endif
753    }
754  }
755#ifdef COIN_HAS_CLP
756  OsiClpSolverInterface *clpSolver = dynamic_cast< OsiClpSolverInterface * >(solver);
757  if (clpSolver && (clpSolver->specialOptions() & 65536) == 0) {
758    // go faster stripes
759    if (clpSolver->getNumRows() < 300 && clpSolver->getNumCols() < 500) {
760      clpSolver->setupForRepeatedUse(2, 0);
761    } else {
762      clpSolver->setupForRepeatedUse(0, 0);
763    }
764    // Turn this off if you get problems
765    // Used to be automatically set
766    clpSolver->setSpecialOptions(clpSolver->specialOptions() | (128 + 64 - 128));
767    ClpSimplex *lpSolver = clpSolver->getModelPtr();
768    lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
769    lpSolver->setSpecialOptions(lpSolver->specialOptions() | (/*16384+*/ 4096 + 512 + 128));
770  }
771#endif
772#ifdef HISTORY_STATISTICS
773  getHistoryStatistics_ = false;
774#endif
775#ifdef COIN_DEVELOP
776  int status = 0;
777#endif
778  int logLevel = model_->logLevel();
779#define LEN_PRINT 250
780  char generalPrint[LEN_PRINT];
781  // Do presolve to see if possible
782  int numberColumns = solver->getNumCols();
783  char *reset = NULL;
784  int returnCode = 1;
785  int saveModelOptions = model_->specialOptions();
786  //assert ((saveModelOptions&2048) == 0);
787  model_->setSpecialOptions(saveModelOptions | 2048);
788  if (fractionSmall < 1.0) {
789    int saveLogLevel = solver->messageHandler()->logLevel();
790    if (saveLogLevel == 1)
791      solver->messageHandler()->setLogLevel(0);
792    OsiPresolve *pinfo = new OsiPresolve();
793    int presolveActions = 0;
794    // Allow dual stuff on integers
795    presolveActions = 1;
796    // Do not allow all +1 to be tampered with
797    //if (allPlusOnes)
798    //presolveActions |= 2;
799    // allow transfer of costs
800    // presolveActions |= 4;
801    pinfo->setPresolveActions(presolveActions);
802    OsiSolverInterface *presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
803    delete pinfo;
804    // see if too big
805
806    if (presolvedModel) {
807      int afterRows = presolvedModel->getNumRows();
808      int afterCols = presolvedModel->getNumCols();
809      //#define COIN_DEVELOP
810#ifdef COIN_DEVELOP_z
811      if (numberNodes < 0) {
812        solver->writeMpsNative("before.mps", NULL, NULL, 2, 1);
813        presolvedModel->writeMpsNative("after1.mps", NULL, NULL, 2, 1);
814      }
815#endif
816      delete presolvedModel;
817      double ratio = sizeRatio(afterRows - shiftRows, afterCols,
818        numberRowsStart, numberColumnsStart);
819      double after = 2 * afterRows + afterCols;
820      if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
821        // Need code to try again to compress further using used
822        const int *used = model_->usedInSolution();
823        int maxUsed = 0;
824        int iColumn;
825        const double *lower = solver->getColLower();
826        const double *upper = solver->getColUpper();
827        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
828          if (upper[iColumn] > lower[iColumn]) {
829            if (solver->isBinary(iColumn))
830              maxUsed = CoinMax(maxUsed, used[iColumn]);
831          }
832        }
833        if (maxUsed) {
834          reset = new char[numberColumns];
835          int nFix = 0;
836          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
837            reset[iColumn] = 0;
838            if (upper[iColumn] > lower[iColumn]) {
839              if (solver->isBinary(iColumn) && used[iColumn] == maxUsed) {
840                bool setValue = true;
841                if (maxUsed == 1) {
842                  double randomNumber = randomNumberGenerator_.randomDouble();
843                  if (randomNumber > 0.3)
844                    setValue = false;
845                }
846                if (setValue) {
847                  reset[iColumn] = 1;
848                  solver->setColLower(iColumn, 1.0);
849                  nFix++;
850                }
851              }
852            }
853          }
854          pinfo = new OsiPresolve();
855          presolveActions = 0;
856          // Allow dual stuff on integers
857          presolveActions = 1;
858          // Do not allow all +1 to be tampered with
859          //if (allPlusOnes)
860          //presolveActions |= 2;
861          // allow transfer of costs
862          // presolveActions |= 4;
863          pinfo->setPresolveActions(presolveActions);
864          presolvedModel = pinfo->presolvedModel(*solver, 1.0e-8, true, 2);
865          delete pinfo;
866          if (presolvedModel) {
867            // see if too big
868            int afterRows2 = presolvedModel->getNumRows();
869            int afterCols2 = presolvedModel->getNumCols();
870            delete presolvedModel;
871            double ratio = sizeRatio(afterRows2 - shiftRows, afterCols2,
872              numberRowsStart, numberColumnsStart);
873            double after = 2 * afterRows2 + afterCols2;
874            if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
875              sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
876                solver->getNumRows(), solver->getNumCols(),
877                afterRows, afterCols, nFix, afterRows2, afterCols2);
878              // If much too big - give up
879              if (ratio > 0.75)
880                returnCode = -1;
881            } else {
882              sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
883                solver->getNumRows(), solver->getNumCols(),
884                afterRows, afterCols, nFix, afterRows2, afterCols2);
885            }
886            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
887              << generalPrint
888              << CoinMessageEol;
889          } else {
890            returnCode = 2; // infeasible
891          }
892        }
893      } else if (ratio > fractionSmall && after > 300 && numberNodes >= 0) {
894        returnCode = -1;
895      }
896    } else {
897      returnCode = 2; // infeasible
898    }
899    solver->messageHandler()->setLogLevel(saveLogLevel);
900  }
901  if (returnCode == 2 || returnCode == -1) {
902    model_->setSpecialOptions(saveModelOptions);
903    delete[] reset;
904#ifdef HISTORY_STATISTICS
905    getHistoryStatistics_ = true;
906#endif
907    //printf("small no good\n");
908    return returnCode;
909  }
910  // Reduce printout
911  bool takeHint;
912  OsiHintStrength strength;
913  solver->getHintParam(OsiDoReducePrint, takeHint, strength);
914  solver->setHintParam(OsiDoReducePrint, true, OsiHintTry);
915  solver->setHintParam(OsiDoPresolveInInitial, false, OsiHintTry);
916  double signedCutoff = cutoff * solver->getObjSense();
917  solver->setDblParam(OsiDualObjectiveLimit, signedCutoff);
918  solver->initialSolve();
919  if (solver->isProvenOptimal()) {
920    CglPreProcess process;
921    OsiSolverInterface *solver2 = NULL;
922    if ((model_->moreSpecialOptions() & 65536) != 0)
923      process.setOptions(2 + 4 + 8 + 16); // no cuts
924    else
925      process.setOptions(16); // no complicated dupcol stuff
926    /* Do not try and produce equality cliques and
927           do up to 2 passes (normally) 5 if restart */
928    int numberPasses = 2;
929    if ((model_->moreSpecialOptions2() & 16) != 0) {
930      // quick
931      process.setOptions(2 + 4 + 8 + 16); // no cuts
932      numberPasses = 1;
933    }
934    if (numberNodes < 0) {
935      numberPasses = 5;
936      // Say some rows cuts
937      int numberRows = solver->getNumRows();
938      if (numberNodes_ < numberRows && true /* think */) {
939        char *type = new char[numberRows];
940        memset(type, 0, numberNodes_);
941        memset(type + numberNodes_, 1, numberRows - numberNodes_);
942        process.passInRowTypes(type, numberRows);
943        delete[] type;
944      }
945    }
946    if (logLevel <= 1)
947      process.messageHandler()->setLogLevel(0);
948    if (!solver->defaultHandler() && solver->messageHandler()->logLevel(0) != -1000)
949      process.passInMessageHandler(solver->messageHandler());
950#ifdef CGL_DEBUG
951    /*
952          We're debugging. (specialOptions 1)
953        */
954    if ((model_->specialOptions() & 1) != 0) {
955      const OsiRowCutDebugger *debugger = solver->getRowCutDebugger();
956      if (debugger) {
957        process.setApplicationData(const_cast< double * >(debugger->optimalSolution()));
958      }
959    }
960#endif
961#ifdef COIN_HAS_CLP
962    OsiClpSolverInterface *clpSolver = dynamic_cast< OsiClpSolverInterface * >(solver);
963    // See if SOS
964    if (clpSolver && clpSolver->numberSOS()) {
965      // SOS
966      int numberSOS = clpSolver->numberSOS();
967      const CoinSet *setInfo = clpSolver->setInfo();
968      int *sosStart = new int[numberSOS + 1];
969      char *sosType = new char[numberSOS];
970      int i;
971      int nTotal = 0;
972      sosStart[0] = 0;
973      for (i = 0; i < numberSOS; i++) {
974        int type = setInfo[i].setType();
975        int n = setInfo[i].numberEntries();
976        sosType[i] = static_cast< char >(type);
977        nTotal += n;
978        sosStart[i + 1] = nTotal;
979      }
980      int *sosIndices = new int[nTotal];
981      double *sosReference = new double[nTotal];
982      for (i = 0; i < numberSOS; i++) {
983        int n = setInfo[i].numberEntries();
984        const int *which = setInfo[i].which();
985        const double *weights = setInfo[i].weights();
986        int base = sosStart[i];
987        for (int j = 0; j < n; j++) {
988          int k = which[j];
989          sosIndices[j + base] = k;
990          sosReference[j + base] = weights ? weights[j] : static_cast< double >(j);
991        }
992      }
993      int numberColumns = solver->getNumCols();
994      char *prohibited = new char[numberColumns];
995      memset(prohibited, 0, numberColumns);
996      int n = sosStart[numberSOS];
997      for (int i = 0; i < n; i++) {
998        int iColumn = sosIndices[i];
999        prohibited[iColumn] = 1;
1000      }
1001      delete[] sosIndices;
1002      delete[] sosReference;
1003      delete[] sosStart;
1004      delete[] sosType;
1005      process.passInProhibited(prohibited, numberColumns);
1006      delete[] prohibited;
1007    }
1008#endif
1009    solver2 = process.preProcessNonDefault(*solver, false,
1010      numberPasses);
1011    if (!solver2) {
1012      if (logLevel > 1)
1013        printf("Pre-processing says infeasible\n");
1014      returnCode = 2; // so will be infeasible
1015    } else {
1016#ifdef COIN_DEVELOP_z
1017      if (numberNodes < 0) {
1018        solver2->writeMpsNative("after2.mps", NULL, NULL, 2, 1);
1019      }
1020#endif
1021      // see if too big
1022      double ratio = sizeRatio(solver2->getNumRows() - shiftRows, solver2->getNumCols(),
1023        numberRowsStart, numberColumnsStart);
1024      double after = 2 * solver2->getNumRows() + solver2->getNumCols();
1025      if (ratio > fractionSmall && (after > 300 || numberNodes < 0)) {
1026        sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
1027          solver->getNumRows(), solver->getNumCols(),
1028          solver2->getNumRows(), solver2->getNumCols());
1029        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1030          << generalPrint
1031          << CoinMessageEol;
1032        returnCode = -1;
1033        //printf("small no good2\n");
1034      } else {
1035        sprintf(generalPrint, "Full problem %d rows %d columns, reduced to %d rows %d columns",
1036          solver->getNumRows(), solver->getNumCols(),
1037          solver2->getNumRows(), solver2->getNumCols());
1038        model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1039          << generalPrint
1040          << CoinMessageEol;
1041      }
1042#ifdef CGL_DEBUG
1043      if ((model_->specialOptions() & 1) != 0) {
1044        const OsiRowCutDebugger *debugger = solver2->getRowCutDebugger();
1045        if (debugger) {
1046          printf("On optimal path after preprocessing\n");
1047        }
1048      }
1049#endif
1050      if (returnCode == 1) {
1051        solver2->resolve();
1052        CbcModel model(*solver2);
1053        double startTime = model_->getDblParam(CbcModel::CbcStartSeconds);
1054        model.setDblParam(CbcModel::CbcStartSeconds, startTime);
1055        // move seed across
1056        model.randomNumberGenerator()->setSeed(model_->randomNumberGenerator()->getSeed());
1057#ifdef COIN_HAS_CLP
1058        // redo SOS
1059        OsiClpSolverInterface *clpSolver
1060          = dynamic_cast< OsiClpSolverInterface * >(model.solver());
1061        if (clpSolver && clpSolver->numberSOS()) {
1062          int numberColumns = clpSolver->getNumCols();
1063          const int *originalColumns = process.originalColumns();
1064          CoinSet *setInfo = const_cast< CoinSet * >(clpSolver->setInfo());
1065          int numberSOS = clpSolver->numberSOS();
1066          for (int iSOS = 0; iSOS < numberSOS; iSOS++) {
1067            //int type = setInfo[iSOS].setType();
1068            int n = setInfo[iSOS].numberEntries();
1069            int *which = setInfo[iSOS].modifiableWhich();
1070            double *weights = setInfo[iSOS].modifiableWeights();
1071            int n2 = 0;
1072            for (int j = 0; j < n; j++) {
1073              int iColumn = which[j];
1074              int i;
1075              for (i = 0; i < numberColumns; i++) {
1076                if (originalColumns[i] == iColumn)
1077                  break;
1078              }
1079              if (i < numberColumns) {
1080                which[n2] = i;
1081                weights[n2++] = weights[j];
1082              }
1083            }
1084            setInfo[iSOS].setNumberEntries(n2);
1085          }
1086        }
1087#endif
1088        if (numberNodes >= 0) {
1089          // normal
1090          model.setSpecialOptions(saveModelOptions | 2048);
1091          if (logLevel <= 1 && feasibilityPumpOptions_ != -3)
1092            model.setLogLevel(0);
1093          else
1094            model.setLogLevel(logLevel);
1095          // No small fathoming
1096          model.setFastNodeDepth(-1);
1097          model.setCutoff(signedCutoff);
1098          model.setStrongStrategy(0);
1099          // Don't do if original fraction > 1.0 and too large
1100          if (fractionSmall_ > 1.0 && fractionSmall_ < 1000000.0) {
1101            /* 1.4 means -1 nodes if >.4
1102                         2.4 means -1 nodes if >.5 and 0 otherwise
1103                         3.4 means -1 nodes if >.6 and 0 or 5
1104                         4.4 means -1 nodes if >.7 and 0, 5 or 10
1105                      */
1106            double fraction = fractionSmall_ - floor(fractionSmall_);
1107            if (ratio > fraction) {
1108              int type = static_cast< int >(floor(fractionSmall_ * 0.1));
1109              int over = static_cast< int >(ceil(ratio - fraction));
1110              int maxNodes[] = { -1, 0, 5, 10 };
1111              if (type > over)
1112                numberNodes = maxNodes[type - over];
1113              else
1114                numberNodes = -1;
1115            }
1116          }
1117          model.setMaximumNodes(numberNodes);
1118          model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
1119          if ((saveModelOptions & 2048) == 0)
1120            model.setMoreSpecialOptions(model_->moreSpecialOptions());
1121          model.setMoreSpecialOptions2(model_->moreSpecialOptions2());
1122          // off conflict analysis
1123          model.setMoreSpecialOptions(model.moreSpecialOptions() & ~4194304);
1124
1125          // Lightweight
1126          CbcStrategyDefaultSubTree strategy(model_, 1, 5, 1, 0);
1127          model.setStrategy(strategy);
1128          model.solver()->setIntParam(OsiMaxNumIterationHotStart, 10);
1129          model.setMaximumCutPassesAtRoot(CoinMin(20, CoinAbs(model_->getMaximumCutPassesAtRoot())));
1130          model.setMaximumCutPasses(CoinMin(10, model_->getMaximumCutPasses()));
1131          // Set best solution (even if bad for this submodel)
1132          if (model_->bestSolution()) {
1133            const double *bestSolution = model_->bestSolution();
1134            int numberColumns2 = model.solver()->getNumCols();
1135            double *bestSolution2 = new double[numberColumns2];
1136            const int *originalColumns = process.originalColumns();
1137            for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1138              int jColumn = originalColumns[iColumn];
1139              bestSolution2[iColumn] = bestSolution[jColumn];
1140            }
1141            model.setBestSolution(bestSolution2, numberColumns2,
1142              1.0e50,
1143              false);
1144            model.setSolutionCount(1);
1145            maximumSolutions++;
1146            delete[] bestSolution2;
1147          }
1148        } else {
1149          // modify for event handler
1150          model.setSpecialOptions(saveModelOptions);
1151          model_->messageHandler()->message(CBC_RESTART, model_->messages())
1152            << solver2->getNumRows() << solver2->getNumCols()
1153            << CoinMessageEol;
1154          // going for full search and copy across more stuff
1155          model.gutsOfCopy(*model_, 2);
1156         
1157          if (model.solver()->getCGraph()) {
1158            model.solver()->buildCGraph();
1159          }
1160
1161#ifdef CGL_DEBUG
1162          if ((model_->specialOptions() & 1) != 0) {
1163            const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger();
1164            if (debugger) {
1165              printf("On optimal path BB\n");
1166            }
1167          }
1168#endif
1169          assert(!model_->heuristicModel());
1170          model_->setHeuristicModel(&model);
1171          for (int i = 0; i < model.numberCutGenerators(); i++) {
1172            CbcCutGenerator *generator = model.cutGenerator(i);
1173            CglGomory *gomory = dynamic_cast< CglGomory * >(generator->generator());
1174            if (gomory && gomory->originalSolver())
1175              gomory->passInOriginalSolver(model.solver());
1176            generator->setTiming(true);
1177            // Turn on if was turned on
1178            int iOften = model_->cutGenerator(i)->howOften();
1179#ifdef CLP_INVESTIGATE
1180            printf("Gen %d often %d %d\n",
1181              i, generator->howOften(),
1182              iOften);
1183#endif
1184            if (iOften > 0)
1185              generator->setHowOften(iOften % 1000000);
1186            if (model_->cutGenerator(i)->howOftenInSub() == -200)
1187              generator->setHowOften(-100);
1188          }
1189          model.setCutoff(signedCutoff);
1190          // make sure can't do nested search! but allow heuristics
1191          model.setSpecialOptions((model.specialOptions() & (~(512 + 2048))) | 1024);
1192          // but say we are doing full search
1193          model.setSpecialOptions(model.specialOptions() | 67108864);
1194          bool takeHint;
1195          OsiHintStrength strength;
1196          // Switch off printing if asked to
1197          model_->solver()->getHintParam(OsiDoReducePrint, takeHint, strength);
1198          model.solver()->setHintParam(OsiDoReducePrint, takeHint, strength);
1199          // no cut generators if none in parent
1200          CbcStrategyDefault
1201            strategy(model_->numberCutGenerators() ? 1 : -1,
1202              model_->numberStrong(),
1203              model_->numberBeforeTrust());
1204          // Set up pre-processing - no
1205          strategy.setupPreProcessing(0); // was (4);
1206          model.setStrategy(strategy);
1207          //model.solver()->writeMps("crunched");
1208          int numberCuts = process.cuts().sizeRowCuts();
1209          if (numberCuts) {
1210            // add in cuts
1211            CglStored cuts = process.cuts();
1212            model.addCutGenerator(&cuts, 1, "Stored from first");
1213            model.cutGenerator(model.numberCutGenerators() - 1)->setGlobalCuts(true);
1214          }
1215        }
1216        // Do search
1217        if (logLevel > 1)
1218          model_->messageHandler()->message(CBC_START_SUB, model_->messages())
1219            << name
1220            << model.getMaximumNodes()
1221            << CoinMessageEol;
1222        // probably faster to use a basis to get integer solutions
1223        model.setSpecialOptions(model.specialOptions() | 2);
1224#ifdef CBC_THREAD
1225        if (model_->getNumberThreads() > 0 && (model_->getThreadMode() & 4) != 0) {
1226          // See if at root node
1227          bool atRoot = model_->getNodeCount() == 0;
1228          int passNumber = model_->getCurrentPassNumber();
1229          if (atRoot && passNumber == 1)
1230            model.setNumberThreads(model_->getNumberThreads());
1231        }
1232#endif
1233        model.setParentModel(*model_);
1234        model.setMaximumSolutions(maximumSolutions);
1235        model.setOriginalColumns(process.originalColumns());
1236        model.setSearchStrategy(-1);
1237        // If no feasibility pump then insert a lightweight one
1238        if (feasibilityPumpOptions_ >= 0 || feasibilityPumpOptions_ == -2) {
1239          CbcHeuristicFPump *fpump = NULL;
1240          for (int i = 0; i < model.numberHeuristics(); i++) {
1241            CbcHeuristicFPump *pump = dynamic_cast< CbcHeuristicFPump * >(model.heuristic(i));
1242            if (pump) {
1243              fpump = pump;
1244              break;
1245            }
1246          }
1247          if (!fpump) {
1248            CbcHeuristicFPump heuristic4;
1249            // use any cutoff
1250            heuristic4.setFakeCutoff(0.5 * COIN_DBL_MAX);
1251            if (fractionSmall_ <= 1.0)
1252              heuristic4.setMaximumPasses(10);
1253            int pumpTune = feasibilityPumpOptions_;
1254            if (pumpTune == -2)
1255              pumpTune = 4; // proximity
1256            if (pumpTune > 0) {
1257              /*
1258                            >=10000000 for using obj
1259                            >=1000000 use as accumulate switch
1260                            >=1000 use index+1 as number of large loops
1261                            >=100 use 0.05 objvalue as increment
1262                            %100 == 10,20 etc for experimentation
1263                            1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
1264                            4 and static continuous, 5 as 3 but no internal integers
1265                            6 as 3 but all slack basis!
1266                            */
1267              double value = solver2->getObjSense() * solver2->getObjValue();
1268              int w = pumpTune / 10;
1269              int ix = w % 10;
1270              w /= 10;
1271              int c = w % 10;
1272              w /= 10;
1273              int r = w;
1274              int accumulate = r / 1000;
1275              r -= 1000 * accumulate;
1276              if (accumulate >= 10) {
1277                int which = accumulate / 10;
1278                accumulate -= 10 * which;
1279                which--;
1280                // weights and factors
1281                double weight[] = { 0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0 };
1282                double factor[] = { 0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5 };
1283                heuristic4.setInitialWeight(weight[which]);
1284                heuristic4.setWeightFactor(factor[which]);
1285              }
1286              // fake cutoff
1287              if (c) {
1288                double cutoff;
1289                solver2->getDblParam(OsiDualObjectiveLimit, cutoff);
1290                cutoff = CoinMin(cutoff, value + 0.1 * fabs(value) * c);
1291                heuristic4.setFakeCutoff(cutoff);
1292              }
1293              if (r) {
1294                // also set increment
1295                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
1296                double increment = 0.0;
1297                heuristic4.setAbsoluteIncrement(increment);
1298                heuristic4.setAccumulate(accumulate);
1299                heuristic4.setMaximumRetries(r + 1);
1300              }
1301              pumpTune = pumpTune % 100;
1302              if (pumpTune == 6)
1303                pumpTune = 13;
1304              if (pumpTune != 13)
1305                pumpTune = pumpTune % 10;
1306              heuristic4.setWhen(pumpTune);
1307              if (ix) {
1308                heuristic4.setFeasibilityPumpOptions(ix * 10);
1309              }
1310            }
1311            model.addHeuristic(&heuristic4, "feasibility pump", 0);
1312          }
1313        } else if (feasibilityPumpOptions_ == -3) {
1314          // add all (except this)
1315          for (int i = 0; i < model_->numberHeuristics(); i++) {
1316            if (strcmp(heuristicName(), model_->heuristic(i)->heuristicName()))
1317              model.addHeuristic(model_->heuristic(i));
1318          }
1319        }
1320        // modify heuristics
1321        for (int i = 0; i < model.numberHeuristics(); i++) {
1322          // reset lastNode
1323          CbcHeuristicRINS *rins = dynamic_cast< CbcHeuristicRINS * >(model.heuristic(i));
1324          if (rins) {
1325            rins->setLastNode(-1000);
1326            rins->setSolutionCount(0);
1327          }
1328        }
1329        //printf("sol %x\n",inputSolution_);
1330#ifdef CGL_DEBUG
1331        if ((model_->specialOptions() & 1) != 0) {
1332          const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger();
1333          if (debugger) {
1334            printf("On optimal path CC\n");
1335          }
1336        }
1337#endif
1338        if (inputSolution_) {
1339          // translate and add a serendipity heuristic
1340          int numberColumns = solver2->getNumCols();
1341          const int *which = process.originalColumns();
1342          OsiSolverInterface *solver3 = solver2->clone();
1343          for (int i = 0; i < numberColumns; i++) {
1344            if (isHeuristicInteger(solver3, i)) {
1345              int k = which[i];
1346              double value = inputSolution_[k];
1347              //if (value)
1348              //printf("orig col %d now %d val %g\n",
1349              //       k,i,value);
1350              solver3->setColLower(i, value);
1351              solver3->setColUpper(i, value);
1352            }
1353          }
1354          solver3->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1355          solver3->resolve();
1356          if (!solver3->isProvenOptimal()) {
1357            // Try just setting nonzeros
1358            OsiSolverInterface *solver4 = solver2->clone();
1359            for (int i = 0; i < numberColumns; i++) {
1360              if (isHeuristicInteger(solver4, i)) {
1361                int k = which[i];
1362                double value = floor(inputSolution_[k] + 0.5);
1363                if (value) {
1364                  solver3->setColLower(i, value);
1365                  solver3->setColUpper(i, value);
1366                }
1367              }
1368            }
1369            solver4->setDblParam(OsiDualObjectiveLimit, COIN_DBL_MAX);
1370            solver4->resolve();
1371            int nBad = -1;
1372            if (solver4->isProvenOptimal()) {
1373              nBad = 0;
1374              const double *solution = solver4->getColSolution();
1375              for (int i = 0; i < numberColumns; i++) {
1376                if (isHeuristicInteger(solver4, i)) {
1377                  double value = floor(solution[i] + 0.5);
1378                  if (fabs(value - solution[i]) > 1.0e-6)
1379                    nBad++;
1380                }
1381              }
1382            }
1383            if (nBad) {
1384              delete solver4;
1385            } else {
1386              delete solver3;
1387              solver3 = solver4;
1388            }
1389          }
1390          if (solver3->isProvenOptimal()) {
1391            // good
1392            CbcSerendipity heuristic(model);
1393            double value = solver3->getObjSense() * solver3->getObjValue();
1394            heuristic.setInputSolution(solver3->getColSolution(), value);
1395            value = value + 1.0e-7 * (1.0 + fabs(value));
1396            value *= solver3->getObjSense();
1397            model.setCutoff(value);
1398            model.addHeuristic(&heuristic, "Previous solution", 0);
1399            //printf("added seren\n");
1400          } else {
1401            double value = model_->getMinimizationObjValue();
1402            value = value + 1.0e-7 * (1.0 + fabs(value));
1403            value *= solver3->getObjSense();
1404            model.setCutoff(value);
1405            sprintf(generalPrint, "Unable to insert previous solution - using cutoff of %g",
1406              value);
1407            model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1408              << generalPrint
1409              << CoinMessageEol;
1410#ifdef CLP_INVESTIGATE
1411            printf("NOT added seren\n");
1412            solver3->writeMps("bad_seren");
1413            solver->writeMps("orig_seren");
1414#endif
1415          }
1416          delete solver3;
1417        }
1418        if (model_->searchStrategy() == 2) {
1419          model.setNumberStrong(5);
1420          model.setNumberBeforeTrust(5);
1421        }
1422        if (model.getNumCols()) {
1423          if (numberNodes >= 0) {
1424            setCutAndHeuristicOptions(model);
1425            // not too many iterations
1426            model.setMaximumNumberIterations(iterationMultiplier * (numberNodes + 10));
1427            // Not fast stuff
1428            model.setFastNodeDepth(-1);
1429            //model.solver()->writeMps("before");
1430          } else if (model.fastNodeDepth() >= 1000000) {
1431            // already set
1432            model.setFastNodeDepth(model.fastNodeDepth() - 1000000);
1433          }
1434          model.setWhenCuts(999998);
1435#define ALWAYS_DUAL
1436#ifdef ALWAYS_DUAL
1437          OsiSolverInterface *solverD = model.solver();
1438          bool takeHint;
1439          OsiHintStrength strength;
1440          solverD->getHintParam(OsiDoDualInResolve, takeHint, strength);
1441          solverD->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
1442#endif
1443          model.passInEventHandler(model_->getEventHandler());
1444          // say model_ is sitting there
1445          int saveOptions = model_->specialOptions();
1446          model_->setSpecialOptions(saveOptions | 1048576);
1447          // and switch off debugger
1448          model.setSpecialOptions(model.specialOptions() & (~1));
1449#if 0 //def COIN_HAS_CLP
1450                    OsiClpSolverInterface * clpSolver
1451                      = dynamic_cast<OsiClpSolverInterface *> (model.solver());
1452                    if (clpSolver)
1453                      clpSolver->zapDebugger();
1454#endif
1455#ifdef CONFLICT_CUTS
1456          if ((model_->moreSpecialOptions() & 4194304) != 0)
1457            model.zapGlobalCuts();
1458#endif
1459#ifdef CGL_DEBUG
1460          if ((model_->specialOptions() & 1) != 0) {
1461            const OsiRowCutDebugger *debugger = model.solver()->getRowCutDebugger();
1462            if (debugger) {
1463              printf("On optimal path DD\n");
1464            }
1465          }
1466#endif
1467          model.setPreProcess(&process);
1468          model.branchAndBound();
1469          model_->setHeuristicModel(NULL);
1470          model_->setSpecialOptions(saveOptions);
1471#ifdef ALWAYS_DUAL
1472          solverD = model.solver();
1473          solverD->setHintParam(OsiDoDualInResolve, takeHint, strength);
1474#endif
1475          numberNodesDone_ = model.getNodeCount();
1476#ifdef COIN_DEVELOP
1477          printf("sub branch %d nodes, %d iterations - max %d\n",
1478            model.getNodeCount(), model.getIterationCount(),
1479            100 * (numberNodes + 10));
1480#endif
1481          if (numberNodes < 0) {
1482            model_->incrementIterationCount(model.getIterationCount());
1483            model_->incrementNodeCount(model.getNodeCount());
1484            // update best solution (in case ctrl-c)
1485            // !!! not a good idea - think a bit harder
1486            //model_->setMinimizationObjValue(model.getMinimizationObjValue());
1487            for (int iGenerator = 0; iGenerator < model.numberCutGenerators(); iGenerator++) {
1488              CbcCutGenerator *generator = model.cutGenerator(iGenerator);
1489              sprintf(generalPrint,
1490                "%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts (%.3f seconds)",
1491                generator->cutGeneratorName(),
1492                generator->numberTimesEntered(),
1493                generator->numberCutsInTotal() + generator->numberColumnCuts(),
1494                generator->numberCutsActive(),
1495                generator->timeInCutGenerator());
1496              CglStored *stored = dynamic_cast< CglStored * >(generator->generator());
1497              if (stored && !generator->numberCutsInTotal())
1498                continue;
1499#ifndef CLP_INVESTIGATE
1500              CglImplication *implication = dynamic_cast< CglImplication * >(generator->generator());
1501              if (implication)
1502                continue;
1503#endif
1504              model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
1505                << generalPrint
1506                << CoinMessageEol;
1507            }
1508          }
1509        } else {
1510          // empty model
1511          model.setMinimizationObjValue(model.solver()->getObjSense() * model.solver()->getObjValue());
1512        }
1513        if (logLevel > 1)
1514          model_->messageHandler()->message(CBC_END_SUB, model_->messages())
1515            << name
1516            << CoinMessageEol;
1517        if (model.getMinimizationObjValue() < CoinMin(cutoff, 1.0e30)) {
1518          // solution
1519          if (model.getNumCols())
1520            returnCode = model.isProvenOptimal() ? 3 : 1;
1521          else
1522            returnCode = 3;
1523            // post process
1524#ifdef COIN_HAS_CLP
1525          OsiClpSolverInterface *clpSolver = dynamic_cast< OsiClpSolverInterface * >(model.solver());
1526          if (clpSolver) {
1527            ClpSimplex *lpSolver = clpSolver->getModelPtr();
1528            lpSolver->setSpecialOptions(lpSolver->specialOptions() | 0x01000000); // say is Cbc (and in branch and bound)
1529          }
1530#endif
1531          //if (fractionSmall_ < 1000000.0)
1532          process.postProcess(*model.solver());
1533          if (solver->isProvenOptimal() && solver->getObjValue() * solver->getObjSense() < cutoff) {
1534            // Solution now back in solver
1535            int numberColumns = solver->getNumCols();
1536            memcpy(newSolution, solver->getColSolution(),
1537              numberColumns * sizeof(double));
1538            newSolutionValue = model.getMinimizationObjValue();
1539#ifdef COIN_HAS_CLP
1540            if (clpSolver) {
1541              if (clpSolver && clpSolver->numberSOS()) {
1542                // SOS
1543                int numberSOS = clpSolver->numberSOS();
1544                const CoinSet *setInfo = clpSolver->setInfo();
1545                int i;
1546                for (i = 0; i < numberSOS; i++) {
1547#ifndef NDEBUG
1548                  int type = setInfo[i].setType();
1549#endif
1550                  int n = setInfo[i].numberEntries();
1551                  const int *which = setInfo[i].which();
1552                  int first = -1;
1553                  int last = -1;
1554                  for (int j = 0; j < n; j++) {
1555                    int iColumn = which[j];
1556                    if (fabs(newSolution[iColumn]) > 1.0e-7) {
1557                      last = j;
1558                      if (first < 0)
1559                        first = j;
1560                    }
1561                  }
1562                  assert(last - first < type);
1563                  for (int j = 0; j < n; j++) {
1564                    if (j < first || j > last) {
1565                      int iColumn = which[j];
1566                      // do I want to fix??
1567                      solver->setColLower(iColumn, 0.0);
1568                      solver->setColUpper(iColumn, 0.0);
1569                      newSolution[iColumn] = 0.0;
1570                    }
1571                  }
1572                }
1573              }
1574            }
1575#endif
1576          } else {
1577            // odd - but no good
1578            returnCode = 0; // so will be infeasible
1579          }
1580        } else {
1581          // no good
1582          returnCode = model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
1583        }
1584        int totalNumberIterations = model.getIterationCount() + process.numberIterationsPre() + process.numberIterationsPost();
1585        if (totalNumberIterations > 100 * (numberNodes + 10)
1586          && fractionSmall_ < 1000000.0) {
1587          // only allow smaller problems
1588          fractionSmall = fractionSmall_;
1589          fractionSmall_ *= 0.9;
1590#ifdef CLP_INVESTIGATE
1591          printf("changing fractionSmall from %g to %g for %s as %d iterations\n",
1592            fractionSmall, fractionSmall_, name.c_str(), totalNumberIterations);
1593#endif
1594        }
1595        if (model.status() == 5)
1596          model_->sayEventHappened();
1597#ifdef COIN_DEVELOP
1598        if (model.isProvenInfeasible())
1599          status = 1;
1600        else if (model.isProvenOptimal())
1601          status = 2;
1602#endif
1603      }
1604    }
1605  } else {
1606    returnCode = 2; // infeasible finished
1607    if (logLevel > 1) {
1608      printf("Infeasible on initial solve\n");
1609    }
1610  }
1611  model_->setSpecialOptions(saveModelOptions);
1612  model_->setLogLevel(logLevel);
1613  if (returnCode == 1 || returnCode == 2) {
1614    OsiSolverInterface *solverC = model_->continuousSolver();
1615    if (false && solverC) {
1616      const double *lower = solver->getColLower();
1617      const double *upper = solver->getColUpper();
1618      const double *lowerC = solverC->getColLower();
1619      const double *upperC = solverC->getColUpper();
1620      bool good = true;
1621      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1622        if (isHeuristicInteger(solverC, iColumn)) {
1623          if (lower[iColumn] > lowerC[iColumn] && upper[iColumn] < upperC[iColumn]) {
1624            good = false;
1625            printf("CUT - can't add\n");
1626            break;
1627          }
1628        }
1629      }
1630      if (good) {
1631        double *cut = new double[numberColumns];
1632        int *which = new int[numberColumns];
1633        double rhs = -1.0;
1634        int n = 0;
1635        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1636          if (isHeuristicInteger(solverC, iColumn)) {
1637            if (lower[iColumn] == upperC[iColumn]) {
1638              rhs += lower[iColumn];
1639              cut[n] = 1.0;
1640              which[n++] = iColumn;
1641            } else if (upper[iColumn] == lowerC[iColumn]) {
1642              rhs -= upper[iColumn];
1643              cut[n] = -1.0;
1644              which[n++] = iColumn;
1645            }
1646          }
1647        }
1648        printf("CUT has %d entries\n", n);
1649        OsiRowCut newCut;
1650        newCut.setLb(-COIN_DBL_MAX);
1651        newCut.setUb(rhs);
1652        newCut.setRow(n, which, cut, false);
1653        model_->makeGlobalCut(newCut);
1654        delete[] cut;
1655        delete[] which;
1656      }
1657    }
1658#ifdef COIN_DEVELOP
1659    if (status == 1)
1660      printf("heuristic could add cut because infeasible (%s)\n", heuristicName_.c_str());
1661    else if (status == 2)
1662      printf("heuristic could add cut because optimal (%s)\n", heuristicName_.c_str());
1663#endif
1664  }
1665  if (reset) {
1666    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1667      if (reset[iColumn])
1668        solver->setColLower(iColumn, 0.0);
1669    }
1670    delete[] reset;
1671  }
1672#ifdef HISTORY_STATISTICS
1673  getHistoryStatistics_ = true;
1674#endif
1675  solver->setHintParam(OsiDoReducePrint, takeHint, strength);
1676  return returnCode;
1677}
1678// Set input solution
1679void CbcHeuristic::setInputSolution(const double *solution, double objValue)
1680{
1681  delete[] inputSolution_;
1682  inputSolution_ = NULL;
1683  if (model_ && solution) {
1684    int numberColumns = model_->getNumCols();
1685    inputSolution_ = new double[numberColumns + 1];
1686    memcpy(inputSolution_, solution, numberColumns * sizeof(double));
1687    inputSolution_[numberColumns] = objValue;
1688  }
1689}
1690
1691//##############################################################################
1692
1693inline int compare3BranchingObjects(const CbcBranchingObject *br0,
1694  const CbcBranchingObject *br1)
1695{
1696  const int t0 = br0->type();
1697  const int t1 = br1->type();
1698  if (t0 < t1) {
1699    return -1;
1700  }
1701  if (t0 > t1) {
1702    return 1;
1703  }
1704  return br0->compareOriginalObject(br1);
1705}
1706
1707//==============================================================================
1708
1709inline bool compareBranchingObjects(const CbcBranchingObject *br0,
1710  const CbcBranchingObject *br1)
1711{
1712  return compare3BranchingObjects(br0, br1) < 0;
1713}
1714
1715//==============================================================================
1716
1717void CbcHeuristicNode::gutsOfConstructor(CbcModel &model)
1718{
1719  //  CbcHeurDebugNodes(&model);
1720  CbcNode *node = model.currentNode();
1721  brObj_ = new CbcBranchingObject *[node->depth()];
1722  CbcNodeInfo *nodeInfo = node->nodeInfo();
1723  int cnt = 0;
1724  while (nodeInfo->parentBranch() != NULL) {
1725    const OsiBranchingObject *br = nodeInfo->parentBranch();
1726    const CbcBranchingObject *cbcbr = dynamic_cast< const CbcBranchingObject * >(br);
1727    if (!cbcbr) {
1728      throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
1729        "gutsOfConstructor",
1730        "CbcHeuristicNode",
1731        __FILE__, __LINE__);
1732    }
1733    brObj_[cnt] = cbcbr->clone();
1734    brObj_[cnt]->previousBranch();
1735    ++cnt;
1736    nodeInfo = nodeInfo->parent();
1737  }
1738  std::sort(brObj_, brObj_ + cnt, compareBranchingObjects);
1739  if (cnt <= 1) {
1740    numObjects_ = cnt;
1741  } else {
1742    numObjects_ = 0;
1743    CbcBranchingObject *br = NULL; // What should this be?
1744    for (int i = 1; i < cnt; ++i) {
1745      if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
1746        int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], br != 0);
1747        switch (comp) {
1748        case CbcRangeSame: // the same range
1749        case CbcRangeDisjoint: // disjoint decisions
1750          // should not happen! we are on a chain!
1751          abort();
1752        case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
1753          delete brObj_[i];
1754          break;
1755        case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
1756          delete brObj_[numObjects_];
1757          brObj_[numObjects_] = brObj_[i];
1758          break;
1759        case CbcRangeOverlap: // overlap
1760          delete brObj_[i];
1761          delete brObj_[numObjects_];
1762          brObj_[numObjects_] = br;
1763          break;
1764        }
1765        continue;
1766      } else {
1767        brObj_[++numObjects_] = brObj_[i];
1768      }
1769    }
1770    ++numObjects_;
1771  }
1772}
1773
1774//==============================================================================
1775
1776CbcHeuristicNode::CbcHeuristicNode(CbcModel &model)
1777{
1778  gutsOfConstructor(model);
1779}
1780
1781//==============================================================================
1782
1783double
1784CbcHeuristicNode::distance(const CbcHeuristicNode *node) const
1785{
1786
1787  const double disjointWeight = 1;
1788  const double overlapWeight = 0.4;
1789  const double subsetWeight = 0.2;
1790  int countDisjointWeight = 0;
1791  int countOverlapWeight = 0;
1792  int countSubsetWeight = 0;
1793  int i = 0;
1794  int j = 0;
1795  double dist = 0.0;
1796#ifdef PRINT_DEBUG
1797  printf(" numObjects_ = %i, node->numObjects_ = %i\n",
1798    numObjects_, node->numObjects_);
1799#endif
1800  while (i < numObjects_ && j < node->numObjects_) {
1801    CbcBranchingObject *br0 = brObj_[i];
1802    const CbcBranchingObject *br1 = node->brObj_[j];
1803#ifdef PRINT_DEBUG
1804    const CbcIntegerBranchingObject *brPrint0 = dynamic_cast< const CbcIntegerBranchingObject * >(br0);
1805    const double *downBounds = brPrint0->downBounds();
1806    const double *upBounds = brPrint0->upBounds();
1807    int variable = brPrint0->variable();
1808    int way = brPrint0->way();
1809    printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1810      variable, static_cast< int >(downBounds[0]), static_cast< int >(downBounds[1]),
1811      static_cast< int >(upBounds[0]), static_cast< int >(upBounds[1]), way);
1812    const CbcIntegerBranchingObject *brPrint1 = dynamic_cast< const CbcIntegerBranchingObject * >(br1);
1813    downBounds = brPrint1->downBounds();
1814    upBounds = brPrint1->upBounds();
1815    variable = brPrint1->variable();
1816    way = brPrint1->way();
1817    printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1818      variable, static_cast< int >(downBounds[0]), static_cast< int >(downBounds[1]),
1819      static_cast< int >(upBounds[0]), static_cast< int >(upBounds[1]), way);
1820#endif
1821    const int brComp = compare3BranchingObjects(br0, br1);
1822    if (brComp < 0) {
1823      dist += subsetWeight;
1824      countSubsetWeight++;
1825      ++i;
1826    } else if (brComp > 0) {
1827      dist += subsetWeight;
1828      countSubsetWeight++;
1829      ++j;
1830    } else {
1831      const int comp = br0->compareBranchingObject(br1, false);
1832      switch (comp) {
1833      case CbcRangeSame:
1834        // do nothing
1835        break;
1836      case CbcRangeDisjoint: // disjoint decisions
1837        dist += disjointWeight;
1838        countDisjointWeight++;
1839        break;
1840      case CbcRangeSubset: // subset one way or another
1841      case CbcRangeSuperset:
1842        dist += subsetWeight;
1843        countSubsetWeight++;
1844        break;
1845      case CbcRangeOverlap: // overlap
1846        dist += overlapWeight;
1847        countOverlapWeight++;
1848        break;
1849      }
1850      ++i;
1851      ++j;
1852    }
1853  }
1854  dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
1855  countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
1856  COIN_DETAIL_PRINT(printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
1857    countOverlapWeight, countDisjointWeight));
1858  return dist;
1859}
1860
1861//==============================================================================
1862
1863CbcHeuristicNode::~CbcHeuristicNode()
1864{
1865  for (int i = 0; i < numObjects_; ++i) {
1866    delete brObj_[i];
1867  }
1868  delete[] brObj_;
1869}
1870
1871//==============================================================================
1872
1873double
1874CbcHeuristicNode::minDistance(const CbcHeuristicNodeList &nodeList) const
1875{
1876  double minDist = COIN_DBL_MAX;
1877  for (int i = nodeList.size() - 1; i >= 0; --i) {
1878    minDist = CoinMin(minDist, distance(nodeList.node(i)));
1879  }
1880  return minDist;
1881}
1882
1883//==============================================================================
1884
1885bool CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList &nodeList,
1886  const double threshold) const
1887{
1888  for (int i = nodeList.size() - 1; i >= 0; --i) {
1889    if (distance(nodeList.node(i)) >= threshold) {
1890      continue;
1891    } else {
1892      return true;
1893    }
1894  }
1895  return false;
1896}
1897
1898//==============================================================================
1899
1900double
1901CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList &nodeList) const
1902{
1903  if (nodeList.size() == 0) {
1904    return COIN_DBL_MAX;
1905  }
1906  double sumDist = 0;
1907  for (int i = nodeList.size() - 1; i >= 0; --i) {
1908    sumDist += distance(nodeList.node(i));
1909  }
1910  return sumDist / nodeList.size();
1911}
1912
1913//##############################################################################
1914
1915// Default Constructor
1916CbcRounding::CbcRounding()
1917  : CbcHeuristic()
1918{
1919  // matrix and row copy will automatically be empty
1920  seed_ = 7654321;
1921  down_ = NULL;
1922  up_ = NULL;
1923  equal_ = NULL;
1924  //whereFrom_ |= 16*(1+256); // allow more often
1925}
1926
1927// Constructor from model
1928CbcRounding::CbcRounding(CbcModel &model)
1929  : CbcHeuristic(model)
1930{
1931  // Get a copy of original matrix (and by row for rounding);
1932  assert(model.solver());
1933  if (model.solver()->getNumRows()) {
1934    matrix_ = *model.solver()->getMatrixByCol();
1935    matrixByRow_ = *model.solver()->getMatrixByRow();
1936    validate();
1937  }
1938  down_ = NULL;
1939  up_ = NULL;
1940  equal_ = NULL;
1941  seed_ = 7654321;
1942  //whereFrom_ |= 16*(1+256); // allow more often
1943}
1944
1945// Destructor
1946CbcRounding::~CbcRounding()
1947{
1948  delete[] down_;
1949  delete[] up_;
1950  delete[] equal_;
1951}
1952
1953// Clone
1954CbcHeuristic *
1955CbcRounding::clone() const
1956{
1957  return new CbcRounding(*this);
1958}
1959// Create C++ lines to get to current state
1960void CbcRounding::generateCpp(FILE *fp)
1961{
1962  CbcRounding other;
1963  fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
1964  fprintf(fp, "3  CbcRounding rounding(*cbcModel);\n");
1965  CbcHeuristic::generateCpp(fp, "rounding");
1966  if (seed_ != other.seed_)
1967    fprintf(fp, "3  rounding.setSeed(%d);\n", seed_);
1968  else
1969    fprintf(fp, "4  rounding.setSeed(%d);\n", seed_);
1970  fprintf(fp, "3  cbcModel->addHeuristic(&rounding);\n");
1971}
1972//#define NEW_ROUNDING
1973// Copy constructor
1974CbcRounding::CbcRounding(const CbcRounding &rhs)
1975  : CbcHeuristic(rhs)
1976  , matrix_(rhs.matrix_)
1977  , matrixByRow_(rhs.matrixByRow_)
1978  , seed_(rhs.seed_)
1979{
1980#ifdef NEW_ROUNDING
1981  int numberColumns = matrix_.getNumCols();
1982  down_ = CoinCopyOfArray(rhs.down_, numberColumns);
1983  up_ = CoinCopyOfArray(rhs.up_, numberColumns);
1984  equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
1985#else
1986  down_ = NULL;
1987  up_ = NULL;
1988  equal_ = NULL;
1989#endif
1990}
1991
1992// Assignment operator
1993CbcRounding &
1994CbcRounding::operator=(const CbcRounding &rhs)
1995{
1996  if (this != &rhs) {
1997    CbcHeuristic::operator=(rhs);
1998    matrix_ = rhs.matrix_;
1999    matrixByRow_ = rhs.matrixByRow_;
2000#ifdef NEW_ROUNDING
2001    delete[] down_;
2002    delete[] up_;
2003    delete[] equal_;
2004    int numberColumns = matrix_.getNumCols();
2005    down_ = CoinCopyOfArray(rhs.down_, numberColumns);
2006    up_ = CoinCopyOfArray(rhs.up_, numberColumns);
2007    equal_ = CoinCopyOfArray(rhs.equal_, numberColumns);
2008#else
2009    down_ = NULL;
2010    up_ = NULL;
2011    equal_ = NULL;
2012#endif
2013    seed_ = rhs.seed_;
2014  }
2015  return *this;
2016}
2017
2018// Resets stuff if model changes
2019void CbcRounding::resetModel(CbcModel *model)
2020{
2021  model_ = model;
2022  // Get a copy of original matrix (and by row for rounding);
2023  assert(model_->solver());
2024  matrix_ = *model_->solver()->getMatrixByCol();
2025  matrixByRow_ = *model_->solver()->getMatrixByRow();
2026  validate();
2027}
2028/* Check whether the heuristic should run at all
2029   0 - before cuts at root node (or from doHeuristics)
2030   1 - during cuts at root
2031   2 - after root node cuts
2032   3 - after cuts at other nodes
2033   4 - during cuts at other nodes
2034   8 added if previous heuristic in loop found solution
2035*/
2036bool CbcRounding::shouldHeurRun(int whereFrom)
2037{
2038  if (whereFrom != 4) {
2039    return CbcHeuristic::shouldHeurRun(whereFrom);
2040  } else {
2041    numCouldRun_++;
2042    return shouldHeurRun_randomChoice();
2043  }
2044}
2045// See if rounding will give solution
2046// Sets value of solution
2047// Assumes rhs for original matrix still okay
2048// At present only works with integers
2049// Fix values if asked for
2050// Returns 1 if solution, 0 if not
2051int CbcRounding::solution(double &solutionValue,
2052  double *betterSolution)
2053{
2054
2055  numCouldRun_++;
2056  // See if to do
2057  if (!when() || (when() % 10 == 1 && model_->phase() != 1) || (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
2058    return 0; // switched off
2059  numRuns_++;
2060#ifdef HEURISTIC_INFORM
2061  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
2062    heuristicName(), numRuns_, numCouldRun_, when_);
2063#endif
2064  OsiSolverInterface *solver = model_->solver();
2065  double direction = solver->getObjSense();
2066  double newSolutionValue = direction * solver->getObjValue();
2067  return solution(solutionValue, betterSolution, newSolutionValue);
2068}
2069// See if rounding will give solution
2070// Sets value of solution
2071// Assumes rhs for original matrix still okay
2072// At present only works with integers
2073// Fix values if asked for
2074// Returns 1 if solution, 0 if not
2075int CbcRounding::solution(double &solutionValue,
2076  double *betterSolution,
2077  double newSolutionValue)
2078{
2079
2080  // See if to do
2081  if (!when() || (when() % 10 == 1 && model_->phase() != 1) || (when() % 10 == 2 && (model_->phase() != 2 && model_->phase() != 3)))
2082    return 0; // switched off
2083  OsiSolverInterface *solver = model_->solver();
2084  const double *lower = solver->getColLower();
2085  const double *upper = solver->getColUpper();
2086  const double *rowLower = solver->getRowLower();
2087  const double *rowUpper = solver->getRowUpper();
2088  const double *solution = solver->getColSolution();
2089  const double *objective = solver->getObjCoefficients();
2090  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
2091  double primalTolerance;
2092  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
2093  double useTolerance = primalTolerance;
2094
2095  int numberRows = matrix_.getNumRows();
2096  assert(numberRows <= solver->getNumRows());
2097  if (numberRows == 0) {
2098    return 0;
2099  }
2100  int numberIntegers = model_->numberIntegers();
2101  const int *integerVariable = model_->integerVariable();
2102  int i;
2103  double direction = solver->getObjSense();
2104  //double newSolutionValue = direction*solver->getObjValue();
2105  int returnCode = 0;
2106  // Column copy
2107  const double *element = matrix_.getElements();
2108  const int *row = matrix_.getIndices();
2109  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
2110  const int *columnLength = matrix_.getVectorLengths();
2111  // Row copy
2112  const double *elementByRow = matrixByRow_.getElements();
2113  const int *column = matrixByRow_.getIndices();
2114  const CoinBigIndex *rowStart = matrixByRow_.getVectorStarts();
2115  const int *rowLength = matrixByRow_.getVectorLengths();
2116
2117  // Get solution array for heuristic solution
2118  int numberColumns = solver->getNumCols();
2119  double *newSolution = new double[numberColumns];
2120  memcpy(newSolution, solution, numberColumns * sizeof(double));
2121
2122  double *rowActivity = new double[numberRows];
2123  memset(rowActivity, 0, numberRows * sizeof(double));
2124  for (i = 0; i < numberColumns; i++) {
2125    CoinBigIndex j;
2126    double value = newSolution[i];
2127    if (value < lower[i]) {
2128      value = lower[i];
2129      newSolution[i] = value;
2130    } else if (value > upper[i]) {
2131      value = upper[i];
2132      newSolution[i] = value;
2133    }
2134    if (value) {
2135      for (j = columnStart[i];
2136           j < columnStart[i] + columnLength[i]; j++) {
2137        int iRow = row[j];
2138        rowActivity[iRow] += value * element[j];
2139      }
2140    }
2141  }
2142  // check was feasible - if not adjust (cleaning may move)
2143  for (i = 0; i < numberRows; i++) {
2144    if (rowActivity[i] < rowLower[i]) {
2145      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
2146      rowActivity[i] = rowLower[i];
2147    } else if (rowActivity[i] > rowUpper[i]) {
2148      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
2149      rowActivity[i] = rowUpper[i];
2150    }
2151  }
2152  for (i = 0; i < numberIntegers; i++) {
2153    int iColumn = integerVariable[i];
2154    double value = newSolution[iColumn];
2155    //double thisTolerance = integerTolerance;
2156    if (fabs(floor(value + 0.5) - value) > integerTolerance) {
2157      double below = floor(value);
2158      double newValue = newSolution[iColumn];
2159      double cost = direction * objective[iColumn];
2160      double move;
2161      if (cost > 0.0) {
2162        // try up
2163        move = 1.0 - (value - below);
2164      } else if (cost < 0.0) {
2165        // try down
2166        move = below - value;
2167      } else {
2168        // won't be able to move unless we can grab another variable
2169        double randomNumber = randomNumberGenerator_.randomDouble();
2170        // which way?
2171        if (randomNumber < 0.5)
2172          move = below - value;
2173        else
2174          move = 1.0 - (value - below);
2175      }
2176      newValue += move;
2177      newSolution[iColumn] = newValue;
2178      newSolutionValue += move * cost;
2179      CoinBigIndex j;
2180      for (j = columnStart[iColumn];
2181           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2182        int iRow = row[j];
2183        rowActivity[iRow] += move * element[j];
2184      }
2185    }
2186  }
2187
2188  double penalty = 0.0;
2189  // see if feasible - just using singletons
2190  for (i = 0; i < numberRows; i++) {
2191    double value = rowActivity[i];
2192    double thisInfeasibility = 0.0;
2193    if (value < rowLower[i] - primalTolerance)
2194      thisInfeasibility = value - rowLower[i];
2195    else if (value > rowUpper[i] + primalTolerance)
2196      thisInfeasibility = value - rowUpper[i];
2197    if (thisInfeasibility) {
2198      // See if there are any slacks I can use to fix up
2199      // maybe put in coding for multiple slacks?
2200      double bestCost = 1.0e50;
2201      CoinBigIndex k;
2202      int iBest = -1;
2203      double addCost = 0.0;
2204      double newValue = 0.0;
2205      double changeRowActivity = 0.0;
2206      double absInfeasibility = fabs(thisInfeasibility);
2207      for (k = rowStart[i]; k < rowStart[i] + rowLength[i]; k++) {
2208        int iColumn = column[k];
2209        // See if all elements help
2210        if (columnLength[iColumn] == 1) {
2211          double currentValue = newSolution[iColumn];
2212          double elementValue = elementByRow[k];
2213          double lowerValue = lower[iColumn];
2214          double upperValue = upper[iColumn];
2215          double gap = rowUpper[i] - rowLower[i];
2216          double absElement = fabs(elementValue);
2217          if (thisInfeasibility * elementValue > 0.0) {
2218            // we want to reduce
2219            if ((currentValue - lowerValue) * absElement >= absInfeasibility) {
2220              // possible - check if integer
2221              double distance = absInfeasibility / absElement;
2222              double thisCost = -direction * objective[iColumn] * distance;
2223              if (isHeuristicInteger(solver, iColumn)) {
2224                distance = ceil(distance - useTolerance);
2225                if (currentValue - distance >= lowerValue - useTolerance) {
2226                  if (absInfeasibility - distance * absElement < -gap - useTolerance)
2227                    thisCost = 1.0e100; // no good
2228                  else
2229                    thisCost = -direction * objective[iColumn] * distance;
2230                } else {
2231                  thisCost = 1.0e100; // no good
2232                }
2233              }
2234              if (thisCost < bestCost) {
2235                bestCost = thisCost;
2236                iBest = iColumn;
2237                addCost = thisCost;
2238                newValue = currentValue - distance;
2239                changeRowActivity = -distance * elementValue;
2240              }
2241            }
2242          } else {
2243            // we want to increase
2244            if ((upperValue - currentValue) * absElement >= absInfeasibility) {
2245              // possible - check if integer
2246              double distance = absInfeasibility / absElement;
2247              double thisCost = direction * objective[iColumn] * distance;
2248              if (isHeuristicInteger(solver, iColumn)) {
2249                distance = ceil(distance - 1.0e-7);
2250                assert(currentValue - distance <= upperValue + useTolerance);
2251                if (absInfeasibility - distance * absElement < -gap - useTolerance)
2252                  thisCost = 1.0e100; // no good
2253                else
2254                  thisCost = direction * objective[iColumn] * distance;
2255              }
2256              if (thisCost < bestCost) {
2257                bestCost = thisCost;
2258                iBest = iColumn;
2259                addCost = thisCost;
2260                newValue = currentValue + distance;
2261                changeRowActivity = distance * elementValue;
2262              }
2263            }
2264          }
2265        }
2266      }
2267      if (iBest >= 0) {
2268        /*printf("Infeasibility of %g on row %d cost %g\n",
2269                  thisInfeasibility,i,addCost);*/
2270        newSolution[iBest] = newValue;
2271        thisInfeasibility = 0.0;
2272        newSolutionValue += addCost;
2273        rowActivity[i] += changeRowActivity;
2274      }
2275      penalty += fabs(thisInfeasibility);
2276    }
2277  }
2278  if (penalty) {
2279    // see if feasible using any
2280    // first continuous
2281    double penaltyChange = 0.0;
2282    int iColumn;
2283    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2284      if (isHeuristicInteger(solver, iColumn))
2285        continue;
2286      double currentValue = newSolution[iColumn];
2287      double lowerValue = lower[iColumn];
2288      double upperValue = upper[iColumn];
2289      CoinBigIndex j;
2290      int anyBadDown = 0;
2291      int anyBadUp = 0;
2292      double upImprovement = 0.0;
2293      double downImprovement = 0.0;
2294      for (j = columnStart[iColumn];
2295           j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2296        int iRow = row[j];
2297        if (rowUpper[iRow] > rowLower[iRow]) {
2298          double value = element[j];
2299          if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2300            // infeasible above
2301            downImprovement += value;
2302            upImprovement -= value;
2303            if (value > 0.0)
2304              anyBadUp++;
2305            else
2306              anyBadDown++;
2307          } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2308            // feasible at ub
2309            if (value > 0.0) {
2310              upImprovement -= value;
2311              anyBadUp++;
2312            } else {
2313              downImprovement += value;
2314              anyBadDown++;
2315            }
2316          } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2317            // feasible in interior
2318          } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2319            // feasible at lb
2320            if (value < 0.0) {
2321              upImprovement += value;
2322              anyBadUp++;
2323            } else {
2324              downImprovement -= value;
2325              anyBadDown++;
2326            }
2327          } else {
2328            // infeasible below
2329            downImprovement -= value;
2330            upImprovement += value;
2331            if (value < 0.0)
2332              anyBadUp++;
2333            else
2334              anyBadDown++;
2335          }
2336        } else {
2337          // equality row
2338          double value = element[j];
2339          if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2340            // infeasible above
2341            downImprovement += value;
2342            upImprovement -= value;
2343            if (value > 0.0)
2344              anyBadUp++;
2345            else
2346              anyBadDown++;
2347          } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2348            // infeasible below
2349            downImprovement -= value;
2350            upImprovement += value;
2351            if (value < 0.0)
2352              anyBadUp++;
2353            else
2354              anyBadDown++;
2355          } else {
2356            // feasible - no good
2357            anyBadUp = -1;
2358            anyBadDown = -1;
2359            break;
2360          }
2361        }
2362      }
2363      // could change tests for anyBad
2364      if (anyBadUp)
2365        upImprovement = 0.0;
2366      if (anyBadDown)
2367        downImprovement = 0.0;
2368      double way = 0.0;
2369      double improvement = 0.0;
2370      if (downImprovement > 0.0 && currentValue > lowerValue) {
2371        way = -1.0;
2372        improvement = downImprovement;
2373      } else if (upImprovement > 0.0 && currentValue < upperValue) {
2374        way = 1.0;
2375        improvement = upImprovement;
2376      }
2377      if (way) {
2378        // can improve
2379        double distance;
2380        if (way > 0.0)
2381          distance = upperValue - currentValue;
2382        else
2383          distance = currentValue - lowerValue;
2384        for (j = columnStart[iColumn];
2385             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2386          int iRow = row[j];
2387          double value = element[j] * way;
2388          if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2389            // infeasible above
2390            assert(value < 0.0);
2391            double gap = rowActivity[iRow] - rowUpper[iRow];
2392            if (gap + value * distance < 0.0)
2393              distance = -gap / value;
2394          } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2395            // infeasible below
2396            assert(value > 0.0);
2397            double gap = rowActivity[iRow] - rowLower[iRow];
2398            if (gap + value * distance > 0.0)
2399              distance = -gap / value;
2400          } else {
2401            // feasible
2402            if (value > 0) {
2403              double gap = rowActivity[iRow] - rowUpper[iRow];
2404              if (gap + value * distance > 0.0)
2405                distance = -gap / value;
2406            } else {
2407              double gap = rowActivity[iRow] - rowLower[iRow];
2408              if (gap + value * distance < 0.0)
2409                distance = -gap / value;
2410            }
2411          }
2412        }
2413        //move
2414        penaltyChange += improvement * distance;
2415        distance *= way;
2416        newSolution[iColumn] += distance;
2417        newSolutionValue += direction * objective[iColumn] * distance;
2418        for (j = columnStart[iColumn];
2419             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2420          int iRow = row[j];
2421          double value = element[j];
2422          rowActivity[iRow] += distance * value;
2423        }
2424      }
2425    }
2426    // and now all if improving
2427    double lastChange = penaltyChange ? 1.0 : 0.0;
2428    int numberPasses = 0;
2429    while (lastChange > 1.0e-2 && numberPasses < 1000) {
2430      lastChange = 0;
2431      numberPasses++;
2432      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2433        bool isInteger = isHeuristicInteger(solver, iColumn);
2434        double currentValue = newSolution[iColumn];
2435        double lowerValue = lower[iColumn];
2436        double upperValue = upper[iColumn];
2437        CoinBigIndex j;
2438        int anyBadDown = 0;
2439        int anyBadUp = 0;
2440        double upImprovement = 0.0;
2441        double downImprovement = 0.0;
2442        for (j = columnStart[iColumn];
2443             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2444          int iRow = row[j];
2445          double value = element[j];
2446          if (isInteger) {
2447            if (value > 0.0) {
2448              if (rowActivity[iRow] + value > rowUpper[iRow] + primalTolerance)
2449                anyBadUp++;
2450              if (rowActivity[iRow] - value < rowLower[iRow] - primalTolerance)
2451                anyBadDown++;
2452            } else {
2453              if (rowActivity[iRow] - value > rowUpper[iRow] + primalTolerance)
2454                anyBadDown++;
2455              if (rowActivity[iRow] + value < rowLower[iRow] - primalTolerance)
2456                anyBadUp++;
2457            }
2458          }
2459          if (rowUpper[iRow] > rowLower[iRow]) {
2460            if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2461              // infeasible above
2462              downImprovement += value;
2463              upImprovement -= value;
2464              if (value > 0.0)
2465                anyBadUp++;
2466              else
2467                anyBadDown++;
2468            } else if (rowActivity[iRow] > rowUpper[iRow] - primalTolerance) {
2469              // feasible at ub
2470              if (value > 0.0) {
2471                upImprovement -= value;
2472                anyBadUp++;
2473              } else {
2474                downImprovement += value;
2475                anyBadDown++;
2476              }
2477            } else if (rowActivity[iRow] > rowLower[iRow] + primalTolerance) {
2478              // feasible in interior
2479            } else if (rowActivity[iRow] > rowLower[iRow] - primalTolerance) {
2480              // feasible at lb
2481              if (value < 0.0) {
2482                upImprovement += value;
2483                anyBadUp++;
2484              } else {
2485                downImprovement -= value;
2486                anyBadDown++;
2487              }
2488            } else {
2489              // infeasible below
2490              downImprovement -= value;
2491              upImprovement += value;
2492              if (value < 0.0)
2493                anyBadUp++;
2494              else
2495                anyBadDown++;
2496            }
2497          } else {
2498            // equality row
2499            if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2500              // infeasible above
2501              downImprovement += value;
2502              upImprovement -= value;
2503              if (value > 0.0)
2504                anyBadUp++;
2505              else
2506                anyBadDown++;
2507            } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2508              // infeasible below
2509              downImprovement -= value;
2510              upImprovement += value;
2511              if (value < 0.0)
2512                anyBadUp++;
2513              else
2514                anyBadDown++;
2515            } else {
2516              // feasible - no good
2517              anyBadUp = -1;
2518              anyBadDown = -1;
2519              break;
2520            }
2521          }
2522        }
2523        // could change tests for anyBad
2524        if (anyBadUp)
2525          upImprovement = 0.0;
2526        if (anyBadDown)
2527          downImprovement = 0.0;
2528        double way = 0.0;
2529        double improvement = 0.0;
2530        if (downImprovement > 0.0 && currentValue > lowerValue) {
2531          way = -1.0;
2532          improvement = downImprovement;
2533          if (isInteger && currentValue < lowerValue + 0.99)
2534            continue; // no good
2535        } else if (upImprovement > 0.0 && currentValue < upperValue) {
2536          way = 1.0;
2537          improvement = upImprovement;
2538          if (isInteger && currentValue > upperValue - 0.99)
2539            continue; // no good
2540        }
2541        if (way) {
2542          // can improve
2543          double distance = COIN_DBL_MAX;
2544          for (j = columnStart[iColumn];
2545               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2546            int iRow = row[j];
2547            double value = element[j] * way;
2548            if (rowActivity[iRow] > rowUpper[iRow] + primalTolerance) {
2549              // infeasible above
2550              assert(value < 0.0);
2551              double gap = rowActivity[iRow] - rowUpper[iRow];
2552              if (gap + value * distance < 0.0) {
2553                // If integer then has to move by 1
2554                if (!isInteger)
2555                  distance = -gap / value;
2556                else
2557                  distance = CoinMax(-gap / value, 1.0);
2558              }
2559            } else if (rowActivity[iRow] < rowLower[iRow] - primalTolerance) {
2560              // infeasible below
2561              assert(value > 0.0);
2562              double gap = rowActivity[iRow] - rowLower[iRow];
2563              if (gap + value * distance > 0.0) {
2564                // If integer then has to move by 1
2565                if (!isInteger)
2566                  distance = -gap / value;
2567                else
2568                  distance = CoinMax(-gap / value, 1.0);
2569              }
2570            } else {
2571              // feasible
2572              if (value > 0) {
2573                double gap = rowActivity[iRow] - rowUpper[iRow];
2574                if (gap + value * distance > 0.0)
2575                  distance = -gap / value;
2576              } else {
2577                double gap = rowActivity[iRow] - rowLower[iRow];
2578                if (gap + value * distance < 0.0)
2579                  distance = -gap / value;
2580              }
2581            }
2582          }
2583          if (isInteger)
2584            distance = floor(distance + 1.05e-8);
2585          if (!distance) {
2586            // should never happen
2587            //printf("zero distance in CbcRounding - debug\n");
2588          }
2589          //move
2590          lastChange += improvement * distance;
2591          distance *= way;
2592          newSolution[iColumn] += distance;
2593          newSolutionValue += direction * objective[iColumn] * distance;
2594          for (j = columnStart[iColumn];
2595               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2596            int iRow = row[j];
2597            double value = element[j];
2598            rowActivity[iRow] += distance * value;
2599          }
2600        }
2601      }
2602      penaltyChange += lastChange;
2603    }
2604    penalty -= penaltyChange;
2605    if (penalty < 1.0e-5 * fabs(penaltyChange)) {
2606      // recompute
2607      penalty = 0.0;
2608      for (i = 0; i < numberRows; i++) {
2609        double value = rowActivity[i];
2610        if (value < rowLower[i] - primalTolerance)
2611          penalty += rowLower[i] - value;
2612        else if (value > rowUpper[i] + primalTolerance)
2613          penalty += value - rowUpper[i];
2614      }
2615    }
2616    // but integer variables may have wandered
2617    for (i = 0; i < numberIntegers; i++) {
2618      int iColumn = integerVariable[i];
2619      double value = newSolution[iColumn];
2620      if (fabs(floor(value + 0.5) - value) > integerTolerance)
2621        penalty += fabs(floor(value + 0.5) - value);
2622    }
2623  }
2624
2625  // Could also set SOS (using random) and repeat
2626  if (!penalty) {
2627    // See if we can do better
2628    //seed_++;
2629    //CoinSeedRandom(seed_);
2630    // Random number between 0 and 1.
2631    double randomNumber = randomNumberGenerator_.randomDouble();
2632    int iPass;
2633    int start[2];
2634    int end[2];
2635    int iRandom = static_cast< int >(randomNumber * (static_cast< double >(numberIntegers)));
2636    start[0] = iRandom;
2637    end[0] = numberIntegers;
2638    start[1] = 0;
2639    end[1] = iRandom;
2640    for (iPass = 0; iPass < 2; iPass++) {
2641      int i;
2642      for (i = start[iPass]; i < end[iPass]; i++) {
2643        int iColumn = integerVariable[i];
2644#ifndef NDEBUG
2645        double value = newSolution[iColumn];
2646        assert(fabs(floor(value + 0.5) - value) <= integerTolerance);
2647#endif
2648        double cost = direction * objective[iColumn];
2649        double move = 0.0;
2650        if (cost > 0.0)
2651          move = -1.0;
2652        else if (cost < 0.0)
2653          move = 1.0;
2654        while (move) {
2655          bool good = true;
2656          double newValue = newSolution[iColumn] + move;
2657          if (newValue < lower[iColumn] - useTolerance || newValue > upper[iColumn] + useTolerance) {
2658            move = 0.0;
2659          } else {
2660            // see if we can move
2661            CoinBigIndex j;
2662            for (j = columnStart[iColumn];
2663                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2664              int iRow = row[j];
2665              double newActivity = rowActivity[iRow] + move * element[j];
2666              if (newActivity < rowLower[iRow] - primalTolerance || newActivity > rowUpper[iRow] + primalTolerance) {
2667                good = false;
2668                break;
2669              }
2670            }
2671            if (good) {
2672              newSolution[iColumn] = newValue;
2673              newSolutionValue += move * cost;
2674              CoinBigIndex j;
2675              for (j = columnStart[iColumn];
2676                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2677                int iRow = row[j];
2678                rowActivity[iRow] += move * element[j];
2679              }
2680            } else {
2681              move = 0.0;
2682            }
2683          }
2684        }
2685      }
2686    }
2687    // Just in case of some stupidity
2688    double objOffset = 0.0;
2689    solver->getDblParam(OsiObjOffset, objOffset);
2690    newSolutionValue = -objOffset;
2691    for (i = 0; i < numberColumns; i++)
2692      newSolutionValue += objective[i] * newSolution[i];
2693    newSolutionValue *= direction;
2694    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
2695    if (newSolutionValue < solutionValue) {
2696      // paranoid check
2697      memset(rowActivity, 0, numberRows * sizeof(double));
2698      for (i = 0; i < numberColumns; i++) {
2699        CoinBigIndex j;
2700        double value = newSolution[i];
2701        if (value) {
2702          for (j = columnStart[i];
2703               j < columnStart[i] + columnLength[i]; j++) {
2704            int iRow = row[j];
2705            rowActivity[iRow] += value * element[j];
2706          }
2707        }
2708      }
2709      // check was approximately feasible
2710      bool feasible = true;
2711      for (i = 0; i < numberRows; i++) {
2712        if (rowActivity[i] < rowLower[i]) {
2713          if (rowActivity[i] < rowLower[i] - 1000.0 * primalTolerance)
2714            feasible = false;
2715        } else if (rowActivity[i] > rowUpper[i]) {
2716          if (rowActivity[i] > rowUpper[i] + 1000.0 * primalTolerance)
2717            feasible = false;
2718        }
2719      }
2720      if (feasible) {
2721        // new solution
2722        memcpy(betterSolution, newSolution, numberColumns * sizeof(double));
2723        solutionValue = newSolutionValue;
2724        //printf("** Solution of %g found by rounding\n",newSolutionValue);
2725        returnCode = 1;
2726      } else {
2727        // Can easily happen
2728        //printf("Debug CbcRounding giving bad solution\n");
2729      }
2730    }
2731  }
2732#ifdef NEW_ROUNDING
2733  if (!returnCode) {
2734#ifdef JJF_ZERO
2735    // back to starting point
2736    memcpy(newSolution, solution, numberColumns * sizeof(double));
2737    memset(rowActivity, 0, numberRows * sizeof(double));
2738    for (i = 0; i < numberColumns; i++) {
2739      int j;
2740      double value = newSolution[i];
2741      if (value < lower[i]) {
2742        value = lower[i];
2743        newSolution[i] = value;
2744      } else if (value > upper[i]) {
2745        value = upper[i];
2746        newSolution[i] = value;
2747      }
2748      if (value) {
2749        for (j = columnStart[i];
2750             j < columnStart[i] + columnLength[i]; j++) {
2751          int iRow = row[j];
2752          rowActivity[iRow] += value * element[j];
2753        }
2754      }
2755    }
2756    // check was feasible - if not adjust (cleaning may move)
2757    for (i = 0; i < numberRows; i++) {
2758      if (rowActivity[i] < rowLower[i]) {
2759        //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
2760        rowActivity[i] = rowLower[i];
2761      } else if (rowActivity[i] > rowUpper[i]) {
2762        //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
2763        rowActivity[i] = rowUpper[i];
2764      }
2765    }
2766#endif
2767    int *candidate = new int[numberColumns];
2768    int nCandidate = 0;
2769    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2770      bool isInteger = isHeuristicInteger(solver, iColumn);
2771      if (isInteger) {
2772        double currentValue = newSolution[iColumn];
2773        if (fabs(currentValue - floor(currentValue + 0.5)) > 1.0e-8)
2774          candidate[nCandidate++] = iColumn;
2775      }
2776    }
2777    if (true) {
2778      // Rounding as in Berthold
2779      while (nCandidate) {
2780        double infeasibility = 1.0e-7;
2781        int iRow = -1;
2782        for (i = 0; i < numberRows; i++) {
2783          double value = 0.0;
2784          if (rowActivity[i] < rowLower[i]) {
2785            value = rowLower[i] - rowActivity[i];
2786          } else if (rowActivity[i] > rowUpper[i]) {
2787            value = rowActivity[i] - rowUpper[i];
2788          }
2789          if (value > infeasibility) {
2790            infeasibility = value;
2791            iRow = i;
2792          }
2793        }
2794        if (iRow >= 0) {
2795          // infeasible
2796        } else {
2797          // feasible
2798        }
2799      }
2800    } else {
2801      // Shifting as in Berthold
2802    }
2803    delete[] candidate;
2804  }
2805#endif
2806  delete[] newSolution;
2807  delete[] rowActivity;
2808  return returnCode;
2809}
2810// update model
2811void CbcRounding::setModel(CbcModel *model)
2812{
2813  model_ = model;
2814  // Get a copy of original matrix (and by row for rounding);
2815  assert(model_->solver());
2816  if (model_->solver()->getNumRows()) {
2817    matrix_ = *model_->solver()->getMatrixByCol();
2818    matrixByRow_ = *model_->solver()->getMatrixByRow();
2819    // make sure model okay for heuristic
2820    validate();
2821  }
2822}
2823// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2824void CbcRounding::validate()
2825{
2826  if (model_ && (when() % 100) < 10) {
2827    if (model_->numberIntegers() != model_->numberObjects() && (model_->numberObjects() || (model_->specialOptions() & 1024) == 0)) {
2828      int numberOdd = 0;
2829      for (int i = 0; i < model_->numberObjects(); i++) {
2830        if (!model_->object(i)->canDoHeuristics())
2831          numberOdd++;
2832      }
2833      if (numberOdd)
2834        setWhen(0);
2835    }
2836  }
2837#ifdef NEW_ROUNDING
2838  int numberColumns = matrix_.getNumCols();
2839  down_ = new unsigned short[numberColumns];
2840  up_ = new unsigned short[numberColumns];
2841  equal_ = new unsigned short[numberColumns];
2842  // Column copy
2843  const double *element = matrix_.getElements();
2844  const int *row = matrix_.getIndices();
2845  const CoinBigIndex *columnStart = matrix_.getVectorStarts();
2846  const int *columnLength = matrix_.getVectorLengths();
2847  const double *rowLower = model.solver()->getRowLower();
2848  const double *rowUpper = model.solver()->getRowUpper();
2849  for (int i = 0; i < numberColumns; i++) {
2850    int down = 0;
2851    int up = 0;
2852    int equal = 0;
2853    if (columnLength[i] > 65535) {
2854      equal[0] = 65535;
2855      break; // unlikely to work
2856    }
2857    for (CoinBigIndex j = columnStart[i];
2858         j < columnStart[i] + columnLength[i]; j++) {
2859      int iRow = row[j];
2860      if (rowLower[iRow] > -1.0e20 && rowUpper[iRow] < 1.0e20) {
2861        equal++;
2862      } else if (element[j] > 0.0) {
2863        if (rowUpper[iRow] < 1.0e20)
2864          up++;
2865        else
2866          down--;
2867      } else {
2868        if (rowLower[iRow] > -1.0e20)
2869          up++;
2870        else
2871          down--;
2872      }
2873    }
2874    down_[i] = (unsigned short)down;
2875    up_[i] = (unsigned short)up;
2876    equal_[i] = (unsigned short)equal;
2877  }
2878#else
2879  down_ = NULL;
2880  up_ = NULL;
2881  equal_ = NULL;
2882#endif
2883}
2884
2885// Default Constructor
2886CbcHeuristicPartial::CbcHeuristicPartial()
2887  : CbcHeuristic()
2888{
2889  fixPriority_ = 10000;
2890}
2891
2892// Constructor from model
2893CbcHeuristicPartial::CbcHeuristicPartial(CbcModel &model, int fixPriority, int numberNodes)
2894  : CbcHeuristic(model)
2895{
2896  fixPriority_ = fixPriority;
2897  setNumberNodes(numberNodes);
2898  validate();
2899}
2900
2901// Destructor
2902CbcHeuristicPartial::~CbcHeuristicPartial()
2903{
2904}
2905
2906// Clone
2907CbcHeuristic *
2908CbcHeuristicPartial::clone() const
2909{
2910  return new CbcHeuristicPartial(*this);
2911}
2912// Create C++ lines to get to current state
2913void CbcHeuristicPartial::generateCpp(FILE *fp)
2914{
2915  CbcHeuristicPartial other;
2916  fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
2917  fprintf(fp, "3  CbcHeuristicPartial partial(*cbcModel);\n");
2918  CbcHeuristic::generateCpp(fp, "partial");
2919  if (fixPriority_ != other.fixPriority_)
2920    fprintf(fp, "3  partial.setFixPriority(%d);\n", fixPriority_);
2921  else
2922    fprintf(fp, "4  partial.setFixPriority(%d);\n", fixPriority_);
2923  fprintf(fp, "3  cbcModel->addHeuristic(&partial);\n");
2924}
2925//#define NEW_PARTIAL
2926// Copy constructor
2927CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial &rhs)
2928  : CbcHeuristic(rhs)
2929  , fixPriority_(rhs.fixPriority_)
2930{
2931}
2932
2933// Assignment operator
2934CbcHeuristicPartial &
2935CbcHeuristicPartial::operator=(const CbcHeuristicPartial &rhs)
2936{
2937  if (this != &rhs) {
2938    CbcHeuristic::operator=(rhs);
2939    fixPriority_ = rhs.fixPriority_;
2940  }
2941  return *this;
2942}
2943
2944// Resets stuff if model changes
2945void CbcHeuristicPartial::resetModel(CbcModel *model)
2946{
2947  model_ = model;
2948  // Get a copy of original matrix (and by row for partial);
2949  assert(model_->solver());
2950  validate();
2951}
2952// See if partial will give solution
2953// Sets value of solution
2954// Assumes rhs for original matrix still okay
2955// At present only works with integers
2956// Fix values if asked for
2957// Returns 1 if solution, 0 if not
2958int CbcHeuristicPartial::solution(double &solutionValue,
2959  double *betterSolution)
2960{
2961  // Return if already done
2962  if (fixPriority_ < 0)
2963    return 0; // switched off
2964#ifdef HEURISTIC_INFORM
2965  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
2966    heuristicName(), numRuns_, numCouldRun_, when_);
2967#endif
2968  const double *hotstartSolution = model_->hotstartSolution();
2969  const int *hotstartPriorities = model_->hotstartPriorities();
2970  if (!hotstartSolution)
2971    return 0;
2972  OsiSolverInterface *solver = model_->solver();
2973
2974  int numberIntegers = model_->numberIntegers();
2975  const int *integerVariable = model_->integerVariable();
2976
2977  OsiSolverInterface *newSolver = model_->continuousSolver()->clone();
2978  const double *colLower = newSolver->getColLower();
2979  const double *colUpper = newSolver->getColUpper();
2980
2981  double primalTolerance;
2982  solver->getDblParam(OsiPrimalTolerance, primalTolerance);
2983
2984  int i;
2985  int numberFixed = 0;
2986  int returnCode = 0;
2987
2988  for (i = 0; i < numberIntegers; i++) {
2989    int iColumn = integerVariable[i];
2990    if (abs(hotstartPriorities[iColumn]) <= fixPriority_) {
2991      double value = hotstartSolution[iColumn];
2992      double lower = colLower[iColumn];
2993      double upper = colUpper[iColumn];
2994      value = CoinMax(value, lower);
2995      value = CoinMin(value, upper);
2996      if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
2997        value = floor(value + 0.5);
2998        newSolver->setColLower(iColumn, value);
2999        newSolver->setColUpper(iColumn, value);
3000        numberFixed++;
3001      }
3002    }
3003  }
3004  if (numberFixed > numberIntegers / 5 - 100000000) {
3005#ifdef COIN_DEVELOP
3006    printf("%d integers fixed\n", numberFixed);
3007#endif
3008    returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
3009      model_->getCutoff(), "CbcHeuristicPartial");
3010    if (returnCode < 0)
3011      returnCode = 0; // returned on size
3012    //printf("return code %d",returnCode);
3013    if ((returnCode & 2) != 0) {
3014      // could add cut
3015      returnCode &= ~2;
3016      //printf("could add cut with %d elements (if all 0-1)\n",nFix);
3017    } else {
3018      //printf("\n");
3019    }
3020  }
3021  fixPriority_ = -1; // switch off
3022
3023  delete newSolver;
3024  return returnCode;
3025}
3026// update model
3027void CbcHeuristicPartial::setModel(CbcModel *model)
3028{
3029  model_ = model;
3030  assert(model_->solver());
3031  // make sure model okay for heuristic
3032  validate();
3033}
3034// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
3035void CbcHeuristicPartial::validate()
3036{
3037  if (model_ && (when() % 100) < 10) {
3038    if (model_->numberIntegers() != model_->numberObjects())
3039      setWhen(0);
3040  }
3041}
3042bool CbcHeuristicPartial::shouldHeurRun(int /*whereFrom*/)
3043{
3044  return true;
3045}
3046
3047// Default Constructor
3048CbcSerendipity::CbcSerendipity()
3049  : CbcHeuristic()
3050{
3051}
3052
3053// Constructor from model
3054CbcSerendipity::CbcSerendipity(CbcModel &model)
3055  : CbcHeuristic(model)
3056{
3057}
3058
3059// Destructor
3060CbcSerendipity::~CbcSerendipity()
3061{
3062}
3063
3064// Clone
3065CbcHeuristic *
3066CbcSerendipity::clone() const
3067{
3068  return new CbcSerendipity(*this);
3069}
3070// Create C++ lines to get to current state
3071void CbcSerendipity::generateCpp(FILE *fp)
3072{
3073  fprintf(fp, "0#include \"CbcHeuristic.hpp\"\n");
3074  fprintf(fp, "3  CbcSerendipity serendipity(*cbcModel);\n");
3075  CbcHeuristic::generateCpp(fp, "serendipity");
3076  fprintf(fp, "3  cbcModel->addHeuristic(&serendipity);\n");
3077}
3078
3079// Copy constructor
3080CbcSerendipity::CbcSerendipity(const CbcSerendipity &rhs)
3081  : CbcHeuristic(rhs)
3082{
3083}
3084
3085// Assignment operator
3086CbcSerendipity &
3087CbcSerendipity::operator=(const CbcSerendipity &rhs)
3088{
3089  if (this != &rhs) {
3090    CbcHeuristic::operator=(rhs);
3091  }
3092  return *this;
3093}
3094
3095// Returns 1 if solution, 0 if not
3096int CbcSerendipity::solution(double &solutionValue,
3097  double *betterSolution)
3098{
3099  if (!model_)
3100    return 0;
3101#ifdef HEURISTIC_INFORM
3102  printf("Entering heuristic %s - nRuns %d numCould %d when %d\n",
3103    heuristicName(), numRuns_, numCouldRun_, when_);
3104#endif
3105  if (!inputSolution_) {
3106    // get information on solver type
3107    OsiAuxInfo *auxInfo = model_->solver()->getAuxiliaryInfo();
3108    OsiBabSolver *auxiliaryInfo = dynamic_cast< OsiBabSolver * >(auxInfo);
3109    if (auxiliaryInfo) {
3110      return auxiliaryInfo->solution(solutionValue, betterSolution, model_->solver()->getNumCols());
3111    } else {
3112      return 0;
3113    }
3114  } else {
3115    int numberColumns = model_->getNumCols();
3116    double value = inputSolution_[numberColumns];
3117    int returnCode = 0;
3118    if (value < solutionValue) {
3119      solutionValue = value;
3120      memcpy(betterSolution, inputSolution_, numberColumns * sizeof(double));
3121      returnCode = 1;
3122    }
3123    delete[] inputSolution_;
3124    inputSolution_ = NULL;
3125    model_ = NULL; // switch off
3126    return returnCode;
3127  }
3128}
3129// update model
3130void CbcSerendipity::setModel(CbcModel *model)
3131{
3132  model_ = model;
3133}
3134// Resets stuff if model changes
3135void CbcSerendipity::resetModel(CbcModel *model)
3136{
3137  model_ = model;
3138}
3139
3140// Default Constructor
3141CbcHeuristicJustOne::CbcHeuristicJustOne()
3142  : CbcHeuristic()
3143  , probabilities_(NULL)
3144  , heuristic_(NULL)
3145  , numberHeuristics_(0)
3146{
3147}
3148
3149// Constructor from model
3150CbcHeuristicJustOne::CbcHeuristicJustOne(CbcModel &model)
3151  : CbcHeuristic(model)
3152  , probabilities_(NULL)
3153  , heuristic_(NULL)
3154  , numberHeuristics_(0)
3155{
3156}
3157
3158// Destructor
3159CbcHeuristicJustOne::~CbcHeuristicJustOne()
3160{
3161  for (int i = 0; i < numberHeuristics_; i++)
3162    delete heuristic_[i];
3163  delete[] heuristic_;
3164  delete[] probabilities_;
3165}
3166
3167// Clone
3168CbcHeuristicJustOne *
3169CbcHeuristicJustOne::clone() const
3170{
3171  return new CbcHeuristicJustOne(*this);
3172}
3173
3174// Create C++ lines to get to current state
3175void CbcHeuristicJustOne::generateCpp(FILE *fp)
3176{
3177  CbcHeuristicJustOne other;
3178  fprintf(fp, "0#include \"CbcHeuristicJustOne.hpp\"\n");
3179  fprintf(fp, "3  CbcHeuristicJustOne heuristicJustOne(*cbcModel);\n");
3180  CbcHeuristic::generateCpp(fp, "heuristicJustOne");
3181  fprintf(fp, "3  cbcModel->addHeuristic(&heuristicJustOne);\n");
3182}
3183
3184// Copy constructor
3185CbcHeuristicJustOne::CbcHeuristicJustOne(const CbcHeuristicJustOne &rhs)
3186  : CbcHeuristic(rhs)
3187  , probabilities_(NULL)
3188  , heuristic_(NULL)
3189  , numberHeuristics_(rhs.numberHeuristics_)
3190{
3191  if (numberHeuristics_) {
3192    probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
3193    heuristic_ = new CbcHeuristic *[numberHeuristics_];
3194    for (int i = 0; i < numberHeuristics_; i++)
3195      heuristic_[i] = rhs.heuristic_[i]->clone();
3196  }
3197}
3198
3199// Assignment operator
3200CbcHeuristicJustOne &
3201CbcHeuristicJustOne::operator=(const CbcHeuristicJustOne &rhs)
3202{
3203  if (this != &rhs) {
3204    CbcHeuristic::operator=(rhs);
3205    for (int i = 0; i < numberHeuristics_; i++)
3206      delete heuristic_[i];
3207    delete[] heuristic_;
3208    delete[] probabilities_;
3209    probabilities_ = NULL;
3210    heuristic_ = NULL;
3211    numberHeuristics_ = rhs.numberHeuristics_;
3212    if (numberHeuristics_) {
3213      probabilities_ = CoinCopyOfArray(rhs.probabilities_, numberHeuristics_);
3214      heuristic_ = new CbcHeuristic *[numberHeuristics_];
3215      for (int i = 0; i < numberHeuristics_; i++)
3216        heuristic_[i] = rhs.heuristic_[i]->clone();
3217    }
3218  }
3219  return *this;
3220}
3221// Sets value of solution
3222// Returns 1 if solution, 0 if not
3223int CbcHeuristicJustOne::solution(double &solutionValue,
3224  double *betterSolution)
3225{
3226#ifdef DIVE_DEBUG
3227  std::cout << "solutionValue = " << solutionValue << std::endl;
3228#endif
3229  ++numCouldRun_;
3230
3231  // test if the heuristic can run
3232  if (!shouldHeurRun_randomChoice() || !numberHeuristics_)
3233    return 0;
3234  double randomNumber = randomNumberGenerator_.randomDouble();
3235  int i;
3236  for (i = 0; i < numberHeuristics_; i++) {
3237    if (randomNumber < probabilities_[i])
3238      break;
3239  }
3240  assert(i < numberHeuristics_);
3241  int returnCode;
3242  //model_->unsetDivingHasRun();
3243#ifdef COIN_DEVELOP
3244  printf("JustOne running %s\n",
3245    heuristic_[i]->heuristicName());
3246#endif
3247  returnCode = heuristic_[i]->solution(solutionValue, betterSolution);
3248#ifdef COIN_DEVELOP
3249  if (returnCode)
3250    printf("JustOne running %s found solution\n",
3251      heuristic_[i]->heuristicName());
3252#endif
3253  return returnCode;
3254}
3255// Resets stuff if model changes
3256void CbcHeuristicJustOne::resetModel(CbcModel *model)
3257{
3258  CbcHeuristic::resetModel(model);
3259  for (int i = 0; i < numberHeuristics_; i++)
3260    heuristic_[i]->resetModel(model);
3261}
3262// update model (This is needed if cliques update matrix etc)
3263void CbcHeuristicJustOne::setModel(CbcModel *model)
3264{
3265  CbcHeuristic::setModel(model);
3266  for (int i = 0; i < numberHeuristics_; i++)
3267    heuristic_[i]->setModel(model);
3268}
3269// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
3270void CbcHeuristicJustOne::validate()
3271{
3272  CbcHeuristic::validate();
3273  for (int i = 0; i < numberHeuristics_; i++)
3274    heuristic_[i]->validate();
3275}
3276// Adds an heuristic with probability
3277void CbcHeuristicJustOne::addHeuristic(const CbcHeuristic *heuristic, double probability)
3278{
3279  CbcHeuristic *thisOne = heuristic->clone();
3280  thisOne->setWhen(-999);
3281  CbcHeuristic **tempH = CoinCopyOfArrayPartial(heuristic_, numberHeuristics_ + 1,
3282    numberHeuristics_);
3283  delete[] heuristic_;
3284  heuristic_ = tempH;
3285  heuristic_[numberHeuristics_] = thisOne;
3286  double *tempP = CoinCopyOfArrayPartial(probabilities_, numberHeuristics_ + 1,
3287    numberHeuristics_);
3288  delete[] probabilities_;
3289  probabilities_ = tempP;
3290  probabilities_[numberHeuristics_] = probability;
3291  numberHeuristics_++;
3292}
3293// Normalize probabilities
3294void CbcHeuristicJustOne::normalizeProbabilities()
3295{
3296  double sum = 0.0;
3297  for (int i = 0; i < numberHeuristics_; i++)
3298    sum += probabilities_[i];
3299  double multiplier = 1.0 / sum;
3300  sum = 0.0;
3301  for (int i = 0; i < numberHeuristics_; i++) {
3302    sum += probabilities_[i];
3303    probabilities_[i] = sum * multiplier;
3304  }
3305  assert(fabs(probabilities_[numberHeuristics_ - 1] - 1.0) < 1.0e-5);
3306  probabilities_[numberHeuristics_ - 1] = 1.000001;
3307}
3308
3309/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
3310*/
Note: See TracBrowser for help on using the repository browser.