decluster.w,v

Lin-Kernighan heuristic for the TSP and minimum weight perfect matching

{
int i;
for (i=0;i<n;i++) {
	digest[i].cost=0;
}
}

@@
@@<Other cleanup code@@>=
free_mem(digest);
mem_deduct((n+n-1)*sizeof(digest_t));

@@
@@<Copy |cost| fields from |T_prime|@@>=
{
int i;
for(i=n;i<n+n-1;i++) digest[i].cost=T_prime->edge[i].cost;
}

@@ Entry |digest[i].level| is just the depth of vertex |i|.   
We compute all of level numbers in linear time with a breadth first search.
The search begins at the root of the tree, |T_prime->edge[n+n-2]|.

Array |queue| serves as a last-in first-out queue.  It is indexed by 
a read cursor |r| and a write cursor |w|.

@@<Compute |level| numbers@@>= 
{ int *queue = new_arr_of(int,n+n-1), r, w, i, ch;
digest[n+n-2].level=0;
queue[0]=n+n-2;
for ( r=0, w=1; r < w ; r++ ) {
	const int here = queue[r], cur_level = digest[here].level+1;
	for(i=0;i<2;i++) if ( (ch=T_prime->edge[here].child[i]) != NO_CHILD ) {
		digest[ch].level=cur_level;
		queue[w++]=ch;
	}
}
free_mem(queue);mem_deduct((n+n-1)*sizeof(int));
}


@@ The inlabel number of a vertex is a function of its preorder number
and the size of subtree rooted there.    These are stored in arrays
|preorder| and |size|, naturally.  
Value |preorder[i]| is one more than the number of vertices visited before
vertex |i| in a preorder traversal of the tree.
Value |size[i]| is the number of vertices in the subtree rooted at |i|.


@@<Compute |inlabel| numbers@@>=
{
int *preorder=new_arr_of(int,n+n-1), *size=new_arr_of(int,n+n-1);
int preorder_number=0;
#define DFS_INLABEL
@@<Do a depth-first search of |T_prime|@@>@@;
#undef DFS_INLABEL

free_mem(preorder);free_mem(size);
mem_deduct(2*(n+n-1)*sizeof(int));
}


@@
Both |preorder| and |size| statistics 
 are computed via a depth-first search of the tree.
I'll separate the bookkeeping from the real action by using subsections.
This way I can use the same \CWEB/ bookkeeping code twice.

Array |in| stores the roots of the subtrees we're currently exploring.
It behaves as a stack.  Entry |cur_child[i]| tells
us  which 
child (0 or 1) of vertex |in[i]| is currently 
being explored.

We know the graph is acyclic, so we don't need to mark nodes.
Also, the graph is a binary tree.

We begin at the root, vertex |n+n-2|.

@@<Do a depth-first search of |T_prime|@@>=
{
int *in=new_arr_of(int,n+n-1), *cur_child=new_arr_of(int,n+n-1);
int top, here, next_child;

top=0;
in[top]=n+n-2;
cur_child[top]=-1;

while(top>=0) {
	here=in[top];
	switch(cur_child[top]) {
	case -1: /* Visit myself, then fall through. */
		@@<Visit |here| the first time@@>@@;
	case 0: /* Visit the children. */
		next_child=T_prime->edge[here].child[++cur_child[top]];
		if ( next_child != NO_CHILD ) {
			/* Push that child on the stack. */
			in[++top] = next_child;
			cur_child[top] = -1;
		}
		break;
	default:
		@@<Visit |here| the last time@@>@@;
		/* $\ldots$ and pop the stack. */
		top--;
	}
}
free_mem(in);free_mem(cur_child);
mem_deduct(2*(n+n-1)*sizeof(int));
}

@@ For |inlabel| numbers, we need to update the preorder number as we
enter a node.

@@<Visit |here| the first time@@>=
#if defined(DFS_INLABEL)
preorder[here] = ++preorder_number; 
#endif 

@@ For |inlabel| numbers, we need to update the |size| array as we leave
a node.

I'll define |inlabel[here]| but not explain its use.  
The preorder numbers of all the nodes in the subtree of |T_prime| rooted
at |here| range from |preorder[here]| through
|preorder[here]+size[here]-1|.  In fact, together they form precisely
that interval, no more and no less.
Value |inlabel[here]| is the value in that interval having most trailing
0 bits when expressed in binary.

