cutset.c

生成直角Steiner树的程序包

				  vert_mask,
				  edge_mask,
				  cip);

	cstack [i] = TRUE;
	cutlist = enumerate_cuts (i + 1,
				  ntaken + 1,
				  cstack,
				  cutlist,
				  comps,
				  ncomps,
				  cut_terms,
				  x,
				  vert_mask,
				  edge_mask,
				  cip);

	return (cutlist);
}

/*
 * This routine generates a single cutset constraint for each
 * component (and merges appropriately similar cuts).  This method
 * can be MUCH faster than the enumeration of all cuts!
 */

	static
	struct constraint *
simple_cuts (

struct constraint *	cutlist,	/* IN - list of cutsets */
bitmap_t *		comps,		/* IN - connected components */
int			ncomps,		/* IN - number of components */
double *		x,		/* IN - current LP solution */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip		/* IN - compatibility info */
)
{
int			i;
int			kmasks;
bitmap_t *		cut_terms;

	kmasks = cip -> num_vert_masks;

	for (i = 0; i < ncomps; i++) {
		/* The terminals in this component. */
		cut_terms = &comps [kmasks * i];

		/* Add this cutset to the list. */
		cutlist = add_cutset_to_list (cut_terms,
					      cutlist,
					      x,
					      vert_mask,
					      edge_mask,
					      cip);
	}

	return (cutlist);
}

/*
 * This routine formulates the Linear Program used to find violated
 * cutset constraints.
 */

	void
build_cutset_separation_formulation (

bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip,		/* IN - compatibility info. */
struct cs_info *	csip		/* OUT - cutset flow formulation. */
)
{
int			i;
int			j;
int			k;
int			n;
int			nedges;
int			nverts;
int			sumsizes;
int			num_used_trees;
int			num_nodes;
int			num_arcs;
struct flow_prob *	prob;
int			p0;
int			q0;
int			arc_num;
struct pset *		pts;
int *			outlist;
int *			inlist;
int *			counts;
int **			ptrs;
int *			vp1;
int *			vp2;
int *			ip;

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

	/* Compute the sum of all hyperedge cardinalities. */
	sumsizes = 0;
	num_used_trees = 0;
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		sumsizes += cip -> edge_size [i];
		++num_used_trees;
	}

#if 0
tracef (" %% nedges = %d, nverts = %d, sumsizes = %d, num_used_trees = %d\n",
	nedges, nverts, sumsizes, num_used_trees);
#endif

	/* The network we construct contains two additional	*/
	/* nodes Pi and Qi per full set Fi (in addition to the	*/
	/* terminals), and the following set of edges for each	*/
	/* full set Fi:						*/
	/*							*/
	/*	- 1 edge from Pi to Qi.				*/
	/*	- 2 other edges per terminal Tj in Fi:		*/
	/*		- from Tj to Pi				*/
	/*		- from Qi to Tj				*/
	/*							*/
	/* NOTE: we have nodes Pi and Qi regardless of whether	*/
	/*	 or not full set Fi is in the problem -- this	*/
	/*	 makes it easier to compute the proper node	*/
	/*	 numbers for Pi and Qi...  We also have the	*/
	/*	 corresponding Pi -> Qi arc regardless of	*/
	/*	 whether or not full set Fi is in the problem.	*/

	num_nodes = nverts + 2 * nedges;
	num_arcs = nedges + 2 * sumsizes;

	prob = &(csip -> prob);

	prob -> num_nodes	= num_nodes;
	prob -> num_arcs	= num_arcs;
	prob -> source		= 0;		/* specified dynamically */
	prob -> sink		= 0;		/* specified dynamically */
	prob -> arc_src		= NEWA (num_arcs, int);
	prob -> arc_dst		= NEWA (num_arcs, int);
	prob -> capacity	= NEWA (num_arcs, double);

	csip -> arc_to_fset	= NEWA (num_arcs, int);

	/* The terminals become nodes 0 through nverts-1,	*/
	/* The "Pi" nodes start at nverts, and the "Qi" nodes	*/
	/* follow.						*/
	p0	= nverts;
	q0	= p0 + nedges;

	/* Generate the proper graph widget for each full set	*/
	/* that is properly a part of the problem...		*/
	/* The arcs numbered 0..nedges-1 are the Pi -> Qi arcs,	*/
	/* giving the total flow through each full set.		*/
	arc_num = 0;
	for (i = 0; i < nedges; i++) {
		/* Generate Pi -> Qi arc. */
		prob -> arc_src [arc_num] = p0 + i;
		prob -> arc_dst [arc_num] = q0 + i;
		csip -> arc_to_fset [arc_num] = i;
		++arc_num;
	}

	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;

		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;

			/* Generate Tj -> Pi arc. */
			prob -> arc_src [arc_num] = j;
			prob -> arc_dst [arc_num] = p0 + i;
			csip -> arc_to_fset [arc_num] = i;
			++arc_num;

			/* Generate Qi -> Tj arc. */
			prob -> arc_src [arc_num] = q0 + i;
			prob -> arc_dst [arc_num] = j;
			csip -> arc_to_fset [arc_num] = i;
			++arc_num;
		}
	}

	if (arc_num NE num_arcs) {
		fatal ("build_cutset_separation_formulation: Bug 1.");
	}

	/* 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...					*/

	outlist	= NEWA (num_arcs, int);
	inlist	= NEWA (num_arcs, int);

	prob -> out = NEWA (num_nodes + 1, int *);
	prob -> in  = NEWA (num_nodes + 1, int *);

	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, &(csip -> soln));
	create_flow_temp_data (prob, &(csip -> temp));
}

