limp-cbc-0.3.2.0: cbits/coin/CglMixedIntegerRounding2.cpp
// LAST EDIT:
//-----------------------------------------------------------------------------
// name: Mixed Integer Rounding Cut Generator
// authors: Joao Goncalves (jog7@lehigh.edu)
// Laszlo Ladanyi (ladanyi@us.ibm.com)
// date: August 11, 2004
//-----------------------------------------------------------------------------
// Copyright (C) 2004, International Business Machines Corporation and others.
// All Rights Reserved.
// This code is published under the Eclipse Public License.
//#include <cmath>
//#include <cstdlib>
#include <cassert>
#include "CoinPragma.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinPackedVector.hpp"
#include "CglMixedIntegerRounding2.hpp"
//-----------------------------------------------------------------------------
// Generate Mixed Integer Rounding inequality
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::generateCuts(const OsiSolverInterface& si,
OsiCuts& cs,
const CglTreeInfo info)
{
// If the LP or integer presolve is used, then need to redo preprocessing
// everytime this function is called. Otherwise, just do once.
bool preInit = false;
bool preReso = false;
si.getHintParam(OsiDoPresolveInInitial, preInit);
si.getHintParam(OsiDoPresolveInResolve, preReso);
if (preInit == false && preReso == false && doPreproc_ == -1 ) { // Do once
if (doneInitPre_ == false) {
mixIntRoundPreprocess(si);
doneInitPre_ = true;
}
}
else {
if(doPreproc_ == 1){ // Do everytime
mixIntRoundPreprocess(si);
doneInitPre_ = true;
}
else {
if (doneInitPre_ == false) {
mixIntRoundPreprocess(si);
doneInitPre_ = true;
}
}
}
int numberRowCutsBefore = cs.sizeRowCuts();
const double* xlp = si.getColSolution(); // LP solution
const double* colUpperBound = si.getColUpper(); // vector of upper bounds
const double* colLowerBound = si.getColLower(); // vector of lower bounds
// get matrix by row
const CoinPackedMatrix & tempMatrixByRow = *si.getMatrixByRow();
CoinPackedMatrix matrixByRow(false,0.0,0.0);
// There are no duplicates but this is faster
matrixByRow.submatrixOfWithDuplicates(tempMatrixByRow, numRows_, indRows_);
CoinPackedMatrix matrixByCol(matrixByRow,0,0,true);
//matrixByCol.reverseOrdering();
//const CoinPackedMatrix & matrixByRow = *si.getMatrixByRow();
const double* LHS = si.getRowActivity();
//const double* coefByRow = matrixByRow.getElements();
//const int* colInds = matrixByRow.getIndices();
//const int* rowStarts = matrixByRow.getVectorStarts();
// get matrix by column
//const CoinPackedMatrix & matrixByCol = *si.getMatrixByCol();
const double* coefByCol = matrixByCol.getElements();
const int* rowInds = matrixByCol.getIndices();
const int* colStarts = matrixByCol.getVectorStarts();
generateMirCuts(si, xlp, colUpperBound, colLowerBound,
matrixByRow, LHS, //coefByRow,
//colInds, rowStarts, //matrixByCol,
coefByCol, rowInds, colStarts,
cs);
if (!info.inTree&&((info.options&4)==4||((info.options&8)&&!info.pass))) {
int numberRowCutsAfter = cs.sizeRowCuts();
for (int i=numberRowCutsBefore;i<numberRowCutsAfter;i++)
cs.rowCutPtr(i)->setGloballyValid();
}
}
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
CglMixedIntegerRounding2::CglMixedIntegerRounding2 ()
:
CglCutGenerator()
{
gutsOfConstruct(1, true, 1, -1);
}
//-------------------------------------------------------------------
// Alternate Constructor
//-------------------------------------------------------------------
CglMixedIntegerRounding2::CglMixedIntegerRounding2 (const int maxaggr,
const bool multiply,
const int criterion,
const int preproc)
:
CglCutGenerator()
{
gutsOfConstruct(maxaggr, multiply, criterion, preproc);
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
CglMixedIntegerRounding2::CglMixedIntegerRounding2 (
const CglMixedIntegerRounding2 & rhs)
:
CglCutGenerator(rhs)
{
gutsOfCopy(rhs);
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
CglCutGenerator *
CglMixedIntegerRounding2::clone() const
{
return new CglMixedIntegerRounding2(*this);
}
//------------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
CglMixedIntegerRounding2 &
CglMixedIntegerRounding2::operator=(const CglMixedIntegerRounding2& rhs)
{
if (this != &rhs) {
gutsOfDelete();
CglCutGenerator::operator=(rhs);
gutsOfCopy(rhs);
}
return *this;
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
CglMixedIntegerRounding2::~CglMixedIntegerRounding2 ()
{
gutsOfDelete();
}
//-------------------------------------------------------------------
// Construct
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::gutsOfConstruct (const int maxaggr,
const bool multiply,
const int criterion,
const int preproc)
{
if (maxaggr > 0) {
MAXAGGR_ = maxaggr;
}
else {
throw CoinError("Unallowable value. maxaggr must be > 0",
"gutsOfConstruct","CglMixedIntegerRounding2");
}
MULTIPLY_ = multiply;
if ((criterion >= 1) && (criterion <= 3)) {
CRITERION_ = criterion;
}
else {
throw CoinError("Unallowable value. criterion must be 1, 2 or 3",
"gutsOfConstruct","CglMixedIntegerRounding2");
}
if ((preproc >= -1) && (preproc <= 2)) {
doPreproc_ = preproc;
}
else {
throw CoinError("Unallowable value. preproc must be -1, 0 or 1",
"gutsOfConstruct","CglMixedIntegerRounding");
}
EPSILON_ = 1.0e-6;
UNDEFINED_ = -1;
TOLERANCE_ = 1.0e-4;
numRows_ = 0;
numCols_ = 0;
doneInitPre_ = false;
vubs_ = 0;
vlbs_ = 0;
rowTypes_ = 0;
indRows_ = 0;
numRowMix_ = 0;
indRowMix_ = 0;
numRowCont_ = 0;
indRowCont_ = 0;
numRowInt_ = 0;
indRowInt_ = 0;
numRowContVB_ = 0;
indRowContVB_ = 0;
integerType_ = NULL;
sense_=NULL;
RHS_=NULL;
}
//-------------------------------------------------------------------
// Delete
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::gutsOfDelete ()
{
if (vubs_ != 0) { delete [] vubs_; vubs_ = 0; }
if (vlbs_ != 0) { delete [] vlbs_; vlbs_ = 0; }
if (rowTypes_ != 0) { delete [] rowTypes_; rowTypes_ = 0; }
if (indRows_ != 0) { delete [] indRows_; indRows_ = 0; }
if (indRowMix_ != 0) { delete [] indRowMix_; indRowMix_ = 0; }
if (indRowCont_ != 0) { delete [] indRowCont_; indRowCont_ = 0; }
if (indRowInt_ != 0) { delete [] indRowInt_; indRowInt_ = 0; }
if (indRowContVB_ != 0) { delete [] indRowContVB_; indRowContVB_ = 0; }
if (integerType_ !=NULL) { delete [] integerType_; integerType_=NULL;}
if (sense_ !=NULL) { delete [] sense_; sense_=NULL;}
if (RHS_ !=NULL) { delete [] RHS_; RHS_=NULL;}
}
//-------------------------------------------------------------------
// Copy
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::gutsOfCopy (const CglMixedIntegerRounding2& rhs)
{
MAXAGGR_ = rhs.MAXAGGR_;
MULTIPLY_ = rhs.MULTIPLY_;
CRITERION_ = rhs.CRITERION_;
EPSILON_ = rhs.EPSILON_;
UNDEFINED_ = rhs.UNDEFINED_;
TOLERANCE_ = rhs.TOLERANCE_;
doPreproc_ = rhs.doPreproc_;
numRows_ = rhs.numRows_;
numCols_ = rhs.numCols_;
doneInitPre_ = rhs.doneInitPre_;
numRowMix_ = rhs.numRowMix_;
numRowCont_ = rhs.numRowCont_;
numRowInt_ = rhs.numRowInt_;
numRowContVB_ = rhs.numRowContVB_;
if (numCols_ > 0) {
vubs_ = new CglMixIntRoundVUB2 [numCols_];
vlbs_ = new CglMixIntRoundVLB2 [numCols_];
CoinDisjointCopyN(rhs.vubs_, numCols_, vubs_);
CoinDisjointCopyN(rhs.vlbs_, numCols_, vlbs_);
integerType_ = CoinCopyOfArray(rhs.integerType_,numCols_);
}
else {
vubs_ = 0;
vlbs_ = 0;
integerType_ = NULL;
}
if (numRows_ > 0) {
rowTypes_ = new RowType [numRows_];
CoinDisjointCopyN(rhs.rowTypes_, numRows_, rowTypes_);
indRows_ = new int [numRows_];
CoinDisjointCopyN(rhs.indRows_, numRows_, indRows_);
sense_ = CoinCopyOfArray(rhs.sense_,numRows_);
RHS_ = CoinCopyOfArray(rhs.RHS_,numRows_);
}
else {
rowTypes_ = 0;
indRows_ = 0;
sense_=NULL;
RHS_=NULL;
}
if (numRowMix_ > 0) {
indRowMix_ = new int [numRowMix_];
CoinDisjointCopyN(rhs.indRowMix_, numRowMix_, indRowMix_);
}
else {
indRowMix_ = 0;
}
if (numRowCont_ > 0) {
indRowCont_ = new int [numRowCont_];
CoinDisjointCopyN(rhs.indRowCont_, numRowCont_, indRowCont_);
indRowContVB_ = new int [numRowCont_];
CoinDisjointCopyN(rhs.indRowContVB_, numRowCont_, indRowContVB_);
}
else {
indRowCont_ = 0;
indRowContVB_ = 0;
}
if (numRowInt_ > 0) {
indRowInt_ = new int [numRowInt_];
CoinDisjointCopyN(rhs.indRowInt_, numRowInt_, indRowInt_);
}
else {
indRowInt_ = 0;
}
}
//-------------------------------------------------------------------
// Do preprocessing
// It determines the type of each row. It also identifies the variable
// upper bounds and variable lower bounds.
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::
mixIntRoundPreprocess(const OsiSolverInterface& si)
{
// get matrix stored by row
const CoinPackedMatrix & matrixByRow = *si.getMatrixByRow();
numRows_ = si.getNumRows();
numCols_ = si.getNumCols();
const double* coefByRow = matrixByRow.getElements();
const int* colInds = matrixByRow.getIndices();
const int* rowStarts = matrixByRow.getVectorStarts();
const int* rowLengths = matrixByRow.getVectorLengths();
// Get copies of sense and RHS so we can modify if ranges
if (sense_) {
delete [] sense_;
delete [] RHS_;
}
sense_ = CoinCopyOfArray(si.getRowSense(),numRows_);
RHS_ = CoinCopyOfArray(si.getRightHandSide(),numRows_);
// Save integer type for speed
if (integerType_)
delete [] integerType_;
integerType_ = new char [numCols_];
int iColumn;
for (iColumn=0;iColumn<numCols_;iColumn++) {
if (si.isInteger(iColumn))
integerType_[iColumn]=1;
else
integerType_[iColumn]=0;
}
if (rowTypes_ != 0) {
delete [] rowTypes_; rowTypes_ = 0;
}
rowTypes_ = new RowType [numRows_]; // Destructor will free memory
// Summarize the row type infomation.
int numUNDEFINED = 0;
int numVARUB = 0;
int numVARLB = 0;
int numVAREQ = 0;
int numMIX = 0;
int numCONT = 0;
int numINT = 0;
int numOTHER = 0;
int iRow;
const double* rowActivity = si.getRowActivity();
const double* rowLower = si.getRowLower();
const double* rowUpper = si.getRowUpper();
for (iRow = 0; iRow < numRows_; ++iRow) {
// If range then choose which to use
if (sense_[iRow]=='R') {
if (rowActivity[iRow]-rowLower[iRow]<
rowUpper[iRow]-rowActivity[iRow]) {
// treat as G row
RHS_[iRow]=rowLower[iRow];
sense_[iRow]='G';
} else {
// treat as L row
RHS_[iRow]=rowUpper[iRow];
sense_[iRow]='L';
}
}
// get the type of a row
const RowType rowType =
determineRowType(/*si,*/ rowLengths[iRow], colInds+rowStarts[iRow],
coefByRow+rowStarts[iRow], sense_[iRow], RHS_[iRow]);
// store the type of the current row
rowTypes_[iRow] = rowType;
// Summarize information about row types
switch(rowType) {
case ROW_UNDEFINED:
++numUNDEFINED;
break;
case ROW_VARUB:
++numVARUB;
break;
case ROW_VARLB:
++numVARLB;
break;
case ROW_VAREQ:
++numVAREQ;
break;
case ROW_MIX:
++numMIX;
break;
case ROW_CONT:
++numCONT;
break;
case ROW_INT:
++numINT;
break;
case ROW_OTHER:
++numOTHER;
break;
default:
throw CoinError("Unknown row type", "MixIntRoundPreprocess",
"CglMixedIntegerRounding2");
}
}
// allocate memory for vector of indices of all rows
if (indRows_ != 0) { delete [] indRows_; indRows_ = 0; }
if (numRows_ > 0)
indRows_ = new int [numRows_]; // Destructor will free memory
// allocate memory for vector of indices of rows of type ROW_MIX
numRowMix_ = numMIX;
if (indRowMix_ != 0) { delete [] indRowMix_; indRowMix_ = 0; }
if (numRowMix_ > 0)
indRowMix_ = new int [numRowMix_]; // Destructor will free memory
// allocate memory for vector of indices of rows of type ROW_CONT
numRowCont_ = numCONT;
if (indRowCont_ != 0) { delete [] indRowCont_; indRowCont_ = 0; }
if (numRowCont_ > 0)
indRowCont_ = new int [numRowCont_]; // Destructor will free memory
// allocate memory for vector of indices of rows of type ROW_INT
numRowInt_ = numINT;
if (indRowInt_ != 0) { delete [] indRowInt_; indRowInt_ = 0; }
if (numRowInt_ > 0)
indRowInt_ = new int [numRowInt_]; // Destructor will free memory
#if CGL_DEBUG
std::cout << "The num of rows = " << numRows_ << std::endl;
std::cout << "Summary of Row Type" << std::endl;
std::cout << "numUNDEFINED = " << numUNDEFINED << std::endl;
std::cout << "numVARUB = " << numVARUB << std::endl;
std::cout << "numVARLB = " << numVARLB << std::endl;
std::cout << "numVAREQ = " << numVAREQ << std::endl;
std::cout << "numMIX = " << numMIX << std::endl;
std::cout << "numCONT = " << numCONT << std::endl;
std::cout << "numINT = " << numINT << std::endl;
std::cout << "numOTHER = " << numOTHER << std::endl;
#endif
//---------------------------------------------------------------------------
// Setup vubs_ and vlbs_
if (vubs_ != 0) { delete [] vubs_; vubs_ = 0; }
vubs_ = new CglMixIntRoundVUB2 [numCols_]; // Destructor will free
if (vlbs_ != 0) { delete [] vlbs_; vlbs_ = 0; }
vlbs_ = new CglMixIntRoundVLB2 [numCols_]; // Destructor will free
// Initialization. Altough this has been done in constructor, it is needed
// for the case where the mixIntRoundPreprocess is called more than once
for (int iCol = 0; iCol < numCols_; ++iCol) {
vubs_[iCol].setVar(UNDEFINED_);
vlbs_[iCol].setVar(UNDEFINED_);
}
int countM = 0;
int countC = 0;
int countI = 0;
for ( iRow = 0; iRow < numRows_; ++iRow) {
RowType rowType = rowTypes_[iRow];
// fill the vector indRows_ with the indices of all rows
indRows_[iRow] = iRow;
// fill the vector indRowMix_ with the indices of the rows of type ROW_MIX
if (rowType == ROW_MIX) {
indRowMix_[countM] = iRow;
countM++;
}
// fill the vector indRowCont_ with the indices of rows of type ROW_CONT
else if (rowType == ROW_CONT) {
indRowCont_[countC] = iRow;
countC++;
}
// fill the vector indRowInt_ with the indices of the rows of type ROW_INT
else if (rowType == ROW_INT) {
indRowInt_[countI] = iRow;
countI++;
}
// create vectors with variable lower and upper bounds
else if ( (rowType == ROW_VARUB) ||
(rowType == ROW_VARLB) ||
(rowType == ROW_VAREQ) ) {
int startPos = rowStarts[iRow];
int stopPos = startPos + rowLengths[iRow];
int xInd = 0, yInd = 0; // x is continuous, y is integer
double xCoef = 0.0, yCoef = 0.0;
for (int i = startPos; i < stopPos; ++i) {
if ( fabs(coefByRow[i]) > EPSILON_ ) {
if( integerType_[colInds[i]] ) {
yInd = colInds[i];
yCoef = coefByRow[i];
}
else {
xInd = colInds[i];
xCoef = coefByRow[i];
}
}
}
switch (rowType) {
case ROW_VARUB: // Inequality: x <= ? * y
vubs_[xInd].setVar(yInd);
vubs_[xInd].setVal(-yCoef / xCoef);
break;
case ROW_VARLB: // Inequality: x >= ? * y
vlbs_[xInd].setVar(yInd);
vlbs_[xInd].setVal(-yCoef / xCoef);
break;
case ROW_VAREQ: // Inequality: x >= AND <= ? * y
vubs_[xInd].setVar(yInd);
vubs_[xInd].setVal(-yCoef / xCoef);
vlbs_[xInd].setVar(yInd);
vlbs_[xInd].setVal(-yCoef / xCoef);
break;
default:
// I am getting compiler bug which gets here - I am disabling - JJF
//throw CoinError("Unknown row type: impossible",
// "MixIntRoundPreprocess",
// "CglMixedIntegerRounding2");
break;
}
}
}
// allocate memory for vector of indices of rows of type ROW_CONT
// that have at least one variable with variable upper or lower bound
if (indRowContVB_ != 0) { delete [] indRowContVB_; indRowContVB_ = 0; }
if (numRowCont_ > 0)
indRowContVB_ = new int [numRowCont_]; // Destructor will free memory
// create vector with rows of type ROW_CONT that have at least
// one variable with variable upper or lower bound
countC = 0;
for (int i = 0; i < numRowCont_; ++i) {
int indRow = indRowCont_[i];
int jStart = rowStarts[indRow];
int jStop = jStart + rowLengths[indRow];
for (int j = jStart; j < jStop; ++j) {
int indCol = colInds[j];
CglMixIntRoundVLB2 VLB = vlbs_[indCol];
CglMixIntRoundVUB2 VUB = vubs_[indCol];
if (( VLB.getVar() != UNDEFINED_ ) || ( VUB.getVar() != UNDEFINED_ ) ){
indRowContVB_[countC] = indRow;
countC++;
break;
}
}
}
numRowContVB_ = countC;
}
//-------------------------------------------------------------------
// Determine the type of a given row
//-------------------------------------------------------------------
CglMixedIntegerRounding2::RowType
CglMixedIntegerRounding2::determineRowType(//const OsiSolverInterface& si,
const int rowLen, const int* ind,
const double* coef, const char sense,
const double rhs) const
{
if (rowLen == 0 || fabs(rhs) > 1.0e20)
return ROW_UNDEFINED;
RowType rowType = ROW_UNDEFINED;
int numPosInt = 0; // num of positive integer variables
int numNegInt = 0; // num of negative integer variables
int numInt = 0; // num of integer variables
int numPosCon = 0; // num of positive continuous variables
int numNegCon = 0; // num of negative continuous variables
int numCon = 0; // num of continuous variables
// Summarize the variable types of the given row.
for ( int i = 0; i < rowLen; ++i ) {
if ( coef[i] < -EPSILON_ ) {
if( integerType_[ind[i]] )
++numNegInt;
else
++numNegCon;
}
else if ( coef[i] > EPSILON_ ) {
if( integerType_[ind[i]] )
++numPosInt;
else
++numPosCon;
}
}
numInt = numNegInt + numPosInt;
numCon = numNegCon + numPosCon;
#if CGL_DEBUG
std::cout << "numNegInt = " << numNegInt << std::endl;
std::cout << "numPosInt = " << numPosInt << std::endl;
std::cout << "numInt = " << numInt << std::endl;
std::cout << "numNegCon = " << numNegCon << std::endl;
std::cout << "numPosCon = " << numPosCon << std::endl;
std::cout << "numCon = " << numCon << std::endl;
std::cout << "rowLen = " << rowLen << std::endl;
#endif
//-------------------------------------------------------------------------
// Classify row type based on the types of variables.
if ((numInt > 0) && (numCon > 0)) {
if ((numInt == 1) && (numCon == 1) && (fabs(rhs) <= EPSILON_)) {
// It's a variable bound constraint
switch (sense) {
case 'L':
rowType = numPosCon == 1 ? ROW_VARUB : ROW_VARLB;
break;
case 'G':
rowType = numPosCon == 1 ? ROW_VARLB : ROW_VARUB;
break;
case 'E':
rowType = ROW_VAREQ;
break;
default:
break;
}
}
else {
// It's a constraint with continuous and integer variables;
// The total number of variables is at least 2
rowType = ROW_MIX;
}
}
else if (numInt == 0) {
// It's a constraint with only continuous variables
rowType = ROW_CONT;
}
else if ((numCon == 0) && ((sense == 'L') || (sense == 'G'))) {
// It's a <= or >= constraint with only integer variables
rowType = ROW_INT;
}
else
// It's a constraint that does not fit the above categories
rowType = ROW_OTHER;
return rowType;
}
//-------------------------------------------------------------------
// Generate MIR cuts
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::generateMirCuts(
const OsiSolverInterface& si,
const double* xlp,
const double* colUpperBound,
const double* colLowerBound,
const CoinPackedMatrix& matrixByRow,
const double* LHS,
//const double* coefByRow,
//const int* colInds,
//const int* rowStarts,
//const CoinPackedMatrix& matrixByCol,
const double* coefByCol,
const int* rowInds,
const int* colStarts,
OsiCuts& cs ) const
{
#if CGL_DEBUG
// OPEN FILE
std::ofstream fout("stats.dat");
#endif
// Define upper limit for the loop where the cMIRs are constructed
int upperLimit;
if (MULTIPLY_)
upperLimit = 2;
else
upperLimit = 1;
// create a vector with the columns that were used in the aggregation
int* listColsSelected = new int[MAXAGGR_];
// create a vector with the rows that were aggregated
int* listRowsAggregated = new int[MAXAGGR_];
// create a vector with the LP solutions of the slack variables
double* xlpExtra = new double[MAXAGGR_];
// loop until maximum number of aggregated rows is reached or a
// violated cut is found
int numRowMixAndRowContVB = numRowMix_ + numRowContVB_;
int numRowMixAndRowContVBAndRowInt = numRowMixAndRowContVB + numRowInt_;
// Get large enough vector
CoinIndexedVector rowAggregated(si.getNumCols());
CoinIndexedVector rowToAggregate(si.getNumCols());
CoinIndexedVector mixedKnapsack(si.getNumCols());
CoinIndexedVector contVariablesInS(si.getNumCols());
CoinIndexedVector rowToUse(si.getNumCols());
// And work vectors
CoinIndexedVector workVectors[4];
for (int i=0; i<4; i++)
workVectors[i].reserve(si.getNumCols());
CoinIndexedVector setRowsAggregated(si.getNumRows());
for (int iRow = 0; iRow < numRowMixAndRowContVBAndRowInt; ++iRow) {
int rowSelected; // row selected to be aggregated next
int colSelected; // column selected for pivot in aggregation
rowAggregated.clear();
double rhsAggregated;
// create a set with the indices of rows selected
setRowsAggregated.clear();
// loop until the maximum number of aggregated rows is reached
for (int iAggregate = 0; iAggregate < MAXAGGR_; ++iAggregate) {
if (iAggregate == 0) {
// select row
if (iRow < numRowMix_) {
rowSelected = indRowMix_[iRow];
}
else if (iRow < numRowMixAndRowContVB) {
rowSelected = indRowContVB_[iRow - numRowMix_];
}
else {
rowSelected = indRowInt_[iRow - numRowMixAndRowContVB];
}
copyRowSelected(iAggregate, rowSelected, setRowsAggregated,
listRowsAggregated, xlpExtra, sense_[rowSelected],
RHS_[rowSelected], LHS[rowSelected],
matrixByRow, rowAggregated, rhsAggregated);
}
else {
// search for a row to aggregate
bool foundRowToAggregate = selectRowToAggregate(
/*si,*/ rowAggregated,
colUpperBound, colLowerBound,
setRowsAggregated, xlp,
coefByCol, rowInds, colStarts,
rowSelected, colSelected);
// if finds row to aggregate, compute aggregated row
if (foundRowToAggregate) {
rowToAggregate.clear();
double rhsToAggregate;
listColsSelected[iAggregate] = colSelected;
copyRowSelected(iAggregate, rowSelected, setRowsAggregated,
listRowsAggregated, xlpExtra, sense_[rowSelected],
RHS_[rowSelected], LHS[rowSelected],
matrixByRow, rowToAggregate, rhsToAggregate);
// call aggregate row heuristic
aggregateRow(colSelected, rowToAggregate, rhsToAggregate,
rowAggregated, rhsAggregated);
}
else
break;
}
// construct cMIR with current rowAggregated
// and, if upperLimit=2 construct also a cMIR with
// the current rowAggregated multiplied by -1
for (int i = 0; i < upperLimit; ++i) {
// create vector for mixed knapsack constraint
double rhsMixedKnapsack;
rowToUse.copy(rowAggregated);
if (i==0) {
rhsMixedKnapsack = rhsAggregated;
} else {
rowToUse *= -1.0;
rhsMixedKnapsack = - rhsAggregated;
}
mixedKnapsack.clear();
double sStar = 0.0;
// create vector for the continuous variables in s
contVariablesInS.clear();
// call bound substitution heuristic
bool foundMixedKnapsack = boundSubstitution(
si, rowToUse,
xlp, xlpExtra,
colUpperBound, colLowerBound,
mixedKnapsack, rhsMixedKnapsack,
sStar, contVariablesInS);
// if it did not find a mixed knapsack it is because there is at
// least one integer variable with lower bound different than zero
// or there are no integer or continuous variables.
// In this case, we continue without trying to generate a c-MIR
if (!foundMixedKnapsack) {
#if CGL_DEBUG
std::cout << "couldn't create mixed knapsack" << std::endl;
#endif
continue;
}
OsiRowCut cMirCut;
// Find a c-MIR cut with the current mixed knapsack constraint
bool hasCut = cMirSeparation(si, matrixByRow, rowToUse,
listRowsAggregated, sense_, RHS_,
//coefByRow, colInds, rowStarts,
xlp, sStar, colUpperBound, colLowerBound,
mixedKnapsack,
rhsMixedKnapsack, contVariablesInS,
workVectors,cMirCut);
#if CGL_DEBUG
// PRINT STATISTICS
printStats(fout, hasCut, si, rowAggregated, rhsAggregated, xlp,
xlpExtra, listRowsAggregated, listColsSelected,
iAggregate+1, colUpperBound, colLowerBound );
#endif
// if a cut was found, insert it into cs
if (hasCut) {
#if CGL_DEBUG
std::cout << "MIR cut generated " << std::endl;
#endif
cs.insert(cMirCut);
}
}
}
}
// free memory
delete [] listColsSelected; listColsSelected = 0;
delete [] listRowsAggregated; listRowsAggregated = 0;
delete [] xlpExtra; xlpExtra = 0;
#if CGL_DEBUG
// CLOSE FILE
fout.close();
#endif
return;
}
//-------------------------------------------------------------------
// Copy row selected to CoinIndexedVector
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::copyRowSelected(
const int iAggregate,
const int rowSelected,
CoinIndexedVector& setRowsAggregated,
int* listRowsAggregated,
double* xlpExtra,
const char sen,
const double rhs,
const double lhs,
const CoinPackedMatrix& matrixByRow,
CoinIndexedVector& rowToAggregate,
double& rhsToAggregate) const
{
// copy the row selected to a vector of type CoinIndexedVector
const CoinShallowPackedVector reqdBySunCC = matrixByRow.getVector(rowSelected) ;
rowToAggregate = reqdBySunCC ;
rhsToAggregate = rhs;
// update list of indices of rows selected
setRowsAggregated.insert(rowSelected,1.0);
listRowsAggregated[iAggregate] = rowSelected;
// Add a slack variable if needed and compute its current value
if (sen == 'L') {
rowToAggregate.insert(numCols_ + iAggregate, 1);
xlpExtra[iAggregate] = rhs - lhs;
}
else if (sen == 'G') {
rowToAggregate.insert(numCols_ + iAggregate, -1);
xlpExtra[iAggregate] = lhs - rhs;
}
}
//-------------------------------------------------------------------
// Construct the set P* and select a row to aggregate
//-------------------------------------------------------------------
bool
CglMixedIntegerRounding2::selectRowToAggregate(
//const OsiSolverInterface& si,
const CoinIndexedVector& rowAggregated,
const double* colUpperBound,
const double* colLowerBound,
const CoinIndexedVector& setRowsAggregated,
const double* xlp, const double* coefByCol,
const int* rowInds, const int* colStarts,
int& rowSelected,
int& colSelected ) const
{
bool foundRowToAggregate = false;
double deltaMax = 0.0; // maximum delta
const int numColsAggregated = rowAggregated.getNumElements();
const int *rowAggregatedIndices = rowAggregated.getIndices();
const double *rowAggregatedElements = rowAggregated.denseVector();
for (int j = 0; j < numColsAggregated; ++j) {
// store the index and coefficient of column j
int indCol = rowAggregatedIndices[j];
if (indCol >= numCols_) continue;
double coefCol = rowAggregatedElements[indCol];
// Consider only continuous variables
if ( (integerType_[indCol]) || (fabs(coefCol) < EPSILON_)) continue;
// Compute current lower bound
CglMixIntRoundVLB2 VLB = vlbs_[indCol];
double LB = ( VLB.getVar() != UNDEFINED_ ) ?
VLB.getVal() * xlp[VLB.getVar()] : colLowerBound[indCol];
// Compute current upper bound
CglMixIntRoundVUB2 VUB = vubs_[indCol];
double UB = ( VUB.getVar() != UNDEFINED_ ) ?
VUB.getVal() * xlp[VUB.getVar()] : colUpperBound[indCol];
// Compute distances from current solution to upper and lower bounds
double delta = CoinMin(xlp[indCol] - LB, UB - xlp[indCol]);
// In case this variable is acceptable look for possible rows
if (delta > deltaMax) {
int iStart = colStarts[indCol];
int iStop = colStarts[indCol+1];
// int count = 0;
// std::vector<int> rowPossible;
// find a row to use in aggregation
for (int i = iStart; i < iStop; ++i) {
int rowInd = rowInds[i];
if (!setRowsAggregated.denseVector()[rowInd]) {
// if the row was not already selected, select it
RowType rType = rowTypes_[rowInd];
if ( ((rType == ROW_MIX) || (rType == ROW_CONT))
&& (fabs(coefByCol[i]) > EPSILON_) ) {
// rowPossible.push_back(rowInd);
rowSelected = rowInd;
deltaMax = delta;
colSelected = indCol;
foundRowToAggregate = true;
//count++;
break;
}
}
}
// if (count > 0)
// rowSelected = rowPossible[rand() % count];
// std::cout << count << std::endl;
}
}
return foundRowToAggregate;
}
//-------------------------------------------------------------------
// Aggregate the selected row with the current aggregated row
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::aggregateRow(
const int colSelected,
CoinIndexedVector& rowToAggregate, double rhs,
CoinIndexedVector& rowAggregated,
double& rhsAggregated ) const
{
// quantity to multiply by the coefficients of the row to aggregate
double multiCoef = rowAggregated[colSelected] / rowToAggregate[colSelected];
rowToAggregate *= multiCoef;
rhs *= multiCoef;
rowAggregated = rowAggregated - rowToAggregate;
rhsAggregated -= rhs;
}
//-------------------------------------------------------------------
// Choose the bound substitution based on the criteria defined by the user
//-------------------------------------------------------------------
inline bool
CglMixedIntegerRounding2::isLowerSubst(const double inf,
const double aj,
const double xlp,
const double LB,
const double UB) const
{
if (CRITERION_ == 1) {
// criterion 1 (the same as criterion (a) in the paper)
return xlp - LB < UB - xlp;
}
else {
if (UB == inf || xlp == LB)
return true;
if (LB == -inf || xlp == UB)
return false;
if (CRITERION_ == 2)
// criterion 2 (the same as criterion (b) in the paper)
return aj < 0;
else
// criterion 3 (the same as criterion (c) in the paper)
return aj > 0;
}
}
//-------------------------------------------------------------------
// Bound substitution heuristic
//-------------------------------------------------------------------
bool
CglMixedIntegerRounding2::boundSubstitution(
const OsiSolverInterface& si,
const CoinIndexedVector& rowAggregated,
const double* xlp,
const double* xlpExtra,
const double* colUpperBound,
const double* colLowerBound,
CoinIndexedVector& mixedKnapsack,
double& rhsMixedKnapsack, double& sStar,
CoinIndexedVector& contVariablesInS ) const
{
bool generated = false;
const int numColsAggregated = rowAggregated.getNumElements();
const int *rowAggregatedIndices = rowAggregated.getIndices();
const double *rowAggregatedElements = rowAggregated.denseVector();
// go through all the variables and if it is continuous and delta is
// negative, store variable in the vector contVariablesInS.
// If it is integer, store variable in the vector mixedKnapsack
int numCont = 0;
double infinity = si.getInfinity();
int j;
for ( j = 0; j < numColsAggregated; ++j) {
// get index and coefficient of column j in the aggregated row
const int indCol = rowAggregatedIndices[j];
const double coefCol = rowAggregatedElements[indCol];
// if the lower bound is equal to the upper bound, remove variable
if ( (indCol < numCols_) &&
(colLowerBound[indCol] == colUpperBound[indCol]) ) {
rhsMixedKnapsack -= coefCol * colLowerBound[indCol];
continue;
}
if (fabs(coefCol) < EPSILON_) continue;
// set the coefficients of the integer variables
if ( (indCol < numCols_) && (integerType_[indCol]) ) {
// Copy the integer variable to the vector mixedKnapsack
mixedKnapsack.add(indCol,coefCol);
continue;
}
// Select the continuous variables and copy the ones in s to
// the vector contVariablesInS
if (indCol < numCols_) { // variable is model variable
// Compute lower bound for variable indCol
const CglMixIntRoundVLB2 VLB = vlbs_[indCol];
const double LB = ( VLB.getVar() != UNDEFINED_ ) ?
VLB.getVal() * xlp[VLB.getVar()] : colLowerBound[indCol];
// Compute upper bound for variable indCol
const CglMixIntRoundVUB2 VUB = vubs_[indCol];
const double UB = ( VUB.getVar() != UNDEFINED_ ) ?
VUB.getVal() * xlp[VUB.getVar()] : colUpperBound[indCol];
// if both bounds are infinite, then we cannot form a mixed knapsack
if ( (LB == -1.0 * infinity) &&
(UB == infinity) ) {
#if CGL_DEBUG
std::cout << "continuous var with infinite bounds. " <<
"Cannot form mixed Knapsack = " << std::endl;
#endif
return generated;
}
// Select the bound substitution
if (isLowerSubst(infinity, rowAggregatedElements[indCol],
xlp[indCol], LB, UB)) {
if (VLB.getVar() != UNDEFINED_ ) {
const int indVLB = VLB.getVar();
mixedKnapsack.add(indVLB, coefCol * VLB.getVal());
}
else {
rhsMixedKnapsack -= coefCol * LB;
}
// Update sStar
if (coefCol < -EPSILON_) {
contVariablesInS.insert(indCol, coefCol);
sStar -= coefCol * (xlp[indCol] - LB);
numCont++;
}
}
else {
if (VUB.getVar() != UNDEFINED_ ) {
const int indVUB = VUB.getVar();
mixedKnapsack.add(indVUB, coefCol * VUB.getVal());
}
else {
rhsMixedKnapsack -= coefCol * UB;
}
// Update sStar
if (coefCol > EPSILON_) {
contVariablesInS.insert(indCol, - coefCol);
sStar += coefCol * (UB - xlp[indCol]);
numCont++;
}
}
}
else { // variable is slack variable
// in this case the LB = 0 and the UB = infinity
// Update sStar
const double tLB = xlpExtra[indCol - numCols_];
if (coefCol < -EPSILON_) {
contVariablesInS.insert(indCol, coefCol);
sStar -= coefCol * tLB;
numCont++;
}
}
}
// if there are no continuous variables to form s, then we stop
#if CGL_DEBUG
std::cout << "# of continuous var in mixedKnapsack = " << numCont <<
std::endl;
#endif
if (numCont == 0) return generated;
// check that the integer variables have lower bound equal to zero
const int numInt = mixedKnapsack.getNumElements();
// if there are not integer variables in mixedKnapsack, then we stop
// CAUTION: all the coefficients could be zero
#if CGL_DEBUG
std::cout << "# of integer var in mixedKnapsack = " << numInt <<
std::endl;
#endif
if (numInt == 0) return generated;
const int *knapsackIndices = mixedKnapsack.getIndices();
const double *knapsackElements = mixedKnapsack.denseVector();
for ( j = 0; j < numInt; ++j) {
int indCol = knapsackIndices[j];
// if the coefficient is zero, disregard
if (fabs(knapsackElements[indCol]) < EPSILON_) continue;
// if the lower bound is not zero, then we stop
if (fabs(colLowerBound[indCol]) > EPSILON_) return generated;
}
// if the lower bounds of all integer variables are zero, proceed
generated = true;
return generated;
}
//-------------------------------------------------------------------
// c-MIR separation heuristic
//-------------------------------------------------------------------
bool
CglMixedIntegerRounding2::cMirSeparation(
const OsiSolverInterface& si,
const CoinPackedMatrix& matrixByRow,
const CoinIndexedVector& rowAggregated,
const int* listRowsAggregated,
const char* sense, const double* RHS,
//const double* coefByRow,
//const int* colInds, const int* rowStarts,
const double* xlp, const double sStar,
const double* colUpperBound,
const double* colLowerBound,
const CoinIndexedVector& mixedKnapsack,
const double& rhsMixedKnapsack,
const CoinIndexedVector& contVariablesInS,
CoinIndexedVector * workVectors,
OsiRowCut& cMirCut) const
{
bool generated = false;
double numeratorBeta = rhsMixedKnapsack;
CoinIndexedVector * cMIR = &workVectors[0];
cMIR->copy(mixedKnapsack);
double rhscMIR;
double maxViolation = 0.0;
double bestDelta = 0.0;
CoinIndexedVector * bestCut = &workVectors[1];
double rhsBestCut = 0.0;
double sCoefBestCut = 0.0;
const int numInt = mixedKnapsack.getNumElements();
const int *knapsackIndices = mixedKnapsack.getIndices();
const double *knapsackElements = mixedKnapsack.denseVector();
const int *contVarInSIndices = contVariablesInS.getIndices();
const double *contVarInSElements = contVariablesInS.denseVector();
// Construct set C, T will be the rest.
// Also, for T we construct a CoinIndexedVector named complT which
// contains the vars in T that are strictly between their bounds
CoinIndexedVector & setC = workVectors[3];
setC.clear();
CoinIndexedVector * complT = &workVectors[2];
complT->clear();
double infinity = si.getInfinity();
int j;
for ( j = 0; j < numInt; ++j) {
const int indCol = knapsackIndices[j];
// if the upper bound is infinity, then indCol is in T and cannot
// be in complT
if (colUpperBound[indCol] != infinity) {
if (xlp[indCol] >= colUpperBound[indCol] / 2.0) {
setC.insert(j,1.0);
numeratorBeta -= knapsackElements[indCol] * colUpperBound[indCol];
} else {
if ( (xlp[indCol] <= EPSILON_) ||
(xlp[indCol] >= colUpperBound[indCol] - EPSILON_))
continue;
complT->insert(j, fabs(xlp[indCol] - colUpperBound[indCol]/2));
}
}
}
// Sort the indices in complT by nondecreasing values
// (which are $|y^*_j-u_j/2|$)
if (complT->getNumElements() > 0) {
complT->sortIncrElement();
}
// Construct c-MIR inequalities and take the one with the largest violation
for ( j = 0; j < numInt; ++j) {
int indCol = knapsackIndices[j];
if ( (xlp[indCol] <= EPSILON_) ||
(xlp[indCol] >= colUpperBound[indCol] - EPSILON_))
continue;
double delta = knapsackElements[indCol];
// delta has to be positive
if (delta <= EPSILON_) continue;
double violation = 0.0;
double sCoef = 0.0;
// form a cMIR inequality
cMirInequality(numInt, delta, numeratorBeta, knapsackIndices,
knapsackElements, xlp, sStar, colUpperBound, setC, *cMIR,
rhscMIR, sCoef, violation);
// store cut if it is the best found so far
if (violation > maxViolation + EPSILON_) {
bestCut->copy(*cMIR);
rhsBestCut = rhscMIR;
sCoefBestCut = sCoef;
maxViolation = violation;
bestDelta = delta;
}
}
// if no violated inequality has been found, exit now
if (maxViolation == 0.0) {
bestCut->clear();
return generated;
}
// improve the best violated inequality.
// try to divide delta by 2, 4 or 8 and see if increases the violation
double deltaBase = bestDelta;
for (int multFactor = 2; multFactor <= 8; multFactor *= 2) {
double delta = deltaBase / multFactor;
double violation = 0.0;
double sCoef = 0.0;
// form a cMIR inequality
cMirInequality(numInt, delta, numeratorBeta, knapsackIndices,
knapsackElements, xlp, sStar, colUpperBound, setC, *cMIR,
rhscMIR, sCoef, violation);
// store cut if it is the best found so far
if (violation > maxViolation + EPSILON_) {
bestCut->copy(cMIR);
rhsBestCut = rhscMIR;
sCoefBestCut = sCoef;
maxViolation = violation;
bestDelta = delta;
}
}
// improve cMIR for the best delta
// complT contains indices into mixedKnapsack for the variables
// which may be complemented and they are already appropriately
// sorted.
const int complTSize = complT->getNumElements();
if (complTSize > 0) {
const int *complTIndices = complT->getIndices();
for (int j = 0; j < complTSize; ++j) {
// move variable in set complT from set T to set C
int jIndex = complTIndices[j];
int indCol = knapsackIndices[jIndex];
// do nothing if upper bound is infinity
if (colUpperBound[indCol] >= infinity) continue;
setC.insert(jIndex,1.0);
double violation = 0.0;
double sCoef = 0.0;
double localNumeratorBeta = numeratorBeta -
mixedKnapsack[indCol] * colUpperBound[indCol];
// form a cMIR inequality
cMirInequality(numInt, bestDelta, localNumeratorBeta, knapsackIndices,
knapsackElements, xlp, sStar, colUpperBound, setC, *cMIR,
rhscMIR, sCoef, violation);
// store cut if it is the best found so far; otherwise, move the variable
// that was added to set C back to set T
if (violation > maxViolation + EPSILON_) {
bestCut->copy(cMIR);
rhsBestCut = rhscMIR;
sCoefBestCut = sCoef;
maxViolation = violation;
numeratorBeta = localNumeratorBeta;
}
else
setC.quickAdd(jIndex,-1.0);
}
}
// write the best cut found with the model variables
int numCont = contVariablesInS.getNumElements();
for ( j = 0; j < numCont; ++j) {
int indCol = contVarInSIndices[j];
double coefCol = contVarInSElements[indCol];
if (indCol < numCols_) { // variable is model variable
// Compute lower bound for variable indCol
CglMixIntRoundVLB2 VLB = vlbs_[indCol];
double LB = ( VLB.getVar() != UNDEFINED_ ) ?
VLB.getVal() * xlp[VLB.getVar()] : colLowerBound[indCol];
// Compute upper bound for variable indCol
CglMixIntRoundVUB2 VUB = vubs_[indCol];
double UB = ( VUB.getVar() != UNDEFINED_ ) ?
VUB.getVal() * xlp[VUB.getVar()] : colUpperBound[indCol];
// Select the bound substitution
if (isLowerSubst(infinity, rowAggregated[indCol],
xlp[indCol], LB, UB)) {
if (VLB.getVar() != UNDEFINED_ ) {
int indVLB = VLB.getVar();
bestCut->add(indVLB, - sCoefBestCut * coefCol * VLB.getVal());
bestCut->insert(indCol, sCoefBestCut * coefCol);
}
else {
rhsBestCut += sCoefBestCut * coefCol * colLowerBound[indCol];
bestCut->insert(indCol, sCoefBestCut * coefCol);
}
}
else {
if (VUB.getVar() != UNDEFINED_ ) {
int indVUB = VUB.getVar();
bestCut->add(indVUB, sCoefBestCut * coefCol * VUB.getVal());
bestCut->insert(indCol, - sCoefBestCut * coefCol);
}
else {
rhsBestCut -= sCoefBestCut * coefCol * colUpperBound[indCol];
bestCut->insert(indCol, - sCoefBestCut * coefCol);
}
}
}
else { // variable is slack variable
// in this case the LB = 0 and the UB = infinity
// copy the row selected to a vector of type CoinIndexedVector
const int iRow = listRowsAggregated[indCol - numCols_];
double multiplier;
if (sense[iRow] == 'L') {
// if it is a <= inequality, the coefficient of the slack is 1
multiplier = (- sCoefBestCut * coefCol);
}
else {
// if it is a <= inequality, the coefficient of the slack is -1
multiplier = (sCoefBestCut * coefCol);
}
rhsBestCut += RHS[iRow]*multiplier;
CoinShallowPackedVector row = matrixByRow.getVector(iRow);
int nElements = row.getNumElements();
const int * column = row.getIndices();
const double * element = row.getElements();
for (int i=0;i<nElements;i++)
bestCut->add(column[i],element[i]*multiplier);
}
}
// Check the violation of the cut after it is written with the original
// variables.
int cutLen = bestCut->getNumElements();
int* cutInd = bestCut->getIndices();
double* cutCoef = bestCut->denseVector();
double cutRHS = rhsBestCut;
double violation = 0.0;
double normCut = 0.0;
//double smallest=COIN_DBL_MAX;
double largest=0.0;
// Also weaken by small coefficients
for ( j = 0; j < cutLen; ++j) {
int column = cutInd[j];
double value = cutCoef[column];
//smallest=CoinMin(smallest,fabs(value));
largest=CoinMax(largest,fabs(value));
//normCut += value * value;
}
//normCut=sqrt(normCut);
//printf("smallest %g largest %g norm %g\n",
// smallest,largest,normCut);
double testValue=CoinMax(1.0e-6*largest,1.0e-12);
//normCut=0.0;
int n=0;
for ( j = 0; j < cutLen; ++j) {
int column = cutInd[j];
double value = cutCoef[column];
if (fabs(value)>testValue) {
violation += value * xlp[column];
normCut += value * value;
cutInd[n++]=column;
} else if (value) {
cutCoef[column]=0.0;
// Weaken
if (value>0.0) {
// Allow for at lower bound
cutRHS -= value*colLowerBound[column];
} else {
// Allow for at upper bound
cutRHS -= value*colUpperBound[column];
}
}
}
cutLen=n;
violation -= cutRHS;
violation /= sqrt(normCut);
if ( violation > TOLERANCE_ ) {
// cutCoef is still unpacked
std::sort(cutInd,cutInd+cutLen);
int i;
for ( i=0;i<cutLen;i++) {
int column=cutInd[i];
assert (cutCoef[column]);
double value = cutCoef[column];
cutCoef[column] =0.0;
cutCoef[i]=value;
}
cMirCut.setRow(cutLen, cutInd, cutCoef);
cMirCut.setLb(-1.0 * infinity);
cMirCut.setUb(cutRHS);
cMirCut.setEffectiveness(violation);
#ifdef CGL_DEBUG
{
for (int k=0; k<cutLen; k++){
assert(cutInd[k]>=0);
assert(cutCoef[k]);
assert (fabs(cutCoef[k])>1.0e-12);
}
}
#endif
// Zero bestCut by hand
bestCut->setNumElements(0);
for ( i=0;i<cutLen;i++) {
cutCoef[i]=0.0;
}
generated = true;
} else {
// just clear
bestCut->clear();
}
return generated;
}
//-------------------------------------------------------------------
// construct a c-MIR inequality
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::cMirInequality(
const int numInt,
const double delta,
const double numeratorBeta,
const int *knapsackIndices,
const double* knapsackElements,
const double* xlp,
const double sStar,
const double* colUpperBound,
const CoinIndexedVector& setC,
CoinIndexedVector& cMIR,
double& rhscMIR,
double& sCoef,
double& violation) const
{
// form a cMIR inequality
double beta = numeratorBeta / delta;
double f = beta - floor(beta);
rhscMIR = floor(beta);
double normCut = 0.0;
// coefficients of variables in set T
for (int i = 0; i < numInt; ++i) {
const int iIndex = knapsackIndices[i];
double G = 0.0;
if (setC.denseVector()[i] != 1.0) {
// i is not in setC, i.e., it is in T
G = functionG(knapsackElements[iIndex] / delta, f);
violation += (G * xlp[iIndex]);
normCut += G * G;
cMIR.setElement(i, G);
} else {
G = functionG( - knapsackElements[iIndex] / delta, f);
violation -= (G * xlp[iIndex]);
normCut += G * G;
rhscMIR -= G * colUpperBound[iIndex];
cMIR.setElement(i, -G);
}
}
sCoef = 1.0 / (delta * (1.0 - f));
violation -= (rhscMIR + sCoef * sStar);
normCut += sCoef * sCoef;
violation /= sqrt(normCut);
}
//-------------------------------------------------------------------
// function G for computing coefficients in cMIR inequality
//-------------------------------------------------------------------
inline double
CglMixedIntegerRounding2::functionG( const double d, const double f ) const
{
double delta = d - floor(d) - f;
if (delta > EPSILON_)
return floor(d) + delta / (1 - f);
else
return floor(d);
}
//-------------------------------------------------------------------
// Printing statistics
//-------------------------------------------------------------------
void
CglMixedIntegerRounding2::printStats(
std::ofstream & fout,
const bool hasCut,
const OsiSolverInterface& si,
const CoinIndexedVector& rowAggregated,
const double& rhsAggregated, const double* xlp,
const double* xlpExtra,
const int* listRowsAggregated,
const int* listColsSelected,
const int level,
const double* colUpperBound,
const double* colLowerBound ) const
{
const int numColsAggregated = rowAggregated.getNumElements();
const int *rowAggregatedIndices = rowAggregated.getIndices();
const double *rowAggregatedElements = rowAggregated.denseVector();
fout << "Rows ";
for (int i = 0; i < level; ++i) {
fout << listRowsAggregated[i] << " ";
}
fout << std::endl;
int numColsBack = 0;
// go through all the variables
for (int j = 0; j < numColsAggregated; ++j) {
// get index and coefficient of column j in the aggregated row
int indCol = rowAggregatedIndices[j];
double coefCol = rowAggregatedElements[indCol];
// check if a column used in aggregation is back into the aggregated row
for (int i = 0; i < level-1; ++i) {
if ( (listColsSelected[i] == indCol) && (coefCol != 0) ) {
numColsBack++;
break;
}
}
if (fabs(coefCol) < EPSILON_) {
// print variable number and coefficient
fout << indCol << " " << 0.0 << std::endl;
continue;
}
else {
// print variable number and coefficient
fout << indCol << " " << coefCol << " ";
}
// integer variables
if ( (indCol < numCols_) && (integerType_[indCol]) ) {
// print
fout << "I " << xlp[indCol] << " " << colLowerBound[indCol] <<
" " << colUpperBound[indCol] << std::endl;
continue;
}
// continuous variables
if (indCol < numCols_) { // variable is model variable
// print
fout << "C " << xlp[indCol] << " " << colLowerBound[indCol] <<
" " << colUpperBound[indCol] << " ";
// variable lower bound?
CglMixIntRoundVLB2 VLB = vlbs_[indCol];
if (VLB.getVar() != UNDEFINED_) {
fout << VLB.getVal() << " " << xlp[VLB.getVar()] << " " <<
colLowerBound[VLB.getVar()] << " " <<
colUpperBound[VLB.getVar()] << " ";
}
else {
fout << "-1 -1 -1 -1 ";
}
// variable upper bound?
CglMixIntRoundVUB2 VUB = vubs_[indCol];
if (VUB.getVar() != UNDEFINED_) {
fout << VUB.getVal() << " " << xlp[VUB.getVar()] << " " <<
colLowerBound[VUB.getVar()] << " " <<
colUpperBound[VUB.getVar()] << " ";
}
else {
fout << "-1 -1 -1 -1 ";
}
}
else { // variable is slack variable
// in this case the LB = 0 and the UB = infinity
// print
fout << "C " << xlpExtra[indCol-numCols_] << " " << 0.0 <<
" " << si.getInfinity() << " ";
}
fout << std::endl;
}
fout << "rhs " << rhsAggregated << std::endl;
fout << "numColsBack " << numColsBack << std::endl;
if (hasCut) {
fout << "CUT: YES" << std::endl;
}
else {
fout << "CUT: NO" << std::endl;
}
}
// This can be used to refresh any inforamtion
void
CglMixedIntegerRounding2::refreshSolver(OsiSolverInterface * solver)
{
if (solver->getNumRows()) {
mixIntRoundPreprocess(*solver);
doneInitPre_ = true;
} else {
doneInitPre_ = false;
}
}
// Create C++ lines to get to current state
std::string
CglMixedIntegerRounding2::generateCpp( FILE * fp)
{
CglMixedIntegerRounding2 other;
fprintf(fp,"0#include \"CglMixedIntegerRounding2.hpp\"\n");
fprintf(fp,"3 CglMixedIntegerRounding2 mixedIntegerRounding2;\n");
if (MAXAGGR_!=other.MAXAGGR_)
fprintf(fp,"3 mixedIntegerRounding2.setMAXAGGR_(%d);\n",MAXAGGR_);
else
fprintf(fp,"4 mixedIntegerRounding2.setMAXAGGR_(%d);\n",MAXAGGR_);
if (MULTIPLY_!=other.MULTIPLY_)
fprintf(fp,"3 mixedIntegerRounding2.setMULTIPLY_(%d);\n",MULTIPLY_);
else
fprintf(fp,"4 mixedIntegerRounding2.setMULTIPLY_(%d);\n",MULTIPLY_);
if (CRITERION_!=other.CRITERION_)
fprintf(fp,"3 mixedIntegerRounding2.setCRITERION_(%d);\n",CRITERION_);
if (doPreproc_!=other.doPreproc_)
fprintf(fp,"3 mixedIntegerRounding2.setDoPreproc_(%d);\n", doPreproc_);
if (getAggressiveness()!=other.getAggressiveness())
fprintf(fp,"3 mixedIntegerRounding2.setAggressiveness(%d);\n",getAggressiveness());
else
fprintf(fp,"4 mixedIntegerRounding2.setAggressiveness(%d);\n",getAggressiveness());
return "mixedIntegerRounding2";
}
void CglMixedIntegerRounding2::setDoPreproc(int value)
{
if (value != -1 && value != 0 && value != 1) {
throw CoinError("setDoPrepoc", "invalid value",
"CglMixedIntegerRounding2");
}
else {
doPreproc_ = value;
}
}
bool CglMixedIntegerRounding2::getDoPreproc() const
{
return (doPreproc_!=0);
}