source: stable/2.3/Cbc/src/CbcHeuristic.cpp @ 1226

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

to fix a seg fault and bonmin

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