limp-cbc-0.3.2.0: cbits/coin/CbcHeuristicRENS.cpp
// $Id: CbcHeuristicRENS.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).
// edwin 12/5/09 carved out of CbcHeuristicRINS
#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 "CbcHeuristicRENS.hpp"
#include "CoinWarmStartBasis.hpp"
#include "CoinSort.hpp"
#include "CbcBranchActual.hpp"
#include "CbcStrategy.hpp"
#include "CglPreProcess.hpp"
// Default Constructor
CbcHeuristicRENS::CbcHeuristicRENS()
: CbcHeuristic()
{
numberTries_ = 0;
rensType_ = 0;
whereFrom_ = 256 + 1;
}
// Constructor with model - assumed before cuts
CbcHeuristicRENS::CbcHeuristicRENS(CbcModel & model)
: CbcHeuristic(model)
{
numberTries_ = 0;
rensType_ = 0;
whereFrom_ = 256 + 1;
}
// Destructor
CbcHeuristicRENS::~CbcHeuristicRENS ()
{
}
// Clone
CbcHeuristic *
CbcHeuristicRENS::clone() const
{
return new CbcHeuristicRENS(*this);
}
// Assignment operator
CbcHeuristicRENS &
CbcHeuristicRENS::operator=( const CbcHeuristicRENS & rhs)
{
if (this != &rhs) {
CbcHeuristic::operator=(rhs);
numberTries_ = rhs.numberTries_;
rensType_ = rhs.rensType_;
}
return *this;
}
// Copy constructor
CbcHeuristicRENS::CbcHeuristicRENS(const CbcHeuristicRENS & rhs)
:
CbcHeuristic(rhs),
numberTries_(rhs.numberTries_),
rensType_(rhs.rensType_)
{
}
// Resets stuff if model changes
void
CbcHeuristicRENS::resetModel(CbcModel * )
{
}
int
CbcHeuristicRENS::solution(double & solutionValue,
double * betterSolution)
{
int returnCode = 0;
const double * bestSolution = model_->bestSolution();
if ((numberTries_&&(rensType_&16)==0) || numberTries_>1 || (when() < 2 && bestSolution))
return 0;
numberTries_++;
double saveFractionSmall=fractionSmall_;
OsiSolverInterface * solver = model_->solver();
int numberIntegers = model_->numberIntegers();
const int * integerVariable = model_->integerVariable();
OsiSolverInterface * newSolver = cloneBut(3); // was model_->continuousSolver()->clone();
const double * currentSolution = newSolver->getColSolution();
int type = rensType_&15;
if (type<12)
newSolver->resolve();
double direction = newSolver->getObjSense();
double cutoff=model_->getCutoff();
newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
//cutoff *= direction;
double gap = cutoff - newSolver->getObjValue() * direction ;
double tolerance;
newSolver->getDblParam(OsiDualTolerance, tolerance) ;
if ((gap > 0.0 || !newSolver->isProvenOptimal())&&type<12) {
gap += 100.0 * tolerance;
int nFix = newSolver->reducedCostFix(gap);
if (nFix) {
char line [200];
sprintf(line, "Reduced cost fixing fixed %d variables", nFix);
model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
<< line
<< CoinMessageEol;
}
} else if (type<12) {
return 0; // finished?
}
int numberColumns = solver->getNumCols();
double * dj = CoinCopyOfArray(solver->getReducedCost(),numberColumns);
double djTolerance = (type!=1) ? -1.0e30 : 1.0e-4;
const double * colLower = newSolver->getColLower();
const double * colUpper = newSolver->getColUpper();
double * contribution = NULL;
int numberFixed = 0;
if (type==3) {
double total=0.0;
int n=0;
CoinWarmStartBasis * basis =
dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
if (basis&&basis->getNumArtificial()) {
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]&&
basis->getStructStatus(iColumn) !=
CoinWarmStartBasis::basic) {
n++;
total += fabs(dj[iColumn]);
}
}
if (n)
djTolerance = (0.01*total)/static_cast<double>(n);
delete basis;
}
} else if (type>=5&&type<=12) {
/* 5 fix sets at one
6 fix on dj but leave unfixed SOS slacks
7 fix sets at one but use pi
8 fix all at zero but leave unfixed SOS slacks
9 as 8 but only fix all at zero if just one in set nonzero
10 fix all "stable" ones
11 fix all "stable" ones - approach 2
12 layered approach
*/
// SOS type fixing
bool fixSets = (type==5)||(type==7)||(type==10)||(type==11);
CoinWarmStartBasis * basis =
dynamic_cast<CoinWarmStartBasis *>(solver->getWarmStart()) ;
if (basis&&basis->getNumArtificial()) {
//const double * rowLower = solver->getRowLower();
const double * rowUpper = solver->getRowUpper();
int numberRows = solver->getNumRows();
// Column copy
const CoinPackedMatrix * matrix = solver->getMatrixByCol();
const double * element = matrix->getElements();
const int * row = matrix->getIndices();
const CoinBigIndex * columnStart = matrix->getVectorStarts();
const int * columnLength = matrix->getVectorLengths();
double * bestDj = new double [numberRows];
for (int i=0;i<numberRows;i++) {
if (rowUpper[i]==1.0)
bestDj[i]=1.0e20;
else
bestDj[i]=1.0e30;
}
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
CoinBigIndex j;
if (currentSolution[iColumn]>1.0e-6&&
currentSolution[iColumn]<0.999999) {
for (j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
if (bestDj[iRow]<1.0e30) {
if (element[j] != 1.0)
bestDj[iRow]=1.0e30;
else
bestDj[iRow]=1.0e25;
}
}
} else if ( basis->getStructStatus(iColumn) !=
CoinWarmStartBasis::basic) {
for (j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
if (bestDj[iRow]<1.0e25) {
if (element[j] != 1.0)
bestDj[iRow]=1.0e30;
else
bestDj[iRow]=CoinMin(fabs(dj[iColumn]),bestDj[iRow]);
}
}
}
}
}
// Just leave one slack in each set
{
const double * objective = newSolver->getObjCoefficients();
int * best = new int [numberRows];
double * cheapest = new double[numberRows];
for (int i=0;i<numberRows;i++) {
best[i]=-1;
cheapest[i]=COIN_DBL_MAX;
}
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
if (columnLength[iColumn]==1) {
CoinBigIndex j = columnStart[iColumn];
int iRow = row[j];
if (bestDj[iRow]<1.0e30) {
double obj = direction*objective[iColumn];
if (obj<cheapest[iRow]) {
cheapest[iRow]=obj;
best[iRow]=iColumn;
}
}
}
}
}
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
if (columnLength[iColumn]==1) {
CoinBigIndex j = columnStart[iColumn];
int iRow = row[j];
if (bestDj[iRow]<1.0e30) {
if (best[iRow]!=-1&&iColumn!=best[iRow]) {
newSolver->setColUpper(iColumn,0.0);
}
}
}
}
}
delete [] best;
delete [] cheapest;
}
int nSOS=0;
double * sort = new double [numberRows];
const double * pi = newSolver->getRowPrice();
if (type==12) {
contribution = new double [numberRows];
for (int i=0;i<numberRows;i++) {
if (bestDj[i]<1.0e30)
contribution[i]=0.0;
else
contribution[i]=-1.0;
}
}
for (int i=0;i<numberRows;i++) {
if (bestDj[i]<1.0e30) {
if (type==5)
sort[nSOS++]=bestDj[i];
else if (type==7)
sort[nSOS++]=-fabs(pi[i]);
else
sort[nSOS++]=fabs(pi[i]);
}
}
if (10*nSOS>8*numberRows) {
if (type<10) {
std::sort(sort,sort+nSOS);
int last = static_cast<int>(nSOS*0.9*fractionSmall_);
double tolerance = sort[last];
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
CoinBigIndex j;
if (currentSolution[iColumn]<=1.0e-6||
currentSolution[iColumn]>=0.999999) {
if (fixSets) {
for (j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
double useDj;
if (type==5)
useDj = bestDj[iRow];
else if (type==7)
useDj= -fabs(pi[iRow]);
else
useDj= fabs(pi[iRow]);
if (bestDj[iRow]<1.0e30&&useDj>=tolerance) {
numberFixed++;
if (currentSolution[iColumn]<=1.0e-6)
newSolver->setColUpper(iColumn,0.0);
else if (currentSolution[iColumn]>=0.999999)
newSolver->setColLower(iColumn,1.0);
}
}
} else if (columnLength[iColumn]==1) {
// leave more slacks
int iRow = row[columnStart[iColumn]];
if (bestDj[iRow]<1.0e30) {
// fake dj
dj[iColumn] *= 0.000001;
}
} else if (type==8||type==9) {
if (currentSolution[iColumn]<=1.0e-6) {
if (type==8) {
dj[iColumn] *= 1.0e6;
} else {
bool fix=false;
for (j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
if (bestDj[iRow]<1.0e25) {
fix=true;
break;
}
}
if (fix) {
dj[iColumn] *= 1.0e6;
}
}
} else {
dj[iColumn] *= 0.000001;
}
}
}
}
}
if (fixSets)
djTolerance = 1.0e30;
} else if (type==10) {
double * saveUpper = new double [numberRows];
memset(saveUpper,0,numberRows*sizeof(double));
char * mark = new char [numberColumns];
char * nonzero = new char [numberColumns];
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
CoinBigIndex j;
for (j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
saveUpper[iRow] += element[j];
}
}
}
double sum=0.0;
double sumRhs=0.0;
const double * rowUpper = newSolver->getRowUpper();
for (int i=0;i<numberRows;i++) {
if (bestDj[i]>=1.0e30) {
sum += saveUpper[i];
sumRhs += rowUpper[i];
}
}
double averagePerSet = sum/static_cast<double>(numberRows);
// allow this extra
double factor = averagePerSet*fractionSmall_*numberRows;
factor = 1.0+factor/sumRhs;
fractionSmall_ = 0.5;
memcpy(saveUpper,rowUpper,numberRows*sizeof(double));
// loosen up
for (int i=0;i<numberRows;i++) {
if (bestDj[i]>=1.0e30) {
newSolver->setRowUpper(i,factor*saveUpper[i]);
}
}
newSolver->resolve();
const double * solution = newSolver->getColSolution();
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
mark[iColumn]=0;
nonzero[iColumn]=0;
if (colUpper[iColumn]>colLower[iColumn]&&
solution[iColumn]>0.9999)
mark[iColumn]=1;
else if (solution[iColumn]>0.00001)
nonzero[iColumn]=1;
}
// slightly small
for (int i=0;i<numberRows;i++) {
if (bestDj[i]>=1.0e30) {
newSolver->setRowUpper(i,saveUpper[i]*0.9999);
}
}
newSolver->resolve();
int nCheck=2;
if (newSolver->isProvenOptimal()) {
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]&&
solution[iColumn]>0.9999)
mark[iColumn]++;
else if (solution[iColumn]>0.00001)
nonzero[iColumn]=1;
}
} else {
nCheck=1;
}
// correct values
for (int i=0;i<numberRows;i++) {
if (bestDj[i]>=1.0e30) {
newSolver->setRowUpper(i,saveUpper[i]);
}
}
newSolver->resolve();
int nFixed=0;
int nFixedToZero=0;
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
if (solution[iColumn]>0.9999&&mark[iColumn]==nCheck) {
newSolver->setColLower(iColumn,1.0);
nFixed++;
} else if (!mark[iColumn]&&!nonzero[iColumn]&&
columnLength[iColumn]>1&&solution[iColumn]<0.00001) {
newSolver->setColUpper(iColumn,0.0);
nFixedToZero++;
}
}
}
char line[100];
sprintf(line,"Heuristic %s fixed %d to one (and %d to zero)",
heuristicName(),
nFixed,nFixedToZero);
model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
<< line
<< CoinMessageEol;
delete [] mark;
delete []nonzero;
delete [] saveUpper;
numberFixed=numberColumns;
djTolerance = 1.0e30;
} else if (type==11) {
double * saveUpper = CoinCopyOfArray(newSolver->getRowUpper(),numberRows);
char * mark = new char [numberColumns];
char * nonzero = new char [numberColumns];
// save basis and solution
CoinWarmStartBasis * basis = dynamic_cast<CoinWarmStartBasis*>(newSolver->getWarmStart()) ;
assert(basis != NULL);
double * saveSolution =
CoinCopyOfArray(newSolver->getColSolution(),
numberColumns);
double factors[] = {1.1,1.05,1.01,0.98};
int nPass = (sizeof(factors)/sizeof(double))-1;
double factor=factors[0];
double proportion = fractionSmall_;
fractionSmall_ = 0.5;
// loosen up
for (int i=0;i<numberRows;i++) {
if (bestDj[i]>=1.0e30) {
newSolver->setRowUpper(i,factor*saveUpper[i]);
}
}
bool takeHint;
OsiHintStrength strength;
newSolver->getHintParam(OsiDoDualInResolve, takeHint, strength);
newSolver->setHintParam(OsiDoDualInResolve, false, OsiHintDo);
newSolver->resolve();
newSolver->setHintParam(OsiDoDualInResolve, true, OsiHintDo);
const double * solution = newSolver->getColSolution();
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
mark[iColumn]=0;
nonzero[iColumn]=0;
if (colUpper[iColumn]>colLower[iColumn]&&
solution[iColumn]>0.9999)
mark[iColumn]=1;
else if (solution[iColumn]>0.00001)
nonzero[iColumn]=1;
}
int nCheck=2;
for (int iPass=0;iPass<nPass;iPass++) {
// smaller
factor = factors[iPass+1];
for (int i=0;i<numberRows;i++) {
if (bestDj[i]>=1.0e30) {
newSolver->setRowUpper(i,saveUpper[i]*factor);
}
}
newSolver->resolve();
if (newSolver->isProvenOptimal()) {
nCheck++;
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]&&
solution[iColumn]>0.9999)
mark[iColumn]++;
else if (solution[iColumn]>0.00001)
nonzero[iColumn]++;
}
}
}
// correct values
for (int i=0;i<numberRows;i++) {
if (bestDj[i]>=1.0e30) {
newSolver->setRowUpper(i,saveUpper[i]);
}
}
newSolver->setColSolution(saveSolution);
delete [] saveSolution;
newSolver->setWarmStart(basis);
delete basis ;
newSolver->setHintParam(OsiDoDualInResolve, takeHint, strength);
newSolver->resolve();
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]&&
solution[iColumn]>0.9999)
mark[iColumn]++;
else if (solution[iColumn]>0.00001)
nonzero[iColumn]++;
}
int nFixed=0;
int numberSetsToFix = static_cast<int>(nSOS*(1.0-proportion));
int * mixed = new int[numberRows];
memset(mixed,0,numberRows*sizeof(int));
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
int iSOS=-1;
for (CoinBigIndex j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
if (bestDj[iRow]<1.0e25) {
iSOS=iRow;
break;
}
}
if (iSOS>=0) {
int numberTimesAtOne = mark[iColumn];
int numberTimesNonZero = nonzero[iColumn]+
numberTimesAtOne;
if (numberTimesAtOne<nCheck&&
numberTimesNonZero) {
mixed[iSOS]+=
CoinMin(numberTimesNonZero,
nCheck-numberTimesNonZero);
}
}
}
}
int nFix=0;
for (int i=0;i<numberRows;i++) {
if (bestDj[i]<1.0e25) {
sort[nFix] = -bestDj[i]+1.0e8*mixed[i];
mixed[nFix++]=i;
}
}
CoinSort_2(sort,sort+nFix,mixed);
nFix = CoinMin(nFix,numberSetsToFix);
memset(sort,0,sizeof(double)*numberRows);
for (int i=0;i<nFix;i++)
sort[mixed[i]]=1.0;
delete [] mixed;
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
if (solution[iColumn]>0.9999) {
int iSOS=-1;
for (CoinBigIndex j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
if (bestDj[iRow]<1.0e25) {
iSOS=iRow;
break;
}
}
if (iSOS>=0&&sort[iSOS]) {
newSolver->setColLower(iColumn,1.0);
nFixed++;
}
}
}
}
char line[100];
sprintf(line,"Heuristic %s fixed %d to one (%d sets)",
heuristicName(),
nFixed,nSOS);
model_->messageHandler()->message(CBC_FPUMP1, model_->messages())
<< line
<< CoinMessageEol;
delete [] mark;
delete [] nonzero;
delete [] saveUpper;
numberFixed=numberColumns;
djTolerance = 1.0e30;
}
}
delete basis;
delete [] sort;
delete [] bestDj;
if (10*nSOS<=8*numberRows) {
// give up
delete [] contribution;
delete newSolver;
return 0;
}
}
}
// Do dj to get right number
if (type==4||type==6||(type>7&&type<10)) {
double * sort = new double [numberColumns];
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
sort[iColumn]=1.0e30;
if (colUpper[iColumn]>colLower[iColumn]) {
sort[iColumn] = fabs(dj[iColumn]);
}
}
std::sort(sort,sort+numberColumns);
int last = static_cast<int>(numberColumns*fractionSmall_);
djTolerance = CoinMax(sort[last],1.0e-5);
delete [] sort;
} else if (type==12) {
// Do layered in a different way
int numberRows = solver->getNumRows();
// Column copy
const CoinPackedMatrix * matrix = newSolver->getMatrixByCol();
const double * element = matrix->getElements();
const int * row = matrix->getIndices();
const CoinBigIndex * columnStart = matrix->getVectorStarts();
const int * columnLength = matrix->getVectorLengths();
int * whichRow = new int[numberRows];
int * whichSet = new int [numberColumns];
int nSOS=0;
for (int i=0;i<numberRows;i++) {
whichRow[i]=0;
if (!contribution[i])
nSOS++;
}
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
whichSet[iColumn]=-2;
if (colUpper[iColumn]>colLower[iColumn]) {
CoinBigIndex j;
double sum=0.0;
int iSOS=-1;
int n=0;
for (j = columnStart[iColumn];
j < columnStart[iColumn] + columnLength[iColumn]; j++) {
int iRow = row[j];
if (contribution[iRow]>=0.0) {
iSOS=iRow;
n++;
} else {
sum += fabs(element[j]);
}
}
if (n>1)
COIN_DETAIL_PRINT(printf("Too many SOS entries (%d) for column %d\n",
n,iColumn));
if (sum) {
assert (iSOS>=0);
contribution[iSOS] += sum;
whichRow[iSOS]++;
whichSet[iColumn]=iSOS;
} else {
whichSet[iColumn]=iSOS+numberRows;
}
}
}
int * chunk = new int [numberRows];
for (int i=0;i<numberRows;i++) {
chunk[i]=-1;
if (whichRow[i]) {
contribution[i]= - contribution[i]/static_cast<double>(whichRow[i]);
} else {
contribution[i] = COIN_DBL_MAX;
}
whichRow[i]=i;
}
newSolver->setDblParam(OsiDualObjectiveLimit, 1.0e100);
double * saveLower = CoinCopyOfArray(colLower,numberColumns);
double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
CoinSort_2(contribution,contribution+numberRows,whichRow);
// Set do nothing solution
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if(whichSet[iColumn]>=numberRows)
newSolver->setColLower(iColumn,1.0);
}
newSolver->resolve();
int nChunk = (nSOS+9)/10;
int nPass=0;
int inChunk=0;
for (int i=0;i<nSOS;i++) {
chunk[whichRow[i]]=nPass;
inChunk++;
if (inChunk==nChunk) {
inChunk=0;
// last two together
if (i+nChunk<nSOS)
nPass++;
}
}
// adjust
nPass++;
for (int iPass=0;iPass<nPass;iPass++) {
// fix last chunk and unfix this chunk
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
int iSOS = whichSet[iColumn];
if (iSOS>=0) {
if (iSOS>=numberRows)
iSOS-=numberRows;
if (chunk[iSOS]==iPass-1&&betterSolution[iColumn]>0.9999) {
newSolver->setColLower(iColumn,1.0);
} else if (chunk[iSOS]==iPass) {
newSolver->setColLower(iColumn,saveLower[iColumn]);
newSolver->setColUpper(iColumn,saveUpper[iColumn]);
}
}
}
// solve
returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
model_->getCutoff(), "CbcHeuristicRENS");
if (returnCode < 0) {
returnCode = 0; // returned on size
break;
} else if ((returnCode&1)==0) {
// no good
break;
}
}
if ((returnCode&2) != 0) {
// could add cut
returnCode &= ~2;
}
delete [] chunk;
delete [] saveLower;
delete [] saveUpper;
delete [] whichRow;
delete [] whichSet;
delete [] contribution;
delete newSolver;
return returnCode;
}
double primalTolerance;
solver->getDblParam(OsiPrimalTolerance, primalTolerance);
int i;
int numberTightened = 0;
int numberAtBound = 0;
int numberContinuous = numberColumns - numberIntegers;
for (i = 0; i < numberIntegers; i++) {
int iColumn = integerVariable[i];
double value = currentSolution[iColumn];
double lower = colLower[iColumn];
double upper = colUpper[iColumn];
value = CoinMax(value, lower);
value = CoinMin(value, upper);
double djValue=dj[iColumn]*direction;
#define RENS_FIX_ONLY_LOWER
#ifndef RENS_FIX_ONLY_LOWER
if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
value = floor(value + 0.5);
if (value == lower || value == upper)
numberAtBound++;
newSolver->setColLower(iColumn, value);
newSolver->setColUpper(iColumn, value);
numberFixed++;
} else if (colUpper[iColumn] - colLower[iColumn] >= 2.0) {
numberTightened++;
newSolver->setColLower(iColumn, floor(value));
newSolver->setColUpper(iColumn, ceil(value));
}
#else
if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
floor(value + 0.5) == lower &&
djValue > djTolerance ) {
value = floor(value + 0.5);
numberAtBound++;
newSolver->setColLower(iColumn, value);
newSolver->setColUpper(iColumn, value);
numberFixed++;
} else if (fabs(value - floor(value + 0.5)) < 1.0e-8 &&
floor(value + 0.5) == upper &&
-djValue > djTolerance && (djTolerance > 0.0||type==2)) {
value = floor(value + 0.5);
numberAtBound++;
newSolver->setColLower(iColumn, value);
newSolver->setColUpper(iColumn, value);
numberFixed++;
} else if (colUpper[iColumn] - colLower[iColumn] >= 2.0 &&
djTolerance <0.0) {
numberTightened++;
if (fabs(value - floor(value + 0.5)) < 1.0e-8) {
value = floor(value + 0.5);
if (value < upper) {
newSolver->setColLower(iColumn, CoinMax(value - 1.0, lower));
newSolver->setColUpper(iColumn, CoinMin(value + 1.0, upper));
} else {
newSolver->setColLower(iColumn, upper - 1.0);
}
} else {
newSolver->setColLower(iColumn, floor(value));
newSolver->setColUpper(iColumn, ceil(value));
}
}
#endif
}
delete [] dj;
if (numberFixed > numberIntegers / 5) {
if ( numberFixed < numberColumns / 5) {
#define RENS_FIX_CONTINUOUS
#ifdef RENS_FIX_CONTINUOUS
const double * colLower = newSolver->getColLower();
//const double * colUpper = newSolver->getColUpper();
int nAtLb = 0;
double sumDj = 0.0;
const double * dj = newSolver->getReducedCost();
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (!newSolver->isInteger(iColumn)) {
double value = currentSolution[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 = currentSolution[iColumn];
if (value < colLower[iColumn] + 1.0e-8) {
double djValue = dj[iColumn] * direction;
if (djValue > threshold) {
sort[nFix2] = -djValue;
which[nFix2++] = iColumn;
}
}
}
}
CoinSort_2(sort, sort + nFix2, which);
nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
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 fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
numberFixed, numberTightened, numberAtBound, nFix2);
#endif
}
#endif
}
#ifdef COIN_DEVELOP
printf("%d integers fixed and %d tightened\n", numberFixed, numberTightened);
#endif
returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
model_->getCutoff(), "CbcHeuristicRENS");
if (returnCode < 0 || returnCode == 0) {
#ifdef RENS_FIX_CONTINUOUS
if (numberContinuous > numberIntegers && numberFixed >= numberColumns / 5) {
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 = currentSolution[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 = currentSolution[iColumn];
if (value < colLower[iColumn] + 1.0e-8) {
double djValue = dj[iColumn] * direction;
if (djValue > threshold) {
sort[nFix2] = -djValue;
which[nFix2++] = iColumn;
}
}
}
}
CoinSort_2(sort, sort + nFix2, which);
nFix2 = CoinMin(nFix2, (numberColumns - numberFixed) / 2);
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 fixed (%d tightened) (%d at bound), and %d continuous fixed at lb\n",
numberFixed, numberTightened, numberAtBound, nFix2);
#endif
}
returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
model_->getCutoff(), "CbcHeuristicRENS");
}
#endif
if (returnCode < 0 || returnCode == 0) {
// Do passes fixing up those >0.9 and
// down those < 0.05
#define RENS_PASS 3
//#define KEEP_GOING
#ifdef KEEP_GOING
double * saveLower = CoinCopyOfArray(colLower,numberColumns);
double * saveUpper = CoinCopyOfArray(colUpper,numberColumns);
bool badPass=false;
int nSolved=0;
#endif
for (int iPass=0;iPass<RENS_PASS;iPass++) {
int nFixed=0;
int nFixedAlready=0;
int nFixedContinuous=0;
for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
if (colUpper[iColumn]>colLower[iColumn]) {
if (newSolver->isInteger(iColumn)) {
double value = currentSolution[iColumn];
double fixTo = floor(value+0.1);
if (fixTo>value || value-fixTo < 0.05) {
// above 0.9 or below 0.05
nFixed++;
newSolver->setColLower(iColumn, fixTo);
newSolver->setColUpper(iColumn, fixTo);
}
}
} else if (newSolver->isInteger(iColumn)) {
nFixedAlready++;
} else {
nFixedContinuous++;
}
}
#ifdef CLP_INVESTIGATE2
printf("%d more integers fixed (total %d) plus %d continuous\n",
nFixed,nFixed+nFixedAlready,nFixedContinuous);
#endif
#ifdef KEEP_GOING
if (nFixed) {
newSolver->resolve();
if (!newSolver->isProvenOptimal()) {
badPass=true;
break;
} else {
nSolved++;
memcpy(saveLower,colLower,numberColumns*sizeof(double));
memcpy(saveUpper,colUpper,numberColumns*sizeof(double));
}
} else {
break;
}
#else
if (nFixed) {
newSolver->resolve();
if (!newSolver->isProvenOptimal()) {
returnCode=0;
break;
}
returnCode = smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
model_->getCutoff(), "CbcHeuristicRENS");
} else {
returnCode=0;
}
if (returnCode>=0)
break;
}
if (returnCode < 0)
returnCode = 0; // returned on size
#endif
}
#ifdef KEEP_GOING
if (badPass) {
newSolver->setColLower(saveLower);
newSolver->setColUpper(saveUpper);
newSolver->resolve();
}
delete [] saveLower;
delete [] saveUpper;
if (nSolved)
returnCode =
smallBranchAndBound(newSolver, numberNodes_, betterSolution, solutionValue,
model_->getCutoff(), "CbcHeuristicRENS");
else
returnCode=0;
}
#endif
}
//printf("return code %d",returnCode);
if ((returnCode&2) != 0) {
// could add cut
returnCode &= ~2;
#ifdef COIN_DEVELOP
if (!numberTightened && numberFixed == numberAtBound)
printf("could add cut with %d elements\n", numberFixed);
#endif
} else {
//printf("\n");
}
}
//delete [] whichRow;
//delete [] contribution;
delete newSolver;
fractionSmall_ = saveFractionSmall;
return returnCode;
}
// update model
void CbcHeuristicRENS::setModel(CbcModel * model)
{
model_ = model;
}