construct.w,v

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

@@
@@<Exported subroutines@@>=
length_t
construct(const int n, int *tour, const int heuristic, const long
heur_param, const long random_seed);

@@ We need the interface to the \module{LENGTH} module which 
defines the |length_t| type.
@@<Headers we need to define our own interface@@>=
#include "length.h"

@@ We also need the interface to the \module{ERROR} module, which defines
the error-checking convenience functions.  
@@<Module headers@@>=
#include "error.h"

@@ For a warmup, let's implement the canonical tour, \ie~ordering 0, 1, 2, \etc.
This us usually only good if the points have already been sorted according
to some space-filling curve.

@@<Heuristic alternatives@@>=
case CONSTRUCT_CANONICAL: {
	int i; length_t len;
	for (i=0;i<n;i++)tour[i]=i;

	len = cost(tour[0],tour[n-1]);
	for (i=1;i<n;i++) len += cost(tour[i-1],tour[i]);
	return len;
	break;	/* Redundant, I know$\ldots$ */
}

@@ We need to export the constant |CONSTRUCT_CANONICAL|.
@@<Exported constants@@>=
#define CONSTRUCT_CANONICAL 0




@@ That was too easy.  
Let's implement a random tour heuristic.  Random tours
are usually very bad, but can be computed very quickly.

We use |heur_param| as the random number generator's seed.

@@<Heuristic alternatives@@>=
case CONSTRUCT_RANDOM: {
	int i; length_t len;
	prng_t *random_stream; 

	for (i=0; i< n;i++) {
		tour[i] = i;
	}

	random_stream = prng_new(PRNG_DEFAULT,heur_param);
	for (i=0; i<n;i++) {
		const int next = prng_unif_int(random_stream,n-i);
		const int t = tour[next];@@+
		tour[next] = tour[n-1-i];@@+
		tour[n-1-i] = t;
	}
	prng_free(random_stream);

	len = cost(tour[0],tour[n-1]);
	for (i=1;i<n;i++) len += cost(tour[i-1],tour[i]);

	return len;
	break;	/* Redundant, I know$\ldots$ */
}


@@ We need to export the constant |CONSTRUCT_RANDOM|.
@@<Exported constants@@>=
#define CONSTRUCT_RANDOM 1


@@ And we must import the interface to the random number generators.
The |cost| function is provided by the \module{READ} module.

@@<Module headers@@>=
#include "prng.h"
#include "read.h"


@@*The Greedy heuristic.
Johnson \etal (INSERT REFERENCE), suggest using the Greedy 
heuristic to construct the
starting tour for Lin-Kernighan.   
Bentley (INSERT REFERENCE), calls it the Multiple Fragment heuristic, and
describes a fast algorithm for it on geometric inputs.

``Greedy'' has a few advantages over other
heuristics.  First, it is quickly computed, especially on geometric inputs.
Second, most of the edges it ends up using are very short, so the local
optimization phase doesn't have to fix up too much.  Third, the tours it
constructs have a few very large defects; the Lin-Kernighan cumulative
gain heuristic takes advantage of this to find very good tours.

%For now, I'll just use the best of the upper bounding heuristics that I
%used for the branch and bound procedure.  See the \module{upper} module.

@@ The greedy heuristic may be described as follows.  Begin with an
empty set of edges $T$.  Examine the edges in smallest to largest order.
If the edge in hand doesn't create a cycle in $T$ and doesn't create a vertex 
of
degree three in $T$, then that edge is added to $T$.  
Once $T$ has $n-1$ edges, join
the two ends of $T$ to make it a Hamiltonian cycle.  (The  first $n-1$ edges 
of $T$ form a Hamiltonian path.)

Notice the similarity to Kruskal's algorithm for constructing a minimum
@@^Kruskal's algorithm@@>
@@^minimum spanning tree@@>
spanning tree.  The main difference is that we never allow a vertex to
become of degree three, \ie, we never let our set of paths ``branch''.
Another way to put this is as follows. Let us say that a vertex is
\term{saturated} if it is of degree two.
@@^saturated vertex@@>
Then the Greedy heuristic is just
Kruskal's algorithm with the
proviso that no edge with a saturated endpoint is ever added,
together with
the addition of an $n$'th edge to complete the tour.

@@ To get different outputs on different runs, we can randomize the
above scheme. Instead of always taking the best available alternative,
we choose randomly among the best few alternatives.
JBMR's implementation picks the best alternative with probability 2/3,
and the second-best alternative the rest of the time.

@@ If implemented directly, the description just given for the Greedy heuristic
would be quite inefficient.
The implementation in this module 
follows Bentley's.

