efst.c

生成直角Steiner树的程序包

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

	File:	efst.c
	Rev:	b-2
	Date:	02/28/2001

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

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

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

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

	Modification Log:

	a-1:	12/11/98	martinz
		: Created.  Derived from Pawel and Martin's C++ program
			    using Warme's infrastructure.
	b-1:	11/22/2000	martinz
		: Dynamic allocation of eq-point array using doubling.
		:  (-e option is now the INITIAL number of eq-points
		:  per terminal.)
		: Memory for rectangle structure freed when done.
		: Improved numerical robustness.
		: Translate points such that their mean is at the
		:  origin in order to compute eq-points with
		:  maximum precision.
		: Split off elementary geometric functions to efuncs.h.
		: Improved upper bound heuristic with BSD info.
		: Added new greedy heuristic.
		: Compute eq-points carefully by using displacements.
		: Compute distances between eq-points carefully by
		:  using displacements.
	b-2:	02/28/2001	warme
		: Use GMP, if available, to compute eq-points and
		:  EFST lengths precisely.

************************************************************************/

#include "bsd.h"
#include "config.h"
#include "efst.h"
#include "efuncs.h"
#include "egmp.h"
#include "p1io.h"
#include "steiner.h"


/*
 * Global Routines
 */

int			main (int, char **);

int			Multiple_Precision = 0;


/*
 * Local Routines
 */

static void		add_zero_length_fsts (struct einfo *, int, int **);
static void		build_fst_list (struct einfo *);
static void		build_efst_graph (struct einfo *,
					  struct point *,
					  struct eqp_t *,
					  struct point **,
					  struct edge **,
					  int *,
					  struct point *,
					  int *);
static void		compute_efsts_for_unique_terminals (struct einfo *);
static void		convert_delta_cpu_time (char *);
static void		decode_params (int, char **);
static int		generate_duplicate_terminal_groups (struct pset *,
							    int *,
							    int ***);
static cpu_time_t	get_delta_cpu_time (void);
static int *		heapsort_x (struct pset *);
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 einfo *,
					    struct pset *,
					    int *);
static dist_t		test_and_save_fst (struct einfo *,
					   struct eqp_t *,
					   struct eqp_t *);
static void		usage (void);


/*
 * Local Variables
 */

static char *		description;
static int		InitialEqPointsTerminal = 100;
static int		EpsilonFactor = 32;
static char *		me;
static int		MaxFSTSize = 0;
static int		output_version = CURRENT_P1IO_VERSION;
static bool		Use_Greedy_Heuristic = FALSE;
static bool		Print_Detailed_Timings = FALSE;
static cpu_time_t	T0;
static cpu_time_t	Tn;


/*
 * Local Macros
 */

#define UPDATE_PTR(p,old,new) ((new) + ((p) - (old)))

#define UPDATE_RECTANGLE_BOUNDS(p) \
	{ *minx = MIN(*minx, p.x); *maxx = MAX(*maxx, p.x); \
	  *miny = MIN(*miny, p.y); *maxy = MAX(*maxy, p.y); }

/*
 * The main routine for the "efst" 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 einfo		einfo;
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, EUCLIDEAN, &scale);

	T0 = get_cpu_time ();
	Tn = T0;

	einfo.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,
						  einfo.x_order,
						  &dup_grps);

	einfo.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 = einfo.x_order [i];
		if ((k < 0) OR (pts -> n < k)) {
			fatal ("main: Bug 1.");
		}
		k = fwd_map [k];
		if (k < 0) continue;
		einfo.x_order [j++] = k;
	}

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

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

	einfo.pts = pts2;

	compute_efsts_for_unique_terminals (&einfo);

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

	renumber_terminals (&einfo, pts, rev_map);

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

	/* Add one FST for each duplicate terminal that was removed. */
	if (ndg > 0) {
		add_zero_length_fsts (&einfo, 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			= einfo.ntrees;
	cinfo.num_verts			= einfo.pts -> n;
	cinfo.num_vert_masks		= BMAP_ELTS (cinfo.num_verts);
	cinfo.num_edge_masks		= BMAP_ELTS (cinfo.num_edges);
	cinfo.edge			= NEWA (einfo.ntrees + 1, int *);
	cinfo.edge_size			= NEWA (einfo.ntrees, int);
	cinfo.cost			= NEWA (einfo.ntrees, dist_t);
	cinfo.tflag			= NEWA (cinfo.num_verts, bool);
	cinfo.metric			= EUCLIDEAN;
	cinfo.scale			= scale;
	cinfo.integrality_delta		= 0;
	cinfo.pts			= einfo.pts;
	cinfo.full_trees		= put_trees_in_array (einfo.full_sets,
							      &ntrees);
	cinfo.mst_length		= einfo.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 < einfo.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 < einfo.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 *) pts2);
#if 0
	free ((char *) pts);
#endif

	free_phase_1_data (&cinfo);

	free ((char *) rev_map);
	free ((char *) fwd_map);
	if (dup_grps NE NULL) {
		if (dup_grps [0] NE NULL) {
			free ((char *) (dup_grps [0]));
		}
		free ((char *) dup_grps);
	}
	free ((char *) (einfo.x_order));

	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 'e':
				/* Number of eq-points per terminal */
				if (*ap EQ '\0') {
					if (argc <= 0) {
						usage ();
					}
					ap = *argv++;
					--argc;
				}
				InitialEqPointsTerminal = atoi (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 'g':
				Use_Greedy_Heuristic = TRUE;
				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;

#ifdef HAVE_GMP
			case 'm':
				if (*ap EQ '\0') {
					if (argc <= 0) {
						usage ();
					}
					ap = *argv++;
					--argc;
				}
				if ((*ap < '0') OR (*ap > '9')) {
					usage ();
				}
				Multiple_Precision = atoi (ap);
				ap = "";
				break;
#endif

			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-e K\tInitially allocate K eq-points per terminal",
	"\t\t(default: 100).",
	"\t-f E\tEpsilon multiplication factor for floating point number",
	"\t\tcomparisons (default: 32).",
	"\t-g\tUse greedy heuristic instead of Smith-Lee-Liebman",
	"\t\t(more time consuming but generates fewer eq-points).",
	"\t-k K\tOnly generate FSTs spanning up to K terminals.",
#ifdef HAVE_GMP
	"\t-m N\tUse multiple precision.  Larger N use it more.",
	"\t\t Default is N=0 which disables multiple precision.",
#endif
	"\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 [-gt]"
			" [-d description]"
			" [-e K] [-f E] [-k K]"
#ifdef HAVE_GMP
			" [-m N]"
#endif