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

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

for my experiments

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