prunefst.c

生成直角Steiner树的程序包

			    (NOT BITON (cip -> required_edges, fsave))) {

				if (NOT BITON (cip -> initial_edge_mask, fsave)) {
					/* fatal: FST already removed */
					fatal("prune_fsts: bug 1.");
				}

				SETBIT (cip -> required_edges, fsave);
				required_total++;
				made_req [req_count++] = fsave;
			}
		}

#if 1
		i = req_count;
		req_count = bcc_find_required (cip, made_req, req_count);
		required_total += (req_count - i);
#else
		/* Test if leaving out an FST disconnects the remaining FSTs */
		/* CHANGE! This is rather time-consuming. Should be implemented using bicomps */

		for (i = 0; i < nedges; i++) {
			if ((BITON (cip -> initial_edge_mask, i)) AND
			   (NOT BITON (cip -> required_edges, i))) {

				/* create union-find data structure */
				dsuf_create (&del_comps, cip -> num_verts);
				for (t = 0; t < cip -> num_verts; t++) {
					dsuf_makeset (&del_comps, t);
				}

				del_comps_count = cip -> num_verts;

				for (j = 0; j < nedges; j++) {
					if (BITON (cip -> initial_edge_mask, j)) {

						if (j EQ i) continue; /* skip FST being tested */

						/* Unite vertices spanned */
						vp1 = cip -> edge [j];
						vp2 = cip -> edge [j + 1] - 1;
						while (vp1 < vp2) {
							r1 = dsuf_find (&del_comps, *vp1);
							r2 = dsuf_find (&del_comps, *vp2);

							if (r1 NE r2) {
								dsuf_unite (&del_comps, r1, r2);
								del_comps_count--;
							}
							vp1++;
						}
						if (del_comps_count EQ 1) break;
					}
				}

				if (del_comps_count > 1) {

					/* FST must be included in SMT */

					SETBIT (cip -> required_edges, i);
					required_total++;
					made_req [req_count++] = i;
				}
				dsuf_destroy (&del_comps);
			}
		}
#endif

		/* Now update data structures for new required FSTs */
		for (j = 0; j < req_count; j++) {

			fsave = made_req [j];

			/* Unite vertices spanned */
			vp1 = cip -> edge [fsave];
			vp2 = cip -> edge [fsave + 1] - 1;
			while (vp1 < vp2) {
				r1 = dsuf_find (&comps, *vp1);
				r2 = dsuf_find (&comps, *vp2);

				if (r1 EQ r2) { /* fatal: cycle created */
					fatal("prune_fsts: bug 2.");
				}
				dsuf_unite (&comps, r1, r2);
				vp1++;
			}

			/* Prune all incompatible FSTs */
			ep1 = cip -> inc_edges [fsave];
			ep2 = cip -> inc_edges [fsave + 1];
			while (ep1 < ep2) {
				i = *ep1++;
				if (BITON (cip -> initial_edge_mask, i)) {

					if (BITON (cip -> required_edges, i)) {
						/* fatal: FST already required */
						fatal("prune_fsts: bug 3.");
					}
					CLRBIT (cip -> initial_edge_mask, i);
					pruned_total++;
				}
			}
		}

		/* Remove FSTs making cycles among required FSTs */
		for (i = 0; i < nedges; i++) {
			if ((BITON (cip -> initial_edge_mask, i)) AND
			    (NOT BITON (cip -> required_edges, i))) {

				/* Check if a pair of vertices span the same component */
				vp1 = cip -> edge [i];
				vp2 = cip -> edge [i + 1];
				while (vp1 < vp2) {
					vp = vp1 + 1;
					while (vp < vp2) {
						r1 = dsuf_find (&comps, *vp1);
						r2 = dsuf_find (&comps, *vp);

						if (r1 EQ r2) { /* cycle created - remove FST */
							CLRBIT (cip -> initial_edge_mask, i);
							pruned_total++;
							vp = vp1 = vp2;
							break;
						}
						vp++;
					}
					vp1++;
				}
			}
		}

		if (req_count > 0) goto try_again;

		if (old_pruned_total EQ pruned_total) break;

		if (Print_Detailed_Timings) {
			fprintf(stderr, "- scan %d finished. %6d FSTs pruned\n",
				scan, pruned_total - old_pruned_total);
		}
	}

	if (Print_Detailed_Timings) {
		fprintf(stderr, "- pruning finished: before: %d  after: %d  required: %d\n",
			nedges, nedges - pruned_total, required_total);
	}

	/* Free pruning data information... */

	for (i = 0; i < nedges; i++) {
		if (pinfo.clt_info[i] NE NULL) {
			free ((char *) pinfo.clt_info [i]);
		}
	}
	free ((char *) pinfo.clt_count);
	free ((char *) pinfo.clt_info);
	free ((char *) pinfo.steiner_index);
	free ((char *) pinfo.pg_edges);

	dsuf_destroy (&comps);

	free ((char *) made_req);
	free ((char *) comps_edge);

	free ((char *) lfmask);
	free ((char *) lflist);
	free ((char *) ltmask);
	free ((char *) ltlist);

	free ((char *) pinfo.compat_mask);
}

