sec_heur.c

生成直角Steiner树的程序包

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

	File:	sec_heur.c
	Rev:	b-1
	Date:	02/28/2001

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

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

	Heuristic separation procedure for SEC's (Subtour
	Elimination Constraints).  The heuristic is that we reduce
	the problem from weighted hypergraphs down to standard
	weighted graphs by computing a geometric MST for each
	hyperedge.

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

	Modification Log:

	a-1:	05/13/96	warme
		: Created.
	b-1:	02/28/2001	warme
		: Numerous changes for 3.1 release.
		: Completely rewrote subtour enumeration code to
		:  use fast incremental update of subtour constraint
		:  LHS and RHS values.

************************************************************************/

#include "bb.h"
#include "constrnt.h"
#include "flow.h"
#include "sec_comp.h"
#include "sec_heur.h"
#include "steiner.h"


/*
 * Global Routines
 */

struct constraint *	check_subtour (bitmap_t *,
				       struct constraint *,
				       double *,
				       bitmap_t *,
				       struct cinfo *);
struct constraint *	check_unique_subtour (bitmap_t *,
					      int,
					      struct constraint *);
struct constraint *	enumerate_all_subtours (struct comp *,
						struct constraint *,
						struct bbinfo *);
struct constraint *	find_integer_cycles (double *,
					     bitmap_t *,
					     bitmap_t *,
					     struct constraint *,
					     struct cinfo *);
struct constraint *	find_small_subtours (struct comp *,
					     struct constraint *,
					     struct bbinfo *);
bool			is_equal (bitmap_t *, bitmap_t *, int);
struct constraint *	sec_flow_heuristic (struct comp *,
					    double *,
					    bitmap_t *,
					    struct cinfo *,
					    struct constraint *);


/*
 * External References
 */

	/* None */


/*
 * Local Equates
 */

#define CYCLE_LIMIT	250


/*
 * Local Types
 */

struct sec_heur_info {

	/* Data used by the flow solver... */
	struct flow_prob	prob;	/* The network flow formulation */
	struct flow_soln	soln;	/* The network flow solution */
	struct flow_temp	temp;	/* Temporary data structures */

	/* Data used to set the arc capacities and modify the	*/
	/* flow network during SEC separation. */
	int		src_arc1;	/* First source arc */
	int		dst_arc1;	/* First destination arc */
	int **		edge_lists;	/* edges for each arc */
};



/*
 * Local Routines
 */

static void			build_heuristic_SEC_formulation (
					struct comp *		comp,
					struct cinfo *		cip,
					struct sec_heur_info *	hsecp);
static struct pset *		compute_vertex_coords (struct comp *,
							 struct cinfo *);
static struct constraint *	cwalk (int,
				       int,
				       bitmap_t *,
				       bitmap_t *,
				       bitmap_t *,
				       bitmap_t *,
				       bitmap_t *,
				       bitmap_t *,
				       bitmap_t *,
				       int *,
				       int *,
				       int *,
				       struct constraint *,
				       struct cinfo *);
#ifdef OLD_ENUMERATION
static struct constraint *	enumerate_subtours (int,
						    int,
						    double,
						    bitmap_t *,
						    int *,
						    int *,
						    struct comp *,
						    struct constraint *,
						    struct bbinfo *);
#endif
static struct constraint *	find_almost_integer_cycles (int,
							    int *,
							    bitmap_t *,
							    bitmap_t *,
							    bitmap_t *,
							    struct cinfo *);
static void			free_heuristic_SEC_formulation (
					struct sec_heur_info *	hsecp);
static int			full_set_mst (struct comp *,
					      int,
					      struct pset *,
					      struct edge *);
static bool			fwalk (int,
				       int,
				       bitmap_t *,
				       bitmap_t *,
				       struct cinfo *);
static struct constraint *	heuristically_find_violated_secs (
					struct comp *		comp,
					bitmap_t *		fset_mask,
					struct cinfo *		cip,
					struct sec_heur_info *	hsecp,
					double *		x,
					struct constraint *	clist);
