limp-cbc-0.3.2.0: cbits/coin/CbcHeuristicRINS.cpp
/* $Id: CbcHeuristicRINS.cpp 1902 2013-04-10 16:58:16Z stefan $ */
// Copyright (C) 2006, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
#if defined(_MSC_VER)
// Turn off compiler warning about long names
# pragma warning(disable:4786)
#endif
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <cfloat>
#include "OsiSolverInterface.hpp"
#include "CbcModel.hpp"
#include "CbcMessage.hpp"
#include "CbcHeuristicRINS.hpp"
#include "CbcBranchActual.hpp"
#include "CbcStrategy.hpp"
#include "CglPreProcess.hpp"
// Default Constructor
CbcHeuristicRINS::CbcHeuristicRINS()
: CbcHeuristic()
{
numberSolutions_ = 0;
numberSuccesses_ = 0;
numberTries_ = 0;
stateOfFixing_ = 0;
shallowDepth_ = 0;
lastNode_ = -999999;
howOften_ = 100;
decayFactor_ = 0.5;
used_ = NULL;
whereFrom_ = 1 + 8 + 16 + 255 * 256;
whereFrom_ = 1 + 8 + 255 * 256;
}
// Constructor with model - assumed before cuts
CbcHeuristicRINS::CbcHeuristicRINS(CbcModel & model)
: CbcHeuristic(model)
{
numberSolutions_ = 0;
numberSuccesses_ = 0;
numberTries_ = 0;
stateOfFixing_ = 0;
shallowDepth_ = 0;
lastNode_ = -999999;
howOften_ = 100;
decayFactor_ = 0.5;
assert(model.solver());
int numberColumns = model.solver()->getNumCols();
used_ = new char[numberColumns];
memset(used_, 0, numberColumns);
whereFrom_ = 1 + 8 + 16 + 255 * 256;
whereFrom_ = 1 + 8 + 255 * 256;
}
// Destructor
CbcHeuristicRINS::~CbcHeuristicRINS ()
{
delete [] used_;
}
// Clone
CbcHeuristic *
CbcHeuristicRINS::clone() const
{
return new CbcHeuristicRINS(*this);
}
// Assignment operator
CbcHeuristicRINS &
CbcHeuristicRINS::operator=( const CbcHeuristicRINS & rhs)
{
if (this != &rhs) {
CbcHeuristic::operator=(rhs);
numberSolutions_ = rhs.numberSolutions_;
howOften_ = rhs.howOften_;
numberSuccesses_ = rhs.numberSuccesses_;
numberTries_ = rhs.numberTries_;
stateOfFixing_ = rhs.stateOfFixing_;
lastNode_ = rhs.lastNode_;
delete [] used_;
if (model_ && rhs.used_) {
int numberColumns = model_->solver()->getNumCols();
used_ = new char[numberColumns];
memcpy(used_, rhs.used_, numberColumns);
} else {
used_ = NULL;
}
}
return *this;
}
// Create C++ lines to get to current state
void
CbcHeuristicRINS::generateCpp( FILE * fp)
{
CbcHeuristicRINS other;
fprintf(fp, "0#include \"CbcHeuristicRINS.hpp\"\n");
fprintf(fp, "3 CbcHeuristicRINS heuristicRINS(*cbcModel);\n");
CbcHeuristic::generateCpp(fp, "heuristicRINS");
if (howOften_ != other.howOften_)
fprintf(fp, "3 heuristicRINS.setHowOften(%d);\n", howOften_);
else
fprintf(fp, "4 heuristicRINS.setHowOften(%d);\n", howOften_);
fprintf(fp, "3 cbcModel->addHeuristic(&heuristicRINS);\n");
}
// Copy constructor
CbcHeuristicRINS::CbcHeuristicRINS(const CbcHeuristicRINS & rhs)
:
CbcHeuristic(rhs),
numberSolutions_(rhs.numberSolutions_),
howOften_(rhs.howOften_),
numberSuccesses_(rhs.numberSuccesses_),
numberTries_(rhs.numberTries_),
stateOfFixing_(rhs.stateOfFixing_),
lastNode_(rhs.lastNode_)
{
if (model_ && rhs.used_) {
int numberColumns = model_->solver()->getNumCols();
used_ = new char[numberColumns];
memcpy(used_, rhs.used_, numberColumns);
} else {
used_ = NULL;
}
}
// Resets stuff if model changes
void
CbcHeuristicRINS::resetModel(CbcModel * /*model*/)
{
//CbcHeuristic::resetModel(model);
delete [] used_;
stateOfFixing_ = 0;
if (model_ && used_) {
int numberColumns = model_->solver()->getNumCols();
used_ = new char[numberColumns];
memset(used_, 0, numberColumns);
} else {
used_ = NULL;
}
}
/*
First tries setting a variable to better value. If feasible then
tries setting others. If not feasible then tries swaps
Returns 1 if solution, 0 if not */
int
CbcHeuristicRINS::solution(double & solutionValue,
double * betterSolution)
{
numCouldRun_++;
int returnCode = 0;
const double * bestSolution = model_->bestSolution();
if (!bestSolution)
return 0; // No solution found yet
if (numberSolutions_ < model_->getSolutionCount()) {
// new solution - add info
numberSolutions_ = model_->getSolutionCount();
int numberIntegers = model_->numberIntegers();
const int * integerVariable = model_->integerVariable();
int i;
for (i = 0; i < numberIntegers; i++) {
int iColumn = integerVariable[i];
const OsiObject * object = model_->object(i);
// get original bounds
double originalLower;
double originalUpper;
getIntegerInformation( object, originalLower, originalUpper);
double value = bestSolution[iColumn];
if (value < originalLower) {
value = originalLower;
} else if (value > originalUpper) {
value = originalUpper;
}
double nearest = floor(value + 0.5);
// if away from lower bound mark that fact
if (nearest > originalLower) {
used_[iColumn] = 1;
}
}
}
int numberNodes = model_->getNodeCount();
if (howOften_ == 100) {
if (numberNodes < lastNode_ + 12)
return 0;
// Do at 50 and 100
if ((numberNodes > 40 && numberNodes <= 50) || (numberNodes > 90 && numberNodes < 100))
numberNodes = howOften_;
}
// Allow for infeasible nodes - so do anyway after a bit
if (howOften_ >= 100 && numberNodes >= lastNode_ + 2*howOften_) {
numberNodes = howOften_;
}
if ((numberNodes % howOften_) == 0 && (model_->getCurrentPassNumber() == 1 ||
model_->getCurrentPassNumber() == 999999)) {
lastNode_ = model_->getNodeCount();
OsiSolverInterface * solver = model_->solver();
int numberIntegers = model_->numberIntegers();
const int * integerVariable = model_->integerVariable();
const double * currentSolution = solver->getColSolution();
const int * used = model_->usedInSolution();
OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
int numberColumns = newSolver->getNumCols();
int numberContinuous = numberColumns - numberIntegers;
double primalTolerance;
solver->getDblParam(OsiPrimalTolerance, primalTolerance);
int i;
int nFix = 0;
for (i = 0; i < numberIntegers; i++) {
int iColumn = integerVariable[i];
const OsiObject * object = model_->object(i);
// get original bounds
double originalLower;
double originalUpper;
getIntegerInformation( object, originalLower, originalUpper);
double valueInt = bestSolution[iColumn];
if (valueInt < originalLower) {
valueInt = originalLower;
} else if (valueInt > originalUpper) {
valueInt = originalUpper;
}
if (fabs(currentSolution[iColumn] - valueInt) < 10.0*primalTolerance) {
double nearest = floor(valueInt + 0.5);
/*
shallowDepth_
0 - normal
1 - only fix if at lb
2 - only fix if not at lb
3 - only fix if at lb and !used
*/
bool fix=false;
switch (shallowDepth_) {
case 0:
fix = true;
break;
case 1:
if (nearest==originalLower)
fix = true;
break;
case 2:
if (nearest!=originalLower)
fix = true;
break;
case 3:
if (nearest==originalLower && !used[iColumn])
fix = true;
break;
}
if (fix) {
newSolver->setColLower(iColumn, nearest);
newSolver->setColUpper(iColumn, nearest);
nFix++;
}
}
}
int divisor = 0;
if (5*nFix > numberIntegers) {
if (numberContinuous > 2*numberIntegers && ((nFix*10 < numberColumns &&
!numRuns_ && numberTries_ > 2) ||
stateOfFixing_)) {
#define RINS_FIX_CONTINUOUS
#ifdef RINS_FIX_CONTINUOUS
const double * colLower = newSolver->getColLower();
//const double * colUpper = newSolver->getColUpper();
int nAtLb = 0;
//double sumDj=0.0;
const double * dj = newSolver->getReducedCost();
double direction = newSolver->getObjSense();
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (!newSolver->isInteger(iColumn)) {
double value = bestSolution[iColumn];
if (value < colLower[iColumn] + 1.0e-8) {
//double djValue = dj[iColumn]*direction;
nAtLb++;
//sumDj += djValue;
}
}
}
if (nAtLb) {
// fix some continuous
double * sort = new double[nAtLb];
int * which = new int [nAtLb];
//double threshold = CoinMax((0.01*sumDj)/static_cast<double>(nAtLb),1.0e-6);
int nFix2 = 0;
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (!newSolver->isInteger(iColumn)) {
double value = bestSolution[iColumn];
if (value < colLower[iColumn] + 1.0e-8) {
double djValue = dj[iColumn] * direction;
if (djValue > 1.0e-6) {
sort[nFix2] = -djValue;
which[nFix2++] = iColumn;
}
}
}
}
CoinSort_2(sort, sort + nFix2, which);
divisor = 4;
if (stateOfFixing_ > 0)
divisor = stateOfFixing_;
else if (stateOfFixing_ < -1)
divisor = (-stateOfFixing_) - 1;
nFix2 = CoinMin(nFix2, (numberColumns - nFix) / divisor);
for (int i = 0; i < nFix2; i++) {
int iColumn = which[i];
newSolver->setColUpper(iColumn, colLower[iColumn]);
}
delete [] sort;
delete [] which;
#ifdef CLP_INVESTIGATE2
printf("%d integers have same value, and %d continuous fixed at lb\n",
nFix, nFix2);
#endif
}
#endif
}
//printf("%d integers have same value\n",nFix);
returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
model_->getCutoff(), "CbcHeuristicRINS");
if (returnCode < 0) {
returnCode = 0; // returned on size
if (divisor) {
stateOfFixing_ = - divisor; // say failed
} else if (numberContinuous > 2*numberIntegers &&
!numRuns_ && numberTries_ > 2) {
stateOfFixing_ = -4; //start fixing
}
} else {
numRuns_++;
if (divisor)
stateOfFixing_ = divisor; // say small enough
}
if ((returnCode&1) != 0)
numberSuccesses_++;
//printf("return code %d",returnCode);
if ((returnCode&2) != 0) {
// could add cut
returnCode &= ~2;
//printf("could add cut with %d elements (if all 0-1)\n",nFix);
} else {
//printf("\n");
}
}
numberTries_++;
if ((numberTries_ % 10) == 0 && numberSuccesses_*3 < numberTries_)
howOften_ += static_cast<int> (howOften_ * decayFactor_);
delete newSolver;
}
return returnCode;
}
// update model
void CbcHeuristicRINS::setModel(CbcModel * model)
{
model_ = model;
// Get a copy of original matrix
assert(model_->solver());
delete [] used_;
int numberColumns = model->solver()->getNumCols();
used_ = new char[numberColumns];
memset(used_, 0, numberColumns);
}