bb.c

生成直角Steiner树的程序包

/***********************************************************************

	File:	bb.c
	Rev:	b-2
	Date:	02/28/2001

	Copyright (c) 1995, 2001 by David M. Warme

************************************************************************

	The guts of the branch-and-cut.

************************************************************************

	Modification Log:

	a-1:	11/17/95	warme
		: Created.
	b-1:	11/14/96	warme
		: Renamed this program.
		: Split off the cutset stuff into another file.
		: Other cleanups for release.
	b-2:	02/28/2001	warme
		: Numerous changes for 3.1 release.
		: Split the main() routine off into bbmain.c.
		: Split certain utility routines off into bbsubs.c.
		: Major improvements to branch var selection.
		: New create_bbinfo routine does setup for
		:  branch_and_cut.
		: Enable variable fixing code.

************************************************************************/

#include "bb.h"
#include "bbsubs.h"
#include "config.h"
#include "constrnt.h"
#include "cra.h"
#include "cutset.h"
#include "genps.h"
#include "p1io.h"
#include "sec2.h"
#include "sec_comp.h"
#include "sec_heur.h"
#include "steiner.h"
#include "ub.h"

#include <signal.h>


/*
 * Global Routines
 */

dist_t			branch_and_cut (struct bbinfo *);
bool			check_for_better_IFS (double *,
					      struct bbinfo *,
					      double *);
struct bbinfo *		create_bbinfo (struct cinfo *);
void			new_upper_bound (double, struct bbinfo *);


int		Branch_Var_Policy = 1;
int		Check_Branch_Vars_Thoroughly = 1;
bool		Check_Root_Constraints = FALSE;
bool		Choose_Branch_Vars_Carefully = TRUE;
double		Initial_Upper_Bound = DBL_MAX;
bool		Print_Root_LP = FALSE;

volatile bool	force_branch_flag;


/*
 * Local Equates
 */


	/* Lower bound outcomes... */
#define	LB_INFEASIBLE		1
#define	LB_CUTOFF		2
#define	LB_INTEGRAL		3
#define	LB_FRACTIONAL		4
#define LB_PREEMPTED		5

	/* Variable fixing outcomes... */
#define	VFIX_NOTHING_FIXED	0
#define	VFIX_VARIABLES_FIXED	1
#define	VFIX_FIXED_FRACTIONAL	2
#define	VFIX_INFEASIBLE		3

#define	UP_FIRST	TRUE


/*
 * Local Types
 */

struct bvar {			/* Info about best branch var seen */
	int	var;		/* Best branch var, or -1 */
	double	z0;		/* Objective when Xj=0 */
	double	z1;		/* Objective when Xj=1 */
	double	test_2nd_val;	/* Only check 2nd branch if 1st > this. */
};

#ifdef CPLEX

struct basis_save {	/* Structure to save basis state for CPLEX */
	int *		cstat;
	int *		rstat;
};

#endif


/*
 * Local Routines
 */

static int		carefully_choose_branching_variable (struct bbinfo *,
							     double *,
							     double *);
static void		change_var_bounds (LP_t *, int, double, double);
static void		check_root_constraints (struct bbinfo *);
static int		choose_branching_variable (struct bbinfo *,
						   double *,
						   double *);
static bool		compare_branch_vars (struct bbinfo *,
					     int,
					     struct bvar *);
static int		compute_good_lower_bound (struct bbinfo *);
static void		cut_off_existing_nodes (double, struct bbtree *);
static struct constraint * do_separations (struct bbinfo *,
					   cpu_time_t **);
static bool		eval_branch_var (struct bbinfo *,
					 int,
					 int,
					 struct basis_save *,
					 double);
static int		fix_variables (struct bbinfo *,
				       int *, int,
				       int *, int);
static RETSIGTYPE	force_branch_signal_handler (int);
static bool		integer_feasible_solution (double *,
						   bitmap_t *,
						   bitmap_t *,
						   struct cinfo *,
						   int *);
static void		new_lower_bound (double, struct bbinfo *);
static void		print_root_lp (struct bbinfo *);
static int		reduced_cost_var_fixing (struct bbinfo *);
static struct bbnode *	select_next_node (struct bbtree *);
static void		sort_branching_vars (int *, int, double *);
static void		trace_node (struct bbinfo *, char, char *);
static void		update_node_preempt_value (struct bbinfo *);

