prunefst.c

生成直角Steiner树的程序包

	static
	void
test_close_terminal (

struct pinfo *		pip,	/* IN	  - pruning data structure  */
int			fst,	/* IN	  - FST index */
int			term,	/* IN	  - Terminal */
struct clt_info**	clip	/* IN/OUT - store close terminal info here */
)
{
int			j;
int			i1;
int			i2;
dist_t			d;
dist_t			d1;
dist_t			d2;
dist_t			pg_longest;
struct point *		pt;
struct point *		p1;
struct point *		p2;
struct point		clp;
struct full_set *	fsp;
struct pset *		fsp_terms;
struct pset *		fsp_steins;
struct cinfo *		cip;
struct clt_info		clt;

	cip	   = pip -> cip;
	pt	   = &(cip -> pts -> a [term]);
	fsp	   = cip -> full_trees [fst];
	fsp_terms  = cip -> full_trees [fst] -> terminals;
	fsp_steins = cip -> full_trees [fst] -> steiners;
	pg_longest = pip -> pg_edges [pip -> num_pg_edges-1].len;

	/* First a rough test to eliminate the terminal */
	if (EDIST(&(fsp_terms -> a[0]), pt) >
	    (fsp_terms -> n - 1) * pg_longest)
		return; /* this terminal is too far away */

	/* Now find the edge that is closest edge to this terminal */
	memset (&clt, 0, sizeof (clt));
	clt.term = term;
	clt.dist = INF_DISTANCE;
	for (j = 0; j < fsp -> nedges; j++) {

		/* Get coordinates and pruning graph indices for endpoints */
		if (fsp -> edges[j].p1 < fsp_terms -> n) {
			p1 = &(fsp_terms  -> a[ fsp -> edges[j].p1 ]);
			i1 = fsp_terms -> a[ fsp -> edges[j].p1 ].pnum;
		}
		else {
			p1 = &(fsp_steins -> a[ fsp -> edges[j].p1 -
						fsp_terms -> n ]);
			i1 = pip -> steiner_index [fst] +
				fsp -> edges[j].p1 - fsp_terms -> n;
		}

		if (fsp -> edges[j].p2 < fsp_terms -> n) {
			p2 = &(fsp_terms  -> a[ fsp -> edges[j].p2 ]);
			i2 = fsp_terms -> a[ fsp -> edges[j].p2 ].pnum;
		}
		else {
			p2 = &(fsp_steins -> a[ fsp -> edges[j].p2 -
						fsp_terms -> n ]);
			i2 = pip -> steiner_index [fst] +
				fsp -> edges[j].p2 - fsp_terms -> n;
		}

		/* Compute closest distance from terminal to edge */
		d = terminal_edge_distance(cip, pt, p1, p2, &clp, &d1, &d2) *
		    (1.0 + pip -> eps * ((double) cip -> edge_size [fst]));
		if (d < clt.dist) {
			clt.dist = d;
			clt.aterm1 = (d1 <= d) ? i1 : pip -> num_pg_verts;
			clt.aterm2 = (d2 <= d) ? i2 : pip -> num_pg_verts;
		}
	}

	/* Is distance smaller than longest edge in pruning graph? */
	if (clt.dist < pg_longest) {
		*((*clip)++) = clt; /* save this terminal */
	}
}

/*
 * Distance from point to edge
 */

	static
	dist_t
