/* maxFlow.c */
#include "../BPG.h"
/*--------------------------------------------------------------------*/
/*
----------------------------------------
a Cell structure represent an edge (x,y)
rxy -- residual for edge
fxy -- flow for edge
xnext -- next Cell in list for x
ynext -- next Cell in list for y
----------------------------------------
*/
typedef struct _Cell Cell ;
struct _Cell {
int x ;
int y ;
int rxy ;
int fxy ;
Cell *xnext ;
Cell *ynext ;
} ;
/*--------------------------------------------------------------------*/
/*
-----------------
static prototypes
-----------------
*/
static void
setupCells(BPG *bpg, Cell ***pheads, Cell **pcells,
int msglvl, FILE *msgFile) ;
static int
findFlowAugmentingPath(BPG *bpg, Cell *heads[], int xexp, int exp[],
int par[], int labels[], int list[], int mark[],
int tag, int msglvl, FILE *msgFile) ;
static void
augmentPath(BPG *bpg, int xexp, int yexp, int nvexp[], Cell *heads[],
int labels[], int par[], int msglvl, FILE *msgFile) ;
static void
setDMflags(BPG *bpg, Cell *heads[], int nvexp[], int list[], int mark[],
int dmflags[], int stats[], int msglvl, FILE *msgFile) ;
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------
find the DM decomposition of a weighted bipartite graph
using a simple max flow algorithm
bpg -- BPG bipartite graph object
dmflags -- flags vector
/ 0 if x in X_R
dmflags[x] = | 1 if x in X_I
\ 2 if x in X_E
/ 0 if y in Y_R
dmflags[y] = | 1 if y in Y_I
\ 2 if y in Y_E
stats -- statistics vector
stats[0] -- weight of X_I
stats[1] -- weight of X_E
stats[2] -- weight of X_R
stats[3] -- weight of Y_I
stats[4] -- weight of Y_E
stats[5] -- weight of Y_R
created -- 96mar08, cca
-------------------------------------------------------
*/
void
BPG_DMviaMaxFlow (
BPG *bpg,
int dmflags[],
int stats[],
int msglvl,
FILE *msgFile
) {
int ierr, nX, nY, tag, x, xexp, y, yexp ;
int *labels, *list, *mark, *nvexp, *par, *vwghts ;
Cell *cell, *cells ;
Cell **heads ;
/*
---------------
check the input
---------------
*/
if ( bpg == NULL || dmflags == NULL || stats == NULL
|| (msglvl > 0 && msgFile == NULL) ) {
fprintf(stderr, "\n fatal error in BPG_DMviaMaxFlow(%p,%p,%p,%d,%p)"
"\n bad input\n", bpg, dmflags, stats, msglvl, msgFile) ;
exit(-1) ;
}
nX = bpg->nX ;
nY = bpg->nY ;
/*
---------------------------------------------
set up the Cell structures and working arrays
---------------------------------------------
*/
setupCells(bpg, &heads, &cells, msglvl, msgFile) ;
nvexp = IVinit(nX + nY, 1) ;
if ( (vwghts = bpg->graph->vwghts) != NULL ) {
IVcopy(nX + nY, nvexp, vwghts) ;
}
list = IVinit(nX + nY, -1) ;
mark = IVinit(nX + nY, -1) ;
par = IVinit(nX + nY, -1) ;
labels = IVinit(nX + nY, -1) ;
/*
------------------------
loop over the x vertices
------------------------
*/
tag = 0 ;
for ( xexp = 0 ; xexp < nX ; xexp++ ) {
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n checking out x node %d, exposure %d",
xexp, nvexp[xexp]) ;
}
while ( nvexp[xexp] > 0 ) {
yexp = findFlowAugmentingPath(bpg, heads, xexp, nvexp,
par, labels, list, mark, tag,
msglvl, msgFile) ;
tag++ ;
if ( yexp == -1 ) {
/*
-------------------------------------------------------
no flow augmenting path has been found, node is exposed
-------------------------------------------------------
*/
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n node %d has exposure %d",
xexp, nvexp[xexp]) ;
}
break ;
} else {
/*
----------------------------------------
flow augmenting path found, augment path
----------------------------------------
*/
augmentPath(bpg, xexp, yexp, nvexp, heads, labels, par,
msglvl, msgFile) ;
}
}
}
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n after maximum flow found") ;
for ( x = 0 ; x < nX ; x++ ) {
fprintf(msgFile, "\n X node %d : ", x) ;
for ( cell = heads[x] ; cell != NULL ; cell = cell->xnext ) {
fprintf(msgFile, "\n (%3d,%3d,%3d)",
cell->y, cell->fxy, cell->rxy) ;
}
}
for ( y = nX ; y < nX + nY ; y++ ) {
fprintf(msgFile, "\n Y node %d : ", y) ;
for ( cell = heads[y] ; cell != NULL ; cell = cell->ynext ) {
fprintf(msgFile, "\n (%3d,%3d,%3d)",
cell->x, cell->fxy, cell->rxy) ;
}
}
}
/*
----------------
set the DM flags
----------------
*/
setDMflags(bpg, heads, nvexp, list, mark,
dmflags, stats, msglvl, msgFile) ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n X dmflags: ") ;
IVfp80(msgFile, nX, dmflags, 12, &ierr) ;
fprintf(msgFile, "\n Y dmflags: ") ;
IVfp80(msgFile, nY, dmflags + nX, 12, &ierr) ;
fflush(msgFile) ;
}
/*
------------------------
free the working storage
------------------------
*/
IVfree(nvexp) ;
IVfree(list) ;
IVfree(mark) ;
IVfree(par) ;
IVfree(labels) ;
FREE(heads) ;
FREE(cells) ;
return ; }
/*--------------------------------------------------------------------*/
/*
------------
set up cells
------------
*/
static void
setupCells (
BPG *bpg,
Cell ***pheads,
Cell **pcells,
int msglvl,
FILE *msgFile
) {
Cell *cell, *cells ;
Cell **heads ;
Graph *graph ;
int ii, ncell, nX, nY, x, xsize, y ;
int *vwghts, *xadj ;
/*
-------------------------
count the number of cells
-------------------------
*/
nX = bpg->nX ;
nY = bpg->nY ;
graph = bpg->graph ;
ncell = 0 ;
for ( x = 0 ; x < nX ; x++ ) {
Graph_adjAndSize(graph, x, &xsize, &xadj) ;
ncell += xsize ;
}
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n %d cells required", ncell) ;
fflush(msgFile) ;
}
/*
------------------------------------------------
allocate the storage and set pointers for return
------------------------------------------------
*/
ALLOCATE(cells, struct _Cell, ncell) ;
ALLOCATE(heads, struct _Cell *, nX + nY) ;
*pcells = cells ;
*pheads = heads ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n cells = %p", cells) ;
fflush(msgFile) ;
}
for ( ii = 0 ; ii < nX + nY ; ii++ ) {
heads[ii] = NULL ;
}
/*
----------------------------------
set the vertex fields of the cells
----------------------------------
*/
for ( x = 0, cell = cells ; x < nX ; x++ ) {
Graph_adjAndSize(graph, x, &xsize, &xadj) ;
for ( ii = 0 ; ii < xsize ; ii++, cell++ ) {
y = xadj[ii] ;
cell->x = x ;
cell->y = y ;
}
}
/*
--------------------------------
set the link fields of the cells
--------------------------------
*/
for ( ii = 0, cell = cells ; ii < ncell ; ii++, cell++ ) {
x = cell->x ;
cell->xnext = heads[x] ;
heads[x] = cell ;
y = cell->y ;
cell->ynext = heads[y] ;
heads[y] = cell ;
}
/*
---------------------------------------------
set the residual and flow fields of the cells
---------------------------------------------
*/
if ( (vwghts = graph->vwghts) == NULL ) {
for ( ii = 0, cell = cells ; ii < ncell ; ii++, cell++ ) {
cell->rxy = 1 ;
cell->fxy = 0 ;
}
} else {
for ( ii = 0, cell = cells ; ii < ncell ; ii++, cell++ ) {
x = cell->x ;
y = cell->y ;
cell->rxy = (vwghts[x] < vwghts[y]) ? vwghts[x] : vwghts[y] ;
cell->fxy = 0 ;
}
}
if ( msglvl > 2 && msgFile != NULL ) {
for ( ii = 0, cell = cells ; ii < ncell ; ii++, cell++ ) {
fprintf(msgFile,
"\n cell %d, address %p : (x,y) = (%3d,%3d), r = %d, f = %d",
ii, cell, cell->x, cell->y, cell->rxy, cell->fxy) ;
}
fflush(msgFile) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------
search for a flow augmenting path starting at xexp
created -- 96mar08, cca
--------------------------------------------------
*/
static int
findFlowAugmentingPath (
BPG *bpg,
Cell *heads[],
int xexp,
int exp[],
int par[],
int labels[],
int list[],
int mark[],
int tag,
int msglvl,
FILE *msgFile
) {
Cell *cell ;
Graph *graph ;
int last, now, nX, x, y ;
/*
---------------------------
get dimensions and pointers
---------------------------
*/
nX = bpg->nX ;
graph = bpg->graph ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n FIND AUGMENTING PATH FOR NODE %d", xexp) ;
fflush(msgFile) ;
}
/*
-----------------------------------
mark the exposed node with the tag,
label with its exposure,
and put into the list for traversal
-----------------------------------
*/
mark[xexp] = tag ;
list[0] = xexp ;
labels[xexp] = exp[xexp] ;
par[xexp] = -1 ;
now = last = 0 ;
while ( now <= last ) {
/*
-------------------------------
check out next node on the list
-------------------------------
*/
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n list[%d] = %d", now, list[now]) ;
fflush(msgFile) ;
}
if ( list[now] < nX ) {
x = list[now++] ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", X node %d", x) ;
fflush(msgFile) ;
}
for ( cell = heads[x] ; cell != NULL ; cell = cell->xnext ) {
y = cell->y ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n adjacent Y node %d, mark %d, r %d",
y, mark[y], cell->rxy) ;
fflush(msgFile) ;
}
if ( mark[y] != tag && cell->rxy > 0 ) {
/*
----------------------------------------------------------
y has not yet been visited and edge (x,y) is not saturated
----------------------------------------------------------
*/
mark[y] = tag ;
par[y] = x ;
if ( labels[x] < cell->rxy ) {
labels[y] = labels[x] ;
} else {
labels[y] = cell->rxy ;
}
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", label = %d", labels[y]) ;
fflush(msgFile) ;
}
if ( exp[y] > 0 ) {
/*
----------------------------------
y is exposed, set label and return
----------------------------------
*/
if ( labels[y] > exp[y] ) {
labels[y] = exp[y] ;
}
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile,
"\n Y node %d is exposed, label = %d",
y,labels[y]) ;
fflush(msgFile) ;
}
return(y) ;
}
list[++last] = y ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", adding to list") ;
fflush(msgFile) ;
}
}
}
} else {
y = list[now++] ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", Y node %d", y) ;
fflush(msgFile) ;
}
for ( cell = heads[y] ; cell != NULL ; cell = cell->ynext ) {
x = cell->x ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n adjacent X node %d, mark %d, f %d",
x, mark[x], cell->fxy) ;
fflush(msgFile) ;
}
if ( mark[x] != tag && cell->fxy > 0 ) {
/*
--------------------------------------------------------
x has not yet been visited and there is flow from x to y
--------------------------------------------------------
*/
mark[x] = tag ;
par[x] = y ;
if ( labels[y] < cell->fxy ) {
labels[x] = labels[y] ;
} else {
labels[x] = cell->fxy ;
}
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", labels = %d", labels[x]) ;
fflush(msgFile) ;
}
list[++last] = x ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", adding to list") ;
fflush(msgFile) ;
}
}
}
}
}
return(-1) ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------------------------
given an augmenting path starting with xexp and ending with yexp,
augment the path by the amount in labels[yexp]
96mar08, cca
-----------------------------------------------------------------
*/
static void
augmentPath (
BPG *bpg,
int xexp,
int yexp,
int nvexp[],
Cell *heads[],
int labels[],
int par[],
int msglvl,
FILE *msgFile
) {
Cell *cell ;
int delta, x, y ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n AUGMENT PATH FOR NODE %d", xexp) ;
fflush(msgFile) ;
}
/*
---------------------------
start at the exposed y node
---------------------------
*/
delta = labels[yexp] ;
y = yexp ;
while ( 1 ) {
x = par[y] ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n (x,y) = (%3d, %3d)", x, y) ;
fflush(msgFile) ;
}
for ( cell = heads[y] ; cell != NULL ; cell = cell->ynext ) {
if ( cell->x == x ) {
cell->rxy -= delta ;
cell->fxy += delta ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", r = %d, f = %d", cell->rxy, cell->fxy) ;
fflush(msgFile) ;
}
break ;
}
}
if ( cell == NULL ) {
fprintf(stderr, "\n 1. error, x = %d, y = %d", x, y) ;
exit(-1) ;
}
if ( x == xexp ) {
break ;
}
y = par[x] ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n (x,y) = (%3d, %3d)", x, y) ;
fflush(msgFile) ;
}
for ( cell = heads[x] ; cell != NULL ; cell = cell->xnext ) {
if ( cell->y == y ) {
cell->rxy += delta ;
cell->fxy -= delta ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", r = %d, f = %d", cell->rxy, cell->fxy) ;
fflush(msgFile) ;
}
break ;
}
}
if ( cell == NULL ) {
fprintf(stderr, "\n 2. error, x = %d, y = %d", x, y) ;
exit(-1) ;
}
}
nvexp[xexp] -= delta ;
nvexp[yexp] -= delta ;
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------
set the flags for the DM decomposition
created -- 96mar08, cca
--------------------------------------
*/
static void
setDMflags (
BPG *bpg,
Cell *heads[],
int nvexp[],
int list[],
int mark[],
int dmflags[],
int stats[],
int msglvl,
FILE *msgFile
) {
Cell *cell ;
int ierr, last, now, nX, nY, x, xwght, y, ywght ;
int *vwghts ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n SET DM FLAGS") ;
fflush(msgFile) ;
}
nX = bpg->nX ;
nY = bpg->nY ;
IVfill(nX + nY, dmflags, 0) ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n nvexp : ") ;
IVfp80(msgFile, nX + nY, nvexp, 20, &ierr) ;
fflush(msgFile) ;
}
/*
-----------------------
load exposed nodes in X
-----------------------
*/
IVzero(nX + nY, mark) ;
last = -1 ;
for ( x = 0 ; x < nX ; x++ ) {
if ( nvexp[x] > 0 ) {
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n X node %d is exposed", x) ;
fflush(msgFile) ;
}
mark[x] = 1 ;
list[++last] = x ;
dmflags[x] = 1 ;
}
}
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n %d exposed X nodes :", last + 1) ;
IVfp80(msgFile, 1 + last, list, 20, &ierr) ;
fflush(msgFile) ;
}
/*
---------------------------------------------------------------
drop an alternating level structure from the exposed nodes in X
---------------------------------------------------------------
*/
now = 0 ;
while ( now <= last ) {
if ( list[now] < nX ) {
x = list[now++] ;
dmflags[x] = 1 ;
for ( cell = heads[x] ; cell != NULL ; cell = cell->xnext ) {
y = cell->y ;
/*
if ( mark[y] != 1 && cell->rxy > 0 ) {
*/
if ( mark[y] != 1 ) {
mark[y] = 1 ;
list[++last] = y ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n adding Y node %d", y) ;
fflush(msgFile) ;
}
}
}
} else {
y = list[now++] ;
dmflags[y] = 2 ;
for ( cell = heads[y] ; cell != NULL ; cell = cell->ynext ) {
x = cell->x ;
if ( mark[x] != 1 && cell->fxy > 0 ) {
mark[x] = 1 ;
list[++last] = x ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n adding X node %d", x) ;
fflush(msgFile) ;
}
}
}
}
}
/*
-----------------------
load exposed nodes in Y
-----------------------
*/
IVzero(nX + nY, mark) ;
last = -1 ;
for ( y = nX ; y < nX + nY ; y++ ) {
if ( nvexp[y] > 0 ) {
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n Y node %d is exposed", y) ;
fflush(msgFile) ;
}
mark[y] = 1 ;
list[++last] = y ;
dmflags[y] = 1 ;
}
}
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n %d exposed Y nodes :", last + 1) ;
IVfp80(msgFile, 1 + last, list, 20, &ierr) ;
fflush(msgFile) ;
}
/*
---------------------------------------------------------------
drop an alternating level structure from the exposed nodes in Y
---------------------------------------------------------------
*/
now = 0 ;
while ( now <= last ) {
if ( list[now] < nX ) {
x = list[now++] ;
dmflags[x] = 2 ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n checking out X node %d", x) ;
fflush(msgFile) ;
}
for ( cell = heads[x] ; cell != NULL ; cell = cell->xnext ) {
y = cell->y ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n adjacent Y node %d, mark %d, f %d",
y, mark[y], cell->fxy) ;
fflush(msgFile) ;
}
if ( mark[y] != 1 && cell->fxy > 0 ) {
mark[y] = 1 ;
list[++last] = y ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", adding ") ;
fflush(msgFile) ;
}
}
}
} else {
y = list[now++] ;
dmflags[y] = 1 ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n checking out Y node %d", y) ;
fflush(msgFile) ;
}
for ( cell = heads[y] ; cell != NULL ; cell = cell->ynext ) {
x = cell->x ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, "\n adjacent X node %d, mark %d, r %d",
x, mark[x], cell->rxy) ;
fflush(msgFile) ;
}
/*
if ( mark[x] != 1 && cell->rxy > 0 ) {
*/
if ( mark[x] != 1 ) {
mark[x] = 1 ;
list[++last] = x ;
if ( msglvl > 2 && msgFile != NULL ) {
fprintf(msgFile, ", adding") ;
fflush(msgFile) ;
}
}
}
}
}
/*
--------------------------
fill the statistics vector
--------------------------
*/
IVzero(6, stats) ;
vwghts = bpg->graph->vwghts ;
for ( x = 0 ; x < nX ; x++ ) {
xwght = (vwghts != NULL) ? vwghts[x] : 1 ;
switch ( dmflags[x] ) {
case 0 : stats[2] += xwght ; break ;
case 1 : stats[0] += xwght ; break ;
case 2 : stats[1] += xwght ; break ;
default :
fprintf(stderr, "\n fatal error in BPG_DMviaMaxFlow"
"\n dmflags[%d] = %d\n", x, dmflags[x]) ;
exit(-1) ;
}
}
for ( y = nX ; y < nX + nY ; y++ ) {
ywght = (vwghts != NULL) ? vwghts[y] : 1 ;
switch ( dmflags[y] ) {
case 0 : stats[5] += ywght ; break ;
case 1 : stats[3] += ywght ; break ;
case 2 : stats[4] += ywght ; break ;
default :
fprintf(stderr, "\n fatal error in BPG_DMviaMaxFlow"
"\n dmflags[%d] = %d\n", y, dmflags[y]) ;
exit(-1) ;
}
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1