rfst.c

生成直角Steiner树的程序包

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

	File:	rfst.c
	Rev:	b-1
	Date:	12/22/2000

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

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

	The main routine for the rectilinear FST generator.  It
	reads a point set from standard input, generates the FSTs,
	and outputs them to standard output.

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

	Modification Log:

	a-1:	08/31/98	warme
		: Created.  Derived from Martin's prototype.
	b-1:	12/22/2000	martinz
		: Added -k option to denote the maximum size of FSTs
		: that should be generated.
		: FSTs for which corner-flipped versions split
		: into two or more FSTs are removed.

************************************************************************/

#include "bsd.h"
#include "emptyr.h"
#include "p1io.h"
#include "rfst.h"
#include "steiner.h"


/*
 * Global Routines
 */

int			main (int, char **);


/*
 * Local Constants
 */

#define EMPTY_DIAMOND_PROPERTY	1
#define UB_SHORTLEG		1
#define NOSPLIT_CORNER_FLIPPED	1
#define KAHNG_ROBINS_HEURISTIC	0
#define DO_STATISTICS		0


/*
 * Local Types
 */

struct ibuf {
	struct ibuf *	next;
	int		count;
	int		buf [1];
};


/*
 * Local Routines
 */

static void		add_zero_length_fsts (struct rinfo *, int, int **);
static void		build_fst_list (struct rinfo *);
static void		build_rfst_graph (struct rinfo *,
					  int,
					  int,
					  int,
					  struct pset *,
					  struct pset *,
					  struct edge *,
					  int);
static void		compute_rfsts_for_unique_terminals (struct rinfo *);
static void		compute_successors (struct rinfo *);
static void		compute_ub0 (struct rinfo *);
static void		compute_ub1 (struct rinfo *);
static void		compute_zt (struct rinfo *);
static void		convert_delta_cpu_time (char *);
static void		decode_params (int, char **);
static dist_t		delta_x_func (struct point *, struct point *);
static dist_t		delta_y_func (struct point *, struct point *);
static bool		diamond_empty (struct rinfo *,
				       struct point *,
				       struct point *,
				       int,
				       int);
static int		generate_duplicate_terminal_groups (struct pset *,
							    int *,
							    int ***);
static cpu_time_t	get_delta_cpu_time (void);
static void		grow_RFST (struct rinfo *	rip,
				   int			size,
				   dist_t		length,
				   int			dir,
				   dist_t		ub_length,
				   dist_t *		ub_shortleg,
				   int			longindex);
static int *		heapsort_x (struct pset *);
static int *		heapsort_y (struct pset *);
static int		lrindex_dir_0 (struct point *, struct point *);
static int		lrindex_dir_1 (struct point *, struct point *);
static int		lrindex_dir_2 (struct point *, struct point *);
static int		lrindex_dir_3 (struct point *, struct point *);
static struct full_set ** put_trees_in_array (struct full_set *, int *);
static struct pset *	remove_duplicates (struct pset * pts,
					   int		 ndg,
					   int **	 list,
					   int **	 fwd_map,
					   int **	 rev_map);
static void		renumber_terminals (struct rinfo *,
					    struct pset *,
					    int *);
static dist_t		test_and_save_fst (struct rinfo *,
					   int,
					   dist_t,
					   int,
					   int);
static void		usage (void);


/*
 * Local Variables
 */

static char *		description;
static char *		me;
static int		MaxFSTSize = 0;
static int		output_version = CURRENT_P1IO_VERSION;
static bool		Print_Detailed_Timings = FALSE;
static cpu_time_t	T0;
static cpu_time_t	Tn;

/*
 * Data tables for implementing the DSTDIR and DSTDIRP macros without
 * conditional branches.  This should be faster on machines for which
 * pipeline flushes are an issue.
 * We could really use the pointer-to-member feature of C++ here!
 */

#define X	offsetof(struct point, x)
#define Y	offsetof(struct point, y)

