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

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

some changes e.g. DINS added and a bit of printing

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