Bentley's innovations are threefold.  First,
he sorts and selects the edges lazily, using a priority queue containing
the shortest edge from each unsaturated vertex.   
Second, he uses $k$-dimensional trees to perform nearest neighbour queries
(where one is allowed to hide vertices).  See module \module{KDTREE} for
an implementation of this data structure.  Third, he uses an efficient
data structure to determine when an edge would create a subcycle.

Unfortunately, not all instances have the benefit of geometry, so we
have to add some cleverness in searching for next-nearest neighbours.
In the non-Euclidean case, this amounts to maintaining
extra bookkeeping about 
which cities remain unsaturated, and buffering multiple nearest neighbours
in the queue to save on repeated exhaustive searches of neigbours of a
single city.

Finally, much of this code for the TSP can be shared with the weighted
perfect matching code.  When there is a difference, the TSP-specific code is
delimited with |DO_TOUR|, and the matching code with |DO_MATCHING|.



@@<Heuristic alternatives@@>=
case CONSTRUCT_GREEDY: 
{
	length_t len=0;
	const int E2_case = E2_supports(tsp_instance);
#define DO_TOUR
	@@<Declare the data structures for Greedy@@>@@;
	@@<Build the data structures for Greedy.@@>@@;
	@@<Build the Greedy tour@@>@@;
	@@<Destroy the data structures for Greedy@@>@@;
#undef DO_TOUR
	return len;
break;
}

@@ The randomized case differs slightly.  We control its code by defining
|DO_RANDOM|.

@@<Heuristic alternatives@@>=
case CONSTRUCT_GREEDY_RANDOM: 
{
	length_t len=0;
	const int E2_case = E2_supports(tsp_instance);
#define DO_RANDOM
#define DO_TOUR
	@@<Declare the data structures for Greedy@@>@@;
	@@<Build the data structures for Greedy.@@>@@;
	@@<Build the Greedy tour@@>@@;
	@@<Destroy the data structures for Greedy@@>@@;
#undef DO_TOUR
#undef DO_RANDOM
	return len;
break;
}

@@ We need to export the constants |CONSTRUCT_GREEDY| and
|CONSTRUCT_GREEDY_RANDOM|.
@@<Exported constants@@>=
#define CONSTRUCT_GREEDY 		2
#define CONSTRUCT_GREEDY_RANDOM 3

@@ The first data structure we'll need is a priority queue of edges.  
%We'll
%record at most one edge per  vertex in the graph.  
%It will be convenient
%to store the edges in an array indexed by city number.  
Each entry records
the name of the other end of the edge, together with the length of that
edge.

I've put the length entry first because we will access it much more often
than the city entry.  This happens because the length entry is the
key in the priority queue.  The comparison function |cmp_pq_edge| 
defined below 
accesses the length entry through a pointer; if that field comes first,
then the offset from the pointer will be zero.  I assume that the \CEE/ 
compiler will be smart enough to elide the pointer arithmetic, making
the program that little bit faster.

@@<Module type definitions@@>=
typedef struct {
	length_t len;
	int this_end;
	int other_end;
} pq_edge_t;

@@ Here's a printing function that might be useful for debugging.
@@<Module subroutines@@>=
static void pq_edge_print(FILE *out, void *edgep);
static void
pq_edge_print(FILE *out, void *edgep)
{
	pq_edge_t *edge = edgep;
	if ( edge==NULL ) { fprintf(out,"(null)"); }
	else {
		fprintf(out,"{%p,%d,%d," length_t_spec "}",
			edge,edge->this_end,edge->other_end, length_t_pcast(edge->len)); 
	}
}

@@ In the Euclidean case, we keep in the queue at most one 
neighbour for each city.  An array is convenient for that, and
that is how Bentley defines it.
He says  ``|nn_link[i]| contains
the index of the nearest neighbour to $i$ that (when originally computed) is
not in the same fragment that contains |i| (though subsequent operations might
invalidate that condition).  These links are the edges that are eventually 
inserted into the tour.''  That is, he city from which the edge is
emanating is implicitly defined by pointer arithmetic.

In the non-Euclidean case, we cache multiple neighbours at a time,
so the array structure doesn't work well for us.  To keep the code
simpler, we use a dynamically allocated pooled of edges instead of just
an array.  This pool is called |nn_link_pool|.

The entries of |nn_link_pool| serve as the domain over which
the priority queue |edge_pq| is defined.
We will use the priority queue routines from module \module{PQ}.

@@<Declare the data structures for Greedy@@>=
pool_t *nn_link_pool = pool_create(sizeof(pq_edge_t),n);
pq_t *pq_edge = pq_create_size(cmp_pq_edge,n);

@@ We need to do a bit of error checking.
@@<Build the data structures for Greedy.@@>=
errorif(pq_edge==NULL,"Couldn't create the priority queue!");
pq_set_print_func(pq_edge,pq_edge_print);

