jbmr.w,v

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

#include "read.h"
#include "nn.h"
#include "lk.h"


@@*Setup and cleanup.
Procedure |jbmr_setup| just allocates the resources required by this module,
and sets up some module-level convenience variables.

@@<Subroutines@@>=
void
jbmr_setup(int the_n) {
	n = the_n;
	@@<Other setup code@@>@@;
}

@@
@@<Exported subroutines@@>=
void jbmr_setup(int the_n);

@@ We should declare |n|.
@@<Module variables@@>=
static int n;	/* The number of cities. */

@@ Deallocation is simple (especially with named sections!).
@@<Subroutines@@>=
void 
jbmr_cleanup(void)
{
	n=0;
	@@<Other cleanup code@@>@@;
}

@@
@@<Exported subroutines@@>=
void jbmr_cleanup(void);

@@
We'll maintain a set of ``dirty'' cities, \ie, those cities 
whose edges have been modified since the last time they were declared
locally optimal.  

Bentley invented this notion (as near as I can tell)
in the paper ``Fast algorithms for the geometric TSP'' (check title).
However, he couched the idea in terms of a ``don't-look'' bit on each
city.
Johnson, Bentley, McGeoch, and Rothberg adopt this notion,
with the same name.  Some people call the set of dirty cities the ``active''
list.  That's probably best, but is too late for this code.

Notice the ``don't look'' bit 
is in the reverse sense in which I'm thinking about this
problem.  That is, Bentley's bit is {\it off} when there is work to
be done.  I prefer to think of it in the positive terms.  That is,
there is an entry in the dictionary for every city for which there
is work to be done.  This is inspired by the notion of 
page table ``dirty'' bits
from operating systems work.


@@ I originally used a dictionary for this structure.  
The other option, a last-in-first-out queue, is described below.

@@<|jbmr_run| variables@@>=
dirty_set_t *dirty_set;

@@
@@<Module headers@@>=
#include "dirty.h"

@@ It will be useful to have a macro for making a city dirty.

@@d mark_dirty(CITY) (dirty_add(dirty_set,(CITY)))


@@ Initially the dirty set contains all the cities.
Some implementations fill the set in random order, so we pass on a random 
number generator seed.  
The seed is a next value from the random stream |random_stream|, 
which itself is passed into |jbmr_run|.

@@<Create the dirty set@@>=
dirty_set = dirty_create(n,1,prng_unif_int(random_stream,n),__FILE__,__LINE__);

@@ We must not forget to clean up the set.
@@<Deallocate |jbmr_run| sets and arrays@@>=
dirty_destroy(dirty_set);

@@ We'll also need to maintain a set of added edges.  JBMR never delete
an added edge.  In particular, this leads to the guarantee that we'll
never explore more than |n| moves for a particular improvement sequence.

Papadimitriou's variation (INSERT REFERENCE) 
dispenses with the added edges list and
instead never inserts a previously deleted edge.  This extends the
possible depth to roughly all |n(n-1)/2| (undirected) edges.  
This 
version of the Lin-Kernighan algorithm is the basis of a PLS-complete
problem.  I would like to implement this variation at some point in the
future.  The machinery is in place now.

I'll use  a dictionary to maintain the set of added edges.  We'll
allocate these as we need them, as this will likely be a small list.
So we don't allocate anything for this purpose on a module-level basis.

An alternative way of doing this is to forget messing around with 
the whole dictionary and to just compare with pairs of entries in the |t| array.
This not only saves space, but it may be faster in the common case, that is,
when the search is shallow.


@@*The Lin-Kernighan algorithm.
Ok.  Here's the big enchilada.

This is a local search routine.  It makes improvements until it can find
no more.  There are possible improvements to be made as long as there
are dirty cities.

We also implement the Iterated Lin-Kernighan scheme of Johnson
(INSERT REFERENCE) which was inspired by Martin-Otto-Felten's use of
chained local optimization.

@@<Subroutines@@>=
void 
jbmr_run (const int iterations, prng_t *random_stream) 
{
	@@<|jbmr_run| variables@@>@@;
	int iteration;

	@@<One-shot initialize@@>@@;
	@@<Allocate |jbmr_run| sets and arrays@@>@@;
	@@<Create the dirty set@@>@@;

	for ( iteration=0; iteration < iterations ; iteration++ ) {
		int dirty;
		@@<Per-iteration initialization@@>@@;
		while ( (dirty = dirty_remove(dirty_set)) >= 0 ) {
			@@<Search for an improving sequence starting at |dirty|@@>@@;
		}
		@@<Verbose: report end of LK step@@>@@;
		@@<Revert to the previous solution if it was better@@>@@;
		@@<If doing another iteration, perturb with a double-bridge@@>@@;
	}
	@@<Verbose: report termination of LK phase@@>@@;
	@@<Show final milestone@@>@@;
	@@<Deallocate |jbmr_run| sets and arrays@@>@@;
	@@<Verbose: print statistics@@>@@;
}

