prunefst.c

生成直角Steiner树的程序包

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

	File:	prunefst.c
	Rev:	a-1
	Date:	11/22/2000

	Copyright (c) 2000, 2001 by Martin Zachariasen & David M. Warme

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

	Pruning of Euclidean and rectilinear FSTs using method
	originally proposed by Fossmeier & Kaufmann.
	Use implementation similar to the one suggested by Althaus,
	but with significant improvements in the compatibility tests.

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

	Modification Log:

	a-1:	11/22/2000	martinz
		: Created.

************************************************************************/

#include "bb.h"
#include "bsd.h"
#include "btsearch.h"
#include "dsuf.h"
#include "efuncs.h"
#include "emptyr.h"
#include "p1io.h"
#include "steiner.h"

/*
 * Global Routines
 */

int			main (int, char **);

/*
 * External References
 */

	/* none */

/*
 * Local Macros
 */

#ifndef MIN
 #define MIN(a,b)	(((a) < (b)) ? (a) : (b))
#endif
#ifndef MAX
 #define MAX(a,b)	(((a) > (b)) ? (a) : (b))
#endif


/*
 * Local Types
 */

struct pg_edge {
	dist_t				len; /* Length of edge */
	int				fst; /* Corresponding FST */
	int				p1;  /* First vertex */
	int				p2;  /* Second vertex */
};

struct clt_info {
	dist_t				dist;
	int				term;
	int				aterm1;
	int				aterm2;
};

struct pinfo {
	struct cinfo *			cip;
	int				num_pg_edges;
	int				num_pg_verts;
	struct pg_edge *		pg_edges;
	int *				steiner_index;
	struct clt_info **		clt_info;
	int *				clt_count;
	bitmap_t *			compat_mask;
	double				eps;
};

struct incompat {
	struct incompat *		next;
	int				fst;
};

struct bc3 {
	struct cinfo *	cip;		/* problem data */
	int		kmasks;		/* size of vert_mask */
	int		nmasks;		/* size of edge_mask */
	int *		dfs;		/* DFS number of each vertex */
	int *		low;		/* lowest DFS num in component */
	int *		parent;		/* parents of vertices in DFS tree */
	int		max_stack;	/* size of stack */
	int *		stack;		/* base-address of edge stack */
	int *		sp;		/* current stack pointer */
	int		counter;	/* DFS number generator */
	bitmap_t *	bcc_vmask;	/* scratch buffer for new BCCs */
	bitmap_t *	bcc_emask;	/* scratch buffer for new BCCs */
	bitmap_t *	edges_seen;	/* edges already pushed */
	int *		bcc_vlist;	/* scratch vertex list buffer */
	int *		degree;		/* temp vertex degree counter */
	int *		made_req;	/* caller's list of required edges */
	int		req_count;	/* cur index into made-req */
};


/*
 * Local Routines
 */

static void		add_incompat (struct incompat **, int, int, int *);
static int		bcc_find_required (struct cinfo *, int *, int);
static void		bcc3 (struct bc3 *, int);
static int		comp_ints (const void *, const void *);
static int		comp_pg_edges (const void *, const void *);
static int		comp_clt (const void *, const void *);
static void		compute_incompatibility_info (struct pinfo *, struct bsd *);
static void		compute_pruning_info (struct cinfo *, struct pinfo *);
static void		convert_delta_cpu_time (char *);
static void		decode_params (int, char**);
static cpu_time_t	get_delta_cpu_time (void);
static bool		passes_upper_bound_tests (struct pinfo *, struct bsd *, int, int,
						  struct pset *, bitmap_t *, int *, bitmap_t *);
static void		process_bcc3 (struct bc3 *, int *, int *);
static void		prune_fsts(struct cinfo *, struct bsd *, double);
static bool		prune_this_fst (struct cinfo *, struct pinfo *, int);
static dist_t		terminal_edge_distance(struct cinfo *, struct point *,
						struct point *, struct point *,
						struct point *,
						dist_t *, dist_t *);