#ifndef OLD_ENUMERATION
static struct constraint *	recurse_enum (int	limit,
					      int	navail,
					      int *	avail,
					      int	nchosen,
					      int *	chosen,
					      int	nexcl,
					      int *	excluded,
					      int *	vstat,
					      double	lhs,
					      double	rhs,
					      bool	supersets,
					      struct comp *	comp,
					      struct constraint * cp,
					      struct bbinfo *	bbip);
#endif
static void			set_arc_capacities (struct comp *,
						    struct sec_heur_info *);
static void			sort_edges (int, struct edge *, int *);


/*
 * Local Variables
 */

	/* none */

/*
 * This routine heuristically reduces the weighted hypergraph of FSTs
 * to a conventional weighted graph by computing a minimum spanning
 * tree for each FST.  We then use the flow formulation of
 * Padberg ("Trees and Cuts" paper) to find SEC violations in the
 * weighted graph.  This method is heuristic since violations will be
 * discovered or not depending on the particular reduction from
 * hypergraph to graph that is chosen.
 */

	struct constraint *
sec_flow_heuristic (

struct comp *		comp,		/* IN - component to separate. */
double *		x,		/* IN - current LP solution. */
bitmap_t *		edge_mask,	/* IN - subset of edges. */
struct cinfo *		cip,		/* IN - compatability info. */
struct constraint *	cp		/* IN - existing constraints. */
)
{
struct sec_heur_info	sec_heur_info;

	if (cip -> pts EQ NULL) {
		/* This heuristic needs coordinates for each vertex! */
		return (cp);
	}

	if (comp -> num_verts > 1) {
		/* Build an SEC formulation for this component... */
		build_heuristic_SEC_formulation (comp,
						 cip,
						 &sec_heur_info);

		/* Do the heuristic separation... */
		cp = heuristically_find_violated_secs (comp,
						       edge_mask,
						       cip,
						       &sec_heur_info,
						       x,
						       cp);

		free_heuristic_SEC_formulation (&sec_heur_info);
	}

	return (cp);
}

/*
 * This routine builds the initial formulation for the heuristic
 * separation procedure.
 */

	static
	void
