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

Last change on this file since 931 was 931, checked in by forrest, 14 years ago

changes to try and improve performance

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