static void		test_close_terminal (struct pinfo *,
					     int, int, struct clt_info **);
static void		usage ();
static void		zap_deleted_fsts(struct cinfo *);


/*
 * Local Variables
 */

static char *		description = NULL;
static int		EpsilonFactor = 32;
static char *		me;
static int		output_version = CURRENT_P1IO_VERSION;
static bool		Print_Detailed_Timings = FALSE;
static cpu_time_t	T0;
static cpu_time_t	Tn;


/*
 * Function Prototypes.
 */

extern int		bmst (struct pset *, struct bsd *, struct edge *);
extern dist_t		smith_lee_liebman(struct pset *);
extern dist_t		greedy_heuristic(struct pset *, struct bsd *);

/*
 * The main routine for the "prunefst" program.	 It takes the output from
 * the Euclidean or rectilinear FST generator, prunes FSTs that cannot
 * appear in an optimal solution and outputs the pruned set of FSTs
 * to the backend (FST concatenation procedure).
 */

	int
main (

int		argc,
char **		argv
)
{
int			nedges;
int			fpsave;
char			buf1 [32];
struct cinfo		cinfo;
struct bsd *		bsd;
struct edge *		mst_edges;
cpu_time_t		Tzap;
bitmap_t *		empty_rect;
double			eps;

	fpsave = set_floating_point_double_precision ();

	setbuf (stdout, NULL);

	decode_params (argc, argv);

	init_tables ();

	/* Read FST data */
	read_phase_1_data (&cinfo);

	if ((cinfo.metric NE RECTILINEAR) AND (cinfo.metric NE EUCLIDEAN)) {
		fprintf (stderr, "Can only prune geometric (Euclidean or rectilinear) FSTs\n");
		exit (1);
	}

	startup_lp_solver ();

	initialize_btsearch ();

	T0 = get_cpu_time ();
	Tn = T0;

	/* Compute minimum spanning tree */

	mst_edges = NEWA (cinfo.pts -> n, struct edge);
	eps = 0.0;
	nedges	  = 0;
	if (cinfo.metric EQ EUCLIDEAN) {
		eps	= ((double) EpsilonFactor) * DBL_EPSILON;
		nedges	= euclidean_mst (cinfo.pts, mst_edges);
	}
	else if (cinfo.metric EQ RECTILINEAR) {
		empty_rect = init_empty_rectangles (cinfo.pts, NULL);
		nedges	   = rect_mst (cinfo.pts, mst_edges, empty_rect);
		free ((char *) empty_rect);
	}
	if (nedges NE cinfo.pts -> n - 1) {
		fatal ("prunefst: bug 1.");
	}

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute MST:            %s\n", buf1);
	}
	/* Compute bottleneck Steiner distances */

	bsd = compute_bsd (nedges, mst_edges, 0);
	free ((char *) mst_edges);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Compute BSD:            %s\n", buf1);
	}

	/* Prune FSTs */
	prune_fsts (&cinfo, bsd, eps);

	if (Print_Detailed_Timings) {
		convert_delta_cpu_time (buf1);
		fprintf (stderr, "Pruning FSTs:           %s\n", buf1);
	}

	/* Remove deleted FSTs permanently */
	zap_deleted_fsts (&cinfo);

	/* Measure zap time.  This also sets Tn so that Tn-T0 is   */
	/* the total processing time.				   */
	Tzap = get_delta_cpu_time ();

	if (Print_Detailed_Timings) {
		convert_cpu_time (Tzap, buf1);
		fprintf (stderr, "Zap deleted FSTs:       %s\n", buf1);
		convert_cpu_time (Tn - T0, buf1);
		fprintf (stderr, "Total:                  %s\n", buf1);
	}
	cinfo.p1time += (Tn - T0);
	if (description NE NULL) {
		free ((char *) cinfo.description);
		cinfo.description = gst_strdup (description);
	}

	/* Print pruned FST data */
	print_phase_1_data (&cinfo, output_version);

	shutdown_bsd (bsd);

	shutdown_lp_solver ();

	free_phase_1_data (&cinfo);

	restore_floating_point_precision (fpsave);

	exit (0);
}

