sec_heur.c

生成直角Steiner树的程序包

			pts -> a [i].x = x / k;
			pts -> a [i].y = y / k;
		}
		else {
			pts -> a [i].x = 0;
			pts -> a [i].y = 0;
		}
		pts -> a [i].pnum	= i;
	}

	return (pts);
}

/*
 * This routine computes a Minimum Spanning Tree for the given FST
 * from the given congested component.
 */

	static
	int
full_set_mst (

struct comp *		comp,		/* IN - congested component */
int			edge,		/* IN - edge number in component */
struct pset *		pts,		/* IN - X,Y coords of vertices */
struct edge *		elist		/* OUT - list of MST edges */
)
{
int			n;
int			t1;
int			t2;
int			nedges;
int			num_mst_edges;
int			num_mst_verts;
struct point *		p1;
struct point *		p2;
struct edge *		ep;
int *			ip1;
int *			ip2;
int *			endp;
struct edge *		edges;

	ip1  = comp -> everts [edge];
	endp = comp -> everts [edge + 1];

	n = endp - ip1;

	n = n * (n - 1) / 2;

	edges = NEWA (n, struct edge);

	/* Generate list of all edges, and compute the length	*/
	/* of each.  Note that we list the endpoints of each	*/
	/* edge in increasing order so that edges may be easily	*/
	/* compared and sorted...				*/
	nedges = 0;
	ep = &edges [0];
	num_mst_verts = endp - ip1;
	while (ip1 < endp) {
		t1 = *ip1++;
		p1 = &(pts -> a [t1]);
		ip2 = ip1;
		while (ip2 < endp) {
			t2 = *ip2++;
			p2 = &(pts -> a [t2]);
#if 1
			/* Rectilinear distance. */
			ep -> len = RDIST (p1, p2);
#else
			/* Euclidean distance. */
			ep -> len = EDIST (p1, p2);
#endif
			if (t1 < t2) {
				ep -> p1 = t1;
				ep -> p2 = t2;
			}
			else {
				ep -> p1 = t2;
				ep -> p2 = t1;
			}
			++ep;
			++nedges;
		}
	}

	if (nedges > 0) {
		num_mst_edges = mst_edge_list (num_mst_verts,
					       nedges,
					       edges,
					       elist);
	}
	else {
		num_mst_edges = 0;
	}

	free ((char *) edges);

	return (num_mst_edges);
}

/*
 * This routine sorts the given list of edges (and corresponding edge
 * numbers) on three keys -- first and second endpoint vertex numbers,
 * and then edge number.  Because p1 < p2 for every edge, this is a
 * well-defined ordering of edges for the undirected graph.
 */

	static
	void
sort_edges (

int			n,	/* IN - number of edges to sort. */
struct edge *		edges,	/* IN/OUT - edges to sort. */
int *			e_num	/* IN/OUT - edge number for each edge. */
)
{
int		h;
struct edge	etmp;
struct edge *	p1;
struct edge *	p2;
struct edge *	p3;
struct edge *	p4;
struct edge *	endp;
int		itmp;
int *		ip1;
int *		ip2;
int *		ip3;

	endp = &edges [n];

	for (h = 1; h <= n; h = 3*h+1) {
	}

	do {
		h = h / 3;
		p4 = &edges [h];
		p1 = p4;
		ip1 = &e_num [h];
		while (p1 < endp) {
			etmp = *p1;
			itmp = *ip1;
			p2 = p1;
			ip2 = ip1;
			while (TRUE) {
				p3 = (p2 - h);
				ip3 = (ip2 - h);
				if (p3 -> p1 < etmp.p1) break;
				if (p3 -> p1 EQ etmp.p1) {
					if (p3 -> p2 < etmp.p2) break;
					if (p3 -> p2 EQ etmp.p2) {
						if (*ip3 <= itmp) {
							break;
						}
					}
				}
				/* *p3 > *p2 */
				*p2 = *p3;
				p2 = p3;
				*ip2 = *ip3;
				ip2 = ip3;
				if (p2 < p4) break;
			}
			*p2 = etmp;
			*ip2 = itmp;
			++p1;
			++ip1;
		}
	} while (h > 1);
}

/*
 * This routine heuristically identifies violated Subtour Elimination
 * Constraints by finding max-flows (i.e. min-cuts) in the previously
 * constructed directed graph.  The given LP solution is used to set
 * the arc capacities before doing the max-flow.
 *
 * A list of violated SEC constraints is returned.  This may be NULL,
 * in which case no violations were detected.
 */

	static
	struct constraint *