terminal_edge_distance (

struct cinfo *		cip,	/* IN  - compatibility info */
struct point *		pt,	/* IN  - point */
struct point *		p1,	/* IN  - first edge point */
struct point *		p2,	/* IN  - second edge point */
struct point *		clp,	/* OUT - closest point */
dist_t *		d1,	/* OUT - distance from clp to p1 */
dist_t *		d2	/* OUT - distance to clp to p2 */
)
{
dist_t			d;
dist_t			l1;
dist_t			l2;
dist_t			l;
struct point *		a;
struct point *		b;
struct point		e;
struct point		ctr;

	d = 0.0;
	if (cip -> metric EQ RECTILINEAR) {

		/* Compute distance to rectangle given by p1 and p2.
		   Also identify (one) closest point. */

		l1 = pt -> x - p1 -> x;
		l2 = pt -> x - p2 -> x;
		clp -> x = pt -> x;
		if ((l1 < 0.0) AND (l2 < 0.0)) {
			clp -> x = MIN(p1 -> x, p2 -> x);
		}
		if ((l1 > 0.0) AND (l2 > 0.0)) {
			clp -> x = MAX(p1 -> x, p2 -> x);
		}

		l1 = pt -> y - p1 -> y;
		l2 = pt -> y - p2 -> y;
		clp -> y = pt -> y;
		if ((l1 < 0.0) AND (l2 < 0.0)) {
			clp -> y = MIN(p1 -> y, p2 -> y);
		}
		if ((l1 > 0.0) AND (l2 > 0.0)) {
			clp -> y = MAX(p1 -> y, p2 -> y);
		}
		d   = RDIST(pt, clp);
		*d1 = RDIST(p1, clp);
		*d2 = RDIST(p2, clp);

	}

	if (cip -> metric EQ EUCLIDEAN) {

		/* Compute Steiner point for the three points and then
		   the distance from pt to this Steiner point */

		/* Make sure points are oriented nicely */
		if (left_turn(p1, p2, pt)) {
			a = p1; b = p2;
		}
		else {
			a = p2; b = p1;
		}

		eq_point (a, b, &e);
		eq_circle_center (a, b, &e, &ctr);

		if (right_turn (&e, a, pt)) {
			if (left_turn (&e, b, pt)) {
				if (sqr_dist(&ctr, pt) > sqr_dist(&ctr, &e)) {
					project_point(&e, &ctr, pt, clp);
					l = EDIST(&e, pt);
				}
				else {
					*clp = *pt;
					l = EDIST(a, pt) + EDIST(b, pt);
				}
			}
			else {
				*clp = *b;
				l = EDIST(a, b) + EDIST(pt, b);
			}
		}
		else {
			*clp = *a;
			l = EDIST(b, a) + EDIST(pt, a);
		}

		d   = l - EDIST(p1, p2);
		*d1 = EDIST(p1, clp);
		*d2 = EDIST(p2, clp);
	}


	return d;
}

/*
 * Find all FSTs whose removal splits the problem.  Make these
 * all required.
 */

	static
	int
bcc_find_required (

struct cinfo *	cip,		/* IN - compatibility info */
int *		made_req,	/* IN/OUT - list of FSTs made REQUIRED */
int		req_count	/* IN - existing required count */
)
{
int			i;
int			j;
int			nverts;
int			nedges;
int			kmasks;
int			nmasks;
int			max_stack;
struct bc3		bc;

	nverts	= cip -> num_verts;
	nedges	= cip -> num_edges;

	if (nverts < 2) {
		/* No BCC's. */
		return (req_count);
	}

	kmasks	= BMAP_ELTS (nverts);
	nmasks	= BMAP_ELTS (nedges);

	bc.cip		= cip;
	bc.kmasks	= kmasks;
	bc.nmasks	= nmasks;
	bc.dfs		= NEWA (nverts, int);
	bc.low		= NEWA (nverts, int);
	bc.parent	= NEWA (nverts, int);

	for (i = 0; i < nverts; i++) {
		bc.dfs	[i] = 0;
		bc.low [i] = 0;
		bc.parent [i] = -1;
	}

	j = (nedges > nverts) ? nedges : nverts;
	bc.max_stack	= j;
	bc.stack	= NEWA (j, int);
	bc.sp		= bc.stack;
	bc.counter	= 0;

	bc.bcc_vmask	= NEWA (kmasks, bitmap_t);
	bc.bcc_emask	= NEWA (nmasks, bitmap_t);
	bc.edges_seen	= NEWA (nmasks, bitmap_t);
	bc.bcc_vlist	= NEWA (nverts, int);
	bc.degree	= NEWA (nverts, int);
	bc.made_req	= made_req;
	bc.req_count	= req_count;

	for (i = 0; i < bc.kmasks; i++) {
		bc.bcc_vmask [i] = 0;
	}
	for (i = 0; i < bc.nmasks; i++) {
		bc.bcc_emask [i] = 0;
		bc.edges_seen [i] = 0;
	}

	/* Traverse each connected component, identifying its BCC's as	*/
	/* we go.							*/
	for (i = 0; i < nverts; i++) {
		if (NOT BITON (cip -> initial_vert_mask, i)) continue;
		if (bc.dfs [i] > 0) continue;

		/* Traverse one connected component, finding	*/
		/* each of its BCCs as we go.			*/

		bcc3 (&bc, i);
	}

	free ((char *) bc.degree);
	free ((char *) bc.bcc_vlist);
	free ((char *) bc.edges_seen);
	free ((char *) bc.bcc_emask);
	free ((char *) bc.bcc_vmask);
	free ((char *) bc.stack);
	free ((char *) bc.parent);
	free ((char *) bc.low);
	free ((char *) bc.dfs);

	return (bc.req_count);
}