Collectively the |inlabel| fields form  an injective mapping from
the |N := n+n-1| vertices of |T_prime| into a complete binary tree $B$ 
of depth $\lfloor \log_2 N\rfloor +1$.

See below for how |inlabel|s are used.

With care I can elide the taking of logarithms, as I've done below 
in the \LCA\ computation.  Maybe later.  Or maybe never, since
this is just an initialization phase, so timing is not so critical.

@@<Visit |here| the last time@@>=
#if defined(DFS_INLABEL)
{ const int child0=T_prime->edge[here].child[0],
			child1=T_prime->edge[here].child[1];
size[here] = 1 + (child0==NO_CHILD?0:size[child0])
			   + (child1==NO_CHILD?0:size[child1]);
}
{
const unsigned int last=preorder[here]+size[here]-1, 
	i=floor_log_2((preorder[here]-1)^last);
	digest[here].inlabel=hi_mask(i)&last; /* All but the bottom $i$ bits of |last|. */
}
#endif

@@ The |ascendant| numbers may also be computed during a depth-first search.
Again, see Schieber.  The |inlabel| numbers must be computed first.

@@<Compute |ascendant| numbers@@>=
#define DFS_ASCENDANT
@@<Do a depth-first search of |T_prime|@@>@@;
#undef DFS_ASCENDANT

@@ Assuming |x| is positive, value |(x^(x-1))&x| is the rightmost 1 bit
of |x|, but isolated.


@@<Visit |here| the first time@@>=
#if defined(DFS_ASCENDANT)
if (top==0) { /* The root. */
	digest[here].ascendant=digest[here].inlabel;
} else {
	const int p=in[top-1], ap=digest[p].ascendant, ip=digest[p].inlabel; /* Parent data. */
	if ( digest[here].inlabel==ip ) digest[here].ascendant=ap;
	else {
		const int ih=digest[here].inlabel;
		digest[here].ascendant=ap+((ih^(ih-1))&ih);
	}
}
#endif

@@
Schieber also prescribes that we compute a |head| table 
telling us the head of each inlabel chain. That is |head[i]| is
the shallowest vertex in |T_prime| with inlabel |i|.

We only ever want to know the |parent| of a head, so instead
I compute  a 
|parent_of_head| with the following defining equation:
$$\hbox{|parent_of_head[inlabel_w]==parent(head[inlabel_w])|}.$$
This saves an indirection.

All the |parent_of_head|
fields can be defined in four easy passes
over the array.  

But we do need the |head| array during initialization.
It is indexed by inlabel number, so it ranges from 1 through |n+n|.
We also compute a |parent| array, indexed by vertex (in |T_prime|).

@@<Compute |parent_of_head| numbers@@>=
{
int i, j, *head=new_arr_of(int,n+n), *parent=new_arr_of(int,n+n-1);
for (i=0;i<n+n;i++) head[i]=-1; /* Empty each bin. */
for (i=0;i<n+n-1;i++) { /* Compute heads of inlabel chains */
	const int ii=digest[i].inlabel, hii=head[ii];
	if ( hii==-1 || digest[i].level < digest[hii].level ) head[ii]=i;
}
for (i=0;i<n+n-1;i++) parent[i]=NO_PARENT,parent_of_head[i]=NO_PARENT; 
	/* Make parents null. */
for (i=0;i<n+n-1;i++)  /* Compute parents. */
	for ( j=0;j<2;j++) {
		const int child=T_prime->edge[i].child[j];
		if (child!=NO_CHILD) parent[child]=i;
	}
for (i=0;i<n+n-1;i++) /* Record parents of heads of chains. */
	parent_of_head[digest[i].inlabel] = parent[head[digest[i].inlabel]];
@@<Verbose: print |parent|@@>@@;
@@<Verbose: print |head|@@>@@;

free_mem(head);mem_deduct((n+n)*sizeof(int));
free_mem(parent);mem_deduct((n+n-1)*sizeof(int));
}

@@ We need to define |parent_of_head|.
@@<Module variables@@>=
int *parent_of_head;

