limp-cbc-0.3.2.0: cbits/coin/CbcSolverHeuristics.cpp
/* $Id: CbcSolverHeuristics.cpp 1902 2013-04-10 16:58:16Z stefan $ */
// Copyright (C) 2007, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
/*! \file CbcSolverHeuristics.cpp
\brief Second level routines for the cbc stand-alone solver.
*/
#include "CbcConfig.h"
#include "CoinPragma.hpp"
#include "CoinTime.hpp"
#include "OsiClpSolverInterface.hpp"
#include "ClpPresolve.hpp"
#include "CbcOrClpParam.hpp"
#include "CbcModel.hpp"
#include "CbcHeuristicLocal.hpp"
#include "CbcHeuristicPivotAndFix.hpp"
//#include "CbcHeuristicPivotAndComplement.hpp"
#include "CbcHeuristicRandRound.hpp"
#include "CbcHeuristicGreedy.hpp"
#include "CbcHeuristicFPump.hpp"
#include "CbcHeuristicRINS.hpp"
#include "CbcHeuristicDiveCoefficient.hpp"
#include "CbcHeuristicDiveFractional.hpp"
#include "CbcHeuristicDiveGuided.hpp"
#include "CbcHeuristicDiveVectorLength.hpp"
#include "CbcHeuristicDivePseudoCost.hpp"
#include "CbcHeuristicDiveLineSearch.hpp"
#include "CbcStrategy.hpp"
#include "OsiAuxInfo.hpp"
#include "ClpSimplexOther.hpp"
// Crunch down model
void
crunchIt(ClpSimplex * model)
{
#ifdef JJF_ZERO
model->dual();
#else
int numberColumns = model->numberColumns();
int numberRows = model->numberRows();
// Use dual region
double * rhs = model->dualRowSolution();
int * whichRow = new int[3*numberRows];
int * whichColumn = new int[2*numberColumns];
int nBound;
ClpSimplex * small = static_cast<ClpSimplexOther *> (model)->crunch(rhs, whichRow, whichColumn,
nBound, false, false);
if (small) {
small->dual();
if (small->problemStatus() == 0) {
model->setProblemStatus(0);
static_cast<ClpSimplexOther *> (model)->afterCrunch(*small, whichRow, whichColumn, nBound);
} else if (small->problemStatus() != 3) {
model->setProblemStatus(1);
} else {
if (small->problemStatus() == 3) {
// may be problems
small->computeObjectiveValue();
model->setObjectiveValue(small->objectiveValue());
model->setProblemStatus(3);
} else {
model->setProblemStatus(3);
}
}
delete small;
} else {
model->setProblemStatus(1);
}
delete [] whichRow;
delete [] whichColumn;
#endif
}
/*
On input
doAction - 0 just fix in original and return NULL
1 return fixed non-presolved solver
2 as one but use presolve Inside this
3 use presolve and fix ones with large cost
? do heuristics and set best solution
? do BAB and just set best solution
10+ then use lastSolution and relax a few
-2 cleanup afterwards if using 2
On output - number fixed
*/
OsiClpSolverInterface *
fixVubs(CbcModel & model, int skipZero2,
int & doAction,
CoinMessageHandler * /*generalMessageHandler*/,
const double * lastSolution, double dextra[6],
int extra[5])
{
if (doAction == 11 && !lastSolution)
lastSolution = model.bestSolution();
assert (((doAction >= 0 && doAction <= 3) && !lastSolution) || (doAction == 11 && lastSolution));
double fractionIntFixed = dextra[3];
double fractionFixed = dextra[4];
double fixAbove = dextra[2];
double fixAboveValue = (dextra[5] > 0.0) ? dextra[5] : 1.0;
#ifdef COIN_DETAIL
double time1 = CoinCpuTime();
#endif
int leaveIntFree = extra[1];
OsiSolverInterface * originalSolver = model.solver();
OsiClpSolverInterface * originalClpSolver = dynamic_cast< OsiClpSolverInterface*> (originalSolver);
ClpSimplex * originalLpSolver = originalClpSolver->getModelPtr();
int * originalColumns = NULL;
OsiClpSolverInterface * clpSolver;
ClpSimplex * lpSolver;
ClpPresolve pinfo;
assert(originalSolver->getObjSense() > 0);
if (doAction == 2 || doAction == 3) {
double * saveLB = NULL;
double * saveUB = NULL;
int numberColumns = originalLpSolver->numberColumns();
if (fixAbove > 0.0) {
#ifdef COIN_DETAIL
double time1 = CoinCpuTime();
#endif
originalClpSolver->initialSolve();
COIN_DETAIL_PRINT(printf("first solve took %g seconds\n", CoinCpuTime() - time1));
double * columnLower = originalLpSolver->columnLower() ;
double * columnUpper = originalLpSolver->columnUpper() ;
const double * solution = originalLpSolver->primalColumnSolution();
saveLB = CoinCopyOfArray(columnLower, numberColumns);
saveUB = CoinCopyOfArray(columnUpper, numberColumns);
const double * objective = originalLpSolver->getObjCoefficients() ;
int iColumn;
int nFix = 0;
int nArt = 0;
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
if (objective[iColumn] > fixAbove) {
if (solution[iColumn] < columnLower[iColumn] + 1.0e-8) {
columnUpper[iColumn] = columnLower[iColumn];
nFix++;
} else {
nArt++;
}
} else if (objective[iColumn] < -fixAbove) {
if (solution[iColumn] > columnUpper[iColumn] - 1.0e-8) {
columnLower[iColumn] = columnUpper[iColumn];
nFix++;
} else {
nArt++;
}
}
}
COIN_DETAIL_PRINT(printf("%d artificials fixed, %d left as in solution\n", nFix, nArt));
lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e-8, true, 10);
if (!lpSolver || doAction == 2) {
// take off fixing in original
memcpy(columnLower, saveLB, numberColumns*sizeof(double));
memcpy(columnUpper, saveUB, numberColumns*sizeof(double));
}
delete [] saveLB;
delete [] saveUB;
if (!lpSolver) {
// try again
pinfo.destroyPresolve();
lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e-8, true, 10);
assert (lpSolver);
}
} else {
lpSolver = pinfo.presolvedModel(*originalLpSolver, 1.0e-8, true, 10);
assert (lpSolver);
}
clpSolver = new OsiClpSolverInterface(lpSolver, true);
assert(lpSolver == clpSolver->getModelPtr());
numberColumns = lpSolver->numberColumns();
originalColumns = CoinCopyOfArray(pinfo.originalColumns(), numberColumns);
doAction = 1;
} else {
OsiSolverInterface * solver = originalSolver->clone();
clpSolver = dynamic_cast< OsiClpSolverInterface*> (solver);
lpSolver = clpSolver->getModelPtr();
}
// Tighten bounds
lpSolver->tightenPrimalBounds(0.0, 11, true);
int numberColumns = clpSolver->getNumCols() ;
double * saveColumnLower = CoinCopyOfArray(lpSolver->columnLower(), numberColumns);
double * saveColumnUpper = CoinCopyOfArray(lpSolver->columnUpper(), numberColumns);
//char generalPrint[200];
const double *objective = lpSolver->getObjCoefficients() ;
double *columnLower = lpSolver->columnLower() ;
double *columnUpper = lpSolver->columnUpper() ;
int numberRows = clpSolver->getNumRows();
int iRow, iColumn;
// Row copy
CoinPackedMatrix matrixByRow(*clpSolver->getMatrixByRow());
const double * elementByRow = matrixByRow.getElements();
const int * column = matrixByRow.getIndices();
const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
const int * rowLength = matrixByRow.getVectorLengths();
// Column copy
CoinPackedMatrix matrixByCol(*clpSolver->getMatrixByCol());
//const double * element = matrixByCol.getElements();
const int * row = matrixByCol.getIndices();
const CoinBigIndex * columnStart = matrixByCol.getVectorStarts();
const int * columnLength = matrixByCol.getVectorLengths();
const double * rowLower = clpSolver->getRowLower();
const double * rowUpper = clpSolver->getRowUpper();
// Get maximum size of VUB tree
// otherColumn is one fixed to 0 if this one zero
int nEl = matrixByCol.getNumElements();
CoinBigIndex * fixColumn = new CoinBigIndex [numberColumns+1];
int * otherColumn = new int [nEl];
int * fix = new int[numberColumns];
char * mark = new char [numberColumns];
memset(mark, 0, numberColumns);
int numberInteger = 0;
int numberOther = 0;
fixColumn[0] = 0;
double large = lpSolver->largeValue(); // treat bounds > this as infinite
#ifndef NDEBUG
double large2 = 1.0e10 * large;
#endif
double tolerance = lpSolver->primalTolerance();
int * check = new int[numberRows];
for (iRow = 0; iRow < numberRows; iRow++) {
check[iRow] = -2; // don't check
if (rowLower[iRow] < -1.0e6 && rowUpper[iRow] > 1.0e6)
continue;// unlikely
// possible row
int numberPositive = 0;
int iPositive = -1;
int numberNegative = 0;
int iNegative = -1;
CoinBigIndex rStart = rowStart[iRow];
CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
CoinBigIndex j;
int kColumn;
for (j = rStart; j < rEnd; ++j) {
double value = elementByRow[j];
kColumn = column[j];
if (columnUpper[kColumn] > columnLower[kColumn]) {
if (value > 0.0) {
numberPositive++;
iPositive = kColumn;
} else {
numberNegative++;
iNegative = kColumn;
}
}
}
if (numberPositive == 1 && numberNegative == 1)
check[iRow] = -1; // try both
if (numberPositive == 1 && rowLower[iRow] > -1.0e20)
check[iRow] = iPositive;
else if (numberNegative == 1 && rowUpper[iRow] < 1.0e20)
check[iRow] = iNegative;
}
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
fix[iColumn] = -1;
if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e-8) {
if (clpSolver->isInteger(iColumn))
numberInteger++;
if (columnLower[iColumn] == 0.0) {
bool infeasible = false;
fix[iColumn] = 0;
// fake upper bound
double saveUpper = columnUpper[iColumn];
columnUpper[iColumn] = 0.0;
for (CoinBigIndex i = columnStart[iColumn];
i < columnStart[iColumn] + columnLength[iColumn]; i++) {
iRow = row[i];
if (check[iRow] != -1 && check[iRow] != iColumn)
continue; // unlikely
// possible row
int infiniteUpper = 0;
int infiniteLower = 0;
double maximumUp = 0.0;
double maximumDown = 0.0;
double newBound;
CoinBigIndex rStart = rowStart[iRow];
CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
CoinBigIndex j;
int kColumn;
// Compute possible lower and upper ranges
for (j = rStart; j < rEnd; ++j) {
double value = elementByRow[j];
kColumn = column[j];
if (value > 0.0) {
if (columnUpper[kColumn] >= large) {
++infiniteUpper;
} else {
maximumUp += columnUpper[kColumn] * value;
}
if (columnLower[kColumn] <= -large) {
++infiniteLower;
} else {
maximumDown += columnLower[kColumn] * value;
}
} else if (value < 0.0) {
if (columnUpper[kColumn] >= large) {
++infiniteLower;
} else {
maximumDown += columnUpper[kColumn] * value;
}
if (columnLower[kColumn] <= -large) {
++infiniteUpper;
} else {
maximumUp += columnLower[kColumn] * value;
}
}
}
// Build in a margin of error
maximumUp += 1.0e-8 * fabs(maximumUp);
maximumDown -= 1.0e-8 * fabs(maximumDown);
double maxUp = maximumUp + infiniteUpper * 1.0e31;
double maxDown = maximumDown - infiniteLower * 1.0e31;
if (maxUp <= rowUpper[iRow] + tolerance &&
maxDown >= rowLower[iRow] - tolerance) {
//printf("Redundant row in vubs %d\n",iRow);
} else {
if (maxUp < rowLower[iRow] - 100.0*tolerance ||
maxDown > rowUpper[iRow] + 100.0*tolerance) {
infeasible = true;
break;
}
double lower = rowLower[iRow];
double upper = rowUpper[iRow];
for (j = rStart; j < rEnd; ++j) {
double value = elementByRow[j];
kColumn = column[j];
double nowLower = columnLower[kColumn];
double nowUpper = columnUpper[kColumn];
if (value > 0.0) {
// positive value
if (lower > -large) {
if (!infiniteUpper) {
assert(nowUpper < large2);
newBound = nowUpper +
(lower - maximumUp) / value;
// relax if original was large
if (fabs(maximumUp) > 1.0e8)
newBound -= 1.0e-12 * fabs(maximumUp);
} else if (infiniteUpper == 1 && nowUpper > large) {
newBound = (lower - maximumUp) / value;
// relax if original was large
if (fabs(maximumUp) > 1.0e8)
newBound -= 1.0e-12 * fabs(maximumUp);
} else {
newBound = -COIN_DBL_MAX;
}
if (newBound > nowLower + 1.0e-12 && newBound > -large) {
// Tighten the lower bound
// check infeasible (relaxed)
if (nowUpper < newBound) {
if (nowUpper - newBound <
-100.0*tolerance) {
infeasible = true;
break;
}
}
}
}
if (upper < large) {
if (!infiniteLower) {
assert(nowLower > - large2);
newBound = nowLower +
(upper - maximumDown) / value;
// relax if original was large
if (fabs(maximumDown) > 1.0e8)
newBound += 1.0e-12 * fabs(maximumDown);
} else if (infiniteLower == 1 && nowLower < -large) {
newBound = (upper - maximumDown) / value;
// relax if original was large
if (fabs(maximumDown) > 1.0e8)
newBound += 1.0e-12 * fabs(maximumDown);
} else {
newBound = COIN_DBL_MAX;
}
if (newBound < nowUpper - 1.0e-12 && newBound < large) {
// Tighten the upper bound
// check infeasible (relaxed)
if (nowLower > newBound) {
if (newBound - nowLower <
-100.0*tolerance) {
infeasible = true;
break;
} else {
newBound = nowLower;
}
}
if (!newBound || (clpSolver->isInteger(kColumn) && newBound < 0.999)) {
// fix to zero
if (!mark[kColumn]) {
otherColumn[numberOther++] = kColumn;
mark[kColumn] = 1;
if (check[iRow] == -1)
check[iRow] = iColumn;
else
assert(check[iRow] == iColumn);
}
}
}
}
} else {
// negative value
if (lower > -large) {
if (!infiniteUpper) {
assert(nowLower < large2);
newBound = nowLower +
(lower - maximumUp) / value;
// relax if original was large
if (fabs(maximumUp) > 1.0e8)
newBound += 1.0e-12 * fabs(maximumUp);
} else if (infiniteUpper == 1 && nowLower < -large) {
newBound = (lower - maximumUp) / value;
// relax if original was large
if (fabs(maximumUp) > 1.0e8)
newBound += 1.0e-12 * fabs(maximumUp);
} else {
newBound = COIN_DBL_MAX;
}
if (newBound < nowUpper - 1.0e-12 && newBound < large) {
// Tighten the upper bound
// check infeasible (relaxed)
if (nowLower > newBound) {
if (newBound - nowLower <
-100.0*tolerance) {
infeasible = true;
break;
} else {
newBound = nowLower;
}
}
if (!newBound || (clpSolver->isInteger(kColumn) && newBound < 0.999)) {
// fix to zero
if (!mark[kColumn]) {
otherColumn[numberOther++] = kColumn;
mark[kColumn] = 1;
if (check[iRow] == -1)
check[iRow] = iColumn;
else
assert(check[iRow] == iColumn);
}
}
}
}
if (upper < large) {
if (!infiniteLower) {
assert(nowUpper < large2);
newBound = nowUpper +
(upper - maximumDown) / value;
// relax if original was large
if (fabs(maximumDown) > 1.0e8)
newBound -= 1.0e-12 * fabs(maximumDown);
} else if (infiniteLower == 1 && nowUpper > large) {
newBound = (upper - maximumDown) / value;
// relax if original was large
if (fabs(maximumDown) > 1.0e8)
newBound -= 1.0e-12 * fabs(maximumDown);
} else {
newBound = -COIN_DBL_MAX;
}
if (newBound > nowLower + 1.0e-12 && newBound > -large) {
// Tighten the lower bound
// check infeasible (relaxed)
if (nowUpper < newBound) {
if (nowUpper - newBound <
-100.0*tolerance) {
infeasible = true;
break;
}
}
}
}
}
}
}
}
for (int i = fixColumn[iColumn]; i < numberOther; i++)
mark[otherColumn[i]] = 0;
// reset bound unless infeasible
if (!infeasible || !clpSolver->isInteger(iColumn))
columnUpper[iColumn] = saveUpper;
else if (clpSolver->isInteger(iColumn))
columnLower[iColumn] = 1.0;
}
}
fixColumn[iColumn+1] = numberOther;
}
delete [] check;
delete [] mark;
// Now do reverse way
int * counts = new int [numberColumns];
CoinZeroN(counts, numberColumns);
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++)
counts[otherColumn[i]]++;
}
numberOther = 0;
CoinBigIndex * fixColumn2 = new CoinBigIndex [numberColumns+1];
int * otherColumn2 = new int [fixColumn[numberColumns]];
fixColumn2[0] = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
numberOther += counts[iColumn];
counts[iColumn] = 0;
fixColumn2[iColumn+1] = numberOther;
}
// Create other way
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
int jColumn = otherColumn[i];
CoinBigIndex put = fixColumn2[jColumn] + counts[jColumn];
counts[jColumn]++;
otherColumn2[put] = iColumn;
}
}
// get top layer i.e. those which are not fixed by any other
int kLayer = 0;
while (true) {
int numberLayered = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (fix[iColumn] == kLayer) {
for (int i = fixColumn2[iColumn]; i < fixColumn2[iColumn+1]; i++) {
int jColumn = otherColumn2[i];
if (fix[jColumn] == kLayer) {
fix[iColumn] = kLayer + 100;
}
}
}
if (fix[iColumn] == kLayer) {
numberLayered++;
}
}
if (numberLayered) {
kLayer += 100;
} else {
break;
}
}
for (int iPass = 0; iPass < 2; iPass++) {
for (int jLayer = 0; jLayer < kLayer; jLayer++) {
int check[] = { -1, 0, 1, 2, 3, 4, 5, 10, 50, 100, 500, 1000, 5000, 10000, COIN_INT_MAX};
int nCheck = static_cast<int> (sizeof(check) / sizeof(int));
int countsI[20];
int countsC[20];
assert (nCheck <= 20);
memset(countsI, 0, nCheck*sizeof(int));
memset(countsC, 0, nCheck*sizeof(int));
check[nCheck-1] = numberColumns;
int numberLayered = 0;
int numberInteger = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (fix[iColumn] == jLayer) {
numberLayered++;
int nFix = fixColumn[iColumn+1] - fixColumn[iColumn];
if (iPass) {
// just integers
nFix = 0;
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
int jColumn = otherColumn[i];
if (clpSolver->isInteger(jColumn))
nFix++;
}
}
int iFix;
for (iFix = 0; iFix < nCheck; iFix++) {
if (nFix <= check[iFix])
break;
}
assert (iFix < nCheck);
if (clpSolver->isInteger(iColumn)) {
numberInteger++;
countsI[iFix]++;
} else {
countsC[iFix]++;
}
}
}
#ifdef COIN_DETAIL
if (numberLayered) {
printf("%d (%d integer) at priority %d\n", numberLayered, numberInteger, 1 + (jLayer / 100));
char buffer[50];
for (int i = 1; i < nCheck; i++) {
if (countsI[i] || countsC[i]) {
if (i == 1)
sprintf(buffer, " == zero ");
else if (i < nCheck - 1)
sprintf(buffer, "> %6d and <= %6d ", check[i-1], check[i]);
else
sprintf(buffer, "> %6d ", check[i-1]);
printf("%s %8d integers and %8d continuous\n", buffer, countsI[i], countsC[i]);
}
}
}
#endif
}
}
delete [] counts;
// Now do fixing
{
// switch off presolve and up weight
ClpSolve solveOptions;
//solveOptions.setPresolveType(ClpSolve::presolveOff,0);
solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
//solveOptions.setSpecialOption(1,3,30); // sprint
int numberColumns = lpSolver->numberColumns();
int iColumn;
bool allSlack = true;
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
if (lpSolver->getColumnStatus(iColumn) == ClpSimplex::basic) {
allSlack = false;
break;
}
}
if (allSlack)
solveOptions.setSpecialOption(1, 2, 50); // idiot
lpSolver->setInfeasibilityCost(1.0e11);
lpSolver->defaultFactorizationFrequency();
if (doAction != 11)
lpSolver->initialSolve(solveOptions);
double * columnLower = lpSolver->columnLower();
double * columnUpper = lpSolver->columnUpper();
double * fullSolution = lpSolver->primalColumnSolution();
const double * dj = lpSolver->dualColumnSolution();
int iPass = 0;
#define MAXPROB 2
ClpSimplex models[MAXPROB];
int kPass = -1;
int kLayer = 0;
int skipZero = 0;
if (skipZero2 == -1)
skipZero2 = 40; //-1;
/* 0 fixed to 0 by choice
1 lb of 1 by choice
2 fixed to 0 by another
3 as 2 but this go
-1 free
*/
char * state = new char [numberColumns];
for (iColumn = 0; iColumn < numberColumns; iColumn++)
state[iColumn] = -1;
while (true) {
double largest = -0.1;
double smallest = 1.1;
int iLargest = -1;
int iSmallest = -1;
int atZero = 0;
int atOne = 0;
int toZero = 0;
int toOne = 0;
int numberFree = 0;
int numberGreater = 0;
columnLower = lpSolver->columnLower();
columnUpper = lpSolver->columnUpper();
fullSolution = lpSolver->primalColumnSolution();
if (doAction == 11) {
{
double * columnLower = lpSolver->columnLower();
double * columnUpper = lpSolver->columnUpper();
// lpSolver->dual();
memcpy(columnLower, saveColumnLower, numberColumns*sizeof(double));
memcpy(columnUpper, saveColumnUpper, numberColumns*sizeof(double));
// lpSolver->dual();
int iColumn;
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e-8) {
if (clpSolver->isInteger(iColumn)) {
double value = lastSolution[iColumn];
int iValue = static_cast<int> (value + 0.5);
assert (fabs(value - static_cast<double> (iValue)) < 1.0e-3);
assert (iValue >= columnLower[iColumn] &&
iValue <= columnUpper[iColumn]);
columnLower[iColumn] = iValue;
columnUpper[iColumn] = iValue;
}
}
}
lpSolver->initialSolve(solveOptions);
memcpy(columnLower, saveColumnLower, numberColumns*sizeof(double));
memcpy(columnUpper, saveColumnUpper, numberColumns*sizeof(double));
}
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnUpper[iColumn] > columnLower[iColumn] + 1.0e-8) {
if (clpSolver->isInteger(iColumn)) {
double value = lastSolution[iColumn];
int iValue = static_cast<int> (value + 0.5);
assert (fabs(value - static_cast<double> (iValue)) < 1.0e-3);
assert (iValue >= columnLower[iColumn] &&
iValue <= columnUpper[iColumn]);
if (!fix[iColumn]) {
if (iValue == 0) {
state[iColumn] = 0;
assert (!columnLower[iColumn]);
columnUpper[iColumn] = 0.0;
} else if (iValue == 1) {
state[iColumn] = 1;
columnLower[iColumn] = 1.0;
} else {
// leave fixed
columnLower[iColumn] = iValue;
columnUpper[iColumn] = iValue;
}
} else if (iValue == 0) {
state[iColumn] = 10;
columnUpper[iColumn] = 0.0;
} else {
// leave fixed
columnLower[iColumn] = iValue;
columnUpper[iColumn] = iValue;
}
}
}
}
int jLayer = 0;
int nFixed = -1;
int nTotalFixed = 0;
while (nFixed) {
nFixed = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
int jColumn = otherColumn[i];
if (columnUpper[jColumn]) {
bool canFix = true;
for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
int kColumn = otherColumn2[k];
if (state[kColumn] == 1) {
canFix = false;
break;
}
}
if (canFix) {
columnUpper[jColumn] = 0.0;
nFixed++;
}
}
}
}
}
nTotalFixed += nFixed;
jLayer += 100;
}
COIN_DETAIL_PRINT(printf("This fixes %d variables in lower priorities\n", nTotalFixed));
break;
}
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (!clpSolver->isInteger(iColumn) || fix[iColumn] > kLayer)
continue;
// skip if fixes nothing
if (fixColumn[iColumn+1] - fixColumn[iColumn] <= skipZero2)
continue;
double value = fullSolution[iColumn];
if (value > 1.00001) {
numberGreater++;
continue;
}
double lower = columnLower[iColumn];
double upper = columnUpper[iColumn];
if (lower == upper) {
if (lower)
atOne++;
else
atZero++;
continue;
}
if (value < 1.0e-7) {
toZero++;
columnUpper[iColumn] = 0.0;
state[iColumn] = 10;
continue;
}
if (value > 1.0 - 1.0e-7) {
toOne++;
columnLower[iColumn] = 1.0;
state[iColumn] = 1;
continue;
}
numberFree++;
// skip if fixes nothing
if (fixColumn[iColumn+1] - fixColumn[iColumn] <= skipZero)
continue;
if (value < smallest) {
smallest = value;
iSmallest = iColumn;
}
if (value > largest) {
largest = value;
iLargest = iColumn;
}
}
if (toZero || toOne)
COIN_DETAIL_PRINT(printf("%d at 0 fixed and %d at one fixed\n", toZero, toOne));
COIN_DETAIL_PRINT(printf("%d variables free, %d fixed to 0, %d to 1 - smallest %g, largest %g\n",
numberFree, atZero, atOne, smallest, largest));
if (numberGreater && !iPass)
COIN_DETAIL_PRINT(printf("%d variables have value > 1.0\n", numberGreater));
//skipZero2=0; // leave 0 fixing
int jLayer = 0;
int nFixed = -1;
int nTotalFixed = 0;
while (nFixed) {
nFixed = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
int jColumn = otherColumn[i];
if (columnUpper[jColumn]) {
bool canFix = true;
for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
int kColumn = otherColumn2[k];
if (state[kColumn] == 1) {
canFix = false;
break;
}
}
if (canFix) {
columnUpper[jColumn] = 0.0;
nFixed++;
}
}
}
}
}
nTotalFixed += nFixed;
jLayer += 100;
}
COIN_DETAIL_PRINT(printf("This fixes %d variables in lower priorities\n", nTotalFixed));
if (iLargest < 0 || numberFree <= leaveIntFree)
break;
double movement;
int way;
if (smallest <= 1.0 - largest && smallest < 0.2 && largest < fixAboveValue) {
columnUpper[iSmallest] = 0.0;
state[iSmallest] = 0;
movement = smallest;
way = -1;
} else {
columnLower[iLargest] = 1.0;
state[iLargest] = 1;
movement = 1.0 - largest;
way = 1;
}
double saveObj = lpSolver->objectiveValue();
iPass++;
kPass = iPass % MAXPROB;
models[kPass] = *lpSolver;
if (way == -1) {
// fix others
for (int i = fixColumn[iSmallest]; i < fixColumn[iSmallest+1]; i++) {
int jColumn = otherColumn[i];
if (state[jColumn] == -1) {
columnUpper[jColumn] = 0.0;
state[jColumn] = 3;
}
}
}
double maxCostUp = COIN_DBL_MAX;
objective = lpSolver->getObjCoefficients() ;
if (way == -1)
maxCostUp = (1.0 - movement) * objective[iSmallest];
lpSolver->setDualObjectiveLimit(saveObj + maxCostUp);
crunchIt(lpSolver);
double moveObj = lpSolver->objectiveValue() - saveObj;
COIN_DETAIL_PRINT(printf("movement %s was %g costing %g\n",
(way == -1) ? "down" : "up", movement, moveObj));
if (way == -1 && (moveObj >= maxCostUp || lpSolver->status())) {
// go up
columnLower = models[kPass].columnLower();
columnUpper = models[kPass].columnUpper();
columnLower[iSmallest] = 1.0;
columnUpper[iSmallest] = saveColumnUpper[iSmallest];
*lpSolver = models[kPass];
columnLower = lpSolver->columnLower();
columnUpper = lpSolver->columnUpper();
fullSolution = lpSolver->primalColumnSolution();
dj = lpSolver->dualColumnSolution();
columnLower[iSmallest] = 1.0;
columnUpper[iSmallest] = saveColumnUpper[iSmallest];
state[iSmallest] = 1;
// unfix others
for (int i = fixColumn[iSmallest]; i < fixColumn[iSmallest+1]; i++) {
int jColumn = otherColumn[i];
if (state[jColumn] == 3) {
columnUpper[jColumn] = saveColumnUpper[jColumn];
state[jColumn] = -1;
}
}
crunchIt(lpSolver);
}
models[kPass] = *lpSolver;
}
lpSolver->dual();
COIN_DETAIL_PRINT(printf("Fixing took %g seconds\n", CoinCpuTime() - time1));
columnLower = lpSolver->columnLower();
columnUpper = lpSolver->columnUpper();
fullSolution = lpSolver->primalColumnSolution();
dj = lpSolver->dualColumnSolution();
int * sort = new int[numberColumns];
double * dsort = new double[numberColumns];
int chunk = 20;
int iRelax = 0;
//double fractionFixed=6.0/8.0;
// relax while lots fixed
while (true) {
if (skipZero2 > 10 && doAction < 10)
break;
iRelax++;
int n = 0;
double sum0 = 0.0;
double sum00 = 0.0;
double sum1 = 0.0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (!clpSolver->isInteger(iColumn) || fix[iColumn] > kLayer)
continue;
// skip if fixes nothing
if (fixColumn[iColumn+1] - fixColumn[iColumn] == 0 && doAction < 10)
continue;
double djValue = dj[iColumn];
if (state[iColumn] == 1) {
assert (columnLower[iColumn]);
assert (fullSolution[iColumn] > 0.1);
if (djValue > 0.0) {
//printf("YY dj of %d at %g is %g\n",iColumn,value,djValue);
sum1 += djValue;
sort[n] = iColumn;
dsort[n++] = -djValue;
} else {
//printf("dj of %d at %g is %g\n",iColumn,value,djValue);
}
} else if (state[iColumn] == 0 || state[iColumn] == 10) {
assert (fullSolution[iColumn] < 0.1);
assert (!columnUpper[iColumn]);
double otherValue = 0.0;
int nn = 0;
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
int jColumn = otherColumn[i];
if (columnUpper[jColumn] == 0.0) {
if (dj[jColumn] < -1.0e-5) {
nn++;
otherValue += dj[jColumn]; // really need to look at rest
}
}
}
if (djValue < -1.0e-2 || otherValue < -1.0e-2) {
//printf("XX dj of %d at %g is %g - %d out of %d contribute %g\n",iColumn,value,djValue,
// nn,fixColumn[iColumn+1]-fixColumn[iColumn],otherValue);
if (djValue < 1.0e-8) {
sum0 -= djValue;
sum00 -= otherValue;
sort[n] = iColumn;
if (djValue < -1.0e-2)
dsort[n++] = djValue + otherValue;
else
dsort[n++] = djValue + 0.001 * otherValue;
}
} else {
//printf("dj of %d at %g is %g - no contribution from %d\n",iColumn,value,djValue,
// fixColumn[iColumn+1]-fixColumn[iColumn]);
}
}
}
CoinSort_2(dsort, dsort + n, sort);
double * originalColumnLower = saveColumnLower;
double * originalColumnUpper = saveColumnUpper;
double * lo = CoinCopyOfArray(columnLower, numberColumns);
double * up = CoinCopyOfArray(columnUpper, numberColumns);
for (int k = 0; k < CoinMin(chunk, n); k++) {
iColumn = sort[k];
state[iColumn] = -2;
}
memcpy(columnLower, originalColumnLower, numberColumns*sizeof(double));
memcpy(columnUpper, originalColumnUpper, numberColumns*sizeof(double));
int nFixed = 0;
int nFixed0 = 0;
int nFixed1 = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (state[iColumn] == 0 || state[iColumn] == 10) {
columnUpper[iColumn] = 0.0;
assert (lo[iColumn] == 0.0);
nFixed++;
nFixed0++;
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
int jColumn = otherColumn[i];
if (columnUpper[jColumn]) {
bool canFix = true;
for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
int kColumn = otherColumn2[k];
if (state[kColumn] == 1 || state[kColumn] == -2) {
canFix = false;
break;
}
}
if (canFix) {
columnUpper[jColumn] = 0.0;
assert (lo[jColumn] == 0.0);
nFixed++;
}
}
}
} else if (state[iColumn] == 1) {
columnLower[iColumn] = 1.0;
nFixed1++;
}
}
COIN_DETAIL_PRINT(printf("%d fixed %d orig 0 %d 1\n", nFixed, nFixed0, nFixed1));
int jLayer = 0;
nFixed = -1;
int nTotalFixed = 0;
while (nFixed) {
nFixed = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnUpper[iColumn] == 0.0 && fix[iColumn] == jLayer) {
for (int i = fixColumn[iColumn]; i < fixColumn[iColumn+1]; i++) {
int jColumn = otherColumn[i];
if (columnUpper[jColumn]) {
bool canFix = true;
for (int k = fixColumn2[jColumn]; k < fixColumn2[jColumn+1]; k++) {
int kColumn = otherColumn2[k];
if (state[kColumn] == 1 || state[kColumn] == -2) {
canFix = false;
break;
}
}
if (canFix) {
columnUpper[jColumn] = 0.0;
assert (lo[jColumn] == 0.0);
nFixed++;
}
}
}
}
}
nTotalFixed += nFixed;
jLayer += 100;
}
nFixed = 0;
int nFixedI = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (columnLower[iColumn] == columnUpper[iColumn]) {
if (clpSolver->isInteger(iColumn))
nFixedI++;
nFixed++;
}
}
COIN_DETAIL_PRINT(printf("This fixes %d variables in lower priorities - total %d (%d integer) - all target %d, int target %d\n",
nTotalFixed, nFixed, nFixedI, static_cast<int>(fractionFixed*numberColumns), static_cast<int> (fractionIntFixed*numberInteger)));
int nBad = 0;
int nRelax = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
if (lo[iColumn] < columnLower[iColumn] ||
up[iColumn] > columnUpper[iColumn]) {
COIN_DETAIL_PRINT(printf("bad %d old %g %g, new %g %g\n", iColumn, lo[iColumn], up[iColumn],
columnLower[iColumn], columnUpper[iColumn]));
nBad++;
}
if (lo[iColumn] > columnLower[iColumn] ||
up[iColumn] < columnUpper[iColumn]) {
nRelax++;
}
}
COIN_DETAIL_PRINT(printf("%d relaxed\n", nRelax));
if (iRelax > 20 && nRelax == chunk)
nRelax = 0;
if (iRelax > 50)
nRelax = 0;
assert (!nBad);
delete [] lo;
delete [] up;
lpSolver->primal(1);
if (nFixed < fractionFixed*numberColumns || nFixedI < fractionIntFixed*numberInteger || !nRelax)
break;
}
delete [] state;
delete [] sort;
delete [] dsort;
}
delete [] fix;
delete [] fixColumn;
delete [] otherColumn;
delete [] otherColumn2;
delete [] fixColumn2;
// See if was presolved
if (originalColumns) {
columnLower = lpSolver->columnLower();
columnUpper = lpSolver->columnUpper();
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
saveColumnLower[iColumn] = columnLower[iColumn];
saveColumnUpper[iColumn] = columnUpper[iColumn];
}
pinfo.postsolve(true);
columnLower = originalLpSolver->columnLower();
columnUpper = originalLpSolver->columnUpper();
double * newColumnLower = lpSolver->columnLower();
double * newColumnUpper = lpSolver->columnUpper();
for (iColumn = 0; iColumn < numberColumns; iColumn++) {
int jColumn = originalColumns[iColumn];
columnLower[jColumn] = CoinMax(columnLower[jColumn], newColumnLower[iColumn]);
columnUpper[jColumn] = CoinMin(columnUpper[jColumn], newColumnUpper[iColumn]);
}
numberColumns = originalLpSolver->numberColumns();
delete [] originalColumns;
}
delete [] saveColumnLower;
delete [] saveColumnUpper;
if (!originalColumns) {
// Basis
memcpy(originalLpSolver->statusArray(), lpSolver->statusArray(), numberRows + numberColumns);
memcpy(originalLpSolver->primalColumnSolution(), lpSolver->primalColumnSolution(), numberColumns*sizeof(double));
memcpy(originalLpSolver->primalRowSolution(), lpSolver->primalRowSolution(), numberRows*sizeof(double));
// Fix in solver
columnLower = lpSolver->columnLower();
columnUpper = lpSolver->columnUpper();
}
double * originalColumnLower = originalLpSolver->columnLower();
double * originalColumnUpper = originalLpSolver->columnUpper();
// number fixed
doAction = 0;
for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
originalColumnLower[iColumn] = columnLower[iColumn];
originalColumnUpper[iColumn] = columnUpper[iColumn];
if (columnLower[iColumn] == columnUpper[iColumn])
doAction++;
}
COIN_DETAIL_PRINT(printf("%d fixed by vub preprocessing\n", doAction));
if (originalColumns) {
originalLpSolver->initialSolve();
}
delete clpSolver;
return NULL;
}
int doHeuristics(CbcModel * model, int type, CbcOrClpParam* parameters_,
int numberParameters_,int noPrinting_,int initialPumpTune)
{
#ifdef JJF_ZERO //NEW_STYLE_SOLVER==0
CbcOrClpParam * parameters_ = parameters;
int numberParameters_ = numberParameters;
bool noPrinting_ = noPrinting_;
#endif
char generalPrint[10000];
CoinMessages generalMessages = model->messages();
CoinMessageHandler * generalMessageHandler = model->messageHandler();
//generalMessageHandler->setPrefix(false);
bool anyToDo = false;
int logLevel = parameters_[whichParam(CLP_PARAM_INT_LOGLEVEL, numberParameters_, parameters_)].intValue();
int useFpump = parameters_[whichParam(CBC_PARAM_STR_FPUMP, numberParameters_, parameters_)].currentOptionAsInteger();
int useRounding = parameters_[whichParam(CBC_PARAM_STR_ROUNDING, numberParameters_, parameters_)].currentOptionAsInteger();
int useGreedy = parameters_[whichParam(CBC_PARAM_STR_GREEDY, numberParameters_, parameters_)].currentOptionAsInteger();
int useCombine = parameters_[whichParam(CBC_PARAM_STR_COMBINE, numberParameters_, parameters_)].currentOptionAsInteger();
int useProximity = parameters_[whichParam(CBC_PARAM_STR_PROXIMITY, numberParameters_, parameters_)].currentOptionAsInteger();
int useCrossover = parameters_[whichParam(CBC_PARAM_STR_CROSSOVER2, numberParameters_, parameters_)].currentOptionAsInteger();
//int usePivotC = parameters_[whichParam(CBC_PARAM_STR_PIVOTANDCOMPLEMENT, numberParameters_, parameters_)].currentOptionAsInteger();
int usePivotF = parameters_[whichParam(CBC_PARAM_STR_PIVOTANDFIX, numberParameters_, parameters_)].currentOptionAsInteger();
int useRand = parameters_[whichParam(CBC_PARAM_STR_RANDROUND, numberParameters_, parameters_)].currentOptionAsInteger();
int useRINS = parameters_[whichParam(CBC_PARAM_STR_RINS, numberParameters_, parameters_)].currentOptionAsInteger();
int useRENS = parameters_[whichParam(CBC_PARAM_STR_RENS, numberParameters_, parameters_)].currentOptionAsInteger();
int useDINS = parameters_[whichParam(CBC_PARAM_STR_DINS, numberParameters_, parameters_)].currentOptionAsInteger();
int useDIVING2 = parameters_[whichParam(CBC_PARAM_STR_DIVINGS, numberParameters_, parameters_)].currentOptionAsInteger();
int useNaive = parameters_[whichParam(CBC_PARAM_STR_NAIVE, numberParameters_, parameters_)].currentOptionAsInteger();
int kType = (type < 10) ? type : 1;
assert (kType == 1 || kType == 2);
// FPump done first as it only works if no solution
if (useFpump >= kType && useFpump <= kType + 1) {
anyToDo = true;
CbcHeuristicFPump heuristic4(*model);
double dextra3 = parameters_[whichParam(CBC_PARAM_DBL_SMALLBAB, numberParameters_, parameters_)].doubleValue();
heuristic4.setFractionSmall(dextra3);
double dextra1 = parameters_[whichParam(CBC_PARAM_DBL_ARTIFICIALCOST, numberParameters_, parameters_)].doubleValue();
if (dextra1)
heuristic4.setArtificialCost(dextra1);
heuristic4.setMaximumPasses(parameters_[whichParam(CBC_PARAM_INT_FPUMPITS, numberParameters_, parameters_)].intValue());
if (parameters_[whichParam(CBC_PARAM_INT_FPUMPITS, numberParameters_, parameters_)].intValue() == 21)
heuristic4.setIterationRatio(1.0);
int pumpTune = parameters_[whichParam(CBC_PARAM_INT_FPUMPTUNE, numberParameters_, parameters_)].intValue();
int pumpTune2 = parameters_[whichParam(CBC_PARAM_INT_FPUMPTUNE2, numberParameters_, parameters_)].intValue();
if (pumpTune > 0) {
bool printStuff = (pumpTune != initialPumpTune || logLevel > 1 || pumpTune2 > 0)
&& !noPrinting_;
if (printStuff) {
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< "Options for feasibility pump - "
<< CoinMessageEol;
}
/*
>=10000000 for using obj
>=1000000 use as accumulate switch
>=1000 use index+1 as number of large loops
>=100 use dextra1 as cutoff
%100 == 10,20 etc for experimentation
1 == fix ints at bounds, 2 fix all integral ints, 3 and continuous at bounds
4 and static continuous, 5 as 3 but no internal integers
6 as 3 but all slack basis!
*/
double value = model->solver()->getObjSense() * model->solver()->getObjValue();
int w = pumpTune / 10;
int i = w % 10;
w /= 10;
int c = w % 10;
w /= 10;
int r = w;
int accumulate = r / 1000;
r -= 1000 * accumulate;
if (accumulate >= 10) {
int which = accumulate / 10;
accumulate -= 10 * which;
which--;
// weights and factors
double weight[] = {0.01, 0.01, 0.1, 0.1, 0.5, 0.5, 1.0, 1.0, 5.0, 5.0};
double factor[] = {0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5, 0.1, 0.5};
heuristic4.setInitialWeight(weight[which]);
heuristic4.setWeightFactor(factor[which]);
if (printStuff) {
sprintf(generalPrint, "Initial weight for objective %g, decay factor %g",
weight[which], factor[which]);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
}
// fake cutoff
if (c) {
double cutoff;
model->solver()->getDblParam(OsiDualObjectiveLimit, cutoff);
cutoff = CoinMin(cutoff, value + 0.05 * fabs(value) * c);
double fakeCutoff = parameters_[whichParam(CBC_PARAM_DBL_FAKECUTOFF, numberParameters_, parameters_)].doubleValue();
if (fakeCutoff)
cutoff = fakeCutoff;
heuristic4.setFakeCutoff(cutoff);
if (printStuff) {
sprintf(generalPrint, "Fake cutoff of %g", cutoff);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
}
int offRandomEtc = 0;
if (pumpTune2) {
if ((pumpTune2 / 1000) != 0) {
offRandomEtc = 1000000 * (pumpTune2 / 1000);
if (printStuff) {
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< "Feasibility pump may run twice"
<< CoinMessageEol;
}
pumpTune2 = pumpTune2 % 1000;
}
if ((pumpTune2 / 100) != 0) {
offRandomEtc += 100 * (pumpTune2 / 100);
if (printStuff) {
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< "Not using randomized objective"
<< CoinMessageEol;
}
}
int maxAllowed = pumpTune2 % 100;
if (maxAllowed) {
offRandomEtc += 1000 * maxAllowed;
if (printStuff) {
sprintf(generalPrint, "Fixing if same for %d passes",
maxAllowed);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
}
}
if (accumulate) {
heuristic4.setAccumulate(accumulate);
if (printStuff) {
if (accumulate) {
sprintf(generalPrint, "Accumulate of %d", accumulate);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
}
}
if (r) {
// also set increment
//double increment = (0.01*i+0.005)*(fabs(value)+1.0e-12);
double increment = 0.0;
double fakeIncrement = parameters_[whichParam(CBC_PARAM_DBL_FAKEINCREMENT, numberParameters_, parameters_)].doubleValue();
if (fakeIncrement)
increment = fakeIncrement;
heuristic4.setAbsoluteIncrement(increment);
heuristic4.setMaximumRetries(r + 1);
if (printStuff) {
if (increment) {
sprintf(generalPrint, "Increment of %g", increment);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
sprintf(generalPrint, "%d retries", r + 1);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
}
if (i + offRandomEtc) {
heuristic4.setFeasibilityPumpOptions(i*10 + offRandomEtc);
if (printStuff) {
sprintf(generalPrint, "Feasibility pump options of %d",
i*10 + offRandomEtc);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
}
pumpTune = pumpTune % 100;
if (pumpTune == 6)
pumpTune = 13;
heuristic4.setWhen((pumpTune % 10) + 10);
if (printStuff) {
sprintf(generalPrint, "Tuning (fixing) %d", pumpTune % 10);
generalMessageHandler->message(CBC_GENERAL, generalMessages)
<< generalPrint
<< CoinMessageEol;
}
}
heuristic4.setHeuristicName("feasibility pump");
//#define ROLF
#ifdef ROLF
CbcHeuristicFPump pump(*model);
pump.setMaximumTime(60);
pump.setMaximumPasses(100);
pump.setMaximumRetries(1);
pump.setFixOnReducedCosts(0);
pump.setHeuristicName("Feasibility pump");
pump.setFractionSmall(1.0);
pump.setWhen(13);
model->addHeuristic(&pump);
#else
model->addHeuristic(&heuristic4);
#endif
}
if (useRounding >= type && useRounding >= kType && useRounding <= kType + 1) {
CbcRounding heuristic1(*model);
heuristic1.setHeuristicName("rounding");
model->addHeuristic(&heuristic1) ;
anyToDo = true;
}
if (useGreedy >= type && useGreedy >= kType && useGreedy <= kType + 1) {
CbcHeuristicGreedyCover heuristic3(*model);
heuristic3.setHeuristicName("greedy cover");
CbcHeuristicGreedyEquality heuristic3a(*model);
heuristic3a.setHeuristicName("greedy equality");
model->addHeuristic(&heuristic3);
model->addHeuristic(&heuristic3a);
anyToDo = true;
}
if ((useRENS==7 && kType==1) || (useRENS==8 && kType==2)) {
useRENS=1+2*(useRENS-7);
CbcHeuristicRENS heuristic6a(*model);
heuristic6a.setHeuristicName("RENSdj");
heuristic6a.setFractionSmall(0.6/*3.4*/);
heuristic6a.setFeasibilityPumpOptions(3);
heuristic6a.setNumberNodes(10);
heuristic6a.setWhereFrom(4*256+4*1);
heuristic6a.setWhen(2);
heuristic6a.setRensType(1+16);
model->addHeuristic(&heuristic6a) ;
heuristic6a.setHeuristicName("RENSub");
heuristic6a.setFractionSmall(0.4);
heuristic6a.setFeasibilityPumpOptions(1008003);
heuristic6a.setNumberNodes(50);
heuristic6a.setWhereFrom(4*256+4*1);
heuristic6a.setWhen(2);
heuristic6a.setRensType(2+16);
model->addHeuristic(&heuristic6a) ;
}
if (useRENS >= kType && useRENS <= kType + 1) {
#ifndef JJF_ONE
CbcHeuristicRENS heuristic6(*model);
heuristic6.setHeuristicName("RENS");
heuristic6.setFractionSmall(0.4);
heuristic6.setFeasibilityPumpOptions(1008003);
int nodes [] = { -2, 50, 50, 50, 200, 1000, 10000};
heuristic6.setNumberNodes(nodes[useRENS]);
#else
CbcHeuristicVND heuristic6(*model);
heuristic6.setHeuristicName("VND");
heuristic5.setFractionSmall(0.5);
heuristic5.setDecayFactor(5.0);
#endif
model->addHeuristic(&heuristic6) ;
anyToDo = true;
}
if (useNaive >= kType && useNaive <= kType + 1) {
CbcHeuristicNaive heuristic5b(*model);
heuristic5b.setHeuristicName("Naive");
heuristic5b.setFractionSmall(0.4);
heuristic5b.setNumberNodes(50);
model->addHeuristic(&heuristic5b) ;
anyToDo = true;
}
int useDIVING = 0;
{
int useD;
useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGV, numberParameters_, parameters_)].currentOptionAsInteger();
useDIVING |= 1 * ((useD >= kType) ? 1 : 0);
useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGG, numberParameters_, parameters_)].currentOptionAsInteger();
useDIVING |= 2 * ((useD >= kType) ? 1 : 0);
useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGF, numberParameters_, parameters_)].currentOptionAsInteger();
useDIVING |= 4 * ((useD >= kType) ? 1 : 0);
useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGC, numberParameters_, parameters_)].currentOptionAsInteger();
useDIVING |= 8 * ((useD >= kType) ? 1 : 0);
useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGL, numberParameters_, parameters_)].currentOptionAsInteger();
useDIVING |= 16 * ((useD >= kType) ? 1 : 0);
useD = parameters_[whichParam(CBC_PARAM_STR_DIVINGP, numberParameters_, parameters_)].currentOptionAsInteger();
useDIVING |= 32 * ((useD >= kType) ? 1 : 0);
}
if (useDIVING2 >= kType && useDIVING2 <= kType + 1) {
int diveOptions = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();
if (diveOptions < 0 || diveOptions > 10)
diveOptions = 2;
CbcHeuristicJustOne heuristicJustOne(*model);
heuristicJustOne.setHeuristicName("DiveAny");
heuristicJustOne.setWhen(diveOptions);
// add in others
CbcHeuristicDiveCoefficient heuristicDC(*model);
heuristicDC.setHeuristicName("DiveCoefficient");
heuristicJustOne.addHeuristic(&heuristicDC, 1.0) ;
CbcHeuristicDiveFractional heuristicDF(*model);
heuristicDF.setHeuristicName("DiveFractional");
heuristicJustOne.addHeuristic(&heuristicDF, 1.0) ;
CbcHeuristicDiveGuided heuristicDG(*model);
heuristicDG.setHeuristicName("DiveGuided");
heuristicJustOne.addHeuristic(&heuristicDG, 1.0) ;
CbcHeuristicDiveLineSearch heuristicDL(*model);
heuristicDL.setHeuristicName("DiveLineSearch");
heuristicJustOne.addHeuristic(&heuristicDL, 1.0) ;
CbcHeuristicDivePseudoCost heuristicDP(*model);
heuristicDP.setHeuristicName("DivePseudoCost");
heuristicJustOne.addHeuristic(&heuristicDP, 1.0) ;
CbcHeuristicDiveVectorLength heuristicDV(*model);
heuristicDV.setHeuristicName("DiveVectorLength");
heuristicJustOne.addHeuristic(&heuristicDV, 1.0) ;
// Now normalize probabilities
heuristicJustOne.normalizeProbabilities();
model->addHeuristic(&heuristicJustOne) ;
}
if (useDIVING > 0) {
int majorIterations=64;
int diveOptions2 = parameters_[whichParam(CBC_PARAM_INT_DIVEOPT, numberParameters_, parameters_)].intValue();
int diveOptions;
if (diveOptions2 > 99) {
// switch on various active set stuff
diveOptions = diveOptions2%100;
diveOptions2 /= 100;
} else {
diveOptions = diveOptions2;
diveOptions2 = 0;
}
if (diveOptions < 0 || diveOptions > 9)
diveOptions = 2;
if ((useDIVING&1) != 0) {
CbcHeuristicDiveVectorLength heuristicDV(*model);
heuristicDV.setHeuristicName("DiveVectorLength");
heuristicDV.setWhen(diveOptions);
if (diveOptions2) {
heuristicDV.setMaxIterations(majorIterations);
heuristicDV.setPercentageToFix(0.0);
heuristicDV.setMaxSimplexIterations(COIN_INT_MAX);
heuristicDV.setMaxSimplexIterationsAtRoot(COIN_INT_MAX-(diveOptions2-1));
}
model->addHeuristic(&heuristicDV) ;
}
if ((useDIVING&2) != 0) {
CbcHeuristicDiveGuided heuristicDG(*model);
heuristicDG.setHeuristicName("DiveGuided");
heuristicDG.setWhen(diveOptions);
if (diveOptions2) {
heuristicDG.setMaxIterations(majorIterations);
heuristicDG.setPercentageToFix(0.0);
heuristicDG.setMaxSimplexIterations(COIN_INT_MAX);
heuristicDG.setMaxSimplexIterationsAtRoot(COIN_INT_MAX-(diveOptions2-1));
}
model->addHeuristic(&heuristicDG) ;
}
if ((useDIVING&4) != 0) {
CbcHeuristicDiveFractional heuristicDF(*model);
heuristicDF.setHeuristicName("DiveFractional");
heuristicDF.setWhen(diveOptions);
if (diveOptions2) {
heuristicDF.setMaxIterations(majorIterations);
heuristicDF.setPercentageToFix(0.0);
heuristicDF.setMaxSimplexIterations(COIN_INT_MAX);
heuristicDF.setMaxSimplexIterationsAtRoot(COIN_INT_MAX-(diveOptions2-1));
}
model->addHeuristic(&heuristicDF) ;
}
if ((useDIVING&8) != 0) {
CbcHeuristicDiveCoefficient heuristicDC(*model);
heuristicDC.setHeuristicName("DiveCoefficient");
heuristicDC.setWhen(diveOptions);
if (diveOptions2) {
heuristicDC.setMaxIterations(majorIterations);
heuristicDC.setPercentageToFix(0.0);
heuristicDC.setMaxSimplexIterations(COIN_INT_MAX);
heuristicDC.setMaxSimplexIterationsAtRoot(COIN_INT_MAX-(diveOptions2-1));
}
model->addHeuristic(&heuristicDC) ;
}
if ((useDIVING&16) != 0) {
CbcHeuristicDiveLineSearch heuristicDL(*model);
heuristicDL.setHeuristicName("DiveLineSearch");
heuristicDL.setWhen(diveOptions);
if (diveOptions2) {
heuristicDL.setMaxIterations(majorIterations);
heuristicDL.setPercentageToFix(0.0);
heuristicDL.setMaxSimplexIterations(COIN_INT_MAX);
heuristicDL.setMaxSimplexIterationsAtRoot(COIN_INT_MAX-(diveOptions2-1));
}
model->addHeuristic(&heuristicDL) ;
}
if ((useDIVING&32) != 0) {
CbcHeuristicDivePseudoCost heuristicDP(*model);
heuristicDP.setHeuristicName("DivePseudoCost");
heuristicDP.setWhen(diveOptions /*+ diveOptions2*/);
if (diveOptions2) {
heuristicDP.setMaxIterations(majorIterations);
heuristicDP.setPercentageToFix(0.0);
heuristicDP.setMaxSimplexIterations(COIN_INT_MAX);
heuristicDP.setMaxSimplexIterationsAtRoot(COIN_INT_MAX-(diveOptions2-1));
}
model->addHeuristic(&heuristicDP) ;
}
anyToDo = true;
}
#ifdef JJF_ZERO
if (usePivotC >= type && usePivotC <= kType + 1) {
CbcHeuristicPivotAndComplement heuristic(*model);
heuristic.setHeuristicName("pivot and complement");
heuristic.setFractionSmall(10.0); // normally 0.5
model->addHeuristic(&heuristic);
anyToDo = true;
}
#endif
if (usePivotF >= type && usePivotF <= kType + 1) {
CbcHeuristicPivotAndFix heuristic(*model);
heuristic.setHeuristicName("pivot and fix");
heuristic.setFractionSmall(10.0); // normally 0.5
model->addHeuristic(&heuristic);
anyToDo = true;
}
if (useRand >= type && useRand <= kType + 1) {
CbcHeuristicRandRound heuristic(*model);
heuristic.setHeuristicName("randomized rounding");
heuristic.setFractionSmall(10.0); // normally 0.5
model->addHeuristic(&heuristic);
anyToDo = true;
}
if (useDINS >= kType && useDINS <= kType + 1) {
CbcHeuristicDINS heuristic5a(*model);
heuristic5a.setHeuristicName("DINS");
heuristic5a.setFractionSmall(0.6);
if (useDINS < 4)
heuristic5a.setDecayFactor(5.0);
else
heuristic5a.setDecayFactor(1.5);
heuristic5a.setNumberNodes(1000);
model->addHeuristic(&heuristic5a) ;
anyToDo = true;
}
if (useRINS >= kType && useRINS <= kType + 1) {
CbcHeuristicRINS heuristic5(*model);
heuristic5.setHeuristicName("RINS");
if (useRINS < 4) {
heuristic5.setFractionSmall(0.5);
heuristic5.setDecayFactor(5.0);
} else {
heuristic5.setFractionSmall(0.6);
heuristic5.setDecayFactor(1.5);
}
model->addHeuristic(&heuristic5) ;
anyToDo = true;
}
if (useCombine >= kType && useCombine <= kType + 1) {
CbcHeuristicLocal heuristic2(*model);
heuristic2.setHeuristicName("combine solutions");
heuristic2.setFractionSmall(0.5);
heuristic2.setSearchType(1);
model->addHeuristic(&heuristic2);
anyToDo = true;
}
if ((useProximity >= kType && useProximity <= kType + 1)||
(kType == 1 && useProximity >= 4)) {
CbcHeuristicProximity heuristic2a(*model);
heuristic2a.setHeuristicName("Proximity Search");
heuristic2a.setFractionSmall(9999999.0);
heuristic2a.setNumberNodes(30);
heuristic2a.setFeasibilityPumpOptions(-2);
if (useProximity>=4) {
const int nodes[]={10,100,300};
heuristic2a.setNumberNodes(nodes[useProximity-4]);
// more print out and stronger feasibility pump
if (useProximity==6)
heuristic2a.setFeasibilityPumpOptions(-3);
}
model->addHeuristic(&heuristic2a);
anyToDo = true;
}
if (useCrossover >= kType && useCrossover <= kType + 1) {
CbcHeuristicCrossover heuristic2a(*model);
heuristic2a.setHeuristicName("crossover");
heuristic2a.setFractionSmall(0.3);
// just fix at lower
heuristic2a.setWhen(11);
model->addHeuristic(&heuristic2a);
model->setMaximumSavedSolutions(5);
anyToDo = true;
}
int heurSwitches = parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() % 100;
if (heurSwitches) {
for (int iHeur = 0; iHeur < model->numberHeuristics(); iHeur++) {
CbcHeuristic * heuristic = model->heuristic(iHeur);
heuristic->setSwitches(heurSwitches);
}
}
if (type == 2 && anyToDo) {
// Do heuristics
#ifndef JJF_ONE
// clean copy
CbcModel model2(*model);
// But get rid of heuristics in model
model->doHeuristicsAtRoot(2);
if (logLevel <= 1)
model2.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
OsiBabSolver defaultC;
//solver_->setAuxiliaryInfo(&defaultC);
model2.passInSolverCharacteristics(&defaultC);
// Save bounds
int numberColumns = model2.solver()->getNumCols();
model2.createContinuousSolver();
bool cleanModel = !model2.numberIntegers() && !model2.numberObjects();
model2.findIntegers(false);
int heurOptions = (parameters_[whichParam(CBC_PARAM_INT_HOPTIONS, numberParameters_, parameters_)].intValue() / 100) % 100;
if (heurOptions == 0 || heurOptions == 2) {
model2.doHeuristicsAtRoot(1);
} else if (heurOptions == 1 || heurOptions == 3) {
model2.setMaximumNodes(-1);
CbcStrategyDefault strategy(0, 5, 5);
strategy.setupPreProcessing(1, 0);
model2.setStrategy(strategy);
model2.branchAndBound();
}
if (cleanModel)
model2.zapIntegerInformation(false);
if (model2.bestSolution()) {
double value = model2.getMinimizationObjValue();
model->setCutoff(value);
model->setBestSolution(model2.bestSolution(), numberColumns, value);
model->setSolutionCount(1);
model->setNumberHeuristicSolutions(1);
}
#else
if (logLevel <= 1)
model->solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
OsiBabSolver defaultC;
//solver_->setAuxiliaryInfo(&defaultC);
model->passInSolverCharacteristics(&defaultC);
// Save bounds
int numberColumns = model->solver()->getNumCols();
model->createContinuousSolver();
bool cleanModel = !model->numberIntegers() && !model->numberObjects();
model->findIntegers(false);
model->doHeuristicsAtRoot(1);
if (cleanModel)
model->zapIntegerInformation(false);
#endif
return 0;
} else {
return 0;
}
}