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

Last change on this file since 1053 was 1053, checked in by forrest, 12 years ago

trying to make go faster

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