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

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

out some printing

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