twolevel.w,v

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

@@*Initializing the two-level tree.
The procedure |twolevel_set(int *t)| takes a tour represented
as an array of $n$ integers
and changes the two-level tree to
match it.  It balances the segment lengths as best as it can.

@@ We have a choice to make here.  We must decide how to map
cities to city nodes; this might have consequences on how 
well the cache will perform over the duration of the local optimization
phase of the program.
Our choices are as follows.
We can store the data for city $i$
in the $i$'th city node;
or, we could store the data for the city that appears in position $i$ of the
tour
in the $i$'th city node.

Storing the data for city $i$ in the $i$'th node carries any locality
properties of the initial city ordering over to this module.
For instance, if we preprocess the data using a space-filling curve to try 
to achieve better locality, then this first scheme carries that mapping into
the manipulations of the two-level tree.  In this way, this scheme
supports the efforts of that data rearrangement.

On the other hand, storing the data in the order that it appears in 
the tour might achieve better locality initially.  However, starting
tours such as the Greedy/Multiple Fragment tours have some very large
defects.  These defects both increase the tour size, but more importantly
for our current
discussion,  also likely decrease locality for the duration of the 
local optimization phase.

For this reason, I've chosen to map the data for city $i$ onto element
$i$ of the |city_node| array.

% Another point is that sequence numbers for cities are initialized to
%be balanced around zero.  That is, they range roughly over the
%interval from $-k$ to $k$, where $k$ is half the number of cities
%in the group.  This is done in the hope that sequence numbers will 
%not go near the overflow values of |INT_MAX| or |INT_MIN|, where
%we'd be in danger of making integer comparisons go wrong.
%
%Hang on: as long as indices are consecutive, we don't have to
%worry about integer comparisons going wrong anyway!  Unless, of course,
%there are more than |UINT_MAX| cities in each segment --- that's $2^{32}-1$
%on a 32-bit machine.  But we're only targeting for up to a million cities, so
%it shouldn't become a problem.

On another front, it turns out that \CEE/ doesn't guarantee that remainders
will be non-negative when the dividend is negative.  
But the pointer arithmetic for initializing sibling pointers
relies on this property.
So we must use
|(group-1+num_groups)%num_groups| in place of the plainer
|(group-1)%num_groups| in place of the plainer.  Similarly,
|(i-1+n)%n| replaces
|(i-1)%n|.

@@<Subroutines@@>=
void
twolevel_set(int const *tour) {
	int i, j, group, num_big_groups = n % num_groups, base_group_size = n/num_groups;
	for ( i=0, group=0; i<n; group++ ) {
		const int this_group_size = base_group_size + (group<num_big_groups);
		parent_node[group].head = city_node+tour[i];
		parent_node[group].tail = city_node+tour[i+this_group_size-1];
		parent_node[group].reverse = 0;
		parent_node[group].seq = group;
		parent_node[group].prev = parent_node + ((group-1+num_groups) % num_groups);
		parent_node[group].next = parent_node + ((group+1) % num_groups);
		for ( j=0; j<this_group_size ; j++, i++ ) {
			city_node[tour[i]].parent = parent_node+group;
			city_node[tour[i]].seq = j;
			city_node[tour[i]].prev = city_node+tour[(i-1+n)%n];
			city_node[tour[i]].next = city_node+tour[(i+1)%n];
		}
	}
	errorif( i != n || group != num_groups, "Bug in my 'rithmetic");
}

@@ We must export this routine.
@@<Exported subroutines@@>=
void twolevel_set(int const *tour);


@@*Variations in the data structure.  
There are few other potential variations on the data structure.  I hope
to implement a few of them, to study their tradeoffs.

@@ First, Fredman \etal.~maintian links between cities at the ends of 
segments.  An alternative is to keep distinct linked lists, one 
per segment.  That is, don't link ends of segments together at the
city node level.

This variation is  designed to  make |flip|s faster.  However, 
|prev| and |next| queries might be slower because they would have
to check for a |NULL| pointer.  The |between| query ought to be unaffected.

@@ Second, we might choose to never move a city from its node.  We would
change pointers instead.  When this is done, we don't have to explicitly
store the city number in each city node.  This should make things faster
because the cache will be able to store a larger portion of the
city node array.

I have to think about whether this change involves extra pointer 
manipulation.  I ought to measure the difference from the standard 
data structure.

@@ Third, perhaps we can borrow use some approximate counting techniques
such as are found in good implementations of the 
$<INSERT,DELETE,RANK>$ ADT.  See Dietz (INSERT REFERENCE).

