jbmr.w,v

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

@@*Table-driven flips.
Now, if you look at the case analysis of April 22, 1996, you will see that there
is no discernable pattern to the flips required to effect the changes
in the initial 2/3/4-change.  Because of this, I use a table-driven 
approach to encoding these flips.  Instead of using in-line code implementing
the case analysis, I'm just going to encode my handiwork.

Encoding my handiwork in a table is simpler, smaller, and maybe even faster.
(It might be faster because the branches are easier to predict, either by
the compiler, or by the processor.)

@@ I should mention that there are in fact many possible ways to effect the
required 2/3/4-change.  I am only recording one.  If it is infeasible, then
this 2/3/4-change is declared infeasible.  Note that this is a conservative
approximation, and that it is conceivable that an alternative implementation
of the same 2/3/4-change might be feasible.  

However, I have chosen not to look for alternative implementations for
the following reasons.  

First, this search may take up too much time, and
is not likely to be of much benefit.  A change is infeasible because it
repeats at least one city, and perhaps repeats many.  If we fail once, then
we are making changes in a tight corner, and it is likely that many alternative
flip sequences encoding that change 
will also fail.  Why go through the bother of enumerating
all these alternatives when a simpler change can probably effect just
as good a gain?

Second, it was difficult enough to find the effecting sequences of flips
in the first place.  In fact, in my March analysis I didn't even see
that the 2/3/4-change required special handling (\ie, different from
the deeper changes in the $\lambda$-change).  Writing a program to 
enumerate all possible flips implementing a given 2/3/4-change would be a 
task comparable to writing a program enumerating 
all the possible solutions to the Rubik's
cube.  I don't want to do that.

@@  Here is the encoding of the flips required to effect the various
cases.

Array
|scheme[s]| holds the sequence of indices into the |t| array that implement
scheme |s|.  This sequence is read from left to right in groups of four.
A group of indices $a,b,c,d$ means apply |tour_flip(t[a],t[b],t[c],t[d])|.

The last entry of scheme |s| is in |scheme_max[s]-1|.
Array |scheme_num_cities| records the number of cities involved in each of
the schemes, \ie, the largest index into |t|.

The translation between case numbers (from my April 22, 1996 notes) and scheme
numbers is given in the comments.

@@<Module variables@@>=
static int scheme[14][16] =@@+{@@;
{
1,2,5,6,
4,3,2,5 },
/* Scheme 0, case 1.1.1: \epsffile{scheme0.eps} */

{
1,2,6,5,
2,6,4,3,
1,5,4,6 },
/* Scheme 1, case 1.1.2: \epsffile{scheme1.eps} */ 

{
5,6,3,4,
1,2,6,3,
6,2,8,7,
1,3,2,8 },
/* Scheme 2, case 1.2.1: \epsffile{scheme2.eps}*/ 

{
5,6,3,4,
8,7,6,3,
1,2,3,8 },
/* Scheme 3, case 1.2.2: \epsffile{scheme3.eps}*/ 

{
1,2,3,4 },
/* Scheme 4, case 2: \epsffile{scheme4.eps}*/ 

{
1,2,3,4,
1,4,6,5,
6,4,8,7,
1,5,4,8 },
/* Scheme 5, case 2.1.1.1: \epsffile{scheme5.eps}*/ 

{
1,2,3,4,
6,5,8,7,
1,4,5,8 },
/* Scheme 6, case 2.1.1.2: \epsffile{scheme6.eps}*/ 

{
1,2,3,4,
1,4,5,6 },
/* Scheme 7, case 2.1.2: \epsffile{scheme7.eps}*/ 

{
1,2,5,6,
5,2,3,4 },
/* Scheme 8, case 2.2.1: \epsffile{scheme8.eps}*/ 

{
6,5,8,7,
4,3,8,5,
1,2,3,8 },
/* Scheme 9, case 2.2.2.1.1: \epsffile{scheme9.eps}*/ 

{
1,2,8,7,
1,7,6,5,
1,5,2,8,
4,3,2,5 },
/* Scheme 10, case 2.2.2.1.2: \epsffile{scheme10.eps}*/ 

{
6,5,4,3,
6,3,8,7,
1,2,3,8 },
/* Scheme 11, case 2.2.2.2.1: \epsffile{scheme11.eps}*/ 

{
6,5,8,7,
1,2,5,8,
5,2,3,4 },
/* Scheme 12, case 2.2.2.2.2: \epsffile{scheme12.eps}*/ 

{-1}	/* Scheme 13, a generic change, so this entry is unused. */
}@@+;

static int scheme_max[14] = { 8,12,16,12,4,16,12,8,8,12,16,12,12,0 };
static int scheme_num_cities[14] = { 6,6,8,8,4,8,8,6,6,8,8,8,8,0 };

