construct.w,v

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

		e=pool_alloc(nn_link_pool);
		*e = nn_work[i];
		pq_insert(pq_edge,e);
		if ( farthest_len < e->len ) {
			farthest_city = e->other_end;
			farthest_len = e->len; 
		}
	}
	farthest_in_queue[x]=farthest_city;
}

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


@@ We must declare the set |nn_work|.
@@<Declare the data structures for Greedy@@>=
pq_edge_t *nn_work=NULL;


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


@@ We need a place to record the tour as its being built up.  We use
an array |adj| indexed by city number.  Entry |adj[i]| records the set of
cities adjacent to |i| in the tour so far.    Because at any time
there are either zero, one or two cities adjacent to a given city, it
is sufficient for each entry to have a |count| field and a two-entry
array |city| of neighbours.

@@<Module type definitions@@>=
typedef struct {
	int count;	/* 0, 1, or 2 */
	int city[2];
} adj_entry_t;

@@ We allocate the |adj| when we declare it.

@@<Declare the data structures for Greedy@@>=
#if defined(DO_TOUR)
adj_entry_t *adj = new_arr_of(adj_entry_t,n);
#endif

@@ Deallocation is easy.
@@<Destroy the data structures for Greedy@@>=
#if defined(DO_TOUR)
free_mem(adj);
#endif

@@ Each vertex begins with no neighbours.
@@<Other initializations for city |i|@@>=
#if defined(DO_TOUR)
adj[i].count = 0;
#endif

@@ 
While
building the tour, we maintain a set of framgents.
Each fragment is
a (non-branching) path.  In general, each path has two tails, and these
are the cities to which we add links.  We avoid creating subcycles by
avoiding joining one endpoint of a fragment to its other end.  
It is therefore convenient to maintain an array, |tail|. When
city |i| is still a tail, entry |tail[i]|
records the index of the city at the other end of the fragment containing
|i|.

@@<Declare the data structures for Greedy@@>=
#if defined(DO_TOUR)
int *tail=new_arr_of(int,n);
#endif

@@
@@<Destroy the data structures for Greedy@@>=
#if defined(DO_TOUR)
free_mem(tail);mem_deduct(sizeof(int)*n);
#endif

@@ Each city begins in a fragment by itself.
@@<Other initializations for city |i|@@>=
#if defined(DO_TOUR)
tail[i] = i;
#endif

@@ Now we're ready for the main part of the Greedy heuristic, which builds
the Hamiltonian path.

Initially, each city is alone in its own fragment.  
We repeatedly take the cheapest remaining
valid link and add it to the tour.  Most of the code here is just
the bookkeeping involved in adding the chosen edge.

Once $n-1$ edges have been added, we connect the two remaining tails.
By sheer luck, these tails are stored in |tx| and |ty|.

@@<Build the Greedy tour@@>=
{@@+ int i, x, y, tx, ty;
errorif(n<3,"Only %d cities.  Can't build a tour.",n);
tx=ty=-1;	/* Satisfy GCC's dataflow analysis. */
for ( i=0;i<n-1;i++ ) {
	pq_edge_t *e;
	@@<Find a short valid edge $e=(x,y)$@@>@@;
	adj[x].city[adj[x].count++] = y;
	adj[y].city[adj[y].count++] = x;
	if ( adj[y].count == 2) { if (E2_case) E2_hide(y); else {const int c=y;
@@<Mark city |c| as saturated@@>}}
	if ( adj[x].count == 2) { if (E2_case) E2_hide(x); else {const int c=x;
@@<Mark city |c| as saturated@@>}}
	else { @@<Insert a new link for |x|, using |e|@@> }
	tx = tail[x];
	ty = tail[y];
	tail[tx]=ty;
	tail[ty]=tx;
}
adj[tx].city[adj[tx].count++] = ty; /* Complete the Hamiltonian cycle */
adj[ty].city[adj[ty].count++] = tx;
if (E2_case) E2_unhide_all();
len += cost(tx,ty);
}

@@ Here's where we pick a short remaining valid edge.
If we're in random mode, then get two eligible candidates
and choose randomly between them.

