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

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

changes for heuristic clone

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