@@ We also need a scheme selection variable |scheme_id|.
It is declared local to the |jbmr_run| procedure, because it should not
be shared with other threads.

We determine the final value of |scheme_id| by performing the case
analysis as in my notes, all the while keeping |scheme_id| at value |-1|.  
During this determination, we use an 
array |base_scheme| indexed in the same way as |t|.
The value of |base_scheme[k]| is as 
tight a lower bound on the
final value of |scheme_id| as we can manage given the selections made
up to and including city |t[k]|.  Since schemes involve at most eight
cities, the values of index $k$ are bounded from above by 8.

When the analysis is finished, 
we set |scheme_id| and implement the changes required
for that scheme.  We maintain the invariant that |scheme_id == -1| if and
only if no intial 2/3/4-change has been implemented.

@@<|jbmr_run| variables@@>=
int scheme_id, base_scheme[9];


@@*Performing the search.
At the very least, this code assumes that we have allocated at least
nine (8+1) entries in the |t| array.

We examine both neighbours of |t[1]| as possibilities for |t[2]|.  
Because it may lead to a better initial gain, we first try the farther
tour neighbour of |t[1]|.  This is a greedy criterion; Lin and Kernighan
use this kind of criterion to prefer one neighbour of |t[7]| over the other.

@@<Search for an improving sequence starting at |dirty|@@>=
{
	int t1_n[2], t1_i;
	length_t t1_l[2];
	@@<Verbose: announce start of search at |dirty|@@>@@;
	t[1] = dirty;
	t1_n[0] = tour_prev(t[1]);
	t1_n[1] = tour_next(t[1]);
	t1_l[0] = cost(t1_n[0],t[1]);
	t1_l[1] = cost(t1_n[1],t[1]);
#if JBMR_FARTHER_T1_FIRST
	if ( t1_l[0] < t1_l[1] ) {
		int tmp;@@+
		length_t tmp_l;
		tmp = t1_n[0];@@+t1_n[0] = t1_n[1];@@+t1_n[1] = tmp;
		tmp_l = t1_l[0];@@+t1_l[0]=t1_l[1];@@+t1_l[1]=tmp_l;
	}
#endif

	@@<Initialize the bookkeeping variables@@>@@;
	put_city(1);
	for ( t1_i = 0 ; t1_i < 2 && more_backtracking; t1_i++ ) {
		t[2] = t1_n[t1_i];
		@@<Verbose: new |(t1,t2)|@@>@@;
#if !SPLIT_GAIN_VAR
		cum_gain = t1_l[t1_i];
#else
		cum_gain_pos = t1_l[t1_i];
		cum_gain_neg = 0;
#endif
		put_city(2);
		@@<Search from |t[2]|@@>@@;
	}
	if ( best_gain > 0 ) {
		incumbent_len -= best_gain;
		@@<Verbose: announce improvement by |best_gain|@@>@@;
		@@<Check milestone@@>@@;
		@@<Set the |instance_epsilon| slop value@@>@@;
	}
}

@@  We initialize the cumulative gain components |cum_gain_pos| and
|cum_gain_neg| so that GCC's dataflow analysis doesn't complain when
the debugging code is allowed.

@@<Initialize the bookkeeping variables@@>=
#if !SPLIT_GAIN_VAR
cum_gain = 0;
#else
cum_gain_pos = cum_gain_neg = 0;
#endif
best_gain = 0; best_two_i = 0; best_exit_a = best_exit_b = -1;
more_backtracking = 1; scheme_id = best_scheme_id = -1;

@@*Searching at level 0.  
As you may have guessed already (for instance from the presence of the
variable |more_backtracking|), the Lin-Kernighan algorithms performs some
backtracking.  Now, it's an odd beast in that it neither performs arbitrary
depth backtracking---for instance as every branch and bound algorithm does---
nor does it perform zero backtracking, as most pure greedy algorithms do.
If it chose either of these strategies, then our code 
would be simplified because of the regularity of these search paradigms.

Instead Lin-Kernighan is a hybrid, performing backtracking on the 
first few levels, and
probing in a purely greedy manner at deeper levels.
Fortunately, ``first few'' really means
``is bounded by a constant''.  So it will be convenient for us to build
and iterate through 
a fixed amount of state using ordinary nested loops.

Just to make things more interesting, 
the combinatorics of $k$-changes forces
each level of the search to have its own oddities, so each must be coded
differently in places.
Yet much of the code has the same overall structure, and we would like
to reflect that single architecture by using a single collection of 
named sections.  

Now, sit down, because I'm about to tell you how I managed to reconcile
these conflicting coding design goals.  Believe me, you won't like it.

