cutset.c

生成直角Steiner树的程序包

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

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

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

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

	Separation routines for cutset constraints.

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

	Modification Log:

	b-1:	11/14/96	warme
		: Split off from bb.c.
	b-2:	02/28/2001	warme
		: Changes for 3.1 release.
		: Split add_cutset_to_list, and its subroutines
		:  off to new cutsubs.c file.

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


/*
 * Global Routines
 */

void		build_cutset_separation_formulation (bitmap_t *,
						     bitmap_t *,
						     struct cinfo *,
						     struct cs_info *);
struct constraint * find_cutset_constraints (double *,
					     bitmap_t *,
					     bitmap_t *,
					     struct cinfo *);
struct constraint * find_fractional_cutsets (double *,
					     struct cs_info *,
					     bitmap_t *,
					     bitmap_t *,
					     struct cinfo *);


/*
 * Local Routines
 */

static struct constraint * enumerate_cuts (int,
					   int,
					   bool *,
					   struct constraint *,
					   bitmap_t *,
					   int,
					   bitmap_t *,
					   double *,
					   bitmap_t *,
					   bitmap_t *,
					   struct cinfo *);
static int		find_comps (bitmap_t *, bitmap_t *, struct cinfo *);
static struct constraint * simple_cuts (struct constraint *,
					bitmap_t *,
					int,
					double *,
					bitmap_t *,
					bitmap_t *,
					struct cinfo *);

/*
 * This routine quickly finds cutsets of zero weight -- totally
 * disconnected solutions.  It first computes the connected components
 * of the solution and then uses either a combinatorially thorough
 * method to generate a complete set of constraints, or a quicker
 * method to generate one constraint per component.  The method used
 * depends upon the number of connected components found.
 */

	struct constraint *
find_cutset_constraints (

double *		x,		/* IN - 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			nverts;
int			nedges;
int			kmasks;
int			nmasks;
int			ncomps;
bitmap_t *		sol;
bitmap_t *		cut_terms;
struct constraint *	cutlist;
bitmap_t *		comps;
bool *			cstack;

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

	comps	= NEWA (nverts * kmasks, bitmap_t);
	sol	= NEWA (nmasks, bitmap_t);

	cutlist = NULL;

	/* Get bit-mask of all full-sets even partially present in soln. */
	memset (sol, 0, nmasks * sizeof (*sol));
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		if (x [i] > FUZZ) {
			SETBIT (sol, i);
		}
	}

	/* Determine the number of connected components it has.	*/
	ncomps = find_comps (sol, comps, cip);
	if (ncomps <= 0) {
		fatal ("find_cutset_constraints: Bug 2.");
	}
	if (ncomps EQ 1) {
		free ((char *) sol);
		free ((char *) comps);
		return (NULL);
	}

#if 1
	tracef (" %% @cutset: %d connected components.\n", ncomps);
#endif

	if (ncomps >= 12) {
		/* Too many components to bother trying all combinations */
		/* of them -- just cut around each component...		 */
		cutlist = simple_cuts (cutlist,
				       comps,
				       ncomps,
				       x,
				       vert_mask,
				       edge_mask,
				       cip);
	}
	else {
		/* Enumerate all cuts induced by the connected	*/
		/* components.  For each cut, we may have to	*/
		/* generate a constraint (if it has not been	*/
		/* generated before).				*/
		cut_terms  = NEWA (kmasks, bitmap_t);
		cstack = NEWA (nverts, bool);
		cstack [0] = TRUE;	/* Always take 1st component... */
		cutlist = enumerate_cuts (1,
					  1,
					  cstack,
					  cutlist,
					  comps,
					  ncomps,
					  cut_terms,
					  x,
					  vert_mask,
					  edge_mask,
					  cip);
		free ((char *) cstack);
		free ((char *) cut_terms);
	}

	free ((char *) sol);
	free ((char *) comps);

	return (cutlist);
}

