read.w,v

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

even casting to |length_t|.  This is done by function 
|cost_from_euc2d_raw|.

The declarations of numerical functions
|floor|, |sqrt|, and |hypot| are already included from \file{math.h}.

@@<Subroutines@@>=
double
cost_from_euc2d_raw(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 my_hypot(xd,yd);
}

@@
@@<Exported subroutines@@>=
double cost_from_euc2d_raw(const int i, const int j);


@@ For comparisons between edges, we can strip off the rounding and
square root computations.  This is about three times faster
on an Intel 486 DX2-66.

@@<Module subroutines@@>=
length_t pseudo_cost_from_euc2d(const int i, const int j);
length_t
pseudo_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)(xd*xd+yd*yd);	/* Worry about rounding? */
}

@@ The GEO type is a geographic great circle distance between 
two longitude and latitude pairs.


Henry Baker suggests using verisine and related functions 
in order to keep better precision for arccos near zero.
See his ftp site.  But we keep the standard algorithm as specified
in the TSPLIB paper.

We use the macro |truncate_to_zero| to replace what is written as
|nint| in the original TSPLIB paper.

GCC is of the opinion that strict
ANSI C doesn't define |M_PI| in \file{math.h}.
No matter, at least $\pi$ is a constant.

@@d truncate_to_zero(x) ((double)(int)(x))
@@d PI (3.14159265358979323846)

@@<Module subroutines@@>=
length_t cost_from_geo(const int i, const int j);
length_t
cost_from_geo(const int i, const int j)
{
	coord_2d *coord_array = p->coord;
	const double x1 = coord_array[i].x[0], x2 = coord_array[j].x[0];
	const double y1 = coord_array[i].x[1], y2 = coord_array[j].x[1];
    double degrees, minutes;
	double latitude1, latitude2, longitude1, longitude2;
    double q1, q2, q3;
	const double radius=6378.388;

    degrees = truncate_to_zero(x1);
    minutes = x1 - degrees;
    latitude1 = (PI/180) * (degrees + minutes * (5.0/3.0));

    degrees = truncate_to_zero(x2);
    minutes = x2 - degrees;
    latitude2 = (PI/180) * (degrees + minutes * (5.0/3.0));

    degrees = truncate_to_zero (y1);
    minutes = y1 - degrees;
    longitude1 = (PI/180) * (degrees + minutes * (5.0/3.0));

    degrees = truncate_to_zero (y2);
    minutes = y2 - degrees;
    longitude2 = (PI/180) * (degrees + minutes * (5.0/3.0));

    q1 = cos (longitude1 - longitude2);
    q2 = cos (latitude1 - latitude2);
    q3 = cos (latitude1 + latitude2);
    return (length_t) 
		truncate_to_zero(radius * acos(0.5*((1+q1)*q2 - (1-q1) * q3)) + 1);
}


@@ Cost function ATT is just a bit funky.  
It is much like a scaled Euclidean metric.

That funky constant is $1/\sqrt{10}$.  (Thanks, \type{bc}.)

@@<Module subroutines@@>=
length_t cost_from_att(const int i, const int j);
length_t
cost_from_att(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)*0.31622776601683793319988935);
}

@@ Cost function |DSJ_RANDOM| generates random 
scaled weights based only on a scale, a seed, and the two indices $i$
and $j$.  

Remember that our vertex numbers are zero-based, so we have to add 1
before generating the salt.

@@<Module subroutines@@>=

length_t cost_from_dsj_random(const int i, const int j);
length_t
cost_from_dsj_random(const int ii, const int jj)
{
	const int32 i=ii, j=jj;
	const int32 salt1 = 0x12345672*(i+1) + 1;
	const int32 salt2 = 0x12345672*(j+1) + 1;
	int32 x,y,z;

    x = salt1 & salt2;
    y = salt1 | salt2;
    z = dsj_random_param;

    x *= z;
    y *= x;
    z *= y;

    z ^= dsj_random_param;

    x *= z;
    y *= x;
    z *= y;

    x = ((salt1 + salt2) ^ z) &0x7fffffff;
    return (length_t) (x*dsj_random_factor);
}

@@ The two parameters are also part of the cost function context.
@@<Set the |cost|...@@>=
dsj_random_factor = p->dsj_random_factor;
dsj_random_param = p->dsj_random_param;

@@ For consistent results, we need a 32 bit integer type.
@@<Module definitions@@>=
#if SIZEOF_INT==4
typedef int int32;
#elif SIZEOF_SHORT==4
typedef short int int32;
#elif SIZEOF_LONG==4
typedef long int32;
#else
#error "I need a 32 bit integer for consistent results with DSJ_RANDOM"
#endif