#ifdef CPLEX
static void		destroy_LP_basis (struct basis_save *);
static double		try_branch (LP_t *,
				    int,
				    int,
				    double *,
				    double,
				    struct basis_save *);
static void		save_LP_basis (LP_t *, struct basis_save *);
#endif


/*
 * Local Variables
 */

static FILE *		rcfile;

/*
 * Set up the initial branch-and-bound problem, including the root
 * node and the initial constraint pool.
 */

	struct bbinfo *
create_bbinfo (

struct cinfo *		cip		/* IN - compatibility info */
)
{
struct bbinfo *		bbip;
int			i;
int			j;
int			k;
int			n;
int			nmasks;
int			nedges;
LP_t *			lp;
int			status;
bitmap_t *		vert_mask;
bitmap_t *		edge_mask;
bitmap_t *		req_edges;
bitmap_t *		fixed;
bitmap_t *		value;
bitmap_t *		smt;
struct cpool *		cpool;
struct bbstats *	statp;
struct bbtree *		bbtree;
struct bbnode *		root;
struct lpmem *		lpmem;
struct rcon *		rcp;

	/* Create the global branch-and-bound info structure... */
	bbip = NEW (struct bbinfo);
	memset (bbip, 0, sizeof (*bbip));

	bbip -> t0 = get_cpu_time ();

	nmasks = cip -> num_edge_masks;
	nedges = cip -> num_edges;

	vert_mask	= cip -> initial_vert_mask;
	edge_mask	= cip -> initial_edge_mask;
	req_edges	= cip -> required_edges;

	/* initialize global pool of constraints. */
	cpool = NEW (struct cpool);
	initialize_constraint_pool (cpool, vert_mask, edge_mask, cip);

	/* Build initial formulation. */
	lpmem = NEW (struct lpmem);
	lp = build_initial_formulation (cpool,
					vert_mask,
					edge_mask,
					cip,
					lpmem);

	/* Initialize the branch-and-bound tree... */
	bbtree = create_bbtree (nmasks);

	/* Create vectors to describe the current problem... */
	fixed	= NEWA (nmasks, bitmap_t);
	value	= NEWA (nmasks, bitmap_t);
	smt	= NEWA (nmasks, bitmap_t);

	for (i = 0; i < nmasks; i++) {
		fixed [i] = 0;
		value [i] = 0;
		smt [i]	  = 0;
	}
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) {
			/* variables that are outside of the problem */
			/* are fixed at zero... */
			SETBIT (fixed, i);
			change_var_bounds (lp, i, 0.0, 0.0);
		}
		else if (BITON (req_edges, i)) {
			/* Front-end has determined that this hyperedge	*/
			/* MUST be present in an optimal solution!	*/
			SETBIT (fixed, i);
			SETBIT (value, i);
			change_var_bounds (lp, i, 1.0, 1.0);
		}
	}

	/* Create the root node... */
	root = NEW (struct bbnode);
	memset (root, 0, sizeof (*root));

	root -> z	= -DBL_MAX;
	root -> optimal	= FALSE;
	root -> num	= (bbtree -> snum)++;
	root -> iter	= 0;
	root -> parent	= -1;
	for (i = 0; i < NUM_BB_HEAPS; i++) {
		root -> index [i] = -1;
	}
	root -> var	= -1;
	root -> dir	= 0;
	root -> depth	= 0;
	root -> br1cnt	= 0;
	root -> x	= NEWA (nedges, double);
	root -> cpiter	= -1;		/* x is not current. */
	root -> zlb	= NEWA (2 * nedges, double);
	root -> fixed	= fixed;
	root -> value	= value;
	root -> n_uids	= 0;
	root -> bc_uids	= NULL;
	root -> bc_row	= NULL;
	root -> rstat	= NULL;
	root -> cstat	= NULL;
	root -> bheur	= NEWA (nedges, double);
	root -> next	= NULL;
	root -> prev	= NULL;

	for (i = 0; i < nedges; i++) {
		root -> bheur [i] = 0.0;
	}

	for (i = 0; i < 2*nedges; i++) {
		root -> zlb [i] = -DBL_MAX;
	}

	/* Create the branch-and-bound statistics structure... */
	statp = NEW (struct bbstats);
	memset (statp, 0, sizeof (*statp));

	statp -> n		= cip -> num_verts;
	statp -> m		= nedges;
	statp -> p1time		= cip -> p1time;
	statp -> num_nodes	= 0;
	statp -> num_lps	= 0;

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

	/* Fill in the global branch-and-bound info structure... */
	bbip -> cip		= cip;
	bbip -> tmap		= vert_mask;
	bbip -> fset_mask	= edge_mask;
	bbip -> lp		= lp;
	bbip -> lpmem		= lpmem;
	bbip -> cpool		= cpool;
	bbip -> bbtree		= bbtree;
	bbip -> csip		= NULL;
	bbip -> preempt_z	= Initial_Upper_Bound;
	bbip -> best_z		= Initial_Upper_Bound;
	bbip -> smt		= smt;
	bbip -> node		= root;
	bbip -> slack_size	= 0;
	bbip -> slack		= NULL;
	bbip -> dj		= NEWA (nedges, double);
	bbip -> fixed		= NULL;
	bbip -> value		= NULL;
	bbip -> statp		= statp;
	bbip -> prevlb		= -DBL_MAX;
	bbip -> ubip		= NULL;

	/* Make the root node inactive by putting it in the bbtree... */

	n = cpool -> nlprows;
	root -> n_uids	= n;
	root -> bc_uids	= NEWA (n, int);
	root -> bc_row	= NEWA (n, int);

	j = 0;
	rcp = &(cpool -> rows [0]);
	for (i = 0; i < cpool -> nrows; i++, rcp++) {
		k = rcp -> lprow;
		if (k < 0) continue;
		++(rcp -> refc);
		root -> bc_uids [j] = rcp -> uid;
		root -> bc_row [j]  = k;
		++j;
	}
	if (j NE n) {
		fatal ("create_bbinfo: Bug 1.");
	}

	root -> next = NULL;
	root -> prev = NULL;
	bbtree -> first = root;

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

	return (bbip);
}

