p1read.c

生成直角Steiner树的程序包

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

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

	Copyright (c) 1993, 2001 by David M. Warme

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

	Routines for reading the data that is output from phase 1.

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

	Modification Log:

	b-1:	01/12/97	warme
		: Split p1io.c into reading and writing parts.
	b-2:	02/28/2001	warme
		: Numerous changes for 3.1 release.
		: New input scaling stuff.
		: Fixed uninitialized incompat array bug.
		: Use common sort_ints routine.
		: Free *all* cinfo data structures, and be more
		:  careful about freeing term_trees.

************************************************************************/

#include "p1io.h"
#include "steiner.h"


/*
 * Global Routines
 */

int **		compute_basic_incompat (struct cinfo *);
void		free_phase_1_data (struct cinfo *);
void		init_term_trees (struct cinfo *);
void		read_phase_1_data (struct cinfo *);


/*
 * External References
 */

	/* none */


/*
 * Local Routines
 */

static void		double_to_hex (double, char *);
static int		get_d (void);
static double		get_dec_double (void);
static double		get_hex_double (void);
static void		get_line (char *, int);
static int		get_version (void);
static bitmap_t		get_x (void);
static double		hex_to_double (char *);
static int		hexdig (char);
static void		init_inc_edges (struct cinfo *, int **, int *);
static void		read_duplicate_terminal_groups (struct cinfo *, int);
static void		read_version_0 (struct cinfo *, int);
static void		read_version_2 (struct cinfo *, int);
static void		remove_duplicates (int, int **, bitmap_t *);
static void		skip (void);
static void		verify_symmetric (int **, int);

/*
 * This routine reads in the data as printed by the print-phase-1-data
 * routine.
 */

	void
read_phase_1_data (

struct cinfo *	cip		/* OUT - compatibility info. */
)
{
int		version;

	/* Zero everything initially... */
	(void) memset (cip, 0, sizeof (*cip));

	version = get_version ();

	switch (version) {
	case P1IO_VERSION_0:
		read_version_0 (cip, version);
		break;

	case P1IO_VERSION_2:
	case P1IO_VERSION_3:
		read_version_2 (cip, version);
		break;

	default:
		fatal ("read_phase_1_data: Bug 1.");
		break;
	}

	/* Set up the output conversion stuff. */
	init_output_conversion (cip -> pts, cip -> metric, &(cip -> scale));
}

/*
 * This routine reads in the extended OR-library format.
 */

	static
	void
