source: branches/heur/Cbc/src/CbcHeuristic.cpp @ 893

Last change on this file since 893 was 893, checked in by jpgoncal, 13 years ago

Added function to print info.

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