limp-cbc-0.3.2.0: cbits/coin/CoinPresolveDupcol.cpp
/* $Id: CoinPresolveDupcol.cpp 1585 2013-04-06 20:42:02Z stefan $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
#include <stdio.h>
#include <math.h>
//#define PRESOLVE_DEBUG 1
// Debugging macros/functions
//#define PRESOLVE_DETAIL 1
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolveDupcol.hpp"
#include "CoinSort.hpp"
#include "CoinFinite.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPresolveUseless.hpp"
#include "CoinMessage.hpp"
#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
#include "CoinPresolvePsdebug.hpp"
#endif
#define DSEED2 2147483647.0
// Can be used from anywhere
void coin_init_random_vec(double *work, int n)
{
double deseed = 12345678.0;
for (int i = 0; i < n; ++i) {
deseed *= 16807.;
int jseed = static_cast<int> (deseed / DSEED2);
deseed -= static_cast<double> (jseed) * DSEED2;
double random = deseed / DSEED2;
work[i]=random;
}
}
namespace { // begin unnamed file-local namespace
/*
For each candidate major-dimension vector in majcands, calculate the sum
over the vector, with each minor dimension weighted by a random amount.
(E.g., calculate column sums with each row weighted by a random amount.)
The sums are returned in the corresponding entries of majsums.
*/
void compute_sums (int /*n*/, const int *majlens, const CoinBigIndex *majstrts,
int *minndxs, double *elems, const double *minmuls,
int *majcands, double *majsums, int nlook)
{ for (int cndx = 0 ; cndx < nlook ; ++cndx)
{ int i = majcands[cndx] ;
PRESOLVEASSERT(majlens[i] > 0) ;
CoinBigIndex kcs = majstrts[i] ;
CoinBigIndex kce = kcs + majlens[i] ;
double value = 0.0 ;
for (CoinBigIndex k = kcs ; k < kce ; k++)
{ int irow = minndxs[k] ;
value += minmuls[irow]*elems[k] ; }
majsums[cndx] = value ; }
return ; }
void create_col (int col, int n, double *els,
CoinBigIndex *mcstrt, double *colels, int *hrow, int *link,
CoinBigIndex *free_listp)
{
int *rows = reinterpret_cast<int *>(els+n) ;
CoinBigIndex free_list = *free_listp;
int xstart = NO_LINK;
for (int i=0; i<n; ++i) {
CoinBigIndex k = free_list;
assert(k >= 0) ;
free_list = link[free_list];
hrow[k] = rows[i];
colels[k] = els[i];
link[k] = xstart;
xstart = k;
}
mcstrt[col] = xstart;
*free_listp = free_list;
}
} // end unnamed file-local namespace
const char *dupcol_action::name () const
{
return ("dupcol_action");
}
/*
Original comment: This is just ekkredc5, adapted into the new framework.
The datasets scorpion.mps and allgrade.mps have duplicate columns.
In case you don't have your OSL manual handy, a somewhat more informative
explanation: We're looking for an easy-to-detect special case of linearly
dependent columns, where the coefficients of the duplicate columns are
exactly equal. The idea for locating such columns is to generate pseudo-
random weights for each row and then calculate the weighted sum of
coefficients of each column. Columns with equal sums are checked more
thoroughly.
Analysis of the situation says there are two major cases:
* If the columns have equal objective coefficients, we can combine
them.
* If the columns have unequal objective coefficients, we may be able to
fix one at bound. If the required bound doesn't exist, we have dual
infeasibility (hence one of primal infeasibility or unboundedness).
In the comments below are a few fragments of code from the original
routine. I don't think they make sense, but I've left them for the nonce in
case someone else recognises the purpose. -- lh, 040909 --
*/
const CoinPresolveAction
*dupcol_action::presolve (CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
# if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
# if PRESOLVE_DEBUG > 0
std::cout
<< "Entering dupcol_action::presolve." << std::endl ;
# endif
presolve_consistent(prob) ;
presolve_links_ok(prob) ;
presolve_check_sol(prob) ;
presolve_check_nbasic(prob) ;
# endif
# if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
int startEmptyRows = 0 ;
int startEmptyColumns = 0 ;
startEmptyRows = prob->countEmptyRows() ;
startEmptyColumns = prob->countEmptyCols() ;
# if COIN_PRESOLVE_TUNING > 0
double startTime = 0.0;
if (prob->tuning_) startTime = CoinCpuTime() ;
# endif
# endif
double maxmin = prob->maxmin_ ;
double *colels = prob->colels_ ;
int *hrow = prob->hrow_ ;
CoinBigIndex *mcstrt = prob->mcstrt_ ;
int *hincol = prob->hincol_ ;
int ncols = prob->ncols_ ;
int nrows = prob->nrows_ ;
double *clo = prob->clo_ ;
double *cup = prob->cup_ ;
double *sol = prob->sol_ ;
double *rlo = prob->rlo_ ;
double *rup = prob->rup_ ;
// If all coefficients positive do more simply
bool allPositive=true;
double * rhs = prob->usefulRowDouble_; //new double[nrows];
CoinMemcpyN(rup,nrows,rhs);
/*
Scan the columns for candidates, and write the indices into sort. We're not
interested in columns that are empty, prohibited, or integral.
Question: Should we exclude singletons, which are useful in other transforms?
Question: Why are we excluding integral columns?
*/
// allow integral columns if asked for
bool allowIntegers = ((prob->presolveOptions_&0x01) != 0) ;
int *sort = prob->usefulColumnInt_; //new int[ncols] ;
int nlook = 0 ;
for (int j = 0 ; j < ncols ; j++) {
if (hincol[j] == 0) continue ;
// sort
CoinSort_2(hrow+mcstrt[j],hrow+mcstrt[j]+hincol[j],
colels+mcstrt[j]);
// check all positive and adjust rhs
if (allPositive) {
double lower = clo[j];
if (lower<cup[j]) {
for (int k=mcstrt[j];k<mcstrt[j]+hincol[j];k++) {
double value=colels[k];
if (value<0.0)
allPositive=false;
else
rhs[hrow[k]] -= lower*value;
}
} else {
for (int k=mcstrt[j];k<mcstrt[j]+hincol[j];k++) {
double value=colels[k];
rhs[hrow[k]] -= lower*value;
}
}
}
if (prob->colProhibited2(j)) continue ;
//#define PRESOLVE_INTEGER_DUPCOL
#ifndef PRESOLVE_INTEGER_DUPCOL
if (prob->isInteger(j)&&!allowIntegers) continue ;
#endif
sort[nlook++] = j ; }
if (nlook == 0)
{ //delete[] sort ;
//delete [] rhs;
return (next) ; }
/*
Prep: add the coefficients of each candidate column. To reduce false
positives, multiply each row by a `random' multiplier when forming the
sums. On return from compute_sums, sort and colsum are loaded with the
indices and column sums, respectively, of candidate columns. The pair of
arrays are then sorted by sum so that equal sums are adjacent.
*/
double *colsum = prob->usefulColumnDouble_; //new double[ncols] ;
double *rowmul;
if (!prob->randomNumber_) {
rowmul = new double[nrows] ;
coin_init_random_vec(rowmul,nrows) ;
} else {
rowmul = prob->randomNumber_;
}
compute_sums(ncols,hincol,mcstrt,hrow,colels,rowmul,sort,colsum,nlook) ;
CoinSort_2(colsum,colsum+nlook,sort) ;
/*
General prep --- unpack the various vectors we'll need, and allocate arrays
to record the results.
*/
presolvehlink *clink = prob->clink_ ;
double *rowels = prob->rowels_ ;
int *hcol = prob->hcol_ ;
const CoinBigIndex *mrstrt = prob->mrstrt_ ;
int *hinrow = prob->hinrow_ ;
double *dcost = prob->cost_ ;
action *actions = new action [nlook] ;
int nactions = 0 ;
# ifdef ZEROFAULT
memset(actions,0,nlook*sizeof(action)) ;
# endif
int *fixed_down = new int[nlook] ;
int nfixed_down = 0 ;
int *fixed_up = new int[nlook] ;
int nfixed_up = 0 ;
#if 0
// Excluded in the original routine. I'm guessing it's excluded because
// it's just not cost effective to worry about this. -- lh, 040908 --
// It may be the case that several columns are duplicate.
// If not all have the same cost, then we have to make sure
// that we set the most expensive one to its minimum
// now sort in each class by cost
{
double dval = colsum[0] ;
int first = 0 ;
for (int jj = 1; jj < nlook; jj++) {
while (colsum[jj]==dval)
jj++ ;
if (first + 1 < jj) {
double buf[jj - first] ;
for (int i=first; i<jj; ++i)
buf[i-first] = dcost[sort[i]]*maxmin ;
CoinSort_2(buf,buf+jj-first,sort+first) ;
//ekk_sortonDouble(buf,&sort[first],jj-first) ;
}
}
}
#endif
// We will get all min/max but only if needed
bool gotStuff=false;
/*
Original comment: It appears to be the case that this loop is finished,
there may still be duplicate cols left. I haven't done anything
about that yet.
Open the main loop to compare column pairs. We'll compare sort[jj] to
sort[tgt]. This allows us to accumulate multiple columns into one. But
we don't manage all-pairs comparison when we can't combine columns.
We can quickly dismiss pairs which have unequal sums or lengths.
*/
int isorted = -1 ;
int tgt = 0 ;
for (int jj = 1 ; jj < nlook ; jj++)
{ if (colsum[jj] != colsum[jj-1]) {
tgt = jj; // Must update before continuing
continue ;
}
int j2 = sort[jj] ;
int j1 = sort[tgt] ;
int len2 = hincol[j2] ;
int len1 = hincol[j1] ;
if (len2 != len1)
{ tgt = jj ;
continue ; }
/*
The final test: sort the columns by row index and compare each index and
coefficient.
*/
CoinBigIndex kcs = mcstrt[j2] ;
CoinBigIndex kce = kcs+hincol[j2] ;
int ishift = mcstrt[j1]-kcs ;
if (len1 > 1 && isorted < j1)
{ CoinSort_2(hrow+mcstrt[j1],hrow+mcstrt[j1]+len1,
colels+mcstrt[j1]) ;
isorted = j1 ; }
if (len2 > 1 && isorted < j2)
{ CoinSort_2(hrow+kcs,hrow+kcs+len2,colels+kcs) ;
isorted = j2 ; }
CoinBigIndex k ;
for (k = kcs ; k < kce ; k++)
{ if (hrow[k] != hrow[k+ishift] || colels[k] != colels[k+ishift])
{ break ; } }
if (k != kce)
{ tgt = jj ;
continue ; }
/*
These really are duplicate columns. Grab values for convenient reference.
Convert the objective coefficients for minimization.
*/
double clo1 = clo[j1] ;
double cup1 = cup[j1] ;
double clo2 = clo[j2] ;
double cup2 = cup[j2] ;
double c1 = dcost[j1]*maxmin ;
double c2 = dcost[j2]*maxmin ;
PRESOLVEASSERT(!(clo1 == cup1 || clo2 == cup2)) ;
// Get reasonable bounds on sum of two variables
double lowerBound=-COIN_DBL_MAX;
double upperBound=COIN_DBL_MAX;
// For now only if lower bounds are zero
if (!clo1&&!clo2) {
// Only need bounds if c1 != c2
if (c1!=c2) {
if (!allPositive) {
#if 0
for (k=kcs;k<kce;k++) {
int iRow = hrow[k];
bool posinf = false;
bool neginf = false;
double maxup = 0.0;
double maxdown = 0.0;
// compute sum of all bounds except for j1,j2
CoinBigIndex kk;
CoinBigIndex kre = mrstrt[iRow]+hinrow[iRow];
double value1=0.0;
for (kk=mrstrt[iRow]; kk<kre; kk++) {
int col = hcol[kk];
if (col == j1||col==j2) {
value1=rowels[kk];
continue;
}
double coeff = rowels[kk];
double lb = clo[col];
double ub = cup[col];
if (coeff > 0.0) {
if (PRESOLVE_INF <= ub) {
posinf = true;
if (neginf)
break; // pointless
} else {
maxup += ub * coeff;
}
if (lb <= -PRESOLVE_INF) {
neginf = true;
if (posinf)
break; // pointless
} else {
maxdown += lb * coeff;
}
} else {
if (PRESOLVE_INF <= ub) {
neginf = true;
if (posinf)
break; // pointless
} else {
maxdown += ub * coeff;
}
if (lb <= -PRESOLVE_INF) {
posinf = true;
if (neginf)
break; // pointless
} else {
maxup += lb * coeff;
}
}
}
if (kk==kre) {
assert (value1);
if (value1>1.0e-5) {
if (!neginf&&rup[iRow]<1.0e10)
if (upperBound*value1>rup[iRow]-maxdown)
upperBound = (rup[iRow]-maxdown)/value1;
if (!posinf&&rlo[iRow]>-1.0e10)
if (lowerBound*value1<rlo[iRow]-maxup)
lowerBound = (rlo[iRow]-maxup)/value1;
} else if (value1<-1.0e-5) {
if (!neginf&&rup[iRow]<1.0e10)
if (lowerBound*value1>rup[iRow]-maxdown) {
#ifndef NDEBUG
double x=lowerBound;
#endif
lowerBound = (rup[iRow]-maxdown)/value1;
assert (lowerBound == CoinMax(x,(rup[iRow]-maxdown)/value1));
}
if (!posinf&&rlo[iRow]>-1.0e10)
if (upperBound*value1<rlo[iRow]-maxup) {
#ifndef NDEBUG
double x=upperBound;
#endif
upperBound = (rlo[iRow]-maxup)/value1;
assert(upperBound == CoinMin(x,(rlo[iRow]-maxup)/value1));
}
}
}
}
double l=lowerBound;
double u=upperBound;
#endif
if (!gotStuff) {
prob->recomputeSums(-1); // get min max
gotStuff=true;
}
int positiveInf=0;
int negativeInf=0;
double lo=0;
double up=0.0;
if (clo1<-PRESOLVE_INF)
negativeInf++;
else
lo+=clo1;
if (clo2<-PRESOLVE_INF)
negativeInf++;
else
lo+=clo2;
if (cup1>PRESOLVE_INF)
positiveInf++;
else
up+=cup1;
if (cup2>PRESOLVE_INF)
positiveInf++;
else
up+=cup2;
for (k=kcs;k<kce;k++) {
int iRow = hrow[k];
double value = colels[k];
int pInf = (value>0.0) ? positiveInf : negativeInf;
int nInf = (value>0.0) ? negativeInf : positiveInf;
int posinf = prob->infiniteUp_[iRow]-pInf;
int neginf = prob->infiniteDown_[iRow]-nInf;
if (posinf>0&&neginf>0)
continue; // this row can't bound
double maxup = prob->sumUp_[iRow];
double maxdown = prob->sumDown_[iRow];
if (value>0.0) {
maxdown -= value*lo;
maxup -= value*up;
} else {
maxdown -= value*up;
maxup -= value*lo;
}
if (value>1.0e-5) {
if (!neginf&&rup[iRow]<1.0e10)
if (upperBound*value>rup[iRow]-maxdown)
upperBound = (rup[iRow]-maxdown)/value;
if (!posinf&&rlo[iRow]>-1.0e10)
if (lowerBound*value<rlo[iRow]-maxup)
lowerBound = (rlo[iRow]-maxup)/value;
} else if (value<-1.0e-5) {
if (!neginf&&rup[iRow]<1.0e10)
if (lowerBound*value>rup[iRow]-maxdown) {
lowerBound = (rup[iRow]-maxdown)/value;
}
if (!posinf&&rlo[iRow]>-1.0e10)
if (upperBound*value<rlo[iRow]-maxup) {
upperBound = (rlo[iRow]-maxup)/value;
}
}
}
//assert (fabs(l-lowerBound)<1.0e-5&&fabs(u-upperBound)<1.0e-5);
} else {
// can do faster
lowerBound=0.0;
for (k=kcs;k<kce;k++) {
int iRow = hrow[k];
double value=colels[k];
if (upperBound*value>rhs[iRow])
upperBound = rhs[iRow]/value;
}
}
}
// relax a bit
upperBound -= 1.0e-9;
} else {
// Not sure what to do so give up
continue;
}
/*
There are two main cases: The objective coefficients are equal or unequal.
For equal objective coefficients c1 == c2, we can combine the columns by
making the substitution x<j1> = x'<j1> - x<j2>. This will eliminate column
sort[jj] = j2 and leave the new variable x' in column sort[tgt] = j1. tgt
doesn't move.
*/
if (c1 == c2)
{
#ifdef PRESOLVE_INTEGER_DUPCOL
if (!allowIntegers) {
if (prob->isInteger(j1)) {
if (!prob->isInteger(j2)) {
if (cup2 < upperBound) //if (!prob->colInfinite(j2))
continue;
else
cup2 = COIN_DBL_MAX;
}
} else if (prob->isInteger(j2)) {
if (cup1 < upperBound) //if (!prob->colInfinite(j1))
continue;
else
cup1 = COIN_DBL_MAX;
}
//printf("TakingINTeq\n");
}
#endif
/*
As far as the presolved lp, there's no difference between columns. But we
need this relation to hold in order to guarantee that we can split the
value of the combined column during postsolve without damaging the basis.
(The relevant case is when the combined column is basic --- we need to be
able to retain column j1 in the basis and make column j2 nonbasic.)
*/
if (!(clo2+cup1 <= clo1+cup2))
{ CoinSwap(j1,j2) ;
CoinSwap(clo1,clo2) ;
CoinSwap(cup1,cup2) ;
tgt = jj ; }
/*
Create the postsolve action before we start to modify the columns.
*/
# if PRESOLVE_DEBUG > 1
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d += %d\n",j1,j2,j1,j2)) ;
PRESOLVE_DETAIL_PRINT(printf("pre_dupcol %dC %dC E\n",j2,j1));
# endif
action *s = &actions[nactions++] ;
s->thislo = clo[j2] ;
s->thisup = cup[j2] ;
s->lastlo = clo[j1] ;
s->lastup = cup[j1] ;
s->ithis = j2 ;
s->ilast = j1 ;
s->nincol = hincol[j2] ;
s->colels = presolve_dupmajor(colels,hrow,hincol[j2],mcstrt[j2]) ;
/*
Combine the columns into column j1. Upper and lower bounds and solution
simply add, and the coefficients are unchanged.
I'm skeptical of pushing a bound to infinity like this, but leave it for now.
-- lh, 040908 --
*/
clo1 += clo2 ;
if (clo1 < -1.0e20)
{ clo1 = -PRESOLVE_INF ; }
clo[j1] = clo1 ;
cup1 += cup2 ;
if (cup1 > 1.0e20)
{ cup1 = PRESOLVE_INF ; }
cup[j1] = cup1 ;
if (sol)
{ sol[j1] += sol[j2] ; }
if (prob->colstat_)
{ if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
{ prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic); } }
/*
Empty column j2.
*/
dcost[j2] = 0.0 ;
if (sol)
{ sol[j2] = clo2 ; }
CoinBigIndex k2cs = mcstrt[j2] ;
CoinBigIndex k2ce = k2cs + hincol[j2] ;
for (CoinBigIndex k = k2cs ; k < k2ce ; ++k)
{ presolve_delete_from_row(hrow[k],j2,mrstrt,hinrow,hcol,rowels) ; }
hincol[j2] = 0 ;
PRESOLVE_REMOVE_LINK(clink,j2) ;
continue ; }
/*
Unequal reduced costs. In this case, we may be able to fix one of the columns
or prove dual infeasibility. Given column a_k, duals y, objective
coefficient c_k, the reduced cost cbar_k = c_k - dot(y,a_k). Given
a_1 = a_2, substitute for dot(y,a_1) in the relation for cbar_2 to get
cbar_2 = (c_2 - c_1) + cbar_1
Independent elements here are variable bounds l_k, u_k, and difference
(c_2 - c_1). Infinite bounds for l_k, u_k will constrain the sign of cbar_k.
Assume minimization. If you do the case analysis, you find these cases of
interest:
l_1 u_1 l_2 u_2 cbar_1 c_2-c_1 cbar_2 result
A any finite -inf any <= 0 > 0 <= 0 x_1 -> NBUB
B -inf any any finite <= 0 < 0 < 0 x_2 -> NBUB
C finite any any +inf >= 0 < 0 >= 0 x_1 -> NBLB
D any +inf finite any >= 0 > 0 >= 0 x_2 -> NBLB
E -inf any any +inf <= 0 < 0 >= 0 dual infeas
F any inf -inf any >= 0 > 0 <= 0 dual infeas
G any finite finite any > 0 no inference
H finite any any finite < 0 no inference
The cases labelled dual infeasible are primal unbounded.
To keep the code compact, we'll always aim to take x_2 to bound. In the cases
where x_1 should go to bound, we'll swap. The implementation is boolean
algebra. Define bits for infinite bounds and (c_2 > c_1), then look for the
correct patterns.
*/
else
{ int minterm = 0 ;
#ifdef PRESOLVE_INTEGER_DUPCOL
if (!allowIntegers) {
if (c2 > c1) {
if (cup1 < upperBound/*!prob->colInfinite(j1)*/ && (prob->isInteger(j1)||prob->isInteger(j2)))
continue ;
} else {
if (cup2 < upperBound/*!prob->colInfinite(j2)*/ && (prob->isInteger(j1)||prob->isInteger(j2)))
continue ;
}
//printf("TakingINTne\n");
}
#endif
bool swapped = false ;
#if PRESOLVE_DEBUG > 1
printf("bounds %g %g\n",lowerBound,upperBound);
#endif
if (c2 > c1) minterm |= 1<<0 ;
if (cup2 >= PRESOLVE_INF/*prob->colInfinite(j2)*/) minterm |= 1<<1 ;
if (clo2 <= -PRESOLVE_INF) minterm |= 1<<2 ;
if (cup1 >= PRESOLVE_INF/*prob->colInfinite(j1)*/) minterm |= 1<<3 ;
if (clo1 <= -PRESOLVE_INF) minterm |= 1<<4 ;
// for now be careful - just one special case
if (!clo1&&!clo2) {
if (c2 > c1 && cup1 >= upperBound)
minterm |= 1<<3;
else if (c2 < c1 && cup2 >= upperBound)
minterm |= 1<<1;
}
/*
The most common case in a well-formed system should be no inference. We're
looking for x00x1 (case G) and 0xx00 (case H). This is where we have the
potential to miss inferences: If there are three or more columns with the
same sum, sort[tgt] == j1 will only be compared to the second in the
group.
*/
if ((minterm&0x0d) == 0x1 || (minterm&0x13) == 0)
{ tgt = jj ;
continue ; }
/*
Next remove the unbounded cases, 1xx10 and x11x1.
*/
if ((minterm&0x13) == 0x12 || (minterm&0x0d) == 0x0d)
{ prob->setStatus(2) ;
# if PRESOLVE_DEBUG > 1
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) Unbounded\n",j1,j2)) ;
# endif
break ; }
/*
With the no inference and unbounded cases removed, all that's left are the
cases where we can push a variable to bound. Swap if necessary (x01x1 or
0xx10) so that we're always fixing index j2. This means that column
sort[tgt] = j1 will be fixed. Unswapped, we fix column sort[jj] = j2.
*/
if ((minterm&0x0d) == 0x05 || (minterm&0x13) == 0x02)
{ CoinSwap(j1, j2) ;
CoinSwap(clo1, clo2) ;
CoinSwap(cup1, cup2) ;
CoinSwap(c1, c2) ;
int tmp1 = minterm&0x18 ;
int tmp2 = minterm&0x06 ;
int tmp3 = minterm&0x01 ;
minterm = (tmp1>>2)|(tmp2<<2)|(tmp3^0x01) ;
swapped = true ; }
/*
Force x_2 to upper bound? (Case B, boolean 1X100, where X == don't care.)
*/
if ((minterm&0x13) == 0x10)
{ fixed_up[nfixed_up++] = j2 ;
# if PRESOLVE_DEBUG > 1
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBUB\n",j1,j2,j2)) ;
# endif
if (prob->colstat_)
{ if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
{ prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic) ; }
prob->setColumnStatus(j2,CoinPrePostsolveMatrix::atUpperBound) ; }
if (sol)
{ double delta2 = cup2-sol[j2] ;
sol[j2] = cup2 ;
sol[j1] -= delta2 ; }
if (swapped)
{ tgt = jj ; }
continue ; }
/*
Force x_2 to lower bound? (Case C, boolean X1011.)
*/
if ((minterm&0x0d) == 0x09)
{ fixed_down[nfixed_down++] = j2 ;
# if PRESOLVE_DEBUG > 1
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBLB\n",j1,j2,j2)) ;
# endif
if (prob->colstat_)
{ if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
{ prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic) ; }
prob->setColumnStatus(j2,CoinPrePostsolveMatrix::atLowerBound) ; }
if (sol)
{ double delta2 = clo2-sol[j2] ;
sol[j2] = clo2 ;
sol[j1] -= delta2 ; }
if (swapped)
{ tgt = jj ; }
continue ; } }
/*
We should never reach this point in the loop --- all cases force a new
iteration or loop termination. If we get here, something happened that we
didn't anticipate.
*/
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) UNEXPECTED!\n",j1,j2)) ; }
/*
What's left? Deallocate vectors, and call make_fixed_action to handle any
variables that were fixed to bound.
*/
if (rowmul != prob->randomNumber_)
delete[] rowmul ;
//delete[] colsum ;
//delete[] sort ;
//delete [] rhs;
# if PRESOLVE_SUMMARY || PRESOLVE_DEBUG
if (nactions+nfixed_down+nfixed_up > 0)
{ printf("DUPLICATE COLS: %d combined, %d lb, %d ub\n",
nactions,nfixed_down,nfixed_up) ; }
# endif
if (nactions)
{ next = new dupcol_action(nactions,CoinCopyOfArray(actions,nactions),next) ;
// we can't go round again in integer
prob->presolveOptions_ |= 0x80000000;
}
deleteAction(actions,action*) ;
if (nfixed_down)
{ next =
make_fixed_action::presolve(prob,fixed_down,nfixed_down,true,next) ; }
delete[]fixed_down ;
if (nfixed_up)
{ next =
make_fixed_action::presolve(prob,fixed_up,nfixed_up,false,next) ; }
delete[]fixed_up ;
# if COIN_PRESOLVE_TUNING > 0
if (prob->tuning_) double thisTime = CoinCpuTime() ;
# endif
# if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
presolve_check_sol(prob) ;
# endif
# if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
int droppedRows = prob->countEmptyRows()-startEmptyRows ;
int droppedColumns = prob->countEmptyCols()-startEmptyColumns ;
std::cout
<< "Leaving dupcol_action::presolve, "
<< droppedRows << " rows, " << droppedColumns << " columns dropped" ;
# if COIN_PRESOLVE_TUNING > 0
std::cout << " in " << thisTime-startTime << "s" ;
# endif
std::cout << "." << std::endl ;
# endif
return (next) ;
}
void dupcol_action::postsolve(CoinPostsolveMatrix *prob) const
{
const action *const actions = actions_;
const int nactions = nactions_;
double *clo = prob->clo_;
double *cup = prob->cup_;
double *sol = prob->sol_;
double *dcost = prob->cost_;
double *colels = prob->colels_;
int *hrow = prob->hrow_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int *link = prob->link_;
double *rcosts = prob->rcosts_;
double tolerance = prob->ztolzb_;
for (const action *f = &actions[nactions-1]; actions<=f; f--) {
int icol = f->ithis; // was fixed
int icol2 = f->ilast; // was kept
dcost[icol] = dcost[icol2];
clo[icol] = f->thislo;
cup[icol] = f->thisup;
clo[icol2] = f->lastlo;
cup[icol2] = f->lastup;
create_col(icol,f->nincol,f->colels,mcstrt,colels,hrow,link,
&prob->free_list_) ;
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
// hincol[icol] = hincol[icol2]; // right? - no - has to match number in create_col
hincol[icol] = f->nincol;
double l_j = f->thislo;
double u_j = f->thisup;
double l_k = f->lastlo;
double u_k = f->lastup;
double x_k_sol = sol[icol2];
PRESOLVE_DETAIL_PRINT(printf("post icol %d %g %g %g icol2 %d %g %g %g\n",
icol,clo[icol],sol[icol],cup[icol],
icol2,clo[icol2],sol[icol2],cup[icol2]));
if (l_j>-PRESOLVE_INF&& x_k_sol-l_j>=l_k-tolerance&&x_k_sol-l_j<=u_k+tolerance) {
// j at lb, leave k
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atLowerBound);
sol[icol] = l_j;
sol[icol2] = x_k_sol - sol[icol];
} else if (u_j<PRESOLVE_INF&& x_k_sol-u_j>=l_k-tolerance&&x_k_sol-u_j<=u_k+tolerance) {
// j at ub, leave k
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atUpperBound);
sol[icol] = u_j;
sol[icol2] = x_k_sol - sol[icol];
} else if (l_k>-PRESOLVE_INF&& x_k_sol-l_k>=l_j-tolerance&&x_k_sol-l_k<=u_j+tolerance) {
// k at lb make j basic
prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
sol[icol2] = l_k;
sol[icol] = x_k_sol - l_k;
prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atLowerBound);
} else if (u_k<PRESOLVE_INF&& x_k_sol-u_k>=l_j-tolerance&&x_k_sol-u_k<=u_j+tolerance) {
// k at ub make j basic
prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
sol[icol2] = u_k;
sol[icol] = x_k_sol - u_k;
prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atUpperBound);
} else {
// both free! superbasic time
sol[icol] = 0.0; // doesn't matter
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::isFree);
}
PRESOLVE_DETAIL_PRINT(printf("post2 icol %d %g icol2 %d %g\n",
icol,sol[icol],
icol2,sol[icol2]));
// row activity doesn't change
// dj of both variables is the same
rcosts[icol] = rcosts[icol2];
// leave until destructor
// deleteAction(f->colels,double *);
# if PRESOLVE_DEBUG > 0
const double ztolzb = prob->ztolzb_;
if (! (clo[icol] - ztolzb <= sol[icol] && sol[icol] <= cup[icol] + ztolzb))
printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol], sol[icol], cup[icol]);
if (! (clo[icol2] - ztolzb <= sol[icol2] && sol[icol2] <= cup[icol2] + ztolzb))
printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol2], sol[icol2], cup[icol2]);
# endif
}
// leave until desctructor
// deleteAction(actions_,action *);
}
dupcol_action::~dupcol_action()
{
for (int i = nactions_-1; i >= 0; --i) {
deleteAction(actions_[i].colels, double *);
}
deleteAction(actions_, action*);
}
/*
Routines for duplicate rows. This is definitely unfinished --- there's no
postsolve action.
*/
const char *duprow_action::name () const
{
return ("duprow_action");
}
// This is just ekkredc4, adapted into the new framework.
/*
I've made minimal changes for compatibility with dupcol: An initial scan to
accumulate rows of interest in sort.
-- lh, 040909 --
*/
const CoinPresolveAction
*duprow_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 *rowels = prob->rowels_;
int *hcol = prob->hcol_;
CoinBigIndex *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
int ncols = prob->ncols_;
int nrows = prob->nrows_;
/*
Scan the rows for candidates, and write the indices into sort. We're not
interested in rows that are empty or prohibited.
Question: Should we exclude singletons, which are useful in other transforms?
Question: Why are we excluding integral columns?
*/
int *sort = new int[nrows] ;
int nlook = 0 ;
for (int i = 0 ; i < nrows ; i++)
{ if (hinrow[i] == 0) continue ;
if (prob->rowProhibited2(i)) continue ;
// sort
CoinSort_2(hcol+mrstrt[i],hcol+mrstrt[i]+hinrow[i],
rowels+mrstrt[i]);
sort[nlook++] = i ; }
if (nlook == 0)
{ delete[] sort ;
return (next) ; }
double * workrow = new double[nrows+1];
double * workcol;
if (!prob->randomNumber_) {
workcol = new double[ncols+1];
coin_init_random_vec(workcol, ncols);
} else {
workcol = prob->randomNumber_;
}
compute_sums(nrows,hinrow,mrstrt,hcol,rowels,workcol,sort,workrow,nlook);
CoinSort_2(workrow,workrow+nlook,sort);
double *rlo = prob->rlo_;
double *rup = prob->rup_;
int nuseless_rows = 0;
bool fixInfeasibility = ((prob->presolveOptions_&0x4000) != 0) ;
bool allowIntersection = ((prob->presolveOptions_&0x10) != 0) ;
double tolerance = prob->feasibilityTolerance_;
double dval = workrow[0];
for (int jj = 1; jj < nlook; jj++) {
if (workrow[jj]==dval) {
int ithis=sort[jj];
int ilast=sort[jj-1];
CoinBigIndex krs = mrstrt[ithis];
CoinBigIndex kre = krs + hinrow[ithis];
if (hinrow[ithis] == hinrow[ilast]) {
int ishift = mrstrt[ilast] - krs;
CoinBigIndex k;
for (k=krs;k<kre;k++) {
if (hcol[k] != hcol[k+ishift] ||
rowels[k] != rowels[k+ishift]) {
break;
}
}
if (k == kre) {
/* now check rhs to see what is what */
double rlo1=rlo[ilast];
double rup1=rup[ilast];
double rlo2=rlo[ithis];
double rup2=rup[ithis];
int idelete=-1;
if (rlo1<=rlo2) {
if (rup2<=rup1) {
/* this is strictly tighter than last */
idelete=ilast;
PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ilast,ithis));
} else if (fabs(rlo1-rlo2)<1.0e-12) {
/* last is strictly tighter than this */
idelete=ithis;
PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ithis,ilast));
// swap so can carry on deleting
sort[jj-1]=ithis;
sort[jj]=ilast;
} else {
if (rup1<rlo2-tolerance&&!fixInfeasibility) {
// infeasible
prob->status_|= 1;
// wrong message - correct if works
prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
prob->messages())
<<ithis
<<rlo[ithis]
<<rup[ithis]
<<CoinMessageEol;
break;
} else if (allowIntersection/*||fabs(rup1-rlo2)<tolerance*/) {
/* overlapping - could merge */
#ifdef CLP_INVESTIGATE7
printf("overlapping duplicate row %g %g, %g %g\n",
rlo1,rup1,rlo2,rup2);
# endif
// pretend this is stricter than last
idelete=ilast;
PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ilast,ithis));
rup[ithis]=rup1;
}
}
} else {
// rlo1>rlo2
if (rup1<=rup2) {
/* last is strictly tighter than this */
idelete=ithis;
PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ithis,ilast));
// swap so can carry on deleting
sort[jj-1]=ithis;
sort[jj]=ilast;
} else {
/* overlapping - could merge */
// rlo1>rlo2
// rup1>rup2
if (rup2<rlo1-tolerance&&!fixInfeasibility) {
// infeasible
prob->status_|= 1;
// wrong message - correct if works
prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
prob->messages())
<<ithis
<<rlo[ithis]
<<rup[ithis]
<<CoinMessageEol;
break;
} else if (allowIntersection/*||fabs(rup2-rlo1)<tolerance*/) {
#ifdef CLP_INVESTIGATE7
printf("overlapping duplicate row %g %g, %g %g\n",
rlo1,rup1,rlo2,rup2);
# endif
// pretend this is stricter than last
idelete=ilast;
PRESOLVE_DETAIL_PRINT(printf("pre_duprow %dR %dR E\n",ilast,ithis));
rlo[ithis]=rlo1;
}
}
}
if (idelete>=0)
sort[nuseless_rows++]=idelete;
}
}
}
dval=workrow[jj];
}
delete[]workrow;
if(workcol != prob->randomNumber_)
delete[]workcol;
if (nuseless_rows) {
# if PRESOLVE_SUMMARY
printf("DUPLICATE ROWS: %d\n", nuseless_rows);
# endif
next = useless_constraint_action::presolve(prob,
sort, nuseless_rows,
next);
}
delete[]sort;
if (prob->tuning_) {
double thisTime=CoinCpuTime();
int droppedRows = prob->countEmptyRows() - startEmptyRows ;
int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
printf("CoinPresolveDuprow(256) - %d rows, %d columns dropped in time %g, total %g\n",
droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
}
return (next);
}
void duprow_action::postsolve(CoinPostsolveMatrix *) const
{
printf("STILL NO POSTSOLVE FOR DUPROW!\n");
abort();
}
/*
Routines for gub rows. This is definitely unfinished --- there's no
postsolve action.
This is potentially called from ClpPresolve and OsiPresolve. Unclear that
it can be backed out --- there's no postsolve.
*/
const char *gubrow_action::name () const
{
return ("gubrow_action");
}
const CoinPresolveAction
*gubrow_action::presolve (CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
double startTime = 0.0;
int droppedElements=0;
int affectedRows=0;
if (prob->tuning_) {
startTime = CoinCpuTime();
}
double *rowels = prob->rowels_;
int *hcol = prob->hcol_;
CoinBigIndex *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
double *colels = prob->colels_ ;
int *hrow = prob->hrow_ ;
CoinBigIndex *mcstrt = prob->mcstrt_ ;
int *hincol = prob->hincol_ ;
int ncols = prob->ncols_;
int nrows = prob->nrows_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
/*
Scan the rows. We're not
interested in rows that are empty or prohibited.
*/
int *which = prob->usefulRowInt_;
int * number = which + nrows;
double * els = prob->usefulRowDouble_;
char * markCol = reinterpret_cast<char *> (prob->usefulColumnInt_);
memset(markCol,0,ncols);
CoinZeroN(els,nrows);
for (int i = 0 ; i < nrows ; i++) {
int nInRow = hinrow[i];
if (nInRow>1 &&!prob->rowProhibited2(i)&&rlo[i]==rup[i]) {
CoinBigIndex rStart = mrstrt[i];
CoinBigIndex k = rStart;
CoinBigIndex rEnd = rStart+nInRow;
double value1=rowels[k];
k++;
for (;k<rEnd;k++) {
if (rowels[k]!=value1)
break;
}
if (k==rEnd) {
// Gub row
int nLook = 0 ;
for (k=rStart;k<rEnd;k++) {
int iColumn = hcol[k];
markCol[iColumn]=1;
CoinBigIndex kk = mcstrt[iColumn];
CoinBigIndex cEnd = kk+hincol[iColumn];
for (;kk<cEnd;kk++) {
int iRow = hrow[kk];
double value = colels[kk];
if (iRow!=i) {
double value2 = els[iRow];
if (value2) {
if (value==value2)
number[iRow]++;
} else {
// first
els[iRow]=value;
number[iRow]=1;
which[nLook++]=iRow;
}
}
}
}
// Now see if any promising
for (int j=0;j<nLook;j++) {
int iRow = which[j];
if (number[iRow]==nInRow) {
// can delete elements and adjust rhs
affectedRows++;
droppedElements += nInRow;
for (CoinBigIndex kk=rStart; kk<rEnd; kk++)
presolve_delete_from_col(iRow,hcol[kk],mcstrt,hincol,hrow,colels) ;
int nInRow2 = hinrow[iRow];
CoinBigIndex rStart2 = mrstrt[iRow];
CoinBigIndex rEnd2 = rStart2+nInRow2;
for (CoinBigIndex kk=rStart2; kk<rEnd2; kk++) {
int iColumn = hcol[kk];
if (markCol[iColumn]==0) {
hcol[rStart2]=iColumn;
rowels[rStart2++]=rowels[kk];
}
}
hinrow[iRow] = nInRow2-nInRow;
if (!hinrow[iRow])
PRESOLVE_REMOVE_LINK(prob->rlink_,iRow) ;
double value =(rlo[i]/value1)*els[iRow];
// correct rhs
if (rlo[iRow]>-1.0e20)
rlo[iRow] -= value;
if (rup[iRow]<1.0e20)
rup[iRow] -= value;
}
els[iRow]=0.0;
}
for (k=rStart;k<rEnd;k++) {
int iColumn = hcol[k];
markCol[iColumn]=0;
}
}
}
}
if (prob->tuning_) {
double thisTime=CoinCpuTime();
printf("CoinPresolveGubrow(1024) - %d elements dropped (%d rows) in time %g, total %g\n",
droppedElements,affectedRows,thisTime-startTime,thisTime-prob->startTime_);
} else if (droppedElements) {
#ifdef CLP_INVESTIGATE
printf("CoinPresolveGubrow(1024) - %d elements dropped (%d rows)\n",
droppedElements,affectedRows);
#endif
}
return (next);
}
void gubrow_action::postsolve(CoinPostsolveMatrix *) const
{
printf("STILL NO POSTSOLVE FOR GUBROW!\n");
abort();
}
/*
Routines for two by two blocks. This is definitely unfinished --- there's no
postsolve action.
*/
const char *twoxtwo_action::name () const
{
return ("twoxtwo_action");
}
const CoinPresolveAction
*twoxtwo_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();
}
// maximum number of records
action * boundRecords = new action[(prob->nrows_+1)>>1];
int nactions=0;
// column-major representation
const int ncols = prob->ncols_ ;
const CoinBigIndex *const mcstrt = prob->mcstrt_ ;
const int *const hincol = prob->hincol_ ;
const int *const hrow = prob->hrow_ ;
const double * colels = prob->colels_ ;
double * cost = prob->cost_ ;
// column type, bounds, solution, and status
const unsigned char *const integerType = prob->integerType_ ;
double * clo = prob->clo_ ;
double * cup = prob->cup_ ;
// row-major representation
//const int nrows = prob->nrows_ ;
const CoinBigIndex *const mrstrt = prob->mrstrt_ ;
const int *const hinrow = prob->hinrow_ ;
const int *const hcol = prob->hcol_ ;
const double * rowels = prob->rowels_ ;
// row bounds
double * rlo = prob->rlo_ ;
double * rup = prob->rup_ ;
// tolerances
//const double ekkinf2 = PRESOLVE_SMALL_INF ;
//const double ekkinf = ekkinf2*1.0e8 ;
//const double ztolcbarj = prob->ztoldj_ ;
//const CoinRelFltEq relEq(prob->ztolzb_) ;
double bound[2];
double alpha[2]={0.0,0.0};
double offset=0.0;
for (int icol=0;icol<ncols;icol++) {
if (hincol[icol]==2) {
CoinBigIndex start=mcstrt[icol];
int row0 = hrow[start];
if (hinrow[row0]!=2)
continue;
int row1 = hrow[start+1];
if (hinrow[row1]!=2)
continue;
double element0 = colels[start];
double rowUpper0=rup[row0];
bool swapSigns0=false;
if (rlo[row0]>-1.0e30) {
if (rup[row0]>1.0e30) {
swapSigns0=true;
rowUpper0=-rlo[row0];
element0=-element0;
} else {
// range or equality
continue;
}
} else if (rup[row0]>1.0e30) {
// free
continue;
}
#if 0
// skip here for speed
// skip if no cost (should be able to get rid of)
if (!cost[icol]) {
PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n",icol));
continue;
}
// skip if negative cost for now
if (cost[icol]<0.0) {
PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
continue;
}
#endif
double element1 = colels[start+1];
double rowUpper1=rup[row1];
bool swapSigns1=false;
if (rlo[row1]>-1.0e30) {
if (rup[row1]>1.0e30) {
swapSigns1=true;
rowUpper1=-rlo[row1];
element1=-element1;
} else {
// range or equality
continue;
}
} else if (rup[row1]>1.0e30) {
// free
continue;
}
double lowerX=clo[icol];
double upperX=cup[icol];
int otherCol=-1;
CoinBigIndex startRow=mrstrt[row0];
for (CoinBigIndex j=startRow;j<startRow+2;j++) {
int jcol=hcol[j];
if (jcol!=icol) {
alpha[0]=swapSigns0 ? -rowels[j] :rowels[j];
otherCol=jcol;
}
}
startRow=mrstrt[row1];
bool possible=true;
for (CoinBigIndex j=startRow;j<startRow+2;j++) {
int jcol=hcol[j];
if (jcol!=icol) {
if (jcol==otherCol) {
alpha[1]=swapSigns1 ? -rowels[j] :rowels[j];
} else {
possible=false;
}
}
}
if (possible) {
// skip if no cost (should be able to get rid of)
if (!cost[icol]) {
PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n",icol));
continue;
}
// skip if negative cost for now
if (cost[icol]<0.0) {
PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
continue;
}
bound[0]=clo[otherCol];
bound[1]=cup[otherCol];
double lowestLowest=COIN_DBL_MAX;
double highestLowest=-COIN_DBL_MAX;
double lowestHighest=COIN_DBL_MAX;
double highestHighest=-COIN_DBL_MAX;
int binding0=0;
int binding1=0;
for (int k=0;k<2;k++) {
bool infLow0=false;
bool infLow1=false;
double sum0=0.0;
double sum1=0.0;
double value=bound[k];
if (fabs(value)<1.0e30) {
sum0+=alpha[0]*value;
sum1+=alpha[1]*value;
} else {
if (alpha[0]>0.0) {
if (value<0.0)
infLow0 =true;
} else if (alpha[0]<0.0) {
if (value>0.0)
infLow0 =true;
}
if (alpha[1]>0.0) {
if (value<0.0)
infLow1 =true;
} else if (alpha[1]<0.0) {
if (value>0.0)
infLow1 =true;
}
}
/* Got sums
*/
double thisLowest0=-COIN_DBL_MAX;
double thisHighest0=COIN_DBL_MAX;
if (element0>0.0) {
// upper bound unless inf&2 !=0
if (!infLow0)
thisHighest0 = (rowUpper0-sum0)/element0;
} else {
// lower bound unless inf&2 !=0
if (!infLow0)
thisLowest0 = (rowUpper0-sum0)/element0;
}
double thisLowest1=-COIN_DBL_MAX;
double thisHighest1=COIN_DBL_MAX;
if (element1>0.0) {
// upper bound unless inf&2 !=0
if (!infLow1)
thisHighest1 = (rowUpper1-sum1)/element1;
} else {
// lower bound unless inf&2 !=0
if (!infLow1)
thisLowest1 = (rowUpper1-sum1)/element1;
}
if (thisLowest0>thisLowest1+1.0e-12) {
if (thisLowest0>lowerX+1.0e-12)
binding0|= 1<<k;
} else if (thisLowest1>thisLowest0+1.0e-12) {
if (thisLowest1>lowerX+1.0e-12)
binding1|= 1<<k;
thisLowest0=thisLowest1;
}
if (thisHighest0<thisHighest1-1.0e-12) {
if (thisHighest0<upperX-1.0e-12)
binding0|= 1<<k;
} else if (thisHighest1<thisHighest0-1.0e-12) {
if (thisHighest1<upperX-1.0e-12)
binding1|= 1<<k;
thisHighest0=thisHighest1;
}
lowestLowest=CoinMin(lowestLowest,thisLowest0);
highestHighest=CoinMax(highestHighest,thisHighest0);
lowestHighest=CoinMin(lowestHighest,thisHighest0);
highestLowest=CoinMax(highestLowest,thisLowest0);
}
// see if any good
//#define PRINT_VALUES
if (!binding0||!binding1) {
PRESOLVE_DETAIL_PRINT(printf("Row redundant for column %d\n",icol));
} else {
PRESOLVE_DETAIL_PRINT(printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
icol,lowerX,upperX,lowestLowest,lowestHighest,
highestLowest,highestHighest));
// if integer adjust
if (integerType[icol]) {
lowestLowest=ceil(lowestLowest-1.0e-5);
highestLowest=ceil(highestLowest-1.0e-5);
lowestHighest=floor(lowestHighest+1.0e-5);
highestHighest=floor(highestHighest+1.0e-5);
}
// if costed may be able to adjust
if (cost[icol]>=0.0) {
if (highestLowest<upperX&&highestLowest>=lowerX&&highestHighest<1.0e30) {
highestHighest=CoinMin(highestHighest,highestLowest);
}
}
if (cost[icol]<=0.0) {
if (lowestHighest>lowerX&&lowestHighest<=upperX&&lowestHighest>-1.0e30) {
lowestLowest=CoinMax(lowestLowest,lowestHighest);
}
}
#if 1
if (lowestLowest>lowerX+1.0e-8) {
PRESOLVE_DETAIL_PRINT(printf("Can increase lower bound on %d from %g to %g\n",
icol,lowerX,lowestLowest));
lowerX=lowestLowest;
}
if (highestHighest<upperX-1.0e-8) {
PRESOLVE_DETAIL_PRINT(printf("Can decrease upper bound on %d from %g to %g\n",
icol,upperX,highestHighest));
upperX=highestHighest;
}
#endif
// see if we can move costs
double xValue;
double yValue0;
double yValue1;
double newLower=COIN_DBL_MAX;
double newUpper=-COIN_DBL_MAX;
double costEqual;
double slope[2];
assert (binding0+binding1==3);
// get where equal
xValue=(rowUpper0*element1-rowUpper1*element0)/(alpha[0]*element1-alpha[1]*element0);
yValue0=(rowUpper0-xValue*alpha[0])/element0;
yValue1=(rowUpper1-xValue*alpha[1])/element1;
newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
double xValueEqual=xValue;
double yValueEqual=yValue0;
costEqual = xValue*cost[otherCol]+yValueEqual*cost[icol];
if (binding0==1) {
// take x 1.0 down
double x=xValue-1.0;
double y=(rowUpper0-x*alpha[0])/element0;
double costTotal = x*cost[otherCol]+y*cost[icol];
slope[0] = costEqual-costTotal;
// take x 1.0 up
x=xValue+1.0;
y=(rowUpper1-x*alpha[1])/element0;
costTotal = x*cost[otherCol]+y*cost[icol];
slope[1] = costTotal-costEqual;
} else {
// take x 1.0 down
double x=xValue-1.0;
double y=(rowUpper1-x*alpha[1])/element0;
double costTotal = x*cost[otherCol]+y*cost[icol];
slope[1] = costEqual-costTotal;
// take x 1.0 up
x=xValue+1.0;
y=(rowUpper0-x*alpha[0])/element0;
costTotal = x*cost[otherCol]+y*cost[icol];
slope[0] = costTotal-costEqual;
}
PRESOLVE_DETAIL_PRINT(printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1)));
PRESOLVE_DETAIL_PRINT(printf("Cost at equality %g for constraint 0 slope %g for constraint 1 slope %g\n",
costEqual,slope[0],slope[1]));
xValue=bound[0];
yValue0=(rowUpper0-xValue*alpha[0])/element0;
yValue1=(rowUpper1-xValue*alpha[1])/element1;
PRESOLVE_DETAIL_PRINT(printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1)));
newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
// cost>0 so will be at lower
//double yValueAtBound0=newLower;
newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
xValue=bound[1];
yValue0=(rowUpper0-xValue*alpha[0])/element0;
yValue1=(rowUpper1-xValue*alpha[1])/element1;
PRESOLVE_DETAIL_PRINT(printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1)));
newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
// cost>0 so will be at lower
//double yValueAtBound1=newLower;
newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
lowerX=CoinMax(lowerX,newLower-1.0e-12*fabs(newLower));
upperX=CoinMin(upperX,newUpper+1.0e-12*fabs(newUpper));
// Now make duplicate row
// keep row 0 so need to adjust costs so same
PRESOLVE_DETAIL_PRINT(printf("Costs for x %g,%g,%g are %g,%g,%g\n",
xValueEqual-1.0,xValueEqual,xValueEqual+1.0,
costEqual-slope[0],costEqual,costEqual+slope[1]));
double costOther=cost[otherCol]+slope[1];
double costThis=cost[icol]+slope[1]*(element0/alpha[0]);
xValue=xValueEqual;
yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
double thisOffset=costEqual-(costOther*xValue+costThis*yValue0);
offset += thisOffset;
PRESOLVE_DETAIL_PRINT(printf("new cost at equal %g\n",costOther*xValue+costThis*yValue0+thisOffset));
xValue=xValueEqual-1.0;
yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
PRESOLVE_DETAIL_PRINT(printf("new cost at -1 %g\n",costOther*xValue+costThis*yValue0+thisOffset));
assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual-slope[0]))<1.0e-5);
xValue=xValueEqual+1.0;
yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
PRESOLVE_DETAIL_PRINT(printf("new cost at +1 %g\n",costOther*xValue+costThis*yValue0+thisOffset));
assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual+slope[1]))<1.0e-5);
action & boundRecord = boundRecords[nactions++] ;
boundRecord.row=row1;
boundRecord.col=icol;
boundRecord.othercol=otherCol;
boundRecord.lbound_row=rlo[row1];
boundRecord.ubound_row=rup[row1];
boundRecord.lbound_col=clo[icol];
boundRecord.ubound_col=cup[icol];
boundRecord.cost_col=cost[icol];
boundRecord.cost_othercol=cost[otherCol];
cost[otherCol] = costOther;
cost[icol] = costThis;
clo[icol]=lowerX;
cup[icol]=upperX;
// make row useless
rlo[row1]=-COIN_DBL_MAX;
rup[row1]=COIN_DBL_MAX;
}
}
}
}
if (nactions) {
# if PRESOLVE_SUMMARY
printf("Cost offset %g - from %d blocks\n",offset,nactions);
printf("TWO by TWO blocks: %d - offset %g\n", nactions,offset);
# endif
action * actions = new action[nactions];
memcpy(actions,boundRecords,nactions*sizeof(action));
next = new twoxtwo_action(nactions,actions,next);
int *sort = prob->usefulColumnInt_;
for (int i=0;i<nactions;i++)
sort[i]=boundRecords[i].row;
next = useless_constraint_action::presolve(prob,
sort, nactions,
next);
// adjust offset
prob->change_bias(offset);
}
delete [] boundRecords;
if (prob->tuning_) {
double thisTime=CoinCpuTime();
int droppedRows = prob->countEmptyRows() - startEmptyRows ;
int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
printf("CoinPresolveTwoxtwo(2048) - %d rows, %d columns dropped in time %g, total %g\n",
droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
}
return (next);
}
void twoxtwo_action::postsolve(CoinPostsolveMatrix * prob) const
{
const CoinBigIndex *const mcstrt = prob->mcstrt_ ;
const int *const hincol = prob->hincol_ ;
const int *const hrow = prob->hrow_ ;
const double * colels = prob->colels_ ;
int *link = prob->link_;
double * cost = prob->cost_ ;
// column type, bounds, solution, and status
//const unsigned char *const integerType = prob->integerType_ ;
double * clo = prob->clo_ ;
double * cup = prob->cup_ ;
// row bounds
double * rlo = prob->rlo_ ;
double * rup = prob->rup_ ;
double *sol = prob->sol_;
double *rcosts = prob->rcosts_;
double * rowacts = prob->acts_;
double * dual = prob->rowduals_;
double tolerance = prob->ztolzb_;
const double maxmin = prob->maxmin_;
for (int iAction=0;iAction<nactions_;iAction++) {
const action & boundRecord = actions_[iAction];
int row1=boundRecord.row;
int icol=boundRecord.col;
int otherCol=boundRecord.othercol;
CoinBigIndex start=mcstrt[icol];
CoinBigIndex nextEl = link[start];
int row0;
// first is otherCol
double els0[2]={0.0,0.0};
double els1[2]={0.0,0.0};
if (hrow[start]==row1) {
row0=hrow[nextEl];
els0[1]=colels[nextEl];
els1[1]=colels[start];
} else {
row0=hrow[start];
els0[1]=colels[start];
els1[1]=colels[nextEl];
}
nextEl=mcstrt[otherCol];
for (CoinBigIndex j=0;j<hincol[otherCol];j++) {
if (hrow[nextEl]==row0)
els0[0]=colels[nextEl];
else if (hrow[nextEl]==row1)
els1[0]=colels[nextEl];
nextEl=link[nextEl];
}
prob->setRowStatus(row1,CoinPrePostsolveMatrix::basic);
// put stuff back
rlo[row1]=boundRecord.lbound_row;
rup[row1]=boundRecord.ubound_row;
clo[icol]=boundRecord.lbound_col;
cup[icol]=boundRecord.ubound_col;
double oldCost=cost[icol];
//double oldOtherCost=cost[otherCol];
cost[icol]=boundRecord.cost_col;
cost[otherCol]=boundRecord.cost_othercol;
double els0real[2];
double els1real[2];
els0real[0]=els0[0];
els0real[1]=els0[1];
els1real[0]=els1[0];
els1real[1]=els1[1];
// make <= rows
double rowUpper0=rup[row0];
if (rlo[row0]>-1.0e30) {
rowUpper0=-rlo[row0];
els0[0]=-els0[0];
els0[1]=-els0[1];
}
double rowUpper1=rup[row1];
bool swapSigns1=false;
if (rlo[row1]>-1.0e30) {
swapSigns1=true;
rowUpper1=-rlo[row1];
els1[0]=-els1[0];
els1[1]=-els1[1];
}
// compute feasible value for icol
double valueOther=sol[otherCol];
double value;
// first see if at bound is OK
bool lowerBoundPossible = clo[icol]>-1.0e30;
value=clo[icol];
if (lowerBoundPossible) {
double value0 = els0[0]*valueOther + els0[1]*value;
if (value0>rowUpper0+tolerance)
lowerBoundPossible=false;
double value1 = els1[0]*valueOther + els1[1]*value;
if (value1>rowUpper1+tolerance)
lowerBoundPossible=false;
}
bool upperBoundPossible = cup[icol]<1.0e30;
value=cup[icol];
if (upperBoundPossible) {
double value0 = els0[0]*valueOther + els0[1]*value;
if (value0>rowUpper0+tolerance)
upperBoundPossible=false;
double value1 = els1[0]*valueOther + els1[1]*value;
if (value1>rowUpper1+tolerance)
upperBoundPossible=false;
}
if (lowerBoundPossible&&cost[icol]>=0.0) {
// set to lower bound
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atLowerBound);
sol[icol]=clo[icol];
rcosts[icol]=maxmin*cost[icol]-dual[row0]*els0real[1];
} else if (upperBoundPossible&&cost[icol]<=0.0) {
// set to upper bound
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atUpperBound);
sol[icol]=cup[icol];
rcosts[icol]=maxmin*cost[icol]-dual[row0]*els0real[1];
} else {
// need to make basic
// we shouldn't get here (at present) if zero cost
assert (cost[icol]);
double value0 = (rowUpper0-els0[0]*valueOther)/els0[1];
double value1 = (rowUpper1-els1[0]*valueOther)/els1[1];
//bool binding0=true;
double value;
if (cost[icol]>0) {
if (value0>value1) {
value=value0;
} else {
value=value1;
//binding0=false;
}
} else {
if (value0<value1) {
value=value0;
} else {
value=value1;
//binding0=false;
}
}
sol[icol]=value;
#if 0
printf("row %d status %d, row %d status %d, col %d status %d, col %d status %d - binding0 %c\n",
row0,prob->getRowStatus(row0),
row1,prob->getRowStatus(row1),
otherCol,prob->getColumnStatus(otherCol),
icol,prob->getColumnStatus(icol),binding0 ? 'T' : 'F');
#endif
if (prob->getColumnStatus(icol)==CoinPrePostsolveMatrix::basic) {
//printf("col %d above was basic\n",icol);
if (prob->getRowStatus(row0)!=CoinPrePostsolveMatrix::basic) {
// adjust dual
dual[row0]=maxmin*((cost[icol]-oldCost)/els0real[1]);
}
continue;
}
//if (binding0)
//printf("Says row0 %d binding?\n",row0);
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::basic);
rcosts[icol]=0.0;
//printf("row1 %d taken out of basis\n",row1);
if (!swapSigns1) {
prob->setRowStatus(row1,CoinPrePostsolveMatrix::atUpperBound);
rowacts[row1]=rup[row1];
} else {
prob->setRowStatus(row1,CoinPrePostsolveMatrix::atLowerBound);
rowacts[row1]=rlo[row1];
}
dual[row1]=maxmin*((cost[icol]-oldCost)/els1real[1]);
if (iAction==-1)
abort();
}
}
}