sec_heur.c

生成直角Steiner树的程序包

{
int			i;
int			j;
int			nverts;
int			nedges;
int			kmasks;
int			limit;
int *			ip1;
int *			ip2;
int *			ip3;
bitmap_t *		bp1;
bitmap_t *		bp2;
struct constraint *	newp;
double			sum;

	nverts = comp -> num_verts;

	if (k <= 0) {
		/* We have chosen a complete subset of the desired	*/
		/* cardinality.  Check it for violations...		*/
		nedges = comp -> num_edges;
		sum = 0.0;
		for (i = 0; i < nedges; i++) {
			ip1 = comp -> everts [i];
			ip2 = comp -> everts [i + 1];
			k = -1;
			while (ip1 < ip2) {
				t = *ip1++;
				if (BITON (chosen, t)) {
					++k;
				}
			}
			if (k > 0) {
				sum += (k * (comp -> x [i]));
			}
		}
		ip1 = start;
		while (ip1 < end) {
			t = *ip1++;
			sum += comp -> tviol [t];
		}

		if (sum > rhs) {
			/* We have a violation!  Assemble a mask of the	*/
			/* actual vertices in the subtour from the	*/
			/* chosen vertices of the congested component. */
			kmasks = bbip -> cip -> num_vert_masks;
			bp1 = NEWA (kmasks, bitmap_t);

			for (i = 0; i < kmasks; i++) {
				bp1 [i] = 0;
			}
			ip1 = start;
			while (ip1 < end) {
				i = *ip1++;
				ip2 = comp -> rverts [i];
				ip3 = comp -> rverts [i + 1];
				while (ip2 < ip3) {
					j = *ip2++;
					SETBIT (bp1, j);
				}
			}

			newp = NEW (struct constraint);
			newp -> next		= cp;
			newp -> iteration	= 0;
			newp -> type		= CT_SUBTOUR;
			newp -> mask		= bp1;
			cp = newp;
#if 0
			print_mask (" %% subtour:", bp1, kmasks * BPW);
#endif
		}

		return (cp);
	}

	limit = nverts - k;

	while (t <= limit) {
		/* Try taking the current vertex... */
		SETBIT (chosen, t);
		*end = t;
		cp = enumerate_subtours (k - 1,
					 t + 1,
					 rhs,
					 chosen,
					 start,
					 end + 1,
					 comp,
					 cp,
					 bbip);
		CLRBIT (chosen, t);
		++t;
	}

	return (cp);
}

/*
 * This routine finds violated SEC's by exhaustively enumerating all
 * subsets of size 2 or more.  We return a list of the constraints that
 * were found.
 */

#else

	struct constraint *
enumerate_all_subtours (

struct comp *		comp,		/* IN - component to check. */
struct constraint *	cp,		/* IN - existing list of constraints */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			nverts;
int *			avail;
int *			chosen;
int *			excluded;
int *			vstat;

	nverts = comp -> num_verts;

#if 0
	tracef (" %% Exhaustively enumerating %d congested vertices.\n", nverts);
#endif

	avail	 = NEWA (4 * nverts, int);
	chosen	 = avail + nverts;
	excluded = chosen + nverts;
	vstat	 = excluded + nverts;

	/* Each vertex is initially free... */
	for (i = 0; i < nverts; i++) {
		vstat [i] = -1;
	}

	for (i = 0; i < nverts; i++) {

		chosen [0]	= i;
		vstat [i]	= 1;

		cp = recurse_enum (nverts,	/* limit */
				   0,		/* navail */
				   avail,
				   1,		/* nchosen */
				   chosen,
				   i,		/* nexcl */
				   excluded,
				   vstat,
				   0.0,		/* LHS */
				   FUZZ,	/* RHS */
				   TRUE,	/* do_supersets */
				   comp,
				   cp,
				   bbip);
		excluded [i] = i;
		vstat [i] = 2;
	}

	free ((char *) avail);

	return (cp);
}

/*
 * This routine finds violated SEC's by explicitly enumerating subsets
 * of increasing cardinality.  Only the vertices of the given congested
 * component are examined.  We return a list of the constraints that
 * were found.
 */

	struct constraint *