@@ It must be allocated$\ldots$
@@<Other setup code@@>=
parent_of_head=new_arr_of(int,n+n);

@@ $\ldots$ and freed.
@@<Other cleanup code@@>=
free_mem(parent_of_head);mem_deduct((n+n)*sizeof(int));


@@*Least common ancestor queries.
We're now ready to process least common ancestor queries.

@@<Module subroutines@@>=
static inline int
decluster_lca(int x, int y)
{
	const digest_t xd=digest[x], yd=digest[y];
	if ( xd.inlabel == yd.inlabel ) return (xd.level<= yd.level)?x:y;
	else { /* The \LCA\ is neither |x| nor |y|. */
		const int xil = xd.inlabel, yil=yd.inlabel;
		int b, inlabel_z, common, jmask, xpp, ypp;
		@@<Set |b| to the \LCA\ of $C_x$ and $C_y$ in the complete tree $B$@@>@@;
		@@<Find |inlabel_z|, |common|, and |jmask|@@>@@;
		@@<Find |xpp| and |ypp|@@>@@;
		return (digest[xpp].level <= digest[ypp].level) ? xpp : ypp;
	}
}


@@ 
The digestion of |T_prime| decomposes it into (contiguous) chains of vertices.
Chain $C_u$ is the chain containing vertex $u$ of |T_prime|.
The inlabel numbers form an injective mapping from these chains of
|T_prime| into  a complete binary tree |B|.  
Tree $B$ is not explicitly stored, it's existence is implicit in the inlabel
numbers.  Futhermore, |inlabel[u]| is the ordinal position of the
vertex $C_u$ of $B$ in an in-order traversal of $B$, counting from 1.
@@^inlabel numbers@@>
@@^in-order traversal@@>
@@^tree $B$@@>

Finding least common ancestors in binary trees is easy, given in-order
numbers.  

For clarity,  Schieber describes
things in terms of logarithms, but in the interest of speed
I implemented the computation without logarithms.

The expression of a 1 followed by many zeroes, namely,
|t^(t>>1)|, used to be |(t+1)>>1|.  The old code is unsafe when |t| 
is the largest positive |int|, \eg, $2^{31}-1$, because adding 1 to
this value makes a large negative number.  Then the right shift might
extend the sign down and preserve the negative sign bit.  

The new
\def\AC{{\it AC}}%
\def\NC{{\it NC}}%
code may even be faster because it replaces an $\AC^0$ operation (addition)
and an $\NC^0$ operation (a single bit shift)
with
$NC^0$ operations (exclusive-or, logical and, and
a single bit shift).  On today's machines, who knows?


%$(t+1)>>1$ is unsafe when $t=2^{31}-1$ because of possible sign extension
%on the right shift.
%Replace it with $t\XOR(t>>1)$. (Here t is a low bits mask.)  Also, this
%might be faster because it replaces an $AC^0$ operation (t+1) with an $NC^0$
%operation (exclusive-or).

% at d rightmost_1(x) ((x^(x-1))&x)

@@d min(A,B) ((A)<(B)?(A):(B))
@@d max(A,B) ((A)>(B)?(A):(B))

@@<Set |b| to the \LCA\ of $C_x$ and $C_y$ in the complete tree $B$@@>=
{
const unsigned int xfuzz = xil^(xil-1), yfuzz = yil^(yil-1);
const unsigned int lomask=max(xfuzz,yfuzz)>>1;
int hi_diff = (xil^yil)&(~lomask);
if ( hi_diff ) { /* Neither ancestor nor descendant. */
	int t=hi_diff;  
	copy_1_down(t); 
		/* Sort of ``sign extend'' the top 1 bit of |t| down through the
lower bits. */
	b=(~t)&xil; /* Keep only the top few bits common to both |xil| and |yil|. */
	b|=t^(t>>1); /* But append a 1 followed by many zeros. */
} else { /* One is an ancestor of the other. */
	b=(xfuzz>=yfuzz)?xil:yil;
}
@@<Verbose: print |b|@@>@@;
}

@@ Let $z=\LCA(x,y)$. This section computes |inlabel_z|$:=$|inlabel[z]| 
without directly
computing |z| itself.  We use the ascendant numbers and some bit masking.
Again Schieber describes things in terms of logarithms, but I have elided
the logarithms.

