prunefst.c

生成直角Steiner树的程序包

	free ((char *) ltlist);
	free ((char *) tmask);
	free ((char *) fmask);
	free ((char *) flist);
	free ((char *) incompat);
	free ((char *) counts);
}

/*
 * Compute pruning information. For every FST identify a list
 * of "close" terminals and find their distance to the FST.
 */

	static
	void
compute_pruning_info (

struct cinfo *		cip,	/* IN	  - compatibility info */
struct pinfo *		pip	/* IN/OUT - pruning data structure  */
)
{
int			i;
int			j;
int			k;
int			t;
int			total;
int			steiner_index;
int			nedges;
int			nverts;
int			kmasks;
int *			vp1;
int *			vp2;
bitmap_t *		tmask;
struct full_set *	fsp;
struct pset *		terms;
struct clt_info*	cli;
struct clt_info*	clip;
struct point *		p1;
struct point *		p2;
dist_t			l;

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

	/* Terminal mask */
	tmask = NEWA (kmasks, bitmap_t);
	for (i = 0; i < kmasks; i++) {
		tmask [i] = 0;
	}

	/* Generate list of all edges in all FSTs */
	total = 0;
	for (i = 0; i < nedges; i++) {
		total += cip -> full_trees [i] -> nedges;
	}

	steiner_index = nverts;
	pip -> num_pg_edges	= total;
	pip -> pg_edges		= NEWA (total, struct pg_edge);
	memset (pip -> pg_edges, 0, total * sizeof (pip -> pg_edges [0]));
	pip -> steiner_index	= NEWA (nedges, int);
	k = 0;
	l = 0.0;
	for (i = 0; i < nedges; i++) {
		fsp = cip -> full_trees [i];
		terms = fsp -> terminals;
		pip -> steiner_index [i] = steiner_index;
		for (j = 0; j < fsp -> nedges; j++) {

			/* Compute length of this edge */
			p1 = (fsp -> edges[j].p1 < terms -> n)
			     ? &(terms -> a[ fsp -> edges[j].p1 ])
			     : &(fsp -> steiners -> a [fsp -> edges[j].p1 - terms -> n]);
			p2 = (fsp -> edges[j].p2 < terms -> n)
			     ? &(terms -> a[ fsp -> edges[j].p2 ])
			     : &(fsp -> steiners -> a [fsp -> edges[j].p2 - terms -> n]);
			if (cip -> metric EQ EUCLIDEAN) {
				l = EDIST(p1, p2) *
				    (1.0 - pip -> eps * ((double) cip -> edge_size [i]));
			}
			if (cip -> metric EQ RECTILINEAR) {
				l = RDIST(p1, p2);
			}

			pip -> pg_edges[k].fst = i;
			pip -> pg_edges[k].len = l;
			pip -> pg_edges[k].p1  =
				(fsp -> edges[j].p1 < terms -> n)
				? terms -> a[ fsp -> edges[j].p1 ].pnum
				: fsp -> edges[j].p1 - terms -> n + steiner_index;
			pip -> pg_edges[k].p2  =
				(fsp -> edges[j].p2 < terms -> n)
				? terms -> a[ fsp -> edges[j].p2 ].pnum
				: fsp -> edges[j].p2 - terms -> n + steiner_index;
			k++;
		}
		if (fsp -> steiners NE NULL) {
			steiner_index += fsp -> steiners -> n;
		}
	}
	pip -> num_pg_verts = steiner_index;

	/* Sort edge list */
	qsort(pip -> pg_edges, pip -> num_pg_edges,
	      sizeof(struct pg_edge), comp_pg_edges);

	/* Pruning graph has been constructed.
	   Now construct close terminal lists */

	pip -> clt_info	 = NEWA (nedges, struct clt_info *);
	pip -> clt_count = NEWA (nedges, int);
	cli = NEWA (nverts, struct clt_info);
	memset (cli, 0, nverts * sizeof (cli [0]));
	if (Print_Detailed_Timings) {
		fprintf(stderr, "- constructing close terminal list for each FST\n");
	}
	for (i = 0; i < nedges; i++) {

		/* Mark terminals in current FST */
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			t = *vp1++;
			SETBIT (tmask, t);
		}

		/* Find all close terminals and add information
		   to close terminal data structure */
		clip = cli;
		for (t = 0; t < nverts; t++) {
			if (NOT BITON (tmask, t)) {
				test_close_terminal(pip, i, t, &clip);
			}
		}

		/* Unmark terminals in current FST */
		vp1 = cip -> edge [i];
		while (vp1 < vp2) {
			t = *vp1++;
			CLRBIT (tmask, t);
		}

		/* Sort close terminals */
		pip -> clt_count [i] = clip - cli;
		pip -> clt_info [i]  = NULL;
		if (pip -> clt_count [i] > 0) {
			qsort(cli, pip -> clt_count [i],
			      sizeof(struct clt_info), comp_clt);
			pip -> clt_info [i] = NEWA (pip -> clt_count [i],
						    struct clt_info);
			for (j = 0; j < pip -> clt_count [i]; j++)
				pip -> clt_info [i][j] = cli[j];
		}
	}

	free ((char *) cli);
	free ((char *) tmask);
}
/*
 * Computes upper bounds for single FST or pair of FSTs.
 * Returns TRUE if all upper bound tests are passed
 * and FALSE otherwise.
 */

	static
	bool
