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

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

Cleaned up a little bit.

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