ub.c

生成直角Steiner树的程序包

					      ranking);

		/* Sort the zero-weight FSTs by rank only. */
		sort_fst_list_by_rank (edges_0, count_0, ranking);

		/* Try several greedy trees using this ordering of FSTs. */
		try_trees (edge_list, nedges, bbip);

		/* Create a second ordering by deleting all of the MST	*/
		/* edges and placing them last.				*/
		ep1 = edge_list;
		ep2 = edge_list + nedges;
		ep3 = temp_edges;
		while (ep1 < ep2) {
			e = *ep1++;
			if (cip -> edge_size [e] <= 2) continue;
			*ep3++ = e;
		}
		ep1 = ubip -> mst_edges;
		ep2 = ep1 + (cip -> num_verts - 1);
		while (ep1 < ep2) {
			*ep3++ = *ep1++;
		}

		/* Try a few greedy trees using this ordering. */
		try_trees (temp_edges, ep3 - temp_edges, bbip);
	}

	free ((char *) temp_edges);
	free ((char *) edge_list);

	if (ubip -> best_z < bbip -> best_z) {
		new_upper_bound (ubip -> best_z, bbip);
	}
}

/*
 * This routine sorts the full sets of integral weight (both 0.0 and 1.0).
 * In this case, only the given ranking is used.
 */

	static
	void
sort_fst_list_by_rank (

int *		array,		/* IN - list of full set numbers to sort */
int		n,		/* IN - number of full sets to sort */
int *		ranking		/* IN - rank ordering of FSTs */
)
{
int		h;
int		t1;
int		t2;
int		key1;
int *		p1;
int *		p2;
int *		p3;
int *		p4;
int *		endp;

	endp = &array [n];

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

	do {
		h = h / 3;
		p4 = &array [h];
		p1 = p4;
		while (p1 < endp) {
			t1 = *p1;
			key1 = ranking [t1];
			p2 = p1;
			while (TRUE) {
				p3 = (p2 - h);
				t2 = *p3;
				if (ranking [t2] <= key1) break;
				*p2 = t2;
				p2 = p3;
				if (p2 < p4) break;
			}
			*p2 = t1;
			++p1;
		}
	} while (h > 1);
}

/*
 * This routine sorts the fractional full sets.  The primary key is LP
 * solution weight, ties are broken using the given ranking.
 */

	static
	void
sort_fst_list_by_lp_and_rank (

int *		array,		/* IN - list of full set numbers to sort */
int		n,		/* IN - number of full sets to sort */
double *	x,		/* IN - LP solution weights */
int *		ranking		/* IN - rank ordering of FSTs */
)
{
int		h;
int		t1;
int		t2;
double		key1;
int *		p1;
int *		p2;
int *		p3;
int *		p4;
int *		endp;

	endp = &array [n];

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

	do {
		h = h / 3;
		p4 = &array [h];
		p1 = p4;
		while (p1 < endp) {
			t1 = *p1;
			key1 = x [t1];
			p2 = p1;
			while (TRUE) {
				p3 = (p2 - h);
				t2 = *p3;
				if (x [t2] > key1) break;
				if ((x [t2] EQ key1) AND
				    (ranking [t2] <= ranking [t1])) break;
				*p2 = t2;
				p2 = p3;
				if (p2 < p4) break;
			}
			*p2 = t1;
			++p1;
		}
	} while (h > 1);
}

/*
 * This routine constructs several trees according to the given sorted
 * list of FSTs.  It uses the greedy Kruskal algorithm to build an
 * initial tree.  It then constructs a sequence of DIFFERENT trees as
 * follows:
 *
 *	1. Find an edge that has NOT been used in any of the trees
 *	   constructed previously by this routine.
 *	2. Build a new tree using the same list, but putting this
 *	   previously unused edge in FIRST.
 *
 * Given M edges, we repeat the edge-forcing operation until we have
 * tried O(log(M)) unused edges, or until we run out of edges to force.
 */

	static
	void