read_version_0 (

struct cinfo *	cip,		/* OUT - compatibility info. */
int		version		/* IN - version number. */
)
{
int			i;
int			j;
int			k;
char			c;
int			nverts;
int			nedges;
int			kmasks;
int			nmasks;
int			scale_factor;
int			total_edge_cardinality;
int			nt;
char *			s;
struct numlist *	line;
struct numlist *	p;
struct numlist **	hookp;
struct numlist **	prev_hookp;
struct numlist *	costs;
struct numlist **	cost_hookp;
int *			ip1;
int *			ip2;

	nverts = get_d ();
	nedges = get_d ();

	kmasks = BMAP_ELTS (nverts);
	nmasks = BMAP_ELTS (nedges);

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

	cip -> edge		= NEWA (nedges + 1, int *);
	cip -> edge_size	= NEWA (nedges, int);
	cip -> cost		= NEWA (nedges, dist_t);
	cip -> tflag		= NEWA (nverts, bool);
	cip -> metric		= PURE_GRAPH;

	cip -> initial_vert_mask = NEWA (kmasks, bitmap_t);
	cip -> initial_edge_mask = NEWA (nmasks, bitmap_t);
	cip -> required_edges	 = NEWA (nmasks, bitmap_t);

	memset (cip -> initial_vert_mask, 0, kmasks * sizeof (bitmap_t));
	memset (cip -> initial_edge_mask, 0, nmasks * sizeof (bitmap_t));
	memset (cip -> required_edges,	  0, nmasks * sizeof (bitmap_t));

	for (i = 0; i < nverts; i++) {
		SETBIT (cip -> initial_vert_mask, i);
	}
	for (i = 0; i < nedges; i++) {
		SETBIT (cip -> initial_edge_mask, i);
	}

	costs = NULL;
	cost_hookp = &costs;

	total_edge_cardinality = 0;

	for (i = 0; i < nedges; i++) {
		line = parse_line_of_numbers (stdin);
		if (line EQ NULL) {
			fprintf (stderr, "Unexpected EOF!\n");
			exit (1);
		}

		/* Count the numbers and detach the last one */
		j = 0;
		prev_hookp = NULL;
		hookp = &line;
		for (;;) {
			p = *hookp;
			if (p EQ NULL) break;
			prev_hookp = hookp;
			hookp = &(p -> next);
			++j;
		}
		if (j < 3) {
			fprintf (stderr,
				 "Hyperedge with less than 2 vertices!\n");
			exit (1);
		}

		/* Transfer last number onto cost list. */
		p = *prev_hookp;
		*prev_hookp = NULL;
		--j;
		cip -> edge_size [i] = j;
		total_edge_cardinality += j;

		*cost_hookp = p;
		cost_hookp = &(p -> next);

		ip1 = NEWA (j, int);
		cip -> edge [i] = ip1;

		/* Verify remaining numbers are integers and convert. */
		j = 0;
		for (p = line; p NE NULL; p = p -> next) {
			if (p -> expon < 0) goto nonint;
			if (p -> expon > 0) {
				do {
					p -> mantissa *= 10.0;
				} while (--(p -> expon) > 0);
			}
			if ((p -> mantissa < 1) OR
			    (nverts < p -> mantissa) OR
			    (p -> mantissa NE floor (p -> mantissa))) {
				fprintf (stderr,
					 "Vertex number %g out of range.\n",
					 p -> mantissa);
				exit (1);
			}
			*ip1++ = p -> mantissa - 1;
		}
		if (cip -> edge [i] + cip -> edge_size [i] NE ip1) {
			fatal ("read_version_0: Bug 1.");
		}
	}

	/* Now convert the hyperedge costs. */
	scale_factor = compute_scaling_factor (costs);

	set_scale_info (&(cip -> scale), scale_factor);

	i = 0;
	while (costs NE NULL) {
		p = costs;
		costs = p -> next;
		p -> next = NULL;
		cip -> cost [i++] = read_numlist (p, &(cip -> scale));
	}
	if (i NE nedges) {
		fatal ("read_version_0: Bug 2.");
	}

	/* Read in the list of terminal vertices... */
	if (scanf (" %d", &nt) NE 1) {
		/* EOF -- assume all vertices are terminals. */
		for (i = 0; i < nverts; i++) {
			cip -> tflag [i] = TRUE;
		}
	}
	else {
		/* All are Steiner vertices unless listed! */
		for (i = 0; i < nverts; i++) {
			cip -> tflag [i] = FALSE;
		}
		for (i = 0; i < nt; i++) {
			j = get_d ();
			if ((j < 1) OR (nverts < j)) {
				fatal ("read_version_0: Bug 3.");
			}
			cip -> tflag [j - 1] = TRUE;
		}
	}

	/* Copy hyperedges into contiguous form. */
	ip1 = NEWA (total_edge_cardinality, int);
	for (i = 0; i < nedges; i++) {
		ip2 = cip -> edge [i];
		cip -> edge [i] = ip1;
		k = cip -> edge_size [i];
		for (j = 0; j < k; j++) {
			*ip1++ = ip2 [j];
		}
		free ((char *) ip2);
	}
	cip -> edge [i] = ip1;

	init_term_trees (cip);
	cip -> inc_edges = compute_basic_incompat (cip);
	return;

	/* Error exit. */
nonint:
	fprintf (stderr, "Expected integer!\n");
	exit (1);
}

/*
 * This routine reads in the new data formats -- versions 2 and 3.
 */

	static
	void