Early on, most edges will appear twice in the queue, and will come
out in succession.  For example, a shortest edge $(u,v)$ will be
quickly followed by its reverse-ordered reprentation $(v,u)$.
But things aren't very randomized if we flip a coin between
two represntations of the same thing!  So if we do get a reversed copy
of the first edge we withdraw, then we hold it off to the side while
getting a third edge (which can't be the same as the first), and
then re-insert the second edge.

@@<Find a short valid edge $e=(x,y)$@@>=
#if defined(DO_RANDOM)
{
	pq_edge_t *candidate[2];
	@@<Get one short valid edge or |NULL|@@>@@; @@+ candidate[0]=e;
	@@<Get one short valid edge or |NULL|@@>@@; 
	if ( e && e->this_end == candidate[0]->other_end @@| 
			&& e->other_end == candidate[0]->this_end ) {
		pq_edge_t *push_back = e;
		@@<Get one short valid edge or |NULL|@@>@@;  /* Get another. */
		pq_insert(pq_edge,push_back);  /* But put the second one back. */
	}
	candidate[1]=e;
	e = candidate[0];
	if ( candidate[1] ) { /* There is a choice to be made. */
		const int chosen = (prng_unif_int(random_stream,3) == 0); 
			/* Prob|[chosen==1]== 1/3| */
		e=candidate[chosen];
		pq_insert(pq_edge,candidate[1-chosen]); /* Put the other one back. */
	}
}
#endif

@@ In the deterministic case, just choose one.  Now isn't that much simpler!
@@<Find a short valid edge $e=(x,y)$@@>=
#if !defined(DO_RANDOM)
@@<Get one short valid edge or |NULL|@@>@@;
#endif

@@ Of course, we need to define a random number stream for the randomized
choice$\ldots$

@@<Declare the data structures for Greedy@@>=
#if defined(DO_RANDOM)
prng_t *random_stream=NULL;
#endif

@@ $\ldots$and initialize it too.
@@^Commodore PET@@>
@@^Commodore 64@@>
@@^C64@@>
@@^6502@@>
@@^6510@@>

@@<Build the data structures for Greedy.@@>=
#if defined(DO_RANDOM)
random_stream = prng_new(PRNG_DEFAULT,random_seed^(6502*6510));
#endif

@@ And when we're done, we should get rid of it too.
@@<Destroy the data structures for Greedy@@>=
#if defined(DO_RANDOM)
prng_free(random_stream);
#endif

@@
An edge is eligible if both ends are
not yet saturated (\ie, if both ends have degree less than two), and
if the two ends belong to different fragments (\ie, are not the two
tails of the
same fragment).
	
First we find the shortest remaining edge, $(x,y)$.  If it satisfies both 
of the above criteria, then we break out of the loop, and pass it on.
Otherwise, it is a stale link, \ie, edges added since its insertion
into the priority queue have made it obsolete.  In that case, we must
find a new valid link and insert it into the priority queue.

@@<Get one short valid edge or |NULL|@@>=
@@<Verbose: PQ before getting an edge@@>@@;
while (1) {
	e = pq_delete_min(pq_edge);
	if ( e==NULL ) break;
	@@<Verbose: extracted |e|@@>@@;
	x = e->this_end; y=e->other_end;
	if ( adj[x].count == 2 ) { /* |x| is saturated */
		continue; 
	}
	if ( adj[y].count <2 && y!=tail[x] )  break;	/* We found a valid edge */
	@@<Insert a new link for |x|, using |e|@@>@@;
}
@@<Verbose: got a valid edge@@>@@;

@@ Once we decide upon an edge, we have some accounting to do.
First we set |x| and |y|.
Then we must add the length of the edge upon which we decided.

@@<Find a short valid edge $e=(x,y)$@@>=
errorif(e==NULL,"Exhausted the priority queue of links.");
@@<Verbose: accept this edge@@>@@;
len += e->len;
x=e->this_end;
y=e->other_end;


@@ This section finds a new candidate link for city |x|.  The new link
is written into the space to which |e| points.

In the non-Euclidean case, we might have many other neighbours of |x| 
already
in the queue.  So we only refresh if we have just used up the furthest
one that was in the queue.

In the Euclidean case, 
we must first hide the
city at the other end of its fragment.

@@<Insert a new link for |x|, using |e|@@>=
if ( E2_case ) {
	E2_hide(tail[x]);
	e->other_end = E2_nn(x);
	e->len = cost(x,e->other_end);
	E2_unhide(tail[x]);
	pq_insert(pq_edge,e);
} else {
	pool_free(nn_link_pool,e);
	if ( farthest_in_queue[x] == y ) {
		const int not_me = tail[x];
		@@<Verbose: need to refresh@@>@@;
		@@<Get fresh neighbours for |x|, excepting |not_me|@@>@@;
	}
}

@@ Once the tour is built in |adj|, we write it in simpler form to |tour|.
This is a simple matter of following indices around. 

@@<Build the Greedy tour@@>=
{@@+ int i=0, prev=-1, here=0;
	do {
		tour[i++]=here;
		if ( adj[here].city[0] != prev ) {@@+prev=here;@@+here=adj[here].city[0];@@+}
		else {@@+prev=here;@@+here=adj[here].city[1];@@+}
	} while( here != 0 );
	errorif(i!=n,"Not a tour.");
}


@@* Weighted perfect matchings.

Procedure |construct_matching| builds a weighted perfect matching.
It is analogous to procedure |construct|.

The first parameter is the number of vertices in the graph.

The second parameter is a buffer in which to write the matching.
After we're through, edge $(u,v)$ is in the matching if and only if
|mate[u]=v| and |mate[v]=u|.

The third parameter specifies which algorithm should be used to construct
the matching, as with procedure |construct|.

The fourth parameter can be specified on the command line
and may be used (or ignored) by the algorithm for whatever purpose
the algorithm desires.

The routine returns the weight of the constructed perfect matching.

@@<Subroutines@@>=
length_t 
construct_matching(int n, int *mate, int alg,long alg_param, const long
random_seed) 
{
	int i;
	length_t weight=0;
	errorif(n%2,
		"Perfect matchings need an even number of vertices; given %d", n);
	errorif(mate==NULL,
			"Tried to construct a matching before space is allocated");
	switch(alg) {
	@@<Matching cases@@>@@;
	default:
		errorif(1,"Unrecognized matching construction algorithm %d",alg);
	}
	return weight;
}

@@ 
@@<Exported subroutines@@>=
length_t construct_matching(int n, int *mate, int alg,long alg_param, const
long random_seed);


@@ The first algorithm is rather trivial.  It just pairs each vertex
with its positional neighbour.

Now, this might seem like an awful matching.  But the user has the
option of pre-sorting the points so that they are in the order
determined by Moore's space-filling curve (a variant of Hilbert's curve).
In that case, this canonical perfect matching  
might be an excellent starting point.