const int8u	dstdir_offsets [] = {
	X, Y, X, Y, X,
};

#undef X
#undef Y

/* Tables for the array-of-function-pointers implementation. */

typedef dist_t	(* dstdir_funcp) (struct point *, struct point *);

const dstdir_funcp	dstdir_funcs [] = {
	delta_x_func,
	delta_y_func,
	delta_x_func,
	delta_y_func,
	delta_x_func,
};


/* Tables for the LRINDEX implementation using an array of functions. */

const lrindex_funcp	lrindex_funcs [] = {
	lrindex_dir_0,
	lrindex_dir_1,
	lrindex_dir_2,
	lrindex_dir_3,
};

/*
 * The main routine for the "rfst" program.  It reads a point set
 * from standard input, generates all of the rectilinear FSTs,
 * and outputs them to standard output in our special "phase 1 I/O"
 * format.
 */

	int
main (

int		argc,
char **		argv
)
{
int			i;
int			j;
int			k;
int			ndg;
int			ntrees;
int			count;
int			fpsave;
int **			dup_grps;
int *			fwd_map;
int *			rev_map;
int *			ip1;
struct full_set *	fsp;
struct pset *		terms;
struct rinfo		rinfo;
struct cinfo		cinfo;
struct pset *		pts;
struct pset *		pts2;
cpu_time_t		Trenum;
struct scale_info	scale;
char			buf1 [32];

	fpsave = set_floating_point_double_precision ();

	setbuf (stdout, NULL);

	decode_params (argc, argv);

	pts = get_points (stdin, &scale);

	init_output_conversion (pts, RECTILINEAR, &scale);

	T0 = get_cpu_time ();
	Tn = T0;

	rinfo.x_order = heapsort_x (pts);

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

	/* Find all duplicate terminals in the input. */
	ndg = generate_duplicate_terminal_groups (pts,
						  rinfo.x_order,
						  &dup_grps);

	rinfo.num_term_masks = BMAP_ELTS (pts -> n);

	/* Remove all but the first of each duplicate terminal. */
	/* Compute forward and reverse maps to renumber the terminals. */
	pts2 = remove_duplicates (pts, ndg, dup_grps, &fwd_map, &rev_map);

	/* Renumber the x_order list -- instead of re-sorting pts2. */
	j = 0;
	for (i = 0; i < pts -> n; i++) {
		k = rinfo.x_order [i];
		if ((k < 0) OR (pts -> n < k)) {
			fatal ("main: Bug 1.");
		}
		k = fwd_map [k];
		if (k < 0) continue;
		rinfo.x_order [j++] = k;
	}

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

	rinfo.y_order = heapsort_y (pts2);

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

	/* From now on, we work only with the reduced terminal set, and */
	/* we assume that all terminals are unique. */

	rinfo.pts = pts2;

	compute_rfsts_for_unique_terminals (&rinfo);

	/* Now put the terminal numbers back the way they were, */
	/* renumber the terminals within each RFST, etc. */

	renumber_terminals (&rinfo, pts, rev_map);

	/* Link the FSTs together into one long list, and number them. */
	build_fst_list (&rinfo);

	/* Add one FST for each duplicate terminal that was removed. */
	if (ndg > 0) {
		add_zero_length_fsts (&rinfo, ndg, dup_grps);
	}

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

	if (Print_Detailed_Timings) {
		convert_cpu_time (Trenum, buf1);
		fprintf (stderr, "Renumber Terminals:     %s\n", buf1);
		convert_cpu_time (Tn - T0, buf1);
		fprintf (stderr, "Total:                  %s\n", buf1);
	}

	/* Zero out the compatibility info. */
	memset (&cinfo, 0, sizeof (cinfo));

	cinfo.num_edges			= rinfo.ntrees;
	cinfo.num_verts			= rinfo.pts -> n;
	cinfo.num_vert_masks		= BMAP_ELTS (cinfo.num_verts);
	cinfo.num_edge_masks		= BMAP_ELTS (cinfo.num_edges);
	cinfo.edge			= NEWA (rinfo.ntrees + 1, int *);
	cinfo.edge_size			= NEWA (rinfo.ntrees, int);
	cinfo.cost			= NEWA (rinfo.ntrees, dist_t);
	cinfo.tflag			= NEWA (cinfo.num_verts, bool);
	cinfo.metric			= RECTILINEAR;
	cinfo.scale			= scale;
	cinfo.integrality_delta		= 1.0;
	cinfo.pts			= rinfo.pts;
	cinfo.full_trees		= put_trees_in_array (rinfo.full_sets,
							      &ntrees);
	cinfo.mst_length		= rinfo.mst_length;
	cinfo.description		= gst_strdup (description);
	cinfo.p1time			= Tn - T0;

	for (i = 0; i < cinfo.num_verts; i++) {
		cinfo.tflag [i] = TRUE;
	}

	count = 0;
	for (i = 0; i < rinfo.ntrees; i++) {
		fsp = cinfo.full_trees [i];
		k = fsp -> terminals -> n;
		cinfo.edge_size [i]	= k;
		cinfo.cost [i]		= fsp -> tree_len;
		count += k;
	}
	ip1 = NEWA (count, int);
	for (i = 0; i < rinfo.ntrees; i++) {
		cinfo.edge [i] = ip1;
		terms = cinfo.full_trees [i] -> terminals;
		k = terms -> n;
		for (j = 0; j < k; j++) {
			*ip1++ = terms -> a [j].pnum;
		}
	}
	cinfo.edge [i] = ip1;

	/* Print all of the data we have generated... */
	print_phase_1_data (&cinfo, output_version);

	/* Clean up. */

	free ((char *) (rinfo.y_order));
	free ((char *) fwd_map);
	free ((char *) rev_map);
	free ((char *) (rinfo.x_order));

	free ((char *) pts2);
#if 0
	free ((char *) pts);
#endif

	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;
int		v;

	--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 'k':
				/* Max number of terminals in FSTs */
				if (*ap EQ '\0') {
					if (argc <= 0) {
						usage ();
					}
					ap = *argv++;
					--argc;
				}
				MaxFSTSize = atoi (ap);
				ap = "";
				break;

			case 't':
				Print_Detailed_Timings = TRUE;
				break;

			case 'v':
				if (*ap EQ '\0') {
					if (argc <= 0) {
						usage ();
					}
					ap = *argv++;
					--argc;
				}
				if ((*ap < '0') OR (*ap > '9')) {
					usage ();
				}
				v = atoi (ap);
				if ((v < 0) OR (v > CURRENT_P1IO_VERSION)) {
					fprintf (stderr,
						 "%s: Bad version `%s'."
						 "  Valid versions range"
						 " from 0 to %d.\n",
						 me,
						 ap,
						 CURRENT_P1IO_VERSION);
					usage ();
				}
				output_version = v;
				ap = "";
				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-k K\tOnly generate FSTs spanning up to K terminals.",
	"\t-t\tPrint detailed timings on stderr.",
	"\t-v N\tGenerates version N output data format.",
	"",
	NULL
};

	static
	void
