limp-cbc-0.3.2.0: cbits/coin/CglRedSplit2.cpp
// Last edit: 04/03/10
//
// Name: CglRedSplit2.cpp
// Author: Giacomo Nannicini
// Singapore University of Technology and Design
// Singapore
// email: nannicini@sutd.edu.sg
// based on CglRedSplit by Francois Margot
// Date: 03/09/09
//-----------------------------------------------------------------------------
// Copyright (C) 2010, Giacomo Nannicini and others. All Rights Reserved.
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cassert>
#include <iostream>
/* Debug output */
//#define RS2_TRACE
/* Print optimal tableau and basis */
//#define RS2_TRACETAB
/* Use LAPACK instead of own code for solving linear systems */
//#define RS2_USE_LAPACK
// Sparse for intNonBasic 0,1,2
#define RS_FAST_INT 2
// Sparse for contNonBasic 0,1,2
#define RS_FAST_CONT 2
// Sparse for workNonBasic 0,1,2
#define RS_FAST_WORK 2
#include "CoinPragma.hpp"
#include "OsiSolverInterface.hpp"
#include "CglRedSplit2.hpp"
#include "CoinPackedVector.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinIndexedVector.hpp"
#include "OsiRowCutDebugger.hpp"
#include "CoinFactorization.hpp"
#include "CoinFinite.hpp"
#define TINY 1e-20
#define rs2round(x) (floor((x)+0.5))
//-------------------------------------------------------------------
// Generate Reduce-and-Split cuts
//-------------------------------------------------------------------
#ifdef RS2_USE_LAPACK
extern "C" void dgesv_( const int * , const int * , double * , const int * , int * , double * , const int * , int * );
#endif
/***************************************************************************/
// Utility functions and definitions for sorting
struct sortElement{
int index;
double cost;
};
// Return -1 if firstE has lower cost or same cost but lower index
int rs2_compareElements(const void* firstE, const void* secondE){
const struct sortElement* a = static_cast<const struct sortElement*>(firstE);
const struct sortElement* b = static_cast<const struct sortElement*>(secondE);
if (a->cost < b->cost)
return -1;
else if (a->cost > b->cost)
return 1;
else if (a->index < b->index)
return -1;
else if (a->index > b->index)
return 1;
return 0;
}
/***************************************************************************/
// Returns (value - floor) but allowing for small errors;
// taken from the CglGomory cut generator
inline double CglRedSplit2::rs_above_integer(const double value) const
{
double value2=floor(value);
double value3=rs2round(value);
if (fabs(value3-value)< param.getEPS() * (fabs(value3)+1.0))
return 0.0;
return value-value2;
} /* rs_above_integer */
/**********************************************************/
void CglRedSplit2::rs_allocmatINT(int ***v, int m, int n)
{
*v = reinterpret_cast<int **> (calloc (m, sizeof(int *)));
if (*v == NULL) {
printf("###ERROR: INTEGER matrix allocation failed\n");
exit(1);
}
for (int i = 0; i < m; ++i) {
(*v)[i] = reinterpret_cast<int *> (calloc (n, sizeof(int)));
if ((*v)[i] == NULL) {
printf("###ERROR: INTEGER matrix allocation failed\n");
exit(1);
}
}
} /* rs_allocmatINT */
/**********************************************************/
void CglRedSplit2::rs_deallocmatINT(int ***v, int m)
{
for (int i = 0; i < m; ++i) {
free(reinterpret_cast<void *> ((*v)[i]));
}
free(reinterpret_cast<void *> (*v));
} /* rs_deallocmatINT */
/**********************************************************/
void CglRedSplit2::rs_allocmatDBL(double ***v, int m, int n)
{
*v = reinterpret_cast<double **> (calloc (m, sizeof(double *)));
if (*v == NULL) {
printf("###ERROR: DOUBLE matrix allocation failed\n");
exit(1);
}
for (int i = 0; i < m; ++i){
(*v)[i] = reinterpret_cast<double *> (calloc (n, sizeof(double)));
if ((*v)[i] == NULL) {
printf("###ERROR: DOUBLE matrix allocation failed\n");
exit(1);
}
}
} /* rs_allocmatDBL */
/**********************************************************/
void CglRedSplit2::rs_deallocmatDBL(double ***v, int m)
{
for (int i = 0; i < m; ++i){
free(reinterpret_cast<void *> ((*v)[i]));
}
free(reinterpret_cast<void *> (*v));
} /* rs_deallocmatDBL */
/**********************************************************/
void CglRedSplit2::rs_printvecINT(const char *vecstr, const int *x, int n) const
{
int num, fromto, upto, j, i;
num = (n/10) + 1;
printf("%s :\n", vecstr);
for(j=0; j<num; j++){
fromto = 10*j;
upto = 10 * (j+1);
if(n <= upto) upto = n;
for(i=fromto; i<upto; i++)
printf(" %4d", x[i]);
printf("\n");
}
printf("\n");
} /* rs_printvecINT */
/**********************************************************/
void CglRedSplit2::rs_printvecDBL(const char *vecstr,
const double *x, int n) const
{
int num, fromto, upto, j, i;
num = (n/10) + 1;
printf("%s :\n", vecstr);
for(j=0; j<num; j++){
fromto = 10*j;
upto = 10 * (j+1);
if(n <= upto) upto = n;
for(i=fromto; i<upto; i++)
printf(" %7.5f", x[i]);
printf("\n");
}
printf("\n");
} /* rs_printvecDBL */
/**********************************************************/
void CglRedSplit2::rs_printmatINT(const char *vecstr, const int * const *x,
int m, int n) const
{
printf("%s :\n", vecstr);
for (int i = 0; i < m; ++i){
for (int j = 0; j < n; ++j){
printf(" %4d", x[i][j]);
}
printf("\n");
}
printf("\n");
} /* rs_printmatINT */
/**********************************************************/
void CglRedSplit2::rs_printmatDBL(const char *vecstr, const double * const *x, int m, int n) const
{
printf("%s :\n", vecstr);
for (int i = 0; i < m; ++i){
for (int j = 0; j < n; ++j){
printf(" %7.3f", x[i][j]);
}
printf("\n");
}
printf("\n");
} /* rs_printmatDBL */
/***************************************************************************/
double CglRedSplit2::rs_dotProd(const double *u, const double *v, int dim) const
{
double result = 0;
for (int i = 0; i < dim; ++i){
result += u[i] * v[i];
}
return(result);
} /* rs_dotProd */
/***************************************************************************/
static double rs_sparseDotProd(const double *u, const double *v,
const int * uInd, const int * vInd)
{
double result = 0;
if (uInd[0]<vInd[0]) {
int n=uInd[0];
uInd++;
for (int i = 0; i < n; ++i) {
int j=uInd[i];
result += u[j] * v[j];
}
} else {
int n=vInd[0];
vInd++;
for (int i = 0; i < n; ++i) {
int j=vInd[i];
result += u[j] * v[j];
}
}
return(result);
} /* rs_sparseDotProd */
/***************************************************************************/
double CglRedSplit2::rs_dotProd(const int *u, const double *v, int dim) const {
double result = 0;
for (int i=0; i < dim; ++i){
result += u[i] * v[i];
}
return(result);
} /* rs_dotProd */
/***************************************************************************/
// From Numerical Recipes in C: LU decomposition
int CglRedSplit2::ludcmp(double **a, int n, int *indx,
double *d, double* vv) const {
int ret = true;
int i, imax = 0, j, k;
double big, dum, sum, temp;
*d=1.0;
for (i = 1; i <= n; ++i){
big = 0.0;
for(j = 1; j <= n; j++) {
temp = fabs(a[i-1][j-1]);
if (temp > big) {
big = temp;
}
}
if (big == 0.0) {
return false;
}
vv[i-1]=1.0/big;
}
for (j = 1; j <= n; ++j){
for (i=1;i<j;i++) {
sum = a[i-1][j-1];
for(k = 1; k < i; k++) {
sum -= a[i-1][k-1] * a[k-1][j-1];
}
a[i-1][j-1]=sum;
}
big=0.0;
for (i = j; i <= n; i++) {
sum = a[i-1][j-1];
for (k = 1; k < j; k++) {
sum -= a[i-1][k-1]*a[k-1][j-1];
}
a[i-1][j-1] = sum;
dum = vv[i-1]*fabs(sum);
if (dum >= big) {
big = dum;
imax = i;
}
}
if (j != imax) {
for (k=1;k<=n;k++) {
dum=a[imax-1][k-1];
a[imax-1][k-1]=a[j-1][k-1];
a[j-1][k-1]=dum;
}
*d = -(*d);
vv[imax-1]=vv[j-1];
}
indx[j-1]=imax;
if (a[j-1][j-1] == 0.0) a[j-1][j-1]=TINY;
if (j != n) {
dum=1.0/(a[j-1][j-1]);
for (i=j+1;i<=n;i++) a[i-1][j-1] *= dum;
}
}
return ret;
}
/***************************************************************************/
// from Numerical Recipes in C: backward substitution
void CglRedSplit2::lubksb(double **a, int n, int *indx, double *b) const {
int i,ii=0,ip,j;
double sum;
for (i = 1; i <= n; ++i){
ip=indx[i-1];
sum=b[ip-1];
b[ip-1]=b[i-1];
if (ii)
for (j=ii;j<=i-1;j++) sum -= a[i-1][j-1]*b[j-1];
else if (sum) ii=i;
b[i-1]=sum;
}
for (i=n;i>=1;i--) {
sum=b[i-1];
for (j=i+1;j<=n;j++) sum -= a[i-1][j-1]*b[j-1];
b[i-1]=sum/a[i-1][i-1];
}
}
/***************************************************************************/
double CglRedSplit2::compute_norm_change(double oldnorm, const int* list,
int numElemList,
const double* multipliers) const {
double newnorm = 0;
double accumulator;
for (int j = 0; j < nTab; ++j){
accumulator = 0;
for (int i = 0; i < numElemList; ++i){
accumulator += multipliers[i]*workNonBasicTab[list[i]][j];
}
newnorm += accumulator*accumulator;
}
return (newnorm-oldnorm);
}
/***************************************************************************/
int CglRedSplit2::sort_rows_by_nonzeroes(struct sortElement* array,
int rowIndex, int maxRows,
int whichTab) const {
// Note: this function only takes into account rows which share at least
// a nonzero column with the row that we want to reduce (rowIndex).
int counter = 0;
int numZeroCost = 0;
for (int i = 0; i < mTab; ++i){
if (!checkTime()){
break;
}
if (numZeroCost == maxRows){
// We found enough rows with cost equal to zero, i.e. "perfect" rows
// No need to continue!
counter = numZeroCost;
break;
}
if ((i != rowIndex) && (norm[i] > param.getNormIsZero())){
// sort rows by number of nonzeros on nonbasic columns.
// check if they have at least one nonzero in the same place,
// otherwise skip
bool match = false;
for (int j = 0; j < nTab; ++j){
if ((fabs(workNonBasicTab[rowIndex][j]) > param.getEPS_COEFF()) &&
(fabs(workNonBasicTab[i][j]) > param.getEPS_COEFF())){
match = true;
break;
}
}
if (!match)
continue;
array[counter].index = i;
array[counter].cost = 0;
// whichTab = 0 means only intNonBasicTab, 1 means only
// workNonBasicTab, 2 means both
if (whichTab == 0 || whichTab == 2){
for (int j = 0; j < card_intNonBasicVar; ++j){
if ((fabs(intNonBasicTab[rowIndex][j]) <= param.getEPS_COEFF()) &&
(fabs(intNonBasicTab[i][j]) > param.getEPS_COEFF())){
array[counter].cost += 1;
}
}
}
if (whichTab == 1 || whichTab == 2){
for (int j = 0; j < nTab; ++j){
if ((fabs(workNonBasicTab[rowIndex][j]) <= param.getEPS_COEFF()) &&
(fabs(workNonBasicTab[i][j]) > param.getEPS_COEFF())){
array[counter].cost += 1;
}
}
}
if (array[counter].cost == 0){
array[counter] = array[numZeroCost];
array[numZeroCost].index = i;
array[numZeroCost].cost = 0;
numZeroCost++;
}
counter++;
}
}
if (counter > maxRows){
qsort(array, counter, sizeof(struct sortElement), rs2_compareElements);
}
return counter;
}
/***************************************************************************/
int CglRedSplit2::sort_rows_by_nonzeroes_greedy(struct sortElement* array,
int rowIndex, int maxRows,
int whichTab) const{
int counter = 0;
counter = sort_rows_by_nonzeroes(array, rowIndex, maxRows, whichTab);
if (counter > maxRows){
int i, j;
// vector of positions of zero elements
int* z_int = NULL;
int* z_cont = NULL;
// whichTab = 0 means only intNonBasicTab, 1 means only
// workNonBasicTab, 2 means both
if (whichTab == 0 || whichTab == 2)
z_int = new int[card_intNonBasicVar];
if (whichTab == 1 || whichTab == 2)
z_cont = new int[nTab];
// number of elements in the vectors above
int numz_int = 0;
int numz_cont = 0;
// compute initial vector of nonzeroes
if (whichTab == 0 || whichTab == 2){
for (j = 0; j < card_intNonBasicVar; ++j){
if (fabs(intNonBasicTab[rowIndex][j]) <= param.getEPS_COEFF()){
z_int[numz_int++] = j;
}
}
}
if (whichTab == 1 || whichTab == 2){
for (j = 0; j < nTab; ++j){
if (fabs(workNonBasicTab[rowIndex][j]) <= param.getEPS_COEFF()){
z_cont[numz_cont++] = j;
}
}
}
int numSorted = counter;
counter = 1;
// store the minimum number of nz introduced by one row
int minNz;
// and the position of that row
int rowMinNz;
// maximum number of nonzero that can be introduced
double maxNz;
int currentNz;
while (counter < maxRows && counter < numSorted){
if (!checkTime()){
break;
}
// initialize minNz to a large value
minNz = numz_int + numz_cont;;
maxNz = array[counter].cost + array[counter-1].cost;
rowMinNz = counter;
for (i = counter; i < numSorted && array[i].cost < maxNz; ++i){
int ii = array[i].index;
// pick the row which introduces the smallest
// number of nonzeros on the nonbasic integer columns
currentNz = 0;
for (j = 0; j < numz_int; ++j){
if ((fabs(intNonBasicTab[ii][z_int[j]]) > param.getEPS_COEFF())){
currentNz++;
}
}
for (j = 0; j < numz_cont; ++j){
if ((fabs(workNonBasicTab[ii][z_cont[j]]) > param.getEPS_COEFF())){
currentNz++;
}
}
array[i].cost = currentNz;
if (currentNz < minNz){
rowMinNz = i;
minNz = currentNz;
}
if (currentNz == 0){
// perfect matching, no need to look further
break;
}
}
// update introduced nonzeroes by the selected row;
// first, put the row in the correct position in the array
int swap1 = array[rowMinNz].index;
double swap2 = array[rowMinNz].cost;
array[rowMinNz] = array[counter];
array[counter].index = swap1;
array[counter].cost = swap2;
// now count how many nonzeroes were introduces, and update
// the vectors z_int and z_cont accordingly
for (j = 0; j < numz_int; ++j){
if (fabs(intNonBasicTab[swap1][z_int[j]]) > param.getEPS_COEFF()){
z_int[j] = z_int[numz_int-1];
numz_int--;
}
}
for (j = 0; j < numz_cont; ++j){
if (fabs(workNonBasicTab[swap1][z_cont[j]]) > param.getEPS_COEFF()){
z_cont[j] = z_cont[numz_cont-1];
numz_cont--;
}
}
counter++;
}
if (z_int)
delete[] z_int;
if (z_cont)
delete[] z_cont;
}
return counter;
}
/***************************************************************************/
int CglRedSplit2::sort_rows_by_cosine(struct sortElement* array,
int rowIndex, int maxRows,
int whichTab) const {
// whichTab == 0 means only integer nonbasic variables, whichTab ==
// 1 means only workTab, whichTab == 2 means both
int counter = 0;
double initnorm = 0.0;
if (whichTab == 0 || whichTab == 2){
#if RS_FAST_INT == 0
initnorm += rs_dotProd(intNonBasicTab[rowIndex],
intNonBasicTab[rowIndex],
card_intNonBasicVar);
#elif RS_FAST_INT == 1
initnorm += rs_dotProd(intNonBasicTab[rowIndex],
intNonBasicTab[rowIndex],
card_intNonBasicVar);
double value = rs_sparseDotProd(intNonBasicTab[rowIndex],
intNonBasicTab[rowIndex],
pi_mat[rowIndex]+mTab,
pi_mat[rowIndex]+mTab);
assert (value==initnorm);
#else
initnorm = rs_sparseDotProd(intNonBasicTab[rowIndex],
intNonBasicTab[rowIndex],
pi_mat[rowIndex]+mTab,
pi_mat[rowIndex]+mTab);
#endif
}
if (whichTab == 1 || whichTab == 2){
initnorm += norm[rowIndex];
}
#if RS_FAST_WORK
int workOffset = mTab+card_intNonBasicVar+card_contNonBasicVar+2;
#endif
for (int i = 0; i < mTab; ++i){
if ((i != rowIndex) && (norm[i] > param.getNormIsZero())){
if (!checkTime()){
break;
}
array[counter].index = i;
array[counter].cost = 0;
if (whichTab == 0 || whichTab == 2){
#if RS_FAST_INT == 0
array[counter].cost -= fabs(rs_dotProd(intNonBasicTab[rowIndex],
intNonBasicTab[i],
card_intNonBasicVar));
#elif RS_FAST == 1
array[counter].cost -= fabs(rs_dotProd(intNonBasicTab[rowIndex],
intNonBasicTab[i],
card_intNonBasicVar));
double value = -fabs(rs_sparseDotProd(intNonBasicTab[rowIndex],
intNonBasicTab[i],
pi_mat[rowIndex]+mTab,
pi_mat[i]+mTab));
assert (value==array[counter].cost);
#else
array[counter].cost = -fabs(rs_sparseDotProd(intNonBasicTab[rowIndex],
intNonBasicTab[i],
pi_mat[rowIndex]+mTab,
pi_mat[i]+mTab));
#endif
}
if (whichTab == 1 || whichTab == 2){
#if RS_FAST_WORK == 0
array[counter].cost -= fabs(rs_dotProd(workNonBasicTab[rowIndex],
workNonBasicTab[i], nTab));
#elif RS_FAST_WORK == 1
double value = array[counter].cost;
array[counter].cost -= fabs(rs_dotProd(workNonBasicTab[rowIndex],
workNonBasicTab[i], nTab));
value -= fabs(rs_sparseDotProd(workNonBasicTab[rowIndex],
workNonBasicTab[i],
pi_mat[rowIndex]+workOffset,
pi_mat[i]+workOffset));
assert (value==array[counter].cost);
#else
array[counter].cost -= fabs(rs_sparseDotProd(workNonBasicTab[rowIndex],
workNonBasicTab[i],
pi_mat[rowIndex]+workOffset,
pi_mat[i]+workOffset));
#endif
}
// Now divide by the two norms
double denom = 0.0;
if (whichTab == 0 || whichTab == 2){
#if RS_FAST_INT == 0
denom += initnorm*rs_dotProd(intNonBasicTab[i], intNonBasicTab[i],
card_intNonBasicVar);
#elif RS_FAST == 1
denom += initnorm*rs_dotProd(intNonBasicTab[i], intNonBasicTab[i],
card_intNonBasicVar);
double value = initnorm*rs_sparseDotProd(intNonBasicTab[i],
intNonBasicTab[i],
pi_mat[i]+mTab,
pi_mat[i]+mTab);
assert (value==denom);
#else
denom = initnorm*rs_sparseDotProd(intNonBasicTab[i],
intNonBasicTab[i],
pi_mat[i]+mTab,
pi_mat[i]+mTab);
#endif
}
if (whichTab == 1 || whichTab == 2){
denom += initnorm*norm[i];
}
array[counter].cost /= sqrt(denom);
if (array[counter].cost != 0.0){
counter++;
}
}
}
if (counter >= maxRows){
qsort(array, counter, sizeof(struct sortElement), rs2_compareElements);
}
return counter;
}
/***************************************************************************/
int CglRedSplit2::get_list_rows_reduction(int rowIndex, int maxRowsReduction,
int* list, const double* norm,
CglRedSplit2Param::RowSelectionStrategy selectionStrategy) const{
#ifdef RS2_TRACE
printf("Obtaining list of rows with column strategy %d\n", selectionStrategy);
#endif
struct sortElement* array = new struct sortElement[mTab];
int counter = 0;
int i, j;
// Look in CglRedSplit2Param for a description of each strategy
if (selectionStrategy == CglRedSplit2Param::RS1){
counter = sort_rows_by_nonzeroes(array, rowIndex, maxRowsReduction-1, 0);
}
else if (selectionStrategy == CglRedSplit2Param::RS2){
counter = sort_rows_by_nonzeroes(array, rowIndex, maxRowsReduction-1, 1);
}
else if (selectionStrategy == CglRedSplit2Param::RS3){
counter = sort_rows_by_nonzeroes(array, rowIndex, maxRowsReduction-1, 2);
}
else if (selectionStrategy == CglRedSplit2Param::RS4){
// compute the *true* number of nonzeroes that we introduce
// on the integer nonbasics
counter = sort_rows_by_nonzeroes_greedy(array, rowIndex, maxRowsReduction-1, 0);
}
else if (selectionStrategy == CglRedSplit2Param::RS5){
// compute the *true* number of nonzeroes that we introduce
// on the working set (continuous nonbasics)
counter = sort_rows_by_nonzeroes_greedy(array, rowIndex, maxRowsReduction-1, 1);
}
else if (selectionStrategy == CglRedSplit2Param::RS6){
// compute the *true* number of nonzeroes that we introduce
// on all the nonbasics: intNonBasicVar and workNonBasicVar
counter = sort_rows_by_nonzeroes_greedy(array, rowIndex, maxRowsReduction-1, 2);
}
else if (selectionStrategy == CglRedSplit2Param::RS7){
// sort rows by cosine of the angle in the space of the nonbasic integer
// and nonbasic continuous (working set) columns, wrt row rowIndex
counter = sort_rows_by_cosine(array, rowIndex, maxRowsReduction-1, 2);
}
else if (selectionStrategy == CglRedSplit2Param::RS8){
// sort rows by cosine of the angle in the space of the
// nonbasic continuous (working set) columns, wrt row rowIndex
counter = sort_rows_by_cosine(array, rowIndex, maxRowsReduction-1, 1);
}
list[0] = rowIndex;
j = 1;
for (i = 0; i < counter && j < maxRowsReduction; ++i){
list[j] = array[i].index;
j++;
}
delete[] array;
return j;
}
/***************************************************************************/
void CglRedSplit2::fill_workNonBasicTab(CglRedSplit2Param::ColumnSelectionStrategy strategy, const int* ignore_list){
// Sort continuous non basic variables by reduced cost and choose the
// columns we will work with.
// We can choose among several strategies; these are described in
// CglRedSplit2Param.hpp.
int i, j;
if (strategy == CglRedSplit2Param::CS_ALLCONT){
// select all continuous nonbasic columns
for (i = 0; i < mTab; ++i){
memcpy(workNonBasicTab[i], contNonBasicTab[i],
card_contNonBasicVar*sizeof(double));
}
nTab = card_contNonBasicVar;
}
else{
// Begin by sorting all continuous nonbasic columns by increasing
// reduced cost (except those in the ingore_list)
struct sortElement* array = new struct sortElement[card_contNonBasicVar];
int pos = 0, iter = 0;
for (i = 0; i < card_contNonBasicVar; ++i){
if (ignore_list != NULL){
iter = 0;
while (ignore_list[iter] >= 0 &&
ignore_list[iter] != contNonBasicVar[i]){
iter++;
}
if (ignore_list[iter] == contNonBasicVar[i]){
continue;
}
}
array[pos].index = i;
// We take absolute value of reduced costs because sign could be
// positive or negative depending on minimization/maximization,
// and which bound the variable is at (upper/lower).
if (contNonBasicVar[i] < ncol)
array[pos].cost = fabs(reducedCost[contNonBasicVar[i]]);
else
array[pos].cost = fabs(rowPrice[contNonBasicVar[i]-ncol]);
pos++;
}
qsort(array,pos,sizeof(struct sortElement),rs2_compareElements);
nTab = 0;
int card_contNonBasicVar = pos;
// Now choose the variables. These strategies are described in
// CglRedSplit2Param.hpp.
// CS1, CS2 and CS3 are the 3 partitions of C-3P: each one
// contains one third of the variables, sorted by reduced costs
if (strategy == CglRedSplit2Param::CS1){
for (i = 0; i < card_contNonBasicVar/3; ++i){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
else if (strategy == CglRedSplit2Param::CS2){
for (i = card_contNonBasicVar/3; i < 2*card_contNonBasicVar/3; ++i){
if (!checkTime()){
break;
}
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
else if (strategy == CglRedSplit2Param::CS3){
for (i = 2*card_contNonBasicVar/3; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
// CS4, CS5, CS6, CS7 and CS8 are the 5 partitions of C-5P; each
// one contains 1/5 of the variables, sorted by increasing reduced
// cost: first set has smallest reduced costs, and so on.
else if (strategy == CglRedSplit2Param::CS4){
for (i = 0; i < card_contNonBasicVar/5; ++i){
if (!checkTime()){
break;
}
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
else if (strategy == CglRedSplit2Param::CS5){
for (i = card_contNonBasicVar/5; i < 2*card_contNonBasicVar/5; ++i){
if (!checkTime()){
break;
}
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
else if (strategy == CglRedSplit2Param::CS6){
for (i = 2*card_contNonBasicVar/5; i < 3*card_contNonBasicVar/5; ++i){
if (!checkTime()){
break;
}
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
else if (strategy == CglRedSplit2Param::CS7){
for (i = 3*card_contNonBasicVar/5; i < 4*card_contNonBasicVar/5; ++i){
if (!checkTime()){
break;
}
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
else if (strategy == CglRedSplit2Param::CS8){
for (i = 4*card_contNonBasicVar/5; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
// CS9 and CS10 are the 2 partitions of I-2P-2/3; this is a
// partition with interleaving pattern XX--X- on the 2/3 of the
// variables with smallest reduced cost
else if (strategy == CglRedSplit2Param::CS9){
int pat;
for (i = 0; i < 2*card_contNonBasicVar/3; ++i){
if (!checkTime()){
break;
}
pat = i%6;
if (pat == 0 || pat == 1 || pat == 4){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS10){
int pat;
for (i = 0; i < 2*card_contNonBasicVar/3; ++i){
if (!checkTime()){
break;
}
pat = i%6;
if (pat == 2 || pat == 3 || pat == 5){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
// CS11 and CS12 are the 2 partitions of I-2P-4/5; partition with
// interleaving pattern X---XXX- on 4/5 of the vars with smallest
// reduced cost
else if (strategy == CglRedSplit2Param::CS11){
int pat;
for (i = 0; i < 4*card_contNonBasicVar/5; ++i){
if (!checkTime()){
break;
}
pat = i%8;
if (pat == 0 || pat == 4 || pat == 5 || pat == 6){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS12){
int pat;
for (i = 0; i < 4*card_contNonBasicVar/5; ++i){
if (!checkTime()){
break;
}
pat = i%8;
if (pat == 1 || pat == 2 || pat == 3 || pat == 7){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
// CS13, and CS14 are the 2 partitions of I-2P-1/2; interleaving
// partition with pattern X--X on 1/2 of the variables
else if (strategy == CglRedSplit2Param::CS13){
for (i = 0; i < card_contNonBasicVar/2; ++i){
if (!checkTime()){
break;
}
if (i%4 == 0 || i%4 == 3){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS14){
for (i = 0; i < card_contNonBasicVar/2; ++i){
if (!checkTime()){
break;
}
if (i%4 == 1 || i%4 == 2){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
// CS15, CS16 and CS17 are the 3 partitions of I-3P; interleaving
// partition with pattern X-/ on all the variables
else if (strategy == CglRedSplit2Param::CS15){
for (i = 0; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
if (i%3 == 0){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS16){
for (i = 0; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
if (i%3 == 1){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS17){
for (i = 0; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
if (i%3 == 2){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
// CS18, CS19, CS20, CS21 are the 4 partitions of I-4P;
// interleaving partition with pattern X-X/\\-/ on all the
// variables
else if (strategy == CglRedSplit2Param::CS18){
int pat;
for (i = 0; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
pat = i%8;
if (pat == 0 || pat == 2){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS19){
int pat;
for (i = 0; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
pat = i%8;
if (pat == 1 || pat == 6){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS20){
int pat;
for (i = 0; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
pat = i%8;
if (pat == 3 || pat == 7){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
else if (strategy == CglRedSplit2Param::CS21){
int pat;
for (i = 0; i < card_contNonBasicVar; ++i){
if (!checkTime()){
break;
}
pat = i%8;
if (pat == 4 || pat == 5){
for (j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = contNonBasicTab[j][array[i].index];
}
nTab++;
}
}
}
delete[] array;
}
// create sparse work stuff
for(i=0; i<mTab; i++) {
#if RS_FAST_WORK == 0
norm[i] = rs_dotProd(workNonBasicTab[i], workNonBasicTab[i], nTab);
#elif RS_FAST_WORK == 1
int * info = pi_mat[i];
int * workSparse=info+mTab+2+card_intNonBasicVar+card_contNonBasicVar;
int n = 0;
const double * values = workNonBasicTab[i];
double checkNorm=0.0;
for (int j=0;j<nTab;j++) {
if (values[j]) {
workSparse[++n]=j;
checkNorm += values[j]*values[j];
}
}
workSparse[0]=n;
norm[i] = rs_dotProd(workNonBasicTab[i], workNonBasicTab[i], nTab);
assert (checkNorm==norm[i]);
#else
int * info = pi_mat[i];
int * workSparse=info+mTab+2+card_intNonBasicVar+card_contNonBasicVar;
int n = 0;
const double * values = workNonBasicTab[i];
double checkNorm=0.0;
for (int j=0;j<nTab;j++) {
if (values[j]) {
workSparse[++n]=j;
checkNorm += values[j]*values[j];
}
}
workSparse[0]=n;
norm[i] = checkNorm;
#endif
}
}
/***************************************************************************/
void CglRedSplit2::reduce_workNonBasicTab(int numRowsReduction,
CglRedSplit2Param::RowSelectionStrategy rowSelectionStrategy,
int maxIterations) {
// Use at most this number of rows
int maxRowsReduction = CoinMin(numRowsReduction, mTab);
if (maxRowsReduction == 1){
return;
}
int i, j, k, h;
// Allocate space to store the matrix for the linear system
#ifdef RS2_USE_LAPACK
double* A = new double[maxRowsReduction*maxRowsReduction];
#else
double** A;
rs_allocmatDBL(&A, maxRowsReduction, maxRowsReduction);
#endif
// Right hand side
double* b = new double[maxRowsReduction];
// Data for LU decomposition
int* indexlu = new int[maxRowsReduction];
double tmpnumlu = 0.0;
double* tmpveclu = new double[maxRowsReduction];
// List of rows involved in the linear combination
int* list = new int[maxRowsReduction];
// Number of rows actually used
int numUsedRows;
bool resolveWithNormalization = false;
for (k = 0; k < mTab && k < maxIterations; ++k) {
if (!checkTime()){
break;
}
if (norm[k] > param.getNormIsZero()) {
// Obtain the list of rows that should be used to reduce row k
numUsedRows = get_list_rows_reduction(k, maxRowsReduction, list, norm,
rowSelectionStrategy);
#ifdef RS2_TRACE
rs_printvecINT("rows used for reduction", list, numUsedRows);
#endif
if (numUsedRows <= 1){
// This means that the list only contains the current row;
// Thus, no rows were selected for reduction. Skip.
continue;
}
// Note: the list must have size maxRowsReduction.
// Now prepare the linear system according to the paper
for (i = 0; i < numUsedRows; ++i){
for (j = 0; j < numUsedRows; ++j){
#ifdef RS2_USE_LAPACK
A[i*numUsedRows+j] = 0;
#else
A[i][j] = 0;
#endif
if (list[i] != k && list[j] != k){
for (h = 0; h < nTab; ++h){
#ifdef RS2_USE_LAPACK
A[i*numUsedRows+j]+=workNonBasicTab[list[i]][h]*workNonBasicTab[list[j]][h];
#else
A[i][j]+=workNonBasicTab[list[i]][h]*workNonBasicTab[list[j]][h];
#endif
}
if (resolveWithNormalization && i == j){
// Penalize the norm of lambda, i.e. the solution
#ifdef RS2_USE_LAPACK
A[i*numUsedRows+j] += norm[k]*param.getNormalization();
#else
A[i][j] += norm[k]*param.getNormalization();
#endif
}
}
}
if (list[i] == k){
b[i] = 1;
#ifdef RS2_USE_LAPACK
A[i*numUsedRows+i] = 1;
#else
A[i][i] = 1;
#endif
}
else{
b[i] = 0;
for (h = 0; h < nTab; ++h){
b[i] -= workNonBasicTab[list[i]][h]*workNonBasicTab[k][h];
}
}
}
// Linear system has been written, now solve it
#ifdef RS2_USE_LAPACK
int info = 0;
int h = 1;
dgesv_(&numUsedRows, &h, A, &numUsedRows, indexlu,
b, &numUsedRows, &info);
if (info)
continue;
#else
// LU decomposition
if (!ludcmp(A, numUsedRows, indexlu, &tmpnumlu, tmpveclu)){
// numerical error: exit
continue;
}
// Backward substitution
lubksb(A, numUsedRows, indexlu, b);
#endif
// Check the 1-norm of the solution
double sumnorm = 0.0;
for (i = 0; i < numUsedRows; ++i){
b[i] = rs2round(b[i]);
sumnorm += fabs(b[i]);
if (sumnorm > param.getMaxSumMultipliers()){
break;
}
}
#ifdef RS2_TRACE
rs_printvecDBL("multipliers", b, numUsedRows);
printf("sumnorm: %f\n", sumnorm);
#endif
if (sumnorm == 1){
// The optimal solution has lambda_k = 1 and all the rest is zero;
// so we cannot reduce any coefficient, let's skip this
continue;
}
else if (!resolveWithNormalization &&
sumnorm > param.getMaxSumMultipliers()){
// The norm of lambda is too large, resolve with normalization
// (note that we do not want to do this more than once)
resolveWithNormalization = true;
k--;
continue;
}
resolveWithNormalization = false;
if (sumnorm > param.getMaxSumMultipliers()){
// If we got this far, even resolving did not help, so we skip
continue;
}
double deltaNorm = compute_norm_change(norm[k], list,
numUsedRows, b);
if (deltaNorm <= -norm[k]*param.getMinNormReduction()){
for (i = 0; i < numUsedRows; ++i){
pi_mat[k][list[i]] = (int)(b[i]);
}
numRedRows++;
}
} /* if (norm[k] > param.getNormIsZero()) */
} /*for (k = 0; k < mTab; ++k) */
delete[] b;
delete[] list;
delete[] indexlu;
delete[] tmpveclu;
#ifdef RS2_USE_LAPACK
delete[] A;
#else
rs_deallocmatDBL(&A, maxRowsReduction);
#endif
#ifdef RS2_TRACE
sum_norms = 0;
for(i=0; i<mTab; i++) {
sum_norms += norm[i];
}
printf("CglRedSplit2::reduce_contNonBasicTab():Final sum of norms: %f\n", sum_norms);
printf("CglRedSplit2::reduce_contNonBasicTab(): %d rows reduced\n",numRedRows);
#endif
} /* reduce_workNonBasicTab */
/************************************************************************/
void CglRedSplit2::generate_row(int index_row, double *row) {
memset(row, 0, (ncol+nrow)*sizeof(double));
// we only deal with nonbasic variables - the cut coefficient will be zero
// on the basic ones anyway
#if RS_FAST_INT == 0 && RS_FAST_CONT == 0 && RS_FAST_WORK == 0
int i;
int locind, j;
for(i=0; i<card_intNonBasicVar; i++) {
locind = intNonBasicVar[i];
row[locind] = 0;
for(j=0; j<mTab; j++) {
row[locind] += pi_mat[index_row][j] * intNonBasicTab[j][i];
}
}
for(i=0; i<card_contNonBasicVar; i++) {
locind = contNonBasicVar[i];
row[locind] = 0;
for(j=0; j<mTab; j++) {
row[locind] += pi_mat[index_row][j] * contNonBasicTab[j][i];
}
}
#elif RS_FAST_INT < 2 || RS_FAST_CONT < 2
double * rowTemp = new double[ncol+nrow];
memset(rowTemp, 0, (ncol+nrow)*sizeof(double));
for(int i=0; i<card_intNonBasicVar; i++) {
int locind = intNonBasicVar[i];
rowTemp[locind] = 0;
for(int j=0; j<mTab; j++) {
rowTemp[locind] += pi_mat[index_row][j] * intNonBasicTab[j][i];
}
}
for(int i=0; i<card_contNonBasicVar; i++) {
int locind = contNonBasicVar[i];
rowTemp[locind] = 0;
for(int j=0; j<mTab; j++) {
rowTemp[locind] += pi_mat[index_row][j] * contNonBasicTab[j][i];
}
}
const int * pi = pi_mat[index_row];
#if RS_FAST_CONT
int contOffset = mTab+card_intNonBasicVar+1;
#endif
for(int j=0; j<mTab; j++) {
double value = pi[j];
if (value) {
const double * tableau = intNonBasicTab[j];
#if RS_FAST_INT < 1
for(int i=0; i<card_intNonBasicVar; i++) {
int locind = intNonBasicVar[i];
row[locind] += value * tableau[i];
}
#else
const int * which = pi_mat[j]+mTab;
int n = which[0];
which++;
for (int k=0;k<n;k++) {
int i=which[k];
int locind = intNonBasicVar[i];
row[locind] += value * tableau[i];
}
#endif
tableau = contNonBasicTab[j];
#if RS_FAST_CONT < 1
for(int i=0; i<card_contNonBasicVar; i++) {
int locind = contNonBasicVar[i];
row[locind] += value * tableau[i];
}
#else
const int * which2 = pi_mat[j]+contOffset;
int n2 = which2[0];
which2++;
for (int k=0;k<n2;k++) {
int i=which2[k];
int locind = contNonBasicVar[i];
row[locind] += value * tableau[i];
}
#endif
}
}
for (int i=0;i<nrow+ncol;i++) {
assert (fabs(row[i]-rowTemp[i])<1.0e-11+1.0e-8*CoinMax(fabs(row[i]),fabs(rowTemp[i])));
}
delete [] rowTemp;
#else
const int * pi = pi_mat[index_row];
int contOffset = mTab+card_intNonBasicVar+1;
for(int j=0; j<mTab; j++) {
double value = pi[j];
if (value) {
const double * tableau = intNonBasicTab[j];
const int * which = pi_mat[j]+mTab;
int n = which[0];
which++;
for (int k=0;k<n;k++) {
int i=which[k];
int locind = intNonBasicVar[i];
row[locind] += value * tableau[i];
}
tableau = contNonBasicTab[j];
which = pi_mat[j]+contOffset;
n = which[0];
which++;
for (int k=0;k<n;k++) {
int i=which[k];
int locind = contNonBasicVar[i];
row[locind] += value * tableau[i];
}
}
}
#endif
} /* generate_row */
/************************************************************************/
int CglRedSplit2::generate_cgcut(double *row, double *rhs) {
#ifdef RS2_TRACE
rs_printvecDBL("CglRedSplit2::generate_cgcut(): starting row",
row, ncol+nrow);
printf("rhs: %f\n", *rhs);
#endif
// Note that basic integer variables with fractional value
// other than basic_ind might have a non zero integer coefficient in the
// "real" row. However, these coefficients become zero anyway in the
// cut. It is thus assumed that all basic integer variables have a
// coefficient of zero in the given row (enforced in generate_row()).
//
// Other integer variables and basic continuous variables
// have a coefficient of zero, as the corresponding row is not selected
// for combination.
double f0 = rs_above_integer(*rhs);
double f0compl = 1 - f0;
// See Wolsey "Integer Programming" (1998), p. 130, fourth line from top
// after correcting typo (Proposition 8.8), flipping all signs to get <=.
if((f0 < param.getAway()) || (f0compl < param.getAway())) {
return(0);
}
int i, locind;
double f;
for(i=0; i<card_intNonBasicVar; i++) {
locind = intNonBasicVar[i];
f = rs_above_integer(row[locind]);
if(f > f0) {
row[locind] = - ((1-f) * f0);
}
else {
row[locind] = -(f*f0compl);
}
}
for(i=0; i<card_contNonBasicVar; i++) {
locind = contNonBasicVar[i];
if(row[locind] < 0) {
row[locind] *= f0;
}
else {
row[locind] = -(row[locind]*f0compl);
}
}
(*rhs) = -f0*f0compl;
#ifdef RS2_TRACE
rs_printvecDBL("CglRedSplit2::generate_cgcut(): row", row, ncol+nrow);
printf("rhs: %f\n", *rhs);
#endif
return(1);
} /* generate_cgcut */
/************************************************************************/
void CglRedSplit2::eliminate_slacks(double *row,
const double *elements,
const int *rowStart,
const int *indices,
const int *rowLength,
const double *rhs, double *tabrowrhs) {
int i, j;
for(i=0; i<nrow; i++) {
if(fabs(row[ncol+i]) > param.getEPS_ELIM()) {
int upto = rowStart[i] + rowLength[i];
for(j=rowStart[i]; j<upto; j++) {
row[indices[j]] -= row[ncol+i] * elements[j];
}
*tabrowrhs -= row[ncol+i] * rhs[i];
}
}
} /* eliminate_slacks */
/************************************************************************/
void CglRedSplit2::flip(double *row) {
int i;
for(i=0; i<card_nonBasicAtUpper; i++) {
row[nonBasicAtUpper[i]] = -row[nonBasicAtUpper[i]];
}
} /* flip */
/************************************************************************/
void CglRedSplit2::unflip(double *row, double *tabrowrhs) {
int i;
for(i=0; i<card_nonBasicAtLower; i++) {
int locind = nonBasicAtLower[i];
if(locind < ncol) {
*tabrowrhs += row[locind] * colLower[locind];
}
}
for(i=0; i<card_nonBasicAtUpper; i++) {
int locind = nonBasicAtUpper[i];
row[locind] = -row[locind];
if(locind < ncol) {
*tabrowrhs += row[locind] * colUpper[locind];
}
}
} /* unflip */
/************************************************************************/
int CglRedSplit2::check_dynamism(double *row) {
int i, nelem = 0;
double val, max_val = 0, min_val = param.getINFINIT();
for(i=0; i<ncol; i++) {
val = fabs(row[i]);
max_val = CoinMax(max_val, val);
if(val > param.getEPS_COEFF()) {
min_val = CoinMin(min_val, val);
nelem++;
}
}
if((max_val < param.getMAXDYN() * min_val) &&
(max_val >= min_val)) {
return 1;
}
#ifdef RS2_TRACE
printf("CglRedSplit2::check_dynamism(): max_val: %6.6f min_val: %6.6f dynamism: %g\n",
max_val, min_val, max_val/min_val);
#endif
return 0;
} /* row_scale_factor */
/************************************************************************/
int CglRedSplit2::generate_packed_row(const double *lclXlp,
double *row,
int *rowind, double *rowelem,
int *card_row, double & rhs) {
int i;
double value;
int max_support = param.getMAX_SUPP_ABS() +
(static_cast<int>(ncol*param.getMAX_SUPP_REL()));
if(!check_dynamism(row)) {
#ifdef RS2_TRACE
printf("CglRedSplit2::generate_packed_row(): Cut discarded (bad dynamism)\n");
#endif
return(0);
}
*card_row = 0;
for(i=0; i<ncol; i++) {
value = row[i];
if(fabs(value) > param.getEPS_COEFF()) {
rowind[*card_row] = i;
rowelem[*card_row] = value;
(*card_row)++;
if(*card_row > max_support) {
#ifdef RS2_TRACE
printf("CglRedSplit2::generate_packed_row(): Cut discarded (too many non zero)\n");
#endif
return(0);
}
} else {
if (value > 0.0) {
rhs -= value * colLower[i];
}
else {
rhs -= value * colUpper[i];
}
}
}
value = 0;
for(i=0; i<(*card_row); i++) {
value += lclXlp[rowind[i]] * rowelem[i];
}
if(value > rhs) {
value = value-rhs;
if(value < param.getMINVIOL()) {
#ifdef RS2_TRACE
printf("CglRedSplit2::generate_packed_row(): Cut discarded: violation: %12.10f\n", value);
#endif
return(0);
}
}
return(1);
} /* generate_packed_row */
/************************************************************************/
bool CglRedSplit2::rs_are_different_vectors(const int *vect1,
const int *vect2,
const int dim) {
int i;
for(i=0; i<dim; i++) {
if(vect1[i] != vect2[i]) {
return(1);
}
}
return(0);
} /* rs_are_different_vectors */
/************************************************************************/
void CglRedSplit2::generateCuts(const OsiSolverInterface &si, OsiCuts & cs,
const CglTreeInfo info)
{
solver = const_cast<OsiSolverInterface *>(&si);
if(solver == NULL) {
printf("### WARNING: CglRedSplit2::generateCuts(): no solver available.\n");
return;
}
if(!solver->optimalBasisIsAvailable()) {
printf("### WARNING: CglRedSplit2::generateCuts(): no optimal basis available.\n");
return;
}
// Reset some members of CglRedSplit2
card_intBasicVar = 0;
card_intBasicVar_frac = 0;
card_intNonBasicVar = 0;
card_contNonBasicVar = 0;
card_nonBasicAtUpper = 0;
card_nonBasicAtLower = 0;
numRedRows = 0;
startTime = CoinCpuTime();
// Get basic problem information from solver
ncol = solver->getNumCols();
nrow = solver->getNumRows();
colLower = solver->getColLower();
colUpper = solver->getColUpper();
rowLower = solver->getRowLower();
rowUpper = solver->getRowUpper();
rowRhs = solver->getRightHandSide();
reducedCost = solver->getReducedCost();
rowPrice = solver->getRowPrice();
objective = solver->getObjCoefficients();
xlp = solver->getColSolution();
rowActivity = solver->getRowActivity();
byRow = solver->getMatrixByRow();
solver->enableFactorization();
generateCuts(&cs, param.getMaxNumCuts());
solver->disableFactorization();
} /* generateCuts */
/************************************************************************/
int CglRedSplit2::generateCuts(OsiCuts* cs, int maxNumCuts, int* lambda)
{
int i;
is_integer = new int[ncol];
compute_is_integer();
int *cstat = new int[ncol];
int *rstat = new int[nrow];
solver->getBasisStatus(cstat, rstat); // 0: free 1: basic
// 2: upper 3: lower
int *basis_index = new int[nrow]; // basis_index[i] =
// index of pivot var in row i
// (slack if number >= ncol)
#ifdef RS2_TRACETAB
rs_printvecINT("cstat", cstat, ncol);
rs_printvecINT("rstat", rstat, nrow);
#endif
solver->getBasics(basis_index);
cv_intBasicVar = new int[ncol];
cv_intBasicVar_frac = new int[ncol];
intBasicVar = new int[ncol];
intNonBasicVar = new int[ncol];
contNonBasicVar = new int[ncol+nrow];
nonBasicAtUpper = new int[ncol+nrow];
nonBasicAtLower = new int[ncol+nrow];
double dist_int;
for(i=0; i<ncol; i++) {
cv_intBasicVar[i] = 0;
cv_intBasicVar_frac[i] = 0;
switch(cstat[i]) {
case 1: // basic variable
dist_int = rs_above_integer(xlp[i]);
if(is_integer[i]){
if ((dist_int > param.getAway()) && (dist_int < 1 - param.getAway())) {
cv_intBasicVar_frac[i] = 1;
card_intBasicVar_frac++;
// intBasicVar_frac computed below,
// as order must be according to selected rows
}
card_intBasicVar++;
cv_intBasicVar[i] = 1;
}
break;
case 2: // Non basic at upper bound: must be flipped and shifted
// so that it becomes non negative with lower bound 0
// It is assumed that bounds for integer variables have been
// tightend so that non basic integer structural variables
// have integer values
nonBasicAtUpper[card_nonBasicAtUpper] = i;
card_nonBasicAtUpper++;
if(is_integer[i]) {
intNonBasicVar[card_intNonBasicVar] = i;
card_intNonBasicVar++;
}
else {
contNonBasicVar[card_contNonBasicVar] = i;
card_contNonBasicVar++;
}
break;
case 3 : // non basic at lower bound: must be shifted so that it becomes
// non negative with lower bound 0
// It is assumed that bounds for integer variables have been
// tightend so that they are integer
nonBasicAtLower[card_nonBasicAtLower] = i;
card_nonBasicAtLower++;
if(is_integer[i]) {
intNonBasicVar[card_intNonBasicVar] = i;
card_intNonBasicVar++;
}
else {
contNonBasicVar[card_contNonBasicVar] = i;
card_contNonBasicVar++;
}
break;
default: // free variable ? Don't know how to handle
printf("### ERROR: CglRedSplit2::generateCuts(): cstat[%d]: %d\n",
i, cstat[i]);
exit(1);
break;
}
}
// Use this instead of rowRhs to allow for ranges
double *effective_rhs = new double[nrow];
for(i=0; i<nrow; i++) {
// effective rhs
effective_rhs[i] = rowRhs[i];
switch(rstat[i]) {
case 1: // basic slack
break;
case 2: // non basic slack at upper; flipped and shifted
effective_rhs[i] = rowLower[i];
nonBasicAtUpper[card_nonBasicAtUpper] = ncol+i;
card_nonBasicAtUpper++;
contNonBasicVar[card_contNonBasicVar] = ncol+i;
card_contNonBasicVar++;
break;
case 3: // non basic slack at lower; shifted
effective_rhs[i] = rowUpper[i];
nonBasicAtLower[card_nonBasicAtLower] = ncol+i;
card_nonBasicAtLower++;
contNonBasicVar[card_contNonBasicVar] = ncol+i;
card_contNonBasicVar++;
break;
default:
printf("### ERROR: CglRedSlpit::generateCuts(): rstat[%d]: %d\n",
i, rstat[i]);
exit(1);
break;
}
assert (fabs(effective_rhs[i])<1.0e100);
}
#ifdef RS2_TRACE
printf("CglRedSplit2()::card_intBasicVar_frac %d %d\n",
card_intBasicVar_frac, card_contNonBasicVar);
#endif
if((card_contNonBasicVar == 0) || (card_intBasicVar_frac == 0)) {
delete[] cstat;
delete[] rstat;
delete[] basis_index;
delete[] cv_intBasicVar;
delete[] cv_intBasicVar_frac;
delete[] intBasicVar;
delete[] intNonBasicVar;
delete[] contNonBasicVar;
delete[] nonBasicAtUpper;
delete[] nonBasicAtLower;
delete[] is_integer;
delete[] effective_rhs;
return 0; // no cuts can be generated
}
double *z = new double[ncol]; // workspace to get row of the tableau
double *slack = new double[nrow]; // workspace to get row of the tableau
#ifdef RS2_TRACETAB
printOptTab(solver);
#endif
mTab = card_intBasicVar;
nTab = card_contNonBasicVar;
rhsTab = new double[mTab];
cv_fracRowsTab = new int[mTab];
memset(cv_fracRowsTab, 0, mTab*sizeof(int));
int card_rowTab = 0;
rs_allocmatDBL(&contNonBasicTab, mTab, card_contNonBasicVar);
rs_allocmatDBL(&workNonBasicTab, mTab, card_contNonBasicVar);
rs_allocmatDBL(&intNonBasicTab, mTab, card_intNonBasicVar);
norm = new double[mTab];
intBasicVar_frac = new int[ncol];
// position of each integer basic variable in the simplex tableau
int* origRow = new int[nrow];
card_intBasicVar = 0; // recompute in pivot order
card_intBasicVar_frac = 0;
for(i=0; i<nrow; i++) {
int ind_row = i;
if(basis_index[ind_row] >= ncol ||
cv_intBasicVar[basis_index[ind_row]] != 1) {
// If basic variable is not an integer variable, skip
continue;
}
origRow[card_intBasicVar] = ind_row;
// row used in generation
intBasicVar[card_intBasicVar] = basis_index[ind_row];
if (cv_intBasicVar_frac[basis_index[ind_row]] == 1){
// row is fractional
intBasicVar_frac[card_intBasicVar_frac] = basis_index[ind_row];
cv_fracRowsTab[card_intBasicVar] = 1;
}
// obtain row of simplex tableau
solver->getBInvARow(ind_row, z, slack);
rhsTab[card_rowTab] = xlp[basis_index[ind_row]];
int nonzeroes = 0;
int ii;
for(ii=0; ii<card_contNonBasicVar; ii++) {
int locind = contNonBasicVar[ii];
if(locind < ncol) {
contNonBasicTab[card_rowTab][ii] = z[locind];
}
else {
contNonBasicTab[card_rowTab][ii] = slack[locind - ncol];
}
if(fabs(contNonBasicTab[card_rowTab][ii]) > TINY)
nonzeroes++;
}
for(ii=0; ii<card_intNonBasicVar; ii++) {
int locind = intNonBasicVar[ii];
if(locind < ncol) {
intNonBasicTab[card_rowTab][ii] = z[locind];
}
else {
printf("### ERROR: CglRedSplit2::generateCuts(): integer slack unexpected\n");
exit(1);
}
if(fabs(intNonBasicTab[card_rowTab][ii]) > TINY)
nonzeroes++;
}
#if RS_FAST_INT > 0 || RS_FAST_CONT > 0 || RS_FAST_WORK > 0
if(nonzeroes < param.getMaxNonzeroesTab()) {
#endif
// The number of nonzeroes is small enough, we keep the row
card_rowTab++;
card_intBasicVar++;
if (cv_intBasicVar_frac[basis_index[ind_row]] == 1)
card_intBasicVar_frac++;
#if RS_FAST_INT > 0 || RS_FAST_CONT > 0 || RS_FAST_WORK > 0
}
#endif
}
#if RS_FAST_INT == 0 && RS_FAST_CONT == 0 && RS_FAST_WORK == 0
rs_allocmatINT(&pi_mat, mTab, mTab);
#else
// Allow for sparse info (and for work)
rs_allocmatINT(&pi_mat, mTab, mTab+
card_intNonBasicVar+2*card_contNonBasicVar+3);
for (int i=0;i<mTab;i++) {
int * info = pi_mat[i];
int * intSparse=info+mTab;
int * contSparse=intSparse+1+card_intNonBasicVar;
int n;
const double * values;
n=0;
values=intNonBasicTab[i];
for (int j=0;j<card_intNonBasicVar;j++) {
if (fabs(values[j]) > TINY)
intSparse[++n]=j;
}
intSparse[0]=n;
n=0;
values=contNonBasicTab[i];
for (int j=0;j<card_contNonBasicVar;j++) {
if (fabs(values[j]))
contSparse[++n]=j;
}
contSparse[0]=n;
}
#endif
#ifdef RS2_TRACE
printf("intBasicVar:\n");
for(i=0; i<card_intBasicVar; i++) {
printf("%d ", intBasicVar[i]);
}
printf("\n");
printf("intBasicVar_frac:\n");
for(i=0; i<card_intBasicVar_frac; i++) {
printf("%d ", intBasicVar_frac[i]);
}
printf("\n");
printf("intNonBasicVar:\n");
for(i=0; i<card_intNonBasicVar; i++) {
printf("%d ", intNonBasicVar[i]);
}
printf("\n");
printf("contNonBasicVar:\n");
for(i=0; i<card_contNonBasicVar; i++) {
printf("%d ", contNonBasicVar[i]);
}
printf("\n");
#endif
int card_row;
double *row = new double[ncol+nrow];
int *rowind = new int[ncol];
double *rowelem = new double[ncol];
const double *elements = byRow->getElements();
const int *rowStart = byRow->getVectorStarts();
const int *indices = byRow->getIndices();
const int *rowLength = byRow->getVectorLengths();
CglRedSplit2Param::ColumnSelectionStrategy columnSelection;
CglRedSplit2Param::RowSelectionStrategy rowSelection;
int numRows;
int maxNumComputedCuts = param.getMaxNumComputedCuts();
int* bufflambda = lambda;
OsiCuts* buffcs = cs;
// quality of cuts
struct sortElement* quality = NULL;
if (maxNumCuts < maxNumComputedCuts){
// user wants the generator to generate several cuts and select only best
if (lambda != NULL){
bufflambda = new int[maxNumComputedCuts*nrow];
}
if (cs != NULL){
buffcs = new OsiCuts();
}
quality = new struct sortElement[maxNumComputedCuts];
}
int initNumCuts = 0;
if (buffcs != NULL){
initNumCuts = buffcs->sizeRowCuts();
}
// reset vector of lambdas
if (bufflambda != NULL){
if (maxNumCuts < maxNumComputedCuts){
memset(bufflambda, 0, maxNumComputedCuts*nrow*sizeof(int));
}
else{
memset(bufflambda, 0, maxNumCuts*nrow*sizeof(int));
}
}
int numCuts = 0;
std::vector<CglRedSplit2Param::ColumnSelectionStrategy> listColSel =
param.getColumnSelectionStrategy();
std::vector<CglRedSplit2Param::RowSelectionStrategy> listRowSel =
param.getRowSelectionStrategy();
std::vector<int> listNumRows = param.getNumRowsReduction();
#if 0
double work= listColSel.size()*listNumRows.size()*listRowSel.size()*nDiag;
work *= mTab*(card_intNonBasicVar+card_contNonBasicVar);
printf("listColSel %d listNumRows %d listRowSel %d mTab %d intNon %d contnon %d - %d generate_rows and dense ops %g\n",
listColSel.size(),listNumRows.size(),listRowSel.size(),mTab,
card_intNonBasicVar,card_contNonBasicVar,
listColSel.size()*listNumRows.size()*listRowSel.size()*nDiag,work);
#endif
for (unsigned int coliter = 0; coliter < listColSel.size(); ++coliter){
if (!checkTime() || numCuts >= maxNumComputedCuts){
break;
}
columnSelection = listColSel[coliter];
#ifdef RS2_TRACE
printf("Filling workNonBasicTab with column strategy %d\n", columnSelection);
#endif
// select the columns
fill_workNonBasicTab(columnSelection);
for (unsigned int nriter = 0; nriter < listNumRows.size(); ++nriter){
if (!checkTime() || numCuts >= maxNumComputedCuts){
break;
}
numRows = listNumRows[nriter];
#ifdef RS2_TRACE
printf("Applying reduction algorithm with %d rows\n", numRows);
#endif
for (unsigned int rowiter = 0; rowiter < listRowSel.size(); ++rowiter){
if (!checkTime() || numCuts >= maxNumComputedCuts){
break;
}
rowSelection = listRowSel[rowiter];
#ifdef RS2_TRACE
printf("Reducing norms with row strategy %d\n", rowSelection);
#endif
// new iteration: reinitialize pi_mat
for(i=0; i<mTab; i++) {
memset(pi_mat[i], 0, mTab*sizeof(int));
if (param.getSkipGomory() == false){
pi_mat[i][i] = 1;
}
}
reduce_workNonBasicTab(numRows, rowSelection, mTab);
// now generate cuts
for(i=0; i<mTab && numCuts < maxNumComputedCuts; i++) {
if (pi_mat[i][i] != 0){
card_row = 0;
generate_row(i, row);
flip(row);
// RHS of equalities after flipping/translating non basic variables
// is given by the current LP solution (re-ordered according to
// basis_index), i.e. what is now in rhsTab.
// RHS of row i of (pi_mat * contNonBasicTab) is then simply
// pi_mat[i] * rhsTab
double tabrowrhs = rs_dotProd(pi_mat[i], rhsTab, mTab);
int got_one = 0;
got_one = generate_cgcut(row, &tabrowrhs);
if (!got_one && param.getSkipGomory() == false){
// Cut improvement failed and we want to generate Gomory
// cuts from single rows; so, revert back to one row and try
// the corresponding Gomory cut. If we still fail, we give up
memset(pi_mat[i], 0, mTab*sizeof(int));
pi_mat[i][i] = 1;
card_row = 0;
generate_row(i, row);
flip(row);
tabrowrhs = rhsTab[i];
got_one = generate_cgcut(row, &tabrowrhs);
}
if (got_one){
bool skip = false;
if (bufflambda != NULL){
for (int k = 0; k < mTab; ++k){
bufflambda[numCuts*nrow + origRow[k]] = pi_mat[i][k];
}
for (int k = 0; k < numCuts; ++k){
if (!rs_are_different_vectors(bufflambda + k*nrow,
bufflambda + numCuts*nrow,
nrow)){
// Cut was previously generated; skip
skip = true;
}
}
}
if (!skip){
// new multipliers found, check if the generated cut is good
unflip(row, &tabrowrhs);
eliminate_slacks(row, elements, rowStart, indices,
rowLength, effective_rhs, &tabrowrhs);
if(generate_packed_row(xlp, row, rowind, rowelem, &card_row,
tabrowrhs)) {
// the cut is good, multipliers should be saved.
// record quality of cut if needed
if (quality){
quality[numCuts].index = numCuts;
quality[numCuts].cost = tabrowrhs;
double rnorm = 0.0;
for (int ii = 0; ii < card_row; ++ii){
rnorm += rowelem[ii]*rowelem[ii];
quality[numCuts].cost -= rowelem[ii]*xlp[rowind[ii]];
}
quality[numCuts].cost /= rnorm;
}
// now, add cut to the collection if needed
if (cs != NULL){
OsiRowCut rc;
rc.setRow(card_row, rowind, rowelem);
rc.setLb(-param.getINFINIT());
double adjust = param.getEPS_RELAX_ABS();
if(param.getEPS_RELAX_REL() > 0.0) {
adjust += fabs(tabrowrhs) * param.getEPS_RELAX_REL();
}
rc.setUb(tabrowrhs + adjust);
// relax the constraint slightly
buffcs->insertIfNotDuplicate(rc, CoinAbsFltEq(param.getEPS()));
numCuts = buffcs->sizeRowCuts() - initNumCuts;
}
else{
numCuts++;
}
}
}
}
}
}
}
}
}
if (maxNumCuts < maxNumComputedCuts){
// select best cuts and copy them back
if (numCuts > maxNumCuts){
qsort(quality, numCuts, sizeof(struct sortElement), rs2_compareElements);
}
// also delete temp data
if (buffcs){
for (int i = 0; i < numCuts && i < maxNumCuts; ++i){
cs->insertIfNotDuplicate(buffcs->rowCut(quality[i].index),
CoinAbsFltEq(param.getEPS_COEFF()));
}
delete buffcs;
}
if (bufflambda){
for (int i = 0; i < numCuts && i < maxNumCuts; ++i){
memcpy(lambda + (i*nrow), bufflambda + (quality[i].index*nrow),
sizeof(int)*nrow);
}
delete[] bufflambda;
}
if (numCuts > maxNumCuts){
numCuts = maxNumCuts;
}
}
delete[] cstat;
delete[] rstat;
delete[] basis_index;
delete[] slack;
delete[] z;
delete[] effective_rhs;
delete[] row;
delete[] rowind;
delete[] rowelem;
delete[] origRow;
delete[] cv_intBasicVar;
delete[] cv_intBasicVar_frac;
delete[] cv_fracRowsTab;
delete[] intBasicVar;
delete[] intBasicVar_frac;
delete[] intNonBasicVar;
delete[] contNonBasicVar;
delete[] nonBasicAtUpper;
delete[] nonBasicAtLower;
delete[] is_integer;
rs_deallocmatDBL(&contNonBasicTab, mTab);
rs_deallocmatDBL(&workNonBasicTab, mTab);
rs_deallocmatDBL(&intNonBasicTab, mTab);
rs_deallocmatINT(&pi_mat, mTab);
delete[] rhsTab;
delete[] norm;
return numCuts;
} /* generateCuts */
/***********************************************************************/
void CglRedSplit2::setParam(const CglRedSplit2Param &source) {
param = source;
} /* setParam */
/***********************************************************************/
void CglRedSplit2::compute_is_integer() {
int i;
for(i=0; i<ncol; i++) {
if(solver->isInteger(i)) {
is_integer[i] = 1;
}
else {
if((colUpper[i] - colLower[i] < param.getEPS()) &&
(rs_above_integer(colUpper[i]) < param.getEPS())) {
// continuous variable fixed to an integer value
is_integer[i] = 1;
}
else {
is_integer[i] = 0;
}
}
}
} /* compute_is_integer */
/***********************************************************************/
void CglRedSplit2::print() const
{
rs_printvecINT("intBasicVar_frac", intBasicVar_frac, card_intBasicVar_frac);
rs_printmatINT("pi_mat", pi_mat, card_intBasicVar_frac,
card_intBasicVar_frac);
rs_printvecINT("intNonBasicVar", intNonBasicVar, card_intNonBasicVar);
rs_printmatDBL("intNonBasicTab", intNonBasicTab, card_intBasicVar_frac,
card_intNonBasicVar);
rs_printvecINT("contNonBasicVar", contNonBasicVar, card_contNonBasicVar);
rs_printmatDBL("contNonBasicTab", contNonBasicTab, card_intBasicVar_frac,
card_contNonBasicVar);
rs_printvecINT("nonBasicAtLower", nonBasicAtLower, card_nonBasicAtLower);
rs_printvecINT("nonBasicAtUpper", nonBasicAtUpper, card_nonBasicAtUpper);
} /* print */
/***********************************************************************/
void CglRedSplit2::printOptTab(OsiSolverInterface *lclSolver) const
{
int i;
int *cstat = new int[ncol];
int *rstat = new int[nrow];
lclSolver->getBasisStatus(cstat, rstat); // 0: free 1: basic
// 2: upper 3: lower
int *basis_index = new int[nrow]; // basis_index[i] =
// index of pivot var in row i
// (slack if number >= ncol)
lclSolver->getBasics(basis_index);
double *z = new double[ncol]; // workspace to get row of the tableau
double *slack = new double[nrow]; // workspace to get row of the tableau
double *slack_val = new double[nrow];
for(i=0; i<nrow; i++) {
slack_val[i] = rowRhs[i] - rowActivity[i];
}
const double *rc = lclSolver->getReducedCost();
const double *dual = lclSolver->getRowPrice();
const double *solution = lclSolver->getColSolution();
rs_printvecINT("cstat", cstat, ncol);
rs_printvecINT("rstat", rstat, nrow);
rs_printvecINT("basis_index", basis_index, nrow);
rs_printvecDBL("solution", solution, ncol);
rs_printvecDBL("slack_val", slack_val, nrow);
rs_printvecDBL("reduced_costs", rc, ncol);
rs_printvecDBL("dual solution", dual, nrow);
printf("Optimal Tableau:\n");
for(i=0; i<nrow; i++) {
lclSolver->getBInvARow(i, z, slack);
int ii;
for(ii=0; ii<ncol; ii++) {
printf("%5.2f ", z[ii]);
}
printf(" | ");
for(ii=0; ii<nrow; ii++) {
printf("%5.2f ", slack[ii]);
}
printf(" | ");
if(basis_index[i] < ncol) {
printf("%5.2f ", solution[basis_index[i]]);
}
else {
printf("%5.2f ", slack_val[basis_index[i]-ncol]);
}
printf("\n");
}
int ii;
for(ii=0; ii<7*(ncol+nrow+1); ii++) {
printf("-");
}
printf("\n");
for(ii=0; ii<ncol; ii++) {
printf("%5.2f ", rc[ii]);
}
printf(" | ");
for(ii=0; ii<nrow; ii++) {
printf("%5.2f ", -dual[ii]);
}
printf(" | ");
printf("%5.2f\n", -lclSolver->getObjValue());
delete[] cstat;
delete[] rstat;
delete[] basis_index;
delete[] slack;
delete[] z;
delete[] slack_val;
} /* printOptTab */
/*********************************************************************/
CglRedSplit2::CglRedSplit2() :
CglCutGenerator(),
nrow(0),
ncol(0),
card_intBasicVar(0),
card_intBasicVar_frac(0),
card_intNonBasicVar(0),
card_contNonBasicVar(0),
card_nonBasicAtUpper(0),
card_nonBasicAtLower(0),
cv_intBasicVar(0),
cv_intBasicVar_frac(0),
cv_fracRowsTab(0),
intBasicVar(0),
intBasicVar_frac(0),
intNonBasicVar(0),
contNonBasicVar(0),
nonBasicAtUpper(0),
nonBasicAtLower(0),
mTab(0),
nTab(0),
pi_mat(0),
contNonBasicTab(0),
intNonBasicTab(0),
rhsTab(0)
{
}
/*********************************************************************/
CglRedSplit2::CglRedSplit2(const CglRedSplit2Param &RS_param) :
CglCutGenerator(),
nrow(0),
ncol(0),
card_intBasicVar(0),
card_intBasicVar_frac(0),
card_intNonBasicVar(0),
card_contNonBasicVar(0),
card_nonBasicAtUpper(0),
card_nonBasicAtLower(0),
cv_intBasicVar(0),
cv_intBasicVar_frac(0),
cv_fracRowsTab(0),
intBasicVar(0),
intBasicVar_frac(0),
intNonBasicVar(0),
contNonBasicVar(0),
nonBasicAtUpper(0),
nonBasicAtLower(0),
mTab(0),
nTab(0),
pi_mat(0),
contNonBasicTab(0),
intNonBasicTab(0),
rhsTab(0)
{
param = RS_param;
}
/*********************************************************************/
CglRedSplit2::CglRedSplit2 (const CglRedSplit2 & source) :
CglCutGenerator(source),
param(source.param),
nrow(0),
ncol(0),
card_intBasicVar(0),
card_intBasicVar_frac(0),
card_intNonBasicVar(0),
card_contNonBasicVar(0),
card_nonBasicAtUpper(0),
card_nonBasicAtLower(0),
cv_intBasicVar(NULL),
cv_intBasicVar_frac(NULL),
cv_fracRowsTab(NULL),
intBasicVar(NULL),
intBasicVar_frac(NULL),
intNonBasicVar(NULL),
contNonBasicVar(NULL),
nonBasicAtUpper(NULL),
nonBasicAtLower(NULL),
mTab(0),
nTab(0),
pi_mat(NULL),
contNonBasicTab(NULL),
intNonBasicTab(NULL),
rhsTab(NULL)
{
}
/*********************************************************************/
CglCutGenerator *
CglRedSplit2::clone() const
{
return new CglRedSplit2(*this);
}
/*********************************************************************/
CglRedSplit2::~CglRedSplit2 ()
{}
/*********************************************************************/
CglRedSplit2 &
CglRedSplit2::operator=(const CglRedSplit2 &source)
{
if (this != &source) {
CglCutGenerator::operator=(source);
param = source.param;
}
return *this;
}
/*********************************************************************/
// Returns true if needs optimal basis to do cuts
bool
CglRedSplit2::needsOptimalBasis() const
{
return true;
}
/************************************************************************/
int CglRedSplit2::generateMultipliers(const OsiSolverInterface &si,
int* lambda,
int maxNumMultipliers,
int* basicVariables,
OsiCuts* cs)
{
solver = const_cast<OsiSolverInterface *>(&si);
if(solver == NULL) {
printf("### WARNING: CglRedSplit2::generateCuts(): no solver available.\n");
return 0;
}
if(!solver->optimalBasisIsAvailable()) {
printf("### WARNING: CglRedSplit2::generateCuts(): no optimal basis available.\n");
return 0;
}
// Reset some members of CglRedSplit2
card_intBasicVar = 0;
card_intBasicVar_frac = 0;
card_intNonBasicVar = 0;
card_contNonBasicVar = 0;
card_nonBasicAtUpper = 0;
card_nonBasicAtLower = 0;
numRedRows = 0;
startTime = CoinCpuTime();
// Get basic problem information from solver
ncol = solver->getNumCols();
nrow = solver->getNumRows();
colLower = solver->getColLower();
colUpper = solver->getColUpper();
rowLower = solver->getRowLower();
rowUpper = solver->getRowUpper();
rowRhs = solver->getRightHandSide();
reducedCost = solver->getReducedCost();
rowPrice = solver->getRowPrice();
objective = solver->getObjCoefficients();
xlp = solver->getColSolution();
rowActivity = solver->getRowActivity();
byRow = solver->getMatrixByRow();
solver->enableFactorization();
if (basicVariables != NULL){
solver->getBasics(basicVariables);
}
int numMultipliers = generateCuts(cs, maxNumMultipliers, lambda);
solver->disableFactorization();
return numMultipliers;
} /* generateMultipliers */
/***********************************************************************/
int CglRedSplit2::tiltLandPcut(const OsiSolverInterface* si,
double* landprow, double landprhs,
int rownumber, const double* xbar,
const int* newnonbasics, OsiRowCut* cs,
int* lambda){
solver = const_cast<OsiSolverInterface*>(si);
if(solver == NULL) {
printf("### WARNING: CglRedSplit2::tiltLandPcut(): no solver available.\n");
return 0;
}
// Reset some members of CglRedSplit2
card_intBasicVar = 0;
card_intBasicVar_frac = 0;
card_intNonBasicVar = 0;
card_contNonBasicVar = 0;
card_nonBasicAtUpper = 0;
card_nonBasicAtLower = 0;
numRedRows = 0;
startTime = CoinCpuTime();
// Get basic problem information from solver
ncol = solver->getNumCols();
nrow = solver->getNumRows();
colLower = solver->getColLower();
colUpper = solver->getColUpper();
rowLower = solver->getRowLower();
rowUpper = solver->getRowUpper();
rowRhs = solver->getRightHandSide();
reducedCost = solver->getReducedCost();
rowPrice = solver->getRowPrice();
objective = solver->getObjCoefficients();
xlp = solver->getColSolution();
rowActivity = solver->getRowActivity();
byRow = solver->getMatrixByRow();
int i;
is_integer = new int[ncol];
compute_is_integer();
int *cstat = new int[ncol];
int *rstat = new int[nrow];
solver->getBasisStatus(cstat, rstat); // 0: free 1: basic
// 2: upper 3: lower
int *basis_index = new int[nrow]; // basis_index[i] =
// index of pivot var in row i
// (slack if number >= ncol)
#ifdef RS2_TRACETAB
rs_printvecINT("cstat", cstat, ncol);
rs_printvecINT("rstat", rstat, nrow);
#endif
solver->getBasics(basis_index);
cv_intBasicVar = new int[ncol];
cv_intBasicVar_frac = new int[ncol];
intBasicVar = new int[ncol];
intNonBasicVar = new int[ncol];
contNonBasicVar = new int[ncol+nrow];
nonBasicAtUpper = new int[ncol+nrow];
nonBasicAtLower = new int[ncol+nrow];
double dist_int;
for(i=0; i<ncol; i++) {
cv_intBasicVar[i] = 0;
cv_intBasicVar_frac[i] = 0;
switch(cstat[i]) {
case 1: // basic variable
dist_int = rs_above_integer(xlp[i]);
if(is_integer[i]){
if ((dist_int > param.getAway()) && (dist_int < 1 - param.getAway())) {
cv_intBasicVar_frac[i] = 1;
card_intBasicVar_frac++;
// intBasicVar_frac computed below,
// as order must be according to selected rows
}
card_intBasicVar++;
cv_intBasicVar[i] = 1;
}
break;
case 2: // Non basic at upper bound: must be flipped and shifted
// so that it becomes non negative with lower bound 0
// It is assumed that bounds for integer variables have been
// tightend so that non basic integer structural variables
// have integer values
nonBasicAtUpper[card_nonBasicAtUpper] = i;
card_nonBasicAtUpper++;
if(is_integer[i]) {
intNonBasicVar[card_intNonBasicVar] = i;
card_intNonBasicVar++;
}
else {
contNonBasicVar[card_contNonBasicVar] = i;
card_contNonBasicVar++;
}
break;
case 3 : // non basic at lower bound: must be shifted so that it becomes
// non negative with lower bound 0
// It is assumed that bounds for integer variables have been
// tightend so that they are integer
nonBasicAtLower[card_nonBasicAtLower] = i;
card_nonBasicAtLower++;
if(is_integer[i]) {
intNonBasicVar[card_intNonBasicVar] = i;
card_intNonBasicVar++;
}
else {
contNonBasicVar[card_contNonBasicVar] = i;
card_contNonBasicVar++;
}
break;
default: // free variable ? Don't know how to handle
printf("### ERROR: CglRedSplit2::generateCuts(): cstat[%d]: %d\n",
i, cstat[i]);
exit(1);
break;
}
}
// Use this instead of rowRhs to allow for ranges
double *effective_rhs = new double[nrow];
for(i=0; i<nrow; i++) {
// effective rhs
effective_rhs[i] = rowRhs[i];
switch(rstat[i]) {
case 1: // basic slack
break;
case 2: // non basic slack at upper; flipped and shifted
effective_rhs[i] = rowLower[i];
nonBasicAtUpper[card_nonBasicAtUpper] = ncol+i;
card_nonBasicAtUpper++;
contNonBasicVar[card_contNonBasicVar] = ncol+i;
card_contNonBasicVar++;
break;
case 3: // non basic slack at lower; shifted
effective_rhs[i] = rowUpper[i];
nonBasicAtLower[card_nonBasicAtLower] = ncol+i;
card_nonBasicAtLower++;
contNonBasicVar[card_contNonBasicVar] = ncol+i;
card_contNonBasicVar++;
break;
default:
printf("### ERROR: CglRedSlpit::generateCuts(): rstat[%d]: %d\n",
i, rstat[i]);
exit(1);
break;
}
assert (fabs(effective_rhs[i])<1.0e100);
}
#ifdef RS2_TRACE
printf("CglRedSplit2()::card_intBasicVar_frac %d %d\n",
card_intBasicVar_frac, card_contNonBasicVar);
#endif
if((card_contNonBasicVar == 0) || (card_intBasicVar == 0)) {
delete[] cstat;
delete[] rstat;
delete[] basis_index;
delete[] cv_intBasicVar;
delete[] cv_intBasicVar_frac;
delete[] intBasicVar;
delete[] intNonBasicVar;
delete[] contNonBasicVar;
delete[] nonBasicAtUpper;
delete[] nonBasicAtLower;
delete[] is_integer;
delete[] effective_rhs;
printf("No vars to generate cut\n");
return 0; // no cuts can be generated
}
double *z = new double[ncol]; // workspace to get row of the tableau
double *slack = new double[nrow]; // workspace to get row of the tableau
#ifdef RS2_TRACETAB
printOptTab(solver);
#endif
if (rownumber >= 0){
// this means that the source row is in the given simplex tableau
mTab = card_intBasicVar;
}
else{
// source row is an extra row, so we need extra space
mTab = card_intBasicVar + 1;
}
nTab = card_contNonBasicVar;
rhsTab = new double[mTab];
cv_fracRowsTab = new int[mTab];
memset(cv_fracRowsTab, 0, mTab*sizeof(int));
int card_rowTab = 0;
int extraColumns = 0;
while (newnonbasics[extraColumns] >= 0){
extraColumns++;
}
rs_allocmatDBL(&contNonBasicTab, mTab, card_contNonBasicVar);
rs_allocmatDBL(&workNonBasicTab, mTab, card_contNonBasicVar + extraColumns);
rs_allocmatDBL(&intNonBasicTab, mTab, card_intNonBasicVar);
norm = new double[mTab];
intBasicVar_frac = new int[ncol];
card_intBasicVar = 0; // recompute in pivot order
card_intBasicVar_frac = 0;
// write L&P row as first row
rhsTab[card_rowTab] = landprhs;
// flip row: L&P flips the nonbasic at upper bound, but we don't want this
for (int ii=0; ii<card_nonBasicAtUpper; ii++){
landprow[nonBasicAtUpper[ii]] = -landprow[nonBasicAtUpper[ii]];
}
for (int ii=0; ii<card_contNonBasicVar; ii++) {
int locind = contNonBasicVar[ii];
contNonBasicTab[card_rowTab][ii] = landprow[locind];
}
for (int ii=0; ii<card_intNonBasicVar; ii++) {
int locind = intNonBasicVar[ii];
intNonBasicTab[card_rowTab][ii] = landprow[locind];
}
// flip back to restore L&P row as it was given
for (int ii=0; ii<card_nonBasicAtUpper; ii++){
landprow[nonBasicAtUpper[ii]] = -landprow[nonBasicAtUpper[ii]];
}
card_rowTab++;
for(i=0; i<nrow; i++) {
int ind_row = i;
if(basis_index[ind_row] >= ncol || ind_row == rownumber) {
continue;
}
if(cv_intBasicVar[basis_index[ind_row]] == 1) {
// row used in generation
intBasicVar[card_intBasicVar] = basis_index[ind_row];
if (cv_intBasicVar_frac[basis_index[ind_row]] == 1){
// row is fractional
intBasicVar_frac[card_intBasicVar_frac] = basis_index[ind_row];
card_intBasicVar_frac++;
cv_fracRowsTab[card_intBasicVar] = 1;
}
card_intBasicVar++;
int ii;
rhsTab[card_rowTab] = xlp[basis_index[ind_row]];
solver->getBInvARow(ind_row, z, slack);
for(ii=0; ii<card_contNonBasicVar; ii++) {
int locind = contNonBasicVar[ii];
if(locind < ncol) {
contNonBasicTab[card_rowTab][ii] = z[locind];
}
else {
contNonBasicTab[card_rowTab][ii] = slack[locind - ncol];
}
}
for(ii=0; ii<card_intNonBasicVar; ii++) {
int locind = intNonBasicVar[ii];
if(locind < ncol) {
intNonBasicTab[card_rowTab][ii] = z[locind];
}
else {
printf("### ERROR: CglRedSplit2::generateCuts(): integer slack unexpected\n");
exit(1);
}
}
card_rowTab++;
}
}
int pi_mat_rows = 1;
rs_allocmatINT(&pi_mat, pi_mat_rows, mTab);
#ifdef RS2_TRACE
printf("intBasicVar:\n");
for(i=0; i<card_intBasicVar; i++) {
printf("%d ", intBasicVar[i]);
}
printf("\n");
printf("intBasicVar_frac:\n");
for(i=0; i<card_intBasicVar_frac; i++) {
printf("%d ", intBasicVar_frac[i]);
}
printf("\n");
printf("intNonBasicVar:\n");
for(i=0; i<card_intNonBasicVar; i++) {
printf("%d ", intNonBasicVar[i]);
}
printf("\n");
printf("contNonBasicVar:\n");
for(i=0; i<card_contNonBasicVar; i++) {
printf("%d ", contNonBasicVar[i]);
}
printf("\n");
#endif
int card_row;
double *row = new double[ncol+nrow];
int *rowind = new int[ncol];
double *rowelem = new double[ncol];
const double *elements = byRow->getElements();
const int *rowStart = byRow->getVectorStarts();
const int *indices = byRow->getIndices();
const int *rowLength = byRow->getVectorLengths();
CglRedSplit2Param::ColumnSelectionStrategy columnSelection;
CglRedSplit2Param::RowSelectionStrategy rowSelection;
int numRows;
int generatedCuts = 0;
std::vector<CglRedSplit2Param::ColumnSelectionStrategy> listColSel
= param.getColumnSelectionStrategyLAP();
std::vector<CglRedSplit2Param::RowSelectionStrategy> listRowSel
= param.getRowSelectionStrategyLAP();
std::vector<int> listNumRows = param.getNumRowsReductionLAP();
for (unsigned int coliter = 0; coliter < listColSel.size(); ++coliter){
if (!checkTime()){
break;
}
columnSelection = listColSel[coliter];
#ifdef RS2_TRACE
printf("Filling workNonBasicTab with column strategy %d\n", columnSelection);
#endif
nTab = 0;
// select the columns
if (columnSelection != CglRedSplit2Param::CS_LAP_NONBASICS){
fill_workNonBasicTab(columnSelection, newnonbasics);
}
fill_workNonBasicTab(newnonbasics, xbar,
param.getColumnScalingStrategyLAP());
for (unsigned int nriter = 0; nriter < listNumRows.size(); ++nriter){
if (!checkTime()){
break;
}
numRows = listNumRows[nriter];
#ifdef RS2_TRACE
printf("Applying reduction algorithm with %d rows\n", numRows);
#endif
for (unsigned int rowiter = 0; rowiter < listRowSel.size(); ++rowiter){
if (!checkTime()){
break;
}
rowSelection = listRowSel[rowiter];
#ifdef RS2_TRACE
printf("Reducing norms with row strategy %d\n", rowSelection);
#endif
// new iteration: reinitialize pi_mat
for(i=0; i<pi_mat_rows; i++) {
memset(pi_mat[i], 0, mTab*sizeof(int));
}
// reduce the coefficients
reduce_workNonBasicTab(numRows, rowSelection, 1);
// now generate cuts
for(i=0; i<pi_mat_rows; i++) {
if (pi_mat[i][i] == 0){
// means we did not generate a cut from this row
continue;
}
card_row = 0;
generate_row(i, row);
flip(row);
// compute RHS
double tabrowrhs = rs_dotProd(pi_mat[i], rhsTab, mTab);
int got_one = 0;
got_one = generate_cgcut(row, &tabrowrhs);
if (got_one){
// new multipliers found, check if the generated cut is good
unflip(row, &tabrowrhs);
eliminate_slacks(row, elements, rowStart, indices,
rowLength, effective_rhs, &tabrowrhs);
if(generate_packed_row(xbar, row, rowind, rowelem, &card_row,
tabrowrhs)) {
// the cut is good, multipliers should be saved;
// first of all, add cut to the collection if needed
cs->setRow(card_row, rowind, rowelem);
cs->setLb(-param.getINFINIT());
double adjust = param.getEPS_RELAX_ABS();
if(param.getEPS_RELAX_REL() > 0.0) {
adjust += fabs(tabrowrhs) * param.getEPS_RELAX_REL();
}
cs->setUb(tabrowrhs + adjust);
if (lambda){
// Modify the initial disjunction by adding the new
// coefficients
for (int k = 1; k < mTab; ++k){
lambda[intBasicVar[k-1]] += pi_mat[i][k];
}
}
generatedCuts++;
}
}
}
}
}
}
delete[] cstat;
delete[] rstat;
delete[] basis_index;
delete[] slack;
delete[] z;
delete[] effective_rhs;
delete[] row;
delete[] rowind;
delete[] rowelem;
delete[] cv_intBasicVar_frac;
delete[] cv_fracRowsTab;
delete[] intBasicVar;
delete[] intBasicVar_frac;
delete[] intNonBasicVar;
delete[] contNonBasicVar;
delete[] nonBasicAtUpper;
delete[] nonBasicAtLower;
delete[] is_integer;
rs_deallocmatDBL(&contNonBasicTab, mTab);
rs_deallocmatDBL(&workNonBasicTab, mTab);
rs_deallocmatDBL(&intNonBasicTab, mTab);
rs_deallocmatINT(&pi_mat, pi_mat_rows);
delete[] rhsTab;
delete[] norm;
return generatedCuts;
} /* tiltLandPcut */
/***************************************************************************/
void CglRedSplit2::fill_workNonBasicTab(const int* newnonbasics,
const double* xbar,
CglRedSplit2Param::ColumnScalingStrategy scaling){
// See CglRedSplit2Param.hpp for a description of the strategies
int i = 0;
int currvar, colpos;
#ifdef RS2_TRACE
while (newnonbasics[i] >= 0){
printf("Newnonbasic[%d]: %d\n", i, newnonbasics[i]);
i++;
}
i = 0;
#endif
// In workNonBasicTab, we write the new nonbasic columns given by the
// Lift & Project algorithm, scaled by the value of the fractional solution.
while (newnonbasics[i] >= 0){
currvar = newnonbasics[i];
if (currvar < ncol && solver->isInteger(currvar)){
for (colpos = 0; colpos < card_intNonBasicVar; ++colpos){
if (intNonBasicVar[colpos] == currvar){
#ifdef RS2_TRACE
printf("Found int nonbasic at pos %d: %d %d\n", colpos, intNonBasicVar[colpos], currvar);
#endif
break;
}
}
double factor = 1.0;
if (scaling == CglRedSplit2Param::SC_LINEAR){
if (fabs(xbar[currvar]) > factor){
factor = fabs(xbar[currvar]);
}
}
else if (scaling == CglRedSplit2Param::SC_LINEAR_BOUNDED){
if (fabs(xbar[currvar]) > factor){
factor = fabs(xbar[currvar]);
}
if (factor < param.getColumnScalingBoundLAP()){
factor = param.getColumnScalingBoundLAP();
}
}
else if (scaling == CglRedSplit2Param::SC_LOG_BOUNDED){
if (log(fabs(xbar[currvar])) > factor){
factor = log(fabs(xbar[currvar]));
}
if (factor < param.getColumnScalingBoundLAP()){
factor = param.getColumnScalingBoundLAP();
}
}
else if (scaling == CglRedSplit2Param::SC_UNIFORM){
factor = param.getColumnScalingBoundLAP();
}
else if (scaling == CglRedSplit2Param::SC_UNIFORM_NZ){
if (fabs(xbar[currvar]) > param.getEPS()){
factor = param.getColumnScalingBoundLAP();
}
}
for (int j = 0; j < mTab; ++j){
workNonBasicTab[j][nTab] = intNonBasicTab[j][colpos]*factor;
}
nTab++;
}
else{
for (colpos = 0; colpos < card_contNonBasicVar; ++colpos){
if (contNonBasicVar[colpos] == currvar){
#ifdef RS2_TRACE
printf("Found cont nonbasic at pos %d: %d %d\n", colpos, contNonBasicVar[colpos], currvar);
#endif
break;
}
}
double factor = 1.0;
if (scaling == CglRedSplit2Param::SC_LINEAR){
if (fabs(xbar[currvar]) > factor){
factor = fabs(xbar[currvar]);
}
}
else if (scaling == CglRedSplit2Param::SC_LINEAR_BOUNDED){
if (fabs(xbar[currvar]) > factor){
factor = fabs(xbar[currvar]);
}
if (factor < param.getColumnScalingBoundLAP()){
factor = param.getColumnScalingBoundLAP();
}
}
else if (scaling == CglRedSplit2Param::SC_LOG_BOUNDED){
if (log(fabs(xbar[currvar])) > factor){
factor = log(fabs(xbar[currvar]));
}
if (factor < param.getColumnScalingBoundLAP()){
factor = param.getColumnScalingBoundLAP();
}
}
else if (scaling == CglRedSplit2Param::SC_UNIFORM){
factor = param.getColumnScalingBoundLAP();
}
else if (scaling == CglRedSplit2Param::SC_UNIFORM_NZ){
if (fabs(xbar[currvar]) > param.getEPS()){
factor = param.getColumnScalingBoundLAP();
}
}
for (int j = 0; j < mTab; ++j){
#ifdef RS2_TRACE
printf("Col %d Row %d: %f, xbar %f\n", colpos, j, contNonBasicTab[j][colpos], xbar[currvar]);
#endif
workNonBasicTab[j][nTab] = contNonBasicTab[j][colpos]*factor;
}
nTab++;
}
i++;
}
#ifdef RS2_TRACE
printf("Printing workNonBasicTab:\n");
#endif
#if RS_FAST_WORK
int workOffset = mTab+card_intNonBasicVar+card_contNonBasicVar+2;
#endif
for (i = 0; i < mTab; ++i) {
#ifdef RS2_TRACE
for (int j = 0; j < nTab; ++j){
printf("%.6f ", workNonBasicTab[i][j]);
}
printf("\n");
#endif
#if RS_FAST_WORK == 0
norm[i] = rs_dotProd(workNonBasicTab[i], workNonBasicTab[i], nTab);
#elif RS_FAST_WORK == 1
norm[i] = rs_dotProd(workNonBasicTab[i], workNonBasicTab[i], nTab);
double value = rs_sparseDotProd(workNonBasicTab[i],
workNonBasicTab[i],
pi_mat[i]+workOffset,
pi_mat[i]+workOffset);
assert (value==norm[i]);
#else
norm[i] = rs_sparseDotProd(workNonBasicTab[i],
workNonBasicTab[i],
pi_mat[i]+workOffset,
pi_mat[i]+workOffset);
#endif
}
} /* fill_workNonBasicTab */
/***********************************************************************/