@@ We choose among the above cost functions based upon the type
of the problem.  In the Euclidean case, we pay attention to the
|noround| variable.

@@<Set the |cost| function according to the current context@@>=
switch( p->edge_weight_type ) {
case EXPLICIT: 
	cost = cost_from_matrix;
	pseudo_cost = cost_from_matrix; 
	break;
case EUC_2D: 
	{
	extern int noround;	/* Defined by module \module{LK}. */
	cost = noround ? cost_from_euc2d_not_rounded : cost_from_euc2d;
	}
	pseudo_cost = pseudo_cost_from_euc2d; 
	break;
case CEIL_2D: 
	{
	extern int noround;	/* Defined by module \module{LK}. */
	cost = noround ? cost_from_euc2d_not_rounded : cost_from_ceil2d;
	}
	pseudo_cost = pseudo_cost_from_euc2d; 
	break;
case GEO: 
	cost = cost_from_geo;
	pseudo_cost = cost_from_geo; 
	break;
case ATT: 
	cost = cost_from_att;
	pseudo_cost = cost_from_att; 
	break;
case DSJ_RANDOM: 
	cost = cost_from_dsj_random;
	pseudo_cost = cost_from_dsj_random; 
	break;
case RANDOM_EDGES: 
	errorif(p->short_edge_weights==NULL,"RANDOM_EDGES specified but no SEED given");
	cost = cost_from_short_matrix;
	pseudo_cost = cost_from_short_matrix; 
	break;
case NO_EDGE_TYPE:
default:
	errorif(1,"Switching to an instance with unknown edge type %d",
		p->edge_weight_type);
}

@@*Forcing even number of vertices.
Some problems require that the number of vertices be even.
For instance, perfect matchings only exist when there are an even number
of vertices.

@@<Possibly remove a vertex@@>=
if ( force_even_num && (p->input_n % 2) ) {
	@@<Remove a vertex@@>@@;
}

@@ There are two cases: geometric and non-geometric.

For a non-geometric instance, we just remove the last vertex.

For a geometric instance, we follow David Applegate and William Cook,
``Solving Large-Scale Matching Problems'', in {\sl Network Flows and
Matching}, First DIMACS Implementation Challenge, David S.~Johnson and
Catherine C.~McGeoch, editors, 1993, pages 557--576.
They remove
the vertex whose coordinates are lexicographically last.
That is, find the vertices with maximum $x$ coordinate.  From
among those, pick one with the maximum $y$ coordinate.  
If there's still a tie, pick the last-occuring one.
(This is a little more spelled-out than what they said: ``we sorted
the $x,y$ coordinates and deleted the last point.'')
@@^Applegate, David@@>
@@^Cook, William@@>
@@^Johnson, David S.@@>
@@^McGeoch, Catherine C.@@>

In the geometric case we {\it swap\/}
the lexicographic-coordinate last point with the
physically-stored last coordinate so that we can always reconstruct an
output that is 
equivalent to the original input.

@@<Remove a vertex@@>=
p->n--;
switch ( p->edge_weight_type ) {
case EUC_2D:
case CEIL_2D:
case GEO:
case ATT:
	{ int i, lex_last=0;
	for (i=1;i<p->input_n;i++) {
		if ( p->coord[i].x[0] > p->coord[lex_last].x[0] @@|
			|| (p->coord[i].x[0] == p->coord[lex_last].x[0]
			   && p->coord[i].x[1] >= p->coord[lex_last].x[1]) )
			lex_last = i;
	}
	{
	coord_2d t = p->coord[p->n]; /* Now swap them. */
	p->coord[p->n] = p->coord[lex_last];
	p->coord[lex_last] = t;
	}
	}
default: break;  /* We already decremented |p->n|. */
}

@@*Checking consistency.
If the arrays are set up right then the cost functions should also be set
correctly.

@@<Check consistency@@>=
if (p->n) {
	switch(p->edge_weight_type) {
	case EUC_2D:
	case CEIL_2D:
	case GEO:
		errorif(p->edge_weight_format != NO_EDGE_FORMAT,
			"Edge weight format should be NO_EDGE_FORMAT, but is %s",
				edge_weight_format_name[p->edge_weight_format]);
		errorif(p->coord == NULL,"No coordinates were read");
		break;
	case EXPLICIT:
		errorif(p->edge_weight_format==NO_EDGE_FORMAT,
			"Edge weight shouldn't be NO_EDGE_FORMAT, but it is");
		errorif(p->edge_weights == NULL,"No edge weights were read");
		break;
	case DSJ_RANDOM:
	default:
		break;
	}
}

@@*Debugging output routines.

%When debugging output is desired, |verbose>0|.

When we're done reading the file, we output some summary statistics
about its contents to a debugging file.  The file handle is |debug|, passed
in to |read_tsp_file| as a parameter.