try_trees (

int *			edge_list,	/* IN - ordered list of edges */
int			nedges,		/* IN - number of edges in list */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			j;
int			nverts;
int			nmasks;
int			limit;
struct cinfo *		cip;
struct ubinfo *		ubip;
int *			temp_edges;
bitmap_t *		used;

	ubip	= bbip -> ubip;
	cip	= bbip -> cip;
	nverts	= cip -> num_verts;
	nmasks	= cip -> num_edge_masks;

	/* None of the edges have been used yet. */
	used = NEWA (nmasks, bitmap_t);
	for (i = 0; i < nmasks; i++) {
		used [i] = 0;
	}

	/* Construct the initial tree.  Note which edges were used. */
	ub_kruskal (edge_list, nedges, used, bbip);

	/* Copy the list so that we can quickly add an edge	*/
	/* onto the front of the list...			*/
	temp_edges = NEWA (nedges + 1, int);
	for (i = 0; i < nedges; i++) {
		temp_edges [i + 1] = edge_list [i];
	}

	/* Determine the limit of edges to try. */
	for (limit = 2; (1 << limit) < nedges; limit++) {
	}

#if 0 /* Randomized rounding. */
	{
	int k;
	for (;;) {
		k = 0;
		for (i = 0; i < nedges; i++) {
			j = edge_list[i];
			if (((BITON (bbip -> fixed, j)) AND
			     (BITON (bbip -> value, j))) OR
			    (cip -> edge_size [j] EQ 2) OR
			    (drand48 () < 0.49*bbip -> x [j] + 0.5)) {
				temp_edges [k++] = j;
			}
		}

		ub_kruskal (temp_edges, k, used, bbip);
		if (--limit <= 0) break;
	}
	}
#endif

#if 1 /* Greedy local search. */
	{
	int e, k, org_limit;
	dist_t l, old_l, small_gap, curr_gap, fraction, quadratic;
	bool have_reset;

	have_reset = FALSE;
	org_limit  = limit;
	old_l      = INF_DISTANCE;
	k = nedges + 1;

	for (i = 0; i < nedges; i++) {
		e = edge_list [i];
		if (BITON (used, e)) continue;

		/* Force edge e to be chosen FIRST. */
		temp_edges [0] = e;

		/* Clear old used map */
		for (j = 0; j < nmasks; j++) {
			used [j] = 0;
		}

		l = ub_kruskal (temp_edges, k, used, bbip);

		/* If improved solution then replace temp_edges */
		if (l < old_l) {
			k = 1;
			for (j = 0; j < nedges; j++) {
				e = edge_list [j];
				if (BITON (used, e) OR
				    (cip -> edge_size [e] EQ 2)) {
					temp_edges [k++] = e;
				}
			}
			old_l = l;

			/* Compute "small" absolute gap value */
			small_gap = 0.0001 * fabs (ubip -> best_z);
			if (small_gap < 1.0) {
				small_gap = 1.0;
			}
		}

		/* If we are really close to optimum then restart	*/
		/* and intensify search.				*/

		curr_gap = l - ubip -> best_z;
		if ((NOT have_reset) AND (curr_gap <= small_gap)) {

#if 0
			tracef (" %% upper bound heuristic intensifying"
				" search (absolute gap is %.3f)\n",
				curr_gap);
#endif

			/* Let the new limit be double, plus an a*x**2	*/
			/* term whose coefficient is linearly dependent	*/
			/* on the gap.					*/
			fraction = (small_gap - curr_gap) / small_gap;
			quadratic = floor (fraction * (org_limit * org_limit));
			limit = 2 * org_limit + (int) quadratic;

			i = 0;
			have_reset = TRUE;
		}

		if (--limit <= 0) break;
	}
	}
#endif

#if 0 /* Round robin. */
	{
	int old_force_edge;
	static int force_edge = 0;
	if (force_edge EQ 0) {
		old_force_edge = nedges - 1;
	}
	else {
		old_force_edge = force_edge - 1;
	}

	for (;;) {
		if (NOT BITON (used, force_edge)) {

			/* Force edge force_edge to be chosen FIRST. */
			temp_edges [0] = force_edge;

			ub_kruskal (temp_edges, nedges + 1, used, bbip);
			--limit;
		}
		if (limit <= 0) break;
		if (force_edge EQ old_force_edge) break;
		++force_edge;
		if (force_edge EQ nedges) {
			force_edge = 0;
		}
	}
	}
#endif

#if 0 /* Very greedy. */
	for (i = 0; i < nedges; i++) {
		j = edge_list [i];
		if (BITON (used, j)) continue;

		/* Force edge j to be chosen FIRST. */
		temp_edges [0] = j;

		ub_kruskal (temp_edges, nedges + 1, used, bbip);

		if (--limit <= 0) break;
	}
#endif

#if 0 /* Less greedy. */
	{
	int iskip, nowskip;
	iskip = 1;
	nowskip = 0;
	limit *= 10;	/* !!!!!!!!! */

	for (i = 0; i < nedges; i++) {
		j = edge_list [i];
		if (BITON (used, j)) continue;
		if (++nowskip < iskip) continue;

		/* Force edge j to be chosen FIRST. */
		temp_edges [0] = j

		ub_kruskal (temp_edges, nedges + 1, used, bbip);

		nowskip = 0;
		iskip++;
		if (--limit <= 0) break;
	}
	}
#endif

	free ((char *) temp_edges);
	free ((char *) used);
}

