grbranching.c

#include	"grhdr.h"

/******************************************************************************
	Chu-Liu-Edmonds algorithm to compute a maximum branching for a
	directed graph with edge weights.  The algorithm alters the edge
	lists to contain just the branching edges.

	Written by Kiem-Phong Vo (02/01/2006)
******************************************************************************/

typedef struct _brcycl_s
{	Gredge_t*	cycl;	/* the cycle of strong components	*/
	Gredge_t*	emin;	/* the minimum weight edge of cycle	*/
	Gredge_t*	entr;	/* the edge that came in to the cycle	*/
} Brcycl_t;

typedef struct _brnode_s
{	Grdata_t	data;
	ssize_t		mark;	/* mark 1 if already searched		*/
} Brnode_t;

typedef struct _bredge_s
{	Grdata_t	data;
	ssize_t		wght;	/* the original edge weight		*/
	ssize_t		wadj;	/* adjusted weight during construction	*/
	Grnode_t*	root;	/* strong component before collapsing	*/
	Gredge_t*	edge;	/* branching edge that this represents	*/
} Bredge_t;

#define BRNODE(n)	((Brnode_t*)grdtnode((n), Grbranching) )
#define BREDGE(e)	((Bredge_t*)grdtedge((e), Grbranching) )

/* compute the representative node in the shadowed union structure */
#define GRLINK(n)	((n)->link == (n) ? (n) : grlink(n))
static Grnode_t* grlink(Grnode_t* n)
{	while(n->link != n)
		n = n->link;
	return n;
}

#ifdef DEBUG
static int	Fd = 2;
static int predge(Gredge_t* ed)
{
	if(!ed)
	{	PRINT(Fd,"Null edge\n",0);
		return 0;
	}
	PRINT(Fd,"%d", (int)ed->tail->label); PRINT(Fd,"(%d) -> ", (int)grfind(ed->tail)->label );
	PRINT(Fd,"%d", (int)ed->head->label); PRINT(Fd,"(%d), ", (int)grfind(ed->head)->label );
	PRINT(Fd,"root=%d, ", BREDGE(ed)->root ? (int)BREDGE(ed)->root->label : -1 );
	PRINT(Fd,"wght=%d, ", (int)BREDGE(ed)->wght); PRINT(Fd,"wadj=%d\n", (int)BREDGE(ed)->wadj);
	return 0;
}
static int prlink(Gredge_t* e)
{	for(; e; e = e->link)
		{ PRINT(Fd, "\t", 0); predge(e); }
	PRINT(Fd,"\n",0);
}
static int prinext(Gredge_t* e)
{	for(; e; e = e->inext)
		{ PRINT(Fd, "\t", 0); predge(e); }
	PRINT(Fd,"\n",0);
}
static int pronext(Gredge_t* e)
{	for(; e; e = e->onext)
		{ PRINT(Fd, "\t", 0); predge(e); }
	PRINT(Fd,"\n",0);
}
static int prnode(Grnode_t* nd)
{	
	if(!nd)
	{	PRINT(Fd,"Null node\n",0);
		return 0;
	}
	PRINT(Fd,"node = %d", (int)nd->label); PRINT(Fd,"(%d), ", grfind(nd)->label);
	PRINT(Fd,"link = %d, ", nd->link ? (int)nd->link->label : -1);
	PRINT(Fd,"mark = %d\n", (int)BRNODE(nd)->mark);
	return 0;
}
static int prid(Graph_t* gr, ssize_t id)
{
	return prnode(grnode(gr, (void*)id, 0));
}
static int prcycl(Brcycl_t* cl)
{
	Gredge_t *e; 
	PRINT(Fd,"Emin: ",0); predge(cl->emin);
	PRINT(Fd,"Entr: ",0); predge(cl->entr);
	prlink(cl->cycl);
	return 0;
}
static int prgraph(Graph_t* gr)
{	Grnode_t	*nd;
	Gredge_t	*ed;
	for(nd = (Grnode_t*)dtflatten(gr->nodes); nd; nd = (Grnode_t*)dtlink(gr->nodes,nd) )
	{	if(nd != grfind(nd))
			continue;
		prnode(nd);
		for(ed = nd->iedge; ed; ed = ed->inext)
			{ PRINT(Fd,"\t",0); predge(ed); }
	}
}
#endif