passes_upper_bound_tests (

struct pinfo *		pip,	/* IN - pruning info */
struct bsd *		bsd,	/* IN - BSD data structure */
int			fst1,	/* IN - first FST */
int			fst2,	/* IN - second FST */
struct pset *		ltlist, /* IN - local terminal list (should just be allocated) */
bitmap_t *		ltmask, /* IN - local terminal mask (should just be cleared) */
int *			lflist, /* IN - local FST list (should just be allocated) */
bitmap_t *		lfmask	/* IN - local FST mask (should just be cleared) */
)
{
int			i;
int			j;
int			t;
int			hid;
int			nverts;
int			nedges;
int			kmasks;
int			nmasks;
int			lnedges;
int			lfs;
int			lfs2;
int			lsmtcount;
int			fstmaxind;
int			isterms;
bool			all_spanned;
bool			shorter_pair_found;
struct bbinfo *		bbip;
struct cinfo		lcinfo;
struct edge *		bsdmst;
struct edge *		ep;
int *			hep;
int *			ep1;
int *			ep2;
int *			vp1;
int *			vp2;
int *			lterm;
struct cinfo *		cip;
bitmap_t		smt_mask;
dist_t			l;
dist_t			bsdl;
dist_t			msthgl;
dist_t			minl;
dist_t			pairl;

	cip = pip -> cip;

	/* Construct list of terminals and set terminal mask */
	i = 0;
	vp1 = cip -> edge [fst1];
	vp2 = cip -> edge [fst1 + 1];
	while (vp1 < vp2) {
		t = *vp1++;
		ltlist -> a[i++] = cip -> pts -> a[t];
		SETBIT (ltmask, t);
	}

	if (fst2 >= 0) { /* negative means not defined */
		vp1 = cip -> edge [fst2];
		vp2 = cip -> edge [fst2 + 1];
		while (vp1 < vp2) {
			t = *vp1++;
			if (BITON (ltmask, t)) continue;
			ltlist -> a[i++] = cip -> pts -> a[t];
			SETBIT (ltmask, t);
		}
	}
	ltlist -> n = i;

	/* Reset terminal masks (in case we quit) */
	for (i = 0; i < ltlist -> n; i++) {
		CLRBIT (ltmask, ltlist -> a[i].pnum);
	}

	/* Lower bound on total length of these two FSTs (or this single FST)*/
	if (fst2 >= 0) {
		l = cip -> cost [fst1] + cip -> cost [fst2];
	}
	else {
		l = cip -> cost [fst1];
	}

	l *= (1.0 - pip -> eps * ((double) ltlist -> n));

	/* Compute BSD-MST */
	bsdmst = NEWA (ltlist -> n - 1, struct edge);
	if (bmst (ltlist, bsd, bsdmst) NE ltlist -> n - 1)
		fatal ("passes_upper_bound_test: bug 1");

	/* BSD-MST test... */
	bsdl = 0.0;
	ep   = bsdmst;
	for (i = 0; i < ltlist -> n - 1; i++, ep++) {
		bsdl += ep -> len;
	}
	if ((ltlist -> n > 3) AND (bsdl <= l)) {
		free ((char *) bsdmst);
		return FALSE;
	}

	/* Heuristic upper bound test (rectilinear) ... */
	if ((cip -> metric EQ RECTILINEAR) AND (ltlist -> n <= 15)) {
		if (kahng_robins_length (ltlist, l) < l) {
			free ((char *) bsdmst);
			return FALSE;
		}
	}

	/* Heuristic upper bound test (Euclidean) ... */
	if (cip -> metric EQ EUCLIDEAN) {
		if (greedy_heuristic (ltlist, bsd) < l) {
			free ((char *) bsdmst);
			return FALSE;
		}
	}

	/* If we are testing a pair of MST edges, we stop here */
	if (ltlist -> n EQ 3) {
		free ((char *) bsdmst);
		return TRUE;
	}

	/* Set terminal masks again */
	for (i = 0; i < ltlist -> n; i++) {
		SETBIT (ltmask, ltlist -> a[i].pnum);
	}

	/* Prepare for calling the branch-and-cut MSTHG procedure.	   */
	/* Construct list of FSTs spaning a subset of the given terminals. */
	/* We only need FSTs of cardinality 3 or larger.		   */

	lnedges = 0;
	for (i = 0; i < ltlist -> n; i++) {
		ep1 = cip -> term_trees [ ltlist -> a[i].pnum ];
		ep2 = cip -> term_trees [ ltlist -> a[i].pnum + 1 ];
		while (ep1 < ep2) {
			lfs = *ep1++;
			if (BITON (lfmask, lfs)) continue;
			if (NOT BITON (cip -> initial_edge_mask, lfs)) continue;
			if (cip -> edge_size [lfs] <= 2) continue;
			if ((lfs EQ fst1) AND (fst2 < 0)) continue; /* skip FST being tested */

			/* Does FST span a subset of given terminals? */
			all_spanned = TRUE;
			vp1 = cip -> edge [lfs];
			vp2 = cip -> edge [lfs + 1];
			while (vp1 < vp2) {
				t = *vp1++;
				if (NOT BITON (ltmask, t)) {
					all_spanned = FALSE;
					break;
				}
			}
			if (all_spanned) {
				SETBIT (lfmask, lfs);
				lflist [lnedges++] = lfs;
			}
		}
	}

	/* Reset terminal and FST masks */
	for (i = 0; i < ltlist -> n; i++) {
		CLRBIT (ltmask, ltlist -> a[i].pnum);
	}
	for (i = 0; i < lnedges; i++) {
		CLRBIT (lfmask, lflist [i]);
	}

	/* If no large FSTs then return (test cannot be performed) */
	if (lnedges EQ 0) {
		free ((char *) bsdmst);
		return TRUE;
	}

	/* Test if there is a PAIR of FSTs that has smaller or equal total length.   */
	/* In case the length is equal we only need to keep the "canonical"	     */
	/* par, i.e., for which the maximum index of (large) FSTs is minimized.	     */

	if (fst2 >= 0) {
		fstmaxind = -1;
		if ((cip -> edge_size [fst1] NE 2) AND (fst1 > fstmaxind)) fstmaxind = fst1;
		if ((cip -> edge_size [fst2] NE 2) AND (fst2 > fstmaxind)) fstmaxind = fst2;
	}
	else {
		fstmaxind = fst1;
	}

	nverts = ltlist -> n;
	shorter_pair_found = FALSE;
	for (i = 0; i < lnedges; i++) {
		lfs = lflist [i];

		/* Mark terminals in this FST */
		vp1 = cip -> edge [lfs];
		vp2 = cip -> edge [lfs + 1];
		while (vp1 < vp2) {
			t = *vp1++;
			SETBIT (ltmask, t);
		}

		if (cip -> edge_size [lfs] + 1 EQ nverts) {

			/* Find shortest BSD-MST edge to remaining terminal */
			ep   = bsdmst;
			minl = INF_DISTANCE;
			for (j = 0; j < nverts - 1; j++, ep++) {
				if ((NOT BITON (ltmask, ltlist -> a[ep -> p1].pnum)) OR
				    (NOT BITON (ltmask, ltlist -> a[ep -> p2].pnum))) {
					if (ep -> len < minl) minl = ep -> len;
				}
			}
			pairl = cip -> cost [lfs] + minl;
			if  (pairl <  l)			shorter_pair_found = TRUE;
			if ((pairl <= l) AND (lfs < fstmaxind)) shorter_pair_found = TRUE;
		}
		else {
			/* Try to combine with another FST */
			for (j = i+1; j < lnedges; j++) {
				lfs2  = lflist [j];
				if (cip -> edge_size [lfs] + cip -> edge_size [lfs2] - 1 NE nverts) continue;
				pairl = cip -> cost [lfs] + cip -> cost [lfs2];
				if (pairl > l) continue;

				/* Count intersecting terminals */
				isterms = 0;
				vp1 = cip -> edge [lfs2];
				vp2 = cip -> edge [lfs2 + 1];
				while (vp1 < vp2) {
					t = *vp1++;
					if (BITON (ltmask, t)) isterms++;
				}
				if (isterms NE 1) continue; /* not spanning all terminals */

				if  (pairl <  l)	    shorter_pair_found = TRUE;
				if ((pairl <= l)       AND
				    (lfs  < fstmaxind) AND
				    (lfs2 < fstmaxind))	    shorter_pair_found = TRUE;
			}
		}

		/* Reset terminal mask */
		for (j = 0; j < nverts; j++) {
			CLRBIT (ltmask, ltlist -> a[j].pnum);
		}
		if (shorter_pair_found) break;
	}

	if (shorter_pair_found) {
		free ((char *) bsdmst);
		return FALSE;
	}

	/* Allocate and initialize local cinfo structure */
	nedges =  lnedges + ltlist -> n - 1; /* we add BSD-MST edges */
	kmasks =  BMAP_ELTS (nverts);
	nmasks =  BMAP_ELTS (nedges);

	memset (&lcinfo, 0, sizeof (lcinfo));
	lcinfo.num_verts		= nverts;
	lcinfo.num_edges		= nedges;
	lcinfo.num_vert_masks		= kmasks;
	lcinfo.initial_vert_mask	= NEWA (kmasks, bitmap_t);
	lcinfo.num_edge_masks		= nmasks;
	lcinfo.initial_edge_mask	= NEWA (nmasks, bitmap_t);
	lcinfo.required_edges		= NEWA (nmasks, bitmap_t);
	lcinfo.edge			= NEWA (nedges + 1, int *);
	lcinfo.edge[0]			= NEWA (cip -> edge[cip -> num_edges] -
						cip -> edge[0], int);
	lcinfo.edge_size		= NEWA (nedges, int);
	lcinfo.cost			= NEWA (nedges, dist_t);
	lcinfo.tflag			= NEWA (nverts, bool);
	lcinfo.metric			= PURE_GRAPH;
	lcinfo.scale.scale		= 0;
	lcinfo.scale.scale_mul		= 1.0;
	lcinfo.scale.scale_div		= 1.0;
	lcinfo.integrality_delta	= 0;
	lcinfo.mst_length		= 0;
	lcinfo.description		= gst_strdup ("Local MSTHG");
	lterm				= NEWA (cip -> num_verts, int);

	for (i = 0; i < kmasks; i++) {
		lcinfo.initial_vert_mask [i] = 0;
	}
	for (i = 0; i < nverts; i++) {
		SETBIT (lcinfo.initial_vert_mask, i);
		lcinfo.tflag [i] = TRUE; /* we only meed to solve MSTHG's */
		lterm [ ltlist -> a[i].pnum ] = i;
	}
	for (i = 0; i < nmasks; i++) {
		lcinfo.initial_edge_mask [i] = 0;
		lcinfo.required_edges [i] = 0;
	}
	for (i = 0; i < nedges; i++) {
		SETBIT (lcinfo.initial_edge_mask, i);
	}

	hid   = 0;
	hep   = lcinfo.edge [0];

	/* Add all identified FSTs to hyperedge list */
	for (i = 0; i < lnedges; i++) {
		lfs = lflist [i];

		/* Make it more attractive to seek an MST with many FSTs */
		lcinfo.cost [hid]      = cip -> cost [lfs] -
					 cip -> integrality_delta / (double) nverts;
		lcinfo.edge_size [hid] = cip -> edge_size [lfs];
		lcinfo.edge [hid]      = hep;
		vp1 = cip -> edge [lfs];
		vp2 = cip -> edge [lfs + 1];
		while (vp1 < vp2) {
			*hep++ = lterm[ *vp1++ ];
		}
		hid++;
	}

	/* Add all BSD-MST edges to hyperedge list */
	ep = bsdmst;
	for (i = 0; i < ltlist -> n - 1; i++, ep++) {
		lcinfo.cost [hid]      = ep -> len -
					 cip -> integrality_delta / (double) nverts;
		lcinfo.edge_size [hid] = 2;
		lcinfo.edge [hid]      = hep;
		*(hep++) = ep -> p1;
		*(hep++) = ep -> p2;
		hid++;
	}
	lcinfo.edge [hid] = hep;

	/* Call branch-and-cut procedure */

	tracef_control.disabled = TRUE;
	init_term_trees (&lcinfo);
	lcinfo.inc_edges = compute_basic_incompat (&lcinfo);

	if ((lcinfo.num_verts <= 24) AND (lcinfo.num_edges <= 32)) {

		/* Small problem -- use backtrack search. */

		smt_mask = 0;
		msthgl = backtrack_search (&lcinfo, &smt_mask);
		lsmtcount = 0;
		for (i = 0; i < lcinfo.num_edges; i++) {
			if ((smt_mask & (1 << i)) NE 0) {
				++lsmtcount;
			}
		}
	}
	else {

		bbip = create_bbinfo (&lcinfo);

#if 0
		/* Print problem */
		{
		int * ep1;
		int * ep2;
		printf("\n%d %d\n", lcinfo.num_verts, lcinfo.num_edges);
		for (i = 0; i < lcinfo.num_edges; i++) {
			ep1 = lcinfo.edge[i];
			ep2 = lcinfo.edge[i+1];
			while (ep1 < ep2)
				printf("%d ", *(ep1++) + 1);
			printf("%f\n", lcinfo.cost[i]);
		}
		printf("%d\n", lcinfo.num_verts);
		for (i = 0; i < lcinfo.num_verts; i++)
			printf("%d\n", i + 1);
		}
#endif

		/* Find MST in hypergraph. */

		msthgl = branch_and_cut (bbip);

		/* Count number of FSTs in SMT */
		lsmtcount = 0;
		for (i = 0; i < lcinfo.num_edges; i++) {
			if (BITON (bbip->smt, i)) lsmtcount++;
		}

		/* Free allocated memory */
		destroy_bbinfo (bbip);

	}

	free_phase_1_data (&lcinfo);
	free ((char *) lterm);
	free ((char *) bsdmst);

	/* Return result of final test */
	if (fst2 >= 0) {
		if (msthgl < l - cip -> integrality_delta) {
			return FALSE;
		}
		/* If equal length then there should be at least three FSTs */
		else if ((msthgl < l) AND (lsmtcount >= 3)) {
			return FALSE;
		}
	}
	else {
		if (msthgl <= l) {
			return FALSE;
		}
	}

	return TRUE; /* all tests passed */

}

/*
 * Find distance and attachment terminals for given terminal
 * to a specific FST
 */