read.w,v

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

@@<Early type definitions@@>@@;
typedef struct {
	char *name;
	char *comment;
	int n;
	int input_n;
	double scale;
	double dsj_random_factor;
	double dsj_random_param;
	@@<Other |tsp_instance_t| fields@@>@@;
} tsp_instance_t;

@@ We specify defaults for these fields so that if there is a problem
reading the file, we can detect it later on.

@@<Set the instance defaults@@>=
p->name = NULL;
p->comment = NULL;
p->n = 0;
p->scale = 1e6;

@@ The different parts of the input file are introduced by various keywords.
So we read a keyword, and then do the appropriate thing.

@@d matchword(STR) (0==strncmp(keyword,STR,strlen(STR)))
	
@@<Read in the instance@@>=
{
#define MAX_KEYWORD_LEN (25)
#define MAX_LINE_LEN (200)
char keyword[MAX_KEYWORD_LEN], rest_of_line[MAX_LINE_LEN];
int more_input = 1, lineno=0;
while( more_input) {
	char *colon; int r;
	keyword[0]=0;
	r = fscanf(in," %s ", keyword); 
	if ( r==EOF )@@+ {@@+more_input=0;@@+break;@@+}
	errorif( r != 1, "%d: Couldn't read the word following \"%s\". (r==%d)",
			 lineno+1,keyword,r);
	if ( NULL!=(colon = strchr(keyword,':')) ) *colon = '\0';
	if ( matchword("EOF") ) {
		more_input = 0;
	}@@+ else if ( matchword("NAME") ) {
		@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
		p->name = dup_string(rest_of_line);@@;
	}@@+ else if ( matchword("COMMENT") ) {
		@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
		p->comment = dup_string(rest_of_line);@@;
	}@@+ else if ( matchword("TYPE") ) {
		@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
		errorif(strcmp(rest_of_line,"TSP"),"Can't read TSPLIB files of type %s.",rest_of_line);@@;
	}@@+ else if ( matchword("DIMENSION") ) {
		@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
		p->input_n = p->n = atoi(rest_of_line);
	}@@+ else if ( matchword("DISPLAY_DATA_TYPE") ) {
		@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
		/* Do nothing. */
	}@@+ else if ( matchword("EDGE_WEIGHT_TYPE") ) {
		@@<Parse the edge weight type@@>@@;
	}@@+ else if ( matchword("EDGE_WEIGHT_FORMAT") ) {
		@@<Parse the edge weight format@@>@@;
	}@@+ else if ( matchword("EDGE_WEIGHT_SECTION") ) {
		@@<Parse the edge weights@@>@@;
	}@@+ else if ( matchword("NODE_COORD_SECTION") ) {
		@@<Parse the node coordinates@@>@@;
	}@@+ else if ( matchword("DISPLAY_DATA_SECTION") ) {
		@@<Parse the display data@@>@@;
	}@@+ else if ( matchword("SEED") ) {
		@@<Set the seed and possibly generate edge lengths@@>@@;
	}@@+ else if ( matchword("SCALE") ) {
		@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
		p->scale = atof(rest_of_line);
		p->dsj_random_factor = p->scale/2147483648.0;
	}@@+ else {
		@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
		errorif(1,"%d: Don't know what the keyword %s is!",lineno,keyword);
	}
}
@@<Set the |cost| function according to the current context@@>@@;
}

@@ We need to include the interface definition for the string handling
functions.
@@<Other includes@@>=
#include <string.h>

@@
We make sure that the colon is actually there, and that the line
lengths aren't violated.

Also, erase trailing whitespace.

@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>=
{ 
int l; char *rcp;
if ( colon == NULL ) {
	int r=fscanf(in," : ");
	errorif( r != 0, "%d: Missed the colon.", lineno);
}
rcp = fgets(rest_of_line,MAX_LINE_LEN,in);
errorif( NULL == rcp,
	@@|"%d: Couldn't read after the colon; truncated file?", lineno);@@;
l = strlen(rest_of_line)-1;@@;
if ( l>=0 && rest_of_line[l] == '\n' ) { rest_of_line[l--]='\0'; lineno++;}
while ( l >= 0 && isspace(rest_of_line[l]) ) rest_of_line[l--] = '\0';@@;
if ( feof(in) )@@+more_input = 0;
}