build_heuristic_SEC_formulation (

struct comp *		comp,		/* IN - congested component. */
struct cinfo *		cip,		/* IN - compatibility info. */
struct sec_heur_info *	hsecp		/* OUT - SEC flow formulation. */
)
{
int			i;
int			j;
int			k;
int			n;
int			nverts;
int			nedges;
int			num_fst_edges;
int			num_unique_edges;
int			total_edges;
int			num_arcs;
int			num_nodes;
int			arc_num;
struct pset *		pts;
struct flow_prob *	prob;
struct full_set *	fsp;
struct edge *		fs_mst_edges;
int *			fs_num;
struct edge *		ep1;
struct edge *		ep2;
struct edge *		ep_endp;
int *			nump;
int *			outlist;
int *			inlist;
int *			fslist;
int *			srcp;
int *			dstp;
int *			fs_ip;
int *			ip;
int *			counts;
int **			ptrs;

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

	/* Compute actual X,Y coordinates for each vertex in	*/
	/* the congested component...				*/
	pts = compute_vertex_coords (comp, cip);

	/* Tally up the total number of edges needed to compute	*/
	/* a minimum spanning tree for the specified vertices	*/
	/* of each edge.  This is the sum of (|e| - 1) for	*/
	/* all edges e...					*/
	total_edges = (comp -> everts [nedges] - comp -> everts [0])
		      - nedges;

	/* Allocate temporary arrays to hold the edges of the	*/
	/* minimum spanning trees of each FST, together with	*/
	/* the FST number for the edge.				*/
	fs_mst_edges	= NEWA (total_edges, struct edge);
	fs_num		= NEWA (total_edges, int);

	/* Compute a Minimum Spanning Tree for each FST.	*/
	/* Record the edges for each, and the FST number	*/
	/* that the edge came from.				*/
	ep1 = fs_mst_edges;
	nump = fs_num;
	for (i = 0; i < nedges; i++) {
		num_fst_edges = full_set_mst (comp, i, pts, ep1);
		ep1 += num_fst_edges;
		for (j = 0; j < num_fst_edges; j++) {
			*nump++ = i;
		}
	}

	if ((ep1 NE (fs_mst_edges + total_edges)) OR
	    (nump NE (fs_num + total_edges))) {
		fatal ("build_heuristic_SEC_formulation: Bug 1.");
	}

	free ((char *) pts);

	/* Sort the list of all MST edges into increasing order	*/
	/* (primary key = low endpoint, secondary key = high	*/
	/* endpoint).						*/
	sort_edges (total_edges, fs_mst_edges, fs_num);

	/* Scan through the sorted edges and count the number	*/
	/* of unique (a,b) pairs.				*/
	if (total_edges <= 0) {
		num_unique_edges = 0;
	}
	else {
		num_unique_edges = 1;
		ep1 = &fs_mst_edges [1];
		for (i = 1; i < total_edges; i++, ep1++) {
			if ((ep1 [-1].p1 NE ep1 [0].p1) OR
			    (ep1 [-1].p2 NE ep1 [0].p2)) {
				++num_unique_edges;
			}
		}
	}

	/* Compute the total number of directed arcs in the	*/
	/* flow graph.  Each undirected edge becomes a pair of	*/
	/* head-to-tail directed arcs, plus 2 extra arcs per	*/
	/* vertex -- one from the source node and one to the	*/
	/* sink node.						*/
	num_nodes		= nverts + 2;
	num_arcs		= 2 * num_unique_edges + 2 * nverts;

	/* Start filling in the flow problem instance... */
	prob = &(hsecp -> prob);
	prob -> num_nodes	= num_nodes;
	prob -> num_arcs	= num_arcs;

	/* Assign node numbers for the source and sink nodes... */
	prob -> source	= nverts;
	prob -> sink	= nverts + 1;

	/* Now that we know how big the directed flow graph is,	*/
	/* allocate storage for the various data structures...	*/

	prob -> out		= NEWA (num_nodes + 1, int *);
	prob -> in		= NEWA (num_nodes + 1, int *);
	hsecp -> edge_lists	= NEWA (2 * num_unique_edges + 1, int *);
	prob -> arc_src		= NEWA (num_arcs, int);
	prob -> arc_dst		= NEWA (num_arcs, int);

	outlist			= NEWA (num_arcs, int);
	inlist			= NEWA (num_arcs, int);
	fslist			= NEWA (2 * total_edges, int);

	prob -> capacity	= NEWA (num_arcs, double);

	/* Loop through all of the undirected edges, generating	*/
	/* the pair of directed arcs for each.  For each arc	*/
	/* generated, generate the list of edges from which	*/
	/* it draws capacity.					*/
	srcp = prob -> arc_src;
	dstp = prob -> arc_dst;

	ep1	= &fs_mst_edges [0];
	ep_endp	= &fs_mst_edges [total_edges];
	fs_ip	= fs_num;
	arc_num	= 0;
	for (i = 0; i < total_edges;) {
		/* Find number of consecutive "copies" of this edge... */
		ep2 = ep1 + 1;
		n = 1;
		for (;;) {
			if (ep2 >= ep_endp) break;
			if (ep1 -> p1 NE ep2 -> p1) break;
			if (ep1 -> p2 NE ep2 -> p2) break;
			++ep2;
			++n;
		}
		/* Record the A -> B arc... */
		*srcp++ = ep1 -> p1;
		*dstp++ = ep1 -> p2;
		hsecp -> edge_lists [arc_num] = fslist;
		for (j = 0; j < n; j++) {
			*fslist++ = fs_ip [j];
		}
		++arc_num;

		/* Record the B -> A arc... */
		*srcp++ = ep1 -> p2;
		*dstp++ = ep1 -> p1;
		hsecp -> edge_lists [arc_num] = fslist;
		for (j = 0; j < n; j++) {
			*fslist++ = fs_ip [j];
		}
		++arc_num;

		/* Advance to the next undirected edge. */
		ep1	+= n;
		fs_ip	+= n;
		i	+= n;
	}
	/* Terminate list of edges comprising each arc. */
	hsecp -> edge_lists [arc_num] = fslist;

	free ((char *) fs_num);
	free ((char *) fs_mst_edges);

	/* Record the source -> vertex arcs... */
	hsecp -> src_arc1 = arc_num;
	for (i = 0; i < nverts; i++) {
		*srcp++ = prob -> source;
		*dstp++ = i;
		++arc_num;
	}

	/* Record the vertex -> sink arcs... */
	hsecp -> dst_arc1 = arc_num;
	for (i = 0; i < nverts; i++) {
		*srcp++ = i;
		*dstp++ = prob -> sink;
		++arc_num;
	}

	if (arc_num NE num_arcs) {
		fatal ("build_heuristic_SEC_formulation: Bug 2.");
	}

	/* We have now specified the directed flow graph as a	*/
	/* list of directed arcs.  Time to construct the	*/
	/* adjacency lists -- for each node we build a list of	*/
	/* outgoing and incoming arc numbers.  Do the outgoing	*/
	/* lists first...					*/

	counts	= NEWA (num_nodes, int);
	ptrs	= NEWA (num_nodes, int *);

	for (i = 0; i < num_nodes; i++) {
		counts [i] = 0;
	}
	for (i = 0; i < num_arcs; i++) {
		++(counts [prob -> arc_src [i]]);
	}
	ip = outlist;
	for (i = 0; i < num_nodes; i++) {
		ptrs [i] = ip;
		prob -> out [i] = ip;
		ip += counts [i];
	}
	prob -> out [i] = ip;
	for (i = 0; i < num_arcs; i++) {
		j = prob -> arc_src [i];
		ip = ptrs [j]++;
		*ip = i;
	}

	/* Now do the incoming arc lists... */
	for (i = 0; i < num_nodes; i++) {
		counts [i] = 0;
	}
	for (i = 0; i < num_arcs; i++) {
		++(counts [prob -> arc_dst [i]]);
	}
	ip = inlist;
	for (i = 0; i < num_nodes; i++) {
		ptrs [i] = ip;
		prob -> in [i] = ip;
		ip += counts [i];
	}
	prob -> in [i] = ip;
	for (i = 0; i < num_arcs; i++) {
		k = prob -> arc_dst [i];
		ip = ptrs [k]++;
		*ip = i;
	}

	/* Free temporary memory used to build things... */
	free ((char *) counts);
	free ((char *) ptrs);

	/* Initialize the buffers used to hold flow solutions */
	/* and temporary data... */
	create_flow_solution_data (prob, &(hsecp -> soln));
	create_flow_temp_data (prob, &(hsecp -> temp));
}

