nn.w,v

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

asked to add: an unqualified (pure) near neighbour, a near neighbour
in a particular quadrant, or a city nearby in the Delauney triangulation.
The list of neighbours for each city $i$ should be:
the |num_pure| closest neighbours of $i$; 
the |num_quadrant| closest neighbours of $i$ in each of the four
surrounding quadrants; 
all cities that are at most |num_delauney| steps from $i$ 
in the Delauney triangulation of the instance.  
Of course, the same city may qualify under more than one of these criteria
for a given seed.  Such cities are included only once for that seed city,
\ie, duplicates are removed.

@@d max(A,B) ((A)>(B) ? (A) : (B))

@@<Subroutines@@>=
void
nn_build(int num_pure, int num_quadrant, int num_delauney) 
{
	@@<|nn_build| variables@@>@@;
	int i;
	n = tsp_instance->n;
	@@<Check build parameters@@>@@;
	
	@@<Set up the array data structures@@>@@;
	nn_max_bound = 0;
	for (i=0;i<n;i++) {
		@@<|nn_build| iteration variables@@>@@;
		@@<Verbose: about to start building list for |i|@@>@@;
		@@<Build a pure work list for city |i|@@>@@;
		@@<Build a quadrant work list for city |i|@@>@@;
		@@<Build a Delauney work list for city |i|@@>@@;
		@@<Compress the work lists and copy them into |list|@@>@@;
		nn_max_bound = max(nn_max_bound, begin[i+1]-begin[i]);
	}
	@@<Clean up the temporary data structures@@>@@;
	@@<Verbose: report number of vertices found@@>@@;
}

@@ We export the building subroutine.
@@<Exported subroutines@@>=
void nn_build(int num_pure, int num_quadrant, int num_delauney);

@@ We check for bad build parameters.  This helps us to spot bugs and 
to avoid misleading the user, \eg, if they specified 50 nearest neighbours
but only get 20 because there are only 21 cities.

@@<Check build parameters@@>=
@@<Verbose: show build parameters@@>@@;
errorif(num_pure<0,
	"Need positive number of nearest neighbours; %d specified",num_pure);
errorif(num_quadrant<0,
	"Need positive number of quadrant neighbours; %d specified",num_quadrant);
errorif(num_delauney<0,
	"Need positive Delauney depth; %d specified",num_delauney);
errorif(num_pure<=0 && num_quadrant<=0 && num_delauney<=0,
	"Must specify some candidates");
errorif(num_pure>= n,
	"%d nearest neighbours specified, but there are only %d cities",num_pure,n);

@@ We need to import the instance |tsp_instance| and its type.
@@<Early module headers@@>=
#include "read.h"
#include "lk.h"


@@ We try to reduce the number of reallocations of 
|list| by guessing the average length of each city's neighbour list.
The following expressions are themselves just guesses.


@@<Guess the average number of neighbours per city@@>=
if ( num_pure )
	guess_avg = num_pure + num_quadrant + num_delauney;
else if ( num_quadrant )
	guess_avg = 4*num_quadrant + num_delauney;
else
	guess_avg = 3*num_delauney*num_delauney;

@@ Once we have a guess at the average length of a list, we initially
allocate |n| times that number of entries for |list|.

@@<Set up the array data structures@@>=
{int guess_avg;
@@<Guess the average number of neighbours per city@@>@@;
list_size = guess_avg * n;
list = new_arr_of(int,list_size);
}

@@
@@<Clean up the array data structures@@>=
if ( list ) { free_mem(begin);mem_deduct(sizeof(int)*list_size); }

@@ Variable |list_size| always contains the number of elements allocated
for |list|.

@@<Module variables@@>=
static int list_size;

@@ We will need work space for each of the three passes through the data.
In the worst (and exceedingly rare) case, 
each city will appear three times in the work list.

@@<Set up the array data structures@@>=
work = new_arr_of(kd_bin_t,3*n);

@@ We don't need the |work| array after we build the lists, so we clean it
up early.  
@@<Clean up the temporary data structures@@>=
free_mem(work);mem_deduct(sizeof(kd_bin_t)*3*n);

@@ But we must declare the work space.
@@<|nn_build| variables@@>=
kd_bin_t *work;

@@ The work list will be built by just appending.  Specifically, the 
working copy of the candidate list for a city is stored in  
|work[0]| through |work[work_next-1]|.  

We sort to find duplicates. After sorting, we don't ever use the
lengths again.

Module \module{LK} sets the |sort| function pointer to a suitable sorting
procedure, \eg, the \CEE/ library |qsort| function.

