decluster.w,v

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

routines.  Procedure |decluster_setup| returns the space into which
a minimum spanning tree can be written.
Procedure |decluster_cleanup_tree| frees that space, given its pointer.  
It is a separate
call because the tree is no longer needed after preprocessing, so
it can be freed early.  In fact, this module throws away its reference
to that storage, so the caller {\it must\/} free it with a
call to |decluster_cleanup_tree|.

Procedure |decluster_mst| builds a minimum spanning tree for the current
graph.  It is passed a pointer to a block of space into which it
should write the tree.  It returns the weight of the tree.

We also export the comparison function |decluster_edge_cmp| so that others
can build priority queues with the same results as here.

Sometimes we want to compute many minimum spanning trees with our own
modified cost function, and without taking advantage of $k$-d trees.
In that case, use |decluster_mst_custom|, which takes 
an $n$-node tree node array into which it should write the answer, two utility
arrays of length $n$, and the custom cost function.  It returns
the length of the resulting tree.

Procedure |decluster_preprocess| digests a given tree so that least
common ancestor queries can be answered quickly by |decluster_d|.  
It overwrites the tree it is given.  Digestion indeed!

Separating the minimum spanning tree building out into a separate call
allows the user to substitute their own minimum spanning tree.  This is
useful in those cases where special knowledge of the input allows them
to build the tree very quickly.

Function |decluster_d| computes cluster distance $d$ efficiently.  
It uses the data structures built by procedure |decluster_preprocess|.

@@ Some parts of the interface supports uses outside just computing cluster
distances.  

Procedure |decluster_topology_tree| returns a pointer to the topology tree
that is preprocessed to answer decluster distance queries.  
It contains much digested information about the structure of the minimum
spanning tree of the instance.

Ordinarily, the topology tree is discarded immediately after preprocessing.
Set variable |decluster_discard_topology_tree| to zero to preserve it
for further use.  In particular, if  |decluster_discard_topology_tree|
is non-zero, 
|decluster_topology_tree| will return |NULL| after |decluster_preprocess|
is called.


Procedure |decluster_print_tree| outputs the structure of an arbitrary
valid |decluster_tree_t|.


@@<Exported subroutines@@>=
decluster_tree_t *decluster_setup(int n);
void decluster_cleanup_tree(decluster_tree_t *T);
void decluster_cleanup(void);
length_t decluster_mst(tsp_instance_t *tsp_instance,decluster_tree_t *T);
length_t decluster_mst_custom(decluster_tree_t *T,int *from,
			length_t *dist,length_t (*cost)(int,int));
int decluster_edge_cmp(const void *a, const void *b);
void decluster_preprocess(decluster_tree_t *T);
length_t decluster_d(int u, int v);

decluster_tree_t *decluster_topology_tree(void);
void decluster_print_tree(FILE *out,decluster_tree_t const *t, const char *name);


@@ We will be using many routines from external libraries.  The interfaces
to those routines are described in the following headers.

@@<System headers@@>=
#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>

@@ The exported interface is contained in the \file{decluster.h} header file,
which has the following form.