/*
 * Use Kruskal's algorithm (extended for hypergraphs) to greedily construct
 * a tree from the given sequence of FSTs.  If this yields a better
 * solution than previously known, record it and update the upper bound.
 */

	static
	dist_t
ub_kruskal (

int *			edge_list,	/* IN - ordered list of edges */
int			nedges,		/* IN - number of edges in list */
bitmap_t *		used,		/* IN - marks FSTs used (ever) */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			j;
int			e;
int			nverts;
int *			ep1;
int *			ep2;
int *			ep3;
int *			vp1;
int *			vp2;
int *			vp3;
int *			vp4;
int *			tree_edges;
struct cinfo *		cip;
struct ubinfo *		ubip;
int *			roots;
bool *			mark;
int			components;
dist_t			length;
struct dsuf		sets;

	ubip		= bbip -> ubip;

	cip		= bbip -> cip;

	nverts		= cip -> num_verts;

	/* Initialize one disjoint subtree for each terminal. */
	dsuf_create (&sets, nverts);
	for (i = 0; i < nverts; i++) {
		dsuf_makeset (&sets, i);
	}

	tree_edges = NEWA (nverts - 1, int);

	mark = NEWA (nverts, bool);
	for (i = 0; i < nverts; i++) {
		mark [i] = FALSE;
	}
	roots = NEWA (nverts, int);

	/* Start the greedy Kruskal procedure... */
	components = nverts;
	length = 0.0;
	ep1 = tree_edges;
	ep2 = edge_list;
	ep3 = edge_list + nedges;
	while (components > 1) {
		if (ep2 >= ep3) {
			/* FSTs ran out before tree constructed! */
			fatal ("kruskal: Bug 1.");
		}
		e = *ep2++;
		vp3 = roots;
		vp1 = cip -> edge [e];
		vp2 = cip -> edge [e + 1];
		for (;;) {
			if (vp1 >= vp2) {
				/* No cycle!  Include e in solution! */
				*ep1++ = e;
				length += cip -> cost [e];
				SETBIT (used, e);
				/* Unite all subtrees joined, and clear	*/
				/* out the mark bits...			*/
				vp4 = roots;
				i = *vp4++;
				mark [i] = FALSE;
				while (vp4 < vp3) {
					j = *vp4++;
					mark [j] = FALSE;
					dsuf_unite (&sets, i, j);
					i = dsuf_find (&sets, i);
					--components;
				}
				break;
			}
			j = dsuf_find (&sets, *vp1++);
			if (mark [j]) {
				/* This FST would form a cycle. */
				while (roots < vp3) {
					mark [*--vp3] = FALSE;
				}
				break;
			}
			/* This FST terminal is in distinct subtree... */
			mark [j] = TRUE;
			*vp3++ = j;
		}
	}

	if (length < ubip -> best_z) {
		/* We have a better feasible integer solution! */
		for (i = 0; i < cip -> num_edge_masks; i++) {
			bbip -> smt [i] = 0;
		}
		while (ep1 > tree_edges) {
			e = *--ep1;
			SETBIT (bbip -> smt, e);
		}

		ubip -> best_z = length;
	}

	free ((char *) roots);
	free ((char *) mark);
	free ((char *) tree_edges);
	dsuf_destroy (&sets);

	return (length);
}