source: branches/heur/Cbc/src/CbcHeuristic.cpp @ 891

Last change on this file since 891 was 891, checked in by jpgoncal, 13 years ago

Added gutsOfConstructor in CbcHeuristicNode?.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.6 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 <cmath>
12#include <cfloat>
13
14#ifdef COIN_HAS_CLP
15#include "OsiClpSolverInterface.hpp"
16#endif
17#include "CbcModel.hpp"
18#include "CbcMessage.hpp"
19#include "CbcHeuristic.hpp"
20#include "CbcHeuristicFPump.hpp"
21#include "CbcStrategy.hpp"
22#include "CglPreProcess.hpp"
23#include "OsiAuxInfo.hpp"
24#include "OsiPresolve.hpp"
25
26//==============================================================================
27
28CbcHeuristicNode::CbcHeuristicNode(const CbcHeuristicNode& rhs)
29{
30  numObjects_ = rhs.numObjects_;
31  brObj_ = new CbcBranchingObject*[numObjects_];
32  for (int i = 0; i < numObjects_; ++i) {
33    brObj_[i] = rhs.brObj_[i]->clone();
34  }
35}
36
37void
38CbcHeuristicNodeList::gutsOfDelete()
39{
40  for (int i = nodes_.size() - 1; i >= 0; --i) {
41    delete nodes_[i];
42  }
43}
44
45void
46CbcHeuristicNodeList::gutsOfCopy(const CbcHeuristicNodeList& rhs)
47{
48  append(rhs);
49}
50
51CbcHeuristicNodeList::CbcHeuristicNodeList(const CbcHeuristicNodeList& rhs)
52{
53  gutsOfCopy(rhs);
54}
55
56CbcHeuristicNodeList& CbcHeuristicNodeList::operator=
57(const CbcHeuristicNodeList& rhs)
58{
59  if (this != &rhs) {
60    gutsOfDelete();
61    gutsOfCopy(rhs);
62  }
63  return *this;
64}
65
66CbcHeuristicNodeList::~CbcHeuristicNodeList()
67{
68  gutsOfDelete();
69}
70
71void
72CbcHeuristicNodeList::append(CbcHeuristicNode*& node)
73{
74  nodes_.push_back(node);
75  node = NULL;
76}
77
78void
79CbcHeuristicNodeList::append(const CbcHeuristicNodeList& nodes)
80{
81  nodes_.reserve(nodes_.size() + nodes.size());
82  for (int i = 0; i < nodes.size(); ++i) {
83    CbcHeuristicNode* node = new CbcHeuristicNode(*nodes.node(i));
84    append(node);
85  }
86}
87
88//==============================================================================
89
90// Default Constructor
91CbcHeuristic::CbcHeuristic() :
92  model_(NULL),
93  when_(2),
94  numberNodes_(200),
95  feasibilityPumpOptions_(-1),
96  fractionSmall_(1.0),
97  heuristicName_("Unknown"),
98  howOften_(100),
99  decayFactor_(0.5),
100  shallowDepth_(0),
101  howOftenShallow_(100),
102  numInvocationsInShallow_(0),
103  numInvocationsInDeep_(0),
104  lastRunDeep_(0),
105  numRuns_(0),
106  minDistanceToRun_(4),
107  runNodes_()
108{
109  // As CbcHeuristic virtual need to modify .cpp if above change
110}
111
112// Constructor from model
113CbcHeuristic::CbcHeuristic(CbcModel & model) :
114  model_(&model),
115  when_(2),
116  numberNodes_(200),
117  feasibilityPumpOptions_(-1),
118  fractionSmall_(1.0),
119  heuristicName_("Unknown"),
120  howOften_(100),
121  decayFactor_(0.5),
122  shallowDepth_(0),
123  howOftenShallow_(100),
124  numInvocationsInShallow_(0),
125  numInvocationsInDeep_(0),
126  lastRunDeep_(0),
127  numRuns_(0),
128  minDistanceToRun_(4),
129  runNodes_()
130{}
131
132void
133CbcHeuristic::gutsOfCopy(const CbcHeuristic & rhs)
134{
135  model_ = rhs.model_;
136  when_ = rhs.when_;
137  numberNodes_ = rhs.numberNodes_;
138  feasibilityPumpOptions_ = rhs.feasibilityPumpOptions_;
139  fractionSmall_ = rhs.fractionSmall_;
140  randomNumberGenerator_ = rhs.randomNumberGenerator_;
141  heuristicName_ = rhs.heuristicName_;
142  howOften_ = rhs.howOften_;
143  decayFactor_ = rhs.howOften_;
144  shallowDepth_= rhs.shallowDepth_;
145  howOftenShallow_= rhs.howOftenShallow_;
146  numInvocationsInShallow_ = rhs.numInvocationsInShallow_;
147  numInvocationsInDeep_ = rhs.numInvocationsInDeep_;
148  lastRunDeep_ = rhs.lastRunDeep_;
149  numRuns_ = rhs.numRuns_;
150  minDistanceToRun_ = rhs.minDistanceToRun_;
151  runNodes_ = rhs.runNodes_;
152}
153// Copy constructor
154CbcHeuristic::CbcHeuristic(const CbcHeuristic & rhs)
155{
156  gutsOfCopy(rhs);
157}
158
159// Assignment operator
160CbcHeuristic & 
161CbcHeuristic::operator=( const CbcHeuristic& rhs)
162{
163  if (this!=&rhs) {
164    gutsOfDelete();
165    gutsOfCopy(rhs);
166  }
167  return *this;
168}
169
170bool
171CbcHeuristic::shouldHeurRun()
172{
173  const CbcNode* currentNode = model_->currentNode();
174  if (currentNode == NULL) {
175    return false;
176  }
177
178  const int depth = currentNode->depth();
179
180  const int nodeCount = model_->getNodeCount();  // FIXME: check that this is
181                                                 // correct in parallel
182
183  if (nodeCount == 0 || depth <= shallowDepth_) {
184    // what to do when we are in the shallow part of the tree
185    if (model_->getCurrentPassNumber() == 1) {
186      // first time in the node...
187      numInvocationsInShallow_ = 0;
188    }
189    ++numInvocationsInShallow_;
190    // Very large howOftenShallow_ will give the original test:
191    // (model_->getCurrentPassNumber() != 1)
192    if ((numInvocationsInShallow_ % howOftenShallow_) != 1) {
193      return false;
194    }
195    // LL: should we save these nodes in the list of nodes where the heur was
196    // LL: run?
197  } else {
198    // deeper in the tree
199    if (model_->getCurrentPassNumber() == 1) {
200      // first time in the node...
201      ++numInvocationsInDeep_;
202    }
203    if (numInvocationsInDeep_ - lastRunDeep_ < howOften_) {
204      return false;
205    }
206    if (model_->getCurrentPassNumber() != 1) {
207      // Run the heuristic only when first entering the node.
208      // LL: I don't think this is right. It should run just before strong
209      // LL: branching, I believe.
210      return false;
211    }
212    // Get where we are and create the appropriate CbcHeuristicNode object
213    CbcHeuristicNode* nodeDesc = new CbcHeuristicNode(*model_);
214    if (nodeDesc->minDistance(runNodes_) < minDistanceToRun_) {
215      delete nodeDesc;
216      return false;
217    }
218    runNodes_.append(nodeDesc);
219    ++lastRunDeep_;
220  }
221  ++numRuns_;
222  return true;
223}
224
225
226// Resets stuff if model changes
227void 
228CbcHeuristic::resetModel(CbcModel * model)
229{
230  model_=model;
231  }
232// Set seed
233void
234CbcHeuristic::setSeed(int value)
235{
236  randomNumberGenerator_.setSeed(value);
237}
238
239// Create C++ lines to get to current state
240void 
241CbcHeuristic::generateCpp( FILE * fp, const char * heuristic) 
242{
243  // hard coded as CbcHeuristic virtual
244  if (when_!=2)
245    fprintf(fp,"3  %s.setWhen(%d);\n",heuristic,when_);
246  else
247    fprintf(fp,"4  %s.setWhen(%d);\n",heuristic,when_);
248  if (numberNodes_!=200)
249    fprintf(fp,"3  %s.setNumberNodes(%d);\n",heuristic,numberNodes_);
250  else
251    fprintf(fp,"4  %s.setNumberNodes(%d);\n",heuristic,numberNodes_);
252  if (fractionSmall_!=1.0)
253    fprintf(fp,"3  %s.setFractionSmall(%g);\n",heuristic,fractionSmall_);
254  else
255    fprintf(fp,"4  %s.setFractionSmall(%g);\n",heuristic,fractionSmall_);
256  if (heuristicName_ != "Unknown")
257    fprintf(fp,"3  %s.setHeuristicName(\"%s\");\n",
258            heuristic,heuristicName_.c_str()) ;
259  else
260    fprintf(fp,"4  %s.setHeuristicName(\"%s\");\n",
261            heuristic,heuristicName_.c_str()) ;
262}
263// Destructor
264CbcHeuristic::~CbcHeuristic ()
265{
266}
267
268// update model
269void CbcHeuristic::setModel(CbcModel * model)
270{
271  model_ = model;
272}
273#ifdef COIN_DEVELOP
274extern bool getHistoryStatistics_;
275#endif
276// Do mini branch and bound (return 1 if solution)
277int 
278CbcHeuristic::smallBranchAndBound(OsiSolverInterface * solver,int numberNodes,
279                                  double * newSolution, double & newSolutionValue,
280                                  double cutoff, std::string name) const
281{
282#ifdef COIN_HAS_CLP
283  OsiClpSolverInterface * osiclp = dynamic_cast< OsiClpSolverInterface*> (solver);
284  if (osiclp&&(osiclp->specialOptions()&65536)==0) {
285    // go faster stripes
286    if (osiclp->getNumRows()<300&&osiclp->getNumCols()<500) {
287      osiclp->setupForRepeatedUse(2,0);
288    } else {
289      osiclp->setupForRepeatedUse(0,0);
290    }
291    // Turn this off if you get problems
292    // Used to be automatically set
293    osiclp->setSpecialOptions(osiclp->specialOptions()|(128+64));
294    ClpSimplex * lpSolver = osiclp->getModelPtr();
295    lpSolver->setSpecialOptions(lpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
296  }
297#endif
298#ifdef COIN_DEVELOP
299  getHistoryStatistics_=false;
300#endif
301  int status=0;
302  int logLevel = model_->logLevel();
303#define LEN_PRINT 250
304  char generalPrint[LEN_PRINT];
305  // Do presolve to see if possible
306  int numberColumns = solver->getNumCols();
307  char * reset = NULL;
308  int returnCode=1;
309  {
310    int saveLogLevel = solver->messageHandler()->logLevel();
311    if (saveLogLevel==1)
312      solver->messageHandler()->setLogLevel(0);
313    OsiPresolve * pinfo = new OsiPresolve();
314    int presolveActions=0;
315    // Allow dual stuff on integers
316    presolveActions=1;
317    // Do not allow all +1 to be tampered with
318    //if (allPlusOnes)
319    //presolveActions |= 2;
320    // allow transfer of costs
321    // presolveActions |= 4;
322    pinfo->setPresolveActions(presolveActions);
323    OsiSolverInterface * presolvedModel = pinfo->presolvedModel(*solver,1.0e-8,true,2);
324    delete pinfo;
325    // see if too big
326    double before = 2*solver->getNumRows()+solver->getNumCols();
327    if (presolvedModel) {
328      int afterRows = presolvedModel->getNumRows();
329      int afterCols = presolvedModel->getNumCols();
330      delete presolvedModel;
331      double after = 2*afterRows+afterCols;
332      if (after>fractionSmall_*before&&after>300) {
333        // Need code to try again to compress further using used
334        const int * used =  model_->usedInSolution();
335        int maxUsed=0;
336        int iColumn;
337        const double * lower = solver->getColLower();
338        const double * upper = solver->getColUpper();
339        for (iColumn=0;iColumn<numberColumns;iColumn++) {
340          if (upper[iColumn]>lower[iColumn]) {
341            if (solver->isBinary(iColumn))
342              maxUsed = CoinMax(maxUsed,used[iColumn]);
343          }
344        }
345        if (maxUsed) {
346          reset = new char [numberColumns];
347          int nFix=0;
348          for (iColumn=0;iColumn<numberColumns;iColumn++) {
349            reset[iColumn]=0;
350            if (upper[iColumn]>lower[iColumn]) {
351              if (solver->isBinary(iColumn)&&used[iColumn]==maxUsed) {
352                bool setValue=true;
353                if (maxUsed==1) {
354                  double randomNumber = randomNumberGenerator_.randomDouble();
355                  if (randomNumber>0.3)
356                    setValue=false;
357                }
358                if (setValue) {
359                  reset[iColumn]=1;
360                  solver->setColLower(iColumn,1.0);
361                  nFix++;
362                }
363              }
364            }
365          }
366          pinfo = new OsiPresolve();
367          presolveActions=0;
368          // Allow dual stuff on integers
369          presolveActions=1;
370          // Do not allow all +1 to be tampered with
371          //if (allPlusOnes)
372          //presolveActions |= 2;
373          // allow transfer of costs
374          // presolveActions |= 4;
375          pinfo->setPresolveActions(presolveActions);
376          presolvedModel = pinfo->presolvedModel(*solver,1.0e-8,true,2);
377          delete pinfo;
378          if(presolvedModel) {
379            // see if too big
380            int afterRows2 = presolvedModel->getNumRows();
381            int afterCols2 = presolvedModel->getNumCols();
382            delete presolvedModel;
383            double after = 2*afterRows2+afterCols2;
384            if (after>fractionSmall_*before&&after>300) {
385              sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - still too large",
386                      solver->getNumRows(),solver->getNumCols(),
387                      afterRows,afterCols,nFix,afterRows2,afterCols2);
388            } else {
389              sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - %d fixed gives %d, %d - ok now",
390                      solver->getNumRows(),solver->getNumCols(),
391                      afterRows,afterCols,nFix,afterRows2,afterCols2);
392            }
393            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
394              << generalPrint
395              <<CoinMessageEol;
396          } else {
397            returnCode=-1; // infeasible
398          }
399        }
400      }
401    } else {
402      returnCode=-1; // infeasible
403    }
404    solver->messageHandler()->setLogLevel(saveLogLevel);
405  }
406  if (returnCode==-1) {
407    delete [] reset;
408#ifdef COIN_DEVELOP
409    getHistoryStatistics_=true;
410#endif
411    return returnCode;
412  }
413  // Reduce printout
414  solver->setHintParam(OsiDoReducePrint,true,OsiHintTry);
415  solver->setHintParam(OsiDoPresolveInInitial,false,OsiHintTry);
416  solver->setDblParam(OsiDualObjectiveLimit,cutoff*solver->getObjSense());
417  solver->initialSolve();
418  if (solver->isProvenOptimal()) {
419    CglPreProcess process;
420    /* Do not try and produce equality cliques and
421       do up to 2 passes */
422    if (logLevel<=1)
423      process.messageHandler()->setLogLevel(0);
424    OsiSolverInterface * solver2= process.preProcessNonDefault(*solver,false,2);
425    if (!solver2) {
426      if (logLevel>1)
427        printf("Pre-processing says infeasible\n");
428      returnCode=2; // so will be infeasible
429    } else {
430      // see if too big
431      double before = 2*solver->getNumRows()+solver->getNumCols();
432      double after = 2*solver2->getNumRows()+solver2->getNumCols();
433      if (after>fractionSmall_*before&&after>300) {
434        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns - too large",
435                solver->getNumRows(),solver->getNumCols(),
436                solver2->getNumRows(),solver2->getNumCols());
437        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
438          << generalPrint
439          <<CoinMessageEol;
440        returnCode = -1;
441      } else {
442        sprintf(generalPrint,"Full problem %d rows %d columns, reduced to %d rows %d columns",
443                solver->getNumRows(),solver->getNumCols(),
444                solver2->getNumRows(),solver2->getNumCols());
445        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
446          << generalPrint
447          <<CoinMessageEol;
448      }
449      if (returnCode==1) {
450        solver2->resolve();
451        CbcModel model(*solver2);
452        if (logLevel<=1)
453          model.setLogLevel(0);
454        else
455          model.setLogLevel(logLevel);
456        if (feasibilityPumpOptions_>=0) {
457          CbcHeuristicFPump heuristic4;
458          int pumpTune=feasibilityPumpOptions_;
459          if (pumpTune>0) {
460            /*
461              >=10000000 for using obj
462              >=1000000 use as accumulate switch
463              >=1000 use index+1 as number of large loops
464              >=100 use 0.05 objvalue as increment
465              >=10 use +0.1 objvalue for cutoff (add)
466              1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
467              4 and static continuous, 5 as 3 but no internal integers
468              6 as 3 but all slack basis!
469            */
470            double value = solver2->getObjSense()*solver2->getObjValue();
471            int w = pumpTune/10;
472            int c = w % 10;
473            w /= 10;
474            int i = w % 10;
475            w /= 10;
476            int r = w;
477            int accumulate = r/1000;
478            r -= 1000*accumulate;
479            if (accumulate>=10) {
480              int which = accumulate/10;
481              accumulate -= 10*which;
482              which--;
483              // weights and factors
484              double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
485              double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
486              heuristic4.setInitialWeight(weight[which]);
487              heuristic4.setWeightFactor(factor[which]);
488            }
489            // fake cutoff
490            if (c) {
491              double cutoff;
492              solver2->getDblParam(OsiDualObjectiveLimit,cutoff);
493              cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
494              heuristic4.setFakeCutoff(cutoff);
495            }
496            if (i||r) {
497              // also set increment
498              //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
499              double increment = 0.0;
500              heuristic4.setAbsoluteIncrement(increment);
501              heuristic4.setAccumulate(accumulate);
502              heuristic4.setMaximumRetries(r+1);
503            }
504            pumpTune = pumpTune%100;
505            if (pumpTune==6)
506              pumpTune =13;
507            heuristic4.setWhen(pumpTune+10);
508          }
509          heuristic4.setHeuristicName("feasibility pump");
510          model.addHeuristic(&heuristic4);
511        }
512        model.setCutoff(cutoff);
513        model.setMaximumNodes(numberNodes);
514        model.solver()->setHintParam(OsiDoReducePrint,true,OsiHintTry);
515        // Lightweight
516        CbcStrategyDefaultSubTree strategy(model_,true,5,1,0);
517        model.setStrategy(strategy);
518        model.solver()->setIntParam(OsiMaxNumIterationHotStart,10);
519        // Do search
520        if (logLevel>1)
521          model_->messageHandler()->message(CBC_START_SUB,model_->messages())
522            << name
523            << model.getMaximumNodes()
524            <<CoinMessageEol;
525        // probably faster to use a basis to get integer solutions
526        model.setSpecialOptions(2);
527#ifdef CBC_THREAD
528        if (model_->getNumberThreads()>0&&(model_->getThreadMode()&1)!=0) {
529          // See if at root node
530          bool atRoot = model_->getNodeCount()==0;
531          int passNumber = model_->getCurrentPassNumber();
532          if (atRoot&&passNumber==1)
533            model.setNumberThreads(model_->getNumberThreads());
534        }
535#endif
536        model.setMaximumCutPassesAtRoot(CoinMin(20,model_->getMaximumCutPassesAtRoot()));
537        model.setParentModel(*model_);
538        model.setOriginalColumns(process.originalColumns());
539        if (model.getNumCols()) {
540          setCutAndHeuristicOptions(model);
541          model.branchAndBound();
542        } else {
543          // empty model
544          model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
545        }
546        if (logLevel>1)
547          model_->messageHandler()->message(CBC_END_SUB,model_->messages())
548            << name
549            <<CoinMessageEol;
550        if (model.getMinimizationObjValue()<CoinMin(cutoff,1.0e30)) {
551          // solution
552          if (model.getNumCols())
553            returnCode=model.isProvenOptimal() ? 3 : 1;
554          else
555            returnCode=3;
556          // post process
557#ifdef COIN_HAS_CLP
558          OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
559          if (clpSolver) {
560            ClpSimplex * lpSolver = clpSolver->getModelPtr();
561            lpSolver->setSpecialOptions(lpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
562          }
563#endif
564          process.postProcess(*model.solver());
565          if (solver->isProvenOptimal()) {
566            // Solution now back in solver
567            int numberColumns = solver->getNumCols();
568            memcpy(newSolution,solver->getColSolution(),
569                   numberColumns*sizeof(double));
570            newSolutionValue = model.getMinimizationObjValue();
571          } else {
572            // odd - but no good
573            returnCode=0; // so will be infeasible
574          }
575        } else {
576        // no good
577          returnCode=model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
578        }
579        if (model.status()==5)
580          returnCode=-2; // stop
581        if (model.isProvenInfeasible())
582          status=1;
583        else if (model.isProvenOptimal())
584          status=2;
585      }
586    }
587  } else {
588    returnCode=2; // infeasible finished
589  }
590  model_->setLogLevel(logLevel);
591  if (reset) {
592    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
593      if (reset[iColumn])
594        solver->setColLower(iColumn,0.0);
595    }
596    delete [] reset;
597  }
598#ifdef COIN_DEVELOP
599  getHistoryStatistics_=true;
600  if (returnCode==1||returnCode==2) {
601    if (status==1)
602      printf("heuristic could add cut because infeasible (%s)\n",heuristicName_.c_str()); 
603    else if (status==2)
604      printf("heuristic could add cut because optimal (%s)\n",heuristicName_.c_str());
605  } 
606#endif
607  return returnCode;
608}
609
610//##############################################################################
611
612inline int compare3BranchingObjects(const CbcBranchingObject* br0,
613                                    const CbcBranchingObject* br1)
614{
615  const int t0 = br0->type();
616  const int t1 = br1->type();
617  if (t0 < t1) {
618    return -1;
619  }
620  if (t0 > t1) {
621    return 1;
622  }
623  return br0->compareOriginalObject(br1);
624}
625
626//==============================================================================
627
628inline bool compareBranchingObjects(const CbcBranchingObject* br0,
629                                    const CbcBranchingObject* br1)
630{
631  return compare3BranchingObjects(br0, br1) < 0;
632}
633
634//==============================================================================
635
636void
637CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
638{
639  CbcNode* node = model.currentNode();
640  int depth = node->depth();
641  numObjects_ = depth; //??
642  brObj_ = new CbcBranchingObject*[numObjects_];
643  CbcNodeInfo* nodeInfo = node->nodeInfo();
644  depth = 0;
645  while (nodeInfo->owner() != NULL) {
646    const OsiBranchingObject* osibr =
647      nodeInfo->owner()->branchingObject();
648    const CbcBranchingObject* cbcbr =
649      dynamic_cast<const CbcBranchingObject*>(osibr);
650    assert(cbcbr);
651    brObj_[depth++] = cbcbr->clone();
652    nodeInfo = nodeInfo->parent();
653  }
654  std::sort(brObj_, brObj_+depth, compareBranchingObjects);
655  int cnt = 0;
656  CbcBranchingObject* br;
657  for (int i = 1; i < depth; ++i) {
658    if (compare3BranchingObjects(brObj_[cnt], brObj_[i]) == 0) {
659      int comp = brObj_[cnt]->compareBranchingObject(brObj_[i], &br);
660      switch (comp) {
661      case CbcRangeSame: // the same range
662      case CbcRangeDisjoint: // disjoint decisions
663        // should not happen! we are on a chain!
664        abort();
665      case CbcRangeSubset: // brObj_[cnt] is a subset of brObj_[i]
666        delete brObj_[i];
667        break;
668      case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[cnt]
669        delete brObj_[cnt];
670        brObj_[cnt] = brObj_[i];
671        break;
672      case CbcRangeOverlap: // overlap
673        delete brObj_[i];
674        delete brObj_[cnt];
675        brObj_[cnt] = br;
676        break;
677      }
678      continue;
679    } else {
680      brObj_[++cnt] = brObj_[i];
681    }
682  }
683  numObjects_ = cnt + 1;
684}
685
686//==============================================================================
687
688CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
689{
690  gutsOfConstructor(model);
691}
692
693//==============================================================================
694
695double
696CbcHeuristicNode::distance(const CbcHeuristicNode* node) const 
697{
698  const double disjointWeight = 1;
699  const double overlapWeight = 0.4;
700  const double subsetWeight = 0.2;
701  int i = 0; 
702  int j = 0;
703  double dist = 0.0;
704  while( i < numObjects_ && j < node->numObjects_) {
705    CbcBranchingObject* br0 = brObj_[i];
706    const CbcBranchingObject* br1 = node->brObj_[j];
707    const int brComp = compare3BranchingObjects(br0, br1);
708    if (brComp < 0) {
709      dist += subsetWeight;
710      ++i;
711    }
712    else if (brComp > 0) {
713      dist += subsetWeight;
714      ++j;
715    }
716    else {
717      const int comp = br0->compareBranchingObject(br1, false);
718      switch (comp) {
719      case CbcRangeSame:
720        // do nothing
721        break;
722      case CbcRangeDisjoint: // disjoint decisions
723        dist += disjointWeight;
724        break;
725      case CbcRangeSubset: // subset one way or another
726      case CbcRangeSuperset:
727        dist += subsetWeight;
728        break;
729      case CbcRangeOverlap: // overlap
730        dist += overlapWeight;
731        break;
732      }
733      ++i;
734      ++j;
735    }
736  }
737  dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
738  return dist;
739}
740
741//==============================================================================
742
743CbcHeuristicNode::~CbcHeuristicNode()
744{
745  for (int i = 0; i < numObjects_; ++i) {
746    delete brObj_[i];
747  }
748  delete [] brObj_;
749}
750
751//==============================================================================
752
753double
754CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList)
755{
756  double minDist = COIN_DBL_MAX;
757  for (int i = nodeList.size() - 1; i >= 0; --i) {
758    minDist = CoinMin(minDist, distance(nodeList.node(i)));
759  }
760  return minDist;
761}
762
763//==============================================================================
764
765double
766CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList)
767{
768  if (nodeList.size() == 0) {
769    return COIN_DBL_MAX;
770  }
771  double sumDist = 0;
772  for (int i = nodeList.size() - 1; i >= 0; --i) {
773    sumDist += distance(nodeList.node(i));
774  }
775  return sumDist/nodeList.size();
776}
777
778//##############################################################################
779
780// Default Constructor
781CbcRounding::CbcRounding() 
782  :CbcHeuristic()
783{
784  // matrix and row copy will automatically be empty
785  seed_=1;
786  down_ = NULL;
787  up_ = NULL;
788  equal_ = NULL;
789}
790
791// Constructor from model
792CbcRounding::CbcRounding(CbcModel & model)
793  :CbcHeuristic(model)
794{
795  // Get a copy of original matrix (and by row for rounding);
796  assert(model.solver());
797  matrix_ = *model.solver()->getMatrixByCol();
798  matrixByRow_ = *model.solver()->getMatrixByRow();
799  validate();
800  seed_=1;
801}
802
803// Destructor
804CbcRounding::~CbcRounding ()
805{
806  delete [] down_;
807  delete [] up_;
808  delete [] equal_;
809}
810
811// Clone
812CbcHeuristic *
813CbcRounding::clone() const
814{
815  return new CbcRounding(*this);
816}
817// Create C++ lines to get to current state
818void 
819CbcRounding::generateCpp( FILE * fp) 
820{
821  CbcRounding other;
822  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
823  fprintf(fp,"3  CbcRounding rounding(*cbcModel);\n");
824  CbcHeuristic::generateCpp(fp,"rounding");
825  if (seed_!=other.seed_)
826    fprintf(fp,"3  rounding.setSeed(%d);\n",seed_);
827  else
828    fprintf(fp,"4  rounding.setSeed(%d);\n",seed_);
829  fprintf(fp,"3  cbcModel->addHeuristic(&rounding);\n");
830}
831//#define NEW_ROUNDING
832// Copy constructor
833CbcRounding::CbcRounding(const CbcRounding & rhs)
834:
835  CbcHeuristic(rhs),
836  matrix_(rhs.matrix_),
837  matrixByRow_(rhs.matrixByRow_),
838  seed_(rhs.seed_)
839{
840#ifdef NEW_ROUNDING
841  int numberColumns = matrix_.getNumCols();
842  down_ = CoinCopyOfArray(rhs.down_,numberColumns);
843  up_ = CoinCopyOfArray(rhs.up_,numberColumns);
844  equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
845#else
846  down_ = NULL;
847  up_ = NULL;
848  equal_ = NULL;
849#endif 
850}
851
852// Assignment operator
853CbcRounding & 
854CbcRounding::operator=( const CbcRounding& rhs)
855{
856  if (this!=&rhs) {
857    CbcHeuristic::operator=(rhs);
858    matrix_ = rhs.matrix_;
859    matrixByRow_ = rhs.matrixByRow_;
860#ifdef NEW_ROUNDING
861    delete [] down_;
862    delete [] up_;
863    delete [] equal_;
864    int numberColumns = matrix_.getNumCols();
865    down_ = CoinCopyOfArray(rhs.down_,numberColumns);
866    up_ = CoinCopyOfArray(rhs.up_,numberColumns);
867    equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
868#else
869    down_ = NULL;
870    up_ = NULL;
871    equal_ = NULL;
872#endif 
873    seed_ = rhs.seed_;
874  }
875  return *this;
876}
877
878// Resets stuff if model changes
879void 
880CbcRounding::resetModel(CbcModel * model)
881{
882  model_=model;
883  // Get a copy of original matrix (and by row for rounding);
884  assert(model_->solver());
885  matrix_ = *model_->solver()->getMatrixByCol();
886  matrixByRow_ = *model_->solver()->getMatrixByRow();
887  validate();
888}
889// See if rounding will give solution
890// Sets value of solution
891// Assumes rhs for original matrix still okay
892// At present only works with integers
893// Fix values if asked for
894// Returns 1 if solution, 0 if not
895int
896CbcRounding::solution(double & solutionValue,
897                      double * betterSolution)
898{
899
900  // See if to do
901  if (!when()||(when()%10==1&&model_->phase()!=1)||
902      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
903    return 0; // switched off
904  OsiSolverInterface * solver = model_->solver();
905  double direction = solver->getObjSense();
906  double newSolutionValue = direction*solver->getObjValue();
907  return solution(solutionValue,betterSolution,newSolutionValue);
908}
909// See if rounding will give solution
910// Sets value of solution
911// Assumes rhs for original matrix still okay
912// At present only works with integers
913// Fix values if asked for
914// Returns 1 if solution, 0 if not
915int
916CbcRounding::solution(double & solutionValue,
917                      double * betterSolution,
918                      double newSolutionValue)
919{
920
921  // See if to do
922  if (!when()||(when()%10==1&&model_->phase()!=1)||
923      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
924    return 0; // switched off
925  OsiSolverInterface * solver = model_->solver();
926  const double * lower = solver->getColLower();
927  const double * upper = solver->getColUpper();
928  const double * rowLower = solver->getRowLower();
929  const double * rowUpper = solver->getRowUpper();
930  const double * solution = solver->getColSolution();
931  const double * objective = solver->getObjCoefficients();
932  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
933  double primalTolerance;
934  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
935
936  int numberRows = matrix_.getNumRows();
937  assert (numberRows<=solver->getNumRows());
938  int numberIntegers = model_->numberIntegers();
939  const int * integerVariable = model_->integerVariable();
940  int i;
941  double direction = solver->getObjSense();
942  //double newSolutionValue = direction*solver->getObjValue();
943  int returnCode = 0;
944  // Column copy
945  const double * element = matrix_.getElements();
946  const int * row = matrix_.getIndices();
947  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
948  const int * columnLength = matrix_.getVectorLengths();
949  // Row copy
950  const double * elementByRow = matrixByRow_.getElements();
951  const int * column = matrixByRow_.getIndices();
952  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
953  const int * rowLength = matrixByRow_.getVectorLengths();
954
955  // Get solution array for heuristic solution
956  int numberColumns = solver->getNumCols();
957  double * newSolution = new double [numberColumns];
958  memcpy(newSolution,solution,numberColumns*sizeof(double));
959
960  double * rowActivity = new double[numberRows];
961  memset(rowActivity,0,numberRows*sizeof(double));
962  for (i=0;i<numberColumns;i++) {
963    int j;
964    double value = newSolution[i];
965    if (value<lower[i]) {
966      value=lower[i];
967      newSolution[i]=value;
968    } else if (value>upper[i]) {
969      value=upper[i];
970      newSolution[i]=value;
971    }
972    if (value) {
973      for (j=columnStart[i];
974           j<columnStart[i]+columnLength[i];j++) {
975        int iRow=row[j];
976        rowActivity[iRow] += value*element[j];
977      }
978    }
979  }
980  // check was feasible - if not adjust (cleaning may move)
981  for (i=0;i<numberRows;i++) {
982    if(rowActivity[i]<rowLower[i]) {
983      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
984      rowActivity[i]=rowLower[i];
985    } else if(rowActivity[i]>rowUpper[i]) {
986      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
987      rowActivity[i]=rowUpper[i];
988    }
989  }
990  for (i=0;i<numberIntegers;i++) {
991    int iColumn = integerVariable[i];
992    double value=newSolution[iColumn];
993    if (fabs(floor(value+0.5)-value)>integerTolerance) {
994      double below = floor(value);
995      double newValue=newSolution[iColumn];
996      double cost = direction * objective[iColumn];
997      double move;
998      if (cost>0.0) {
999        // try up
1000        move = 1.0 -(value-below);
1001      } else if (cost<0.0) {
1002        // try down
1003        move = below-value;
1004      } else {
1005        // won't be able to move unless we can grab another variable
1006        double randomNumber = randomNumberGenerator_.randomDouble();
1007        // which way?
1008        if (randomNumber<0.5) 
1009          move = below-value;
1010        else
1011          move = 1.0 -(value-below);
1012      }
1013      newValue += move;
1014      newSolution[iColumn] = newValue;
1015      newSolutionValue += move*cost;
1016      int j;
1017      for (j=columnStart[iColumn];
1018           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1019        int iRow = row[j];
1020        rowActivity[iRow] += move*element[j];
1021      }
1022    }
1023  }
1024
1025  double penalty=0.0;
1026  const char * integerType = model_->integerType();
1027  // see if feasible - just using singletons
1028  for (i=0;i<numberRows;i++) {
1029    double value = rowActivity[i];
1030    double thisInfeasibility=0.0;
1031    if (value<rowLower[i]-primalTolerance)
1032      thisInfeasibility = value-rowLower[i];
1033    else if (value>rowUpper[i]+primalTolerance)
1034      thisInfeasibility = value-rowUpper[i];
1035    if (thisInfeasibility) {
1036      // See if there are any slacks I can use to fix up
1037      // maybe put in coding for multiple slacks?
1038      double bestCost = 1.0e50;
1039      int k;
1040      int iBest=-1;
1041      double addCost=0.0;
1042      double newValue=0.0;
1043      double changeRowActivity=0.0;
1044      double absInfeasibility = fabs(thisInfeasibility);
1045      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
1046        int iColumn = column[k];
1047        // See if all elements help
1048        if (columnLength[iColumn]==1) {
1049          double currentValue = newSolution[iColumn];
1050          double elementValue = elementByRow[k];
1051          double lowerValue = lower[iColumn];
1052          double upperValue = upper[iColumn];
1053          double gap = rowUpper[i]-rowLower[i];
1054          double absElement=fabs(elementValue);
1055          if (thisInfeasibility*elementValue>0.0) {
1056            // we want to reduce
1057            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
1058              // possible - check if integer
1059              double distance = absInfeasibility/absElement;
1060              double thisCost = -direction*objective[iColumn]*distance;
1061              if (integerType[iColumn]) {
1062                distance = ceil(distance-primalTolerance);
1063                if (currentValue-distance>=lowerValue-primalTolerance) {
1064                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
1065                    thisCost=1.0e100; // no good
1066                  else
1067                    thisCost = -direction*objective[iColumn]*distance;
1068                } else {
1069                  thisCost=1.0e100; // no good
1070                }
1071              }
1072              if (thisCost<bestCost) {
1073                bestCost=thisCost;
1074                iBest=iColumn;
1075                addCost = thisCost;
1076                newValue = currentValue-distance;
1077                changeRowActivity = -distance*elementValue;
1078              }
1079            }
1080          } else {
1081            // we want to increase
1082            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
1083              // possible - check if integer
1084              double distance = absInfeasibility/absElement;
1085              double thisCost = direction*objective[iColumn]*distance;
1086              if (integerType[iColumn]) {
1087                distance = ceil(distance-1.0e-7);
1088                assert (currentValue-distance<=upperValue+primalTolerance);
1089                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
1090                  thisCost=1.0e100; // no good
1091                else
1092                  thisCost = direction*objective[iColumn]*distance;
1093              }
1094              if (thisCost<bestCost) {
1095                bestCost=thisCost;
1096                iBest=iColumn;
1097                addCost = thisCost;
1098                newValue = currentValue+distance;
1099                changeRowActivity = distance*elementValue;
1100              }
1101            }
1102          }
1103        }
1104      }
1105      if (iBest>=0) {
1106        /*printf("Infeasibility of %g on row %d cost %g\n",
1107          thisInfeasibility,i,addCost);*/
1108        newSolution[iBest]=newValue;
1109        thisInfeasibility=0.0;
1110        newSolutionValue += addCost;
1111        rowActivity[i] += changeRowActivity;
1112      }
1113      penalty += fabs(thisInfeasibility);
1114    }
1115  }
1116  if (penalty) {
1117    // see if feasible using any
1118    // first continuous
1119    double penaltyChange=0.0;
1120    int iColumn;
1121    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1122      if (integerType[iColumn])
1123        continue;
1124      double currentValue = newSolution[iColumn];
1125      double lowerValue = lower[iColumn];
1126      double upperValue = upper[iColumn];
1127      int j;
1128      int anyBadDown=0;
1129      int anyBadUp=0;
1130      double upImprovement=0.0;
1131      double downImprovement=0.0;
1132      for (j=columnStart[iColumn];
1133           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1134        int iRow = row[j];
1135        if (rowUpper[iRow]>rowLower[iRow]) {
1136          double value = element[j];
1137          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1138            // infeasible above
1139            downImprovement += value;
1140            upImprovement -= value;
1141            if (value>0.0) 
1142              anyBadUp++;
1143            else 
1144              anyBadDown++;
1145          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
1146            // feasible at ub
1147            if (value>0.0) {
1148              upImprovement -= value;
1149              anyBadUp++;
1150            } else {
1151              downImprovement += value;
1152              anyBadDown++;
1153            }
1154          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
1155            // feasible in interior
1156          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
1157            // feasible at lb
1158            if (value<0.0) {
1159              upImprovement += value;
1160              anyBadUp++;
1161            } else {
1162              downImprovement -= value;
1163              anyBadDown++;
1164            }
1165          } else {
1166            // infeasible below
1167            downImprovement -= value;
1168            upImprovement += value;
1169            if (value<0.0) 
1170              anyBadUp++;
1171            else 
1172              anyBadDown++;
1173          }
1174        } else {
1175          // equality row
1176          double value = element[j];
1177          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1178            // infeasible above
1179            downImprovement += value;
1180            upImprovement -= value;
1181            if (value>0.0) 
1182              anyBadUp++;
1183            else 
1184              anyBadDown++;
1185          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1186            // infeasible below
1187            downImprovement -= value;
1188            upImprovement += value;
1189            if (value<0.0) 
1190              anyBadUp++;
1191            else 
1192              anyBadDown++;
1193          } else {
1194            // feasible - no good
1195            anyBadUp=-1;
1196            anyBadDown=-1;
1197            break;
1198          }
1199        }
1200      }
1201      // could change tests for anyBad
1202      if (anyBadUp)
1203        upImprovement=0.0;
1204      if (anyBadDown)
1205        downImprovement=0.0;
1206      double way=0.0;
1207      double improvement=0.0;
1208      if (downImprovement>0.0&&currentValue>lowerValue) {
1209        way=-1.0;
1210        improvement = downImprovement;
1211      } else if (upImprovement>0.0&&currentValue<upperValue) {
1212        way=1.0;
1213        improvement = upImprovement;
1214      }
1215      if (way) {
1216        // can improve
1217        double distance;
1218        if (way>0.0)
1219          distance = upperValue-currentValue;
1220        else
1221          distance = currentValue-lowerValue;
1222        for (j=columnStart[iColumn];
1223             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1224          int iRow = row[j];
1225          double value = element[j]*way;
1226          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1227            // infeasible above
1228            assert (value<0.0);
1229            double gap = rowActivity[iRow]-rowUpper[iRow];
1230            if (gap+value*distance<0.0) 
1231              distance = -gap/value;
1232          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1233            // infeasible below
1234            assert (value>0.0);
1235            double gap = rowActivity[iRow]-rowLower[iRow];
1236            if (gap+value*distance>0.0) 
1237              distance = -gap/value;
1238          } else {
1239            // feasible
1240            if (value>0) {
1241              double gap = rowActivity[iRow]-rowUpper[iRow];
1242              if (gap+value*distance>0.0) 
1243              distance = -gap/value;
1244            } else {
1245              double gap = rowActivity[iRow]-rowLower[iRow];
1246              if (gap+value*distance<0.0) 
1247                distance = -gap/value;
1248            }
1249          }
1250        }
1251        //move
1252        penaltyChange += improvement*distance;
1253        distance *= way;
1254        newSolution[iColumn] += distance;
1255        newSolutionValue += direction*objective[iColumn]*distance;
1256        for (j=columnStart[iColumn];
1257             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1258          int iRow = row[j];
1259          double value = element[j];
1260          rowActivity[iRow] += distance*value;
1261        }
1262      }
1263    }
1264    // and now all if improving
1265    double lastChange= penaltyChange ? 1.0 : 0.0;
1266    while (lastChange>1.0e-2) {
1267      lastChange=0;
1268      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1269        bool isInteger = (integerType[iColumn]!=0);
1270        double currentValue = newSolution[iColumn];
1271        double lowerValue = lower[iColumn];
1272        double upperValue = upper[iColumn];
1273        int j;
1274        int anyBadDown=0;
1275        int anyBadUp=0;
1276        double upImprovement=0.0;
1277        double downImprovement=0.0;
1278        for (j=columnStart[iColumn];
1279             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1280          int iRow = row[j];
1281          double value = element[j];
1282          if (isInteger) {
1283            if (value>0.0) {
1284              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
1285                anyBadUp++;
1286              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
1287                anyBadDown++;
1288            } else {
1289              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
1290                anyBadDown++;
1291              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
1292                anyBadUp++;
1293            }
1294          }
1295          if (rowUpper[iRow]>rowLower[iRow]) {
1296            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1297              // infeasible above
1298              downImprovement += value;
1299              upImprovement -= value;
1300              if (value>0.0) 
1301                anyBadUp++;
1302              else 
1303                anyBadDown++;
1304            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
1305              // feasible at ub
1306              if (value>0.0) {
1307                upImprovement -= value;
1308                anyBadUp++;
1309              } else {
1310                downImprovement += value;
1311                anyBadDown++;
1312              }
1313            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
1314              // feasible in interior
1315            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
1316              // feasible at lb
1317              if (value<0.0) {
1318                upImprovement += value;
1319                anyBadUp++;
1320              } else {
1321                downImprovement -= value;
1322                anyBadDown++;
1323              }
1324            } else {
1325              // infeasible below
1326              downImprovement -= value;
1327              upImprovement += value;
1328              if (value<0.0) 
1329                anyBadUp++;
1330              else 
1331                anyBadDown++;
1332            }
1333          } else {
1334            // equality row
1335            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1336              // infeasible above
1337              downImprovement += value;
1338              upImprovement -= value;
1339              if (value>0.0) 
1340                anyBadUp++;
1341              else 
1342                anyBadDown++;
1343            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1344              // infeasible below
1345              downImprovement -= value;
1346              upImprovement += value;
1347              if (value<0.0) 
1348                anyBadUp++;
1349              else 
1350                anyBadDown++;
1351            } else {
1352              // feasible - no good
1353              anyBadUp=-1;
1354              anyBadDown=-1;
1355              break;
1356            }
1357          }
1358        }
1359        // could change tests for anyBad
1360        if (anyBadUp)
1361          upImprovement=0.0;
1362        if (anyBadDown)
1363          downImprovement=0.0;
1364        double way=0.0;
1365        double improvement=0.0;
1366        if (downImprovement>0.0&&currentValue>lowerValue) {
1367          way=-1.0;
1368          improvement = downImprovement;
1369        } else if (upImprovement>0.0&&currentValue<upperValue) {
1370          way=1.0;
1371          improvement = upImprovement;
1372        }
1373        if (way) {
1374          // can improve
1375          double distance=COIN_DBL_MAX;
1376          for (j=columnStart[iColumn];
1377               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1378            int iRow = row[j];
1379            double value = element[j]*way;
1380            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1381              // infeasible above
1382              assert (value<0.0);
1383              double gap = rowActivity[iRow]-rowUpper[iRow];
1384              if (gap+value*distance<0.0) {
1385                // If integer then has to move by 1
1386                if (!isInteger)
1387                  distance = -gap/value;
1388                else
1389                  distance = CoinMax(-gap/value,1.0);
1390              }
1391            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1392              // infeasible below
1393              assert (value>0.0);
1394              double gap = rowActivity[iRow]-rowLower[iRow];
1395              if (gap+value*distance>0.0) {
1396                // If integer then has to move by 1
1397                if (!isInteger)
1398                  distance = -gap/value;
1399                else
1400                  distance = CoinMax(-gap/value,1.0);
1401              }
1402            } else {
1403              // feasible
1404              if (value>0) {
1405                double gap = rowActivity[iRow]-rowUpper[iRow];
1406                if (gap+value*distance>0.0) 
1407                  distance = -gap/value;
1408              } else {
1409                double gap = rowActivity[iRow]-rowLower[iRow];
1410                if (gap+value*distance<0.0) 
1411                  distance = -gap/value;
1412              }
1413            }
1414          }
1415          if (isInteger)
1416            distance = floor(distance+1.05e-8);
1417          if (!distance) {
1418            // should never happen
1419            //printf("zero distance in CbcRounding - debug\n");
1420          }
1421          //move
1422          lastChange += improvement*distance;
1423          distance *= way;
1424          newSolution[iColumn] += distance;
1425          newSolutionValue += direction*objective[iColumn]*distance;
1426          for (j=columnStart[iColumn];
1427               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1428            int iRow = row[j];
1429            double value = element[j];
1430            rowActivity[iRow] += distance*value;
1431          }
1432        }
1433      }
1434      penaltyChange += lastChange;
1435    }
1436    penalty -= penaltyChange;
1437    if (penalty<1.0e-5*fabs(penaltyChange)) {
1438      // recompute
1439      penalty=0.0;
1440      for (i=0;i<numberRows;i++) {
1441        double value = rowActivity[i];
1442        if (value<rowLower[i]-primalTolerance)
1443          penalty += rowLower[i]-value;
1444        else if (value>rowUpper[i]+primalTolerance)
1445          penalty += value-rowUpper[i];
1446      }
1447    }
1448  }
1449
1450  // Could also set SOS (using random) and repeat
1451  if (!penalty) {
1452    // See if we can do better
1453    //seed_++;
1454    //CoinSeedRandom(seed_);
1455    // Random number between 0 and 1.
1456    double randomNumber = randomNumberGenerator_.randomDouble();
1457    int iPass;
1458    int start[2];
1459    int end[2];
1460    int iRandom = (int) (randomNumber*((double) numberIntegers));
1461    start[0]=iRandom;
1462    end[0]=numberIntegers;
1463    start[1]=0;
1464    end[1]=iRandom;
1465    for (iPass=0;iPass<2;iPass++) {
1466      int i;
1467      for (i=start[iPass];i<end[iPass];i++) {
1468        int iColumn = integerVariable[i];
1469#ifndef NDEBUG
1470        double value=newSolution[iColumn];
1471        assert (fabs(floor(value+0.5)-value)<integerTolerance);
1472#endif
1473        double cost = direction * objective[iColumn];
1474        double move=0.0;
1475        if (cost>0.0)
1476          move = -1.0;
1477        else if (cost<0.0)
1478          move=1.0;
1479        while (move) {
1480          bool good=true;
1481          double newValue=newSolution[iColumn]+move;
1482          if (newValue<lower[iColumn]-primalTolerance||
1483              newValue>upper[iColumn]+primalTolerance) {
1484            move=0.0;
1485          } else {
1486            // see if we can move
1487            int j;
1488            for (j=columnStart[iColumn];
1489                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1490              int iRow = row[j];
1491              double newActivity = rowActivity[iRow] + move*element[j];
1492              if (newActivity<rowLower[iRow]-primalTolerance||
1493                  newActivity>rowUpper[iRow]+primalTolerance) {
1494                good=false;
1495                break;
1496              }
1497            }
1498            if (good) {
1499              newSolution[iColumn] = newValue;
1500              newSolutionValue += move*cost;
1501              int j;
1502              for (j=columnStart[iColumn];
1503                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
1504                int iRow = row[j];
1505                rowActivity[iRow] += move*element[j];
1506              }
1507            } else {
1508              move=0.0;
1509            }
1510          }
1511        }
1512      }
1513    }
1514    // Just in case of some stupidity
1515    double objOffset=0.0;
1516    solver->getDblParam(OsiObjOffset,objOffset);
1517    newSolutionValue = -objOffset;
1518    for ( i=0 ; i<numberColumns ; i++ )
1519      newSolutionValue += objective[i]*newSolution[i];
1520    newSolutionValue *= direction;
1521    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
1522    if (newSolutionValue<solutionValue) {
1523      // paranoid check
1524      memset(rowActivity,0,numberRows*sizeof(double));
1525      for (i=0;i<numberColumns;i++) {
1526        int j;
1527        double value = newSolution[i];
1528        if (value) {
1529          for (j=columnStart[i];
1530               j<columnStart[i]+columnLength[i];j++) {
1531            int iRow=row[j];
1532            rowActivity[iRow] += value*element[j];
1533          }
1534        }
1535      }
1536      // check was approximately feasible
1537      bool feasible=true;
1538      for (i=0;i<numberRows;i++) {
1539        if(rowActivity[i]<rowLower[i]) {
1540          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
1541            feasible = false;
1542        } else if(rowActivity[i]>rowUpper[i]) {
1543          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
1544            feasible = false;
1545        }
1546      }
1547      if (feasible) {
1548        // new solution
1549        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1550        solutionValue = newSolutionValue;
1551        //printf("** Solution of %g found by rounding\n",newSolutionValue);
1552        returnCode=1;
1553      } else {
1554        // Can easily happen
1555        //printf("Debug CbcRounding giving bad solution\n");
1556      }
1557    }
1558  }
1559#ifdef NEW_ROUNDING
1560  if (!returnCode) {
1561#if 0
1562    // back to starting point
1563    memcpy(newSolution,solution,numberColumns*sizeof(double));
1564    memset(rowActivity,0,numberRows*sizeof(double));
1565    for (i=0;i<numberColumns;i++) {
1566      int j;
1567      double value = newSolution[i];
1568      if (value<lower[i]) {
1569        value=lower[i];
1570        newSolution[i]=value;
1571      } else if (value>upper[i]) {
1572        value=upper[i];
1573        newSolution[i]=value;
1574      }
1575      if (value) {
1576        for (j=columnStart[i];
1577             j<columnStart[i]+columnLength[i];j++) {
1578          int iRow=row[j];
1579          rowActivity[iRow] += value*element[j];
1580        }
1581      }
1582    }
1583    // check was feasible - if not adjust (cleaning may move)
1584    for (i=0;i<numberRows;i++) {
1585      if(rowActivity[i]<rowLower[i]) {
1586        //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1587        rowActivity[i]=rowLower[i];
1588      } else if(rowActivity[i]>rowUpper[i]) {
1589        //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1590        rowActivity[i]=rowUpper[i];
1591      }
1592    }
1593#endif
1594    int * candidate = new int [numberColumns];
1595    int nCandidate=0;
1596    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1597      bool isInteger = (integerType[iColumn]!=0);
1598      if (isInteger) {
1599        double currentValue = newSolution[iColumn];
1600        if (fabs(currentValue-floor(currentValue+0.5))>1.0e-8)
1601          candidate[nCandidate++]=iColumn;
1602      }
1603    }
1604    if (true) {
1605      // Rounding as in Berthold
1606      while (nCandidate) {
1607        double infeasibility =1.0e-7;
1608        int iRow=-1;
1609        for (i=0;i<numberRows;i++) {
1610          double value=0.0;
1611          if(rowActivity[i]<rowLower[i]) {
1612            value = rowLower[i]-rowActivity[i];
1613          } else if(rowActivity[i]>rowUpper[i]) {
1614            value = rowActivity[i]-rowUpper[i];
1615          }
1616          if (value>infeasibility) {
1617            infeasibility = value;
1618            iRow=i;
1619          }
1620        }
1621        if (iRow>=0) {
1622          // infeasible
1623        } else {
1624          // feasible
1625        }
1626      }
1627    } else {
1628      // Shifting as in Berthold
1629    }
1630    delete [] candidate;
1631  }
1632#endif
1633  delete [] newSolution;
1634  delete [] rowActivity;
1635  return returnCode;
1636}
1637// update model
1638void CbcRounding::setModel(CbcModel * model)
1639{
1640  model_ = model;
1641  // Get a copy of original matrix (and by row for rounding);
1642  assert(model_->solver());
1643  matrix_ = *model_->solver()->getMatrixByCol();
1644  matrixByRow_ = *model_->solver()->getMatrixByRow();
1645  // make sure model okay for heuristic
1646  validate();
1647}
1648// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1649void 
1650CbcRounding::validate() 
1651{
1652  if (model_&&when()<10) {
1653    if (model_->numberIntegers()!=
1654        model_->numberObjects())
1655      setWhen(0);
1656  }
1657#ifdef NEW_ROUNDING
1658  int numberColumns = matrix_.getNumCols();
1659  down_ = new unsigned short [numberColumns];
1660  up_ = new unsigned short [numberColumns];
1661  equal_ = new unsigned short [numberColumns];
1662  // Column copy
1663  const double * element = matrix_.getElements();
1664  const int * row = matrix_.getIndices();
1665  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1666  const int * columnLength = matrix_.getVectorLengths();
1667  const double * rowLower = model.solver()->getRowLower();
1668  const double * rowUpper = model.solver()->getRowUpper();
1669  for (int i=0;i<numberColumns;i++) {
1670    int down=0;
1671    int up=0;
1672    int equal=0;
1673    if (columnLength[i]>65535) {
1674      equal[0]=65535; 
1675      break; // unlikely to work
1676    }
1677    for (CoinBigIndex j=columnStart[i];
1678         j<columnStart[i]+columnLength[i];j++) {
1679      int iRow=row[j];
1680      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
1681        equal++;
1682      } else if (element[j]>0.0) {
1683        if (rowUpper[iRow]<1.0e20)
1684          up++;
1685        else
1686          down--;
1687      } else {
1688        if (rowLower[iRow]>-1.0e20)
1689          up++;
1690        else
1691          down--;
1692      }
1693    }
1694    down_[i] = (unsigned short) down;
1695    up_[i] = (unsigned short) up;
1696    equal_[i] = (unsigned short) equal;
1697  }
1698#else
1699  down_ = NULL;
1700  up_ = NULL;
1701  equal_ = NULL;
1702#endif 
1703}
1704
1705// Default Constructor
1706CbcHeuristicPartial::CbcHeuristicPartial() 
1707  :CbcHeuristic()
1708{
1709  fixPriority_ = 10000;
1710}
1711
1712// Constructor from model
1713CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
1714  :CbcHeuristic(model)
1715{
1716  fixPriority_ = fixPriority;
1717  setNumberNodes(numberNodes);
1718  validate();
1719}
1720
1721// Destructor
1722CbcHeuristicPartial::~CbcHeuristicPartial ()
1723{
1724}
1725
1726// Clone
1727CbcHeuristic *
1728CbcHeuristicPartial::clone() const
1729{
1730  return new CbcHeuristicPartial(*this);
1731}
1732// Create C++ lines to get to current state
1733void 
1734CbcHeuristicPartial::generateCpp( FILE * fp) 
1735{
1736  CbcHeuristicPartial other;
1737  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
1738  fprintf(fp,"3  CbcHeuristicPartial partial(*cbcModel);\n");
1739  CbcHeuristic::generateCpp(fp,"partial");
1740  if (fixPriority_!=other.fixPriority_)
1741    fprintf(fp,"3  partial.setFixPriority(%d);\n",fixPriority_);
1742  else
1743    fprintf(fp,"4  partial.setFixPriority(%d);\n",fixPriority_);
1744  fprintf(fp,"3  cbcModel->addHeuristic(&partial);\n");
1745}
1746//#define NEW_PARTIAL
1747// Copy constructor
1748CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
1749:
1750  CbcHeuristic(rhs),
1751  fixPriority_(rhs.fixPriority_)
1752{
1753}
1754
1755// Assignment operator
1756CbcHeuristicPartial & 
1757CbcHeuristicPartial::operator=( const CbcHeuristicPartial& rhs)
1758{
1759  if (this!=&rhs) {
1760    CbcHeuristic::operator=(rhs);
1761    fixPriority_ = rhs.fixPriority_;
1762  }
1763  return *this;
1764}
1765
1766// Resets stuff if model changes
1767void 
1768CbcHeuristicPartial::resetModel(CbcModel * model)
1769{
1770  model_=model;
1771  // Get a copy of original matrix (and by row for partial);
1772  assert(model_->solver());
1773  validate();
1774}
1775// See if partial will give solution
1776// Sets value of solution
1777// Assumes rhs for original matrix still okay
1778// At present only works with integers
1779// Fix values if asked for
1780// Returns 1 if solution, 0 if not
1781int
1782CbcHeuristicPartial::solution(double & solutionValue,
1783                      double * betterSolution)
1784{
1785  // Return if already done
1786  if (fixPriority_<0)
1787    return 0; // switched off
1788  const double * hotstartSolution = model_->hotstartSolution();
1789  const int * hotstartPriorities = model_->hotstartPriorities();
1790  if (!hotstartSolution)
1791    return 0;
1792  OsiSolverInterface * solver = model_->solver();
1793 
1794  int numberIntegers = model_->numberIntegers();
1795  const int * integerVariable = model_->integerVariable();
1796 
1797  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
1798  const double * colLower = newSolver->getColLower();
1799  const double * colUpper = newSolver->getColUpper();
1800
1801  double primalTolerance;
1802  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
1803   
1804  int i;
1805  int numberFixed=0;
1806  int returnCode=0;
1807
1808  for (i=0;i<numberIntegers;i++) {
1809    int iColumn=integerVariable[i];
1810    if (abs(hotstartPriorities[iColumn])<=fixPriority_) {
1811      double value = hotstartSolution[iColumn];
1812      double lower = colLower[iColumn];
1813      double upper = colUpper[iColumn];
1814      value = CoinMax(value,lower);
1815      value = CoinMin(value,upper);
1816      if (fabs(value-floor(value+0.5))<1.0e-8) {
1817        value = floor(value+0.5);
1818        newSolver->setColLower(iColumn,value);
1819        newSolver->setColUpper(iColumn,value);
1820        numberFixed++;
1821      }
1822    }
1823  }
1824  if (numberFixed>numberIntegers/5-100000000) {
1825#ifdef COIN_DEVELOP
1826    printf("%d integers fixed\n",numberFixed);
1827#endif
1828    returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
1829                                     model_->getCutoff(),"CbcHeuristicPartial");
1830    if (returnCode<0)
1831      returnCode=0; // returned on size
1832    //printf("return code %d",returnCode);
1833    if ((returnCode&2)!=0) {
1834      // could add cut
1835      returnCode &= ~2;
1836      //printf("could add cut with %d elements (if all 0-1)\n",nFix);
1837    } else {
1838      //printf("\n");
1839    }
1840  }
1841  fixPriority_=-1; // switch off
1842 
1843  delete newSolver;
1844  return returnCode;
1845}
1846// update model
1847void CbcHeuristicPartial::setModel(CbcModel * model)
1848{
1849  model_ = model;
1850  assert(model_->solver());
1851  // make sure model okay for heuristic
1852  validate();
1853}
1854// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1855void 
1856CbcHeuristicPartial::validate() 
1857{
1858  if (model_&&when()<10) {
1859    if (model_->numberIntegers()!=
1860        model_->numberObjects())
1861      setWhen(0);
1862  }
1863}
1864
1865// Default Constructor
1866CbcSerendipity::CbcSerendipity() 
1867  :CbcHeuristic()
1868{
1869}
1870
1871// Constructor from model
1872CbcSerendipity::CbcSerendipity(CbcModel & model)
1873  :CbcHeuristic(model)
1874{
1875}
1876
1877// Destructor
1878CbcSerendipity::~CbcSerendipity ()
1879{
1880}
1881
1882// Clone
1883CbcHeuristic *
1884CbcSerendipity::clone() const
1885{
1886  return new CbcSerendipity(*this);
1887}
1888// Create C++ lines to get to current state
1889void 
1890CbcSerendipity::generateCpp( FILE * fp) 
1891{
1892  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
1893  fprintf(fp,"3  CbcSerendipity serendipity(*cbcModel);\n");
1894  CbcHeuristic::generateCpp(fp,"serendipity");
1895  fprintf(fp,"3  cbcModel->addHeuristic(&serendipity);\n");
1896}
1897
1898// Copy constructor
1899CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
1900:
1901  CbcHeuristic(rhs)
1902{
1903}
1904
1905// Assignment operator
1906CbcSerendipity & 
1907CbcSerendipity::operator=( const CbcSerendipity& rhs)
1908{
1909  if (this!=&rhs) {
1910    CbcHeuristic::operator=(rhs);
1911  }
1912  return *this;
1913}
1914
1915// Returns 1 if solution, 0 if not
1916int
1917CbcSerendipity::solution(double & solutionValue,
1918                         double * betterSolution)
1919{
1920  if (!model_)
1921    return 0;
1922  // get information on solver type
1923  OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
1924  OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
1925  if (auxiliaryInfo)
1926    return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
1927  else
1928    return 0;
1929}
1930// update model
1931void CbcSerendipity::setModel(CbcModel * model)
1932{
1933  model_ = model;
1934}
1935// Resets stuff if model changes
1936void 
1937CbcSerendipity::resetModel(CbcModel * model)
1938{
1939  model_ = model;
1940}
1941 
Note: See TracBrowser for help on using the repository browser.