/*
 * This routine computes an X,Y coordinate for each vertex in the
 * congested component.  This is tricky since there may be "vertices"
 * in the component that represent the MERGING of numerous "real"
 * vertices from the problem data.  We handle this by computing the
 * AVERAGE x,y coordinate of each real vertex in the component vertex.
 */

	static
	struct pset *
compute_vertex_coords (

struct comp *		comp,		/* IN - congested component */
struct cinfo *		cip		/* IN - compatibility info */
)
{
int			i;
int			j;
int			k;
int			nverts;
int *			ip1;
int *			ip2;
struct pset *		pts;
struct point *		p1;
coord_t			x;
coord_t			y;

	nverts = comp -> num_verts;

	pts = (struct pset *) new (PSET_SIZE (nverts));

	pts -> n = nverts;

	for (i = 0; i < nverts; i++) {
		/* Compute the "average" of the X,Y coordinates of each	*/
		/* "real" vertex embedded in this "fake" congested	*/
		/* component vertex.  We probably don't even really	*/
		/* care if these additions overflow -- we are simply	*/
		/* choosing SOME spanning tree for the FST...		*/
		x = 0;
		y = 0;
		k = 0;
		ip1 = comp -> rverts [i];
		ip2 = comp -> rverts [i + 1];
		while (ip1 < ip2) {
			j = *ip1++;
			p1 = &(cip -> pts -> a [j]);
			x += p1 -> x;
			y += p1 -> y;
			++k;
		}
		if (k > 0) {