sec_comp.c

生成直角Steiner树的程序包

/***********************************************************************

	File:	sec_comp.c
	Rev:	b-1
	Date:	02/28/2001

	Copyright (c) 1996, 2001 by David M. Warme

************************************************************************

	Routines to decompose SEC separation problems into smaller
	and simpler sub-problems.

************************************************************************

	Modification Log:

	a-1:	10/05/96	warme
		: Created.
	b-1:	02/28/2001	warme
		: Changes for 3.1 release.
		: Rename many fields in "struct comp".
		: Fix incorrect BCC routine.
		: Use consistent terminology: prefer "vertex" to
		:  "terminal" and "edge" to "full set".
		: Use rverts instead of tmasks.
		: Completely re-write of merge_chains routine.
		: Fix incomplete reduction bug.

************************************************************************/

#include "bb.h"
#include "constrnt.h"
#include "dsuf.h"
#include "genps.h"
#include "sec_comp.h"
#include "sec_heur.h"
#include "steiner.h"


/*
 * Global Routines
 */

struct constraint *	check_component_subtour (bitmap_t *,
						 struct comp *,
						 struct constraint *,
						 double *,
						 bitmap_t *,
						 struct bbinfo *);
struct comp *	delete_vertex_from_component (int, struct comp *);
struct comp *	find_congested_components (double *,
					   bitmap_t *,
					   bitmap_t *,
					   bool,
					   struct cinfo *);
int		find_least_congested_vertex (bitmap_t *, struct comp *);
void		free_congested_component (struct comp *);


/*
 * External References
 */

	/* None */

/*
 * Local Types
 */

struct bc {
	struct comp *	comp;		/* component being split */
	struct comp *	list;		/* output list of BCC's */
	int *		dfs;		/* DFS number of each vertex */
	int *		low;		/* lowest DFS num in component */
	int *		parent;		/* parents of vertices in DFS tree */
	int *		stack;		/* base-address of edge stack */
	int *		sp;		/* current stack pointer */
	int		counter;	/* DFS number generator */
	bitmap_t *	vert_mask;	/* scratch buffer for new components */
	bitmap_t *	edge_mask;	/* scratch buffer for new components */
	bitmap_t *	edges_seen;	/* edges already pushed */
	int		nvmasks;	/* size of vert_mask */
	int		nemasks;	/* size of edge_mask */
	int		ncomps;		/* number of BCC's split off */
	int		max_stack;	/* size of stack */
};


/*
 * Local Routines
 */

static struct constraint * add_component_subtour (struct comp *,
						  double *,
						  bitmap_t *,
						  struct cinfo *,
						  struct constraint *,
						  bitmap_t *);
static void		bcc (struct bc *, int);
static void		component_verts_to_real_verts (bitmap_t *,
						       struct comp *,
						       bitmap_t *,
						       int);
static struct comp *	copy_subcomponent (struct comp *,
					   bitmap_t *,
					   bitmap_t *);
static void		create_masks (struct comp *);
static struct comp *	find_first_component (double *,
					      bitmap_t *,
					      bitmap_t *,
					      struct cinfo *);
static void		merge_chains (struct comp *);
static void		reduce_component_in_place (struct comp *);
static struct comp *	simplify_one_component (struct comp *);
static struct comp *	split_biconnected_components (struct comp *);
static void		split_connected_components (struct comp *);
static void		strip_uncongested_vertices (struct comp *);
static int		vacuous_component (struct comp *);