@@<Matching cases@@>=
case CONSTRUCT_CANONICAL:
	for(i=0;i<n;i+=2) {
		mate[i]=i+1;
		mate[i+1]=i;
		weight += cost(i,i+1);
	}
	break;

@@ The second option is a random perfect matching.
It's probably very bad, but at least it's fast.

We use |alg_param| as the random number generator's seed.
(We may want to use different pseudo random number generator algorithms
via the facilities provided by module \module{PRNG}, but I'm not too
worried about using \module{GB\_FLIP} in this case.)


@@<Matching cases@@>=
case CONSTRUCT_RANDOM:
	{ @@;
	int *unmated=new_arr_of(int,n), num_unmated=n, u, v, ui, vi;
	prng_t *random_stream = prng_new(PRNG_DEFAULT,alg_param);
	for (i=0;i<n;i++) unmated[i]=i;
	while ( num_unmated > 0 ) {
		ui = prng_unif_int(random_stream,num_unmated);
		u = unmated[ui]; unmated[ui]=unmated[--num_unmated];
		vi = prng_unif_int(random_stream,num_unmated);
		v = unmated[vi]; unmated[vi]=unmated[--num_unmated];
		mate[u]=v;
		mate[v]=u;
		weight+=cost(u,v);
	}
	prng_free(random_stream);
	free_mem(unmated);mem_deduct(sizeof(int)*n);
	}
	break;

@@ The last option is a perfect matching discovered by a greedy algorithm.
Beginning with an empty matching, we repeatedly add the next shortest
edge that is not incident upon any vertex already in the matching.

It has much the same structure as the greedy algorithm for constructing
a greedy tour.  This time around, we don't need to use the |adj| or |tail|
arrays.  

@@<Matching cases@@>=
case CONSTRUCT_GREEDY:
	{
	const int E2_case = E2_supports(tsp_instance);
#define DO_MATCHING
	@@<Declare the data structures for Greedy@@>@@;
	@@<Build the data structures for Greedy.@@>@@;
	@@<Build the Greedy matching@@>@@;
	@@<Destroy the data structures for Greedy@@>@@;
#undef DO_MATCHING
	}
	break;

@@ Here's the randomized case.  It's all the same, but with different
code as marked by |DO_RANDOM|.

@@<Matching cases@@>=
case CONSTRUCT_GREEDY_RANDOM:
	{
	const int E2_case = E2_supports(tsp_instance);
#define DO_RANDOM
#define DO_MATCHING
	@@<Declare the data structures for Greedy@@>@@;
	@@<Build the data structures for Greedy.@@>@@;
	@@<Build the Greedy matching@@>@@;
	@@<Destroy the data structures for Greedy@@>@@;
#undef DO_MATCHING
#undef DO_RANDOM
	}
	break;