@@(decluster.h@@>=
extern const char *decluster_rcs_id;
@@<Exported types@@>@@;
@@<Exported variables@@>@@;
@@<Exported subroutines@@>@@;

@@ To ensure consistency between the interface and the implementation,
we include our own header.
@@<Module headers@@>=
#include "decluster.h"


@@ 
Up front, we know we'll need interfaces to
the error checking and memory allocation
modules, the |length_t| type (from 
\module{length}), and the |cost| function and internal
access to the internals of the graph (from \module{read}).

@@<Early module headers@@>=
#include "error.h"
#include "memory.h"
#include "length.h"
#include "read.h"


@@*Setup and cleanup.
Procedure |decluster_setup| just 
allocates the resources required by this module,
and sets up some module-level convenience variables.  It returns space
into which a minimum spanning tree should be written.

@@<Subroutines@@>=
decluster_tree_t *
decluster_setup(int the_n) 
{
	decluster_tree_t *mst;
	n = the_n;
	@@<Allocate space for a spanning tree |mst|@@>@@;
	@@<Other setup code@@>@@;
	return mst;
}

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

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

@@*Minimum spanning trees.
We'll be manipulating spanning trees throughout this module.  We begin
by defining a spanning tree type.
A spanning tree on $n$ vertices has $n-1$ edges.  Each edge is a pair
of vertices.  
The tree is an array of edges, each with their associated cost.

@@<Exported types@@>=
typedef struct {
	int city[2];
	length_t cost;
} decluster_edge_t;

typedef struct {
	int n;
	decluster_edge_t *edge; /* An array of $n$ edges */
} decluster_tree_t;

@@ Since a tree is just an array of edges, we can allocate it in
one swoop.
@@<Allocate space for a spanning tree |mst|@@>=
mst = new_of(decluster_tree_t);
mst->n = n-1;
mst->edge = new_arr_of(decluster_edge_t,n-1);


@@ Freeing the tree is just as easy.  

@@<Subroutines@@>=
void 
decluster_cleanup_tree(decluster_tree_t *T)
{
	if (T) { /* Be a little defensive. */
		size_t r = (T->n)*sizeof(decluster_edge_t) + sizeof(decluster_tree_t);
		T->n=0;
		free_mem(T->edge);
		free_mem(T);
		mem_deduct(r);
	}
}

@@ 
We can choose among different minimum spanning tree algorithms according
to our knowledge of the underlying structure of the graph.
In
particular, if the graph is embedded in Euclidean space then we can use
Bentley's $k$-d trees for semidynamic point sets, implemented in module
\module{kdtree}, to speed things up
considerably.  Otherwise we say the graph is unstructured and we use
a general minimum spanning tree algorithm.

Our graphs are complete, \ie, there is an edge between each pair of
vertices.  Prim's MST algorithm is efficient for unstructured 
complete graphs because it can be programmed to examine each edge only once.

Prim's algorithm begins with an empty tree $T_0=(V,\{\})$, and designates
a seed vertex $v_0$.  It builds up a single component from $V_0=\{v_0\}$ to
the entire vertex set $V$.
At step $i$ it forms the new version of the component $V_{i+1}$ from 
$V_i$ by adding the shortest edge between a vertex $v\in V_i$ in the current
version of the
component and a vertex $v_{i+1}\in V\setminus V_i$ not in the current
version of the 
component.

Oddly enough, Prim is also good for the Euclidean case because it
builds up a single component from nothing to the entire tree.  We can
filter out many of the edges from consideration by using
the $k$-d tree's ability to hide cities.  We hide cities in the current 
version of the component, \ie, at step $i$ all vertices in $V_i$ are
hidden.
That way nearest neighbour searches will avoid edges that
would
create cycles.


@@<Subroutines@@>=
length_t 
decluster_mst(tsp_instance_t *tsp_instance,decluster_tree_t *T)
{
	extern int verbose;
	length_t total_len=0;
	errorif(T->n!=n-1,"decluster_mst: passed storage for a tree with %d"@@;
		" vertices instead of %d vertices", T->n, n-1);
	if ( E2_supports(tsp_instance) ) {
		@@<Build an MST using $k$-d trees@@>@@;
	} else { /* Use a plain MST algorithm. */
		int *from=new_arr_of(int,n);
		length_t *dist=new_arr_of(length_t,n);
		total_len = 
			decluster_mst_custom(T,from,dist,cost);
		free_mem(from);mem_deduct(n*sizeof(int));
		free_mem(dist);mem_deduct(n*sizeof(length_t));
	}
	return total_len;
}

@@ Let's implement a plain MST algorithm, one which uses only the
given |cost| function and the knowledge that all possible edges are
in the graph.
We use a generic version of Prim's algorithm, outlined in the previous
section.

Value |from| is a pointer to an array of |n| integers, 
and |dist| is a pointer to an array of |n| length values.

If vertex |i| is not in the component, then |from[i]| is the
vertex in the component that is closest to |i|, and
|dist[i]==cost(from[i],i)|.  Otherwise, |from[i]=-1|.

Variable |short_len| is the minimum distance between a vertex in
the component and a vertex outside the component.  Variable |short_to|
is the outside vertex closest to the component.

We use city 0 as the seed city.


@@<Subroutines@@>=
length_t
decluster_mst_custom(
	decluster_tree_t *T,int *from,length_t *dist,length_t (*cost)(int,int))
{
int i, next_edge, short_to, n; /* Override file-scope |n|. */
length_t short_len, total_len = 0;
from[0]=-1;
errorif(T==NULL || T->n<0,"Bug!");
n=T->n+1;
for (i=1,short_len=INFINITY,short_to=-1;i<n;i++) {
	from[i]=0; 
	dist[i]=cost(0,i);
	if ( short_len > dist[i] ) {
		short_len=dist[i];
		short_to=i;
	}
}
for ( next_edge=0; next_edge < n-1; next_edge++ ) {
	@@<Add the edge to |short_to|@@>@@;
	@@<Update distance arrays and |short_to|@@>@@;
}
return total_len;
}


@@
@@<Add the edge to |short_to|@@>=
if (verbose>=1000) printf("decluster_mst_plain: adding edge (%d,%d) "length_t_spec"\n",
	from[short_to],short_to,length_t_pcast(short_len));
T->edge[next_edge].city[0] = short_to;
T->edge[next_edge].city[1] = from[short_to];
T->edge[next_edge].cost = short_len;
total_len += short_len;
from[short_to]=-1;

@@ Now that |short_to| has been added to the component, update 
arrays |from| and |dist|.  That is, see if |short_to| is closer to
an outside vertex than any of the other inside vertices.

We'll also be updating |short_to|, so first we copy the value of 
|short_to| to |new_inside_city|.

Scanning through all the cities wastes effort   % half, in fact
because we touch cities
that are in the component.  We could save time by collecting these cities
in one spot.  But that would entail still other bookkeeping.  
Hopefully this simpler version is faster because of its simplicity.

This section completes the general version of  Prim's algorithm.

@@<Update distance arrays and |short_to|@@>=
{ 
const int new_inside_city=short_to;
short_len=INFINITY;
for (i=1;i<n;i++) {
	length_t d;
	if ( from[i] == -1 ) continue; /* City |i| is already in the component. */
	d=cost(new_inside_city,i);
	if ( d < dist[i] ) {
		dist[i] = d;
		from[i] = new_inside_city;
	}
	if ( dist[i]<short_len ) {
		short_len=dist[i];
		short_to=i;
	}
}
}

@@  Now we're ready to implement Prim's algorithm for geometric inputs.
Specifically, we'll use Bentley's $k$-d trees for semidynamic point sets.
Throughout the execution of the algorithm, cities in the component will
be hidden, and cities outside the component will be visible.
This allows us to eliminate the vast majority of edges from consideration.

(It might be nice to measure this.   I suspect that this scheme will work
very well for uniform instances, because the shortest bridge drawn from
the priority queue will probably be consistently short.
The scheme might not work so well for clustered instances.)

(Hmmm.  Check 
J.~L.~Bentley and J.~H.~Friedman, {\sl Fast Algorithms for Constructing
Minimal Spanning Trees in Coordinate Spaces}, IEEE Transactions on
Computers, C-27:2, pp.~97--105, 1978.   Also, 
J.~L.~Bentley and J.~B.~Saxe {\sl An Analysis of Two Heuristics for
theEuclidean Travelling Salesman Problem}, pp.~41--49, in
{\sl 18th Annual Allerton Conference on Communication, Control, and
Computing}, October 1980.)

We'll keep a priority queue |bridge_pq| of edges that have one endpoint
in the component (recorded in field |city[0]|) 
and the other endpoint outside of the component (recorded in field
|city[1]|).
We'll eventually have an entry for each of $n-1$ vertices, so we
preallocate an array |bridge| of $n$ of them.  (We don't know in advance which
vertex will be the last to enter the component.)

Initially the priority queue contains just the edge from vertex 0 to
its nearest neighbour.  I'm assuming that the $k$-d tree is initialized
to have no cities hidden.  We clean up after ourselves by unhiding all
of them at the end.

@@<Build an MST using $k$-d trees@@>=
{
decluster_edge_t *bridge=new_arr_of(decluster_edge_t,n);
pq_t *bridge_pq=pq_create_size(decluster_edge_cmp,n);
int next_edge; char *is_in_component=new_arr_of_zero(char,n);

errorif(bridge_pq==NULL,"Couldn't allocate a priority queue!");

E2_hide(0);
bridge[0].city[0]=0;
bridge[0].city[1]=E2_nn(0);
bridge[0].cost = cost(0,bridge[0].city[1]);
is_in_component[0]=1;
pq_insert(bridge_pq,&bridge[0]);

for (next_edge=0 ; next_edge < n-1 ; next_edge++ ) {
	int in,out;
	decluster_edge_t *short_bridge;
	while (1) { /* Get a valid bridge. */
		short_bridge = pq_delete_min(bridge_pq);
		in=short_bridge->city[0]; 
			/* The city already in the component, and therefore hidden. */
		out=short_bridge->city[1]; /* The city possibly outside the component. */
		if ( !is_in_component[out] ) break;
		bridge[in].city[1]=E2_nn(in);
		bridge[in].cost=cost(in,bridge[in].city[1]);
		pq_insert(bridge_pq,bridge+in);
	}
	T->edge[next_edge] = *short_bridge;
	total_len += short_bridge->cost;
	E2_hide(out);
	is_in_component[out]=1;

	bridge[in].city[1]=E2_nn(in);
	bridge[in].cost=cost(in,bridge[in].city[1]);
	pq_insert(bridge_pq,bridge+in);
	bridge[out].city[0]=out;
	bridge[out].city[1]=E2_nn(out);
	bridge[out].cost=cost(out,bridge[out].city[1]);
	pq_insert(bridge_pq,bridge+out);
}
pq_destroy(bridge_pq);
free_mem(is_in_component);mem_deduct(n*sizeof(char));
free_mem(bridge);mem_deduct(n*sizeof(decluster_edge_t));
E2_unhide_all();
}

@@ We need to import the interfaces to priority queues and to $k$-d trees.
@@<Module headers@@>=
#include "pq.h"
#include "kdtree.h"

@@ One last thing we need is the comparison function |decluster_edge_cmp| for edges
so the priority queue can give edges to us in the right order.
Function |decluster_edge_cmp| has the same interface as required by the \CEE/
library function |qsort|.  

We force determinacy by comparing pointer values
in case the lengths are the same.   
(I don't think this is necessary in the use in Prim's algorithm above,
but it can't hurt.)

@@<Edge comparison function@@>=
int
decluster_edge_cmp(const void *a, const void *b) 
{
    length_t len_diff = ((const decluster_edge_t *)a)->cost 
						- ((const decluster_edge_t *)b)->cost;
    return len_diff < 0 ? -1 : 
		(len_diff>0 ? 
		 1 :
	     ((int)(((const decluster_edge_t *)a)-((const decluster_edge_t*)b))) );
}

@@  I also use this function in the test routines, but it is local to
both modules.

@@<Subroutines@@>=
@@<Edge comparison function@@>@@;


@@*Topology trees.
Now, given a minimum spanning tree $T$ such as the one constructed in the
previous sections, we'd like to create the related topology tree $T'$ as