/*
 * Remove FSTs permanently that have been marked as not needed.
 * However, MST edges (2-terminal FSTs) are not deleted even if
 * marked as never needed.
 * It should be noted that no arrays are reallocated; data is packed in place.
 */
	static
	void
zap_deleted_fsts (

struct cinfo *	cip		/* IN/OUT - compatibility info */
)
{
int			i;
int			j;
int			k;
int			new_nedges;
int *			ni;
int *			vp;
int *			vp1;
int *			vp2;
int *			ep;
int *			ep1;
int *			ep2;
struct full_set *	fsp;

	/* First we count the number of FSTs that remain */
	/* and set up map from old to new edge index */

	new_nedges = 0;
	ni  = NEWA (cip -> num_edges, int);
	for (i = 0; i < cip -> num_edges; i++) {
		if ((BITON (cip -> initial_edge_mask, i)) OR (cip -> edge_size [i] EQ 2))
			ni[i] = new_nedges++;
		else
			ni[i] = -1;
	}

	/* Pack edge_size and cost arrays */
	for (i = 0; i < cip -> num_edges; i++) {
		if (ni[i] >= 0) {
			cip -> edge_size [ ni[i] ] = cip -> edge_size [i];
			cip -> cost	 [ ni[i] ] = cip -> cost [i];
		}
	}

	/* Pack edge array */
	vp = cip -> edge [0];
	for (i = 0; i < cip -> num_edges; i++) {
		if (ni[i] >= 0) {
			vp1 = cip -> edge [i];
			vp2 = cip -> edge [i + 1];
			cip -> edge [ ni[i] ] = vp;
			while (vp1 < vp2)
				*(vp++) = *(vp1++);
		}
	}
	cip -> edge [ new_nedges ] = vp;

	/* Pack term_trees array */
	ep = cip -> term_trees [0];
	for (j = 0; j < cip -> num_verts; j++) {
		vp1 = cip -> term_trees [j];
		vp2 = cip -> term_trees [j + 1];
		cip -> term_trees [j] = ep;
		while (vp1 < vp2) {
			i = *(vp1++);
			if (ni[i] >= 0)
				*(ep++) = ni[i];
		}
	}
	cip -> term_trees[ cip -> num_verts ] = ep;

	/* Pack inc_edges array */
	ep = cip -> inc_edges [0];
	for (i = 0; i < cip -> num_edges; i++) {
		if (ni[i] >= 0) {
			ep1 = cip -> inc_edges [i];
			ep2 = cip -> inc_edges [i + 1];
			cip -> inc_edges [ ni[i] ] = ep;
			while (ep1 < ep2) {
				k = *(ep1++);
				if (ni[k] >= 0)
					*(ep++) = ni [k];
			}
		}
	}
	cip -> inc_edges [ new_nedges ] = ep;

	/* Pack full_trees array */
	if (cip -> full_trees NE NULL) {
		for (i = 0; i < cip -> num_edges; i++) {
			if (ni[i] >= 0) {
				cip -> full_trees [ ni[i] ] = cip -> full_trees [i];
				cip -> full_trees [ ni[i] ] -> tree_num = ni[i];
			}
			else {

				/* Free FST data */
				fsp = cip -> full_trees [i];
				free ((char *) (fsp -> terminals));
				free ((char *) (fsp -> steiners));
				free ((char *) (fsp -> edges));
				free ((char *) fsp);
			}
		}
	}

	/* Pack bit maps */
	for (i = 0; i < cip -> num_edges; i++) {
		if (ni[i] >= 0) {
			if (BITON (cip -> initial_edge_mask, i))
				SETBIT (cip -> initial_edge_mask, ni[i]);
			else
				CLRBIT (cip -> initial_edge_mask, ni[i]);

			if (BITON (cip -> required_edges, i))
				SETBIT (cip -> required_edges, ni[i]);
			else
				CLRBIT (cip -> required_edges, ni[i]);
		}
	}

	/* Finally set edge count */
	cip -> num_edges      = new_nedges;
	cip -> num_edge_masks = BMAP_ELTS (new_nedges);

	free((char *) ni);
}

