source: trunk/Cbc/src/CbcHeuristic.cpp @ 1132

Last change on this file since 1132 was 1132, checked in by forrest, 10 years ago

chnages to try and make faster

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