read_version_2 (

struct cinfo *	cip,		/* OUT - compatibility info. */
int		version		/* IN - version number. */
)
{
int			i;
int			j;
int			k;
int			nverts;
int			nedges;
int			kmasks;
int			nmasks;
int			nt;
int			ns;
int			scale_factor;
int			fst_edge_count;
int			num_incompat;
int			ndg;
int			count;
int			total_edge_card;
struct pset *		terms;
struct pset *		steins;
struct point *		p1;
struct full_set *	fsp;
struct edge *		ep;
bitmap_t *		bp1;
bitmap_t *		bp2;
bitmap_t *		bp3;
bitmap_t		mask;
int			bitcount;
int *			ip1;
int *			ip2;
int *			icounts;
int **			incompat;
bool			geometric;
dist_t			len;
char			buf1 [128];

	skip ();	/* Skip rest of version line */

	get_line (buf1, sizeof (buf1));
	i = strlen (buf1);
	if (i > 79) {
		fprintf (stderr, "Description field too long.\n");
		exit (1);
	}
	cip -> description = new (i + 1);
	strcpy (cip -> description, buf1);

	cip -> metric = get_d ();
	if ((cip -> metric NE RECTILINEAR) AND
	    (cip -> metric NE EUCLIDEAN) AND
	    (cip -> metric NE PURE_GRAPH)) {
		fprintf (stderr, "Bad metric: %d\n", cip -> metric);
		exit (1);
	}
	geometric = (cip -> metric NE PURE_GRAPH);

	nverts = (int) get_d ();

	if (geometric) {
		(void) get_dec_double ();
		cip -> mst_length = get_hex_double ();
	}
	else {
		cip -> mst_length = 0;
	}

	ndg = 0;
	if ((version <= P1IO_VERSION_2) AND geometric) {
		/* Version 2 has duplicate terminal groups... */
		ndg = (int) get_d ();
	}

	scale_factor = (int) get_d ();

	set_scale_info (&(cip -> scale), scale_factor);

	cip -> integrality_delta = 0;
	if (version >= P1IO_VERSION_3) {
		/* Version 3 has integrality delta... */
		(void) get_dec_double ();
		cip -> integrality_delta = get_hex_double ();
	}

	skip ();	/* skip rest of line before machine description */
	get_line (buf1, sizeof (buf1));	/* read machine description line */

	cip -> p1time = get_d ();

	nedges = (int) get_d ();

	kmasks = BMAP_ELTS (nverts);
	nmasks = BMAP_ELTS (nedges);

	cip -> initial_vert_mask = NEWA (kmasks, bitmap_t);
	for (i = 0; i < kmasks; i++) {
		cip -> initial_vert_mask [i] = 0;
	}
	for (i = 0; i < nverts; i++) {
		SETBIT (cip -> initial_vert_mask, i);
	}
	cip -> initial_edge_mask	= NEWA (nmasks, bitmap_t);
	cip -> required_edges		= NEWA (nmasks, bitmap_t);
	for (i = 0; i < nmasks; i++) {
		cip -> initial_edge_mask [i]	= 0;
		cip -> required_edges [i]	= 0;
	}

	cip -> num_edges	= nedges;
	cip -> num_verts	= nverts;
	cip -> num_vert_masks	= kmasks;
	cip -> num_edge_masks	= nmasks;
	cip -> edge		= NEWA (nedges + 1, int *);
	cip -> edge_size	= NEWA (nedges, int);
	cip -> cost		= NEWA (nedges, dist_t);
	cip -> tflag		= NEWA (nverts, bool);
	if (geometric) {
		cip -> pts	= NEW_PSET (nverts);
		ZERO_PSET (cip -> pts, nverts);
		cip -> pts -> n = nverts;
		cip -> full_trees = NEWA (nedges, struct full_set *);
	}

	incompat		= NEWA (nedges, int *);
	icounts			= NEWA (nedges, int);

	if (geometric) {
		/* Read in the terminals... */
		for (i = 0; i < nverts; i++) {
			p1 = &(cip -> pts -> a [i]);
			(void) get_dec_double ();
			(void) get_dec_double ();
			p1 -> x		 = get_hex_double ();
			p1 -> y		 = get_hex_double ();
			p1 -> pnum	 = i;
		}
	}

	if (version >= P1IO_VERSION_3) {
		/* Version 3 has terminal/Steiner flag for each vertex. */
		for (i = 0; i < nverts; i++) {
			cip -> tflag [i] = get_d ();
		}
	}
	else {
		/* Version 2: assume all vertices are terminals. */
		for (i = 0; i < nverts; i++) {
			cip -> tflag [i] = TRUE;
		}
	}

	if ((version <= P1IO_VERSION_2) AND geometric) {
		read_duplicate_terminal_groups (cip, ndg);
	}

	/* hyperedges... */
	total_edge_card = 0;
	for (i = 0; i < nedges; i++) {
		nt = get_d ();
		total_edge_card += nt;
		cip -> edge_size [i] = nt;
		ip1 = NEWA (nt, int);
		cip -> edge [i] = ip1;
		if (geometric) {
			fsp = NEW (struct full_set);
			(void) memset (fsp, 0, sizeof (*fsp));
			cip -> full_trees [i] = fsp;
			fsp -> tree_num = i;
			terms = NEW_PSET (nt);
			ZERO_PSET (terms, nt);
			terms -> n = nt;
		}
		for (j = 0; j < nt; j++) {
			k = get_d ();
			if ((k < 1) OR (k > nverts)) {
				fprintf (stderr,
					 "Terminal index out of range.\n");
				exit (1);
			}
			--k;
			*ip1++ = k;
			if (geometric) {
				terms -> a [j] = cip -> pts -> a [k];
			}
		}
		(void) get_dec_double ();
		len = get_hex_double ();
		cip -> cost [i] = len;

		if (geometric) {
			fsp -> tree_len = len;
			ns = get_d ();
			steins = NEW_PSET (ns);
			ZERO_PSET (steins, ns);
			steins -> n = ns;
			for (j = 0; j < ns; j++) {
				p1 = &(steins -> a [j]);
				(void) get_dec_double ();
				(void) get_dec_double ();
				p1 -> x = get_hex_double ();
				p1 -> y = get_hex_double ();
				p1 -> pnum = j;
			}
			fsp -> terminals = terms;
			fsp -> steiners = steins;
			fst_edge_count = get_d ();
			fsp -> nedges = fst_edge_count;
			ep = NEWA (fst_edge_count, struct edge);
			fsp -> edges = ep;
			for (j = 0; j < fst_edge_count; j++, ep++) {
				ep -> len = 0;	/* should be unused... */
				k = get_d ();
				if ((-ns <= k) AND (k < 0)) {
					k = nt - k - 1;
				}
				else if ((0 < k) AND (k <= nt)) {
					--k;
				}
				else {
					fprintf (stderr, "Invalid edge endpoint!\n");
					exit (1);
				}
				ep -> p1 = k;
				k = get_d ();
				if ((-ns <= k) AND (k < 0)) {
					k = nt - k - 1;
				}
				else if ((0 < k) AND (k <= nt)) {
					--k;
				}
				else {
					fprintf (stderr, "Invalid edge endpoint!\n");
					exit (1);
				}
				ep -> p2 = k;
			}
		}
		k = get_d ();	/* full set status... */
		switch (k) {
		case 0:
			break;

		case 1:
			SETBIT (cip -> initial_edge_mask, i);
			break;

		case 2:
			SETBIT (cip -> required_edges, i);
			SETBIT (cip -> initial_edge_mask, i);
			break;

		default:
			fprintf (stderr, "Invalid full set status: %d\n", k);
			exit (1);
		}

		num_incompat = get_d ();
		icounts [i] = num_incompat;
		ip1 = NULL;
		if (num_incompat > 0) {
			ip1 = NEWA (num_incompat, int);
			for (j = 0; j < num_incompat; j++) {
				k = get_d ();
				if ((k <= 0) OR (k > nedges)) {
					fprintf (stderr, "Bad incompatible index.\n");
					exit (1);
				}
				ip1 [j] = k - 1;
			}
		}
		incompat [i] = ip1;

		if (version <= P1IO_VERSION_2) {
			/* Version 2 has strongly compatible full sets... */
			k = get_d ();
			for (j = 0; j < k; j++) {
				(void) get_d ();
			}
		}
	}

	/* Copy hyperedges into contiguous form. */
	ip1 = NEWA (total_edge_card, int);
	for (i = 0; i < nedges; i++) {
		ip2 = cip -> edge [i];
		cip -> edge [i] = ip1;
		k = cip -> edge_size [i];
		for (j = 0; j < k; j++) {
			*ip1++ = ip2 [j];
		}
		free ((char *) ip2);
	}
	cip -> edge [i] = ip1;

	init_term_trees (cip);
	init_inc_edges (cip, incompat, icounts);

	free ((char *) icounts);
	free ((char *) incompat);
}