/*
 * Check if a given FST can be pruned
 */
	static
	bool
prune_this_fst (

struct cinfo *	cip,		/* IN - compatibility info */
struct pinfo *	pip,		/* IN - pruning data structure */
int		fst
)
{
int			i;
int			clt_count;
int			curr_clt;
int			root;
int			root1;
int			root2;
bool			prune_fst;
struct dsuf		comps;
struct pg_edge *	pg_edge;
struct clt_info *	clt;

	clt_count = pip -> clt_count [fst];
	if (clt_count == 0) return FALSE;

	/* Create disjoint set */
	dsuf_create (&comps, pip -> num_pg_verts + 1);
	for (i = 0; i < pip -> num_pg_verts + 1; i++)
		dsuf_makeset (&comps, i);

	/* Add pruning graph edges (in sorted order) */

	prune_fst = FALSE;
	curr_clt  = 0;
	clt = &(pip -> clt_info [fst][curr_clt]);
	for (i = 0; i < pip -> num_pg_edges; i++) {
		pg_edge = &(pip -> pg_edges[i]);
		while (pg_edge -> len > clt -> dist) {
			root  = dsuf_find (&comps, clt -> term);
			root1 = dsuf_find (&comps, clt -> aterm1);
			root2 = dsuf_find (&comps, clt -> aterm2);

			if ((root NE root1) AND (root NE root2)) {
				prune_fst = TRUE; /* This FST can be pruned! */
				goto prune_exit;
			}

			curr_clt++;
			if (curr_clt >= clt_count)
				goto prune_exit; /* This FST cannot be pruned... */
			clt = &(pip -> clt_info [fst][curr_clt]);
		}

		/* Add edge if FST is not deleted or incompatible */
		if (BITON (pip -> compat_mask, pg_edge -> fst)) {
			root1 = dsuf_find (&comps, pg_edge -> p1);
			root2 = dsuf_find (&comps, pg_edge -> p2);
			if (root1 NE root2)
				dsuf_unite (&comps, root1, root2);
		}

	}
	prune_exit:

	dsuf_destroy (&comps);
	return prune_fst;
}

/*
 * Add a pair of FSTs as incompatible
 */

	static
	void
add_incompat (

struct incompat **	incompat,	/* IN/OUT - incomp. data structure */
int			fst1,		/* IN	  - first FST */
int			fst2,		/* IN	  - second FST */
int*			counts		/* IN/OUT - incomp. counts */
)
{
struct incompat *	icp;

	icp = incompat [fst1];
	incompat [fst1] = NEW (struct incompat);
	incompat [fst1] -> fst	= fst2;
	incompat [fst1] -> next = icp;
	counts [fst1]++;

	icp = incompat [fst2];
	incompat [fst2] = NEW (struct incompat);
	incompat [fst2] -> fst	= fst1;
	incompat [fst2] -> next = icp;
	counts [fst2]++;
}

/*
 * Computes for each FST, a list of those FSTs which are incompatible,
 * that is, fulfill one of the following conditions:
 * 1. The FSTs have two ore more terminals in common
 * 2. The FSTs have one terminal in common and the BSD of their
 * terminals is shorter than the total length of the FSTs
 * 3. An heuristic tree spanning the terminals is shorter.
 * 4. Another pair of FSTs has shorther or equal length.
 * 5. The MSTHG for the FSTs within the terminals spanned
 *    together with BSD-MST edges is shorter.
 */

	static
	void