/*
 * This routine decodes the various command-line arguments.
 */

	static
	void
decode_params (

int		argc,
char **		argv
)
{
char *		ap;
char		c;

	--argc;
	me = *argv++;
	while (argc > 0) {
		ap = *argv++;
		if (*ap NE '-') {
			usage ();
		}
		++ap;
		while ((c = *ap++) NE '\0') {
			switch (c) {
			case 'd':
				if (*ap EQ '\0') {
					if (argc <= 0) {
						usage ();
					}
					ap = *argv++;
					--argc;
				}
				if (strlen (ap) >= 80) {
					fprintf (stderr,
						"Description must be less"
						" than 80 characters.\n");
					usage ();
				}
				description = ap;
				/* Change newlines to spaces... */
				for (;;) {
					ap = strchr (ap, '\n');
					if (ap EQ NULL) break;
					*ap++ = ' ';
				}
				ap = "";
				break;

			case 'f':
				/* Epsilon multiplication factor */
				if (*ap EQ '\0') {
					if (argc <= 0) {
						usage ();
					}
					ap = *argv++;
					--argc;
				}
				EpsilonFactor = atoi (ap);
				ap = "";
				break;
			case 't':
				Print_Detailed_Timings = TRUE;
				break;
			default:
				usage ();
				break;
			}
		}
		--argc;
	}
}

/*
 * This routine prints out the proper usage and exits.
 */

static char *	arg_doc [] = {
	"",
	"\t-d txt\tDescription of problem instance.",
	"\t-f F\tEpsilon multiplication factor F for floating point number",
	"\t\tcomparisons (default: 32).",
	"\t-t\tPrint detailed timings on stderr.",
	"",
	NULL
};

	static
	void
usage (void)

{
char **		pp;
char *		p;

	(void) fprintf (stderr,
			"\nUsage: %s [-t] [-d description] [-f F]",
			me);
	pp = &arg_doc [0];
	while ((p = *pp++) NE NULL) {
		(void) fprintf (stderr, "%s\n", p);
	}
	exit (1);
}

/*
 * Main pruning procedure. We use a method proposed by Fossmeier and Kaufmann
 * based on thesing whether is advantageous to extend a given FST with a
 * terminal no currently spanned. If so, the FST is discarded.
 */
	static
	void