@@<Compress the work lists and copy them into |list|@@>=
errorif( work_next < 1, "Must have nonempty candidate list"); 
{
int r, w, last_city;
sort(work,(size_t)work_next,sizeof(kd_bin_t),kd_bin_cmp_increasing);
for ( r=w=0, last_city=kd_bin_city(&work[r])-1; r<work_next; r++ ) {
	if ( kd_bin_city(&work[r]) != last_city )
		last_city = kd_bin_city(&work[w++]) = kd_bin_city(&work[r]);
}
@@<Make sure |list| has space for |w| more elements@@>@@;
for (r=0;r<w;r++) list[begin[i]+r] = kd_bin_city(&work[r]);
begin[i+1] = begin[i]+w;
}

@@ We must declare |work_next|.
@@<|nn_build| iteration variables@@>=
int work_next=0;


@@ We use repeated doubling to expand the |list| array to achieve a constant
amortized allocation time per cell.  

It's ok if we end up using horribly 
less than the total allocated space.  It would only end up wasting address
space, which we of course know is infinite.  
There is no performance penalty because we never access the memory we don't
use.  Duh.  Unused memory is never paged in.

@@<Make sure |list| has space for |w| more elements@@>=
if ( begin[i]+w > list_size ) {
	int new_size = list_size;
	while ( begin[i]+w > new_size ) new_size*=2;
	@@<Verbose: resize |list|@@>@@;
	list = mem_realloc(list,sizeof(int)*new_size);
		mem_deduct(sizeof(int)*list_size);
	list_size = new_size;
}

@@*1 Pure lists.
Whenever possible, we use the $k$-d tree to speed things up.

Otherwise, we use a quick and dirty method to calculate
the nearest neighbour list.  For each city, compute the distances to all
the cities, sort, and then pick off the first few entries.
In the future we may want to do smart partitioning, like a truncated
Quicksort, perhaps with randomized pivot sampling.


@@<Build a pure work list for city |i|@@>=
if ( num_pure ) {
	int start_work = work_next;
	if (E2_supports(tsp_instance)) {
		@@<Use the $k$-d tree to find nearest neighbours@@>@@;
	} else {
		@@<Compute all the distances except from $i$@@>@@;
		@@<Move the |num_pure| closest cities to the front@@>@@;
	}
}

@@ The |E2| routines for 2-d trees are in the \module{KDTREE} module.
@@<Early module headers@@>=
#include "kdtree.h"

@@ Let's do the fast way first,
using the $k$-d search structure implemented in \module{KDTREE}.
Here we take advantage of the bulk search routine.
only use the nearest neighbour search and the hiding and unhiding 
routines.  Given those calls, building the list is rather easy.

@@<Use the $k$-d tree to find nearest neighbours@@>=
work_next += E2_nn_bulk(i,num_pure,work+work_next);

@@ Now we're ready to do the brute force method. 
Computing all the distances is easy.

@@<Compute all the distances except from $i$@@>=
{ int j;
	for (j=0;j<i;j++) {
		kd_bin_city(&work[start_work+j]) = j;
		kd_bin_monotonic_len(&work[start_work+j]) = cost(i,j);
	}
	for (j=i+1;j<n;j++) {
		kd_bin_city(&work[start_work+j-1]) = j;
		kd_bin_monotonic_len(&work[start_work+j-1]) = cost(i,j);
	}
}


@@  We use the handy-dandy partial sorting selection function |select_range| from
the \module{DSORT} module to
move the closest |num_pure| cities to the front of the array, and sort them
too.

Recall |num_pure <= n-1|, so we won't run off the end of what we've already
written.

@@<Move the |num_pure| closest cities to the front@@>=
select_range((void *)work,(size_t)n-1,sizeof(kd_bin_t),kd_bin_cmp_increasing,
	0,num_pure,0/*Don't sort them */);
work_next += num_pure;

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


@@*1 Quadrant-balanced lists.
Quadrant balanced lists only make sense when the instance lives in a
coordinate system.  Since the only coordinate system we support is
2-dimensional real space, we just use the 2-d trees.  In fact, quadrant
lists may be used precisely when 2-d trees are supported.

@@<Check build parameters@@>=
errorif(num_quadrant>0 && !E2_supports(tsp_instance),
	"Quadrant lists supported only when 2-d trees supported",num_pure);

@@ Quadrant lists are interesting for two reasons.

First, building them efficiently requires special support in the 2-d tree.
That is, we need to be able to find nearest neighbours while masking
out entire quadrants of space surrounding the seed city.

Second, although the user may ask for a certain number of cities from
each quadrant, there may not be that many available there.  In fact, this
is certain to be the case for cities at the ``corners'' of the instance.
The only consequence of running short is that the caller gets back a
shorter list.

