bb.c

生成直角Steiner树的程序包

 * And besides, it would be a shame to NOT notice something this important
 * -- especially if we don't have ANY upper bound yet!
 */

	bool
check_for_better_IFS (

double *		x,		/* IN - solution to test */
struct bbinfo *		bbip,		/* IN - branch and bound info */
double *		true_z		/* OUT - true Z value, if integral */
)
{
int			i;
int			nedges;
int			nmasks;
struct cinfo *		cip;
double			z;
double			real_z;
int			num_frac;

	cip = bbip -> cip;
	nedges	= cip -> num_edges;

	/* The special try_branch code for lp_solve can leave us with	*/
	/* an LP solution that isn't primal feasible.  In fact, it may	*/
	/* not even satisfy the variable bounds!  Detect this case	*/
	/* right up front.						*/

	for (i = 0; i < nedges; i++) {
		if (x [i] < -FUZZ) return (FALSE);
		if (x [i] > 1.0 + FUZZ) return (FALSE);
	}

	if (NOT integer_feasible_solution (x,
					   bbip -> tmap,
					   bbip -> fset_mask,
					   cip,
					   &num_frac)) {
		/* Not integer feasible -- get out. */
		return (FALSE);
	}

	/* We literally stumbled across an Integer Feasible Solution!	*/

	/* Re-calculate the final objective function (in order to	*/
	/* eliminate some of the LP solver's numerical errors)...	*/

	z = 0.0;
	for (i = 0; i < nedges; i++) {
		if (x [i] + FUZZ < 1.0) continue;
		z += cip -> cost [i];
	}

	/* Give caller the correct Z value... */
	*true_z = z;

	/* Damn Intel FPU is keeping 80 bits for Z... */
	/* Damn compilers keep getting smarter too!  ;^> */
	store_double (&real_z, z);

	if (real_z >= bbip -> best_z) {
		/* No better than what we have... */
		return (FALSE);
	}

	/* We have a new best integer feasible solution!  Record it!	*/
	nmasks = cip -> num_edge_masks;
	for (i = 0; i < nmasks; i++) {
		bbip -> smt [i] = 0;
	}
	for (i = 0; i < nedges; i++) {
		if (x [i] >= 0.5) {
			SETBIT (bbip -> smt, i);
		}
	}

	new_upper_bound (real_z, bbip);

	return (TRUE);
}

/*
 * This routine saves the current basis of the given LP.
 */

#ifdef CPLEX

	static
	void
save_LP_basis (

LP_t *			lp,		/* IN - LP to save basis for */
struct basis_save *	basp		/* OUT - saved basis info */
)
{
int		rows;
int		cols;

	rows = _MYCPX_getnumrows (lp);
	cols = _MYCPX_getnumcols (lp);

	basp -> cstat = NEWA (cols, int);
	basp -> rstat = NEWA (rows, int);

	if (_MYCPX_getbase (lp, basp -> cstat, basp -> rstat) NE 0) {
		fatal ("save_LP_basis: Bug 1.");
	}
}

/*
 * Destroy the saved basis info...
 */

	static
	void
destroy_LP_basis (

struct basis_save *	basp		/* IN - basis info to free up */
)
{
#if CPLEX >= 40
	free ((char *) (basp -> rstat));
	free ((char *) (basp -> cstat));
#endif
}

#endif

/*
 * This routine tries the given branch by solving the LP.  It
 * returns the resulting objective value, or INF_DISTANCE if something
 * goes wrong (like infeasible).
 */

#ifdef CPLEX

	static
	double