find_small_subtours (

struct comp *		comp,		/* IN - component to check. */
struct constraint *	cp,		/* IN - existing list of constraints */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			k;
int			nverts;
int			klimit;
int			ticks_limit;
cpu_time_t		t0;
cpu_time_t		t1;
int *			avail;
int *			chosen;
int *			excluded;
int *			vstat;
double			est;

	nverts = comp -> num_verts;

	avail	 = NEWA (4 * nverts, int);
	chosen	 = avail + nverts;
	excluded = chosen + nverts;
	vstat	 = excluded + nverts;

	tracef (" %% Enumerating %d congested vertices.\n", nverts);

	t0 = get_cpu_time ();

	/* Determine cardinality limit for partial enumeration. */

	if (nverts <= SEC_ENUM_LIMIT) {
		klimit = nverts;
	}
	else if (nverts <= 2 * SEC_ENUM_LIMIT) {
		klimit = 5;
	}
	else if (nverts <= 3 * SEC_ENUM_LIMIT) {
		klimit = 3;
	}
	else {
		klimit = 2;
	}

	/* Establish time limit for each cardinality. */
	ticks_limit = nverts * TICKS_PER_SEC;

	/* Don't need to check 2-vertex subtours if they are	*/
	/* already present in the constraint pool!		*/
	k = Seed_Pool_With_2SECs ? 3 : 2;

	for (; k <= klimit; k++) {
#if 0
		tracef (" %% Checking %d subtours\n", k);
#endif

		/* Each vertex is initially free... */
		for (i = 0; i < nverts; i++) {
			vstat [i] = -1;
		}

		for (i = 0; i < nverts; i++) {

			chosen [0]	= i;
			vstat [i]	= 1;

			cp = recurse_enum (k,		/* limit */
					   0,		/* navail */
					   avail,
					   1,		/* nchosen */
					   chosen,
					   i,		/* nexcl */
					   excluded,
					   vstat,
					   0.0,		/* LHS */
					   FUZZ,	/* RHS */
					   FALSE,	/* do_supersets */
					   comp,
					   cp,
					   bbip);
			excluded [i] = i;
			vstat [i] = 2;
		}

		/* Determine how long this took... */
		t1 = get_cpu_time ();
#if 0
		{ char	tbuf [64];
		convert_cpu_time (t1 - t0, tbuf);
		(void) tracef (" %%	%s seconds\n", tbuf);
		}
#endif

		/* if we found any violations, get out! */
		if (cp NE NULL) break;

		/* Estimate how long the next pass will take... */
		est = ((double) (nverts - k)) / ((double) (k + 1));
		est *= (t1 - t0);
		if (est > ticks_limit) {
			/* Do not take more than specified	*/
			/* number of ticks...			*/
			break;
		}
		t0 = t1;
	}

	free ((char *) avail);

	return (cp);
}

/*
 * This routine recursively enumerates connected subsets of size k of the
 * congested vertices.
 */

	static
	struct constraint *