@@ We need to include the character classifiction macros.
@@<Other includes@@>=
#include <ctype.h>

@@ The |RANDOM_EDGES| edge weight type is an extension to the TSPLIB format.  
It specifies that the edge lengths should be random and uniformly
distributed between 1 and 1000.  The keyword |SEED| must be used
to specify a seed value and to generate the edge lengths.

@@<Parse the edge weight type@@>=
@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
if ( 0==strcmp(rest_of_line,"EUC_2D") ) {
	p->edge_weight_type = EUC_2D;
}@@+else if ( 0==strcmp(rest_of_line,"CEIL_2D") ) {
	p->edge_weight_type = CEIL_2D;
}@@+else if ( 0==strcmp(rest_of_line,"GEO") ) {
	p->edge_weight_type = GEO;
}@@+else if ( 0==strcmp(rest_of_line,"ATT") ) {
	p->edge_weight_type = ATT;
}@@+else if ( 0==strcmp(rest_of_line,"DSJ_RANDOM") ) {
	p->edge_weight_type = DSJ_RANDOM;
}@@+else if ( 0==strcmp(rest_of_line,"EXPLICIT") ) {
	p->edge_weight_type = EXPLICIT;
}@@+else if ( 0==strcmp(rest_of_line,"RANDOM_EDGES") ) {
	p->edge_weight_type = RANDOM_EDGES;
}@@+else {
	errorif(1,"%d: Apology: Unknown edge weight type \"%s\".",lineno,rest_of_line);
}

@@ @@<Early type definitions@@>=
typedef enum 
{ NO_EDGE_TYPE, EUC_2D, CEIL_2D, GEO, ATT, DSJ_RANDOM, EXPLICIT, RANDOM_EDGES } edge_weight_type_t;

@@ 
@@<Module variables@@>=
static const char edge_weight_type_name[] =
{ "NO_EDGE_TYPE", "EUC_2D", "CEIL_2D", "GEO", "ATT", "DSJ_RANDOM", "EXPLICIT", "RANDOM_EDGES"};

@@ @@<Other |tsp_instance_t| fields@@>=
edge_weight_type_t edge_weight_type;

@@ @@<Set the instance defaults@@>=
p->edge_weight_type = NO_EDGE_TYPE;

@@ @@<Other includes@@>=
#include "gb_flip.h"




@@ @@<Parse the edge weight format@@>=
@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
if ( 0==strcmp(rest_of_line,"LOWER_DIAG_ROW") ) {
	p->edge_weight_format = LOWER_DIAG_ROW;
}@@+else if ( 0==strcmp(rest_of_line,"FULL_MATRIX") ) {
	p->edge_weight_format = FULL_MATRIX;
}@@+else if ( 0==strcmp(rest_of_line,"UPPER_ROW") ) {
	p->edge_weight_format = UPPER_ROW;
}@@+else {
	errorif(1,"%d: Unknown edge weight format \"%s\".",lineno,rest_of_line);
}

@@ @@<Early type definitions@@>=
typedef enum { NO_EDGE_FORMAT, LOWER_DIAG_ROW, FULL_MATRIX, UPPER_ROW } 
	edge_weight_format_t;

@@ 
@@<Module variables@@>=
static const char edge_weight_format_name[] =
{ "NO_EDGE_FORMAT", "LOWER_DIAG_ROW", "FULL_MATRIX", "UPPER_ROW"};

@@ @@<Other |tsp_instance_t| fields@@>=
edge_weight_format_t edge_weight_format;

@@ @@<Set the instance defaults@@>=
p->edge_weight_format = NO_EDGE_FORMAT;


@@*Explicit edge weights.
We parse the edge weights according to the format they are given.  The
file format guarantees that we know this information before this section
arrives.

