htswanalysis/MACS/lib/gsl/gsl-1.11/siman/siman_tsp.c

   1 /* siman/siman_tsp.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <math.h>
  22 #include <string.h>
  23 #include <stdio.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_rng.h>
  26 #include <gsl/gsl_siman.h>
  27 #include <gsl/gsl_ieee_utils.h>
  28
  29 /* set up parameters for this simulated annealing run */
  30
  31 #define N_TRIES 200             /* how many points do we try before stepping */
  32 #define ITERS_FIXED_T 2000      /* how many iterations for each T? */
  33 #define STEP_SIZE 1.0           /* max step size in random walk */
  34 #define K 1.0                   /* Boltzmann constant */
  35 #define T_INITIAL 5000.0        /* initial temperature */
  36 #define MU_T 1.002              /* damping factor for temperature */
  37 #define T_MIN 5.0e-1
  38
  39 gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
  40                              K, T_INITIAL, MU_T, T_MIN};
  41
  42 struct s_tsp_city {
  43   const char * name;
  44   double lat, longitude;        /* coordinates */
  45 };
  46 typedef struct s_tsp_city Stsp_city;
  47
  48 void prepare_distance_matrix(void);
  49 void exhaustive_search(void);
  50 void print_distance_matrix(void);
  51 double city_distance(Stsp_city c1, Stsp_city c2);
  52 double Etsp(void *xp);
  53 double Mtsp(void *xp, void *yp);
  54 void Stsp(const gsl_rng * r, void *xp, double step_size);
  55 void Ptsp(void *xp);
  56
  57 /* in this table, latitude and longitude are obtained from the US
  58    Census Bureau, at http://www.census.gov/cgi-bin/gazetteer */
  59
  60 Stsp_city cities[] = {{"Santa Fe",    35.68,   105.95},
  61                       {"Phoenix",     33.54,   112.07},
  62                       {"Albuquerque", 35.12,   106.62},
  63                       {"Clovis",      34.41,   103.20},
  64                       {"Durango",     37.29,   107.87},
  65                       {"Dallas",      32.79,    96.77},
  66                       {"Tesuque",     35.77,   105.92},
  67                       {"Grants",      35.15,   107.84},
  68                       {"Los Alamos",  35.89,   106.28},
  69                       {"Las Cruces",  32.34,   106.76},
  70                       {"Cortez",      37.35,   108.58},
  71                       {"Gallup",      35.52,   108.74}};
  72
  73 #define N_CITIES (sizeof(cities)/sizeof(Stsp_city))
  74
  75 double distance_matrix[N_CITIES][N_CITIES];
  76
  77 /* distance between two cities */
  78 double city_distance(Stsp_city c1, Stsp_city c2)
  79 {
  80   const double earth_radius = 6375.000; /* 6000KM approximately */
  81   /* sin and cos of lat and long; must convert to radians */
  82   double sla1 = sin(c1.lat*M_PI/180), cla1 = cos(c1.lat*M_PI/180),
  83     slo1 = sin(c1.longitude*M_PI/180), clo1 = cos(c1.longitude*M_PI/180);
  84   double sla2 = sin(c2.lat*M_PI/180), cla2 = cos(c2.lat*M_PI/180),
  85     slo2 = sin(c2.longitude*M_PI/180), clo2 = cos(c2.longitude*M_PI/180);
  86
  87   double x1 = cla1*clo1;
  88   double x2 = cla2*clo2;
  89
  90   double y1 = cla1*slo1;
  91   double y2 = cla2*slo2;
  92
  93   double z1 = sla1;
  94   double z2 = sla2;
  95
  96   double dot_product = x1*x2 + y1*y2 + z1*z2;
  97
  98   double angle = acos(dot_product);
  99
 100   /* distance is the angle (in radians) times the earth radius */
 101   return angle*earth_radius;
 102 }
 103
 104 /* energy for the travelling salesman problem */
 105 double Etsp(void *xp)
 106 {
 107   /* an array of N_CITIES integers describing the order */
 108   int *route = (int *) xp;
 109   double E = 0;
 110   unsigned int i;
 111
 112   for (i = 0; i < N_CITIES; ++i) {
 113     /* use the distance_matrix to optimize this calculation; it had
 114        better be allocated!! */
 115     E += distance_matrix[route[i]][route[(i + 1) % N_CITIES]];
 116   }
 117
 118   return E;
 119 }
 120
 121 double Mtsp(void *xp, void *yp)
 122 {
 123   int *route1 = (int *) xp, *route2 = (int *) yp;
 124   double distance = 0;
 125   unsigned int i;
 126
 127   for (i = 0; i < N_CITIES; ++i) {
 128     distance += ((route1[i] == route2[i]) ? 0 : 1);
 129   }
 130
 131   return distance;
 132 }
 133
 134 /* take a step through the TSP space */
 135 void Stsp(const gsl_rng * r, void *xp, double step_size)
 136 {
 137   int x1, x2, dummy;
 138   int *route = (int *) xp;
 139
 140   step_size = 0 ; /* prevent warnings about unused parameter */
 141
 142   /* pick the two cities to swap in the matrix; we leave the first
 143      city fixed */
 144   x1 = (gsl_rng_get (r) % (N_CITIES-1)) + 1;
 145   do {
 146     x2 = (gsl_rng_get (r) % (N_CITIES-1)) + 1;
 147   } while (x2 == x1);
 148
 149   dummy = route[x1];
 150   route[x1] = route[x2];
 151   route[x2] = dummy;
 152 }
 153
 154 void Ptsp(void *xp)
 155 {
 156   unsigned int i;
 157   int *route = (int *) xp;
 158   printf("  [");
 159   for (i = 0; i < N_CITIES; ++i) {
 160     printf(" %d ", route[i]);
 161   }
 162   printf("]  ");
 163 }
 164
 165 int main(void)
 166 {
 167   int x_initial[N_CITIES];
 168   unsigned int i;
 169
 170   const gsl_rng * r = gsl_rng_alloc (gsl_rng_env_setup()) ;
 171
 172   gsl_ieee_env_setup ();
 173
 174   prepare_distance_matrix();
 175
 176   /* set up a trivial initial route */
 177   printf("# initial order of cities:\n");
 178   for (i = 0; i < N_CITIES; ++i) {
 179     printf("# \"%s\"\n", cities[i].name);
 180     x_initial[i] = i;
 181   }
 182
 183   printf("# distance matrix is:\n");
 184   print_distance_matrix();
 185
 186   printf("# initial coordinates of cities (longitude and latitude)\n");
 187   /* this can be plotted with */
 188   /* ./siman_tsp > hhh ; grep city_coord hhh | awk '{print $2 "   " $3}' | xyplot -ps -d "xy" > c.eps */
 189   for (i = 0; i < N_CITIES+1; ++i) {
 190     printf("###initial_city_coord: %g %g \"%s\"\n",
 191            -cities[x_initial[i % N_CITIES]].longitude,
 192            cities[x_initial[i % N_CITIES]].lat,
 193            cities[x_initial[i % N_CITIES]].name);
 194   }
 195
 196 /*   exhaustive_search(); */
 197
 198   gsl_siman_solve(r, x_initial, Etsp, Stsp, Mtsp, Ptsp, NULL, NULL, NULL,
 199                   N_CITIES*sizeof(int), params);
 200
 201   printf("# final order of cities:\n");
 202   for (i = 0; i < N_CITIES; ++i) {
 203     printf("# \"%s\"\n", cities[x_initial[i]].name);
 204   }
 205
 206   printf("# final coordinates of cities (longitude and latitude)\n");
 207   /* this can be plotted with */
 208   /* ./siman_tsp > hhh ; grep city_coord hhh | awk '{print $2 "   " $3}' | xyplot -ps -d "xy" > c.eps */
 209   for (i = 0; i < N_CITIES+1; ++i) {
 210     printf("###final_city_coord: %g %g %s\n",
 211            -cities[x_initial[i % N_CITIES]].longitude,
 212            cities[x_initial[i % N_CITIES]].lat,
 213            cities[x_initial[i % N_CITIES]].name);
 214   }
 215
 216   printf("# ");
 217   fflush(stdout);
 218 #if 0
 219   system("date");
 220 #endif /* 0 */
 221   fflush(stdout);
 222
 223   return 0;
 224 }
 225
 226 void prepare_distance_matrix()
 227 {
 228   unsigned int i, j;
 229   double dist;
 230
 231   for (i = 0; i < N_CITIES; ++i) {
 232     for (j = 0; j < N_CITIES; ++j) {
 233       if (i == j) {
 234         dist = 0;
 235       } else {
 236         dist = city_distance(cities[i], cities[j]);
 237       }
 238       distance_matrix[i][j] = dist;
 239     }
 240   }
 241 }
 242
 243 void print_distance_matrix()
 244 {
 245   unsigned int i, j;
 246
 247   for (i = 0; i < N_CITIES; ++i) {
 248     printf("# ");
 249     for (j = 0; j < N_CITIES; ++j) {
 250       printf("%15.8f   ", distance_matrix[i][j]);
 251     }
 252     printf("\n");
 253   }
 254 }
 255
 256 /* [only works for 12] search the entire space for solutions */
 257 static double best_E = 1.0e100, second_E = 1.0e100, third_E = 1.0e100;
 258 static int best_route[N_CITIES];
 259 static int second_route[N_CITIES];
 260 static int third_route[N_CITIES];
 261 static void do_all_perms(int *route, int n);
 262
 263 void exhaustive_search()
 264 {
 265   static int initial_route[N_CITIES] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11};
 266   printf("\n# ");
 267   fflush(stdout);
 268 #if 0
 269   system("date");
 270 #endif
 271   fflush(stdout);
 272   do_all_perms(initial_route, 1);
 273   printf("\n# ");
 274   fflush(stdout);
 275 #if 0
 276   system("date");
 277 #endif /* 0 */
 278   fflush(stdout);
 279
 280   printf("# exhaustive best route: ");
 281   Ptsp(best_route);
 282   printf("\n# its energy is: %g\n", best_E);
 283
 284   printf("# exhaustive second_best route: ");
 285   Ptsp(second_route);
 286   printf("\n# its energy is: %g\n", second_E);
 287
 288   printf("# exhaustive third_best route: ");
 289   Ptsp(third_route);
 290   printf("\n# its energy is: %g\n", third_E);
 291 }
 292
 293 /* James Theiler's recursive algorithm for generating all routes */
 294 static void do_all_perms(int *route, int n)
 295 {
 296   if (n == (N_CITIES-1)) {
 297     /* do it! calculate the energy/cost for that route */
 298     double E;
 299     E = Etsp(route);            /* TSP energy function */
 300     /* now save the best 3 energies and routes */
 301     if (E < best_E) {
 302       third_E = second_E;
 303       memcpy(third_route, second_route, N_CITIES*sizeof(*route));
 304       second_E = best_E;
 305       memcpy(second_route, best_route, N_CITIES*sizeof(*route));
 306       best_E = E;
 307       memcpy(best_route, route, N_CITIES*sizeof(*route));
 308     } else if (E < second_E) {
 309       third_E = second_E;
 310       memcpy(third_route, second_route, N_CITIES*sizeof(*route));
 311       second_E = E;
 312       memcpy(second_route, route, N_CITIES*sizeof(*route));
 313     } else if (E < third_E) {
 314       third_E = E;
 315       memcpy(route, third_route, N_CITIES*sizeof(*route));
 316     }
 317   } else {
 318     int new_route[N_CITIES];
 319     unsigned int j;
 320     int swap_tmp;
 321     memcpy(new_route, route, N_CITIES*sizeof(*route));
 322     for (j = n; j < N_CITIES; ++j) {
 323       swap_tmp = new_route[j];
 324       new_route[j] = new_route[n];
 325       new_route[n] = swap_tmp;
 326       do_all_perms(new_route, n+1);
 327     }
 328   }
 329 }