// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include <stdio.h>
#include <math.h>
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolveDual.hpp"
#include "CoinMessage.hpp"
#include "CoinHelperFunctions.hpp"
//#define PRESOLVE_TIGHTEN_DUALS 1
//#define PRESOLVE_DEBUG 1
// this looks for "dominated columns"
// following ekkredc
// rdmin == dwrow2
// rdmax == dwrow
// this transformation is applied in: k4.mps, lp22.mps
// make inferences using this constraint on reduced cost:
// at optimality dj>0 ==> var must be at lb
// dj<0 ==> var must be at ub
//
// This implies:
// there is no lb ==> dj<=0 at optimality
// there is no ub ==> dj>=0 at optimality
/*
This routine looks to be something of a work in progres. See the comment
that begins with `Gack!'. And down in the bound propagation loop, why do we
only work with variables with u_j = infty? The corresponding section of code
for l_j = -infty is ifdef'd away. And why exclude the code protected by
PRESOLVE_TIGHTEN_DUALS?
*/
const CoinPresolveAction *remove_dual_action::presolve(CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
double startTime = 0.0;
int startEmptyRows=0;
int startEmptyColumns = 0;
if (prob->tuning_) {
startTime = CoinCpuTime();
startEmptyRows = prob->countEmptyRows();
startEmptyColumns = prob->countEmptyCols();
}
double *colels = prob->colels_;
int *hrow = prob->hrow_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int ncols = prob->ncols_;
double *clo = prob->clo_;
double *cup = prob->cup_;
double *rowels = prob->rowels_;
int *hcol = prob->hcol_;
CoinBigIndex *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
double *csol = prob->sol_;
int nrows = prob->nrows_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
double *dcost = prob->cost_;
const double maxmin = prob->maxmin_;
const double ekkinf = 1e28;
const double ekkinf2 = 1e20;
const double ztoldj = prob->ztoldj_;
double *rdmin = new double[nrows];
double *rdmax = new double[nrows];
// This combines the initialization of rdmin/rdmax to extreme values
// (PRESOLVE_INF/-PRESOLVE_INF) with a version of the next loop specialized
// for row slacks.
// In this case, it is always the case that dprice==0.0 and coeff==1.0.
int i;
for ( i = 0; i < nrows; i++) {
double sup = -rlo[i]; // slack ub; corresponds to cup[j]
double slo = -rup[i]; // slack lb; corresponds to clo[j]
bool no_lb = (slo <= -ekkinf);
bool no_ub = (sup >= ekkinf);
// here, dj==-row dual
// the row slack has no lower bound, but it does have an upper bound,
// then the reduced cost must be <= 0, so the row dual must be >= 0
rdmin[i] = (no_lb && !no_ub) ? 0.0 : -PRESOLVE_INF;
rdmax[i] = (no_ub && !no_lb) ? 0.0 : PRESOLVE_INF;
}
// Look for col singletons and update bounds on dual costs
// Take the min of the maxes and max of the mins
int j;
for ( j = 0; j<ncols; j++) {
bool no_ub = (cup[j] >= ekkinf);
bool no_lb = (clo[j] <= -ekkinf);
if (hincol[j] == 1 &&
// we need singleton cols that have exactly one bound
(no_ub ^ no_lb)) {
int row = hrow[mcstrt[j]];
double coeff = colels[mcstrt[j]];
PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
// I don't think the sign of dcost[j] matters
// this row dual would make this col's reduced cost be 0
double dprice = maxmin * dcost[j] / coeff;
// no_ub == !no_lb
// no_ub ==> !(dj<0)
// no_lb ==> !(dj>0)
// I don't think the case where !no_ub has actually been tested
if ((coeff > 0.0) == no_ub) {
// coeff>0 ==> making the row dual larger would make dj *negative*
// ==> dprice is an upper bound on dj if no *ub*
// coeff<0 ==> making the row dual larger would make dj *positive*
// ==> dprice is an upper bound on dj if no *lb*
if (rdmax[row] > dprice) // reduced cost may be positive
rdmax[row] = dprice;
} else { // no lower bound
if (rdmin[row] < dprice) // reduced cost may be negative
rdmin[row] = dprice;
}
}
}
int *fixup_cols = new int[ncols];
int *fixdown_cols = new int[ncols];
#if PRESOLVE_TIGHTEN_DUALS
double *djmin = new double[ncols];
double *djmax = new double[ncols];
#endif
int nfixup_cols = 0;
int nfixdown_cols = 0;
while (1) {
int tightened = 0;
/* Perform duality tests */
for (int j = 0; j<ncols; j++) {
if (hincol[j] > 0) {
CoinBigIndex kcs = mcstrt[j];
CoinBigIndex kce = kcs + hincol[j];
// Number of infinite rows
int nflagu = 0;
int nflagl = 0;
// Number of ordinary rows
int nordu = 0;
int nordl = 0;
double ddjlo = maxmin * dcost[j];
double ddjhi = ddjlo;
for (CoinBigIndex k = kcs; k < kce; k++) {
int i = hrow[k];
double coeff = colels[k];
if (coeff > 0.0) {
if (rdmin[i] >= -ekkinf2) {
ddjhi -= coeff * rdmin[i];
nordu ++;
} else {
nflagu ++;
}
if (rdmax[i] <= ekkinf2) {
ddjlo -= coeff * rdmax[i];
nordl ++;
} else {
nflagl ++;
}
} else {
if (rdmax[i] <= ekkinf2) {
ddjhi -= coeff * rdmax[i];
nordu ++;
} else {
nflagu ++;
}
if (rdmin[i] >= -ekkinf2) {
ddjlo -= coeff * rdmin[i];
nordl ++;
} else {
nflagl ++;
}
}
}
// See if we may be able to tighten a dual
if (cup[j]>ekkinf) {
// dj can not be negative
if (nflagu==1&&ddjhi<-ztoldj) {
// We can make bound finite one way
for (CoinBigIndex k = kcs; k < kce; k++) {
int i = hrow[k];
double coeff = colels[k];
if (coeff > 0.0&&rdmin[i] < -ekkinf2) {
// rdmax[i] has upper bound
if (ddjhi<rdmax[i]*coeff-ztoldj) {
double newValue = ddjhi/coeff;
// re-compute lo
if (rdmax[i] > ekkinf2 && newValue <= ekkinf2) {
nflagl--;
ddjlo -= coeff * newValue;
} else if (rdmax[i] <= ekkinf2) {
ddjlo -= coeff * (newValue-rdmax[i]);
}
rdmax[i] = newValue;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
#endif
}
} else if (coeff < 0.0 && rdmax[i] > ekkinf2) {
// rdmin[i] has lower bound
if (ddjhi<rdmin[i]*coeff-ztoldj) {
double newValue = ddjhi/coeff;
// re-compute lo
if (rdmin[i] < -ekkinf2 && newValue >= -ekkinf2) {
nflagl--;
ddjlo -= coeff * newValue;
} else if (rdmax[i] >= -ekkinf2) {
ddjlo -= coeff * (newValue-rdmin[i]);
}
rdmin[i] = newValue;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
#endif
ddjlo = 0.0;
}
}
}
} else if (nflagl==0&&nordl==1&&ddjlo<-ztoldj) {
// We may be able to tighten
for (CoinBigIndex k = kcs; k < kce; k++) {
int i = hrow[k];
double coeff = colels[k];
if (coeff > 0.0) {
rdmax[i] += ddjlo/coeff;
ddjlo =0.0;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
#endif
} else if (coeff < 0.0 ) {
rdmin[i] += ddjlo/coeff;
ddjlo =0.0;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
#endif
}
}
}
}
#if 0
if (clo[j]<-ekkinf) {
// dj can not be positive
if (ddjlo>ztoldj&&nflagl==1) {
// We can make bound finite one way
for (CoinBigIndex k = kcs; k < kce; k++) {
int i = hrow[k];
double coeff = colels[k];
if (coeff < 0.0&&rdmin[i] < -ekkinf2) {
// rdmax[i] has upper bound
if (ddjlo>rdmax[i]*coeff+ztoldj) {
double newValue = ddjlo/coeff;
// re-compute hi
if (rdmax[i] > ekkinf2 && newValue <= ekkinf2) {
nflagu--;
ddjhi -= coeff * newValue;
} else if (rdmax[i] <= ekkinf2) {
ddjhi -= coeff * (newValue-rdmax[i]);
}
rdmax[i] = newValue;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
#endif
}
} else if (coeff > 0.0 && rdmax[i] > ekkinf2) {
// rdmin[i] has lower bound
if (ddjlo>rdmin[i]*coeff+ztoldj) {
double newValue = ddjlo/coeff;
// re-compute lo
if (rdmin[i] < -ekkinf2 && newValue >= -ekkinf2) {
nflagu--;
ddjhi -= coeff * newValue;
} else if (rdmax[i] >= -ekkinf2) {
ddjhi -= coeff * (newValue-rdmin[i]);
}
rdmin[i] = newValue;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
#endif
}
}
}
} else if (nflagu==0&&nordu==1&&ddjhi>ztoldj) {
// We may be able to tighten
for (CoinBigIndex k = kcs; k < kce; k++) {
int i = hrow[k];
double coeff = colels[k];
if (coeff < 0.0) {
rdmax[i] += ddjhi/coeff;
ddjhi =0.0;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
#endif
} else if (coeff > 0.0 ) {
rdmin[i] += ddjhi/coeff;
ddjhi =0.0;
tightened++;
#if PRESOLVE_DEBUG
printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
#endif
}
}
}
}
#endif
#if PRESOLVE_TIGHTEN_DUALS
djmin[j] = (nflagl ? -PRESOLVE_INF : ddjlo);
djmax[j] = (nflagu ? PRESOLVE_INF : ddjhi);
#endif
if (ddjlo > ztoldj && nflagl == 0&&!prob->colProhibited2(j)) {
// dj>0 at optimality ==> must be at lower bound
if (clo[j] <= -ekkinf) {
prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDB,
prob->messages())
<<j
<<CoinMessageEol;
prob->status_ |= 2;
break;
} else {
fixdown_cols[nfixdown_cols++] = j;
//if (csol[j]-clo[j]>1.0e-7)
//printf("down %d row %d nincol %d\n",j,hrow[mcstrt[j]],hincol[j]);
// User may have given us feasible solution - move if simple
if (csol&&csol[j]-clo[j]>1.0e-7&&hincol[j]==1) {
double value_j = colels[mcstrt[j]];
double distance_j = csol[j]-clo[j];
int row=hrow[mcstrt[j]];
// See if another column can take value
for (CoinBigIndex kk=mrstrt[row];kk<mrstrt[row]+hinrow[row];kk++) {
int k = hcol[kk];
if (hincol[k]==1&&k!=j) {
double value_k = rowels[kk];
double movement;
if (value_k*value_j>0.0) {
// k needs to increase
double distance_k = cup[k]-csol[k];
movement = CoinMin((distance_j*value_j)/value_k,distance_k);
} else {
// k needs to decrease
double distance_k = clo[k]-csol[k];
movement = CoinMax((distance_j*value_j)/value_k,distance_k);
}
csol[k] += movement;
distance_j -= (movement*value_k)/value_j;
csol[j] -= (movement*value_k)/value_j;
if (distance_j<1.0e-7)
break;
}
}
}
}
} else if (ddjhi < -ztoldj && nflagu == 0&&!prob->colProhibited2(j)) {
// dj<0 at optimality ==> must be at upper bound
if (cup[j] >= ekkinf) {
prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDA,
prob->messages())
<<j
<<CoinMessageEol;
prob->status_ |= 2;
break;
} else {
fixup_cols[nfixup_cols++] = j;
// User may have given us feasible solution - move if simple
//if (cup[j]-csol[j]>1.0e-7)
//printf("up %d row %d nincol %d\n",j,hrow[mcstrt[j]],hincol[j]);
if (csol&&cup[j]-csol[j]>1.0e-7&&hincol[j]==1) {
double value_j = colels[mcstrt[j]];
double distance_j = csol[j]-cup[j];
int row=hrow[mcstrt[j]];
// See if another column can take value
for (CoinBigIndex kk=mrstrt[row];kk<mrstrt[row]+hinrow[row];kk++) {
int k = hcol[kk];
if (hincol[k]==1&&k!=j) {
double value_k = rowels[kk];
double movement;
if (value_k*value_j<0.0) {
// k needs to increase
double distance_k = cup[k]-csol[k];
movement = CoinMin((distance_j*value_j)/value_k,distance_k);
} else {
// k needs to decrease
double distance_k = clo[k]-csol[k];
movement = CoinMax((distance_j*value_j)/value_k,distance_k);
}
csol[k] += movement;
distance_j -= (movement*value_k)/value_j;
csol[j] -= (movement*value_k)/value_j;
if (distance_j>-1.0e-7)
break;
}
}
}
}
}
}
}
// I don't know why I stopped doing this.
#if PRESOLVE_TIGHTEN_DUALS
const double *rowels = prob->rowels_;
const int *hcol = prob->hcol_;
const CoinBigIndex *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
// tighten row dual bounds, as described on p. 229
for (int i = 0; i < nrows; i++) {
bool no_ub = (rup[i] >= ekkinf);
bool no_lb = (rlo[i] <= -ekkinf);
if ((no_ub ^ no_lb) == true) {
CoinBigIndex krs = mrstrt[i];
CoinBigIndex kre = krs + hinrow[i];
double rmax = rdmax[i];
double rmin = rdmin[i];
// all row columns are non-empty
for (CoinBigIndex k=krs; k<kre; k++) {
double coeff = rowels[k];
int icol = hcol[k];
double djmax0 = djmax[icol];
double djmin0 = djmin[icol];
if (no_ub) {
// dj must not be negative
if (coeff > ZTOLDP && djmax0 <PRESOLVE_INF && cup[icol]>=ekkinf) {
double bnd = djmax0 / coeff;
if (rmax > bnd) {
#if PRESOLVE_DEBUG
printf("MAX TIGHT[%d,%d]: %g --> %g\n", i,hrow[k], rdmax[i], bnd);
#endif
rdmax[i] = rmax = bnd;
tightened ++;;
}
} else if (coeff < -ZTOLDP && djmax0 <PRESOLVE_INF && cup[icol] >= ekkinf) {
double bnd = djmax0 / coeff ;
if (rmin < bnd) {
#if PRESOLVE_DEBUG
printf("MIN TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], rdmin[i], bnd);
#endif
rdmin[i] = rmin = bnd;
tightened ++;;
}
}
} else { // no_lb
// dj must not be positive
if (coeff > ZTOLDP && djmin0 > -PRESOLVE_INF && clo[icol]<=-ekkinf) {
double bnd = djmin0 / coeff ;
if (rmin < bnd) {
#if PRESOLVE_DEBUG
printf("MIN1 TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], rdmin[i], bnd);
#endif
rdmin[i] = rmin = bnd;
tightened ++;;
}
} else if (coeff < -ZTOLDP && djmin0 > -PRESOLVE_INF && clo[icol] <= -ekkinf) {
double bnd = djmin0 / coeff ;
if (rmax > bnd) {
#if PRESOLVE_DEBUG
printf("MAX TIGHT1[%d,%d]: %g --> %g\n", i,hrow[k], rdmax[i], bnd);
#endif
rdmax[i] = rmax = bnd;
tightened ++;;
}
}
}
}
}
}
#endif
if (tightened<100||nfixdown_cols||nfixup_cols)
break;
#if PRESOLVE_TIGHTEN_DUALS
else
printf("DUAL TIGHTENED! %d\n", tightened);
#endif
}
if (nfixup_cols) {
#if PRESOLVE_DEBUG
printf("NDUAL: %d\n", nfixup_cols);
#endif
next = make_fixed_action::presolve(prob, fixup_cols, nfixup_cols, false, next);
}
if (nfixdown_cols) {
#if PRESOLVE_DEBUG
printf("NDUAL: %d\n", nfixdown_cols);
#endif
next = make_fixed_action::presolve(prob, fixdown_cols, nfixdown_cols, true, next);
}
// If dual says so then we can make equality row
// Also if cost is in right direction and only one binding row for variable
// We may wish to think about giving preference to rows with 2 or 3 elements
/*
Gack! Ok, I can appreciate the thought here, but I'm seriously sceptical
about writing canFix[0] before reading rdmin[0]. After that, we should be out
of the interference zone for the typical situation where sizeof(double) is
twice sizeof(int).
*/
int * canFix = (int *) rdmin;
for ( i = 0; i < nrows; i++) {
bool no_lb = (rlo[i] <= -ekkinf);
bool no_ub = (rup[i] >= ekkinf);
canFix[i]=0;
if (no_ub && !no_lb ) {
if ( rdmin[i]>0.0)
canFix[i]=-1;
else
canFix[i]=-2;
} else if (no_lb && !no_ub ) {
if (rdmax[i]<0.0)
canFix[i]=1;
else
canFix[i]=2;
}
}
for (j = 0; j<ncols; j++) {
if (hincol[j]<=1)
continue;
CoinBigIndex kcs = mcstrt[j];
CoinBigIndex kce = kcs + hincol[j];
int bindingUp=-1;
int bindingDown=-1;
if (cup[j]<ekkinf)
bindingUp=-2;
if (clo[j]>-ekkinf)
bindingDown=-2;
for (CoinBigIndex k = kcs; k < kce; k++) {
int i = hrow[k];
if (abs(canFix[i])!=2) {
bindingUp=-2;
bindingDown=-2;
break;
}
double coeff = colels[k];
if (coeff>0.0) {
if (canFix[i]==2) {
// binding up
if (bindingUp==-1)
bindingUp = i;
else
bindingUp = -2;
} else {
// binding down
if (bindingDown==-1)
bindingDown = i;
else
bindingDown = -2;
}
} else {
if (canFix[i]==2) {
// binding down
if (bindingDown==-1)
bindingDown = i;
else
bindingDown = -2;
} else {
// binding up
if (bindingUp==-1)
bindingUp = i;
else
bindingUp = -2;
}
}
}
double cost = maxmin * dcost[j];
if (bindingUp>-2&&cost<=0.0) {
// might as well make equality
if (bindingUp>=0) {
canFix[bindingUp] /= 2; //So -2 goes to -1 etc
//printf("fixing row %d to ub by %d\n",bindingUp,j);
} else {
//printf("binding up row by %d\n",j);
}
} else if (bindingDown>-2 &&cost>=0.0) {
// might as well make equality
if (bindingDown>=0) {
canFix[bindingDown] /= 2; //So -2 goes to -1 etc
//printf("fixing row %d to lb by %d\n",bindingDown,j);
} else {
//printf("binding down row by %d\n",j);
}
}
}
// can't fix if integer and non-unit coefficient
//const double *rowels = prob->rowels_;
//const int *hcol = prob->hcol_;
//const CoinBigIndex *mrstrt = prob->mrstrt_;
//int *hinrow = prob->hinrow_;
const unsigned char *integerType = prob->integerType_;
for ( i = 0; i < nrows; i++) {
if (abs(canFix[i])==1) {
CoinBigIndex krs = mrstrt[i];
CoinBigIndex kre = krs + hinrow[i];
for (CoinBigIndex k=krs; k<kre; k++) {
int icol = hcol[k];
if (cup[icol]>clo[icol]&&integerType[icol]) {
canFix[i]=0; // not safe
}
}
}
if (canFix[i]==1) {
rlo[i]=rup[i];
prob->addRow(i);
} else if (canFix[i]==-1) {
rup[i]=rlo[i];
prob->addRow(i);
}
}
delete[]rdmin;
delete[]rdmax;
delete[]fixup_cols;
delete[]fixdown_cols;
#if PRESOLVE_TIGHTEN_DUALS
delete[]djmin;
delete[]djmax;
#endif
if (prob->tuning_) {
double thisTime=CoinCpuTime();
int droppedRows = prob->countEmptyRows() - startEmptyRows;
int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
printf("CoinPresolveDual(1) - %d rows, %d columns dropped in time %g, total %g\n",
droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
}
return (next);
}
syntax highlighted by Code2HTML, v. 0.9.1