source: stable/2.4/Cbc/src/CbcHeuristic.cpp @ 1271

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

Creating new stable branch 2.4 from trunk (rev 1270)

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