usage (void)

{
char **		pp;
char *		p;

	(void) fprintf (stderr,
			"\nUsage: %s [-t] [-d description] [-k K] [-v N]"
			" <points-file\n",
			me);
	pp = &arg_doc [0];
	while ((p = *pp++) NE NULL) {
		(void) fprintf (stderr, "%s\n", p);
	}
	exit (1);
}

/*
 * Use the heapsort algorithm to sort the given terminals in increasing
 * order by the following keys:
 *
 *	1.	X coordinate
 *	2.	Y coordinate
 *	3.	index (i.e., position within input data)
 *
 * Of course, we do not move the points, but rather permute an array
 * of indexes into the points.
 */

	static
	int *
heapsort_x (

struct pset *		pts		/* IN - the terminals to sort */
)
{
int			i, i1, i2, j, k, n;
struct point *		p1;
struct point *		p2;
int *			index;

	n = pts -> n;

	index = NEWA (n, int);
	for (i = 0; i < n; i++) {
		index [i] = i;
	}

	/* Construct the heap via sift-downs, in O(n) time. */
	for (k = n >> 1; k >= 0; k--) {
		j = k;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is greater. */
				i1 = index [i];
				i2 = index [i + 1];
				p1 = &(pts -> a [i1]);
				p2 = &(pts -> a [i2]);
				if ((p2 -> x > p1 -> x) OR
				    ((p2 -> x EQ p1 -> x) AND
				     ((p2 -> y > p1 -> y) OR
				      ((p2 -> y EQ p1 -> y) AND
				       (i2 > i1))))) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}
			i1 = index [j];
			i2 = index [i];
			p1 = &(pts -> a [i1]);
			p2 = &(pts -> a [i2]);
			if ((p1 -> x > p2 -> x) OR
			    ((p1 -> x EQ p2 -> x) AND
			     ((p1 -> y > p2 -> y) OR
			      ((p1 -> y EQ p2 -> y) AND
			       (i1 > i2))))) {
				/* Greatest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			index [j] = i2;
			index [i] = i1;
			j = i;
		}
	}

	/* Now do actual sorting.  Exchange first/last and sift down. */
	while (n > 1) {
		/* Largest is at index [0], swap with index [n-1],	*/
		/* thereby putting it into final position.		*/
		--n;
		i = index [0];
		index [0] = index [n];
		index [n] = i;

		/* Now restore the heap by sifting index [0] down. */
		j = 0;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is greater. */
				i1 = index [i];
				i2 = index [i + 1];
				p1 = &(pts -> a [i1]);
				p2 = &(pts -> a [i2]);
				if ((p2 -> x > p1 -> x) OR
				    ((p2 -> x EQ p1 -> x) AND
				     ((p2 -> y > p1 -> y) OR
				      ((p2 -> y EQ p1 -> y) AND
				       (i2 > i1))))) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}
			i1 = index [j];
			i2 = index [i];
			p1 = &(pts -> a [i1]);
			p2 = &(pts -> a [i2]);
			if ((p1 -> x > p2 -> x) OR
			    ((p1 -> x EQ p2 -> x) AND
			     ((p1 -> y > p2 -> y) OR
			      ((p1 -> y EQ p2 -> y) AND
			       (i1 > i2))))) {
				/* Greatest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			index [j] = i2;
			index [i] = i1;
			j = i;
		}
	}

	return (index);
}

/*
 * Use the heapsort algorithm to sort the given terminals in increasing
 * order by the following keys:
 *
 *	1.	Y coordinate
 *	2.	X coordinate
 *	3.	index (i.e., position within input data)
 *
 * Of course, we do not move the points, but rather permute an array
 * of indexes into the points.
 */

	static
	int *
heapsort_y (

struct pset *		pts		/* IN - the terminals to sort */
)
{
int			i, i1, i2, j, k, n;
struct point *		p1;
struct point *		p2;
int *			index;

	n = pts -> n;

	index = NEWA (n, int);
	for (i = 0; i < n; i++) {
		index [i] = i;
	}

	/* Construct the heap via sift-downs, in O(n) time. */
	for (k = n >> 1; k >= 0; k--) {
		j = k;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is greater. */
				i1 = index [i];
				i2 = index [i + 1];
				p1 = &(pts -> a [i1]);
				p2 = &(pts -> a [i2]);
				if ((p2 -> y > p1 -> y) OR
				    ((p2 -> y EQ p1 -> y) AND
				     ((p2 -> x > p1 -> x) OR
				      ((p2 -> x EQ p1 -> x) AND
				       (i2 > i1))))) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}