ssize_t grbranching(Graph_t* gr)
{
	Grnode_t	*n, *nc;
	Gredge_t	*e, *ec, *ep, *en, *path, *emin;
	Brcycl_t	*clist, *cl, *endcl;
	ssize_t		w;

	for(w = 0, n = (Grnode_t*)dtflatten(gr->nodes); n; w += 1, n = (Grnode_t*)dtlink(gr->nodes,n) )
	{	n->link = n->fold = n; /* union structures: link kept as-is, fold does path-compression */
		n->oedge = NIL(Gredge_t*); /* wipe the out-edges */
		BRNODE(n)->mark = 0;

		/* compute heaviest weight ec and move it to front */
		for(ep = en = NIL(Gredge_t*), ec = e = n->iedge; e; en = e, e = e->inext )
		{	e->link = NIL(Gredge_t*);
			BREDGE(e)->edge = e; /* at start, representing self */
			if((BREDGE(e)->wadj = BREDGE(e)->wght) > BREDGE(ec)->wadj )
				{ ec = e; ep = en; }
		}
		if(ep)
			{ ep->inext = ec->inext; ec->inext = n->iedge; n->iedge = ec; }
	}
	if(w == 0)
		return 0;

	/* space to keep cycle structures */
	if(!(clist = cl = (Brcycl_t*)calloc(w+1,sizeof(Brcycl_t))) )
		return -1;
	endcl = cl+w+1;

	/* search and collapse cycles */
	for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes,n) )
	{	nc = grfind(n);
		if(BRNODE(nc)->mark) /* already searched */
			continue;

		for(path = NIL(Gredge_t*); nc; ) /* depth-first search */
		{	if(!BRNODE(nc)->mark) /* not searched yet */
			{	BRNODE(nc)->mark = 1;

				if(!(ec = nc->iedge) ) /* this path cannot be a cycle */
					break;
				else /* continue to build the path */
				{	ec->link = path; path = ec;
					nc = grfind(ec->tail);
					continue; /* dfs recursion */
				}
			}

			/* potential cycle, check that out and also compute min edge */
			for(emin = NIL(Gredge_t*), ec = path; ec; ec = ec->link)
			{	if(!emin || BREDGE(ec)->wadj < BREDGE(emin)->wadj )
					emin = ec;
				if(grfind(ec->head) == nc) /* end of cycle */
					break;
			}
			if(!ec) /* run up against an old path */
				break;

			if(cl >= endcl) /* hmm, this should not happen */
			{	/**/ASSERT(cl < endcl);
				free(clist);
				return -1;
			}

			cl->cycl = path;
			path = ec->link; /* will restart search from here */
			ec->link = NIL(Gredge_t*);

			/* make list of incoming edges and adjust their weights */
			en = NIL(Gredge_t*);
			for(ec = cl->cycl; ec; ec = ec->link)
			{	BREDGE(ec)->root = grfind(ec->head); /* save node union structure */

				w = BREDGE(ec)->wadj - BREDGE(emin)->wadj;
				for(ep = NIL(Gredge_t*), e = grfind(ec->head)->iedge; e; e = e->inext)
				{	BREDGE(e)->wadj -= w;
					if(!e->inext) /* last of list */
						ep = e;
				}
				ep->inext = en; en = grfind(ec->head)->iedge; /* catenate lists */
				grfind(ec->head)->iedge = NIL(Gredge_t*);
			}

			/* collapsing cycle onto nc */
			for(ec = cl->cycl; ec; ec = ec->link)
			{	if(grfind(ec->head) == nc)
					continue;
				grfind(ec->head)->link = nc; /* union history kept as-is */
				grfind(ec->head)->fold = nc; /* union with path-compression */
			}
			nc->fold = nc->link = nc;

			/* make new edge list, keep heaviest edge in front */
			for(e = en; e; e = en)
			{	en = e->inext;
				if(BREDGE(e)->wadj <= 0 || grfind(e->head) == grfind(e->tail) )
					continue;
				if(!nc->iedge || BREDGE(nc->iedge)->wadj <= BREDGE(e)->wadj )
					{ e->inext = nc->iedge; nc->iedge = e; }
				else	{ e->inext = nc->iedge->inext; nc->iedge->inext = e; }
			}

			cl->emin = emin;
			cl->entr = nc->iedge;
			cl += 1;
			BRNODE(nc)->mark = 0; /* force nc to be searched again */
		}
	}

	/* move the remaining branching edges to their real nodes */
	path = NIL(Gredge_t*);
	for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes,n) )
	{	if(!(e = n->iedge) )
			continue;
		e->link = path; path = e;
	}
	for(ec = path; ec; ec = ec->link)
		ec->head->iedge = ec; 

	/* unroll collapsed cycles in reverse order to construct branching */
	for(cl -= 1; cl >= clist; --cl )
	{	/* restore the union structure to just before this cycle collapsed */
		for(ec = cl->cycl; ec; ec = ec->link)
			BREDGE(ec)->root->link = BREDGE(ec)->root;

		if((en = cl->entr) )
			en = BREDGE(en)->edge;
		for(ec = cl->cycl; ec; ec = ec->link)
		{	if(en && GRLINK(ec->head) == GRLINK(en->head))
				BREDGE(ec)->edge = en;
			else	ec->head->iedge = ec;
		}
		if(!en) /* isolated cycle, remove minimum edge */
		{	cl->emin->head->iedge = NIL(Gredge_t*);
			if(BREDGE(cl->emin)->edge == cl->emin)
				BREDGE(cl->emin)->edge = NIL(Gredge_t*);
		}
	}

	w = 0; /* construct the external branching representation */
	for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes,n) )
	{	if(!(e = n->iedge) )
			continue;
		e->inext = NIL(Gredge_t*);
		e->onext = e->tail->oedge; e->tail->oedge = e;
		w += BREDGE(e)->wght;
	}

	free(clist);

	return w;
}