We build the quadrant city lists last, and can therefore take advantage of
the fact that the nearest neighbour search has already saturated certain
quadrants.  If any of the non-trivial quadrants already have cities,
then quadrant 0 need not be searched.

@@<Build a quadrant work list for city |i|@@>=
if ( num_quadrant ) {
	int quadrant;
	int q_count[5] = {0,0,0,0,0};
	@@<Count the number of cities in each quadrant so far@@>@@;
	if ( 0 == q_count[1]+q_count[2]+q_count[3]+q_count[4] ) {
		@@<Find all cities in quadrant 0@@>@@;
	}
	for ( quadrant = 1; quadrant <= 4; quadrant++ ) {
		if ( q_count[quadrant] < num_quadrant ) {
			@@<Find up to |num_quadrant| neighbours in quadrant |quadrant|@@>@@;
		}
	}
}

@@ 
I borrow from and slightly extend ordinary mathematical 
terminology: if city |i| lies
at $(x,y)$, then 
quadrant 0 is point $(x,y)$;
quadrant 1 is all points $(x',y')$ with $x'>x$ and $y'\ge y$;
quadrant 2 is all points $(x',y')$ with $x'\le x$ and $y'>y$;
quadrant 3 is all points $(x',y')$ with $x'<x$ and $y'\le y$;
quadrant 4 is all points $(x',y')$ with $x'\ge x$ and $y'<y$.

I've been careful to define the boudaries of quadrants 1 through 4
to be balanced: those quadrants are pairwise congruent to each other.

@@<Count the number of cities in each quadrant so far@@>=
{ int j;
coord_2d *coord = tsp_instance->coord;
for ( j=0; j < work_next; j++ ) {
	const double diff_x = coord[i].x[0] - coord[j].x[0];
	const double diff_y = coord[i].x[1] - coord[j].x[1];
	q_count[E2_quadrant(diff_x,diff_y)]++;
}
}


@@ Since all five quadrants are disjoint, we don't need to worry about
duplicates.  So we write directly to the |work| array.

The quadrant parameter passed to the |E2_nn_quadrant| 
is a bit mask of all the quadrants in which to search.

@@<Find up to |num_quadrant| neighbours in quadrant |quadrant|@@>=
work_next += E2_nn_quadrant_bulk(i,num_quadrant,work+work_next,1<<quadrant);

@@ 
Now, we would hope to
not have cities with duplicated coordinates.  But they might be there.
Let's get all of them, no matter what |num_quadrant| is.

We cheat a little by temporarily redefining the |num_quadrant| variable.

@@<Find all cities in quadrant 0@@>=
{ const int num_quadrant = n-1, quadrant=0;
@@<Find up to |num_quadrant| neighbours in quadrant |quadrant|@@>@@;
}


@@*2 Delauney lists.
We don't support Delauney lists.  Alas.  For further pointers, see
either Knuth's {\sl The Stanford Graphbase} or 
Preparata and Shamos' {\sl Computational Geometry}.

@@
@@<Check build parameters@@>=
errorif(num_delauney,"Delauney neighbours not supported (yet)");

@@
@@<Build a Delauney work list for city |i|@@>=

@@*Verbose.

@@<Verbose: resize |list|@@>=
if ( verbose >= 750 ) {
	printf("nn: Resize list from %d elements to %d elements; begin[i]=%d, w=%d\n", 
		list_size, new_size, begin[i], w);
	fflush(stdout);
}

@@
@@<Verbose: show build parameters@@>=
if ( verbose >= 500 ) {
printf("nn: build nn %d nq %d del %d\n", num_pure, num_quadrant, num_delauney);
	fflush(stdout);
}

@@
@@<Verbose: report number of vertices found@@>=
if ( verbose >= 75 ) {
printf("nn: build nn %d nq %d del %d got %d total neighbours\n",
	num_pure, num_quadrant, num_delauney, begin[n] );
	fflush(stdout);
}

@@
@@<Verbose: about to start building list for |i|@@>=
if ( verbose >= 1250 ) {
	printf("nn: about to build for %d; work_next=%d list_size=%d\n", 
		i, work_next, list_size);
	fflush(stdout);
}

@@*Index.

@


1.133
log
@Must always look for quadrant 0 cities...
@
text
@d54 3
d222 1
a222 1
const char *nn_rcs_id="$Id: nn.w,v 1.132 1998/10/09 21:34:52 neto Exp neto $";
d596 1
a596 1
We build the quadrant city last, and can therefore take advantage of
d606 3
a608 1
	@@<Find all cities in quadrant 0@@>@@;
@