@@ First we output the number of cities and create 1-based arrays
for the |x| and |y| coordinates.

Unfortunately, the maximum array size on some PostScript implementations
is  quite small.  For example, the default in the Ghostscript implementation
is 8191.  If we have too many cities, then we can't create the arrays.
In these cases, we must take care later on to not use any PostScript routines
that require them.

@@<Debugging output@@>=
if (debug) {
	fprintf(debug,"/N {%d} def        %% number of nodes\n",p->input_n);
	if ( p->input_n < 8191 )
		fprintf(debug,"/xs N 1 add array def\n/ys N 1 add array def\n");
}


@@ Then we output the bounding positions.

@@<Debugging output@@>=
if (debug) {
const double xdiff = ((double) p->xmax) - ((double)p->xmin);
const double ydiff = ((double) p->ymax) - ((double)p->ymin);
const double maxrange = max(xdiff,ydiff);
const double xoffset = xdiff<ydiff ? (maxrange-xdiff)/2 : 0.0;
const double yoffset = ydiff<xdiff ? (maxrange-ydiff)/2 : 0.0;
fprintf(debug,"/xmin {%f} def \n",p->xmin);
fprintf(debug,"/xmax {%f} def \n",p->xmax);
fprintf(debug,"/ymin {%f} def \n",p->ymin);
fprintf(debug,"/ymax {%f} def \n",p->ymax);
fprintf(debug,"/maxrange {%f} def \n", maxrange );
fprintf(debug,"/xoffset  {%f} def \n", xoffset);
fprintf(debug,"/yoffset  {%f} def \n", yoffset);
}

@@ And the coordinates themselves.
@@<Debugging output@@>=
if (debug) {int i;
for (i=0; i<p->input_n ; i++) {
	fprintf(debug,"%f %f ts\n",p->coord[i].x[0],p->coord[i].x[1]);
}
fflush(debug);
}

@@*Writing TSPLIB files.
Given a |tsp_instance_t| object, we'd like to write it out again.

We have two routines.  Procedure |write_tsp_file|

@@<Subroutines@@>=
void
write_tsp_file(tsp_instance_t *tsp,FILE *out)
{
	write_tsp_file_clip(tsp, out, 0);
}

void write_tsp_file_clip(tsp_instance_t *tsp, FILE *out, int force_even_num)
{
	write_tsp_file_clip_matrix(tsp, out, force_even_num, 0);
}

void write_tsp_file_clip_matrix(
	tsp_instance_t *tsp, FILE *out, int force_even_num, int force_matrix)
{
int n=tsp->input_n;
if ( (n%2) && force_even_num ) n--;

if(out) {
fprintf(out,"NAME: %s\n",tsp->name);
fprintf(out,"TYPE: TSP\n");
fprintf(out,"COMMENT: %s%s\n",
	(tsp->comment?tsp->comment:""),
	(n<tsp->input_n?"| clip":""));
fprintf(out,"DIMENSION: %d\n",n);
if ( force_matrix ) {
	int format=UPPER_ROW; 
	fprintf(out,"EDGE_WEIGHT_TYPE: EXPLICIT\n");
	@@<Write the explicit weights in format |format|@@>@@;
} else {
	switch(tsp->edge_weight_type) {
	case CEIL_2D: fprintf(out,"EDGE_WEIGHT_TYPE: CEIL_2D\n"); 
		@@<Write the coordinates@@>@@;
		break;
	case EUC_2D:  fprintf(out,"EDGE_WEIGHT_TYPE: EUC_2D\n"); 
		@@<Write the coordinates@@>@@;
		break;
	case GEO:  fprintf(out,"EDGE_WEIGHT_TYPE: GEO\n"); 
		@@<Write the coordinates@@>@@;
		break;
	case ATT:  fprintf(out,"EDGE_WEIGHT_TYPE: ATT\n"); 
		@@<Write the coordinates@@>@@;
		break;
	case DSJ_RANDOM:  fprintf(out,"EDGE_WEIGHT_TYPE: DSJ_RANDOM\n"); 
		fprintf(out,"SCALE: %f\nSEED: %ld",tsp->scale,tsp->seed);
		break;
	case EXPLICIT: fprintf(out,"EDGE_WEIGHT_TYPE: EXPLICIT\n");
		{ int format=tsp->edge_weight_format;
		@@<Write the explicit weights in format |format|@@>@@;
		}
		break;
	case RANDOM_EDGES: fprintf(out,"EDGE_WEIGHT_TYPE: EXPLICIT\n");
		{ int format=UPPER_ROW;
		@@<Write the explicit weights in format |format|@@>@@;
		}
		break;
	case NO_EDGE_TYPE:
	default:
		errorif(1,"No edge type specified by instance: have %d instead",
			tsp->edge_weight_type);
	}
}
fprintf(out,"EOF\n");
}
}

