tsp.c

Algorithm based on the outline of (Lin, 1965) for unix

/** @(#)tsp 970101 KLD
 *  Demonstration of minimumtour (solves Travelling Salesman Problem)
 *
 *  This test case ...
 *     node:  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
 *        x:  2  7 13 17 16 15 12  2  2  0  4 10 14 15 13  8  7  5  6 10
 *        y:  2  1  0  3 13 15 16 18 12 10  6  4  5  9 12 13 14 13  8  9
 *  ... should result in path visiting nodes in this order:
 *           16 17 18  8  9 10  1 11 19 20 12  2  3  4 13 14 15  5  6  7
 *
 *  Usage:
 *     tsp x1 y1 x2 y2 x3 y3 ...
 *  There may be up to 20 coordinate pairs, that is, 40 arguments.  The
 *  coordinates are integers with values in the range [0,1023].
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXNODES       (20)
#define TMX            (MAXNODES / 5)
#define COSTMATSIZE    ((MAXNODES * (MAXNODES-1)) / 2)
#define	MAX_XY_VAL     (1023)
#define COSTSCALE      (10.0)
#define SEED           (1027.0)

#define MN1 (MAXNODES+1)	/* Used in declaring tour arrays */

double  rannum;
int     dist[COSTMATSIZE + 1];
int     tour[MN1];
int     nn;
char    *me;

/* Local prototypes */
double  randreal(double *rnx);
void    randomtour(int n);
int     cost(int i, int j);
void    exchangelink(int sense, int k, int j, int n);
void    minimumtour(int n, int *path, int *x, int *y);
void	usage(char *msg);

int     main(int argc, char **argv)
{
  int     x[MN1], y[MN1], p[MN1];
  int     i, coord;

  me = argv[0];
  if ( (argc < 3) || (argc > 41) || !(argc & 1) ) {
    usage("ERROR: Too many, too few, or odd number of arguments");
    exit(1);
  }
  for (i = 2; i <= argc; ++i) {
    if (((coord = atoi(argv[i - 1])) < 0) || (coord > MAX_XY_VAL)) {
      usage("ERROR: Coordinate values must be in [0,1023]");
      exit(1);
    }
    if (i & 1) {
      y[i / 2] = coord;			/* Set up a path ordinate */
    } else {
      x[i / 2] = coord;			/* Set up a path abscissa */
    }
  }
  minimumtour(argc / 2, p, x, y);	/* Get path indices in p array */
  for (i = 1; i <= argc / 2; ++i) {	/* Write out the optimized path */
    printf(" %d", p[i]);
  }
  printf("\n");
  return (0);
}

void	usage(char *msg)
{
  printf("%s\n", msg);
  printf("\nUsage:\n");
  printf("   %s x1 y1 x2 y2 x3 y3 ...\n", me);
  printf("There may be up to 20 coordinate pairs, that is, 40 arguments.\n");
  printf("The coordinates are integers with values in the range [0,1023].\n");
  printf("The output is a list of node numbers identifying the pairs\n");
  printf("in a nearly optimal (minimum cost) order.\n");
  exit(1);
}

/**
  Minimumtour presents in array 'path' a sequence of integers
  in [1,n], which describe an optimal (or nearly optimal)
  path through the nodes given by the arrays x and y.

  Algorithm based on the outline in Appendix B of:
     Computer Solutions of the Traveling Salesman problem
     Shen Lin
     Bell System Technical Journal 44(1965), pp. 2245-2269

  Arguments:
     n:            Number of nodes
     path:         Node index list (returned resulting optimal path)
     x,y:          Coordinate arrays of nodes listed in 'path'

  Called procedures:
     randomtour:   Generates a random path
     exchangelink: Swaps segments of the path
     cost:         Figures cost of a link in the path
     randreal:     Generates a random number in (0,1)
     usage:        Reports a fatal error and exits to op system

  Parameters required:
     MAXNODES:     max no. of nodes to tour
     TMX:          MAXNODES / 5
     COSTMATSIZE:  (MAXNODES * (MAXNODES-1)) / 2
     COSTSCALE:    scale factor for cost analyzer
     SEED:         start value for random number generator
*/

