bb.c

生成直角Steiner树的程序包

			/* Finished the root node... */
			statp -> root_z = node -> z;
			statp -> root_lps = statp -> num_lps;

			/* Slack rows should have already been deleted... */
			statp -> cs_root.num_prows = cpool -> nrows;
			statp -> cs_root.num_lprows =
				GET_LP_NUM_ROWS (bbip -> lp);
			statp -> cs_root.num_pnz = cpool -> num_nz;
			statp -> cs_root.num_lpnz =
				GET_LP_NUM_NZ (bbip -> lp);
			statp -> root_opt = node -> optimal;

			t1 = get_cpu_time ();
			statp -> root_time = t1 - bbip -> t0;

			if ((status EQ LB_FRACTIONAL) AND Print_Root_LP) {
				print_root_lp (bbip);
			}

			if (Check_Root_Constraints) {
				check_root_constraints (bbip);
			}
		}

		node_to_free = NULL;	/* Default is no node to free... */

		switch (status) {
		case LB_INFEASIBLE:
			/* Node is fathomed! */
			trace_node (bbip, ' ', "infeasible");
			node_to_free = node;
			break;

		case LB_CUTOFF:
			/* Node is fathomed! */
			trace_node (bbip, ' ', "cutoff");
			node_to_free = node;
			break;

		case LB_INTEGRAL:
			if (node -> z >= bbip -> best_z) {
				trace_node (bbip, ' ', "cutoff");
				break;
			}
			/* We have a new best feasible integer solution! */
			for (i = 0; i < nmasks; i++) {
				smt [i] = 0;
			}
			x = node -> x;
			for (i = 0; i < nedges; i++) {
				if (x [i] >= 0.5) {
					SETBIT (smt, i);
				}
			}
			new_upper_bound (node -> z, bbip);

			trace_node (bbip, '*', NULL);

			node_to_free = node;
			break;

		case LB_FRACTIONAL:
			j = choose_branching_variable (bbip, &z0, &z1);
			if (j < 0) {
				/* At least one variable was fixed due	*/
				/* to cutoff or infeasibility.  It is	*/
				/* possible that the entire node is now	*/
				/* cutoff...				*/
				if (node -> z >= bbip -> best_z) {
					trace_node (bbip, ' ', "cutoff");
					node_to_free = node;
					break;
				}
				goto suspend;
			}

			/* Create two nodes... */
			if (z0 < z1) {
				add_bbnode (bbip, j, 0, z0);
				add_bbnode (bbip, j, 1, z1);
			}
			else if (z1 < z0) {
				add_bbnode (bbip, j, 1, z1);
				add_bbnode (bbip, j, 0, z0);
			}
			else if (UP_FIRST) {	/* To break ties... */
				add_bbnode (bbip, j, 0, z0);
				add_bbnode (bbip, j, 1, z1);
			}
			else {
				add_bbnode (bbip, j, 1, z1);
				add_bbnode (bbip, j, 0, z0);
			}
			trace_node (bbip, ' ', NULL);

			/* This node is done (became 2 children), free it. */
			node_to_free = node;
			break;

		case LB_PREEMPTED:
suspend:
			tracef ("%% suspending node %d at %24.20f\n",
				node -> num,
				UNSCALE (node -> z, &(cip -> scale)));
			/* This node is no longer the best.  Put it	*/
			/* back into the heap and get another one...	*/
			node2 = bbtree -> first;
			if (node2 NE NULL) {
				node2 -> prev = node;
			}
			node -> next = node2;
			node -> prev = NULL;
			bbtree -> first = node;

			bbheap_insert (node, bbtree, BEST_NODE_HEAP);
			bbheap_insert (node, bbtree, WORST_NODE_HEAP);

			/* Deactivating this node -- remember the basis */
			save_node_basis (node, bbip);

			/* Do NOT free this node! */
			break;
		}

		/* If there is a node to free, do so now... */
		if (node_to_free NE NULL) {
			/* Free up saved basis info and decrement	*/
			/* constraint reference counts before freeing.	*/
			destroy_node_basis (node_to_free, bbip);

			node_to_free -> next = bbtree -> free;
			bbtree -> free = node_to_free;
		}
	}

	statp -> cs_final.num_prows	= cpool -> nrows;
	statp -> cs_final.num_lprows	= GET_LP_NUM_ROWS (bbip -> lp);
	statp -> cs_final.num_pnz	= cpool -> num_nz;
	statp -> cs_final.num_lpnz	= GET_LP_NUM_NZ (bbip -> lp);

	/* New lower bound! */
	new_lower_bound (bbip -> best_z, bbip);

	t1 = get_cpu_time ();
	statp -> p2time = t1 - bbip -> t0;
	statp -> z = bbip -> best_z;

	/* Restore normal handling of SIGTERM... */
	(void) signal (SIGTERM, save_sigterm);