/*
 * This routine performs the separation procedure for cutset constraints.
 * Given an LP solution to the main problem, the task is to find one or
 * more cutsets constraints that are violated (have total weight less
 * than 1) -- or to prove that no such cutsets exist.  We do this by picking
 * the first valid terminal, and computing the max-flow/min-cut to each of
 * the other N-1 terminals.  If the flow is less that 1, we have a
 * violated constraint.  The actual cutset of full-sets consists of those
 * full-sets that span the cut.
 */

	struct constraint *
find_fractional_cutsets (

double *		x,		/* IN - LP solution to separate. */
struct cs_info *	csip,		/* IN - cutset separation info. */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip		/* IN - compatibility info. */
)
{
int			i;
int			kmasks;
int			nverts;
int			num_arcs;
double *		capacity;
int *			arc_to_fset;
int			I;
int			J;
struct constraint *	cutlist;
double			z;

	nverts	= cip -> num_verts;
	kmasks	= cip -> num_vert_masks;

	num_arcs	= csip -> prob.num_arcs;
	capacity	= csip -> prob.capacity;
	arc_to_fset	= csip -> arc_to_fset;

	/* Distribute the LP solution weight for each full set	*/
	/* to every edge in the max-flow network that belongs	*/
	/* to that full set.					*/
	for (i = 0; i < num_arcs; i++) {
		capacity [i] = x [arc_to_fset [i]];
	}

	/* Find first valid terminal. */
	I = 0;
	for (;;) {
		if (I >= nverts) {
			/* No valid terminal found... */
			fatal ("find_fractional_cutsets: Bug 1.");
		}
		if (BITON (vert_mask, I)) break;
		++I;
	}

	cutlist = NULL;

	/* Loop through remaining terminals J.  Do a single max-flow/	*/
	/* min-cut problem for each one...				*/
	for (J = I + 1; J < nverts; J++) {
		if (NOT BITON (vert_mask, J)) continue;

		csip -> prob.source = J;
		csip -> prob.sink   = I;

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

		z = csip -> soln.z;

		if ((1.0 - z) < FUZZ) continue;

		/* We have a violated constraint!  Keep	*/
		/* only bits for valid terminals...	*/
		for (i = 0; i < kmasks; i++) {
			csip -> soln.cut [i] &= vert_mask [i];
		}

		/* Add this cutset to the list! */
		cutlist = add_cutset_to_list (csip -> soln.cut,
					      cutlist,
					      x,
					      vert_mask,
					      edge_mask,
					      cip);
	}

	return (cutlist);
}