@@
@@<Exported subroutines@@>=
void jbmr_run(const int iterations,prng_t*random_stream);



@@ We'll need an array to remember the sequence of cities in the LK change.
I follow the notation introduced by Lin and Kernighan and use a 1-based
array named |t|.

For example, |t[1]| is the anchor to the Hamiltonian path, and |(t[1],t[2])|
is the first edge removed, and  |(t[2],t[3])| is the first edge added.
In general, edge |(t[2i-1],t[2i])| is the $i$'th edge removed, and
|(t[2i],t[2i+1])| is the $i$'th edge added.  
Edge 
|(t[2i],t[1])| is always a ``phantom'' edge, \ie, it is the edge which 
completes the Hamiltonian path to a Hamiltonian cycle, but is never drawn
in the descriptions of the LK algorithm.

This array 1-based and is treated as a growable array.
The array itself is |t|; the number of entries allocated to it is 
|t_max_alloc|.



@@<|jbmr_run| variables@@>=
int *t, t_max_alloc;

@@ 
The common case is for the 
|t| array to be
small. The average number of entries is roughly at most 
$14=8+6$ --- see the chapter
by Johnson and McGeoch.  Allocating 128 entries is a good compromise 
between saving on reallocations and saving space.

@@<Allocate |jbmr_run| sets and arrays@@>=
t_max_alloc = 128;
t = new_arr_of(int,t_max_alloc);

@@
@@<Deallocate |jbmr_run| sets and arrays@@>=
free_mem(t);

@@ We need a code fragment that will grow the array when necessary.

Doubling the size of the array each time we grow it 
requires at most constant amortized work per cell used. 

Under some  circumstances (known at compile-time),
we must notify others that we have moved the array.  
In particular, we must
update the nodes in the tabu splay tree.

@@<Make sure |t[two_i+2]| is valid@@>=
if ( two_i+3 >= t_max_alloc ) {
#if defined(TABU_SPLAY)
	int *old_t = t;
#endif
	do {
		t_max_alloc *= 2;	
	} while ( two_i+3 >= t_max_alloc );
	t = (int *)mem_realloc(t,sizeof(int)*t_max_alloc);
	@@<|t| has moved from |old_t| to |t|@@>@@;
}


@@ An extra tour query.  In my analysis, I found it useful to create the
following oriented tour query: |tour_inorder(a,b,c,d)| assumes that
|(a,b)| is an edge in the current tour, and answers the question ``in
a traversal of the tour beginning at |a| and going in the direction of
|b|, do we reach |c| no later than |d|? (We always count |a| as last.)''.

@@<Module subroutines@@>=
static int tour_inorder(int a, int b, int c, int d);
static int
tour_inorder(int a, int b, int c, int d) {
	if ( tour_next(a) == b ) return tour_between(b,c,d);
	else if ( tour_prev(a) == b ) return tour_between(d,c,b);
	else {
		@@<Debug: print the tour@@>@@;
		errorif(1,"Bad tour_inorder(%d,%d,%d,%d)\n"); 
		return -1;	/* Satisfy {\tt gcc -Wall}. */
	}
}

@@ We will need to record the cumulative improvement in the {\it Hamiltonian
path\/} following 
this sequence of moves.  This is recorded in the variable |cum_gain|.

Positive values mean net gains.

We'll also record, in |best_gain|, the gain made by the best {\it tour\/}
along this sequence (or along any two-city offshoot of this sequence).  
Note that |best_gain| and |cum_gain| are
not directly comparable because |cum_gain| doesn't include the cost of the
closing edge, \ie, the edge that joins the two ends of the Hamiltonian path.


%We use |possible_gain| as a working variable to compute the
%gain to be had by terminating the search now.
% |length_t possible_gain|;

@@ Some experiments that I run use a floating point type for |length_t|.
Unfortunately, sometimes these experiments don't terminate, especially
as problem sizes become larger, \eg, random Euclidean inputs with
$10^5$ or more cities.   (This happened the week of December 11, 1996.)