#if CPLEX
	free ((char *) b_bd);
	free ((char *) b_lu);
	free ((char *) b_index);
#endif

	free ((char *) delta);
	free ((char *) value);
	free ((char *) fixed);

#if CPLEX
	CPXsetdblparam (cplex_env, CPX_PARAM_OBJULIM, save_objlim);
#endif

	/* Return length of final tree... */

	return ((dist_t) bbip -> best_z);
}

/*
 * A handler for the SIGTERM signal.  When we receive this signal we
 * stop generating constraints for the current node and branch instead.
 */

	static
	RETSIGTYPE
force_branch_signal_handler (

int		sig		/* IN - signal being handled */
)
{
	/* Stop generating constraints and force a branch instead. */
	force_branch_flag = TRUE;

	/* Prepare to handle the signal again... */
	(void) signal (SIGTERM, force_branch_signal_handler);
}

/*
 * This routine selects the next node to process from the given
 * branch-and-bound tree.  This is where we implement the specific
 * search policy.
 */

	static
	struct bbnode *
select_next_node (

struct bbtree *		tp	/* IN - branch-and-bound tree */
)
{
struct bbnode *		p;

	if (tp -> first EQ NULL) {
		/* No more nodes! */
		return (NULL);
	}

	switch (tp -> node_policy) {
	case NN_DEPTH_FIRST:
		/* Get node created most recently... */
		p = tp -> first;
		break;

	case NN_BEST_NODE:
		/* Get node with lowest objective function value... */
		p = tp -> heap [BEST_NODE_HEAP].array [0];
		break;

	default:
		fatal ("select_next_node: Bug 1.");
	}

	delete_node_from_bbtree (p, tp);

	return (p);
}

/*
 * This routine traces the result for a given node.
 */

	static
	void
trace_node (

struct bbinfo *		bbip,	/* IN - the branch-and-bound info */
char			c1,	/* IN - space or * char */
char *			msg	/* IN - message OR NULL */
)
{
struct bbnode *		p;
struct bbtree *		tp;
struct bbheap *		hp;
double			lowz;
char *			p1;
char *			p2;
char			c2;
char			buf1 [32];
char			buf2 [32];
char			buf3 [32];
char			buf4 [32];
char			buf5 [32];
char			buf6 [32];
char			line [136];

	p	= bbip -> node;
	tp	= bbip -> bbtree;

	if (msg EQ NULL) {
		(void) sprintf (buf1, "%14.4f", p -> z);
	}
	else {
		(void) sprintf (buf1, "%14s", msg);
	}
	if (bbip -> best_z EQ DBL_MAX) {
		buf2 [0] = '\0';
	}
	else {
		(void) sprintf (buf2, "%14.4f", bbip -> best_z);
	}
	hp = &(tp -> heap [BEST_NODE_HEAP]);
	if (hp -> nheap > 0) {
		(void) sprintf (buf3,
				"%14.4f",
				hp -> array [0] -> z);
	}
	else {
		buf3 [0] = '\0';
	}
	if (p -> var < 0) {
		buf4 [0] = '\0';
		c2 = ' ';
	}
	else {
		(void) sprintf (buf4, "x%d", p -> var);
		c2 = (p -> dir EQ 0) ? 'D' : 'U';
	}
	if (p -> parent < 0) {
		buf5 [0] = '\0';
	}
	else {
		(void) sprintf (buf5, "%d", p -> parent);
	}
	if (p -> depth <= 0) {
		buf6 [0] = '\0';
	}
	else {
		(void) sprintf (buf6, "%d", p -> depth);
	}

	(void) sprintf (line,
			"%c%6d%6d%14s%14s%14s%6s %c%6s%6s",
			c1,		/* space or * */
			p -> num,	/* node number */
			hp -> nheap,	/* nodes left */
			buf1,		/* objective/cutoff/infeas */
			buf2,		/* best integer soln */
			buf3,		/* best node */
			buf4,		/* variable name */
			c2,		/* branch direction */
			buf5,		/* parent node number */
			buf6);		/* node depth */

	p1 = line;
	for (p2 = p1; *p2 NE '\0'; p2++) {
	}
	while ((p2 > p1) AND (p2 [-1] EQ ' ')) {
		--p2;
	}
	*p2 = '\0';
	(void) tracef ("  %% %s\n", line);
}