@@ We use the |mate| array itself to record when a city has been saturated.
We initialize each entry to a negative number, which can never be a city 
label.

@@<Other initializations for city |i|@@>=
#if defined(DO_MATCHING)
mate[i]=-1;
#endif


@@ With the priority queue and the nearest neighbour links all set up, 
the greedy algorithm looks like the following.


@@<Build the Greedy matching@@>=
{
int num_remaining =n/2;@@;
while( num_remaining > 0 ) {@@/
	int u, v;
	pq_edge_t *next_edge;
	@@<Find a short |next_edge| for matching@@>@@;
	errorif ( next_edge==NULL,
		"Priority queue exhausted while we expect 2*%d more", num_remaining);
	u = next_edge->this_end;
	v = next_edge->other_end;
	if( mate[u]>=0 || mate[v]>=0 ) {
		errorif(mate[u]>=0,"Mate[%d] = %d is not -1",u,mate[u]);
		errorif(mate[v]>=0,"Mate[%d] = %d is not -1",v,mate[v]);
	}
	mate[u]=v; /* Insert |(u,v)| into the matching. */
	mate[v]=u;
	if (E2_case) E2_hide(u); else { const int c=u; @@<Mark city |c| as
saturated@@>@@;}
	if (E2_case) E2_hide(v); else { const int c=v; @@<Mark city |c| as
saturated@@>@@;}
	weight += cost(u,v);
	num_remaining--;
}
}
if (E2_case) E2_unhide_all();


@@ Just as with tours, we have the deterministic case  and the randomized
case.

The randomized case looks the same as in the TSP context.  
It's a shame (and my fault) that we couldn't just use the same code sequence.  
Alas.

@@<Find a short |next_edge| for matching@@>=
#if defined(DO_RANDOM)
{
	pq_edge_t *candidate[2];
	@@<Get one short |next_edge| for matching, or |NULL|@@>@@; 
	candidate[0]=next_edge;
	@@<Get one short |next_edge| for matching, or |NULL|@@>@@; 
	if ( next_edge && next_edge->this_end == candidate[0]->other_end @@| 
			&& next_edge->other_end == candidate[0]->this_end ) {
		pq_edge_t *push_back = next_edge;
		@@<Get one short |next_edge| for matching, or |NULL|@@>@@;  /* Get another. */
		pq_insert(pq_edge,push_back);  /* But put the second one back. */
	}
	candidate[1]=next_edge;
	next_edge = candidate[0];
	if ( candidate[1] ) { /* There is a choice to be made. */
		const int chosen = prng_unif_int(random_stream,3) == 0; 
			/* Prob|[chosen==1]== 1/3| */
		next_edge=candidate[chosen];
		pq_insert(pq_edge,candidate[1-chosen]); /* Put the other one back. */
	}
}
#endif

@@ But the deterministic case is simple.
@@<Find a short |next_edge| for matching@@>=
#if !defined(DO_RANDOM)
@@<Get one short |next_edge| for matching, or |NULL|@@>@@;
#endif


@@ 
In the perfect matching case, a node is saturated if it has even a single
edge incident upon it.  We maintain the invariant that a city is hidden
in the $k$-d tree if and only if it is saturated.

If this end of the edge |u| is already saturated, then don't bother
finding another possible mate for it.  If |u| is not saturated,
then we have a hope of adding the current edge, but only if |v| (the
other end of this edge) is not already saturated.

You'll have to agree this is much simpler than the greedy algorithm
for tours.  (You don't?  Oh. OK.)

Note that, as opposed to the tour case, we might add extra elements
to the priority queue while getting just one edge.  But since the
(partial) matching is not changing, the first edge we withdrew
(in the random case) is still valid.  So this is ok.

@@<Get one short |next_edge| for matching, or |NULL|@@>=
while (1) {
	next_edge = pq_delete_min(pq_edge);
	if ( next_edge == NULL ) break; /* Queue is empty. */
	u = next_edge->this_end;
	v = next_edge->other_end;
	if ( mate[u] >= 0 ) continue; /* Node |u| is saturated. */
	if ( mate[v] >= 0 ) {
		@@<Pick a new potential mate for |u|@@>;
	} else break;		/* We got a valid edge. */
}
	

@@ Recall that all mated cities are hidden.  Also, the $k$-tree never
says that a city is its own nearest neighbour.

In the non-Euclidean case, we insert new possible mates only if we've
exhausted
the list of neighbours for |u| already in the queue.