/*
 * This is the recursive part of the bi-connected-components algorithm.	 It
 * is the standard method, with a few tweaks to work on hypergraphs instead.
 * We process each bi-connected component individually.
 */

	static
	void
bcc3 (

struct bc3 *		bcp,		/* IN - global BCC data */
int			v		/* IN - current DFS vertex */
)
{
int			i;
int			e;
int			e2;
int			w;
int *			ep1;
int *			ep2;
int *			vp1;
int *			vp2;
int *			vp3;
int *			vp4;
int *			stack_endp;
int *			sp;
int *			stack;
struct cinfo *		cip;

	cip = bcp -> cip;

	if ((v < 0) OR (v >= cip -> num_verts)) {
		fatal ("bcc3:  Bug 1.");
	}

	++(bcp -> counter);
	bcp -> dfs [v] = bcp -> counter;
	bcp -> low [v] = bcp -> counter;
	ep1 = cip -> term_trees [v];
	ep2 = cip -> term_trees [v + 1];
	while (ep1 < ep2) {
		e = *ep1++;
		if ((e < 0) OR (e >= cip -> num_edges)) {
			fatal ("bcc3: Bug 2.");
		}
		if (NOT BITON (cip -> initial_edge_mask, e)) continue;
		if (NOT BITON (bcp -> edges_seen, e)) {
			/* We haven't seen this edge before.  Push	*/
			/* it onto the stack...				*/
			stack_endp = &(bcp -> stack [bcp -> max_stack]);
			if ((bcp -> sp < bcp -> stack) OR
			    (bcp -> sp >= stack_endp)) {
				fatal ("bcc3: Bug 3.");
			}
			*(bcp -> sp)++ = e;
			SETBIT (bcp -> edges_seen, e);
		}
		/* Scan the vertices and process them... */
		vp1 = cip -> edge [e];
		vp2 = cip -> edge [e + 1];
		while (vp1 < vp2) {
			w = *vp1++;
			if ((w < 0) OR (w >= cip -> num_verts)) {
				fatal ("bcc3: Bug 4.");
			}
			if (bcp -> dfs [w] EQ 0) {
				bcp -> parent [w] = v;
				bcc3 (bcp, w);
				if (bcp -> low [w] >= bcp -> dfs [v]) {
					/* We have a new BCC! */
					stack	= bcp -> stack;
					sp	= bcp -> sp;
					do {
						if (sp <= stack) {
							fatal ("bcc3: Bug 5.");
						}
						e2 = *--sp;
					} while (e2 NE e);

					/* Process the bi-connected comp. */
					process_bcc3 (bcp, sp, bcp -> sp);

					/* Pop BCC edges from stack */
					bcp -> sp = sp;
				}
				if (bcp -> low [w] < bcp -> low [v]) {
					bcp -> low [v] = bcp -> low [w];
				}
			}
			else if ((w NE bcp -> parent [v]) AND
				 (bcp -> dfs [w] < bcp -> low [v])) {
				bcp -> low [v] = bcp -> dfs [w];
			}
		}
	}
}

/*
 * Process a single bi-connected component, specified as a list of edges.
 * We look for vertices having degree 1 within the component.  When this
 * happens, the incident edge is required.
 */

	static
	void