ssize_t grbrweight(Gredge_t* e, ssize_t w)
{
	if(w > 0)
		BREDGE(e)->wght = w;
	return BREDGE(e)->wght;
}


static Grdata_t* brdata(Gralgo_t* algo, int type, Grdata_t* data)
{
	Brnode_t*	brn;
	Bredge_t*	bre;

	if(algo != Grbranching)
		return NIL(Grdata_t*);

	switch(type)
	{ default:
		return NIL(Grdata_t*);

	  case GR_NODE|GR_CLOSING:
	  case GR_EDGE|GR_CLOSING:
		if(data)
			free(data);
		return NIL(Grdata_t*);

	  case GR_NODE|GR_OPENING:
		if(!(brn = (Brnode_t*)calloc(1, sizeof(Brnode_t))) )
			return NIL(Grdata_t*);
		return (Grdata_t*)brn;

	  case GR_EDGE|GR_OPENING:
		if(!(bre = (Bredge_t*)calloc(1, sizeof(Bredge_t))) )
			return NIL(Grdata_t*);
		return (Grdata_t*)bre;
	}
}

static Gralgo_t	_Grbranching = { brdata, 0 };
Gralgo_t*	Grbranching = &_Grbranching;

#ifdef KPVTEST
#include	<stdio.h>
main(int argc, char** argv)
{
	Graph_t		*gr;
	Grnode_t	*n, *hn, *tn;
	Gredge_t	*e;
	int		h, t, w;
	char		buf[1024];
	int		type = 0;

	if(argc > 1 && strcmp(argv[1],"-l") == 0)
		type = -1;
	if(argc > 1 && strcmp(argv[1],"-r") == 0)
		type = 1;

	gr = gropen(0);
	while(fgets(buf,sizeof(buf),stdin))
	{	if(buf[0] == '#')
			continue;
		if(sscanf(buf,"%d,%d,%d", &t, &h, &w) != 3)
			continue;
		tn = grnode(gr, t, 1);
		hn = grnode(gr, h, 1);
		if(type > 0 && t >= h) /* only use edges with t < h */
			continue; 
		if(type < 0 && t <= h) /* only use edges with t > h */
			continue; 
		e  = gredge(gr, tn, hn, 0, 1);
		grbrweight(e, w);
	}

	printf("\nTotal weight=%d\n", grbranching(gr));
	for(n = (Grnode_t*)dtflatten(gr->nodes); n; n = (Grnode_t*)dtlink(gr->nodes, n) )
		if((e = n->iedge) )
			fprintf(stderr, "%d -> %d [%d]\n", e->tail->label, e->head->label, grbrweight(e, 0) );
}
#endif