prune_fsts (

struct cinfo *		cip,	/* IN/OUT - compatibility info */
struct bsd *		bsd,	/* IN	  - BSD data structure	*/
double			eps	/* IN	  - epsilon value */
)
{
int			i, j, t, r1, r2;
int			tcomp;
int			fsave;
int			pruned_total;
int			required_total;
int			old_pruned_total;
int			scan;
int			nverts;
int			kmasks;
int			nedges;
int			nmasks;
int			adj_edge;
int			del_comps_count;
int *			ep1;
int *			ep2;
int *			vp;
int *			vp1;
int *			vp2;
int *			comps_edge;
int *			made_req;
int			req_count;
bool			first_2edge;
bool			this_2edge;
struct dsuf		comps;
struct dsuf		del_comps;
struct pinfo		pinfo;
struct pset *		ltlist;
bitmap_t *		ltmask;
int *			lflist;
bitmap_t *		lfmask;

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

	pinfo.cip = cip;
	pinfo.eps = eps;
	pinfo.compat_mask = NEWA (nmasks, bitmap_t);
	memset (pinfo.compat_mask, 0, nmasks * sizeof (bitmap_t));

	/* The following are only used for calling upper bound procedure */
	ltlist	  = NEW_PSET (nverts);
	ltmask	  = NEWA (kmasks, bitmap_t);
	lflist	  = NEWA (nedges, int);
	lfmask	  = NEWA (nmasks, bitmap_t);

	/* Initialize masks */
	for (i = 0; i < kmasks; i++) {
		ltmask [i] = 0;
	}
	for (i = 0; i < nmasks; i++) {
		lfmask [i] = 0;
	}

	/* Perform thorough upper bound test for every not-yet pruned FST */
	pruned_total = 0;
	for (i = 0; i < nedges; i++) {

		if (BITON (cip -> initial_edge_mask, i)) {
			if (NOT passes_upper_bound_tests (&pinfo, bsd, i, -1,
							  ltlist, ltmask, lflist, lfmask)) {
				CLRBIT (cip -> initial_edge_mask, i);
				pruned_total++;
			}
			else {
				SETBIT (pinfo.compat_mask, i);
			}
		}
		else {
			pruned_total++;
		}
	}

	/* Compute basic compatibility */
	/* CHANGE: USE COMP INFO THAT IS ALREADY THERE? */
	compute_incompatibility_info (&pinfo, bsd);

	/* Compute pruning information */
	compute_pruning_info(cip, &pinfo);

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

	comps_edge = NEWA (cip -> num_verts, int);
	made_req   = NEWA (nedges, int);

	/* Perform actual pruning */
	if (Print_Detailed_Timings) {
		fprintf(stderr, "- scan 0 finished. %6d FSTs pruned\n", pruned_total);
	}
	required_total = 0;
	for (scan = 1; scan < nedges; scan++) {

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

				/* Set up mask of compatible FSTs */
				ep1 = cip -> inc_edges [i];
				ep2 = cip -> inc_edges [i + 1];
				while (ep1 < ep2) {
					j = *ep1++;
					CLRBIT (pinfo.compat_mask, j);
				}

				/* Test if FST can be pruned */
				if (prune_this_fst(cip, &pinfo, i)) {
					CLRBIT (cip -> initial_edge_mask, i);
					CLRBIT (pinfo.compat_mask, i);
					pruned_total++;
				}

				/* Reset mask */
				ep1 = cip -> inc_edges [i];
				while (ep1 < ep2) {
					j = *ep1++;
					if (BITON (cip -> initial_edge_mask, j))
						SETBIT (pinfo.compat_mask, j);
				}
			}
		}

		/* Test if any connected component (initially one
		   for each terminal) only has one adjacent FST */

	try_again:
		req_count = 0;
		for (t = 0; t < cip -> num_verts; t++) {
			comps_edge [t] = -1;
		}

		/* First find the number of adjacent FSTs... */
		for (t = 0; t < cip -> num_verts; t++) {
			tcomp = dsuf_find (&comps, t); /* t's component */
			if (comps_edge [tcomp] EQ -2) continue;
			ep1 = cip -> term_trees [t];
			ep2 = cip -> term_trees [t + 1];
			while (ep1 < ep2) {
				i = *ep1++;
				if ((BITON (cip -> initial_edge_mask, i)) AND
				    (NOT BITON (cip -> required_edges, i))) {
					adj_edge = comps_edge [tcomp];
					if (adj_edge EQ -2) break;
					if (adj_edge EQ -1) {
						/* First FST */
						comps_edge [tcomp] = i;
					}
					else {
						/* Second (or more) FST */
						first_2edge = (cip -> edge_size [adj_edge] EQ 2);
						this_2edge  = (cip -> edge_size [i] EQ 2);

						/* If all adjacent edges have been 2-edges, */
						/* then pick the shortest (if equal then take edge */
						/* with minimum index) */

						if (first_2edge AND this_2edge) {
							if ((cip -> cost [i] < cip -> cost [adj_edge]) OR
							    ((cip -> cost [i] EQ cip -> cost [adj_edge]) AND
							     (i < adj_edge))) {
								comps_edge [tcomp] = i;
							}
						}
						else {
							comps_edge [tcomp] = -2; break;
						}
					}
				}
			}
		}

		/* ... then check the counts */
		for (t = 0; t < cip -> num_verts; t++) {
			fsave = comps_edge [t];
			if ((fsave >= 0) AND