/*
 * Compute the smallest possible subsets of the problem on which to
 * run the SEC separation routines.  We use the following principles
 * to prune down the vertices and edges examined:
 *
 *	- Let S be a set of vertices that define a violated SEC.  Let
 *	  t be a vertex in S with delta(t) <= 1.  Then there is a
 *	  subset S' of S that does not contain t, yet still defines a
 *	  violated SEC.  Therefore, all such vertices t can be
 *	  iteratively removed from consideration.  (Such vertices are
 *	  NOT "congested".)
 *
 *	- Let the hypergraph of vertices and edges be partitioned
 *	  into its connected components C1 = (V1,F1), C2 = (V2,F2), ...,
 *	  Ck = (Vk,Fk).  Let S be any violated SEC.  Then there is at
 *	  least one i in 1..k for which Si (a subset of
 *	  [Vi \intersect S]) is a violated SEC.
 *
 *	- [Same, but with C1, C2, ..., Ck being biconnected components.]
 *
 *	- Let F1 be an edge containing exactly two congested vertices
 *	  s and t.  Let S be an SEC violation containing s and t.  Define
 *	  a new hypergraph C' = (V - {s,t} + {u}, F') where u is a new
 *	  vertex not in V, and F' is derived from F by replacing every
 *	  reference to s or t with u.  [Resulting edges with fewer than
 *	  2 distinct vertices are deleted from F'.]  Let S be a
 *	  violation in C.  Then S' = S - {s,t} + ({u} if s in S or t in S)
 *	  defines a violation in C'.
 *
 * These rules may be applied recursively to arrive at a final set of
 * congested components that are each usually quit small.
 */

	struct comp *
find_congested_components (

double *		x,		/* IN - current LP solution */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
bool			print_flag,	/* IN - TRUE ==> print info */
struct cinfo *		cip		/* IN - compatibility info */
)
{
int			i;
int			j;
int			k;
int			n;
struct comp *		comp;
struct comp *		p;
bitmap_t *		bp1;

	/* Find the first component by iteratively applying the	*/
	/* delta(t) <= 1 rule...				*/

	comp = find_first_component (x, vert_mask, edge_mask, cip);

	if (comp EQ NULL) return (NULL);

	if (print_flag) {
		tracef (" %% initially %d congested vertices:\n",
			comp -> num_verts);
#if 0
		k = 0;
		for (i = 0; i < comp -> num_verts; i++) {
			int * vp1 = comp -> rverts [i];
			int * vp2 = comp -> rverts [i + 1];
			while (vp1 < vp2) {
				j = *vp1++;
				if (k EQ 0) {
					tracef (" %%	");
				}
				tracef (" %3d", j);
				++k;
				if (k >= 16) {
					tracef ("\n");
					k = 0;
				}
			}
		}
		if (k > 0) {
			tracef ("\n");
		}
#endif
	}

	comp = simplify_one_component (comp);

	if (print_flag) {

		n = 0;
		for (p = comp; p NE NULL; p = p -> next) {
			++n;
		}
		tracef (" %% find_congested_components found %d components:\n",
			n);
		n = 0;
		for (p = comp; p NE NULL; p = p -> next) {
			tracef (" %%\tcomponent %d:\t%d verts,\t%d edges\n",
				n++, p -> num_verts, p -> num_edges);
#if 0
			tracef (" %%\t   ");
			for (i = 0; i < p -> num_verts; i++) {
				int * vp1 = p -> rverts [i];
				int * vp2 = p -> rverts [i + 1];
				int k = (vp2 - vp1);
				if (k NE 1) {
					tracef (" {");
				}
				while (vp1 < vp2) {
					j = *vp1++;
					tracef (" %d", j);
				}
				if (k NE 1) {
					tracef ("}");
				}
			}
			tracef ("\n");
#endif
		}
	}

	return (comp);
}

/*
 * This routine finds the initial component by iteratively applying
 * the delta(t) <= 1 rule.  This normally shrinks the problem down
 * substantially.  We then convert the problem into "struct comp"
 * form -- after which we no longer need to refer to the cinfo stuff.
 */

	static
	struct comp *
