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

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

compiler warnings, deterministic parallel and stability

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