@@ Deallocation is easy.
@@<Destroy the data structures for Greedy@@>=
pool_destroy(nn_link_pool);
pq_destroy(pq_edge);

@@ We need the interface to both the memory allocation routines, the
priority queue routines, and to the pool manager.

@@<Module headers@@>=
#include "memory.h"
#include "pq.h"
#include "pool.h"

@@ We also need the definition of |NULL|, which is in one of the following
header files.  Also, the |free_mem| macro needs the declaration of |free|
from \file{stdlib.h}.
@@<System headers@@>=
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define FIXINCLUDES_NEED_NRAND48 /* Interface in \file{prng.h} needs this. */
#include "fixincludes.h"
#undef FIXINCLUDES_NEED_NRAND48

@@ We've had to provide a comparison function for the priority queue routines.
The comparison function has the same interface as the one provided to
the |qsort| standard library function.  This one in particular compares
edge lengths.  We want our priority queue to return the shortest edge
first, so smaller edge lengths compare smaller.

Our function returns an |int|, and we don't know exactly what |length_t| really
is,  so we need to be a bit careful about how we translate the difference
between the lengths into a sign value.  Using an intermediate value should
be safe.  Note that both lengths should be non-negative, so their difference
should be smaller in magnitude than the greater of the two lengths, and
therefore within range of the |length_t| data type.
@@^overflow@@>

When the lengths compare equal, we must break the tie.  Otherwise we'd never
be able to insert two edges of the same length: the priority queue requires
unique keys.  We break ties by returning the difference in the pointers.

@@<Module subroutines@@>=
static int
cmp_pq_edge(const void *a, const void *b) {
	length_t len_diff = ((const pq_edge_t *)a)->len - ((const pq_edge_t *)b)->len;
	return len_diff < 0 ? -1 : (len_diff >0 ? 1 : 
		(int)(((const pq_edge_t *)a)-((const pq_edge_t*)b)) );
}

@@ The priority queue is initialized by filling in each city's entry
in the |nn_link_pool| collection, 
and then inserting that entry into the priority queue.
The entry for a city contains the that city's nearest neighbour and the
length to that neighbour.

%The distance function is assumed to be symmetric.  So if we have $(u,v)$
%in the queue, there is no point in 

We assume that the all the cities in the $k$-d search structure are
visible.

We separate out the non-Euclidean case.

@@<Build the data structures for Greedy.@@>=
if ( E2_case ) {
	int i;
	for (i=0;i<n;i++) {
		pq_edge_t *e = pool_alloc(nn_link_pool);
		e->this_end = i;
		e->other_end = E2_nn(i);
		e->len = cost(i,e->other_end);
		pq_insert(pq_edge,e);
		@@<Other initializations for city |i|@@>@@;
	}
} else {
	@@<Build the data structures for Greedy, non-Euclidean case@@>@@;
}

@@  We need the interface to $k$-d trees and the declaration of |tsp_instance|.
@@<Module headers@@>=
#include "lk.h"
#include "kdtree.h"

@@ In the non-Euclidean case, we can have multiple neighbours in the queue
for each city.  When those are exhausted, we get a fresh batch.  We
need a way of knowing when they are exhausted. For that we use an array
with entry
|farthest_in_queue[i]| being the city at the other end of the heaviest
edge in the queue and which had been added for city |i|.

@@<Declare the data structures for Greedy@@>=
int *farthest_in_queue=NULL;

@@
@@<Build the data structures for Greedy, non-Euclidean case@@>=
farthest_in_queue = new_arr_of(int,n);

@@
@@<Destroy the data structures for Greedy@@>=
if ( !E2_case ) { free_mem(farthest_in_queue); mem_deduct(sizeof(int)*n); }


@@ Furthermore, the non-Euclidean case needs an anolog to the ``hidden''
cities.  It is in fact more convenient to keep track of which cities
remain unsaturated.  Values |unsaturated[0]| through
|unsaturated[num_unsaturated-1]| are the cities which remain, well,
unsaturated.  We also support efficient updates to this array by
maintaining
an inverse mapping of the unsaturated cities: if city |i| 
is unsaturated, then $$\hbox{|unsaturated[inv_unsaturated[i]]==i|}.$$

@@<Declare the data structures for Greedy@@>=
int *unsaturated=NULL, num_unsaturated=0, *inv_unsaturated=NULL;

@@
@@<Build the data structures for Greedy, non-Euclidean case@@>=
unsaturated = new_arr_of(int,n);
inv_unsaturated = new_arr_of(int,n);
num_unsaturated=n;
{ int i;
	for (i=0;i<n;i++) inv_unsaturated[i]=unsaturated[i]=i;
}

@@
@@<Destroy the data structures for Greedy@@>=
if ( !E2_case ) { 
free_mem(unsaturated); mem_deduct(sizeof(int)*n); 
free_mem(inv_unsaturated); mem_deduct(sizeof(int)*n); 
}