/*
 * This routine reads in the duplicate terminal groups.  We do not know
 * how big this is going to be.  Therefore we must plan for the worst
 * and reallocate when we are done...
 */

	static
	void
read_duplicate_terminal_groups (

struct cinfo *		cip,		/* IN/OUT - compatibility info */
int			ndg		/* IN - number of duplicate groups */
)
{
int		i;
int		j;
int		k;
int		n;
int		t;
int		nverts;
int		kmasks;
int *		ip;
int *		real_terms;
bitmap_t *	mark;
int *		index;
int *		terms;
int **		dup_grps;

	if (ndg <= 0) {
		return;				/* nothing to read */
	}

	nverts = cip -> num_verts;

	if (2 * ndg > nverts) {
		/* Too many duplicate terminal groups! */
		fatal ("read_duplicate_terminal_groups: Bug 1.");
	}

	kmasks = cip -> num_vert_masks;

	mark = NEWA (kmasks, bitmap_t);
	index = NEWA (ndg + 1, int);
	terms = NEWA (nverts, int);

	for (i = 0; i < kmasks; i++) {
		mark [i] = 0;
	}

	k = 0;
	for (i = 0; i < ndg; i++) {
		index [i] = k;
		n = (int) get_d ();
		for (j = 0; j < n; j++) {
			t = get_d ();
			if ((t < 1) OR (cip -> num_verts < t)) {
				fatal ("read_duplicate_terminal_groups: Bug 2.");
			}
			--t;
			if (BITON (mark, t)) {
				fatal ("read_duplicate_terminal_groups: Bug 3.");
			}
			SETBIT (mark, t);
			terms [k++] = t;
		}
	}
	index [i] = k;

	real_terms = NEWA (k, int);
	(void) memcpy ((char *) real_terms,
		       (char *) terms,
		       k * sizeof (int));

	dup_grps = NEWA (ndg + 1, int *);
	for (i = 0; i <= ndg; i++) {
		dup_grps [i] = real_terms + index [i];
	}

	free ((char *) dup_grps);
	free ((char *) real_terms);
	free ((char *) terms);
	free ((char *) index);
	free ((char *) mark);

	/* Remove duplicate terminals from the problem... */
	remove_duplicates (ndg, dup_grps, cip -> initial_vert_mask);
}

/*
 * This routine creates the "term_trees" array that is indexed by point
 * number and gives a list of all tree-numbers involving that point.
 */

	void
init_term_trees (

struct cinfo *	cip		/* IN - compatibility info */
)
{
int			i;
int			j;
int			n;
int			nverts;
int			nedges;
struct full_set *	fsp;
struct pset *		terms;
struct point *		p1;
int			total;
int *			ip;
int *			counts;
int **			ptrs;
int *			vp1;
int *			vp2;

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

	counts	= NEWA (nverts, int);

	total = 0;
	for (i = 0; i < nverts; i++) {
		counts [i] = 0;
	}

	for (i = 0; i < nedges; i++) {
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			++counts [*vp1++];
			++total;