@@<Parse the edge weights@@>=
switch( p->edge_weight_format ) {
long int long_dummy;
case LOWER_DIAG_ROW :
	{ 
		int row, col;
		p->edge_weights = new_arr_of(length_t*,p->input_n);
			/* Allocate the space for the edge weights. */
		for ( row = 0 ; row < p->input_n ; row++ ) {
			p->edge_weights[row] = new_arr_of(length_t,p->input_n);
		}
		for ( row = 0 ; row < p->input_n ; row++ ) { /* Now actually read the weights. */
			for ( col = 0 ; col <= row ; col++ ) {
				int r= fscanf(in," %ld ", &long_dummy);
				errorif( 1 != r,
					"Couldn't convert an edge weight: %d to %d.",row+1,col+1);
				p->edge_weights[col][row] = p->edge_weights[row][col] = long_dummy;
			}
		}
	}
	break;
case FULL_MATRIX :
	{
		int row, col;
		p->edge_weights = new_arr_of(length_t*,p->input_n);
			/* Allocate the space for the edge weights. */
		for ( row = 0 ; row < p->input_n ; row++ ) {
			p->edge_weights[row] = new_arr_of(length_t,p->input_n);
		}
		for ( row = 0 ; row < p->input_n ; row++ ) { /* Now actually read the weights. */
			for ( col = 0 ; col < p->input_n ; col++ ) {
				int r= fscanf(in," %ld ", &long_dummy);
				errorif( 1 != r,
					"Couldn't convert an edge weight: %d to %d.",row+1,col+1);
				p->edge_weights[row][col] = long_dummy;
			}
		}
		for ( row = 0; row < p->input_n ; row++ ) {
			for ( col = 0; col < row ; col++ ) {
				errorif(p->edge_weights[row][col] != p->edge_weights[col][row],
					"Asymmetric FULL_MATRIX:  (%d,%d) does not match (%d,%d)",
					row,col,col,row);
			}
		}
	}
	break;
case UPPER_ROW:
	{ 
		int row, col;
		p->edge_weights = new_arr_of(length_t*,p->input_n);
			/* Allocate the space for the edge weights. */
		for ( row = 0 ; row < p->input_n ; row++ ) {
			p->edge_weights[row] = new_arr_of(length_t,p->input_n);
			p->edge_weights[row][row] = 0;
		}
		for ( row = 0 ; row < p->input_n ; row++ ) { 
			/* Now actually read the weights. */
			for ( col = row+1 ; col < p->input_n ; col++ ) {
				int r= fscanf(in," %ld ", &long_dummy);
				errorif( 1 != r,
					"Couldn't convert an edge weight: %d to %d.",row+1,col+1);
				p->edge_weights[col][row] = p->edge_weights[row][col] = long_dummy;
			}
		}
	}
	break;
default: break;
}

@@ We need the |length_t| type.
@@<Other includes@@>=
#include "length.h"

@@ @@<Other |tsp_instance_t| fields@@>=
length_t **edge_weights;

@@ @@<Set the instance defaults@@>=
p->edge_weights = NULL;


@@*Implicit edge weights.
In this case, the weights of the edges are not given explicitly in
the file, and we don't store them either.  Instead, each weight can
be computed from other data given in the file, usually coordinates
of the vertices of the problem.  

Currently the only implicit edge weight scheme we support is the
two-dimensional Euclidean distance.  
Unsupported schemes include
the geographical distances.

Beware of off-by-one errors in input and output:
the input format counts from 1; we count from 0.

For the purposes of graphical output, we also keep track of the bounding
box of the cities.

@@d min(X,Y) ((X)>(Y)?(Y):(X))
@@d max(X,Y) ((X)<(Y)?(Y):(X))
@@d check_num_read(R,K) 
	errorif((R)!=(K),"Error: expected %d arguments, got %d",K,R)

@@<Parse the node coordinates@@>=
{
	int i,j, r;
	p->coord = new_arr_of(coord_2d,p->input_n);	/* Allocate the space for the coordinates. */
	for ( i=0; i<p->input_n; i++ ) {
		r = fscanf(in," %d ",&j);  
		check_num_read(r,1);
		r = fscanf(in," %lf %lf ",&p->coord[j-1].x[0],&p->coord[j-1].x[1]);
		check_num_read(r,2);
		p->xmin = min(p->xmin,p->coord[j-1].x[0]);
		p->ymin = min(p->ymin,p->coord[j-1].x[1]);
		p->xmax = max(p->xmax,p->coord[j-1].x[0]);
		p->ymax = max(p->ymax,p->coord[j-1].x[1]);
		lineno++;
	}
}

@@ Index |0| refers to the $x$ coordinate,
and index |1| refers to the $y$ coordinate.

