twolevel.w,v

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

balancing in case 1.  In particular, we've got code to split off the 
portion of a segment to the left of the pointer |ca|, and code
to split off code to the right of pointer |cc|.

That splitting code takes time proportional to the length of list that is
split off.  So we prefer to split off the shorter portion of the segment.

@@<Split the $a-b$ segment@@>=
@@<Verbose: split the $a-b$ segment@@>@@;
{ city_node_t *l, *r; parent_node_t *p = ca->parent; /* Same as |cb->parent| */
if ( ca->seq < cb->seq ) l=ca, r=cb;
else l=cb, r=ca;
if ( l->seq - p->head->seq < p->tail->seq - r->seq ) {
	city_node_t *ca=r;	/* No relation to the previous |ca| */
	@@<Split off the end to the left of |ca|@@>@@;
} else {
	city_node_t *cc=l;	/* No relation to the previous |cc| */
	@@<Split off the end to the right of |cc|@@>@@;
}
}
@@<Renew the parent sequence numbers@@>@@;

@@ We've moved one of $a$ or $b$.  We need to update its parent sequence number.
It's simpler (and probably faster) to just update both than to figure
out (or remember) which one was moved.

Cities |c| and |d| might have been moved with the $a-b$ split.  We renew
their parent sequence numbers just in case.

The nice symmetry of this code allows us to use it after we split
the segment containing $c$ and $d$.

@@<Renew the parent sequence numbers@@>=
psa=ca->parent->seq; 
psb=cb->parent->seq;
psc=cc->parent->seq; 
psd=cd->parent->seq;

@@ Our aim so far in this third case is to arrange things so that case 2 applies,
\ie, so that both $a-c$ and $b-d$ are each a contiguous sequence of complete 
segments.  However once we've split the $a-b$ segment, the {\it first\/}
case might apply. That is, either $a-c$ or $b-d$ might lie entirely within
one segment.  
If this is the case, then we update the data structure as in
case 1.  Fortunately, \CWEB/ lets us do this quite naturally by reusing
that section.  Recall that if case 1 does apply, then that code
performs its own |return| to complete this procedure execution.

@@<Split the $a-b$ segment@@>=
@@<Handle case 1 of flipping, if it applies@@>@@;

@@ Now lets split the $c-d$ segment.  This is analogous to splitting the
$a-b$ segment.

@@<Split the $c-d$ segment@@>=
@@<Verbose: split the $c-d$ segment@@>@@;
{ city_node_t *l, *r; parent_node_t *p = cc->parent;  /* Same as |cd->parent| */
if ( cc->seq < cd->seq ) l=cc, r=cd;
else l=cd, r=cc;
if ( l->seq - p->head->seq < p->tail->seq - r->seq ) {
	city_node_t *ca=r;	/* No relation to the previous |ca| */
	@@<Split off the end to the left of |ca|@@>@@;
} else {
	city_node_t *cc=l;	/* No relation to the previous |cc| */
	@@<Split off the end to the right of |cc|@@>@@;
}
}
@@<Renew the parent sequence numbers@@>@@;

@@ Just as in the $a-b$ case, $c$ may have been moved into $a$'s segment,
or $d$ may have been moved into $b$'s segment.   If either has occurred,
then we're back to case 1.  So again, we check for case 1 and handle it
if it occurs.

@@<Split the $c-d$ segment@@>=
@@<Handle case 1 of flipping, if it applies@@>@@;

@@ Now we're ready to handle case 2: each of the $a-c$ and the $b-d$ portions
of the tour consists of a contiguous sequence of complete segments.

As with flipping a sequence of city nodes, we first rename cities so that
the portion to be flipped is delimited by the parents of cities |a| and |c|.

@@<Flip a sequence of segments@@>=
@@<Verbose: flip a sequence of segments@@>@@;
errorif((psa==psb||psc==psd),"psa %d==psb %d or psc %d == psd %d",psa,psb,psc,psd);
@@<Rename so that |psa| through |psc| is shorter than |psb| through |psd|@@>@@;
@@<Flip the sequence of segments from $a$ to $c$@@>@@;