/*
 * This routine finds the connected components of the given set of
 * full-sets.  The result is a partition of the given terminals that
 * we represent as an array of terminal bit-masks, one per component.
 */

	static
	int
find_comps (

bitmap_t *	sol,		/* IN - solution to get components of. */
bitmap_t *	comps,		/* OUT - components. */
struct cinfo *	cip		/* IN - compatibility info. */
)
{
int			i;
int			j;
int			e;
int			v;
int			nverts;
int			nedges;
int			kmasks;
int			nmasks;
int			ncomps;
int *			ep1;
int *			ep2;
int *			vp1;
int *			vp2;
bitmap_t *		visited;
int *			sp;
int *			stack;

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

	visited = NEWA (kmasks, bitmap_t);
	stack	= NEWA (nverts, int);

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

	ncomps = 0;
	for (;;) {
		/* Find next remaining full-set... */
		for (e = 0; e < nedges; e++) {
			if (BITON (sol, e)) break;
		}
		if (e >= nedges) break;
		CLRBIT (sol, e);

		/* Identify next component... */
		for (i = 0; i < kmasks; i++) {
			comps [i] = 0;
		}
		sp = &stack [0];
		vp1 = cip -> edge [e];
		vp2 = cip -> edge [e + 1];
		while (vp1 < vp2) {
			v = *vp1++;
			if (BITON (visited, v)) continue;
			SETBIT (visited, v);
			SETBIT (comps, v);
			*sp++ = v;
		}

		while (sp > stack) {
			v = *--sp;
			ep1 = cip -> term_trees [v];
			ep2 = cip -> term_trees [v + 1];
			while (ep1 < ep2) {
				e = *ep1++;
				if (NOT BITON (sol, e)) continue;
				CLRBIT (sol, e);
				vp1 = cip -> edge [e];
				vp2 = cip -> edge [e + 1];
				while (vp1 < vp2) {
					v = *vp1++;
					if (BITON (visited, v)) continue;
					SETBIT (visited, v);
					SETBIT (comps, v);
					*sp++ = v;
				}
			}
		}
		++ncomps;
		comps += kmasks;
	}

	free ((char *) stack);
	free ((char *) visited);

	return (ncomps);
}

/*
 * This routine recursively enumerates all non-trivial partitions of the
 * connected components into two disjoint groups.  Such a partition
 * defines a CUT of the original Steiner tree problem.  For each such
 * cut, we determine the set of all full-sets that span the cut.  If
 * this cutset has not been seen before, we add it to the cutlist.
 */

	static
	struct constraint *
enumerate_cuts (

int			i,		/* IN - current recursion level */
int			ntaken,		/* IN - number of components taken */
bool *			cstack,		/* IN - which components taken */
struct constraint *	cutlist,	/* IN - list of cutsets */
bitmap_t *		comps,		/* IN - connected components */
int			ncomps,		/* IN - number of components */
bitmap_t *		cut_terms,	/* IN - another temporary cut-set */
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			j;
int			kmasks;
bitmap_t *		bp1;
bitmap_t *		bp2;
bitmap_t *		bp3;

	if (i >= ncomps) {
		/* Base case. */
		if ((ntaken <= 0) OR (ntaken >= ncomps)) {
			/* Ignore trivial partitions... */
			return (cutlist);
		}

		kmasks = cip -> num_vert_masks;

		/* Compute the set of terminals in the cut... */
		memset (cut_terms, 0, kmasks * sizeof (*cut_terms));
		for (j = 0; j < ncomps; j++) {
			if (NOT cstack [j]) continue;
			bp1 = cut_terms;
			bp2 = bp1 + kmasks;
			bp3 = &comps [j * kmasks];
			while (bp1 < bp2) {
				*bp1++ |= *bp3++;
			}
		}

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

		return (cutlist);
	}

	/* Recurse here... */
	cstack [i] = FALSE;
	cutlist = enumerate_cuts (i + 1,
				  ntaken,
				  cstack,
				  cutlist,
				  comps,
				  ncomps,
				  cut_terms,
				  x,