/*
 * This routine plots the LP relaxation solution for the root node.  We
 * do this whenever -r is specified and we get an optimal LP relaxation
 * solution that is fractional.
 */

	static
	void
print_root_lp (

struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
struct cinfo *		cip;
struct bbnode *		nodep;
struct bbstats *	statp;
char *			descr;
char			buf1 [32];
char			buf2 [32];
char			title [256];

	cip	= bbip -> cip;
	nodep	= bbip -> node;
	statp	= bbip -> statp;

	convert_cpu_time (statp -> root_time, buf1);
	dist_to_string (buf2, nodep -> z, &(cip -> scale));

	descr = "Steiner Minimal Tree";
	if ((cip -> description NE NULL) AND
	    (cip -> description [0] NE '\0')) {
		descr = cip -> description;
	}
	sprintf (title,
		 "%s:  %lu points,  Root Node, LP %d,  length = %s,  %s seconds",
		 descr,
		 cip -> num_verts,
		 nodep -> iter,
		 buf2,
		 buf1);

	plot_lp_solution (title, bbip -> node -> x, cip, BIG_PLOT);
}

/*
 * This routine chooses the next variable to branch on at this
 * node.
 */

	static
	int
choose_branching_variable (

struct bbinfo *		bbip,		/* IN - branch-and-bound info */
double *		z0,		/* OUT - value to give Xi=0 node */
double *		z1		/* OUT - value to give Xi=1 node */
)
{
int			i;
int			nedges;
int			best_var;
struct cinfo *		cip;
struct bbnode *		nodep;
double *		x;
bitmap_t *		fset_mask;
double			xi;
double			max_infeas;
double			infeas;

	if (Choose_Branch_Vars_Carefully) {
		/* Do it very carefully! */
		return (carefully_choose_branching_variable (bbip, z0, z1));
	}

	/* There are LOTS of things that we MIGHT do here in	*/
	/* the FUTURE!  For right now, just take the variable	*/
	/* that is closest to 1/2 -- take larger variables in	*/
	/* the event of a tie...				*/

	cip		= bbip -> cip;
	nodep		= bbip -> node;
	x		= nodep -> x;
	fset_mask	= bbip -> fset_mask;

	nedges	= cip -> num_edges;

	best_var	= -1;
	max_infeas	= 0.0;

	for (i = nedges - 1; i >= 0; i--) {
		if (NOT BITON (fset_mask, i)) continue;
		xi = x [i];
		if (xi <= FUZZ) continue;
		if (xi + FUZZ >= 1.0) continue;
		infeas = 1.0 - xi;
		if (xi < infeas) {
			infeas = xi;
		}
		if (infeas > max_infeas) {
			best_var = i;
			max_infeas = infeas;
		}
	}
	if (best_var < 0) {
		fatal ("choose_branching_variable: Bug 1.");
	}

	/* Give both nodes the same value... */
	*z0	= nodep -> z;
	*z1	= nodep -> z;

	return (best_var);
}

/*
 * This routine does a very careful job of choosing the next variable
 * to branch on.  For each fractional variable Xi, we solve the LP
 * first with Xi=0 and then Xi=1, yielding objective values Zi0 and Zi1
 * correspondingly.  We then choose variable Xi for which the value
 * min(Zi0, Zi1) is MAXIMIZED.  This provides us with the best possible
 * "one-branch/no-cuts" improvement in the lower bound.
 */

	static
	int