1.132
log
@Report number of vertices found.
@
text
@d54 3
d219 1
a219 1
const char *nn_rcs_id="$Id: nn.w,v 1.131 1998/10/09 19:21:19 neto Exp $";
d603 1
a603 3
	if ( 0 == q_count[1]+q_count[2]+q_count[3]+q_count[4] ) {
		@@<Find all cities in quadrant 0@@>@@;
	}
@


1.131
log
@This works with bulk load
@
text
@d54 3
d216 1
a216 1
const char *nn_rcs_id="$Id: nn.w,v 1.130 1998/10/09 16:47:34 neto Exp neto $";
d373 1
d681 8
@


1.130
log
@Added changes need for bulk neighbour searches.
@
text
@d54 3
d213 1
a213 1
const char *nn_rcs_id="$Id: nn.w,v 1.129 1998/07/16 21:58:55 neto Exp neto $";
a599 1
		printf("%d qc[%d]=%d\n",i,quadrant,q_count[quadrant]);
a601 2
		} else {
			printf("%d skip quadrant %d\n",i,quadrant);
@


1.129
log
@Added the LGPL notice in each file.
@
text
@d54 3
a208 1
@@<Module subroutines@@>@@;
d210 1
a210 1
const char *nn_rcs_id="$Id: nn.w,v 1.128 1998/04/16 15:40:03 neto Exp neto $";
d434 1
a434 1
work = new_arr_of(nn_entry_t,3*n);
d439 1
a439 1
free_mem(work);mem_deduct(sizeof(nn_entry_t)*3*n);
d441 1
a441 1
@@ But we must declare the variable.
d443 1
a443 1
nn_entry_t *work;
d459 4
a462 4
sort(work,(size_t)work_next,sizeof(nn_entry_t),cmp_entry_compress);
for ( r=w=0, last_city=work[r].city-1; r<work_next; r++ ) {
	if ( work[r].city != last_city )
		last_city = work[w++].city = work[r].city;
d465 1
a465 1
for (r=0;r<w;r++) list[begin[i]+r] = work[r].city;
a473 16
@@ We do have to provide a comparison function. 
We promise to the outside world to put nearer
cities first on lists, so that's the first criterion.
Second, we want duplicates to be brought together, so that's the second
criterion.

@@<Module subroutines@@>=
static int 
cmp_entry_compress(const void *a, const void *b) {
	const length_t ad = ((const nn_entry_t *)a)->len;
	const length_t bd = ((const nn_entry_t *)b)->len;
	if ( ad < bd ) return -1;
	else if ( ad > bd ) return 1;
	else return ((const nn_entry_t *)a)->city - ((const nn_entry_t *)b)->city;
}

d489 1
d520 2
a521 1
Here we only use the nearest neighbour search and the hiding and unhiding 
d525 1
a525 12
{
	int j, city;
	for (j=0;j<num_pure;j++) {
		city = E2_nn(i);
		if ( city == -1 ) break; /* No such neighbour. */
		work[work_next].city = city;
		work[work_next].len = cost(i,city);
		work_next++;
		E2_hide(city);
	}
	for (j=start_work;j<work_next;j++) E2_unhide(work[j].city);
}
d533 2
a534 2
		work[start_work+j].city = j;
		work[start_work+j].len = cost(i,j);
d537 2
a538 2
		work[start_work+j-1].city = j;
		work[start_work+j-1].len = cost(i,j);
d552 1
a552 1
select_range((void *)work,(size_t)n-1,sizeof(nn_entry_t),cmp_entry_compress,
d583 5
d591 5
a595 1
	@@<Find all cities in quadrant 0@@>@@;
d597 6
a602 1
		@@<Find up to |num_quadrant| neighbours in quadrant |quadrant|@@>@@;
d619 12
a630 1
Since all five quadrants are disjoint, we don't need to worry about
d637 1
a637 11
{ int j, start_work=work_next;
	for (j=0;j<num_quadrant;j++) {
		const int city = E2_nn_quadrant(i,1<<quadrant);
		if ( city == -1 ) break; /* No such neighbour. */
		work[work_next].city = city;
		work[work_next].len = cost(i,city);
		work_next++;
		E2_hide(city);
	}
	for (j=start_work;j<work_next;j++) E2_unhide(work[j].city);
}
@


1.128
log
@Need to include dsort.h to get the select range prototype.
@
text
@d1 47
a47 1
@@i copyrt.w
d54 3
d208 1
a208 1
const char *nn_rcs_id="$Id: nn.w,v 1.127 1998/04/11 14:34:09 neto Exp neto $";
@


1.127
log
@Duh.  select is a Unix system call.  Renamed my select to select_range