limp-cbc-0.3.2.0: cbits/coin/CoinPresolveEmpty.cpp
/* $Id: CoinPresolveEmpty.cpp 1607 2013-07-16 09:01:29Z 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>
#include "CoinFinite.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveEmpty.hpp"
#include "CoinMessage.hpp"
#if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
#include "CoinPresolvePsdebug.hpp"
#endif
/* \file
Routines to remove/reinsert empty columns and rows.
*/
/*
Physically remove empty columns, compressing mcstrt and hincol. The major
side effect is that columns are renumbered, thus clink_ is no longer valid
and must be rebuilt. And we're totally inconsistent with the row-major
representation.
It's necessary to rebuild clink_ in order to do direct conversion of a
CoinPresolveMatrix to a CoinPostsolveMatrix by transferring the data arrays.
Without clink_, it's impractical to build link_ to match the transferred bulk
storage.
*/
const CoinPresolveAction
*drop_empty_cols_action::presolve (CoinPresolveMatrix *prob,
const int *ecols,
int necols,
const CoinPresolveAction *next)
{
# if PRESOLVE_CONSISTENCY > 0
presolve_links_ok(prob) ;
presolve_consistent(prob) ;
# endif
const int n_orig = prob->ncols_ ;
CoinBigIndex *mcstrt = prob->mcstrt_ ;
int *hincol = prob->hincol_ ;
presolvehlink *clink = prob->clink_ ;
double *clo = prob->clo_ ;
double *cup = prob->cup_ ;
double *cost = prob->cost_ ;
const double ztoldj = prob->ztoldj_ ;
unsigned char *integerType = prob->integerType_ ;
int *originalColumn = prob->originalColumn_ ;
const double maxmin = prob->maxmin_ ;
double *sol = prob->sol_ ;
unsigned char *colstat = prob->colstat_ ;
action *actions = new action[necols] ;
int *colmapping = new int [n_orig+1] ;
CoinZeroN(colmapping,n_orig) ;
// More like `ignore infeasibility'
bool fixInfeasibility = ((prob->presolveOptions_&0x4000) != 0) ;
/*
Open a loop to walk the list of empty columns. Mark them as empty in
colmapping.
*/
for (int ndx = necols-1 ; ndx >= 0 ; ndx--) {
const int j = ecols[ndx] ;
colmapping[j] = -1 ;
/*
Groom bounds on integral variables. Check for previously undetected
infeasibility unless the user wants to ignore it. If we find it, quit.
*/
double &lj = clo[j] ;
double &uj = cup[j] ;
if (integerType[j]) {
lj = ceil(lj-1.0e-9) ;
uj = floor(uj+1.0e-9) ;
if (lj > uj && !fixInfeasibility) {
prob->status_ |= 1 ;
prob->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
prob->messages())
<< j << lj << uj << CoinMessageEol ;
break ;
}
}
/*
Load up the postsolve action with the index and rim vector components.
*/
action &e = actions[ndx] ;
e.jcol = j ;
e.clo = lj ;
e.cup = uj ;
e.cost = cost[j] ;
/*
There are no more constraints on this variable so we had better be able to
compute the answer now. Try to make it nonbasic at bound. If we're
unbounded, say so and quit.
*/
if (fabs(cost[j]) < ztoldj) cost[j] = 0.0 ;
if (cost[j] == 0.0) {
if (-PRESOLVE_INF < lj)
e.sol = lj ;
else if (uj < PRESOLVE_INF)
e.sol = uj ;
else
e.sol = 0.0 ;
} else if (cost[j]*maxmin > 0.0) {
if (-PRESOLVE_INF < lj)
e.sol = lj ;
else {
prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDB,
prob->messages())
<< j << CoinMessageEol ;
prob->status_ |= 2 ;
break ;
}
} else {
if (uj < PRESOLVE_INF)
e.sol = uj ;
else {
prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDA,
prob->messages())
<< j << CoinMessageEol ;
prob->status_ |= 2 ;
break ;
}
}
# if PRESOLVE_DEBUG > 2
if (e.sol*cost[j]) {
std::cout
<< " non-zero cost " << cost[j] << " for empty col " << j << "."
<< std::endl ;
}
# endif
prob->change_bias(e.sol*cost[j]) ;
}
/*
No sense doing the actual work of compression unless we're still feasible
and bounded. If we are, start out by compressing out the entries associated
with empty columns. Empty columns are nonzero in colmapping.
*/
if (prob->status_ == 0) {
int n_compressed = 0 ;
for (int ndx = 0 ; ndx < n_orig ; ndx++) {
if (!colmapping[ndx]) {
mcstrt[n_compressed] = mcstrt[ndx] ;
hincol[n_compressed] = hincol[ndx] ;
clo[n_compressed] = clo[ndx] ;
cup[n_compressed] = cup[ndx] ;
cost[n_compressed] = cost[ndx] ;
if (sol) {
sol[n_compressed] = sol[ndx] ;
colstat[n_compressed] = colstat[ndx] ;
}
integerType[n_compressed] = integerType[ndx] ;
originalColumn[n_compressed] = originalColumn[ndx] ;
colmapping[ndx] = n_compressed++ ;
}
}
mcstrt[n_compressed] = mcstrt[n_orig] ;
colmapping[n_orig] = n_compressed ;
/*
Rebuild clink_. At this point, all empty columns are linked out, so the
only columns left are columns that are to be saved, hence available in
colmapping. All we need to do is walk clink_ and write the new entries
into a new array.
*/
presolvehlink *newclink = new presolvehlink [n_compressed+1] ;
for (int oldj = n_orig ; oldj >= 0 ; oldj = clink[oldj].pre) {
presolvehlink &oldlnk = clink[oldj] ;
int newj = colmapping[oldj] ;
assert(newj >= 0 && newj <= n_compressed) ;
presolvehlink &newlnk = newclink[newj] ;
if (oldlnk.suc >= 0) {
newlnk.suc = colmapping[oldlnk.suc] ;
} else {
newlnk.suc = NO_LINK ;
}
if (oldlnk.pre >= 0) {
newlnk.pre = colmapping[oldlnk.pre] ;
} else {
newlnk.pre = NO_LINK ;
}
}
delete [] clink ;
prob->clink_ = newclink ;
prob->ncols_ = n_compressed ;
}
delete [] colmapping ;
# if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
presolve_links_ok(prob,true,false) ;
presolve_check_sol(prob) ;
presolve_check_nbasic(prob) ;
# endif
return (new drop_empty_cols_action(necols,actions,next)) ;
}
/*
The top-level method scans the matrix for empty columns and calls a worker
routine to do the heavy lifting.
NOTE: At the end of this routine, the column- and row-major representations
are not consistent. Empty columns have been compressed out,
effectively renumbering the columns.
*/
const CoinPresolveAction
*drop_empty_cols_action::presolve (CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
# if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
# if PRESOLVE_DEBUG > 0
std::cout << "Entering drop_empty_cols_action::presolve." << std::endl ;
# endif
presolve_consistent(prob) ;
presolve_links_ok(prob) ;
presolve_check_sol(prob) ;
presolve_check_nbasic(prob) ;
# endif
const int *hincol = prob->hincol_ ;
int ncols = prob->ncols_ ;
int nempty = 0 ;
int *empty = new int [ncols] ;
CoinBigIndex nelems2 = 0 ;
// count empty cols
for (int i = 0 ; i < ncols ; i++) {
nelems2 += hincol[i] ;
if (hincol[i] == 0&&!prob->colProhibited2(i)) {
# if PRESOLVE_DEBUG > 1
if (nempty == 0)
std::cout << "UNUSED COLS:" ;
else
if (i < 100 && nempty%25 == 0)
std::cout << std::endl ;
else
if (i >= 100 && i < 1000 && nempty%19 == 0)
std::cout << std::endl ;
else
if (i >= 1000 && nempty%15 == 0)
std::cout << std::endl ;
std::cout << " " << i ;
# endif
empty[nempty++] = i;
}
}
prob->nelems_ = nelems2 ;
if (nempty)
next = drop_empty_cols_action::presolve(prob,empty,nempty,next) ;
delete [] empty ;
# if PRESOLVE_DEBUG > 0 || COIN_PRESOLVE_TUNING > 0
std::cout << "Leaving drop_empty_cols_action::presolve" ;
if (nempty) std::cout << ", dropped " << nempty << " columns" ;
std::cout << "." << std::endl ;
# endif
return (next);
}
/*
Reintroduce empty columns dropped at the end of presolve.
*/
void drop_empty_cols_action::postsolve(CoinPostsolveMatrix *prob) const
{
const int nactions = nactions_ ;
const action *const actions = actions_ ;
int ncols = prob->ncols_ ;
# if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
# if PRESOLVE_DEBUG > 0
std::cout
<< "Entering drop_empty_cols_action::postsolve, initial system "
<< prob->nrows_ << "x" << ncols << ", " << nactions
<< " columns to restore." << std::endl ;
# endif
char *cdone = prob->cdone_ ;
presolve_check_sol(prob,2,2,2) ;
presolve_check_nbasic(prob) ;
# endif
CoinBigIndex *colStarts = prob->mcstrt_ ;
int *colLengths = prob->hincol_ ;
double *clo = prob->clo_ ;
double *cup = prob->cup_ ;
double *sol = prob->sol_ ;
double *cost = prob->cost_ ;
double *rcosts = prob->rcosts_ ;
unsigned char *colstat = prob->colstat_ ;
const double maxmin = prob->maxmin_ ;
/*
Set up a mapping vector, coded 0 for existing columns, -1 for columns we're
about to reintroduce.
*/
int ncols2 = ncols+nactions ;
int *colmapping = new int [ncols2] ;
CoinZeroN(colmapping,ncols2) ;
for (int ndx = 0 ; ndx < nactions ; ndx++) {
const action *e = &actions[ndx] ;
int j = e->jcol ;
colmapping[j] = -1 ;
}
/*
Working back from the highest index, expand the existing ncols columns over
the full range ncols2, leaving holes for the columns we want to reintroduce.
*/
for (int j = ncols2-1 ; j >= 0 ; j--) {
if (!colmapping[j]) {
ncols-- ;
colStarts[j] = colStarts[ncols] ;
colLengths[j] = colLengths[ncols] ;
clo[j] = clo[ncols] ;
cup[j] = cup[ncols] ;
cost[j] = cost[ncols] ;
if (sol) sol[j] = sol[ncols] ;
if (rcosts) rcosts[j] = rcosts[ncols] ;
if (colstat) colstat[j] = colstat[ncols] ;
# if PRESOLVE_DEBUG > 0
cdone[j] = cdone[ncols] ;
# endif
}
}
assert (!ncols) ;
delete [] colmapping ;
/*
Reintroduce the dropped columns.
*/
for (int ndx = 0 ; ndx < nactions ; ndx++) {
const action *e = &actions[ndx] ;
int j = e->jcol ;
colLengths[j] = 0 ;
colStarts[j] = NO_LINK ;
clo[j] = e->clo ;
cup[j] = e->cup ;
cost[j] = e->cost ;
if (sol) sol[j] = e->sol ;
if (rcosts) rcosts[j] = maxmin*cost[j] ;
if (colstat) prob->setColumnStatusUsingValue(j) ;
# if PRESOLVE_DEBUG > 0
cdone[j] = DROP_COL ;
# if PRESOLVE_DEBUG > 1
std::cout
<< " restoring col " << j << ", lb = " << clo[j] << ", ub = " << cup[j]
<< std::endl ;
# endif
# endif
}
prob->ncols_ += nactions ;
# if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
presolve_check_threads(prob) ;
presolve_check_sol(prob,2,2,2) ;
presolve_check_nbasic(prob) ;
# if PRESOLVE_DEBUG > 0
std::cout << "Leaving drop_empty_cols_action::postsolve, system "
<< prob->nrows_ << "x" << prob->ncols_ << "." << std::endl ;
# endif
# endif
}
const CoinPresolveAction
*drop_empty_rows_action::presolve (CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
# if PRESOLVE_DEBUG > 0
std::cout << "Entering drop_empty_rows_action::presolve." << std::endl ;
# endif
int ncols = prob->ncols_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int *hrow = prob->hrow_;
int nrows = prob->nrows_;
// This is done after row copy needed
//int *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
//int *hcol = prob->hcol_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
unsigned char *rowstat = prob->rowstat_;
double *acts = prob->acts_;
int * originalRow = prob->originalRow_;
//presolvehlink *rlink = prob->rlink_;
bool fixInfeasibility = ((prob->presolveOptions_&0x4000) != 0) ;
// Relax tolerance
double tolerance = 10.0*prob->feasibilityTolerance_;
int i;
int nactions = 0;
for (i=0; i<nrows; i++)
if (hinrow[i] == 0)
nactions++;
/*
Bail out if there's nothing to be done.
*/
if (nactions == 0) {
# if PRESOLVE_DEBUG > 0
std::cout << "Leaving drop_empty_rows_action::presolve." << std::endl ;
# endif
return (next) ;
}
/*
Work to do.
*/
action *actions = new action[nactions];
int * rowmapping = new int [nrows];
nactions = 0;
int nrows2=0;
for (i=0; i<nrows; i++) {
if (hinrow[i] == 0) {
action &e = actions[nactions];
# if PRESOLVE_DEBUG > 1
if (nactions == 0)
std::cout << "UNUSED ROWS:" ;
else
if (i < 100 && nactions%25 == 0)
std::cout << std::endl ;
else
if (i >= 100 && i < 1000 && nactions%19 == 0)
std::cout << std::endl ;
else
if (i >= 1000 && nactions%15 == 0)
std::cout << std::endl ;
std::cout << " " << i ;
# endif
nactions++;
if (rlo[i] > 0.0 || rup[i] < 0.0) {
if ((rlo[i]<=tolerance &&
rup[i]>=-tolerance)||fixInfeasibility) {
rlo[i]=0.0;
rup[i]=0.0;
} else {
prob->status_|= 1;
prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
prob->messages())
<<i
<<rlo[i]
<<rup[i]
<<CoinMessageEol;
break;
}
}
e.row = i;
e.rlo = rlo[i];
e.rup = rup[i];
rowmapping[i]=-1;
} else {
// move down - we want to preserve order
rlo[nrows2]=rlo[i];
rup[nrows2]=rup[i];
originalRow[nrows2]=i;
if (acts) {
acts[nrows2]=acts[i];
rowstat[nrows2]=rowstat[i];
}
rowmapping[i]=nrows2++;
}
}
# if PRESOLVE_DEBUG > 1
std::cout << std::endl ;
# endif
// remap matrix
for (i=0;i<ncols;i++) {
int j;
for (j=mcstrt[i];j<mcstrt[i]+hincol[i];j++)
hrow[j] = rowmapping[hrow[j]];
}
delete [] rowmapping;
prob->nrows_ = nrows2;
next = new drop_empty_rows_action(nactions,actions,next) ;
# if PRESOLVE_DEBUG > 0 || PRESOLVE_CONSISTENCY > 0
presolve_check_nbasic(prob) ;
# if PRESOLVE_DEBUG > 0
std::cout << "Leaving drop_empty_rows_action::presolve" ;
if (nactions) std::cout << ", dropped " << nactions << " rows" ;
std::cout << "." << std::endl ;
# endif
# endif
return (next) ;
}
void drop_empty_rows_action::postsolve(CoinPostsolveMatrix *prob) const
{
const int nactions = nactions_ ;
const action *const actions = actions_ ;
int ncols = prob->ncols_ ;
int nrows0 = prob->nrows0_ ;
int nrows = prob->nrows_ ;
CoinBigIndex *mcstrt = prob->mcstrt_ ;
int *hincol = prob->hincol_ ;
int *hrow = prob->hrow_ ;
double *rlo = prob->rlo_ ;
double *rup = prob->rup_ ;
unsigned char *rowstat = prob->rowstat_ ;
double *rowduals = prob->rowduals_ ;
double *acts = prob->acts_ ;
# if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
# if PRESOLVE_DEBUG > 0
std::cout
<< "Entering drop_empty_rows_action::postsolve, initial system "
<< nrows << "x" << ncols << ", " << nactions
<< " rows to restore." << std::endl ;
# endif
char *rdone = prob->rdone_ ;
presolve_check_sol(prob,2,2,2) ;
presolve_check_nbasic(prob) ;
# endif
/*
Process the array of actions and mark rowmapping[i] if constraint i was
eliminated in presolve.
*/
int *rowmapping = new int [nrows0] ;
CoinZeroN(rowmapping,nrows0) ;
for (int k = 0 ; k < nactions ; k++) {
const action *e = &actions[k] ;
int i = e->row ;
rowmapping[i] = -1 ;
}
/*
Now walk the vectors for row bounds, activity, duals, and status. Expand
the existing entries in 0..(nrows-1) to occupy 0..(nrows0-1), leaving
holes for the rows we're about to reintroduce.
*/
for (int i = nrows0-1 ; i >= 0 ; i--) {
if (!rowmapping[i]) {
nrows-- ;
rlo[i] = rlo[nrows] ;
rup[i] = rup[nrows] ;
acts[i] = acts[nrows] ;
rowduals[i] = rowduals[nrows] ;
if (rowstat)
rowstat[i] = rowstat[nrows] ;
# if PRESOLVE_DEBUG > 0
rdone[i] = rdone[nrows] ;
# endif
}
}
assert (!nrows) ;
/*
Rewrite rowmapping so that it maps presolved row indices to row indices in
the restored matrix.
*/
for (int i = 0 ; i < nrows0 ; i++) {
if (!rowmapping[i])
rowmapping[nrows++] = i ;
}
/*
Now walk the row index array for each column, rewriting the row indices so
they are correct for the restored matrix.
*/
for (int j = 0 ; j < ncols ; j++) {
const CoinBigIndex &start = mcstrt[j] ;
const CoinBigIndex &end = start+hincol[j] ;
for (CoinBigIndex k = start ; k < end ; k++) {
hrow[k] = rowmapping[hrow[k]] ;
}
}
delete [] rowmapping;
/*
And reintroduce the (still empty) rows that were removed in presolve. The
assumption is that an empty row cannot be tight, hence the logical is basic
and the dual is zero.
*/
for (int k = 0 ; k < nactions ; k++) {
const action *e = &actions[k] ;
int i = e->row ;
rlo[i] = e->rlo ;
rup[i] = e->rup ;
acts[i] = 0.0 ;
if (rowstat)
prob->setRowStatus(i,CoinPrePostsolveMatrix::basic) ;
rowduals[i] = 0.0 ;
# if PRESOLVE_DEBUG > 0
rdone[i] = DROP_ROW;
# if PRESOLVE_DEBUG > 1
std::cout
<< " restoring row " << i << ", LB = " << rlo[i] << ", UB = " << rup[i]
<< std::endl ;
# endif
# endif
}
prob->nrows_ += nactions ;
assert(prob->nrows_ == prob->nrows0_) ;
# if PRESOLVE_CONSISTENCY > 0 || PRESOLVE_DEBUG > 0
presolve_check_threads(prob) ;
presolve_check_sol(prob,2,2,2) ;
presolve_check_nbasic(prob) ;
# if PRESOLVE_DEBUG > 0
std::cout << "Leaving drop_empty_rows_action::postsolve, system "
<< prob->nrows_ << "x" << prob->ncols_ << "." << std::endl ;
# endif
# endif
}