Here it is.  
Most of the code is in a single collection of named sections.   The
top level of these is replicated once for each of the four levels of 
backtracking.  However, the oddities of each are separated out via
conditional compilation using the macro |BL|, which stands for {\sl
backtracking level}, and takes on one of the values 0, 1, 2, or 3.  

Actually, this code isn't nearly as ugly as it used to be.  

@@
The search from |t[2]| is a prototype for the search from any city
in an even-numbered position in |t|.    They all share the same style:
first generate the list of all eligible moves from this point, and then
iterate through them.

INSERT MATERIAL ABOUT DECLUSTERING.

@@<Search from |t[2]|@@>=
two_i = 2;
@@<Verbose: update |probe_depth|@@>@@;
#define BL 0
{
#if JBMR_DECLUSTER_IN_ELIGIBILITY_TEST
const length_t cluster_dist=decluster_d(t[1],t[2]);
#else
const length_t cluster_dist=0;
#endif
@@<Update |best_gain| and compose a list of eligible moves@@>@@;
}
@@<Sort |e[BL]|@@>@@;
#undef BL
@@<Backtrack at level 0@@>@@;

@@ We employ lookahead as described in section 2B of the original LK
paper.  
It is also adopted by JBMR.  

This lookahead criterion states that we must choose the next two cities
so that
$$\cost(t[2i],t[2i+1]) - \cost(t[2i+1],t[2i+2])$$
is minimized.  That is, we go to the shortest succeeding Hamiltonian Path.

If we are in an initial special case, then we may be doing backtracking.
So we need to create a list
of such candidate moves; it will be convenient to sort them according to
the above criterion.

In such an array, each move's entry records the proposed choices for
|t[2i+1]| and |t[2i+2]|, 
the scheme ID (if applicable), and the net gain as expressed in 
|gain_pos| and |gain_neg|, the positive and negative parts of the gain.

@@<Module types@@>=
typedef struct {
	length_t gain_for_comparison;
#if JBMR_DECLUSTER_IN_ELIGIBILITY_TEST || JBMR_DECLUSTER_IN_GREEDY 
	length_t cluster_dist;
#endif
#if !SPLIT_GAIN_VAR
	length_t gain;
#else /* |SPLIT_GAIN_VAR| */
	length_t gain_pos, gain_neg;
#endif /* |SPLIT_GAIN_VAR| */
	int t2ip1, t2ip2, scheme_id;
} eligible_t;

@@ Each of the four possible search levels --- 
for searches at |two_i| being equal to
2,
4, 6, or higher --- has its own array of candidates.
These are stored in
the |e| array, in positions 0, 1, 2, and 3, respectively.  
% Only the first entry on the deepest level is used, because there is no
% backtracking at the deepest level.

@@<|jbmr_run| variables@@>=
eligible_t *(e[4]);

@@ Each entry must contain space for
least twice |nn_max_bound| entries because for each choice of |t[2i+1]|, there
are up to two choices for |t[2i+2]|.

@@<Allocate |jbmr_run| sets and arrays@@>=
e[0] = new_arr_of(eligible_t,nn_max_bound*2);
e[1] = new_arr_of(eligible_t,nn_max_bound*2);
e[2] = new_arr_of(eligible_t,nn_max_bound*2);
e[3] = new_arr_of(eligible_t,nn_max_bound*2);

@@
@@<Deallocate |jbmr_run| sets and arrays@@>=
free_mem(e[0]);
free_mem(e[1]);
free_mem(e[2]);
free_mem(e[3]);

@@ Now, we may have varying numbers of eligible sequences at each of the levels.
In particular, we may have fewer than |nn_max_bound| valid entries in the 
0, 1, or 2
entries 
of |e|.  Array |en| holds
these counts, and is indexed in the same way as |e|.

@@<|jbmr_run| variables@@>=
int en[4];

@@ During the actual backtracking search, we also need to know where amongst
the |e| array we are right now.  So we need four cursors, one for each
backtracking level.  These are stored in the integer array |ec|.  

At any point in the search, only entries |ec[0]| through |ec[BL]| are valid.
Furthermore, if |ec[i]| is valid, then |0<=ec[i]<en[i]|.

@@<|jbmr_run| variables@@>=
int ec[4];


@@
Some of my experiments with |length_t==double| did not terminate. 
This happened even
after I split the gain variables (|cum_gain|, |cum_1|, and |cum_2|) into 
positive and negative parts in order to avoid catastrophic cancellation.
Examining a particular run, I noticed that this program eventually
got itself into
an endless loop, bouncing back and forth between 
two 3-changes every 25 tries at a new |t[1]|.  I think each of these
3-changes cancelled the other, but very small rounding errors in the 
caused the program to think that each was an improving 3-change.  Clearly,
the program was wrong about one of them.