@@
@@<Exported subroutines@@>=
void write_tsp_file(tsp_instance_t *tsp,FILE *out);
void write_tsp_file_clip(tsp_instance_t *tsp,FILE *out, int force_even_num);
void write_tsp_file_clip_matrix(
	tsp_instance_t *tsp,FILE *out, int force_even_num, int force_matrix);

@@
@@<Write the coordinates@@>=
fprintf(out,"NODE_COORD_SECTION\n");
{int i;
for (i=0;i<n;i++) {
fprintf(out,"%d %g %g\n",i+1,tsp->coord[i].x[0],tsp->coord[i].x[1]);
}
}


@@
@@<Write the explicit weights in format |format|@@>=
{ int row, col;
	tsp_instance_t *old_p = switch_to(tsp);
	switch (format) {
	case LOWER_DIAG_ROW:
		fprintf(out,"EDGE_WEIGHT_FORMAT: LOWER_DIAG_ROW\n");
		fprintf(out,"EDGE_WEIGHT_SECTION\n");
		for (row=0;row<n;row++) {
			for (col=0;col<=row;col++) {
				fprintf(out," %ld",(long)cost(row,col));
			}
			fprintf(out,"\n");
		}
		break;
	case FULL_MATRIX:
		fprintf(out,"EDGE_WEIGHT_FORMAT: FULL_MATRIX\n");
		fprintf(out,"EDGE_WEIGHT_SECTION:\n");
		for (row=0;row<n;row++) {
			for (col=0;col<n;col++) {
				fprintf(out," %ld",(long)cost(row,col));
			}
			fprintf(out,"\n");
		}
		break;
	case UPPER_ROW:
		fprintf(out,"EDGE_WEIGHT_FORMAT: UPPER_ROW\n");
		fprintf(out,"EDGE_WEIGHT_SECTION:\n");
		for (row=0;row<n;row++) {
			for (col=row+1;col<n;col++) {
				fprintf(out," %ld",(long)cost(row,col));
			}
			fprintf(out,"\n");
		}
		break;
	default: 
		errorif(1,"Unknown explicit format %d\n",format);
	}
	switch_to(old_p);
}
	

@@*Index.
@


1.142
log
@Support DSJ RANDOM, ATT, and GEO
Better error messages.
@
text
@d52 4
d195 1
a195 1
const char *read_rcs_id="$Id: read.w,v 1.141 1998/10/08 15:14:47 neto Exp neto $";
d226 1
a226 1
switch_to(tsp_instance_t *new_problem) 
d1196 6
d1212 3
a1214 18
switch(tsp->edge_weight_type) {
case CEIL_2D: fprintf(out,"EDGE_WEIGHT_TYPE: CEIL_2D\n"); 
	@@<Write the coordinates@@>@@;
	break;
case EUC_2D:  fprintf(out,"EDGE_WEIGHT_TYPE: EUC_2D\n"); 
	@@<Write the coordinates@@>@@;
	break;
case GEO:  fprintf(out,"EDGE_WEIGHT_TYPE: GEO\n"); 
	@@<Write the coordinates@@>@@;
	break;
case ATT:  fprintf(out,"EDGE_WEIGHT_TYPE: GEO\n"); 
	@@<Write the coordinates@@>@@;
	break;
case DSJ_RANDOM:  fprintf(out,"EDGE_WEIGHT_TYPE: DSJ_RANDOM\n"); 
	fprintf(out,"SCALE: %f\nSEED: %ld",tsp->scale,tsp->seed);
	break;
case EXPLICIT: fprintf(out,"EDGE_WEIGHT_TYPE: EXPLICIT\n");
	{ int format=tsp->edge_weight_format;
d1216 31
a1247 10
	break;
case RANDOM_EDGES: fprintf(out,"EDGE_WEIGHT_TYPE: EXPLICIT\n");
	{ int format=UPPER_ROW;
	@@<Write the explicit weights in format |format|@@>@@;
	}
	break;
case NO_EDGE_TYPE:
default:
	errorif(1,"No edge type specified by instance: have %d instead",
		tsp->edge_weight_type);
d1257 2
d1273 1
d1280 1
a1280 1
				fprintf(out," %ld",(long)(tsp->edge_weights[row][col]));
d1290 1
a1290 1
				fprintf(out," %ld",(long)(tsp->edge_weights[row][col]));
d1300 1
a1300 1
				fprintf(out," %ld",(long)(tsp->edge_weights[row][col]));
d1305 2
a1306 1
	default: break;
d1308 1
@


1.141
log
@Allow clipping of last point on output too.
@