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

Last change on this file since 1020 was 1020, checked in by forrest, 11 years ago

tweaks to heuristics

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