@@
We prefer to flip the shorter of the two sequences of segments.  
However, we must know how the sequences are laid out before we
know which is shorter.  Schematically, the parents of the cities $a$, $b$,
$c$, and $d$ may be laid out in forward order, $abdc$, in reverse, $cdba$,
or in any of the three rotational variants of each of these:
$bdca$, $dcab$, $cabd$, or $dbac$, $bacd$, $acdb$.  We'd rather not 
expand the code unnecessarily by writing eight different versions, so
we'll rename cities so that the sequence of parents from |a| through |c|
is to be reversed.

The first step of this renaming task is to identify which  of the above
eight orderings hold.  There are two tricks to doing this easily.
First, we filter out trivial flips, \ie~where
one of the portions consists of a single parent.  
We want the predecessor and successor
of the parent of an end city (one of $a$, $b$, $c$, or $d$) to be 
in the other portion.  
Second, we must know the number of groups, |numgroups|; this makes identifying
successors and predecessors easy by using modular arithmetic.

We will make use of |normp| and |normm|, a pair of 
a fast but restricted  macros for the remainder modulo |num_groups|.
Macro |normp| normalizes numbers from $[0,2\hbox{|num_groups|}-1]$ down
to $[0,\hbox{|num_groups|}-1]$; it is used when adding two numbers (remember
$p$ for `plus').
Macro |normm| normalizes numbers from 
$[-\hbox{|num_groups|},\hbox{|num_groups|}-1]$ down
to $[0,\hbox{|num_groups|}-1]$; it is used when subtracting 
two numbers (remember
$m$ for `minus').

@@d normp(X) ((X)<num_groups ? (X) : (X)-num_groups)	
@@d normm(X) ((X)<0 ? (X)+num_groups : (X))
@@<Rename so that |psa| through |psc| is shorter than |psb| through |psd|@@>=
if ( psb==psd ) {
	SWAP(a,b,ti), @@+ SWAP(c,d,ti),@@+
	SWAP(ca,cb,tcn), @@+ SWAP(cc,cd,tcn),@@+
	SWAP(psa,psb,ti), @@+ SWAP(psc,psd,ti);
} else if ( psa!=psc ) {	/* Nontrivial sequence of parents. Determine the shorter of the two. */
	if ( normp(psa+1)==psb ) { /* Forward order */
		const int dmb = psd-psb, amc=psa-psc;
		if ( normm(dmb) < normm(amc) )
			SWAP(a,b,ti), @@+ SWAP(c,d,ti),@@+
			SWAP(ca,cb,tcn), @@+ SWAP(cc,cd,tcn),@@+
			SWAP(psa,psb,ti), @@+ SWAP(psc,psd,ti);
	} else { /* Reverse order */
		const int bmd = psb-psd, cma=psc-psa;
		if ( normm(bmd) < normm(cma) )
			SWAP(a,b,ti), @@+ SWAP(c,d,ti),@@+
			SWAP(ca,cb,tcn), @@+ SWAP(cc,cd,tcn),@@+
			SWAP(psa,psb,ti), @@+ SWAP(psc,psd,ti);
	}
}

@@ It is convenient to know which parent comes first.  We'll write |u| and
|v| for the lefthand and righthand parents of |ca| and |cb|, respectively.

There are two kinds of tasks to perform.  The first involves
changing pointers at the boundary.  The second involves making changes
at each of the parent nodes in the reversed sequence.



@@<Flip the sequence of segments from $a$ to $c$@@>=
{
parent_node_t *u, *v;
if ( normp(psa+1)==psb ) u=cc->parent, v=ca->parent;
else u=ca->parent, v=cc->parent;
@@<Fix inbound and outbound sibling pointers@@>@@;
@@<Fix per-parent data in the reversal from |u| to |v|@@>@@;
}


@@ It is convenient to do the end-pointer manipulation first, 
while we have a good conceptual 
handle on where pointers point.  

We must fix outbound and inbound city sibling pointers at the boundaries,
and outbound and inbound parent sibling pointers at the boundaries.

This code looks hairy, but only because to fix
the city siblings, 
we have to handle all possible cases of reversal bits.  Also,
we need some double indirection in order to swap properly.

@@<Fix inbound and outbound sibling pointers@@>=
{
parent_node_t *tpn;
int ur, vr;
city_node_t **u_outbound, **v_outbound, **u_inbound, **v_inbound, *u_first, *v_last;
ur=u->reverse;
vr=v->reverse;
u_first = u->city_link[ur^CITY_LINK_HEAD];
v_last  = v->city_link[vr^CITY_LINK_TAIL];
u_outbound = u_first->link + (ur^LINK_PREV);
v_outbound = v_last ->link + (vr^LINK_NEXT);
u_inbound = (*u_outbound)->link + 
	((*u_outbound)->link[LINK_NEXT] == u_first? LINK_NEXT:LINK_PREV);
v_inbound = (*v_outbound)->link + 
	((*v_outbound)->link[LINK_NEXT] == v_last ? LINK_NEXT:LINK_PREV);
	
SWAP(*u_inbound,*v_inbound,tcn);	/* Fix inbound city sibling pointers */
SWAP(*u_outbound,*v_outbound,tcn);	/* Fix outbound city sibling pointers */
u->prev->next = v;	/* Fix inbound parent sibling pointers */
v->next->prev = u;
SWAP(u->prev,v->next,tpn); 	/* Fix outbound parent sibling pointers */
}

@@
This section updates the per-node data on the reversed sequence of parents.

The first job is to fix the parent sibling pointers.  This code is similar
to the code to reverse a sequence of city links a single segment, namely
swap the sibling link pointers: make a former successor into a predecessor
and vice versa.  The outbound parent sibling pointers have already been taken
care of.

The second and third jobs are to update the parent sequence numbers and to
flip the reversal bits.  

Updating sequence numbers is interesting only
because we must do all arithmetic modulo |num_groups|.  In the future,
we might take the extra care
so that we would never overflow the number representation, but this
would only be possible (on a 32-bit or better machine) 
if we had about two {\it billion\/} groups or more.  Instead, I've
opted for the go for simpler and hopefully faster code.

On the other hand, if you end up fixing this problem, beware that you must
be very careful.  For instance, it is not good enough to leave the
code as it is, changing |int|s to |unsigned int|s, because the definition
of |upvn| in effect becomes 
$$upvn = (\hbox{|u->seq + v->seq|}\ \mod\ m)\ \mod\ \hbox{|num_groups|}$$
which is {\it not\/}
in general equal to 
$$(\hbox{|u->seq + v->seq|}\ \mod\ \hbox{|num_groups|})\ \mod\ m.$$
For example, $(4\ \mod\ 3)\ \mod\ 2 = 1 \not=0=(4 \ \mod\ 2)\ \mod 3$.
Here I've taken $m$ to be 1 greater than the maximum value of an unsigned
integer, \eg\ $m=2^{32}$ on a 32-bit machine.


@@<Fix per-parent data in the reversal from |u| to |v|@@>=
{
	const int upv = u->seq+v->seq, upvn = normp(upv);	
	parent_node_t *i, *done=v->next, *tpn;
	errorif(upv<u->seq || upv<v->seq, 
		"We've overflowed the integer representation");
	for ( i=u; i != done; i=i->prev ) {	/* Yes, |prev|: it's the old |next| pointer. */
		const int new_seq = upvn - i->seq;
		i->seq = normm(new_seq);
		i->reverse ^= 1;
		SWAP(i->next,i->prev,tpn);
	}
}

@@*Debugging.
I'm not perfect.  This module didn't work the first time it compiled.
The following routines were created to help me debug the above routines.

Externally, they duplicate the interface using different names.  Internally,
they perform the queries and updates on both the array-based implementation
and the two-level-tree-based implementation, and check
the answers by comparing them.

Here's the interface.

@@<Exported subroutines@@>=
#if defined(TWOLEVEL_DEBUG)
void twolevel_debug_setup(const int num_vertices, const int start_seg_size);
void twolevel_debug_cleanup(void);
void twolevel_debug_set(int const *tour);
int twolevel_debug_next(int a);
int twolevel_debug_prev(int a);
int twolevel_debug_between(int a, int b, int c);
void twolevel_debug_flip(int a, int b, int c, int d);
#endif


@@ The main wrinkle in this debuggging code is that the 
|flip| operation of the ADT does not specify the orientation
of the resulting tour.    So we must remember and adjust for the
possibly different orientations of tours under both implementations.  
Variable |reverse| is non-zero when the orientations differ.  

We will
take the two-level tree as the base implementation and adjust calls
to the array implementation.  That is, the answers to the queries
will always be the answers given by the two-level tree; the answers
from the array queries will be adjusted before comparing.  I've made
this choice because I want to process exactly the same sequence of operations
as when the bare two-level tree implementation is used.  I want to
minimize the entropy during debugging.
@@^entropy@@>


@@<Module variables@@>=
#if defined(TWOLEVEL_DEBUG)
static int reverse;
#endif

@@ The setup, cleanup and tour-setting routines call the appropriate
routines in both modules, then set the reversal variables as appropriate.

@@<Subroutines@@>=
#if defined(TWOLEVEL_DEBUG)
void 
twolevel_debug_setup(const int num_vertices, const int start_seg_size) {
	array_setup(num_vertices);
	twolevel_setup(num_vertices,start_seg_size);
	using_two_representations=1;
}

void 
twolevel_debug_cleanup(void) {
	twolevel_cleanup();
	array_cleanup();
	using_two_representations=0;
}

void 
twolevel_debug_set(int const *tour) {
	if ( verbose >= 100) { printf("set\n"); }
	array_set(tour);
	twolevel_set(tour);
	reverse = array_next(0) != twolevel_next(0);
	if ( verbose >=200) 
		printf("\t\treverse %d == (an0=%d != tn0=%d)\n",reverse,array_next(0),
				twolevel_next(0));
	check_tours_match();
}
#endif

@@ We need the interface to the \module{ARRAY} module.
@@<Module headers@@>=
#if defined(TWOLEVEL_DEBUG)
#include "array.h"
#endif

@@ 
@@<Module variables@@>=
#if defined(TWOLEVEL_DEBUG)
static int using_two_representations;
#endif

@@ Here's a linear-time check on the tours.  
We check most of the internal 
consistency of the data structures, but not 
any of the |between| queries.  Checking all the possible |between|
queries would take cubic time!

I also check sequence numbers; they weren't being properly maintained.

@@<Subroutines@@>=
#if defined(TWOLEVEL_DEBUG)
static int
check_self_consistency(void)
{	int i, c,cnt, an_error=0, cnpt, gs,s,tail_s,lgs,ls, ng=0, ph,pt;
	const int first_city=
		parent_node[0].city_link[CITY_LINK_HEAD^parent_node[0].reverse]
		-city_node;
	parent_node_t *p;
	if (verbose >= 150)
		printf("Checking twolevel tour consistency, reverse==%d\n",reverse);
	tail_s = ls = city_node[twolevel_prev(first_city)].seq;
	lgs = parent_node[0].seq-1;
	for ( i=0, c=first_city; i<n && !an_error ; i++, c=cnt ) {
		if ( c == first_city && i > 0) {an_error=1;printf("Not a tour\n");}
		cnt = twolevel_next(c);
		cnpt = twolevel_prev(cnt);
		if ( cnpt != c ) {
			an_error=1;
			printf("twolevel next/prev inconsistent pos %d city %d next: %d, nextprev: %d\n",
				i,c,cnt,cnpt);
		}
		p = city_node[c].parent;
		if ( lgs != (gs=p->seq) ) { /* New segment. */
			ng++;	/* Number of groups encountered */
			if ( gs != ((lgs+1)%num_groups) ) {
				an_error=1;
				printf("Parent sequence numbers %d to %d not consecutive\n",
					lgs,gs);
			}
			lgs=gs;
			if ( tail_s != ls ) {
				an_error=1;
				printf("Seq of last city in segment %d doesn't seq of \"tail\"%d\n",
					ls,tail_s);
			}
			tail_s = p->city_link[CITY_LINK_TAIL ^ p->reverse]->seq;
			ls = city_node