carefully_choose_branching_variable (

struct bbinfo *		bbip,		/* IN - branch-and-bound info */
double *		node_z0,	/* OUT - value for Xi=0 node */
double *		node_z1		/* OUT - value for Xi=1 node */
)
{
int			i;
int			j;
int			n;
int			nedges;
int			nfrac;
int			limit;
int			new_limit;
int			logn;
int			failure_limit;
int			num_failures;
struct cinfo *		cip;
struct bbnode *		nodep;
double *		x;
double *		rank;
bitmap_t *		edge_mask;
int *			fvars;
bool			fixed;
bool			cur_var_is_better;
double			xi;
double			z0;
double			z1;
double			z;
double			delta;
double			test_2nd_val;
double			num;
double			den;
struct bvar		best;
struct basis_save	bsave;

	cip		= bbip -> cip;
	nodep		= bbip -> node;
	x		= nodep -> x;
	edge_mask	= bbip -> fset_mask;

	nedges	= cip -> num_edges;

	fvars = NEWA (nedges, int);

start_all_over:

	nfrac = 0;
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		xi = x [i];
		if (xi <= FUZZ) continue;
		if (xi + FUZZ >= 1.0) continue;
		fvars [nfrac++] = i;
	}
#if 1
	tracef ("\n %% Carefully choosing branching variable, nfrac = %d\n",
		nfrac);
#endif

	/* Compute final heuristic ranking of the candidate branch	*/
	/* variables.  We multiply the node's "bheur" values by a	*/
	/* "complexity factor" that is 1 for variables sitting at 0.5,	*/
	/* and increases as the variable's "rational" value becomes	*/
	/* more complicated.  Thus, we prefer vars stuck at 1/2, and	*/
	/* then prefer vars stuck at 1/3 or 2/3, vars stuck at 1/4 or	*/
	/* 3/4, etc.							*/

	rank = NEWA (nedges, double);
	memcpy (rank, nodep -> bheur, nedges * sizeof (double));

	for (j = 0; j < nfrac; j++) {
		i = fvars [j];

		/* Compute closest rational approximation. */

		(void) cra (x [i], &num, &den);

		/* Factor is 1 + the denominator's "distance" from 1/2	*/
		/* + the numerator's "distance" from 1/2.		*/

		rank [i] *= ((den - 1.0) + fabs (num - 0.5 * den));
	}

	/* Sort the fractional variables so that good branch choices	*/
	/* appear first with high probability.				*/

	sort_branching_vars (fvars, nfrac, rank);

#if 0
	tracef (" %% Branch Variable Info:\n");
	for (j = 0; j < nfrac; j++) {
		i = fvars [j];
		tracef (" %% %d:\tx%d\t= %.6f,\tbheur = %g,\trank = %g\n",
			j, i, x [i], nodep -> bheur [i], rank [i]);
	}
#endif

	free ((char *) rank);

	/* Snapshot the current basis so that we can quickly	*/
	/* get back to it each time...				*/
	save_LP_basis (bbip -> lp, &bsave);

	/* Compute the non-improvement limit.  When we have tested this	*/
	/* many consecutive variables without finding a better choice,	*/
	/* we punt.  We use 2 * log(N), where N is the number of	*/
	/* fractional vars, and log(x) is the floor of the base-2 log.	*/

	failure_limit = nfrac;
	if (nfrac >= 20) {
		logn = 0;
		for (n = nfrac; n > 1; n >>= 1) {
			++logn;
		}
		failure_limit = 2 * Check_Branch_Vars_Thoroughly * logn;
	}

	best.var	= -1;
	best.z0		= nodep -> z;
	best.z1		= nodep -> z;

	test_2nd_val	= -DBL_MAX;

	/* Do a quick scan without forcing anything to determine the	*/
	/* best initial choice of branch variable.			*/

	for (j = 0; j < nfrac; j++) {
		i = fvars [j];
		if (i < 0) continue;	/* var was fixed! */

		cur_var_is_better = compare_branch_vars (bbip, i, &best);
		if (cur_var_is_better) {
			test_2nd_val = best.test_2nd_val;
		}
	}

#if 1
	if (best.var >= 0) {