/*
 * This routine is the top-level of the branch-and-cut.
 */

	dist_t
branch_and_cut (

struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			j;
int			nmasks;
int			nedges;
int			status;
bitmap_t *		fixed;
bitmap_t *		value;
bitmap_t *		delta;
bitmap_t *		smt;
double *		x;
double			z0;
double			z1;
double			tmpz;
cpu_time_t		t1;
RETSIGTYPE		(*save_sigterm) (int);
struct cinfo *		cip;
LP_t *			lp;
struct cpool *		cpool;
struct bbtree *		bbtree;
struct bbstats *	statp;
struct bbnode *		node;
struct bbnode *		node_to_free;
struct bbnode *		node2;

#ifdef CPLEX
int *			b_index;
char *			b_lu;
double *		b_bd;
double			objlim;
double			save_objlim;
#endif

	cip	= bbip -> cip;
	cpool	= bbip -> cpool;
	lp	= bbip -> lp;
	statp	= bbip -> statp;
	smt	= bbip -> smt;
	bbtree	= bbip -> bbtree;

	nmasks = cip -> num_edge_masks;
	nedges = cip -> num_edges;

#if CPLEX
	/* Save the existing objective limit, and set it to infinity. */
	CPXgetdblparam (cplex_env, CPX_PARAM_OBJULIM, &save_objlim);
	objlim = DBL_MAX;
	CPXsetdblparam (cplex_env, CPX_PARAM_OBJULIM, objlim);
#endif

	if (Check_Root_Constraints) {
		rcfile = fopen ("/tmp/lp.x", "w");
		if (rcfile EQ NULL) {
			fprintf (stderr,
				 "Warning: unable to create /tmp/lp.x\n");
			Check_Root_Constraints = FALSE;
		}
	}

	/* Establish handler for the SIGTERM (kill -15) signal... */
	force_branch_flag = FALSE;
	save_sigterm = signal (SIGTERM, force_branch_signal_handler);

	/* Create vectors to describe the current problem... */
	fixed	= NEWA (nmasks, bitmap_t);
	value	= NEWA (nmasks, bitmap_t);
	delta	= NEWA (nmasks, bitmap_t);

	/* No variables have been fixed yet... */
	for (i = 0; i < nmasks; i++) {
		fixed [i] = 0;
		value [i] = 0;
	}
	bbip -> fixed = fixed;
	bbip -> value = value;

#ifdef CPLEX
	/* Create arrays for changing variable bounds... */
	b_index	= NEWA (2 * nedges, int);
	b_lu	= NEWA (2 * nedges, char);
	b_bd	= NEWA (2 * nedges, double);
#endif

#if 0
	/* Build cutset separation formulation. */
	bbip -> csip = NEW (struct cs_info);
	build_cutset_separation_formulation (vert_mask,
					     edge_mask,
					     cip,
					     bbip -> csip);
#endif

	/* Init the heuristic upper bound. */

	bbip -> ubip = startup_heuristic_upper_bound (bbip -> cip);


	/* At this point, all nodes are inactive. */

	for (;;) {

		/* Select the next node to process. */
		node = select_next_node (bbtree);
		if (node EQ NULL) break;

		/* This is perhaps a new lower bound... */
		new_lower_bound (node -> z, bbip);

		if (node -> z > -DBL_MAX) {
			tracef ("%% Resuming node %d at %24.20f\n",
				node -> num,
				UNSCALE (node -> z, &(cip -> scale)));
		}
		else {
			tracef ("%% Resuming node %d\n", node -> num);
		}

		/* Restore the LP tableaux and basis for this node.	*/
		/* Decrement the reference counts on this node's	*/
		/* binding rows.  Since it is now the ACTIVE node,	*/
		/* there is no reason to continue protecting these rows	*/
		/* from deletion by other nodes.			*/
		restore_node_basis (node, bbip);

		/* Determine new preemption value (i.e. the objective	*/
		/* value of the next-best node).			*/
		update_node_preempt_value (bbip);

		/* Modify LP to represent problem from new node. */
		for (i = 0; i < nmasks; i++) {
			delta [i] =   (fixed [i] ^ node -> fixed [i])
				    | (value [i] ^ node -> value [i]);
		}
#ifdef CPLEX
		j = 0;
		for (i = 0; i < nedges; i++) {
			if (NOT BITON (delta, i)) continue;
			/* Force bounds for variable 'i' to be correct... */
			b_index [j]	= i;	/* variable i, */
			b_lu [j]	= 'L';	/*	lower bound */
			b_index [j+1]	= i;	/* variable i, */
			b_lu [j+1]	= 'U';	/*	upper bound */
			if (NOT BITON (node -> fixed, i)) {
				/* new variable is NOT fixed... */
				b_bd [j]	= 0.0;
				b_bd [j+1]	= 1.0;
			}
			else if (NOT BITON (node -> value, i)) {
				/* new variable is fixed to 0 */
				b_bd [j]	= 0.0;
				b_bd [j+1]	= 0.0;
			}
			else {
				/* new variable is fixed to 1 */
				b_bd [j]	= 1.0;
				b_bd [j+1]	= 1.0;
			}
			j += 2;
		}
		if (j > 0) {
			if (_MYCPX_chgbds (bbip -> lp, j, b_index, b_lu, b_bd) NE 0) {
				fatal ("branch_and_cut: Bug 1.");
			}
#if 0
			++(bbip -> cpool -> uid);
#endif
		}
#endif

#ifdef LPSOLVE
		j = 0;
		for (i = 0; i < nedges; i++) {
			if (NOT BITON (delta, i)) continue;
			++j;
			/* Force bounds on variable 'i' to be correct... */
			if (NOT BITON (node -> fixed, i)) {
				/* variable is NOT fixed... */
				set_bounds (lp, i + 1, 0.0, 1.0);
			}
			else if (NOT BITON (node -> value, i)) {
				/* variable is fixed to 0 */
				set_bounds (lp, i + 1, 0.0, 0.0);
			}
			else {
				/* variable is fixed to 1 -- must set	*/
				/* bounds in this order to avoid	*/
				/* lb > ub condition between calls...	*/
				set_bounds (lp, i + 1, 1.0, 1.0);
			}
		}
		if (j > 0) {
#if 0
			++(bbip -> cpool -> uid);
#endif
		}
#endif

		for (i = 0; i < nmasks; i++) {
			fixed [i] = node -> fixed [i];
			value [i] = node -> value [i];
		}

		/* Set up new node to be processed */
		bbip -> node = node;
		if (node -> iter <= 0) {
			/* Haven't processed this node before	*/
			/* -- tally another node...		*/
			++(statp -> num_nodes);
		}

		/* Process the current node... */
		status = compute_good_lower_bound (bbip);

		if (node -> depth EQ 0) {