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

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

fixes

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