bb.c

生成直角Steiner树的程序包

		tracef (" %% Initial guess is x%d,"
			" Z0 = %-24.15g, Z1 = %-24.15g\n\n",
			best.var,
			best.z0,
			best.z1);
	}
#endif

	/* Now do the expensive part -- testing one or both branches	*/
	/* of good candidate variables.					*/

	new_limit = -1;
	limit = nfrac;

	num_failures = 0;

again:

	for (j = 0; j < limit; j++) {

		if (force_branch_flag) {
			if (best.var >= 0) goto get_out;

			/* Ignore this interrupt! */
			force_branch_flag = FALSE;
		}

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

		xi = x [i];

		z0 = nodep -> zlb [2 * i + 0];
		z1 = nodep -> zlb [2 * i + 1];

		if (((z0 > nodep -> z) AND (z0 < z1)) OR
		    ((z0 EQ z1) AND (xi <= 0.5))) {
			/* Check the Xi=0 branch, and then the Xi=1 branch. */
			fixed = eval_branch_var (bbip,
						 i,
						 0,	/* Xi=0, then Xi=1 */
						 &bsave,
						 test_2nd_val);
		}
		else {
			/* Check the Xi=1 branch, and then the Xi=0 branch. */
			fixed = eval_branch_var (bbip,
						 i,
						 1,	/* Xi=1, then Xi=0 */
						 &bsave,
						 test_2nd_val);
		}

		if (fixed) {
#if 1
			/* Special return code that says to try */
			/* re-solving the LP again.		*/
			destroy_LP_basis (&bsave);
			free ((char *) fvars);
			return (-1);
#elif 0
			goto start_all_over;
#else
			fvars [j] = -1;
			new_limit = j;
			limit = nfrac;
			num_failures = 0;
#endif
		}

		/* Test the current var to see if it is better than the	*/
		/* best seen so var.					*/

		cur_var_is_better = compare_branch_vars (bbip, i, &best);

		if (cur_var_is_better) {
			/* Establish a new threshold for testing 2nd branch. */
			test_2nd_val = best.test_2nd_val;

			z0 = nodep -> zlb [2 * i + 0];
			z1 = nodep -> zlb [2 * i + 1];

			z = z0;
			if (z1 < z) {
				z = z1;
			}

#if 1
			tracef (" %%   New best:  x%d, Z = %-24.15g\n",
				best.var, z);
#endif

			if (z >= bbip -> best_z) {
				/* Nice deal!  This variable	*/
				/* forces a cutoff on both	*/
				/* branches!  No need to look	*/
				/* at the rest...		*/
				new_limit = -1;
				break;
			}
			num_failures = 0;
		}
		else {
			++num_failures;
			if (num_failures >= failure_limit) {
#if 1
				tracef (" %% %d consecutive failures: giving up.\n",
					num_failures);
#endif
				goto get_out;
			}
		}
	}

	if (best.var < 0) {
		fatal ("carefully_choose_branching_variable: Bug 1.");
	}

	if (new_limit >= 0) {
		/* We have fixed some variables.  Must retest. */
		limit = new_limit;
		new_limit = -1;
		goto again;
	}

get_out:

#if 1
	tracef (" %% Best branch is x%d, Z0 = %-24.15g, Z1 = %-24.15g\n\n",
		best.var, best.z0, best.z1);
#endif

	destroy_LP_basis (&bsave);

	free ((char *) fvars);

	*node_z0 = best.z0;
	*node_z1 = best.z1;

	return (best.var);
}

/*
 * Sort the given list of fractional variables so that the best branching
 * variables appear early in the list with high probability.
 */

	static
	void