try_branch (

LP_t *			lp,		/* IN - LP to re-optimize */
int			var,		/* IN - variable to try branching */
int			dir,		/* IN - branch direction, 0 or 1 */
double *		x,		/* OUT - LP solution obtained */
double			ival,		/* IN - value to give if infeasible */
struct basis_save *	basp		/* IN - basis to restore when done */
)
{
int		status;
double		z;
int		b_index [2];
char		b_lu [2];
double		b_bd [2];

	--var;		/* vars are zero-origined in CPLEX... */

	b_index [0] = var;	b_lu [0] = 'L';
	b_index [1] = var;	b_lu [1] = 'U';
	if (dir EQ 0) {
		b_bd [0] = 0.0;
		b_bd [1] = 0.0;
	}
	else {
		b_bd [0] = 1.0;
		b_bd [1] = 1.0;
	}
	if (_MYCPX_chgbds (lp, 2, b_index, b_lu, b_bd) NE 0) {
		fatal ("try_branch: Bug 1.");
	}

	/* Solve the current LP instance... */
	status = _MYCPX_dualopt (lp);
	if (status NE 0) {
		tracef (" %%  WARNING dualopt: status = %d\n", status);
	}

	/* Get current LP solution... */
	if (_MYCPX_solution (lp, &status, &z, x, NULL, NULL, NULL) NE 0) {
		fatal ("try_branch: Bug 2.");
	}

	/* Determine type of LP result... */
	switch (status) {
	case CPX_OPTIMAL:
	case CPX_OPTIMAL_INFEAS:
		break;

	case CPX_INFEASIBLE:
	case CPX_UNBOUNDED:
			/* (CPLEX 3.0 sometimes gives us infeasible!) */
	case CPX_OBJ_LIM:	/* Objective limit exceeded in Phase II. */
		z = ival;
		break;

	default:
		tracef (" %% Status = %d\n", status);
		_MYCPX_lpwrite (lp, "core.lp");
		fatal ("try_branch: Bug 2.");
		break;
	}

	b_bd [0] = 0.0;
	b_bd [1] = 1.0;
	if (_MYCPX_chgbds (lp, 2, b_index, b_lu, b_bd) NE 0) {
		fatal ("try_branch: Bug 3.");
	}

	/* Restore the basis... */
	status = _MYCPX_copybase (lp, basp -> cstat, basp -> rstat);
	if (status NE 0) {
		fprintf (stderr, "try_branch: status = %d\n", status);
		fatal ("try_branch: Bug 4.");
	}

	return (z);
}

#endif

/*
 * This routine computes the lower-bound for the current node, which
 * consists of solving the LP and generating violated constraints
 * until either:
 *
 *	- LP becomes infeasible
 *	- LP objective meets or exceeds cutoff value
 *	- LP solution is integral
 *	- separation finds no more violated constraints
 */

	static
	int
