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

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

changes

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