find_first_component (

double *		x,		/* IN - current LP solution */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip		/* IN - compatibility info */
)
{
int			i;
int			j;
int			k;
int			e;
int			t;
int			nedges;
int			nverts;
int			nmasks;
int			kmasks;
int			ecount;
int			vcount;
int *			vp1;
int *			vp2;
int *			vp3;
int *			ep1;
int *			ep2;
int *			ep3;
double			sum;
bool			changed;
int *			vleft;
int *			sp;
int *			new_vnum;
int *			new_enum;
int *			old_enum;
int *			ip1;
struct comp *		newp;
bitmap_t *		bp1;
bitmap_t *		bp2;
bitmap_t *		bp3;
bitmap_t *		bp4;
bitmap_t *		verts_stacked;
bitmap_t *		cvmask;
int *			stack;
double *		B;

	nedges = cip -> num_edges;
	nverts = cip -> num_verts;
	nmasks = cip -> num_edge_masks;
	kmasks = cip -> num_vert_masks;

	verts_stacked = NEWA (2 * kmasks, bitmap_t);
	cvmask	      = verts_stacked + kmasks;

	for (i = 0; i < kmasks; i++) {
		verts_stacked [i] = 0;
	}

	/* Compute "vleft" for each valid hyperedge having non-zero	*/
	/* weight.  vleft is 1 less than the number of valid vertices	*/
	/* in the hyperedge.  We decrement vleft every time we delete a	*/
	/* vertex from the hyperedge.  When this count goes to zero, we	*/
	/* can delete the hyperedge.  Also count the valid hyperedges	*/
	/* with non-zero weight.					*/

	vleft = NEWA (nedges, int);
	ecount = 0;
	for (i = 0; i < nedges; i++) {
		vleft [i] = 0;
		if (BITON (edge_mask, i) AND (x [i] > FUZZ)) {
			k = 0;
			vp1 = cip -> edge [i];
			vp2 = cip -> edge [i + 1];
			while (vp1 < vp2) {
				j = *vp1++;
				if (BITON (vert_mask, j)) {
					++k;
				}
			}
			if (k >= 2) {
				vleft [i] = k - 1;
				++ecount;
			}
		}
	}

	/* Compute the congestion level Bi of each vertex i. */
	stack = NEWA (nverts, int);
	B = NEWA (nverts, double);
	sp = &stack [0];
	vcount = 0;
	for (i = 0; i < nverts; i++) {
		sum = 0.0;
		if (BITON (vert_mask, i)) {
			++vcount;
			ep1 = cip -> term_trees [i];
			ep2 = cip -> term_trees [i + 1];
			while (ep1 < ep2) {
				e = *ep1++;
				if (vleft [e] > 0) {
					sum += x [e];
				}
			}
		}
		B [i] = sum;
	}

	/* Add every vertex with weight <= 1 to a list of vertices	*/
	/* to be deleted...						*/
	for (i = 0; i < nverts; i++) {
		if (NOT BITON (vert_mask, i)) {
			/* Pretend this vertex has already been	*/
			/* stacked and deleted (already has weight 0).	*/
			SETBIT (verts_stacked, i);
		}
		else if (B [i] <= (1.0 + FUZZ)) {
			/* Schedule this vertex for deletion! */
			*sp++ = i;
			SETBIT (verts_stacked, i);
		}
	}

	/* Iteratively pop and delete vertices until stack is empty.	*/
	while (sp > stack) {
		t = *--sp;

		/* Prune vertex "t" from the remaining structure... */
		--vcount;
		B [t] = 0.0;
		/* Drop one vertex from every remaining edge that	*/
		/* contains vertex "t"...				*/
		ep1 = cip -> term_trees [t];
		ep2 = cip -> term_trees [t + 1];
		while (ep1 < ep2) {
			e = *ep1++;
			if (vleft [e] <= 0) continue;
			--(vleft [e]);
			if (vleft [e] > 0) continue;
			/* Count for this hyperedge has now gone to	*/
			/* zero, so we can delete it.  Find the edge's	*/
			/* one remaining vertex j and subtract edge e's	*/
			/* weight from Bj.				*/
			--ecount;
			vp1 = cip -> edge [e];
			vp2 = cip -> edge [e + 1];
			for (;;) {
				if (vp1 >= vp2) {
					fatal ("find_first_component: Bug 1.");
				}
				j = *vp1++;
				if (B [j] > FUZZ) break;
			}
			B [j] -= x [e];

			if ((B [j] <= 1.0 + FUZZ) AND
			    (NOT BITON (verts_stacked, j))) {
				/* Schedule this vertex for deletion... */
				*sp++ = j;
				SETBIT (verts_stacked, j);
			}
		}
	}

	/* Construct a mask of the vertices that are left (the	*/
	/* congested vertices).					*/
	bp1 = cvmask;
	bp2 = bp1 + kmasks;
	bp3 = vert_mask;
	bp4 = verts_stacked;
	while (bp1 < bp2) {
		*bp1++ = (*bp3++ & ~(*bp4++));
	}

	if ((vcount EQ 0) OR (ecount EQ 0)) {
		/* Nothing congested left... no components! */
		free ((char *) B);
		free ((char *) stack);
		free ((char *) vleft);
		free ((char *) verts_stacked);
		return (NULL);
	}

	/* Now construct the component for what's left... */

	newp = NEW (struct comp);

	newp -> next		= NULL;
	newp -> num_verts	= vcount;
	newp -> num_edges	= ecount;
	newp -> flags		= CFLG_CONG;
	newp -> x		= NEWA (ecount, double);
	newp -> everts		= NEWA (ecount + 1, int *);
	newp -> vedges		= NEWA (vcount + 1, int *);
	newp -> tviol		= NEWA (vcount, double);
	newp -> rverts		= NEWA (vcount + 1, int *);
	newp -> vert_mask	= NULL;
	newp -> edge_mask	= NULL;
	newp -> cp		= NULL;

	new_vnum = NEWA (nverts, int);
	for (i = 0; i < nverts; i++) {
		new_vnum [i] = -1;
	}
	ip1 = NEWA (vcount, int);
	k = 0;
	for (i = 0; i < nverts; i++) {
		if (NOT BITON (cvmask, i)) continue;
		/* This vertex was retained... */
		new_vnum [i] = k;
		newp -> tviol [k] = 0.0;
		newp -> rverts [k] = ip1;
		*ip1++ = i;
		++k;
	}
	newp -> rverts [k] = ip1;
	if (k NE vcount) {
		fatal ("find_first_component: Bug 2.");
	}

	old_enum = NEWA (nedges, int);
	new_enum = NEWA (nedges, int);
	for (i = 0; i < nedges; i++) {
		new_enum [i] = -1;
	}

	/* Construct hyperedge info, including edge -> vertices map. */
	j = 0;
	k = 0;
	for (i = 0; i < nedges; i++) {
		if (vleft [i] < 1) continue;
		new_enum [i] = j;
		old_enum [j] = i;
		newp -> x [j]	= x [i];
		k += (1 + vleft [i]);
		++j;
	}
	if (j NE ecount) {
		fatal ("find_first_component: Bug 3.");
	}
	vp3 = NEWA (k, int);
	ep3 = NEWA (k, int);
	for (j = 0; j < ecount; j++) {
		newp -> everts [j] = vp3;
		i = old_enum [j];
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			t = *vp1++;
			if (NOT BITON (cvmask, t)) continue;
			t = new_vnum [t];
			if (t < 0) {
				fatal ("find_first_component: Bug 4.");
			}
			*vp3++ = t;
		}
	}
	newp -> everts [ecount] = vp3;
	if (vp3 NE newp -> everts [0] + k) {
		fatal ("find_first_component: Bug 5.");
	}
	for (i = 0; i < nverts; i++) {
		j = new_vnum [i];
		if (j < 0) continue;
		newp -> vedges [j] = ep3;
		ep1 = cip -> term_trees [i];
		ep2 = cip -> term_trees [i + 1];
		while (ep1 < ep2) {
			e = new_enum [*ep1++];
			if (e >= 0) {
				*ep3++ = e;
			}
		}
	}
	newp -> vedges [vcount] = ep3;
	if (ep3 NE newp -> vedges [0] + k) {
		fatal ("find_first_component: Bug 6.");
	}

	free ((char *) new_enum);
	free ((char *) old_enum);
	free ((char *) new_vnum);
	free ((char *) B);
	free ((char *) stack);
	free ((char *) vleft);
	free ((char *) verts_stacked);

	return (newp);
}

/*
 * This routine simplifies a SINGLE COMPONENT, producing 0 or more
 * simplified components that are concatenated onto the front of the
 * successors that may follow the given component in the linked list.
 */

	static
	struct comp *
simplify_one_component (

struct comp *		comp		/* IN - component to simplify */
)
{
struct comp *		rest;
struct comp *		p;
struct comp *		temp;
struct comp *		done;
struct comp **		done_hookp;