@@ Here's how to make a city unsaturated.
Once we verify that it really was unsaturated to begin with, 
all we have to do is copy down the last unsaturated city down over its
entry in the |unsaturated| array.

Note that because the values in |inv_unsaturated| are  always non-negative,
we don't have to check for underflow of |num_unsaturated|.
@@<Mark city |c| as saturated@@>=
{ const int inv= inv_unsaturated[c];
	if (inv<num_unsaturated && unsaturated[inv]==c) {
		const int move_city = unsaturated[--num_unsaturated];
		unsaturated[inv]=move_city;
		inv_unsaturated[move_city]=inv;
		@@<Verbose: mark as saturated@@>@@;
	}
}

@@ Initializing the priority queue in the non-Euclidean case is done
city by city.  The initialization for each city is rather complicated,
and we shall reuse it, so we break it off into a separate section.

We'll be needing |sqrt_n| when building the neighbour lists, so we
define it once.

@@<Build the data structures for Greedy, non-Euclidean case@@>=
{
int i;
nn_work = new_arr_of(pq_edge_t,n);
for (i=0;i<n;i++) {
	const int x=i, not_me = -1;
	@@<Get fresh neighbours for |x|, excepting |not_me|@@>@@;
	@@<Other initializations for city |i|@@>@@;
}
}

@@
@@<Declare the data structures for Greedy@@>=
int sqrt_n = (int)sqrt((double)n);

@@  In the non-Euclidean case we don't know anything about the structure
of the distance matrix.  We have to scan all unsaturated cities for
possible neighbours.  However, for a large fraction of cities, we'll
be scanning their rows multiple times.  To save on these scans, 
we insert multiple entries into the queue for every scan we perform.

Of course, there is a tension between the space required to store the
queue, how quickly we can find those nearest neighbours, and how many
of them end up wasted (by the time we draw an edge from the queue,
the other end may have become saturated).

Initially 
we choose to put in up to |init_max_per_city_nn_cache| neighbours at a time
for each city.  We may make this a user paramter.  The default
value is 30, which is a complete guess.
After the number of saturated vertices declines, say to below $\sqrt{n}$ cities,
we can afford space-wise to put in all the rest of the possible edges.

In scanning the list of unsaturated cities, we 
skip over |x| itself and |not_me|.
(The |not_me| is used in the tour algorithm.   Trust me.)

We use the |select_range| procedure from module \module{DSORT}, which
requires the entries to be in an array.  We'll use |nn_work|, defined
below.

@@d init_max_per_city_nn_cache 30
@@d max(X,Y) ((X)<(Y)?(Y):(X))
@@d min(X,Y) ((X)>(Y)?(Y):(X))

@@<Get fresh neighbours for |x|, excepting |not_me|@@>=
if ( num_unsaturated > sqrt_n ) {
	@@<Get fresh neighbours by selecting from |unsaturated|@@>;
} else {
	@@<Put in all unsaturated neighbours of |x|, except for |not_me|@@>@@;
}

@@ 
@@<Put in all unsaturated neighbours of |x|, except for |not_me|@@>=
{ int i;
@@<Verbose: All neighbours@@>@@;
for (i=0;i<num_unsaturated;i++) {
	const int y = unsaturated[i];
	length_t len;
	if ( y==x || y==not_me ) continue;
	len = cost(x,y);
	@@<Put edge |(x,y)| with length |len| into the queue@@>@@;
}
}


@@
@@<Put edge |(x,y)| with length |len| into the queue@@>=
{
	pq_edge_t *e;
	e = pool_alloc(nn_link_pool);
	e->len = len;
	e->this_end = x;
	e->other_end = y;
	pq_insert(pq_edge,e);
	@@<Verbose: allocate new link@@>@@;
}


@@
@@<Get fresh neighbours by selecting from |unsaturated|@@>=
{
	int i,num_chosen, farthest_city;
	size_t w;
	length_t farthest_len;
	@@<Verbose: select fresh neighbours@@>@@;
	for (i=w=0;i<num_unsaturated;i++) {
		const int c=unsaturated[i];
		if ( c==x || c==not_me ) continue;
		nn_work[w].this_end=x;
		nn_work[w].other_end=c;
		nn_work[w].len=cost(x,c);
		w++;
	}
	num_chosen=min(w,init_max_per_city_nn_cache);
	errorif(num_chosen==0,"Bug!");
	(void)select_range(nn_work,w,sizeof(pq_edge_t),cmp_pq_edge,
		0,num_chosen,0 /* Don't sort. */);
	farthest_len = nn_work[num_chosen-1].len;
	farthest_city = nn_work[num_chosen-1].other_end;
	for ( i=0;i<num_chosen; i++) {
		pq_edge_t *e;