compute_good_lower_bound (

struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			j;
int			n;
int			num_const;
int			status;
bitmap_t *		tmap;
bitmap_t *		fset_mask;
struct cinfo *		cip;
LP_t *			lp;
struct bbnode *		nodep;
double *		x;
double			z;
struct constraint *	cp;
struct constraint *	tmp;
int			iteration;
int			fix_status;
int			num_fractional;
cpu_time_t		Tlp;
cpu_time_t *		Tp;
bool			is_int;
char			tbuf1 [16];
char			tbuf2 [256];
char			time_str [256];
char			title [256];
cpu_time_t		Tn [20];

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

	Tp = &Tn [0];

	*Tp++ = get_cpu_time ();

	iteration = 1;
	num_const = 0;

	for (;;) {
		status = solve_LP_over_constraint_pool (bbip);
		z = nodep -> z;

		Tlp = get_cpu_time ();

		++(bbip -> node -> iter);
		++(bbip -> statp -> num_lps);

#if 0
		/* Display LP solution vector in machine-readable form... */
		for (i = 0; i < cip -> num_edges; i++) {
			tracef (" %% %08lx %08lx\n",
				((bitmap_t *) &x [i]) [0],
				((bitmap_t *) &x [i]) [1]);
		}
#endif

		if (Check_Root_Constraints AND
		    (bbip -> node -> depth EQ 0)) {
			fwrite (x,
				1,
				cip -> num_edges * sizeof (*x),
				rcfile);
		}

#if 1
		convert_cpu_time (Tlp - *--Tp, time_str);
		while (Tp > &Tn [0]) {
			--Tp;
			convert_cpu_time (Tp [1] - Tp [0], tbuf1);
			(void) sprintf (tbuf2, "%s/%s", tbuf1, time_str);
			strcpy (time_str, tbuf2);
		}
		(void) sprintf (title,
				"Node %d LP %d Solution, length = %f, %s %d",
				bbip -> node -> num, bbip -> node -> iter,
				z, time_str, num_const);
#if 0
		plot_lp_solution (title, x, cip, BIG_PLOT);
#else
		(void) tracef ("  %% %s\n", title);
#endif
#endif

		switch (status) {
		case BBLP_OPTIMAL:
			if (z >= bbip -> best_z) {
				nodep -> z = bbip -> best_z;
				return (LB_CUTOFF);
			}
			break;

		case BBLP_CUTOFF:
			nodep -> z = bbip -> best_z;
			return (LB_CUTOFF);

		case BBLP_INFEASIBLE:
			nodep -> z = bbip -> best_z;
			return (LB_INFEASIBLE);

		default:
			tracef ("%% solve status = %d\n", status);
			fatal ("compute_good_lower_bound: Bug 3.");
		}

#ifdef CPLEX
		/* Now get rid of any rows that have become	*/
		/* slack.  (We don't lose these constraints:	*/
		/* they're still sitting around in the		*/
		/* constraint pool.)				*/
		delete_slack_rows_from_LP (bbip);
#endif

		/* Solution is feasible, check for integer-feasible... */
		is_int = integer_feasible_solution (x,
						    tmap,
						    fset_mask,
						    cip,
						    &num_fractional);

		tracef (" %% %d fractional variables\n", num_fractional);

		if (is_int) {
			/* All vars are either 0 or 1 and the	*/
			/* solution is connected -- we have a	*/
			/* Steiner tree!			*/

			/* Re-calculate the final objective	*/
			/* function to eliminate numerical	*/
			/* errors in the value of Z...		*/
			z = 0.0;
			for (i = 0; i < cip -> num_edges; i++) {
				if (x [i] + FUZZ < 1.0) continue;
				z += cip -> cost [i];
			}
			nodep -> z = z;
			if (z >= bbip -> best_z) {
				/* probably a repeat performance... */
				nodep -> z = bbip -> best_z;
				return (LB_CUTOFF);
			}
			bbip -> node -> optimal = TRUE;
			return (LB_INTEGRAL);
		}

		/* Check to see if this node's objective value	*/
		/* is now high enough to be preempted...	*/
		if (nodep -> z > bbip -> preempt_z) {
			/* Node's value is no longer the lowest...	*/
			/* Preempt this one in favor of another.	*/
			return (LB_PREEMPTED);
		}

		/* Perhaps we have a new lower bound? */
		new_lower_bound (z, bbip);

		Tp = &Tn [0];
		*Tp++ = get_cpu_time ();

		compute_heuristic_upper_bound (x, bbip);

		/* If we have improved the upper bound, it is possible	*/
		/* that this node can now be cutoff...			*/
		if (nodep -> z >= bbip -> best_z) {
			nodep -> z = bbip -> best_z;
			return (LB_CUTOFF);
		}

		/* Try to fix some variables using reduced costs... */
		fix_status = reduced_cost_var_fixing (bbip);
		if (fix_status EQ VFIX_INFEASIBLE) {
			nodep -> z = bbip -> best_z;
			return (LB_INFEASIBLE);
		}
		if (fix_status EQ VFIX_FIXED_FRACTIONAL) {
			continue;
		}

		if (force_branch_flag) {
			/* User kicked us!  Stop separating and branch! */
			force_branch_flag = FALSE;
			break;
		}

		/* Apply all separation algorithms to solution... */
		cp = do_separations (bbip, &Tp);

		if (cp EQ NULL) {
			/* No more violated constraints found! */
			break;
		}

#ifdef LPSOLVE
		/* Now get rid of any rows that have become	*/
		/* slack.  (We don't lose these constraints:	*/
		/* they're still sitting around in the		*/
		/* constraint pool.)				*/
		delete_slack_rows_from_LP (bbip);
#endif

		/* Add new contraints to the constraint pool. */
		num_const = add_constraints (bbip, cp);

		if (num_const <= 0) {
			/* Separation routines found violations, but	*/
			/* the constraint pool disagrees...		*/
			fatal ("compute_good_lower_bound: Bug 4.");
		}

		while (cp NE NULL) {
			tmp = cp;
			cp = tmp -> next;
			free ((char *) (tmp -> mask));
			free ((char *) tmp);
		}
		++iteration;
	}

#if 1
	/* Print execution times of final iteration... */
	convert_cpu_time (0, time_str);
	--Tp;
	while (Tp > &Tn [0]) {
		--Tp;
		convert_cpu_time (Tp [1] - Tp [0], tbuf1);
		(void) sprintf (tbuf2, "%s/%s", tbuf1, time_str);
		strcpy (time_str, tbuf2);
	}
	(void) tracef ("  %% Final iteration: %s\n", time_str);
#endif

	/* Only get here with fractional solution and	*/
	/* no more violated constraints were found.	*/

	return (LB_FRACTIONAL);
}