Some of the smaller runs that do manage to terminate end with 
a badly underestimated tour length; I check the computed ``incumbent'' length 
in |incumbent_len| against a freshly computed tour length.  This effect
increases with the problem size.

After
a little bit of thought, my best guess is that I run into the problem of
\term{catastrophic cancellation}.  
@@^catastrophic cancellation@@>
@@^precision@@>
@@^numerical analysis@@>
@@^floating point@@>
That is, when we compute the difference
between two nearly-equal numbers, most of the precision is lost. 
In our case, we are alternately adding and then subtracting edge costs
from the variable |cum_gain|.  We may even end up getting the wrong
sign in the end.  This, of course, easily leads to non-termination.

So, I split |cum_gain| into a positive part and a negative part,
|cum_gain_pos| and |cum_gain_neg|.  These three numbers are  are related
by the following equation:
$$
\hbox{|cum_gain|} = \hbox{|cum_gain_pos|} - \hbox{|cum_gain_neg|}.
$$
Actually, this allows a degree of freedom, and I constrain it somewhat
by stipulating
that both |cum_gain_pos| and |cum_gain_neg| should be non-negative numbers.
This is just like the construction of the set of integers as a set
of equivalence classes of pairs of natural numbers.  

In the present situation,
edge costs that increase the cumulative
gain (lengths of removed edges) are added to |cum_gain_pos|, and  edge costs
that decrease the cumulative
gain (lengths of added edges) are added to |cum_gain_neg|.  
I assume here that the |cost| function is non-negative.  In fact, it will
usually be positive-definite.

@@ On the other hand, I split these gain variables before I implemented the
idea of an \term{instance epsilon} (see below), which might alleviate the
same problem.  
@@^instance epsilon@@>
@@^epsilon, instance@@>
When the length type is inexact, the default is to split the gain variable
into positive and negative parts because that is safest.  
But this default may be overridden at compile time by defining the
symbol
|JBMR_REQUIRE_JOINED_GAIN_VAR|, say in file \file{lkconfig.h}.  
This symbol has no effect when the length type is exact.

Compile-time constant |LENGTH_TYPE_IS_EXACT|, provided by module
\module{LENGTH}
is non-zero when |length_t| is a type in which computations
are exact, \eg\ |length_t| is an integer or rational type.  
It is zero
when  there might be a loss of precision, \eg\ when |length_t| is
either |float| or |double|.

@@<Module definitions@@>=
#if LENGTH_TYPE_IS_EXACT || defined(JBMR_REQUIRE_JOINED_GAIN_VAR)
#define SPLIT_GAIN_VAR 0
#else
#define SPLIT_GAIN_VAR 1
#endif

@@ Now we can define the positive and negative cumulative gain variables.
@@<|jbmr_run| variables@@>=
#if SPLIT_GAIN_VAR
length_t cum_gain_pos, cum_gain_neg;
#endif

@@ Variable |best_gain| is always assigned the current value of the gain,
minus an edge cost.  So there is no need to break it up into positive
and negative parts.

@@<|jbmr_run| variables@@>=
length_t best_gain;

@@ Now, when using an integer type for |length_t|, we don't want to
overflow the representation.  So summing the alternating series in
a single variable is in fact decidedly {\it better\/} than splitting
it up.  So we need to retain the integer manipulations of |cum_gain|
and related variables anyway.  

We use the compile-time constant |SPLIT_GAIN_VAR| to distinguish the
two cases.

@@<|jbmr_run| variables@@>=
#if !SPLIT_GAIN_VAR
length_t cum_gain;
#endif

@@
The variables |best_two_i|, |best_exit_a|, and 
|best_exit_b| are defined by the following invariant:
{\medskip\leftskip=0.75in\rightskip=0.75in\noindent
The gain recorded in |best_gain| is put into effect 
by the changes encoded in |t[1]| through |t[best_two_i]|, followed
by the two cities |best_exit_a|, |best_exit_b|.\medskip}


The Lin-Kernighan algorithm includes many special cases for the initial 
2/3/4-change, and this affects the moves required for the exit.  If 
the best tour is found amongst these first few moves, then we must know
which sequence of moves to perform to get this best gain.  The identity
of this scheme is remembered in |best_scheme_id|, which takes on
values between 0 and 12, inclusive, for the initial 2/3/4-changes, and
13 for a generic change beyond that.  
See below
for a detailed discussion and encoding of these schemes.  It takes on the
value |-1| if it hasn't been defined yet.  