compute_incompatibility_info (

struct pinfo *		pip,	/* IN/OUT - pruning info */
struct bsd *		bsd	/* IN	  - BSD data structure	*/
)
{
int			i;
int			j;
int			k;
int			t;
int			fs;
int			common;
int			comterm;
int			nverts;
int			nedges;
int			kmasks;
int			nmasks;
int			total;
bitmap_t *		fmask;
bitmap_t *		lfmask;
bitmap_t *		edge_mask;
bitmap_t *		tmask;
bitmap_t *		ltmask;
struct pset *		ltlist;
struct incompat *	icp;
struct incompat *	icpn;
struct incompat **	incompat;
int *			ep1;
int *			ep2;
int *			vp1;
int *			vp2;
int *			flist;
int *			lflist;
int **			inc_edges;
int *			counts;
struct cinfo *		cip;

	/* Initialize and allocate various variables and arrays */
	cip	  = pip -> cip;
	nverts	  = cip -> num_verts;
	kmasks	  = cip -> num_vert_masks;
	nedges	  = cip -> num_edges;
	nmasks	  = cip -> num_edge_masks;
	edge_mask = cip -> initial_edge_mask;

	inc_edges = NEWA (nedges + 1, int *);
	counts	  = NEWA (nedges, int);
	incompat  = NEWA (nedges, struct incompat *);
	flist	  = NEWA (nedges, int);
	fmask	  = NEWA (nmasks, bitmap_t);
	tmask	  = NEWA (kmasks, bitmap_t);
	ltlist	  = NEW_PSET (nverts);
	ltmask	  = NEWA (kmasks, bitmap_t);
	lflist	  = NEWA (nedges, int);
	lfmask	  = NEWA (nmasks, bitmap_t);

	for (i = 0; i < kmasks; i++) {
		tmask [i] = 0;
		ltmask [i] = 0;
	}
	for (i = 0; i < nmasks; i++) {
		fmask [i] = 0;
		lfmask [i] = 0;
	}
	for (i = 0; i < nedges; i++) {
		incompat [i] = NULL;
		counts [i]   = 0;
	}

	/* Compute the list of (lists of) incomatible FSTs... */
	total = 0;

	if (Print_Detailed_Timings) {
		fprintf(stderr, "- computing incompatible FSTs for each FST\n");
	}
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;

		/* Develop list of all FSTs adjacent to FST i... */
		k = 0; j = 0;
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			t = *vp1++;
			SETBIT (tmask, t);
			ep1 = cip -> term_trees [t];
			ep2 = cip -> term_trees [t + 1];
			while (ep1 < ep2) {
				fs = *ep1++;
				if (BITON (fmask, fs)) continue;
				if (NOT BITON (edge_mask, fs)) continue;
				if (fs <= i) continue; /* test pairs once */
				SETBIT (fmask, fs);
				flist [k++] = fs;
			}
		}

		/* Now loop through all adjacent FSTs */
		ep1 = &flist [0];
		ep2 = &flist [k];
		while (ep1 < ep2) {
			fs = *ep1++;
			CLRBIT (fmask, fs);

			/* Count number of vertices in common. */
			common = 0; comterm = 0;
			vp1 = cip -> edge [fs];
			vp2 = cip -> edge [fs + 1];
			while (vp1 < vp2) {
				t = *vp1++;
				if (BITON (tmask, t)) {
					++common; comterm = t;
				}
			}
			if (common >= 2) {
				/* Too many - retain as incompatible */
				add_incompat(incompat, i, fs, counts); total += 2;
				continue;
			}

			/* One terminal in common. Perform thorough upper tests. */
			if (NOT passes_upper_bound_tests (pip, bsd, i, fs,
							  ltlist, ltmask, lflist, lfmask)) {
				/* Did not pass - retain as incompatible */
				add_incompat(incompat, i, fs, counts); total += 2;
				continue;
			}
		}

		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;
			CLRBIT (tmask, j);
		}
	}

	/* Now allocate and copy into contiguous memory... */
	ep1 = NEWA (total, int);
	for (i = 0; i < nedges; i++) {
		inc_edges [i] = ep1;
		k = counts [i];
		if (k <= 0) continue;
		icp = incompat [i];
		while (icp NE NULL) {
			*ep1++ = icp -> fst;
			icpn = icp -> next;
			free ((char *) icp);
			icp = icpn;
		}
		qsort(inc_edges [i], k, sizeof(int), comp_ints);
	}
	inc_edges [i] = ep1;

	if (ep1 - inc_edges [0] NE total) {
		fatal ("compute_incompatibility_info: bug 1.");
	}
	if (cip -> inc_edges NE NULL) {
		if (cip -> inc_edges [0] NE NULL) {
			free ((char *) (cip -> inc_edges [0]));
		}
		free ((char *) (cip -> inc_edges));
	}
	cip -> inc_edges = inc_edges;

	/* Free allocated memory */

	free ((char *) lfmask);
	free ((char *) lflist);
	free ((char *) ltmask);