recurse_enum (

int			limit,		/* IN - max num verts to choose */
int			navail,		/* IN - num verts available */
int *			avail,		/* IN - unchosen verts available */
int			nchosen,	/* IN - num verts chosen */
int *			chosen,		/* IN - vertices chosen */
int			nexcl,		/* IN - num vertices excluded */
int *			excluded,	/* IN - verts excluded from choice */
int *			vstat,		/* IN - status of each vertex: */
					/* -1=free, 0=avail, 1=chosen, */
					/* 2=excluded. */
double			lhs,		/* IN - LHS of constraint */
double			rhs,		/* IN - RHS of constraint */
bool			do_supersets,	/* IN - do supersets of violations? */
struct comp *		comp,		/* IN - component to enumerate */
struct constraint *	cp,		/* IN - existing constraints */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			j;
int			k;
int			t;
int			e;
int			kmasks;
int			navail_save;
int			new_nexcl;
int *			vp1;
int *			vp2;
int *			ep1;
int *			ep2;
bitmap_t *		bp1;
struct constraint *	newp;
double			sum;

	/* Check if chosen vertices yield a violation: */
	if (lhs > rhs) {
		kmasks = bbip -> cip -> num_vert_masks;
		bp1 = NEWA (kmasks, bitmap_t);

		for (i = 0; i < kmasks; i++) {
			bp1 [i] = 0;
		}
		for (i = 0; i < nchosen; i++) {
			j = chosen [i];
			vp1 = comp -> rverts [j];
			vp2 = comp -> rverts [j + 1];
			while (vp1 < vp2) {
				k = *vp1++;
				SETBIT (bp1, k);
			}
		}

		newp = NEW (struct constraint);
		newp -> next		= cp;
		newp -> iteration	= 0;
		newp -> type		= CT_SUBTOUR;
		newp -> mask		= bp1;

		cp = newp;
#if 0
		print_mask (" %% subtour:", bp1, bbip -> cip -> num_verts);
#endif

		if (NOT do_supersets) {
			/* Don't enumerate supersets of any violation... */
			return (cp);
		}
	}

	if (nchosen >= limit) {
		/* Don't recurse any deeper. */
		return (cp);
	}

	navail_save = navail;

	/* Find newly available vertices.  All free vertices adjacent */
	/* to the most-recently-chosen vertex now become available. */
	t = chosen [nchosen - 1];
	ep1 = comp -> vedges [t];
	ep2 = comp -> vedges [t + 1];
	while (ep1 < ep2) {
		e = *ep1++;
		vp1 = comp -> everts [e];
		vp2 = comp -> everts [e + 1];
		while (vp1 < vp2) {
			i = *vp1++;
			if (vstat [i] < 0) {
				avail [navail++] = i;
				vstat [i] = 0;
			}
		}
	}

	/* Choose each available vertex.  After recursing on	*/
	/* this choice, exclude it from future consideration.	*/

	new_nexcl = nexcl;
	while (navail > 0) {
		/* Choose (pop) last available vertex t. */
		t = avail [--navail];

		/* Compute weight of edges connecting vertex t to all	*/
		/* previously chosen vertices.  This is used to update	*/
		/* the LHS of the constraint.				*/
		sum = 0.0;
		ep1 = comp -> vedges [t];
		ep2 = comp -> vedges [t + 1];
		while (ep1 < ep2) {
			e = *ep1++;
			vp1 = comp -> everts [e];
			vp2 = comp -> everts [e + 1];
			while (vp1 < vp2) {
				j = *vp1++;
				if (vstat [j] EQ 1) {
					sum += comp -> x [e];
					break;
				}
			}
		}

		/* Vertex t is now CHOSEN. */
		chosen [nchosen]	= t;
		vstat [t]		= 1;

		cp = recurse_enum (limit,
				   navail,
				   avail,
				   nchosen + 1,
				   chosen,
				   new_nexcl,
				   excluded,
				   vstat,
				   lhs + sum,
				   rhs + 1.0,
				   do_supersets,
				   comp,
				   cp,
				   bbip);

		/* Now exclude vertex t from further consideration. */
		vstat [t]		= 2;
		excluded [new_nexcl++]	= t;
	}

	/* Restore those vertices that we excluded back to	*/
	/* their available state.				*/
	while (new_nexcl > nexcl) {
		t = excluded [--new_nexcl];
		vstat [t] = 0;
		avail [navail++] = t;
	}

	/* Restore those vertices that we made available back	*/
	/* to their original "free" state.			*/
	while (navail > navail_save) {
		t = avail [--navail];
		vstat [t] = -1;
	}

	return (cp);
}

#endif

/*
 * This routine checks the given subset of vertices to see if it defines
 * a violated subtour constraint.  If so it adds the constraint to the
 * current list of constraints.
 */

	struct constraint *
check_subtour (

bitmap_t *		stour,		/* IN - subset to check */
struct constraint *	clist,		/* IN - existing constraints */
double *		x,		/* IN - current LP solution */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip		/* IN - compatibility info */
)
{
int			i;
int			j;
int			kmasks;
int			nedges;
int			ssize;
int			count;
int *			vp1;
int *			vp2;
bitmap_t *		bp1;
bitmap_t		mask;
double			total;
double			coeff;
struct constraint *	cp;

	nedges = cip -> num_edges;
	kmasks = cip -> num_vert_masks;

	/* Get size of subtour - 1... */
	ssize = -1;
	bp1 = stour;
	for (i = 0; i < kmasks; i++) {
		mask = *bp1++;
		ssize += NBITSON (mask);
	}

	if (ssize <= 0) return (clist);

	total = 0.0;
	for (j = 0; j < nedges; j++) {
		if (NOT BITON (edge_mask, j)) continue;
		count = -1;
		vp1 = cip -> edge [j];
		vp2 = cip -> edge [j + 1];
		while (vp1 < vp2) {
			i = *vp1++;
			if (BITON (stour, i)) {
				++count;
			}
		}
		if (count <= 0) continue;
		coeff = ((double) count);
		total += coeff * x [j];
	}

	if (total <= ((double) ssize) + FUZZ) {
		/* Constraint is not violated... */
		return (clist);
	}

	/* Create new constraint and add it to the list... */

	bp1 = NEWA (kmasks, bitmap_t);
	for (i = 0; i < kmasks; i++) {
		bp1 [i] = stour [i];
	}

	cp = NEW (struct constraint);
	cp -> next	= clist;
	cp -> iteration	= 0;
	cp -> type	= CT_SUBTOUR;
	cp -> mask	= bp1;

#if 0
	print_mask (" %% subtour:", stour, cip -> num_verts);
#endif
#if 0
	debug_print_constraint (" % Subtour:  ",
				" %\t",
				cp,
				x,
				edge_mask,
				cip);
#endif

	return (cp);
}