My proposed fix is to use a positive slop value, |instance_epsilon|, and prune
any sequence of moves with a smaller cumulative gain than this.
Now, I call this value |instance_epsilon| to distinguish it from the 
machine epsilon, |LENGTH_MACHINE_EPSILON|. (|LENGTH_MACHINE_EPSILON| is
is the smallest number that can be added to 1 so that the result
is distinguishable from 1.)  

Actually, we prune any sequence of moves with a cumulative gain less
than the current best tour gain plus the instance epsilon.  Since we're
always testing against |best_gain+instance_epsilon|, I define a new
variable |best_gain_with_slop| that is just this sum.   This saves
us a floating point addition every time we test against this value.

%We initialize |instance_epsilon| so that GCC's dataflow analysis will 
%be silenced on this point.  We use a value of 1, which if kept, would
%result in no
@@<|jbmr_run| variables@@>=
#if !(LENGTH_TYPE_IS_EXACT)
length_t instance_epsilon=incumbent_len * LENGTH_MACHINE_EPSILON;
length_t best_gain_with_slop;
#endif

@@ We must initialize |best_gain_with_slop| at the same time that we
intialize the bookkeeping variables.

@@<Initialize the bookkeeping variables@@>=
#if !(LENGTH_TYPE_IS_EXACT)
best_gain_with_slop = instance_epsilon;
#endif



@@
However, 
|instance_epsilon| and |LENGTH_MACHINE_EPSILON| are related.  
Since we eventually subtract each tour gain from the incumbent length,
we require that the tour gain be large enough so that this subtraction
makes a difference.  (Excuse the pun.)  
@@^pun@@>
So the minimum value of |instance_epsilon| is 
$$\hbox{|instance_epsilon|} = 
\hbox{|incumbent_len|} \times \hbox{|LENGTH_MACHINE_EPSILON|}.$$
This gets set upon entry to |jbmr_run| and whenever we commit to
a net tour length reduction.

@@<Set the |instance_epsilon| slop value@@>=
#if !(LENGTH_TYPE_IS_EXACT)
instance_epsilon = incumbent_len * LENGTH_MACHINE_EPSILON;
#endif

@@ The code looks really awful when we have two quasi-independent
dimensions of parameterization: exact vs.~inexact, and split vs.~joined.
The |CAREFUL_OP(LHS,OP,RHS)| takes out some of the visual clutter.   
It implements numerical comparison |OP| on left-hand side |LHS| and
right-hand side |RHS|.  The twist is that the |LHS| is a variable 
that might be split into positive and negative parts.  We try to avoid
catastrophic cancellation by adding negative parts to the right hand side.
The right hand side must end in a variable to which we may add
|_with_slop|.  Right now the only variable with that ending
is |best_gain_with_slop|.

The double hash |##| is the ANSI C token pasting
operator.  So, for example, if the length type is |double|, an inexact
type, and |SPLIT_GAIN_VAR| is true, then 
$$\hbox{|@@[CAREFUL_OP(cum_gain,<,best_gain)@@]|}$$ translates into
$$\hbox{|((cum_gain_pos)<(best_gain_with_slop)+(cum_gain_neg))|.}$$

@@<Module definitions@@>=
#if LENGTH_TYPE_IS_EXACT
#define CAREFUL_OP(LHS,OP,RHS) ((LHS) OP (RHS))
#elif SPLIT_GAIN_VAR
#define CAREFUL_OP(LHS,OP,RHS) ((LHS##_pos) OP ((RHS##_with_slop)+(LHS##_neg)))
#else
#define CAREFUL_OP(LHS,OP,RHS) ((LHS) OP (RHS##_with_slop))
#endif




@@  We use a few bookkeeping variables when we update |best_gain| and compose 
an eligible list.

The neighbour list for city |i| is computed by |nn_list(i,&nn_bound)|.
Such a list an array
of |nn_bound| city numbers,
sorted in ascending order of
distance from city |i|;
|nn_bound| is bounded above by |nn_max_bound|.

Cities |t2ip1| and |t2ip2| are the candidates for |t[two_i+1]| and 
|t[two_i+2]|, respectively.  Lengths |cum_1| and |cum_2|
are the cumulative gains when we append |t2ip1| and |t2ip2| to the |t| list,
respectively.
When |length_t| is an inexact type, these are represented by their positive
and negative parts |cum_1_pos| and |cum_1_neg|, and 
|cum_2_pos| and |cum_2_neg|, respectively.  

City |t2ip1| is chosen from the nearest neighbour list of |t[two_i]|.
Since nearest neighbour lists are sorted, we stop as soon as 
the cumulative gain d