I use an array because code that indexes can be made shorter than
code that uses different names.  Also, it generalizes more easily
to higher dimensions.  As a bonus, a good compiler ought to emit identical
anyway.

@@<Early type definitions@@>=
typedef struct { double x[2]; } coord_2d;

@@ @@<Other |tsp_instance_t| fields@@>=
coord_2d *coord;
double xmin, xmax, ymin, ymax;

@@ I've used the limits on integer values for the minimum and maximum
coordinate values.  It would be nice to use |double|-valued constants,
but I can't seem to find them under SunOS.  (Header {\tt math.h} has
|HUGE_VAL|, but that appears to be an infinity; it breaks on the KSR.)
In any case, these integer values should be good enough.

@@<Set the instance defaults@@>=
p->coord = NULL;
p->xmin = p->ymin = INT_MAX;
p->xmax = p->ymax = INT_MIN;

@@ @@<Other includes@@>=
#include <limits.h>

@@ Currently we don't do anything with the display coordinates.  
In the future, we may want to output PostScript.
@@<Parse the display data@@>=
{ 	
	int i;
	double dummy;
	for ( i=0; i<p->input_n; i++ ) {
		fscanf(in," %lf %lf %lf ",&dummy,&dummy,&dummy);
	}
}

@@ The random seed can be used in one of two ways.  Either we can 
generate all the edges in advance (in case |RANDOM_EDGES|) or
on the fly (in case |DSJ_RANDOM|).  The latter case saves us 
a whole whack of space ---quadratic in the number of vertices---
and therefore can reasonably used for instances with 
upwards of a hundred thousand
vertices.

@@<Set the seed and possibly generate edge lengths@@>=
@@<Get the rest of the line into |rest_of_line|, but skip the colon@@>@@;
p->seed = atol(rest_of_line);
switch(p->edge_weight_type) {
case DSJ_RANDOM: /* Don't do anything now. */
	p->dsj_random_param = 1+104*p->seed;
	break;
case RANDOM_EDGES:
	@@<Generate random edge lengths using the given seed@@>@@;
	break;
default: 
	errorif(1,"SEED directive used for edge type %s",
			edge_weight_type_name[p->edge_weight_type]);
}

@@ The seed is part of the instance structure.
@@<Other |tsp_instance_t| fields@@>=
long seed;

@@ @@<Set the instance defaults@@>=
p->seed = 1;
p->dsj_random_factor=1.0;
p->dsj_random_param=1;

@@
@@<Module variables@@>=
static double dsj_random_factor=1;
static int32 dsj_random_param=1;

@@ We use |short|s to store the randomly generated edge lengths because
they fit.  This scheme saves space and time because of cache effects.

@@<Generate random edge lengths using the given seed@@>=
gb_init_rand( p->seed );
{ int i,j;
p->short_edge_weights = new_arr_of(short*,p->input_n);
for (i=0;i<p->input_n;i++) {
	p->short_edge_weights[i]  = new_arr_of(short,p->input_n);
	p->short_edge_weights[i][i] = 0;
}
for (i=0;i<p->input_n;i++)
	for (j=0;j<i;j++) 
		p->short_edge_weights[i][j] = p->short_edge_weights[j][i] = 
			(short)(1+gb_unif_rand(1000L));
if ( p->input_n <=10 ) {
	printf("Cost matrix:\n");
	for (i=0;i<p->input_n;i++) {
		for (j=0;j<p->input_n;j++) printf(" %4d",p->short_edge_weights[i][j]);
		printf("\n");
	}
}
}

@@ @@<Other |tsp_instance_t| fields@@>=
short **short_edge_weights;

@@ @@<Set the instance defaults@@>=
p->short_edge_weights = NULL;

@@*Cost functions.
The cost function depends upon the type of the current problem.
The |cost| variable is a function pointer that refers to the appropriate
function, and it is always the case that 
|cost(i,j)| is the cost of the edge from |i| to |j| in the current problem.

We will often only be comparing the lengths of a pair of edges.
In these cases, we can substitute a function which is monotonic
in the real cost function.  
The |pseudo_cost| function should be monotonic in |cost|.
That is,
$${\mit cost}(i,j)\le {\mit cost}(k,l) \iff
{\mit pseudo\_cost}(i,j) \le {\mit pseudo\_cost}(k,l).$$
%$$ |cost(i,j)|\le|cost(k,l)|\iff
%|pseudo_cost(i,j)|\le|pseudo_cost(k,l)|.$$