@@ Fourth, perhaps we can soften the requirement for consecutive
numberings within segments and at the root level.  (Q: Is this
even done at all in the 
`practical' variation?)  Instead of wholesale renumbering, just
take an average of neighbours when inserting a node between two others.
Then ``rebalance'' when we run out of precision.  In the common case,
rebalancing ought to be rare.


@@*Queries.
We'll define the |prev| and |next| queries first.  The  |between| query
will follow later.

The only interesting thing about the |prev| and |next| queries is that
we must decide which link to follow by examining the parent's 
reversal bit.  We take advantage of the fact that the reversal bit
is either 0 or 1 by exclusive-or'ing it with an index into the
city node's link array.  This straightens the code path: the decision is
made 
via arithemtic and indexing, not by branching in the code.  That should
save time on heavily pipelined architectures.
@@^pipelined architectures@@>
In fact, this trick, er, {\it technique}, 
is the main motivation for declaring the link
pointers as an array.
@@^trick@@>
@@^technique@@>

Both queries take constant time.

@@<Subroutines@@>=
int
twolevel_next(int a) {
	const city_node_t *ca = city_node+a;
	return (ca->link[ LINK_NEXT ^ ca->parent->reverse ])-city_node;
}

int
twolevel_prev(int a) {
	const city_node_t *ca = city_node+a;
	return (ca->link[ LINK_PREV ^ ca->parent->reverse ])-city_node;
}

@@ We must export these routines.
@@<Exported subroutines@@>=
int twolevel_next(int a);
int twolevel_prev(int a);

@@ Query |between(a,b,c)| asks: ``In a forward traversal starting at $a$,
do we reach $b$ no later than $c$?''

The |between| query is drudgery in comparison with queries |prev|
and |next|.  It is a case analysis
on the cities' intra-segment number, 
parent numbers, and parent reversal
bits.
Fredman \etal.~(page 444) describe this query quite well, and I take
this description mostly from them.

If all three cities have distinct parents, then we need only compare
their parents' sequence numbers.

If all three cities are in the same segment, then we need only 
compare their intra-segment sequence numbers.  The actual comparison in
this case depends on the common parent's reversal bit.

If exactly two of the cities have the same parent, then the answer
depends only on their intra-segment sequence numbers and their common
parent's reversal bit.
The placement of the third city is irrelevent: it is outside the
segment.

In some of the conditional expressions in the following,
I use bitwise operators $\mid$ and |&|
in place of short-circuiting
operators |||| and |&&|.    
Specifically, I do this when both operands of the logical operators
are comparsions.  I expect that the extra time performing a possibly redundant 
integer comparison is worth the savings on 
pipelined architectures in not having to always take
a branch.
This trick also relies on \CEE/'s semantics
of providing consistent values of 0 and 1 for boolean expressions
evaluating to false and true values, respectively.
@@^pipelined architectures@@>
@@^trick@@>
@@^technique@@>

I also use bitwise exclusive-or, |^|, but there is no short-circuiting
logical exclusive-or operator in \CEE/;
exclusive-or must always evaluate both its arguments anyway.

We're also careful to use as sequence numbers the offsets from 
the beginning of the segment instead of the bare intra-segment sequence number.
This protects us in case some sequence of flips has forced the sequence
numbers to migrate and wrap around at one of the ends of the
integer range, |INT_MIN|$\ldots$|INT_MAX|.

This query takes constant time.

@@<Subroutines@@>=
int
twolevel_between(int a, int b, int c) {
	const city_node_t *ca = city_node+a, *cb= city_node+b, *cc = city_node+c;
	const parent_node_t *pa = ca->parent, *pb = cb->parent, *pc = cc->parent;
	const int sa=ca->seq-pa->head->seq, sb=cb->seq-pb->head->seq, sc=cc->seq-pc->head->seq;

	if ( pa == pb ) 
		if ( pa == pc) 	/* All in the same segment */
			if ( pa->reverse )
				 return (sa>=sb? ((sb>=sc)|(sc>sa)) : ((sb>=sc) & (sc>sa)));
			else return (sa<=sb? ((sb<=sc)|(sc<sa)) : ((sb<=sc) & (sc<sa)));
		else 	/* |a|, |b| in same, |c| elswhere */
		return (sa==sb) | (pa->reverse ^ (sa<sb));
	else 
		if ( pa==pc )	/* |a|, |c| in same, |b| elswhere */
			return (sa!=sc) & (pa->reverse ^ (sa>sc));
		else 
			if (pb==pc) /* |b|, |c| in same, |a| elswhere */
				return (sb==sc) | ((pb->reverse) ^ (sb<sc));
			else {/* All in different; much like all in the same, but use parent numbers. */
				const int psa=pa->seq, psb=pb->seq, psc=pc->seq;
				return (psa<=psb? ((psb<=psc)|(psc<psa)) : ((psb<=psc) & (psc<psa)));
			}
}

@@ We export this function.
@@<Exported subroutines@@>=
int twolevel_between(int a, int b, int c);


@@*Flipping.
Flipping is the most difficult operation, as it changes the data structure,
and must update pointers, reversal bits, and sequence numbers.

With explicit balancing, one can make a flip take $O(\sqrt{n})$ time;
Fredman \etal.~describe how.
However, I've implemented the practical variant they describe. First,
it uses
implicit balancing instead of explicit balancing.
Second, 
when 
the portion of the tour to be flipped
is no more than $3/4\times groupsize$ and lies entirely within one segment,
that segment is split and reversed using its parent's reversal bit.
This practical variant incurs a lower overhead, but its performance
may degenerate to $\Omega(n)$ time in the worst case.  One hopes that
this occurs only vary rarely.

@@ 
An important observation is that we can flip either the $b-d$ portion 
of the tour
or the $a-c$ portion.  We have this freedom because the tour may be 
arbitrarily oriented after
the flip.  
This simplifies the code somewhat.

Flipping a tour segment splits into three cases.  

The first case occurs if 
either the $b-d$ portion or the $a-c$ portion
lies entirely within one city list segment.

The second case occurs when both the $a-c$ and the $b-c$ parts are made
up only of entire segments.  That is, no segment contains both 
$a$ and $d$ nor both $c$ and $b$.

The third case covers the rest of the possibilities.  Fortunately, 
we may always rearrange things so that the second case applies.

(I used to define |SWAP| and |abs| globally using \CWEB/'s {\tt @@@@d}
construct, but that interferes with the IRIX 5.3 definition of |abs|
in \file{stdlib.h}.  So I moved these definitions to after that inclusion.)
@@^system dependencies@@>
@@^IRIX@@>

@@<Subroutines@@>=
#define SWAP(x,y,t)  ((t)=(x),(x)=(y),(y)=(t))
#define abs(x) ((x)<0?-(x):(x))
void
twolevel_flip(int a, int b, int c, int d) {
	city_node_t *ca = city_node+a, *cb= city_node+b, 
		*cc = city_node+c,*cd = city_node+d, *tcn;
	int psa=ca->parent->seq, psb=cb->parent->seq, 
		psc=cc->parent->seq, psd=cd->parent->seq, ti;
#if defined(TWOLEVEL_FLIP_CHECK_PRECONDITION)
	errorif( a != twolevel_next(b), "a != twolevel_next(b)" );
	errorif( d != twolevel_next(c), "d != twolevel_next(c)" );
#endif 

	@@<Handle case 1 of flipping, if it applies@@>@@;
	if (psa==psb) {		/* Arrange for case 2 to apply, part 1. */
		@@<Split the $a-b$ segment@@>@@;
	}
	if (psc==psd) {		/* Arrange for case 2 to apply, part 2. */
		@@<Split the $c-d$ segment@@>@@;
	}
	@@<Flip a sequence of segments@@>@@;
}

@@ We export this routine.
@@<Exported subroutines@@>=
void twolevel_flip(int a, int b, int c, int d);

@@ As I've said, 
the first case occurs if 
either the $b-d$ portion or the $a-c$ portion
lies entirely within one city list segment.

Note that both portions can't both be in single segments unless there
are at most two segments; that's a degenerate case.  I had given some
thought to adding code to flip the shorter portion until I realized
this fact.  So, following Bentley and McIlroy's advice (INSERT REFERENCE), 
I'm keeping
the code simple.

This code gets used in three places, as we'll see.  So I've split it off
and made it a separate section.

@@<Handle case 1 of flipping, if it applies@@>=
@@<Verbose: handle case 1@@>@@;
if ( psb==psd ) {	/* Case 1, but rename so that $a-c$ lies in one segment. */
	SWAP(ca,cb,tcn); @@+ SWAP(cc,cd,tcn);
	SWAP(psa,psb,ti);@@+  SWAP(psc,psd,ti);
	/* |a|, |b|, |c|, |d| are not used again. */
}
if ( psa==psc ) {	/* Case 1, really do it. */
	@@<Flip: |ca| and |cc| are in the same segment@@>@@;
	return;
}

@@ Let's start with the first case,
\ie, when both $a$ and $c$ lie in the same data structure segment.
This means one of two things.  Either the $a-c$ portion lies entirely within
this segment, or the $b-d$ portion lies within this segment and is
sandwiched by cities $a$ and $c$.

First we rename variables so that the list of cities from
|ca| to |cc| lies entirely within one segment and contains no part
of the |cb| to |cd| portion of the tour.
Then we flip the segment.


%Now, it may be the case that the $b-d$ segment lies entirely within this
%segment, and is sandwiched by cities |ca| and |cc|.  It would then be 
%wrong to flip the sequence of cities from |ca| through |cc| within this
%segment: that would have the effect of reversing the $b-d$ segment, {\it and\/}
%swapping cities $a$ and $c$.  So we must check for this condition and
%move |u| and |v| one city `inward' if it is the case.  Note that it is
%sufficient to check that |u->next| ...

@@<Flip: |ca| and |cc| are in the same segment@@>=
@@<Rename so that the |a-c| portion lies entirely in one segment@@>@@;
@@<Flip: |ca| to |cc| lies entirely within one segment@@>@@;

@@  
We rename variables so that variables
|ca| and |cc| point to the portion of the tour 
lying entirely within this segment.
This condition is convenient to check if we first rename variables so that
|ca| is to the left of |cc|.

Of course, if |ca==cc|, then the |ca| to |cc| portion of the tour is just
one city long, so there is no work to do.  We don't check for this condition
because we think it will be common. Rather, we do it so that checking
for the containment of the |cb| to |cd| portion easy: it becomes the
simple test |ca->next == cb|.

It turns out that we won't be needing the identities of cities |b| or |d|