heuristically_find_violated_secs (

struct comp *		comp,		/* IN - congested component */
bitmap_t *		edge_mask,	/* IN - valid edges */
struct cinfo *		cip,		/* IN - compatibility info */
struct sec_heur_info *	hsecp,		/* IN - directed flow graph */
double *		x,		/* IN - LP solution to separate */
struct constraint *	clist		/* IN - existing constraints */
)
{
int			i;
int			j;
int			k;
int			n;
int			kmasks;
int			num_nodes;
int			num_arcs;
int			src_arc1;
int			dst_arc1;
int *			ip1;
int *			ip2;
double *		flow;
double *		c;
double			save;
bitmap_t *		stour;

	num_nodes	= hsecp -> prob.num_nodes;
	num_arcs	= hsecp -> prob.num_arcs;
	c		= hsecp -> prob.capacity;

	flow		= hsecp -> soln.flow;

	src_arc1	= hsecp -> src_arc1;
	dst_arc1	= hsecp -> dst_arc1;


	/* First we must set the arc capacities from the LP solution. */
	set_arc_capacities (comp, hsecp);

	/* Run one max-flow problem for each problem vertex... */

	n = num_nodes - 2;

	kmasks = cip -> num_vert_masks;

	stour  = NEWA (kmasks, bitmap_t);

	for (i = 0; i < n; i++) {
		/* Modify problem to find worst SEC involving	*/
		/* vertex i.  Save capacity (for restore).	*/
		save = c [src_arc1 + i];
		c [src_arc1 + i] = INFINITE_FLOW;

		compute_max_flow (&(hsecp -> prob),
				  &(hsecp -> temp),
				  &(hsecp -> soln));

		/* Construct bitmask of "real" vertices... */
		for (j = 0; j < kmasks; j++) {
			stour [j] = 0;
		}
		for (j = 0; j < n; j++) {
			if (NOT BITON (hsecp -> soln.cut, j)) continue;
			ip1 = comp -> rverts [j];
			ip2 = comp -> rverts [j + 1];
			while (ip1 < ip2) {
				k = *ip1++;
				SETBIT (stour, k);
			}
		}

		/* Generate constraint if a violation is present. */
		clist = check_subtour (stour,
				       clist,
				       x,
				       edge_mask,
				       cip);

		/* Restore the modified capacity... */
		c [src_arc1 + i] = save;

		/* Force node i to the far side of the cut... */
		c [dst_arc1 + i] = INFINITE_FLOW;
	}

	free ((char *) stour);

	return (clist);
}

/*
 * This routine sets the all of the arc capacities from the LP solution.
 * There are 3 kinds of arcs to handle:  the "normal" arcs comprised of
 * flow taken from 1 or more edges, "source" arcs, and "sink" arcs.
 */

	static
	void
set_arc_capacities (

struct comp *		comp,		/* IN - congested component */
struct sec_heur_info *	hsecp		/* IN - directed flow graph */
)
{
int			i;
int			j;
int			num_nodes;
int			nverts;
int			src_arc1;
int			dst_arc1;
double *		c;
int *			ip1;
int *			ip2;
double			sum;
double *		x;

	num_nodes	= hsecp -> prob.num_nodes;
	nverts		= num_nodes - 2;
	c		= hsecp -> prob.capacity;

	src_arc1	= hsecp -> src_arc1;
	dst_arc1	= hsecp -> dst_arc1;

	x = comp -> x;

	/* All arcs before the first "source" arc are "normal" arcs. */
	for (i = 0; i < src_arc1; i++) {
		ip1 = hsecp -> edge_lists [i];
		ip2 = hsecp -> edge_lists [i + 1];
		sum = 0.0;
		while (ip1 < ip2) {
			/* Total weight of all edges that contribute */
			/* to this arc... */
			j = *ip1++;
			sum += x [j];
		}
		c [i] = 0.5 * sum;
	}

	/* Now initialize all "source" and "sink" arcs... */
	for (i = 0; i < nverts; i++) {
		ip1 = comp -> vedges [i];
		ip2 = comp -> vedges [i + 1];
		sum = 0.0;
		while (ip1 < ip2) {
			j = *ip1++;
			sum += x [j];
		}

		if (sum > 2.0) {
			c [src_arc1 + i] = 0.5 * sum - 1.0;
			c [dst_arc1 + i] = 0.0;
		}
		else {
			c [src_arc1 + i] = 0.0;
			c [dst_arc1 + i] = 1.0 - 0.5 * sum;
		}
	}
}

/*
 * This routine frees up the memory allocated by the given SEC max-flow
 * formulation.
 */

	static
	void
free_heuristic_SEC_formulation (

struct sec_heur_info *	hsecp		/* IN - SEC flow formulation */
)
{
	/* Free up the buffers used to hold flow solutions */
	/* and temporary data... */
	free_flow_temp_data (&(hsecp -> temp));
	free_flow_solution_data (&(hsecp -> soln));

	/* Free up the problem formulation... */
	free ((char *) (hsecp -> prob.out [0]));
	free ((char *) (hsecp -> prob.in [0]));
	free ((char *) (hsecp -> edge_lists [0]));

	free ((char *) (hsecp -> prob.out));
	free ((char *) (hsecp -> prob.in));
	free ((char *) (hsecp -> edge_lists));
	free ((char *) (hsecp -> prob.arc_src));
	free ((char *) (hsecp -> prob.arc_dst));
	free ((char *) (hsecp -> prob.capacity));
}