/*
 * Routine to print out an updated lower bound value.
 */

	static
	void
new_lower_bound (

double			lb,		/* IN - new lower bound value */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
struct cinfo *		cip;
struct scale_info *	sip;
double			prev;
double			old_gap;
double			new_gap;
cpu_time_t		t1;
char			buf1 [32];

	cip = bbip -> cip;
	sip = &(cip -> scale);

	prev = bbip -> prevlb;
	if (lb <= prev) {
		/* There has been no improvement - get out. */
		return;
	}

	if (prev <= -DBL_MAX) {
		/* Don't make lower bound jump from initial value... */
		prev = lb;
	}

	/* Print out the old and new lower bounds, with timestamp. */
	t1 = get_cpu_time ();
	convert_cpu_time (t1 - bbip -> t0, buf1);

	if ((bbip -> best_z >= DBL_MAX) OR (bbip -> best_z EQ 0.0)) {
		old_gap = 99.9;
		new_gap = 99.9;
	}
	else {
		old_gap = 100.0 * (bbip -> best_z - prev) / bbip -> best_z;
		new_gap = 100.0 * (bbip -> best_z - lb) / bbip -> best_z;
	}

	tracef (" %% @LO %s %24.20f %2.10f\n",
		buf1, UNSCALE (prev, sip), old_gap);
	tracef (" %% @LN %s %24.20f %2.10f\n",
		buf1, UNSCALE (lb, sip), new_gap);

	bbip -> prevlb = lb;
}

/*
 * Routine to print out an updated upper bound value.
 */

	void
new_upper_bound (

double			ub,		/* IN - new upper bound value */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
struct cinfo *		cip;
struct scale_info *	sip;
double			prev;
int			i;
cpu_time_t		t1;
double			old_gap;
double			new_gap;
char			buf1 [64];
char			buf2 [64];

	cip = bbip -> cip;
	sip = &(cip -> scale);

	prev = bbip -> best_z;
	if (ub >= prev) {
		/* Supposed to be an improvement! */
		fatal ("new_upper_bound: Bug 1.");
	}

	/* We HAVE a new best solution! */
	bbip -> best_z = ub;

#if CPLEX
	{ double toobig, toosmall, ulim;
	  /* Set new cutoff value for future LPs... */
	  ulim = ldexp (ub, -(bbip -> lpmem -> obj_scale));
	  if (_MYCPX_setobjulim (ulim, &toosmall, &toobig) NE 0) {
		fatal ("new_upper_bound: Bug 2.");
	  }
	}
#endif

#if LPSOLVE
	/* Set new cutoff value for future LPs... */
	/* (This may not really work in lp_solve.) */
	bbip -> lp -> obj_bound = ub;
#endif

	cut_off_existing_nodes (ub, bbip -> bbtree);

	/* Might want to do this if all other nodes were cut off. */
	update_node_preempt_value (bbip);

	/* Now print out the trace messages. */
	if (prev >= DBL_MAX) {
		/* Don't make upper bound jump from infinity... */
		prev = ub;