sort_branching_vars (

int *			fvars,	/* IN/OUT - fractional variables to sort */
int			n,	/* IN - number of fractional vars */
double *		rank	/* IN - heuristic rank of each variable */
)
{
int			i, i1, i2, j, k;

	/* Construct the heap via sift-downs, in O(n) time. */
	for (k = n >> 1; k >= 0; k--) {
		j = k;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is larger. */
				i1 = fvars [i];
				i2 = fvars [i + 1];
				if (rank [i2] > rank [i1]) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}
			i1 = fvars [j];
			i2 = fvars [i];
			if (rank [i2] < rank [i1]) {
				/* Largest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			fvars [j] = i2;
			fvars [i] = i1;
			j = i;
		}
	}

	/* Now do actual sorting.  Exchange first/last and sift down. */
	while (n > 1) {
		/* Largest is at fvars [0], swap with fvars [n-1],	*/
		/* thereby putting it into final position.		*/
		--n;
		i = fvars [0];
		fvars [0] = fvars [n];
		fvars [n] = i;

		/* Now restore the heap by sifting fvars [0] down. */
		j = 0;
		for (;;) {
			i = (j << 1) + 1;
			if (i + 1 < n) {
				/* Increment i (to right subchild of j) */
				/* if the right subchild is larger. */
				i1 = fvars [i];
				i2 = fvars [i + 1];
				if (rank [i2] > rank [i1]) {
					++i;
				}
			}
			if (i >= n) {
				/* Hit bottom of heap, sift-down is done. */
				break;
			}
			i1 = fvars [j];
			i2 = fvars [i];
			if (rank [i2] < rank [i1]) {
				/* Largest child is smaller.  Sift-	*/
				/* down is done. */
				break;
			}
			/* Sift down and continue. */
			fvars [j] = i2;
			fvars [i] = i1;
			j = i;
		}
	}
}

/*
 * This routine does the guts of carefully testing a single fractional
 * branching variable.  The caller indicates which direction should
 * be tested first.
 */

	static
	bool
eval_branch_var (

struct bbinfo *		bbip,		/* IN - branch-and-bound info */
int			var,		/* IN - variable to branch */
int			dir1,		/* IN - first branch direction */
struct basis_save *	basp,		/* IN - basis to restore when done */
double			test_2nd_val	/* IN - test 2nd if 1st is > this */
)
{
int			i;
int			nedges;
int			dir2;
struct cinfo *		cip;
LP_t *			lp;
struct bbnode *		nodep;
bool			found;
bool			fixed;
double *		x;
double			z;

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

	nedges = cip -> num_edges;

	x = NEWA (nedges, double);

	dir2 = 1 - dir1;

	fixed = FALSE;

	/* Try the first branch direction... */
	z = try_branch (lp, var + 1, dir1, x, DBL_MAX, basp);

#if CPLEX
	z = ldexp (z, bbip -> lpmem -> obj_scale);
#endif

	/* Check for a better integer feasible solution... */
	found = check_for_better_IFS (x, bbip, &z);

	/* Update per-variable lower bounds. */
	i = 2 * var + dir1;
	if (z > nodep -> zlb [i]) {
		nodep -> zlb [i] = z;
	}
	else {
		z = nodep -> zlb [i];
	}

#if 1
	tracef (" %%%s\tx%d = %d,\tZ%d = %-24.15g\n",
		found ? " !!!" : "",
		var, dir1, dir1, z);
#endif

	/* Try finding a good heuristic solution on the branched solution. */
	compute_heuristic_upper_bound (x, bbip);

#if 1
	if (z >= bbip -> best_z + 1.0e-8 * fabs (bbip -> best_z)) {
		/* Cutoff or infeasible.  Var must be fixed to	*/
		/* other direction.  Reoptimize and get new	*/
		/* basis.					*/
		SETBIT (bbip -> fixed, var);
		SETBIT (bbip -> node -> fixed, var);
		if (dir1) {
			CLRBIT (bbip -> value, var);
			CLRBIT (bbip -> node -> value, var);
		}
		else {
			SETBIT (bbip -> value, var);
			SETBIT (bbip -> node -> value, var);
		}
		change_var_bounds (lp,
				   var,
				   (double) dir2,
				   (double) dir2);

		nodep -> cpiter = -1;	/* force re-solve of LP */

		solve_LP_over_constraint_pool (bbip);
		lp = bbip -> lp;
		save_LP_basis (lp, basp);

		/* Try finding a good heuristic solution on the	*/
		/* new fixed solution...			*/
		compute_heuristic_upper_bound (bbip -> node -> x, bbip);

		/* The variable has been fixed! */
		fixed = TRUE;
		goto all_done;
	}
#endif

	if (z <= test_2nd_val)
	{
		/* No need to test the second branch. */
		goto all_done;
	}

	/* Try the second branch direction... */
	z = try_branch (lp, var + 1, dir2, x, DBL_MAX, basp);

#if CPLEX
	z = ldexp (z, bbip -> lpmem -> obj_scale);
#endif

	/* Check for better integer feasible solution... */
	found = check_for_better_IFS (x, bbip, &z);

	/* Update per-variable lower bounds. */
	i = 2 * var + dir2;
	if (z > nodep -> zlb [i]) {
		nodep -> zlb [i] = z;
	}
	else {
		z = nodep -> zlb [i];
	}

#if 1
	tracef (" %%%s\tx%d = %d,\tZ%d = %-24.15g\n",
		found ? " !!!" : "",
		var, dir2, dir2, z);
#endif

	/* Try finding a good heuristic solution on the branched solution. */
	compute_heuristic_upper_bound (x, bbip);

#if 1
	if (z >= bbip -> best_z + 1.0e-8 * fabs (bbip -> best_z)) {
		/* Cutoff or infeasible.  Var must be fixed to	*/
		/* other direction.  Reoptimize and get new	*/
		/* basis.					*/
		SETBIT (bbip -> fixed, var);
		SETBIT (bbip -> node -> fixed, var);
		if (dir2) {
			CLRBIT (bbip -> value, var);
			CLRBIT (bbip -> node -> value, var);
		}
		else {
			SETBIT (bbip -> value, var);
			SETBIT (bbip -> node -> value, var);
		}
		change_var_bounds (lp,
				   var,
				   (double) dir1,
				   (double) dir1);

		nodep -> cpiter = -1;	/* force re-solve of LP */

		solve_LP_over_constraint_pool (bbip);
		lp = bbip -> lp;
		save_LP_basis (lp, basp);

		/* Try finding a good heuristic solution on the	*/
		/* new fixed solution...			*/
		compute_heuristic_upper_bound (bbip -> node -> x, bbip);

		/* The variable has been fixed! */
		fixed = TRUE;
	}
#endif

all_done:

	free ((char *) x);

	return (fixed);
}

/*
 * See if one candidate branch variable is better than another.
 * We implement various policies here.
 */

	static
	bool
compare_branch_vars (

struct bbinfo *		bbip,		/* IN - branch and bound info */
int			i1,		/* IN - first branch var */
struct bvar *		bvp		/* IN/OUT - current best branch var */
)
{
struct bbnode *		nodep;
bool			cur_var_is_better;
double			z;
double			z0;
double			z1;
double			zmin;
double			zmax;
double			best_z0;
double			best_z1;
double			best_zmin;
double			best_zmax;
double			ub;
double			gap;
double			delta;
double			prod;
double			best_prod;
double			test2;

#define TOLERANCE	1.0e-10

	nodep	= bbip -> node;

	ub	= bbip -> best_z;
	z	= nodep -> z;

	if (i1 < 0) {
		return (FALSE);
	}

	z0	= nodep -> zlb [2 * i1 + 0];
	z1	= nodep -> zlb [2 * i1 + 1];

	if (z0 < z1) {
		zmin = z0;
		zmax = z1;
	}
	else {
		zmin = z1;
		zmax = z0;
	}

	if (bvp -> var < 0) {
		/* New branch var is better.  Still have to provide a	*/
		/* test_2nd_val threshold, which is policy specific.	*/

		switch (Branch_Var_Policy) {
		case 0:
			test2 = zmin;
			break;

		case 1:
			test2 = zmin - TOLERANCE * fabs (zmin);
			break;

		case 2:
			gap = fabs (ub - z);
			if (gap <= 0.0) {
				fatal ("compare_branch_vars: Bug 1.");
			}
			prod = fabs ((z0 - z) * (z1 - z));
			test2 = z + prod / gap;
			break;

		default:
			fatal ("compare_branch_vars: Bug 2.");
			break;
		}

		bvp -> var		= i1;
		bvp -> z0		= z0;
		bvp -> z1		= z1;
		bvp -> test_2nd_val	= test2;

		return (TRUE);
	}

	best_z0	= bvp -> z0;
	best_z1	= bvp -> z1;

	if (best_z0 < best_z1) {
		best_zmin = best_z0;
		best_zmax = best_z1;
	}
	else {
		best_zmin = best_z1;
		best_zmax = best_z0;
	}

	switch (Branch_Var_Policy) {
	case 0:
		/* Naive max of mins. */
		cur_var_is_better = FALSE;
		if (zmin > best_zmin) {
			cur_var_is_better = TRUE;
			test2 = zmin;
		}
		break;

	case 1:
		/* Smarter lexicographic max of mins.  If the mins	*/
		/* are "about equal", use the max values to decide.	*/
fuzzy_lexical:
		cur_var_is_better = FALSE;

		delta = TOLERANCE * fabs (best_zmin);
		if ((zmin - best_zmin) > delta) {
			cur_var_is_better = TRUE;
		}
		else if ((fabs (zmin - best_zmin) <= delta) AND
			 (zmax > best_zmax)) {
			cur_var_is_better = TRUE;
		}
		test2 = zmin - TOLERANCE * fabs (zmin);
		break;

	case 2:
		/* Product of improvements.  Uses method 1 to break	*/
		/* close ties.						*/
		prod = fabs ((z0 - z) * (z1 - z));
		best_prod = fabs ((best_z0 - z) * (best_z1 - z));

		/* Compute tolerance factor that is a fraction of the	*/
		/* current gap.						*/
		gap = fabs (ub - z);
		if (gap <= 0.0) {
			fatal ("compare_branch_vars: Bug 3.");
		}
		delta = 1.0E-5 * gap;
		if (fabs (prod - best_prod) <= delta) {
			/* Products are nearly equal.  Use fuzzy	*/
			/* lexicographic max of mins.			*/
			goto fuzzy_lexical;
		}
		cur_var_is_better = FALSE;
		if (prod > best_prod) {
			cur_var_is_better = TRUE;
			test2 = z + prod / gap;
		}
		break;

	default:
		fatal ("compare_branch_vars: Bug 4.");
		break;
	}

	if (cur_var_is_better) {

		bvp -> var		= i1;
		bvp -> z0		= z0;
		bvp -> z1		= z1;
		bvp -> test_2nd_val	= test2;
	}

	return (cur_var_is_better);

#undef TOLERANCE
}

/*
 * This routine checks to see if the result of doig a "test-branch" on
 * a variable JUST HAPPENS to result in an integer feasible solution
 * that is the best so far.  Although this would be total serendipity,
 * we do it anyway.  The check takes virtually no time because we almost
 * always locate a fractional variable within the first few probes.