process_bcc3 (

struct bc3 *	bcp,		/* IN - global BCC data */
int *		edge_ptr,	/* IN - list of BCC edges */
int *		endp		/* IN - end of BCC edge list */
)
{
int		i;
int		e;
int		v;
int		n;
int *		ep1;
int *		ep2;
int *		ep3;
int *		vp1;
int *		vp2;
int *		vlp;
int *		elist;
struct cslist * cslp;
struct cinfo *	cip;
double		xe;

	cip = bcp -> cip;

	/* Gather a list of all vertices in this BCC.  Compute	*/
	/* their degrees (with respect to the BCC).		*/
	vlp = bcp -> bcc_vlist;

	ep1 = edge_ptr;
	ep2 = endp;
	while (ep1 < ep2) {
		e = *ep1++;
		SETBIT (bcp -> bcc_emask, e);
		vp1 = cip -> edge [e];
		vp2 = cip -> edge [e + 1];
		while (vp1 < vp2) {
			v = *vp1++;
			if (NOT BITON (cip -> initial_vert_mask, v)) continue;
			if (NOT BITON (bcp -> bcc_vmask, v)) {
				*vlp++ = v;
				SETBIT (bcp -> bcc_vmask, v);
				bcp -> degree [v] = 0;
			}
			++(bcp -> degree [v]);
		}
	}

	/* All of the vertices of this BCC are now known, as are their	*/
	/* degrees (relative to the component).	 Now look for vertices	*/
	/* of degree 1.							*/

	vp1 = bcp -> bcc_vlist;
	vp2 = vlp;
	while (vp1 < vp2) {
		v = *vp1++;
		CLRBIT (bcp -> bcc_vmask, v);		/* clean up as we go */

		n = bcp -> degree [v];
		if (n > 1) continue;

		ep1 = cip -> term_trees [v];
		ep2 = cip -> term_trees [v + 1];
		for (;;) {
			if (ep1 >= ep2) {
				fatal ("process_bcc3: Bug 1.");
			}
			e = *ep1++;
			if (BITON (bcp -> bcc_emask, e)) break;
		}
		if (BITON (cip -> required_edges, e)) continue;
		i = (bcp -> req_count)++;
		bcp -> made_req [i] = e;
		SETBIT (cip -> required_edges, e);
	}

	/* Clean up edge mark flags. */
	ep1 = edge_ptr;
	ep2 = endp;
	while (ep1 < ep2) {
		e = *ep1++;
		CLRBIT (bcp -> bcc_emask, e);
	}
}

/*
 * For sorting integers in place
 */

	static
	int
comp_ints (

const void *		p1,
const void *		p2
)
{
int			l1;
int			l2;

	l1 = *((int *) p1);
	l2 = *((int *) p2);

	if (l1 < l2) return (-1);
	if (l1 > l2) return (1);
	return (0);
}


/*
 * For sorting pruning graph in place
 */

	static
	int
comp_pg_edges (

const void *		p1,
const void *		p2
)
{
dist_t			l1;
dist_t			l2;

	l1 = ((struct pg_edge *) p1) -> len;
	l2 = ((struct pg_edge *) p2) -> len;

	if (l1 < l2) return (-1);
	if (l1 > l2) return (1);
	return (0);
}

/*
 * For sorting close terminals in place
 */

	static
	int
comp_clt (

const void *		p1,
const void *		p2
)
{
dist_t			l1;
dist_t			l2;

	l1 = ((struct clt_info *) p1) -> dist;
	l2 = ((struct clt_info *) p2) -> dist;

	if (l1 < l2) return (-1);
	if (l1 > l2) return (1);
	return (0);
}

/*
 * Compute the CPU time used since the last time we called this routine.
 */

	static
	cpu_time_t
get_delta_cpu_time (void)

{
cpu_time_t		now;
cpu_time_t		delta;

	now = get_cpu_time ();

	delta = now - Tn;

	Tn = now;

	return (delta);
}


/*
 * Compute and format the CPU time used since the last time we called
 * this routine.
 */

	static
	void
convert_delta_cpu_time (

char *		buf		/* OUT - ASCII time string, CPU seconds */
)
{
cpu_time_t	delta;

	delta = get_delta_cpu_time ();
	convert_cpu_time (delta, buf);
}