@@<Global variables@@> =
length_t (*cost)(const int,const int);
length_t (*pseudo_cost)(const int,const int);

@@ We make this variable known to the outside world.
@@<Exported variables@@>=
extern length_t (*cost)(const int,const int);
extern length_t (*pseudo_cost)(const int,const int);


@@ When the edge weights are given explicitly, the cost function is merely
a table lookup.  When the edge lengths are generated randomly, they
save space
@@<Module subroutines@@>=
length_t cost_from_matrix(const int i,const int j);

length_t
cost_from_matrix(const int i,const int j)
{
	return (length_t)p->edge_weights[i][j];
}

length_t cost_from_short_matrix(const int i, const int j);

length_t
cost_from_short_matrix(const int i, const int j)
{
	return (length_t)p->short_edge_weights[i][j];
}


@@ We use the following cost function when the weight of an edge weights between
two nodes is the two dimensional Euclidean distance between those nodes.

The original incarnation of the TSPLIB documentation used the |aint| 
function in place of the |floor| 
function.  However, |aint| is not in any ANSI~C standard library.  
For non-negative arguments, |aint| and |floor| behave the same way: each
returns
a |double| value representing the greatest integer not more than its |double|
argument.
(I compared the |aint| manual page under SunOS~4.1 dated 15 October 1987
with the description of |floor| in Plauger's {\sl The Standard C Library}.)

The newer documentation uses |nint| which also isn't in the ANSI~C library.
I believe it is the same ``round to the nearest integer'', which I take
to mean |nint(\alpha)=floor(\alpha)|.


Maybe we should be more intelligent about integer square roots???

@@<Module subroutines@@>=
length_t cost_from_euc2d(const int i, const int j);
length_t
cost_from_euc2d(const int i, const int j)
{
	coord_2d *coord_array = p->coord;
	double xd = coord_array[i].x[0] - coord_array[j].x[0];
	double yd = coord_array[i].x[1] - coord_array[j].x[1];
	return (length_t) floor(0.5+my_hypot(xd,yd));
}

@@ Taking the ceiling instead of just rounding to the nearest integer
is interesting because the Euclidean distance with ceiling is a metric.
In particular, the triangle inequality is preserved under taking ceilings.

@@<Module subroutines@@>=
length_t cost_from_ceil2d(const int i, const int j);
length_t
cost_from_ceil2d(const int i, const int j)
{
	coord_2d *coord_array = p->coord;
	double xd = coord_array[i].x[0] - coord_array[j].x[0];
	double yd = coord_array[i].x[1] - coord_array[j].x[1];
	return (length_t) ceil(my_hypot(xd,yd));
}

@@ I converted over to using the |hypot| function.  But things were running
about 15% slower on an SGI Challenge (150MHz MIPS R4400), 
even using the native C compiler
{\tt cc -xansi}.  The {\tt -xansi} option makes some floating point 
intrinsics inline.  This should be faster, but one should always measure.

So I make the use of |hypot| optional via the |my_hypot| macro.  The
macro definition is changed according to the compile time option
|COST_USE_HYPOT|.

Note that the parameters to this macro should not have side effects.

@@<Module definitions@@>=
#if defined(COST_USE_HYPOT)
#define my_hypot(A,B) (hypot((A),(B)))
#else
#define my_hypot(A,B) (sqrt((A)*(A)+(B)*(B)))
#endif

@@ David Johnson's experiments don't use rounded Euclidean norms.  
If the external variable |noround| is true, then we avoid rounding.

@@<Module subroutines@@>=
length_t cost_from_euc2d_not_rounded(const int i, const int j);
length_t
cost_from_euc2d_not_rounded(const int i, const int j)
{
	coord_2d *coord_array = p->coord;
	double xd = coord_array[i].x[0] - coord_array[j].x[0];
	double yd = coord_array[i].x[1] - coord_array[j].x[1];
	return (length_t) my_hypot(xd,yd);
}

@@
I've had floating point precision problems that I believe are
based on massive cancellation.
So in some experiments I may want to compute the raw distances without