/*
 * This routine locates any cycles formed by hyperedges hyperedges having
 * weight 1.0 in the solution -- in other words, solid integer cycles.
 *
 * Each time we identify a connected component of solid integral edges,
 * we then look for "almost integer" cycles -- integral except for a single
 * fractional edge.  We do this by checking each fractional edge for two
 * or more vertices in common with the integral connected component.
 */

	struct constraint *
find_integer_cycles (

double *		x,		/* IN - LP solution to check */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct constraint *	cp,		/* IN - previous constraints */
struct cinfo *		cip		/* IN - compatibility info */
)
{
int			i;
int			j;
int			nverts;
int			nedges;
int			kmasks;
int			nmasks;
bitmap_t *		integral_edges;
bitmap_t *		edges_left;
bitmap_t *		cc_edges;
bitmap_t *		emark;
int *			stack;
int			num_frac;
int			num_cycles;
int *			frac_edge_nums;
struct constraint *	cp1;
struct constraint *	cp2;
bool			cycles_found;
bitmap_t *		vmark;
bitmap_t *		stour;
bitmap_t *		cc_verts;

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

	stack			= NEWA (2 * nverts, int);
	frac_edge_nums		= NEWA (nedges, int);
	integral_edges		= NEWA (nmasks, bitmap_t);
	edges_left		= NEWA (nmasks, bitmap_t);
	cc_edges		= NEWA (nmasks, bitmap_t);
	emark			= NEWA (nmasks, bitmap_t);
	vmark			= NEWA (kmasks, bitmap_t);
	stour			= NEWA (kmasks, bitmap_t);
	cc_verts		= NEWA (kmasks, bitmap_t);

	for (i = 0; i < nmasks; i++) {
		integral_edges [i] = 0;
	}

	num_frac = 0;
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		if (x [i] + FUZZ >= 1.0) {
			/* Edge is present with weight 1.0 */
			SETBIT (integral_edges, i);
		}
		else if (x [i] > FUZZ) {
			/* Edge is fractionally present. */
			frac_edge_nums [num_frac++] = i;
		}
	}

	for (i = 0; i < nmasks; i++) {
		edges_left [i] = integral_edges [i];
		emark [i] = 0;
	}

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

	num_cycles = 0;
	for (i = 0; i < nedges; i++) {
		/* Find an edge that has not yet been traversed. */
		if (NOT BITON (edges_left, i)) continue;

		/* Prepare a place to hold the vertices and edges	*/
		/* of this integral connected-component...		*/
		for (j = 0; j < kmasks; j++) {
			cc_verts [j] = 0;
		}
		for (j = 0; j < nmasks; j++) {
			cc_edges [j] = 0;
		}

		/* Walk the component containing edge i, looking	*/
		/* for cycles...					*/
		cp1 = NULL;
		cp2 = cwalk (i,
			     -1,
			     integral_edges,
			     edges_left,
			     emark,
			     vmark,
			     stour,
			     cc_verts,
			     cc_edges,
			     stack,
			     stack,
			     &num_cycles,
			     cp1,
			     cip);

#if 0
		print_mask (" % Int CC verts:", cc_verts, nverts);
		print_mask (" % Int CC edges:", cc_edges, nedges);
#endif

		/* Add any cycles we found onto the main list... */
		cycles_found = FALSE;
		while (cp2 NE NULL) {
			cp1 = cp2 -> next;
			cp2 -> next = cp;
			cp = cp2;
			cp2 = cp1;
			cycles_found = TRUE;
		}

		if (num_cycles >= CYCLE_LIMIT) break;
		if (cycles_found) continue;

		/* The current connected-component does not contain any	*/
		/* solid integer cycles.  But now that we know all of	*/
		/* the vertices of this CC we can check for any		*/
		/* fractional edges that have at least two vertices in	*/
		/* common with the CC.  These represent "almost		*/
		/* integer" cycles...  The cycle is unique since the	*/
		/* integer CC is acyclic.				*/
		cp2 = find_almost_integer_cycles (num_frac,
						  frac_edge_nums,
						  cc_verts,
						  cc_edges,
						  stour,
						  cip);

		/* Add any "almost integer" cycles onto main list...	*/
		while (cp2 NE NULL) {
			cp1 = cp2 -> next;
			cp2 -> next = cp;
			cp = cp2;
			cp2 = cp1;
		}
	}

	free ((char *) cc_verts);
	free ((char *) stour);
	free ((char *) vmark);
	free ((char *) emark);
	free ((char *) cc_edges);
	free ((char *) edges_left);
	free ((char *) integral_edges);
	free ((char *) frac_edge_nums);
	free ((char *) stack);

	return (cp);
}

/*
 * This routine recursively walks connected components of hyperedges
 * looking for cycles.  It generates a Subtour Elimination Constraint
 * for every cycle found.
 */

	static
	struct constraint *
cwalk (

int			e,		/* IN - edge number */
int			vertex,		/* IN - vertex by which we entered */
					/*      this hyperedge (or -1) */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
bitmap_t *		edges_left,	/* IN - unvisited hyperedges */