The value of |best_two_i| is always even.    
We maintain the invariant that  |best_scheme_id<0 == (best_gain == 0)|.
@@<|jbmr_run| variables@@>=
int best_two_i, best_exit_a, best_exit_b, best_scheme_id;


@@ I'll also use a boolean variable |more_backtracking| to signal when we should
keep on looking.  This variable will be non-zero when we should continue
looking for starting sequences for the
tabu search, and false otherwise.  In particular, its value is 
the negation of the boolean value of the phrase ``an improving move has
been implemented''.
This variable is used to escape from
various loops in this search.  

Now, |more_backtracking| will almost be synonymous with |best_gain==0|.  
However,
there is one place where they mean something different.  If the initial
segment of the |t| list is an improving (valid) 2-change, then |best_gain>0|
will hold, but we still want to look deeper, at least three levels deep.

So, you ask, why not postpone the setting of  |best_gain| until we enter
the non-backtracking portion of the tabu search?  The problem with this
is that, in case the 2-change is the best improvement we can make, we
don't end up pruning the search as much as we ought to.  That is, all
of the search, including the backtracking search for the initial 3-opt, ought
to be pruned by the 1-tree condition: |best_gain < cum_gain|, which 
is the same as |best_gain + cum_gain_neg < cum_gain_pos|.

This is why we need a seperate boolean variable.

In fact, an improving initial 2-change is not the only situtation where we
need to keep searching in spite of having an improvement already in hand.
{\it Lookahead\/} (see below) also requires us to have this kind of escape.

@@<|jbmr_run| variables@@>=
int more_backtracking;

@@ We will need a few bookkeeping variables.  

The integer |two_i| will be
an index in to the |t| array that points to the active endpoint, \ie, the
city at the end of the Hamiltonian path that is opposite to t[1], the anchor.
To put it another way, |(t[1],t[two_i])| is the phantom edge.

%The eight-entry (though 1-based) array |nn| is a set of read
%cursors into nearest neighbour lists.
%This is used to enumerate candidate neighbours.
% |int nn[9];|

@@<|jbmr_run| variables@@>=
int two_i;

@@*Constraints on the search. 
See my notes of March 7, 1996and April 22, 1996
for an explanation of the case analysis involved
in cycling through all the possible 3-changes (and sometimes 4-changes).

There are four kinds of constraints on the cities that we choose:
\emphpar{%
1.~The LK condition says that the cumulative 1-tree gain must be
greater than the improvement represented by the best tour gain discovered
so far.
In terms of our program variables, this is expressed by |cum_gain > best_gain|,
\ie, |cum_gain_pos > best_gain + cum_gain_neg|.
This is
the main stopping criterion, and it is checked both as we generate the
initial 2/3/4-change, and as we generate the cities on the
deep $\lambda$-change.  I will also call these ``deep'' changes 
{\sl generic}.
}

\emphpar{%
2.~The JBMR tabu edge condition says 
``never delete an added edge''.  This condition
will be checked as we generate both the initial 2/3/4-change 
and the deep portion of the $\lambda$-change.
}

\emphpar{%
3.~The specified edge changes must be feasible.  
That is, there must
always be a legal sequence of flips to get from the current tour to the
tour specified by the changes.  This is where the handling
of the initial 2/3/4-change differs from the rest of the $\lambda$-change.
}

\emphpar{%
The starting 2/3/4-change is checked for feasibility only after the entire
2/3/4-change is specified.  This is because the entire sequence of
flips required to effect the
change depends on the placement of even the last city in the sequence.
For example, see how the flips that implement cases 2.2.2.1.1 and 2.2.2.1.2
differ even in their first flip, even though only the position of the
last city changes between the two cases.  
}

\emphpar{%
Beyond the 2/3/4-change, \ie, the regular ``deep'' $\lambda$-changes
are checked as they are generated, just like for the other two conditions.
}

\emphpar{%
4.~I have added a ``stupid check'': during backtracking, never add the
most recently deleted edge.
This way we avoid covering the same ground twice while generating the
initial  2/3/4-change.  It also quickly eliminates some illegal choices,
specifically, in schemes 1, 8, and 11.
% THINK ABOUT THIS MORE CAREFULLY IN THE t[7] case, because we may want
%  to end up with  a different sequence ending?
%  Well, scheme 11 shows that it is useful at the t[7] level.
}