In declarative form, 
|inlabel[z]|$:=$|inlabel[x]&~lo_mask(j)|
where |j=index_rightmost_1(u)|,
|u=common&~lo_mask(i)|, |common=ascendant[x]&ascendant[y]|,
and 
|i=index_rightmost_1(b)|.

Value |jmask| is |j| low-order 1 bits.

@@<Find |inlabel_z|, |common|, and |jmask|@@>=
{ const int imask=(b^(b-1))>>1; /* Mask of bits corresponding to low 0's in $b$. */
const int u = (common=xd.ascendant & yd.ascendant) & ~imask;
jmask = (u^(u-1))>>1; /* Mask of bits corresponding to low 0's in
$u$. */
inlabel_z = (xd.inlabel & ~jmask)|(jmask+1); /* |yd.inlabel| would also do. */
@@<Verbose: print |inlabel_z|, |common|, and |jmask|@@>@@;
}




@@ Vertex |xpp| is either |x| itself or is the parent of the head
of the inlabel chain at |inlabel_w_x|. Value |inlabel_w_x| is
all but the bottom |k+1| bits of |inlabel[x]|, followed by a 1 and
|k| zeroes.

Vertex |ypp| is similar.

@@<Find |xpp| and |ypp|@@>=
if (xd.inlabel==inlabel_z) xpp=x;
else {
	int inlabel_w_x;
	int kmask=xd.ascendant & jmask;
	copy_1_down(kmask);
	inlabel_w_x =(~kmask)&xil; /* Keep only the top few bits of |inlabel[x]|. */
	inlabel_w_x |=kmask^(kmask>>1); /* But append a 1 followed by many zeros. */
	xpp=parent_of_head[inlabel_w_x];
	@@<Verbose: print |inlabel_w_x|, |xpp|@@>@@;
}
if (yd.inlabel==inlabel_z) ypp=y;
else {
	int inlabel_w_y;
	int kmask=yd.ascendant & jmask;
	copy_1_down(kmask);
	inlabel_w_y =(~kmask)&yil; /* Keep only the top few bits of |inlabel[y]|. */
	inlabel_w_y |=kmask^(kmask>>1); /* But append a 1 followed by many zeros. */
	ypp=parent_of_head[inlabel_w_y];
	@@<Verbose: print |inlabel_w_y|, |ypp|@@>@@;
}


@@ Given the least common ancestor function, computing the cluster distance
|d| is very easy.

We might consider reducing memory requirements by computing the
cost on the fly instead of storing it in the topology tree |T_prime|.
@@<Subroutines@@>=
length_t 
decluster_d(int u, int v)
{
	return digest[decluster_lca(u,v)].cost;
}

@@*Extras.
Here are the extra support routines. That is, these routines are
not required for supporting decluster distance queries.

@@ We may need to see what |T_prime| looks like.
The caller can safely change its contents if want to, because we've
already preprocessed it.

@@<Subroutines@@>=
decluster_tree_t *
decluster_topology_tree(void)
{
	return T_prime;
}

@@ It is be useful to be able to output a tree.
@@<Subroutines@@>=
void
decluster_print_tree(FILE *out,decluster_tree_t const *t, const char *name) 
{
	if ( t )  {
		int n=t->n, i;
		const char *print_name = name ? name : "";
		errorif(t==NULL,"decluster_print_tree: given a NULL tree\n");
		errorif(n<0,"decluster_print_tree: tree %s size %d < 0\n",print_name, n);
		fprintf(out,"%s->n==%d\n",print_name,t->n);
		for (i=0;i<n;i++) {
			fprintf(out," %d (%d,%d) "length_t_spec"\n",
				i,
				t->edge[i].city[0], t->edge[i].city[1],
				length_t_pcast(t->edge[i].cost));
		}
	} else {
		fprintf(out,"Tree %s is NULL\n",name);
		fprintf(out,"For more data, make sure variable decluster_discard_topology_tree is zero)\n");
	}
}

@@*Testing.
Since I don't claim to understand the least common ancestor code
completely, let's test it.

@@(declustertest.c@@>=
#include <config.h>
#include "lkconfig.h"
#include <stdio.h>
#include <stdlib.h>
#include "error.h