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

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

changes to try and improve code

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 66.8 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          model_->messageHandler()->message(CBC_RESTART,model_->messages())
660            <<solver2->getNumRows()<<solver2->getNumCols()
661            <<CoinMessageEol;
662          // going for full search and copy across more stuff
663          model.gutsOfCopy(*model_,2);
664          for (int i=0;i<model.numberCutGenerators();i++)
665            model.cutGenerator(i)->setTiming(true);
666          model.setCutoff(cutoff);
667          bool takeHint;
668          OsiHintStrength strength;
669          // Switch off printing if asked to
670          model_->solver()->getHintParam(OsiDoReducePrint,takeHint,strength);
671          model.solver()->setHintParam(OsiDoReducePrint,takeHint,strength);
672          CbcStrategyDefault strategy(true,model_->numberStrong(),
673                                      model_->numberBeforeTrust());
674          // Set up pre-processing - no
675          strategy.setupPreProcessing(0); // was (4);
676          model.setStrategy(strategy);
677        }
678        if (inputSolution_) {
679          // translate and add a serendipity heuristic
680          int numberColumns = solver2->getNumCols();
681          const int * which = process.originalColumns();
682          OsiSolverInterface * solver3 = solver2->clone();
683          for (int i=0;i<numberColumns;i++) {
684            if (solver3->isInteger(i)) {
685              int k=which[i];
686              double value = inputSolution_[k];
687              solver3->setColLower(i,value);
688              solver3->setColUpper(i,value);
689            }
690          }
691          solver3->setDblParam(OsiDualObjectiveLimit,COIN_DBL_MAX);
692          solver3->resolve();
693          if (solver3->isProvenOptimal()) {
694            // good
695            CbcSerendipity heuristic(model);
696            double value = solver3->getObjSense()*solver3->getObjValue();
697            heuristic.setInputSolution(solver3->getColSolution(),value);
698            model.setCutoff(value+1.0e-7*(1.0+fabs(value)));
699            model.addHeuristic(&heuristic,"previous solution");
700          }
701          delete solver3;
702        }
703        // Do search
704        if (logLevel>1)
705          model_->messageHandler()->message(CBC_START_SUB,model_->messages())
706            << name
707            << model.getMaximumNodes()
708            <<CoinMessageEol;
709        // probably faster to use a basis to get integer solutions
710        model.setSpecialOptions(2);
711#ifdef CBC_THREAD
712        if (model_->getNumberThreads()>0&&(model_->getThreadMode()&1)!=0) {
713          // See if at root node
714          bool atRoot = model_->getNodeCount()==0;
715          int passNumber = model_->getCurrentPassNumber();
716          if (atRoot&&passNumber==1)
717            model.setNumberThreads(model_->getNumberThreads());
718        }
719#endif
720        model.setParentModel(*model_);
721        model.setOriginalColumns(process.originalColumns());
722        model.setSearchStrategy(-1);
723        // If no feasibility pump then insert a lightweight one
724        if (feasibilityPumpOptions_>=0) {
725          bool gotPump=false;
726          for (int i=0;i<model.numberHeuristics();i++) {
727            const CbcHeuristicFPump* pump =
728              dynamic_cast<const CbcHeuristicFPump*>(model_->heuristic(i));
729            if (pump) 
730              gotPump=true;
731          }
732          if (!gotPump) {
733            CbcHeuristicFPump heuristic4;
734            int pumpTune=feasibilityPumpOptions_;
735            if (pumpTune>0) {
736              /*
737                >=10000000 for using obj
738                >=1000000 use as accumulate switch
739                >=1000 use index+1 as number of large loops
740                >=100 use 0.05 objvalue as increment
741                >=10 use +0.1 objvalue for cutoff (add)
742                1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
743                4 and static continuous, 5 as 3 but no internal integers
744                6 as 3 but all slack basis!
745              */
746              double value = solver2->getObjSense()*solver2->getObjValue();
747              int w = pumpTune/10;
748              int c = w % 10;
749              w /= 10;
750              int i = w % 10;
751              w /= 10;
752              int r = w;
753              int accumulate = r/1000;
754              r -= 1000*accumulate;
755              if (accumulate>=10) {
756                int which = accumulate/10;
757                accumulate -= 10*which;
758                which--;
759                // weights and factors
760                double weight[]={0.1,0.1,0.5,0.5,1.0,1.0,5.0,5.0};
761                double factor[] = {0.1,0.5,0.1,0.5,0.1,0.5,0.1,0.5};
762                heuristic4.setInitialWeight(weight[which]);
763                heuristic4.setWeightFactor(factor[which]);
764              }
765              // fake cutoff
766              if (c) {
767                double cutoff;
768                solver2->getDblParam(OsiDualObjectiveLimit,cutoff);
769                cutoff = CoinMin(cutoff,value + 0.1*fabs(value)*c);
770                heuristic4.setFakeCutoff(cutoff);
771              }
772              if (i||r) {
773                // also set increment
774                //double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
775                double increment = 0.0;
776                heuristic4.setAbsoluteIncrement(increment);
777                heuristic4.setAccumulate(accumulate);
778                heuristic4.setMaximumRetries(r+1);
779              }
780              pumpTune = pumpTune%100;
781              if (pumpTune==6)
782                pumpTune =13;
783              heuristic4.setWhen(pumpTune+10);
784            }
785            model.addHeuristic(&heuristic4,"feasibility pump",0);
786          }
787        }
788        if (model_->searchStrategy()==2) {
789          model.setNumberStrong(5);
790          model.setNumberBeforeTrust(5);
791        }
792        if (model.getNumCols()) {
793          if (numberNodes>=0)
794            setCutAndHeuristicOptions(model);
795          model.branchAndBound();
796          if (numberNodes<0) {
797            model_->incrementIterationCount(model.getIterationCount());
798            model_->incrementNodeCount(model.getNodeCount());
799            for (int iGenerator=0;iGenerator<model.numberCutGenerators();iGenerator++) {
800              CbcCutGenerator * generator = model.cutGenerator(iGenerator);
801              printf("%s was tried %d times and created %d cuts of which %d were active after adding rounds of cuts",
802                      generator->cutGeneratorName(),
803                      generator->numberTimesEntered(),
804                      generator->numberCutsInTotal()+
805                      generator->numberColumnCuts(),
806                      generator->numberCutsActive());
807              printf(" (%.3f seconds)\n)",generator->timeInCutGenerator());
808            }
809          }
810        } else {
811          // empty model
812          model.setMinimizationObjValue(model.solver()->getObjSense()*model.solver()->getObjValue());
813        }
814        if (logLevel>1)
815          model_->messageHandler()->message(CBC_END_SUB,model_->messages())
816            << name
817            <<CoinMessageEol;
818        if (model.getMinimizationObjValue()<CoinMin(cutoff,1.0e30)) {
819          // solution
820          if (model.getNumCols())
821            returnCode=model.isProvenOptimal() ? 3 : 1;
822          else
823            returnCode=3;
824          // post process
825#ifdef COIN_HAS_CLP
826          OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model.solver());
827          if (clpSolver) {
828            ClpSimplex * lpSolver = clpSolver->getModelPtr();
829            lpSolver->setSpecialOptions(lpSolver->specialOptions()|0x01000000); // say is Cbc (and in branch and bound)
830          }
831#endif
832          process.postProcess(*model.solver());
833          if (solver->isProvenOptimal()) {
834            // Solution now back in solver
835            int numberColumns = solver->getNumCols();
836            memcpy(newSolution,solver->getColSolution(),
837                   numberColumns*sizeof(double));
838            newSolutionValue = model.getMinimizationObjValue();
839          } else {
840            // odd - but no good
841            returnCode=0; // so will be infeasible
842          }
843        } else {
844        // no good
845          returnCode=model.isProvenInfeasible() ? 2 : 0; // so will be infeasible
846        }
847        if (model.status()==5)
848          returnCode=-2; // stop
849        if (model.isProvenInfeasible())
850          status=1;
851        else if (model.isProvenOptimal())
852          status=2;
853      }
854    }
855  } else {
856    returnCode=2; // infeasible finished
857  }
858  model_->setLogLevel(logLevel);
859  if (reset) {
860    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
861      if (reset[iColumn])
862        solver->setColLower(iColumn,0.0);
863    }
864    delete [] reset;
865  }
866#ifdef COIN_DEVELOP
867  getHistoryStatistics_=true;
868  if (returnCode==1||returnCode==2) {
869    if (status==1)
870      printf("heuristic could add cut because infeasible (%s)\n",heuristicName_.c_str()); 
871    else if (status==2)
872      printf("heuristic could add cut because optimal (%s)\n",heuristicName_.c_str());
873  } 
874#endif
875  solver->setHintParam(OsiDoReducePrint,takeHint,strength);
876  return returnCode;
877}
878// Set input solution
879void 
880CbcHeuristic::setInputSolution(const double * solution, double objValue)
881{
882  delete [] inputSolution_;
883  inputSolution_=NULL;
884  if (model_&&solution) {
885    int numberColumns = model_->getNumCols();
886    inputSolution_ = new double [numberColumns+1];
887    memcpy(inputSolution_,solution,numberColumns*sizeof(double));
888    inputSolution_[numberColumns]=objValue;
889  }
890}
891
892//##############################################################################
893
894inline int compare3BranchingObjects(const CbcBranchingObject* br0,
895                                    const CbcBranchingObject* br1)
896{
897  const int t0 = br0->type();
898  const int t1 = br1->type();
899  if (t0 < t1) {
900    return -1;
901  }
902  if (t0 > t1) {
903    return 1;
904  }
905  return br0->compareOriginalObject(br1);
906}
907
908//==============================================================================
909
910inline bool compareBranchingObjects(const CbcBranchingObject* br0,
911                                    const CbcBranchingObject* br1)
912{
913  return compare3BranchingObjects(br0, br1) < 0;
914}
915
916//==============================================================================
917
918void
919CbcHeuristicNode::gutsOfConstructor(CbcModel& model)
920{
921  //  CbcHeurDebugNodes(&model);
922  CbcNode* node = model.currentNode();
923  brObj_ = new CbcBranchingObject*[node->depth()];
924  CbcNodeInfo* nodeInfo = node->nodeInfo();
925  int cnt = 0;
926  while (nodeInfo->parentBranch() != NULL) {
927    const OsiBranchingObject* br = nodeInfo->parentBranch();
928    const CbcBranchingObject* cbcbr = dynamic_cast<const CbcBranchingObject*>(br);
929    if (! cbcbr) {
930      throw CoinError("CbcHeuristicNode can be used only with CbcBranchingObjects.\n",
931                      "gutsOfConstructor",
932                      "CbcHeuristicNode",
933                      __FILE__, __LINE__);
934    }
935    brObj_[cnt] = cbcbr->clone();
936    brObj_[cnt]->previousBranch();
937    ++cnt;
938    nodeInfo = nodeInfo->parent();
939  }
940  std::sort(brObj_, brObj_+cnt, compareBranchingObjects);
941  if (cnt <= 1) {
942    numObjects_ = cnt;
943  } else {
944    numObjects_ = 0;
945    CbcBranchingObject* br=NULL; // What should this be?
946    for (int i = 1; i < cnt; ++i) {
947      if (compare3BranchingObjects(brObj_[numObjects_], brObj_[i]) == 0) {
948        int comp = brObj_[numObjects_]->compareBranchingObject(brObj_[i], &br);
949        switch (comp) {
950        case CbcRangeSame: // the same range
951        case CbcRangeDisjoint: // disjoint decisions
952          // should not happen! we are on a chain!
953          abort();
954        case CbcRangeSubset: // brObj_[numObjects_] is a subset of brObj_[i]
955          delete brObj_[i];
956          break;
957        case CbcRangeSuperset: // brObj_[i] is a subset of brObj_[numObjects_]
958          delete brObj_[numObjects_];
959          brObj_[numObjects_] = brObj_[i];
960          break;
961        case CbcRangeOverlap: // overlap
962          delete brObj_[i];
963          delete brObj_[numObjects_];
964          brObj_[numObjects_] = br;
965          break;
966      }
967        continue;
968      } else {
969        brObj_[++numObjects_] = brObj_[i];
970      }
971    }
972    ++numObjects_;
973  }
974}
975
976//==============================================================================
977
978CbcHeuristicNode::CbcHeuristicNode(CbcModel& model)
979{
980  gutsOfConstructor(model);
981}
982
983//==============================================================================
984
985double
986CbcHeuristicNode::distance(const CbcHeuristicNode* node) const 
987{
988
989  const double disjointWeight = 1;
990  const double overlapWeight = 0.4;
991  const double subsetWeight = 0.2;
992  int countDisjointWeight = 0;
993  int countOverlapWeight = 0;
994  int countSubsetWeight = 0;
995  int i = 0; 
996  int j = 0;
997  double dist = 0.0;
998#ifdef PRINT_DEBUG
999  printf(" numObjects_ = %i, node->numObjects_ = %i\n",
1000         numObjects_, node->numObjects_);
1001#endif
1002  while( i < numObjects_ && j < node->numObjects_) {
1003    CbcBranchingObject* br0 = brObj_[i];
1004    const CbcBranchingObject* br1 = node->brObj_[j];
1005#ifdef PRINT_DEBUG
1006    const CbcIntegerBranchingObject* brPrint0 =
1007      dynamic_cast<const CbcIntegerBranchingObject*>(br0);
1008    const double* downBounds = brPrint0->downBounds();
1009    const double* upBounds = brPrint0->upBounds();
1010    int variable = brPrint0->variable();
1011    int way = brPrint0->way();
1012    printf("   br0: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1013           variable, (int)downBounds[0], (int)downBounds[1],
1014           (int)upBounds[0], (int)upBounds[1], way);
1015    const CbcIntegerBranchingObject* brPrint1 =
1016      dynamic_cast<const CbcIntegerBranchingObject*>(br1);
1017    downBounds = brPrint1->downBounds();
1018    upBounds = brPrint1->upBounds();
1019    variable = brPrint1->variable();
1020    way = brPrint1->way();
1021    printf("   br1: var %i downBd [%i,%i] upBd [%i,%i] way %i\n",
1022           variable, (int)downBounds[0], (int)downBounds[1],
1023           (int)upBounds[0], (int)upBounds[1], way);
1024#endif
1025    const int brComp = compare3BranchingObjects(br0, br1);
1026    if (brComp < 0) {
1027      dist += subsetWeight;
1028      countSubsetWeight++;
1029      ++i;
1030    }
1031    else if (brComp > 0) {
1032      dist += subsetWeight;
1033      countSubsetWeight++;
1034      ++j;
1035    }
1036    else {
1037      const int comp = br0->compareBranchingObject(br1, false);
1038      switch (comp) {
1039      case CbcRangeSame:
1040        // do nothing
1041        break;
1042      case CbcRangeDisjoint: // disjoint decisions
1043        dist += disjointWeight;
1044        countDisjointWeight++;
1045        break;
1046      case CbcRangeSubset: // subset one way or another
1047      case CbcRangeSuperset:
1048        dist += subsetWeight;
1049        countSubsetWeight++;
1050        break;
1051      case CbcRangeOverlap: // overlap
1052        dist += overlapWeight;
1053        countOverlapWeight++;
1054        break;
1055      }
1056      ++i;
1057      ++j;
1058    }
1059  }
1060  dist += subsetWeight * (numObjects_ - i + node->numObjects_ - j);
1061  countSubsetWeight += (numObjects_ - i + node->numObjects_ - j);
1062  printf("subset = %i, overlap = %i, disjoint = %i\n", countSubsetWeight,
1063         countOverlapWeight, countDisjointWeight);
1064  return dist;
1065}
1066
1067//==============================================================================
1068
1069CbcHeuristicNode::~CbcHeuristicNode()
1070{
1071  for (int i = 0; i < numObjects_; ++i) {
1072    delete brObj_[i];
1073  }
1074  delete [] brObj_;
1075}
1076
1077//==============================================================================
1078
1079double
1080CbcHeuristicNode::minDistance(const CbcHeuristicNodeList& nodeList) const
1081{
1082  double minDist = COIN_DBL_MAX;
1083  for (int i = nodeList.size() - 1; i >= 0; --i) {
1084    minDist = CoinMin(minDist, distance(nodeList.node(i)));
1085  }
1086  return minDist;
1087}
1088
1089//==============================================================================
1090
1091bool
1092CbcHeuristicNode::minDistanceIsSmall(const CbcHeuristicNodeList& nodeList,
1093                                     const double threshold) const
1094{
1095  for (int i = nodeList.size() - 1; i >= 0; --i) {
1096    if (distance(nodeList.node(i)) >= threshold) {
1097      continue;
1098    } else {
1099      return true;
1100    }
1101  }
1102  return false;
1103}
1104
1105//==============================================================================
1106
1107double
1108CbcHeuristicNode::avgDistance(const CbcHeuristicNodeList& nodeList) const
1109{
1110  if (nodeList.size() == 0) {
1111    return COIN_DBL_MAX;
1112  }
1113  double sumDist = 0;
1114  for (int i = nodeList.size() - 1; i >= 0; --i) {
1115    sumDist += distance(nodeList.node(i));
1116  }
1117  return sumDist/nodeList.size();
1118}
1119
1120//##############################################################################
1121
1122// Default Constructor
1123CbcRounding::CbcRounding() 
1124  :CbcHeuristic()
1125{
1126  // matrix and row copy will automatically be empty
1127  seed_=1;
1128  down_ = NULL;
1129  up_ = NULL;
1130  equal_ = NULL;
1131}
1132
1133// Constructor from model
1134CbcRounding::CbcRounding(CbcModel & model)
1135  :CbcHeuristic(model)
1136{
1137  // Get a copy of original matrix (and by row for rounding);
1138  assert(model.solver());
1139  matrix_ = *model.solver()->getMatrixByCol();
1140  matrixByRow_ = *model.solver()->getMatrixByRow();
1141  validate();
1142  seed_=1;
1143}
1144
1145// Destructor
1146CbcRounding::~CbcRounding ()
1147{
1148  delete [] down_;
1149  delete [] up_;
1150  delete [] equal_;
1151}
1152
1153// Clone
1154CbcHeuristic *
1155CbcRounding::clone() const
1156{
1157  return new CbcRounding(*this);
1158}
1159// Create C++ lines to get to current state
1160void 
1161CbcRounding::generateCpp( FILE * fp) 
1162{
1163  CbcRounding other;
1164  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
1165  fprintf(fp,"3  CbcRounding rounding(*cbcModel);\n");
1166  CbcHeuristic::generateCpp(fp,"rounding");
1167  if (seed_!=other.seed_)
1168    fprintf(fp,"3  rounding.setSeed(%d);\n",seed_);
1169  else
1170    fprintf(fp,"4  rounding.setSeed(%d);\n",seed_);
1171  fprintf(fp,"3  cbcModel->addHeuristic(&rounding);\n");
1172}
1173//#define NEW_ROUNDING
1174// Copy constructor
1175CbcRounding::CbcRounding(const CbcRounding & rhs)
1176:
1177  CbcHeuristic(rhs),
1178  matrix_(rhs.matrix_),
1179  matrixByRow_(rhs.matrixByRow_),
1180  seed_(rhs.seed_)
1181{
1182#ifdef NEW_ROUNDING
1183  int numberColumns = matrix_.getNumCols();
1184  down_ = CoinCopyOfArray(rhs.down_,numberColumns);
1185  up_ = CoinCopyOfArray(rhs.up_,numberColumns);
1186  equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
1187#else
1188  down_ = NULL;
1189  up_ = NULL;
1190  equal_ = NULL;
1191#endif 
1192}
1193
1194// Assignment operator
1195CbcRounding & 
1196CbcRounding::operator=( const CbcRounding& rhs)
1197{
1198  if (this!=&rhs) {
1199    CbcHeuristic::operator=(rhs);
1200    matrix_ = rhs.matrix_;
1201    matrixByRow_ = rhs.matrixByRow_;
1202#ifdef NEW_ROUNDING
1203    delete [] down_;
1204    delete [] up_;
1205    delete [] equal_;
1206    int numberColumns = matrix_.getNumCols();
1207    down_ = CoinCopyOfArray(rhs.down_,numberColumns);
1208    up_ = CoinCopyOfArray(rhs.up_,numberColumns);
1209    equal_ = CoinCopyOfArray(rhs.equal_,numberColumns);
1210#else
1211    down_ = NULL;
1212    up_ = NULL;
1213    equal_ = NULL;
1214#endif 
1215    seed_ = rhs.seed_;
1216  }
1217  return *this;
1218}
1219
1220// Resets stuff if model changes
1221void 
1222CbcRounding::resetModel(CbcModel * model)
1223{
1224  model_=model;
1225  // Get a copy of original matrix (and by row for rounding);
1226  assert(model_->solver());
1227  matrix_ = *model_->solver()->getMatrixByCol();
1228  matrixByRow_ = *model_->solver()->getMatrixByRow();
1229  validate();
1230}
1231// See if rounding will give solution
1232// Sets value of solution
1233// Assumes rhs for original matrix still okay
1234// At present only works with integers
1235// Fix values if asked for
1236// Returns 1 if solution, 0 if not
1237int
1238CbcRounding::solution(double & solutionValue,
1239                      double * betterSolution)
1240{
1241
1242  // See if to do
1243  if (!when()||(when()%10==1&&model_->phase()!=1)||
1244      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
1245    return 0; // switched off
1246  OsiSolverInterface * solver = model_->solver();
1247  double direction = solver->getObjSense();
1248  double newSolutionValue = direction*solver->getObjValue();
1249  return solution(solutionValue,betterSolution,newSolutionValue);
1250}
1251// See if rounding will give solution
1252// Sets value of solution
1253// Assumes rhs for original matrix still okay
1254// At present only works with integers
1255// Fix values if asked for
1256// Returns 1 if solution, 0 if not
1257int
1258CbcRounding::solution(double & solutionValue,
1259                      double * betterSolution,
1260                      double newSolutionValue)
1261{
1262
1263  // See if to do
1264  if (!when()||(when()%10==1&&model_->phase()!=1)||
1265      (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
1266    return 0; // switched off
1267  OsiSolverInterface * solver = model_->solver();
1268  const double * lower = solver->getColLower();
1269  const double * upper = solver->getColUpper();
1270  const double * rowLower = solver->getRowLower();
1271  const double * rowUpper = solver->getRowUpper();
1272  const double * solution = solver->getColSolution();
1273  const double * objective = solver->getObjCoefficients();
1274  double integerTolerance = model_->getDblParam(CbcModel::CbcIntegerTolerance);
1275  double primalTolerance;
1276  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
1277
1278  int numberRows = matrix_.getNumRows();
1279  assert (numberRows<=solver->getNumRows());
1280  int numberIntegers = model_->numberIntegers();
1281  const int * integerVariable = model_->integerVariable();
1282  int i;
1283  double direction = solver->getObjSense();
1284  //double newSolutionValue = direction*solver->getObjValue();
1285  int returnCode = 0;
1286  // Column copy
1287  const double * element = matrix_.getElements();
1288  const int * row = matrix_.getIndices();
1289  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
1290  const int * columnLength = matrix_.getVectorLengths();
1291  // Row copy
1292  const double * elementByRow = matrixByRow_.getElements();
1293  const int * column = matrixByRow_.getIndices();
1294  const CoinBigIndex * rowStart = matrixByRow_.getVectorStarts();
1295  const int * rowLength = matrixByRow_.getVectorLengths();
1296
1297  // Get solution array for heuristic solution
1298  int numberColumns = solver->getNumCols();
1299  double * newSolution = new double [numberColumns];
1300  memcpy(newSolution,solution,numberColumns*sizeof(double));
1301
1302  double * rowActivity = new double[numberRows];
1303  memset(rowActivity,0,numberRows*sizeof(double));
1304  for (i=0;i<numberColumns;i++) {
1305    int j;
1306    double value = newSolution[i];
1307    if (value<lower[i]) {
1308      value=lower[i];
1309      newSolution[i]=value;
1310    } else if (value>upper[i]) {
1311      value=upper[i];
1312      newSolution[i]=value;
1313    }
1314    if (value) {
1315      for (j=columnStart[i];
1316           j<columnStart[i]+columnLength[i];j++) {
1317        int iRow=row[j];
1318        rowActivity[iRow] += value*element[j];
1319      }
1320    }
1321  }
1322  // check was feasible - if not adjust (cleaning may move)
1323  for (i=0;i<numberRows;i++) {
1324    if(rowActivity[i]<rowLower[i]) {
1325      //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1326      rowActivity[i]=rowLower[i];
1327    } else if(rowActivity[i]>rowUpper[i]) {
1328      //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1329      rowActivity[i]=rowUpper[i];
1330    }
1331  }
1332  for (i=0;i<numberIntegers;i++) {
1333    int iColumn = integerVariable[i];
1334    double value=newSolution[iColumn];
1335    if (fabs(floor(value+0.5)-value)>integerTolerance) {
1336      double below = floor(value);
1337      double newValue=newSolution[iColumn];
1338      double cost = direction * objective[iColumn];
1339      double move;
1340      if (cost>0.0) {
1341        // try up
1342        move = 1.0 -(value-below);
1343      } else if (cost<0.0) {
1344        // try down
1345        move = below-value;
1346      } else {
1347        // won't be able to move unless we can grab another variable
1348        double randomNumber = randomNumberGenerator_.randomDouble();
1349        // which way?
1350        if (randomNumber<0.5) 
1351          move = below-value;
1352        else
1353          move = 1.0 -(value-below);
1354      }
1355      newValue += move;
1356      newSolution[iColumn] = newValue;
1357      newSolutionValue += move*cost;
1358      int j;
1359      for (j=columnStart[iColumn];
1360           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1361        int iRow = row[j];
1362        rowActivity[iRow] += move*element[j];
1363      }
1364    }
1365  }
1366
1367  double penalty=0.0;
1368  const char * integerType = model_->integerType();
1369  // see if feasible - just using singletons
1370  for (i=0;i<numberRows;i++) {
1371    double value = rowActivity[i];
1372    double thisInfeasibility=0.0;
1373    if (value<rowLower[i]-primalTolerance)
1374      thisInfeasibility = value-rowLower[i];
1375    else if (value>rowUpper[i]+primalTolerance)
1376      thisInfeasibility = value-rowUpper[i];
1377    if (thisInfeasibility) {
1378      // See if there are any slacks I can use to fix up
1379      // maybe put in coding for multiple slacks?
1380      double bestCost = 1.0e50;
1381      int k;
1382      int iBest=-1;
1383      double addCost=0.0;
1384      double newValue=0.0;
1385      double changeRowActivity=0.0;
1386      double absInfeasibility = fabs(thisInfeasibility);
1387      for (k=rowStart[i];k<rowStart[i]+rowLength[i];k++) {
1388        int iColumn = column[k];
1389        // See if all elements help
1390        if (columnLength[iColumn]==1) {
1391          double currentValue = newSolution[iColumn];
1392          double elementValue = elementByRow[k];
1393          double lowerValue = lower[iColumn];
1394          double upperValue = upper[iColumn];
1395          double gap = rowUpper[i]-rowLower[i];
1396          double absElement=fabs(elementValue);
1397          if (thisInfeasibility*elementValue>0.0) {
1398            // we want to reduce
1399            if ((currentValue-lowerValue)*absElement>=absInfeasibility) {
1400              // possible - check if integer
1401              double distance = absInfeasibility/absElement;
1402              double thisCost = -direction*objective[iColumn]*distance;
1403              if (integerType[iColumn]) {
1404                distance = ceil(distance-primalTolerance);
1405                if (currentValue-distance>=lowerValue-primalTolerance) {
1406                  if (absInfeasibility-distance*absElement< -gap-primalTolerance)
1407                    thisCost=1.0e100; // no good
1408                  else
1409                    thisCost = -direction*objective[iColumn]*distance;
1410                } else {
1411                  thisCost=1.0e100; // no good
1412                }
1413              }
1414              if (thisCost<bestCost) {
1415                bestCost=thisCost;
1416                iBest=iColumn;
1417                addCost = thisCost;
1418                newValue = currentValue-distance;
1419                changeRowActivity = -distance*elementValue;
1420              }
1421            }
1422          } else {
1423            // we want to increase
1424            if ((upperValue-currentValue)*absElement>=absInfeasibility) {
1425              // possible - check if integer
1426              double distance = absInfeasibility/absElement;
1427              double thisCost = direction*objective[iColumn]*distance;
1428              if (integerType[iColumn]) {
1429                distance = ceil(distance-1.0e-7);
1430                assert (currentValue-distance<=upperValue+primalTolerance);
1431                if (absInfeasibility-distance*absElement< -gap-primalTolerance)
1432                  thisCost=1.0e100; // no good
1433                else
1434                  thisCost = direction*objective[iColumn]*distance;
1435              }
1436              if (thisCost<bestCost) {
1437                bestCost=thisCost;
1438                iBest=iColumn;
1439                addCost = thisCost;
1440                newValue = currentValue+distance;
1441                changeRowActivity = distance*elementValue;
1442              }
1443            }
1444          }
1445        }
1446      }
1447      if (iBest>=0) {
1448        /*printf("Infeasibility of %g on row %d cost %g\n",
1449          thisInfeasibility,i,addCost);*/
1450        newSolution[iBest]=newValue;
1451        thisInfeasibility=0.0;
1452        newSolutionValue += addCost;
1453        rowActivity[i] += changeRowActivity;
1454      }
1455      penalty += fabs(thisInfeasibility);
1456    }
1457  }
1458  if (penalty) {
1459    // see if feasible using any
1460    // first continuous
1461    double penaltyChange=0.0;
1462    int iColumn;
1463    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1464      if (integerType[iColumn])
1465        continue;
1466      double currentValue = newSolution[iColumn];
1467      double lowerValue = lower[iColumn];
1468      double upperValue = upper[iColumn];
1469      int j;
1470      int anyBadDown=0;
1471      int anyBadUp=0;
1472      double upImprovement=0.0;
1473      double downImprovement=0.0;
1474      for (j=columnStart[iColumn];
1475           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1476        int iRow = row[j];
1477        if (rowUpper[iRow]>rowLower[iRow]) {
1478          double value = element[j];
1479          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1480            // infeasible above
1481            downImprovement += value;
1482            upImprovement -= value;
1483            if (value>0.0) 
1484              anyBadUp++;
1485            else 
1486              anyBadDown++;
1487          } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
1488            // feasible at ub
1489            if (value>0.0) {
1490              upImprovement -= value;
1491              anyBadUp++;
1492            } else {
1493              downImprovement += value;
1494              anyBadDown++;
1495            }
1496          } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
1497            // feasible in interior
1498          } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
1499            // feasible at lb
1500            if (value<0.0) {
1501              upImprovement += value;
1502              anyBadUp++;
1503            } else {
1504              downImprovement -= value;
1505              anyBadDown++;
1506            }
1507          } else {
1508            // infeasible below
1509            downImprovement -= value;
1510            upImprovement += value;
1511            if (value<0.0) 
1512              anyBadUp++;
1513            else 
1514              anyBadDown++;
1515          }
1516        } else {
1517          // equality row
1518          double value = element[j];
1519          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1520            // infeasible above
1521            downImprovement += value;
1522            upImprovement -= value;
1523            if (value>0.0) 
1524              anyBadUp++;
1525            else 
1526              anyBadDown++;
1527          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1528            // infeasible below
1529            downImprovement -= value;
1530            upImprovement += value;
1531            if (value<0.0) 
1532              anyBadUp++;
1533            else 
1534              anyBadDown++;
1535          } else {
1536            // feasible - no good
1537            anyBadUp=-1;
1538            anyBadDown=-1;
1539            break;
1540          }
1541        }
1542      }
1543      // could change tests for anyBad
1544      if (anyBadUp)
1545        upImprovement=0.0;
1546      if (anyBadDown)
1547        downImprovement=0.0;
1548      double way=0.0;
1549      double improvement=0.0;
1550      if (downImprovement>0.0&&currentValue>lowerValue) {
1551        way=-1.0;
1552        improvement = downImprovement;
1553      } else if (upImprovement>0.0&&currentValue<upperValue) {
1554        way=1.0;
1555        improvement = upImprovement;
1556      }
1557      if (way) {
1558        // can improve
1559        double distance;
1560        if (way>0.0)
1561          distance = upperValue-currentValue;
1562        else
1563          distance = currentValue-lowerValue;
1564        for (j=columnStart[iColumn];
1565             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1566          int iRow = row[j];
1567          double value = element[j]*way;
1568          if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1569            // infeasible above
1570            assert (value<0.0);
1571            double gap = rowActivity[iRow]-rowUpper[iRow];
1572            if (gap+value*distance<0.0) 
1573              distance = -gap/value;
1574          } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1575            // infeasible below
1576            assert (value>0.0);
1577            double gap = rowActivity[iRow]-rowLower[iRow];
1578            if (gap+value*distance>0.0) 
1579              distance = -gap/value;
1580          } else {
1581            // feasible
1582            if (value>0) {
1583              double gap = rowActivity[iRow]-rowUpper[iRow];
1584              if (gap+value*distance>0.0) 
1585              distance = -gap/value;
1586            } else {
1587              double gap = rowActivity[iRow]-rowLower[iRow];
1588              if (gap+value*distance<0.0) 
1589                distance = -gap/value;
1590            }
1591          }
1592        }
1593        //move
1594        penaltyChange += improvement*distance;
1595        distance *= way;
1596        newSolution[iColumn] += distance;
1597        newSolutionValue += direction*objective[iColumn]*distance;
1598        for (j=columnStart[iColumn];
1599             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1600          int iRow = row[j];
1601          double value = element[j];
1602          rowActivity[iRow] += distance*value;
1603        }
1604      }
1605    }
1606    // and now all if improving
1607    double lastChange= penaltyChange ? 1.0 : 0.0;
1608    while (lastChange>1.0e-2) {
1609      lastChange=0;
1610      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1611        bool isInteger = (integerType[iColumn]!=0);
1612        double currentValue = newSolution[iColumn];
1613        double lowerValue = lower[iColumn];
1614        double upperValue = upper[iColumn];
1615        int j;
1616        int anyBadDown=0;
1617        int anyBadUp=0;
1618        double upImprovement=0.0;
1619        double downImprovement=0.0;
1620        for (j=columnStart[iColumn];
1621             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1622          int iRow = row[j];
1623          double value = element[j];
1624          if (isInteger) {
1625            if (value>0.0) {
1626              if (rowActivity[iRow]+value>rowUpper[iRow]+primalTolerance)
1627                anyBadUp++;
1628              if (rowActivity[iRow]-value<rowLower[iRow]-primalTolerance)
1629                anyBadDown++;
1630            } else {
1631              if (rowActivity[iRow]-value>rowUpper[iRow]+primalTolerance)
1632                anyBadDown++;
1633              if (rowActivity[iRow]+value<rowLower[iRow]-primalTolerance)
1634                anyBadUp++;
1635            }
1636          }
1637          if (rowUpper[iRow]>rowLower[iRow]) {
1638            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1639              // infeasible above
1640              downImprovement += value;
1641              upImprovement -= value;
1642              if (value>0.0) 
1643                anyBadUp++;
1644              else 
1645                anyBadDown++;
1646            } else if (rowActivity[iRow]>rowUpper[iRow]-primalTolerance) {
1647              // feasible at ub
1648              if (value>0.0) {
1649                upImprovement -= value;
1650                anyBadUp++;
1651              } else {
1652                downImprovement += value;
1653                anyBadDown++;
1654              }
1655            } else if (rowActivity[iRow]>rowLower[iRow]+primalTolerance) {
1656              // feasible in interior
1657            } else if (rowActivity[iRow]>rowLower[iRow]-primalTolerance) {
1658              // feasible at lb
1659              if (value<0.0) {
1660                upImprovement += value;
1661                anyBadUp++;
1662              } else {
1663                downImprovement -= value;
1664                anyBadDown++;
1665              }
1666            } else {
1667              // infeasible below
1668              downImprovement -= value;
1669              upImprovement += value;
1670              if (value<0.0) 
1671                anyBadUp++;
1672              else 
1673                anyBadDown++;
1674            }
1675          } else {
1676            // equality row
1677            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1678              // infeasible above
1679              downImprovement += value;
1680              upImprovement -= value;
1681              if (value>0.0) 
1682                anyBadUp++;
1683              else 
1684                anyBadDown++;
1685            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1686              // infeasible below
1687              downImprovement -= value;
1688              upImprovement += value;
1689              if (value<0.0) 
1690                anyBadUp++;
1691              else 
1692                anyBadDown++;
1693            } else {
1694              // feasible - no good
1695              anyBadUp=-1;
1696              anyBadDown=-1;
1697              break;
1698            }
1699          }
1700        }
1701        // could change tests for anyBad
1702        if (anyBadUp)
1703          upImprovement=0.0;
1704        if (anyBadDown)
1705          downImprovement=0.0;
1706        double way=0.0;
1707        double improvement=0.0;
1708        if (downImprovement>0.0&&currentValue>lowerValue) {
1709          way=-1.0;
1710          improvement = downImprovement;
1711        } else if (upImprovement>0.0&&currentValue<upperValue) {
1712          way=1.0;
1713          improvement = upImprovement;
1714        }
1715        if (way) {
1716          // can improve
1717          double distance=COIN_DBL_MAX;
1718          for (j=columnStart[iColumn];
1719               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1720            int iRow = row[j];
1721            double value = element[j]*way;
1722            if (rowActivity[iRow]>rowUpper[iRow]+primalTolerance) {
1723              // infeasible above
1724              assert (value<0.0);
1725              double gap = rowActivity[iRow]-rowUpper[iRow];
1726              if (gap+value*distance<0.0) {
1727                // If integer then has to move by 1
1728                if (!isInteger)
1729                  distance = -gap/value;
1730                else
1731                  distance = CoinMax(-gap/value,1.0);
1732              }
1733            } else if (rowActivity[iRow]<rowLower[iRow]-primalTolerance) {
1734              // infeasible below
1735              assert (value>0.0);
1736              double gap = rowActivity[iRow]-rowLower[iRow];
1737              if (gap+value*distance>0.0) {
1738                // If integer then has to move by 1
1739                if (!isInteger)
1740                  distance = -gap/value;
1741                else
1742                  distance = CoinMax(-gap/value,1.0);
1743              }
1744            } else {
1745              // feasible
1746              if (value>0) {
1747                double gap = rowActivity[iRow]-rowUpper[iRow];
1748                if (gap+value*distance>0.0) 
1749                  distance = -gap/value;
1750              } else {
1751                double gap = rowActivity[iRow]-rowLower[iRow];
1752                if (gap+value*distance<0.0) 
1753                  distance = -gap/value;
1754              }
1755            }
1756          }
1757          if (isInteger)
1758            distance = floor(distance+1.05e-8);
1759          if (!distance) {
1760            // should never happen
1761            //printf("zero distance in CbcRounding - debug\n");
1762          }
1763          //move
1764          lastChange += improvement*distance;
1765          distance *= way;
1766          newSolution[iColumn] += distance;
1767          newSolutionValue += direction*objective[iColumn]*distance;
1768          for (j=columnStart[iColumn];
1769               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1770            int iRow = row[j];
1771            double value = element[j];
1772            rowActivity[iRow] += distance*value;
1773          }
1774        }
1775      }
1776      penaltyChange += lastChange;
1777    }
1778    penalty -= penaltyChange;
1779    if (penalty<1.0e-5*fabs(penaltyChange)) {
1780      // recompute
1781      penalty=0.0;
1782      for (i=0;i<numberRows;i++) {
1783        double value = rowActivity[i];
1784        if (value<rowLower[i]-primalTolerance)
1785          penalty += rowLower[i]-value;
1786        else if (value>rowUpper[i]+primalTolerance)
1787          penalty += value-rowUpper[i];
1788      }
1789    }
1790  }
1791
1792  // Could also set SOS (using random) and repeat
1793  if (!penalty) {
1794    // See if we can do better
1795    //seed_++;
1796    //CoinSeedRandom(seed_);
1797    // Random number between 0 and 1.
1798    double randomNumber = randomNumberGenerator_.randomDouble();
1799    int iPass;
1800    int start[2];
1801    int end[2];
1802    int iRandom = (int) (randomNumber*((double) numberIntegers));
1803    start[0]=iRandom;
1804    end[0]=numberIntegers;
1805    start[1]=0;
1806    end[1]=iRandom;
1807    for (iPass=0;iPass<2;iPass++) {
1808      int i;
1809      for (i=start[iPass];i<end[iPass];i++) {
1810        int iColumn = integerVariable[i];
1811#ifndef NDEBUG
1812        double value=newSolution[iColumn];
1813        assert (fabs(floor(value+0.5)-value)<integerTolerance);
1814#endif
1815        double cost = direction * objective[iColumn];
1816        double move=0.0;
1817        if (cost>0.0)
1818          move = -1.0;
1819        else if (cost<0.0)
1820          move=1.0;
1821        while (move) {
1822          bool good=true;
1823          double newValue=newSolution[iColumn]+move;
1824          if (newValue<lower[iColumn]-primalTolerance||
1825              newValue>upper[iColumn]+primalTolerance) {
1826            move=0.0;
1827          } else {
1828            // see if we can move
1829            int j;
1830            for (j=columnStart[iColumn];
1831                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
1832              int iRow = row[j];
1833              double newActivity = rowActivity[iRow] + move*element[j];
1834              if (newActivity<rowLower[iRow]-primalTolerance||
1835                  newActivity>rowUpper[iRow]+primalTolerance) {
1836                good=false;
1837                break;
1838              }
1839            }
1840            if (good) {
1841              newSolution[iColumn] = newValue;
1842              newSolutionValue += move*cost;
1843              int j;
1844              for (j=columnStart[iColumn];
1845                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
1846                int iRow = row[j];
1847                rowActivity[iRow] += move*element[j];
1848              }
1849            } else {
1850              move=0.0;
1851            }
1852          }
1853        }
1854      }
1855    }
1856    // Just in case of some stupidity
1857    double objOffset=0.0;
1858    solver->getDblParam(OsiObjOffset,objOffset);
1859    newSolutionValue = -objOffset;
1860    for ( i=0 ; i<numberColumns ; i++ )
1861      newSolutionValue += objective[i]*newSolution[i];
1862    newSolutionValue *= direction;
1863    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
1864    if (newSolutionValue<solutionValue) {
1865      // paranoid check
1866      memset(rowActivity,0,numberRows*sizeof(double));
1867      for (i=0;i<numberColumns;i++) {
1868        int j;
1869        double value = newSolution[i];
1870        if (value) {
1871          for (j=columnStart[i];
1872               j<columnStart[i]+columnLength[i];j++) {
1873            int iRow=row[j];
1874            rowActivity[iRow] += value*element[j];
1875          }
1876        }
1877      }
1878      // check was approximately feasible
1879      bool feasible=true;
1880      for (i=0;i<numberRows;i++) {
1881        if(rowActivity[i]<rowLower[i]) {
1882          if (rowActivity[i]<rowLower[i]-1000.0*primalTolerance)
1883            feasible = false;
1884        } else if(rowActivity[i]>rowUpper[i]) {
1885          if (rowActivity[i]>rowUpper[i]+1000.0*primalTolerance)
1886            feasible = false;
1887        }
1888      }
1889      if (feasible) {
1890        // new solution
1891        memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
1892        solutionValue = newSolutionValue;
1893        //printf("** Solution of %g found by rounding\n",newSolutionValue);
1894        returnCode=1;
1895      } else {
1896        // Can easily happen
1897        //printf("Debug CbcRounding giving bad solution\n");
1898      }
1899    }
1900  }
1901#ifdef NEW_ROUNDING
1902  if (!returnCode) {
1903#if 0
1904    // back to starting point
1905    memcpy(newSolution,solution,numberColumns*sizeof(double));
1906    memset(rowActivity,0,numberRows*sizeof(double));
1907    for (i=0;i<numberColumns;i++) {
1908      int j;
1909      double value = newSolution[i];
1910      if (value<lower[i]) {
1911        value=lower[i];
1912        newSolution[i]=value;
1913      } else if (value>upper[i]) {
1914        value=upper[i];
1915        newSolution[i]=value;
1916      }
1917      if (value) {
1918        for (j=columnStart[i];
1919             j<columnStart[i]+columnLength[i];j++) {
1920          int iRow=row[j];
1921          rowActivity[iRow] += value*element[j];
1922        }
1923      }
1924    }
1925    // check was feasible - if not adjust (cleaning may move)
1926    for (i=0;i<numberRows;i++) {
1927      if(rowActivity[i]<rowLower[i]) {
1928        //assert (rowActivity[i]>rowLower[i]-1000.0*primalTolerance);
1929        rowActivity[i]=rowLower[i];
1930      } else if(rowActivity[i]>rowUpper[i]) {
1931        //assert (rowActivity[i]<rowUpper[i]+1000.0*primalTolerance);
1932        rowActivity[i]=rowUpper[i];
1933      }
1934    }
1935#endif
1936    int * candidate = new int [numberColumns];
1937    int nCandidate=0;
1938    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1939      bool isInteger = (integerType[iColumn]!=0);
1940      if (isInteger) {
1941        double currentValue = newSolution[iColumn];
1942        if (fabs(currentValue-floor(currentValue+0.5))>1.0e-8)
1943          candidate[nCandidate++]=iColumn;
1944      }
1945    }
1946    if (true) {
1947      // Rounding as in Berthold
1948      while (nCandidate) {
1949        double infeasibility =1.0e-7;
1950        int iRow=-1;
1951        for (i=0;i<numberRows;i++) {
1952          double value=0.0;
1953          if(rowActivity[i]<rowLower[i]) {
1954            value = rowLower[i]-rowActivity[i];
1955          } else if(rowActivity[i]>rowUpper[i]) {
1956            value = rowActivity[i]-rowUpper[i];
1957          }
1958          if (value>infeasibility) {
1959            infeasibility = value;
1960            iRow=i;
1961          }
1962        }
1963        if (iRow>=0) {
1964          // infeasible
1965        } else {
1966          // feasible
1967        }
1968      }
1969    } else {
1970      // Shifting as in Berthold
1971    }
1972    delete [] candidate;
1973  }
1974#endif
1975  delete [] newSolution;
1976  delete [] rowActivity;
1977  return returnCode;
1978}
1979// update model
1980void CbcRounding::setModel(CbcModel * model)
1981{
1982  model_ = model;
1983  // Get a copy of original matrix (and by row for rounding);
1984  assert(model_->solver());
1985  matrix_ = *model_->solver()->getMatrixByCol();
1986  matrixByRow_ = *model_->solver()->getMatrixByRow();
1987  // make sure model okay for heuristic
1988  validate();
1989}
1990// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
1991void 
1992CbcRounding::validate() 
1993{
1994  if (model_&&when()<10) {
1995    if (model_->numberIntegers()!=
1996        model_->numberObjects()&&(model_->numberObjects()||
1997                                  (model_->specialOptions()&1024)==0))
1998      setWhen(0);
1999  }
2000#ifdef NEW_ROUNDING
2001  int numberColumns = matrix_.getNumCols();
2002  down_ = new unsigned short [numberColumns];
2003  up_ = new unsigned short [numberColumns];
2004  equal_ = new unsigned short [numberColumns];
2005  // Column copy
2006  const double * element = matrix_.getElements();
2007  const int * row = matrix_.getIndices();
2008  const CoinBigIndex * columnStart = matrix_.getVectorStarts();
2009  const int * columnLength = matrix_.getVectorLengths();
2010  const double * rowLower = model.solver()->getRowLower();
2011  const double * rowUpper = model.solver()->getRowUpper();
2012  for (int i=0;i<numberColumns;i++) {
2013    int down=0;
2014    int up=0;
2015    int equal=0;
2016    if (columnLength[i]>65535) {
2017      equal[0]=65535; 
2018      break; // unlikely to work
2019    }
2020    for (CoinBigIndex j=columnStart[i];
2021         j<columnStart[i]+columnLength[i];j++) {
2022      int iRow=row[j];
2023      if (rowLower[iRow]>-1.0e20&&rowUpper[iRow]<1.0e20) {
2024        equal++;
2025      } else if (element[j]>0.0) {
2026        if (rowUpper[iRow]<1.0e20)
2027          up++;
2028        else
2029          down--;
2030      } else {
2031        if (rowLower[iRow]>-1.0e20)
2032          up++;
2033        else
2034          down--;
2035      }
2036    }
2037    down_[i] = (unsigned short) down;
2038    up_[i] = (unsigned short) up;
2039    equal_[i] = (unsigned short) equal;
2040  }
2041#else
2042  down_ = NULL;
2043  up_ = NULL;
2044  equal_ = NULL;
2045#endif 
2046}
2047
2048// Default Constructor
2049CbcHeuristicPartial::CbcHeuristicPartial() 
2050  :CbcHeuristic()
2051{
2052  fixPriority_ = 10000;
2053}
2054
2055// Constructor from model
2056CbcHeuristicPartial::CbcHeuristicPartial(CbcModel & model, int fixPriority, int numberNodes)
2057  :CbcHeuristic(model)
2058{
2059  fixPriority_ = fixPriority;
2060  setNumberNodes(numberNodes);
2061  validate();
2062}
2063
2064// Destructor
2065CbcHeuristicPartial::~CbcHeuristicPartial ()
2066{
2067}
2068
2069// Clone
2070CbcHeuristic *
2071CbcHeuristicPartial::clone() const
2072{
2073  return new CbcHeuristicPartial(*this);
2074}
2075// Create C++ lines to get to current state
2076void 
2077CbcHeuristicPartial::generateCpp( FILE * fp) 
2078{
2079  CbcHeuristicPartial other;
2080  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
2081  fprintf(fp,"3  CbcHeuristicPartial partial(*cbcModel);\n");
2082  CbcHeuristic::generateCpp(fp,"partial");
2083  if (fixPriority_!=other.fixPriority_)
2084    fprintf(fp,"3  partial.setFixPriority(%d);\n",fixPriority_);
2085  else
2086    fprintf(fp,"4  partial.setFixPriority(%d);\n",fixPriority_);
2087  fprintf(fp,"3  cbcModel->addHeuristic(&partial);\n");
2088}
2089//#define NEW_PARTIAL
2090// Copy constructor
2091CbcHeuristicPartial::CbcHeuristicPartial(const CbcHeuristicPartial & rhs)
2092:
2093  CbcHeuristic(rhs),
2094  fixPriority_(rhs.fixPriority_)
2095{
2096}
2097
2098// Assignment operator
2099CbcHeuristicPartial & 
2100CbcHeuristicPartial::operator=( const CbcHeuristicPartial& rhs)
2101{
2102  if (this!=&rhs) {
2103    CbcHeuristic::operator=(rhs);
2104    fixPriority_ = rhs.fixPriority_;
2105  }
2106  return *this;
2107}
2108
2109// Resets stuff if model changes
2110void 
2111CbcHeuristicPartial::resetModel(CbcModel * model)
2112{
2113  model_=model;
2114  // Get a copy of original matrix (and by row for partial);
2115  assert(model_->solver());
2116  validate();
2117}
2118// See if partial will give solution
2119// Sets value of solution
2120// Assumes rhs for original matrix still okay
2121// At present only works with integers
2122// Fix values if asked for
2123// Returns 1 if solution, 0 if not
2124int
2125CbcHeuristicPartial::solution(double & solutionValue,
2126                      double * betterSolution)
2127{
2128  // Return if already done
2129  if (fixPriority_<0)
2130    return 0; // switched off
2131  const double * hotstartSolution = model_->hotstartSolution();
2132  const int * hotstartPriorities = model_->hotstartPriorities();
2133  if (!hotstartSolution)
2134    return 0;
2135  OsiSolverInterface * solver = model_->solver();
2136 
2137  int numberIntegers = model_->numberIntegers();
2138  const int * integerVariable = model_->integerVariable();
2139 
2140  OsiSolverInterface * newSolver = model_->continuousSolver()->clone();
2141  const double * colLower = newSolver->getColLower();
2142  const double * colUpper = newSolver->getColUpper();
2143
2144  double primalTolerance;
2145  solver->getDblParam(OsiPrimalTolerance,primalTolerance);
2146   
2147  int i;
2148  int numberFixed=0;
2149  int returnCode=0;
2150
2151  for (i=0;i<numberIntegers;i++) {
2152    int iColumn=integerVariable[i];
2153    if (abs(hotstartPriorities[iColumn])<=fixPriority_) {
2154      double value = hotstartSolution[iColumn];
2155      double lower = colLower[iColumn];
2156      double upper = colUpper[iColumn];
2157      value = CoinMax(value,lower);
2158      value = CoinMin(value,upper);
2159      if (fabs(value-floor(value+0.5))<1.0e-8) {
2160        value = floor(value+0.5);
2161        newSolver->setColLower(iColumn,value);
2162        newSolver->setColUpper(iColumn,value);
2163        numberFixed++;
2164      }
2165    }
2166  }
2167  if (numberFixed>numberIntegers/5-100000000) {
2168#ifdef COIN_DEVELOP
2169    printf("%d integers fixed\n",numberFixed);
2170#endif
2171    returnCode = smallBranchAndBound(newSolver,numberNodes_,betterSolution,solutionValue,
2172                                     model_->getCutoff(),"CbcHeuristicPartial");
2173    if (returnCode<0)
2174      returnCode=0; // returned on size
2175    //printf("return code %d",returnCode);
2176    if ((returnCode&2)!=0) {
2177      // could add cut
2178      returnCode &= ~2;
2179      //printf("could add cut with %d elements (if all 0-1)\n",nFix);
2180    } else {
2181      //printf("\n");
2182    }
2183  }
2184  fixPriority_=-1; // switch off
2185 
2186  delete newSolver;
2187  return returnCode;
2188}
2189// update model
2190void CbcHeuristicPartial::setModel(CbcModel * model)
2191{
2192  model_ = model;
2193  assert(model_->solver());
2194  // make sure model okay for heuristic
2195  validate();
2196}
2197// Validate model i.e. sets when_ to 0 if necessary (may be NULL)
2198void 
2199CbcHeuristicPartial::validate() 
2200{
2201  if (model_&&when()<10) {
2202    if (model_->numberIntegers()!=
2203        model_->numberObjects())
2204      setWhen(0);
2205  }
2206}
2207
2208// Default Constructor
2209CbcSerendipity::CbcSerendipity() 
2210  :CbcHeuristic()
2211{
2212}
2213
2214// Constructor from model
2215CbcSerendipity::CbcSerendipity(CbcModel & model)
2216  :CbcHeuristic(model)
2217{
2218}
2219
2220// Destructor
2221CbcSerendipity::~CbcSerendipity ()
2222{
2223}
2224
2225// Clone
2226CbcHeuristic *
2227CbcSerendipity::clone() const
2228{
2229  return new CbcSerendipity(*this);
2230}
2231// Create C++ lines to get to current state
2232void 
2233CbcSerendipity::generateCpp( FILE * fp) 
2234{
2235  fprintf(fp,"0#include \"CbcHeuristic.hpp\"\n");
2236  fprintf(fp,"3  CbcSerendipity serendipity(*cbcModel);\n");
2237  CbcHeuristic::generateCpp(fp,"serendipity");
2238  fprintf(fp,"3  cbcModel->addHeuristic(&serendipity);\n");
2239}
2240
2241// Copy constructor
2242CbcSerendipity::CbcSerendipity(const CbcSerendipity & rhs)
2243:
2244  CbcHeuristic(rhs)
2245{
2246}
2247
2248// Assignment operator
2249CbcSerendipity & 
2250CbcSerendipity::operator=( const CbcSerendipity& rhs)
2251{
2252  if (this!=&rhs) {
2253    CbcHeuristic::operator=(rhs);
2254  }
2255  return *this;
2256}
2257
2258// Returns 1 if solution, 0 if not
2259int
2260CbcSerendipity::solution(double & solutionValue,
2261                         double * betterSolution)
2262{
2263  if (!model_)
2264    return 0;
2265  if (!inputSolution_) {
2266    // get information on solver type
2267    OsiAuxInfo * auxInfo = model_->solver()->getAuxiliaryInfo();
2268    OsiBabSolver * auxiliaryInfo = dynamic_cast< OsiBabSolver *> (auxInfo);
2269    if (auxiliaryInfo) {
2270      return auxiliaryInfo->solution(solutionValue,betterSolution,model_->solver()->getNumCols());
2271    } else {
2272      return 0;
2273    }
2274  } else {
2275    int numberColumns = model_->getNumCols();
2276    double value =inputSolution_[numberColumns];
2277    int returnCode=0;
2278    if (value<solutionValue) {
2279      solutionValue = value;
2280      memcpy(betterSolution,inputSolution_,numberColumns*sizeof(double));
2281      returnCode=1;
2282    }
2283    delete [] inputSolution_;
2284    inputSolution_=NULL;
2285    model_ = NULL; // switch off
2286    return returnCode;
2287  }
2288}
2289// update model
2290void CbcSerendipity::setModel(CbcModel * model)
2291{
2292  model_ = model;
2293}
2294// Resets stuff if model changes
2295void 
2296CbcSerendipity::resetModel(CbcModel * model)
2297{
2298  model_ = model;
2299}
2300 
Note: See TracBrowser for help on using the repository browser.