void    minimumtour(int n, int *path, int *x, int *y)
{
  int     i, j, k, m, r, count;
  int     nobetter, sense;
  int     d1, d2, da, addedcost, removedcost, cp;
  double  xd, yd;
  int     c[MN1];
  int     trial[TMX + 1][MN1];

  if (n > MAXNODES)
    usage("ERROR: Too many nodes");
  if (n < 3) {			/* Handle one-node and two-node cases */
    path[1] = 1;
    path[2] = 2;
  } else {			/* Handle cases having three or more nodes */
    rannum = SEED;
    if (n < 5)
      r = 1;
    else
      r = n / 5;
    nn = n + n;
    for (i = 1; i <= n - 1; ++i) {	/* Set up cost matrix */
      m = ((nn - i) * (i - 1)) / 2 - i;
      for (j = i + 1; j <= n; ++j) {
	xd = x[i] - x[j];
	yd = y[i] - y[j];
	da = sqrt((xd * xd) + (yd * yd)) * COSTSCALE;
	dist[m + j] = da;
	if (da < 2) {
	  printf("ERROR: Node given more than once\n");
	}
      }				/* j */
    }				/* i */

    for (m = 1; m <= r; ++m) {	/* Minimize r random paths */
      randomtour(n);
      if (n > 3)
	do {			/* Paths having <= 3 nodes already minimal */
	  nobetter = 1;
	  for (count = 1; count <= n; ++count) {	/* Try all rotations */
	    k = 1;
	    while ((k <= (n - 3)) && nobetter) {	/* Pick a series of
							 * links */
	      j = k + 1;
	      while ((j <= (n - 1)) && nobetter) {	/* Pick a place to go */
		d2 = cost(1, j + 1) + cost(k, j);
		d1 = cost(k, j + 1) + cost(1, j);
		sense = d1 <= d2;
		if (sense)
		  da = d1;
		else
		  da = d2;
		addedcost = da + cost(k + 1, n);
		removedcost = cost(1, n) + cost(k, k + 1) + cost(j, j + 1);
		nobetter = addedcost >= removedcost;
		if (!nobetter)
		  exchangelink(sense, k, j, n);
		++j;
	      }			/* j */
	      ++k;
	    }			/* k */
	    if (nobetter) {	/* Can't improve it more: rotate path */
	      k = tour[n];
	      for (j = n; j >= 2; --j)
		tour[j] = tour[j - 1];
	      tour[1] = k;
	    }
	  }			/* count */
	} while (!nobetter);	/* if (n > 3) do ... */
      cp = cost(1, n);
      trial[m][1] = tour[1];
      for (i = 2; i <= n; ++i) { /* Save trial path */
	cp += cost(i - 1, i);
	trial[m][i] = tour[i];
      }				/* i */
      c[m] = cp;		/* Cost of mth path */
    }				/* m */
    j = 1;			/* Return the cheapest of the minimized trial
				 * paths */
    for (m = 1; m <= r; ++m)
      if (c[m] > c[j])
	j = m;
    for (i = 1; i <= n; ++i)
      path[i] = trial[j][i];
  }				/* n NOT < 3 */
}				/* minimumtour */

void    screwup(int errinx)
{
  switch (errinx) {
  case 1:
    printf("ERROR: Too many nodes in the list (max = %d)\n", MAXNODES);
    break;
  case 2:
    printf("ERROR: Logical error in randomtour\n");
    break;
  default:
    printf("ERROR: Undefined error!\n");
  }
  exit(1);
}				/* screwup */

#define RANMPY       1021.0
#define RANMODULUS   2048.0

double  randreal(double *rnx)
/** Produces pseudorandom numbers in the interval (0,1).
  Argument must initially be non-zero; will be updated   */
{
  *rnx *= RANMPY;
  *rnx -= RANMODULUS * ((long) (*rnx / RANMODULUS));
  return (*rnx / RANMODULUS);
}				/* randreal */

void    randomtour(int n)
/* Generates a random sequence of the integers 1..n in global array 'tour' */
{
  int     i, j, k, m;

  for (i = 1; i <= n; ++i)
    tour[i] = 0;
  for (i = 1; i <= n; ++i) {
    j = 1 + (randreal(&rannum) * (n + 1 - i));	/* Random j in (1, n + 1 - i) */
    m = 0;
    k = 0;
    while (m < j) {		/* Take the jth zeroed tour[k] next */
      ++k;
      if (k > n)
	usage("ERROR: Logical error in randomtour");
      if (!tour[k])
	++m;
    }
    tour[k] = i;
  }				/* i */
}				/* randomtour */

int     cost(int i, int j)
/* Computes cost of link tour[i] --> tour[j] using array dist */
{
  int     it, jt, tem;

  it = tour[i];
  jt = tour[j];
  if (it > jt) {
    tem = it;
    it = jt;
    jt = tem;
  }
  return (dist[((nn - it) * (it - 1)) / 2 - it + jt]);
}				/* cost */

/**
  Exchangelink moves the series at the beginning of 'tour' to
  a specified new place in 'tour', and rotates the result so it
  will begin at the second node following the moved series.

  Ex: Given tour  = 1  2  3  4  5  6  7  8  9 10, n=10, k=3, j=7:
      Output tour = 9 10  4  5  6  7  1  2  3  8
      Same tour, n, k, and j but with sense==FALSE:
      Output tour = 9 10  4  5  6  7  3  2  1  8

  sense: TRUE:normal, FALSE:reverse travel through moved series
  k:     Last node of series to move (first is 1)
  j:     Series gets moved to between nodes j and j + 1
  n:     Number of nodes in the (closed) tour
*/

void    exchangelink(int sense, int k, int j, int n)
{
  int     temptour[MN1];
  int     i, m;

  for (i = 1; i <= n; ++i)
    temptour[i] = tour[i];
  m = 1;
  if (j <= (n - 2)) {
    for (i = j + 2; i <= n; ++i) {	/* 2nd node after new */
      tour[m] = temptour[i];
      ++m;
    }
  }
  for (i = k + 1; i <= j; ++i) { /* Patch around old */
    tour[m] = temptour[i];
    ++m;
  }
  if (sense) {
    for (i = 1; i <= k; ++i) {	/* Forward - walk new */
      tour[m] = temptour[i];
      ++m;
    }
  } else {
    for (i = k; i >= 1; --i) {	/* Reverse - walk new */
      tour[m] = temptour[i];
      ++m;
    }
  }
  tour[m] = temptour[j + 1];	/* 1st node after new */
}				/* exchangelink */