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

Last change on this file since 884 was 884, checked in by jpgoncal, 11 years ago

Draft code for heuristic to control heuristics.

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