chiark - git - mdw - rocl/blob - graph.c

   1 /* -*-c-*-
   2  *
   3  * Graph theory stuff
   4  *
   5  * (c) 2003 Mark Wooding
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This program is free software; you can redistribute it and/or modify
  11  * it under the terms of the GNU General Public License as published by
  12  * the Free Software Foundation; either version 2 of the License, or
  13  * (at your option) any later version.
  14  *
  15  * This program is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  * GNU General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU General Public License
  21  * along with this program; if not, write to the Free Software Foundation,
  22  * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  23  */
  24
  25 /*----- Header files ------------------------------------------------------*/
  26
  27 #include <assert.h>
  28 #include <math.h>
  29 #include <stdio.h>
  30 #include <stdlib.h>
  31 #include <string.h>
  32
  33 #include <tcl.h>
  34
  35 #include "vec.h"
  36
  37 /*----- Static variables --------------------------------------------------*/
  38
  39 #define INF ((unsigned long)-1)
  40
  41 /*----- Utility functions -------------------------------------------------*/
  42
  43 static int err(Tcl_Interp *ti, /*const*/ char *p)
  44 {
  45   Tcl_SetResult(ti, p, TCL_STATIC);
  46   return (TCL_ERROR);
  47 }
  48
  49 /* --- @import@ --- *
  50  *
  51  * Arguments:   @Tcl_Interp *ti@ = interpreter to leave errors in
  52  *              @vec *v@ = pointer to input adjacency matrix
  53  *              @unsigned long *tt@ = pointer to output adjacency matrix
  54  *              @size_t *nn@ = where to put the table size
  55  *
  56  * Returns:     Tcl return code.
  57  *
  58  * Use:         Imports an adjacency matrix.
  59  */
  60
  61 static int import(Tcl_Interp *ti, vec *v, unsigned long **tt, size_t *nn)
  62 {
  63   size_t i;
  64   unsigned long *t;
  65   size_t n;
  66
  67   /* --- Check the table is well-formed --- */
  68
  69   if (v->ndim != 2)
  70     return (err(ti, "adjacency matrix must be two-dimensional"));
  71   if (v->dim[0].lo != 0 || v->dim[1].lo || v->dim[0].hi != v->dim[1].hi)
  72     return (err(ti, "adjacency matrix must be square and zero-origin"));
  73   n = *nn = v->dim[0].hi;
  74
  75   /* --- Copy the data over --- */
  76
  77   n *= n;
  78   assert(n == v->n);
  79   t = (void *)Tcl_Alloc(n * sizeof(*t));
  80   for (i = 0; i < n; i++) {
  81     long l;
  82     if (Tcl_GetLongFromObj(ti, v->v[i], &l) != TCL_OK) {
  83       Tcl_Free((void *)t);
  84       return (TCL_ERROR);
  85     }
  86     t[i] = l >= 0 ? l : INF;
  87   }
  88   *tt = t;
  89   return (TCL_OK);
  90 }
  91
  92 /* --- @export@ --- *
  93  *
  94  * Arguments:   @Tcl_Interp *ti@ = interpreter to create output vector
  95  *              @unsigned long *t@ = pointer to table
  96  *              @size_t n@ = size of the table
  97  *
  98  * Returns:     A pointer to the vector, or null.
  99  *
 100  * Use:         Exports an adjacency matrix.
 101  */
 102
 103 static vec *export(Tcl_Interp *ti, unsigned long *t, size_t n)
 104 {
 105   vec_bound b[2];
 106   vec *v;
 107   size_t i;
 108   Tcl_Obj *o;
 109
 110   b[0].lo = b[1].lo = 0;
 111   b[0].hi = b[1].hi = n;
 112   if ((v = vec_create(ti, 2, b, 0)) == 0)
 113     return (0);
 114   o = Tcl_NewLongObj(-1);
 115   Tcl_IncrRefCount(o);
 116   for (i = 0; i < v->n; i++) {
 117     v->v[i] = t[i] == INF ? o : Tcl_NewLongObj(t[i]);
 118     Tcl_IncrRefCount(v->v[i]);
 119   }
 120   Tcl_DecrRefCount(o);
 121   return (v);
 122 }
 123
 124 /*----- Floyd-Warshall all-points shortest path ---------------------------*/
 125
 126 /* --- @graph-shortest-path VEC@ --- *
 127  *
 128  * Returns a pair of vectors containing, respectively, the shortest path
 129  * length and the successor element in the shortest path.  If you say
 130  *
 131  *   destructure {len path} [graph-shortest-path $v]
 132  *
 133  * then [$len get I J] is the shortest path length from node I to node J, and
 134  * [$path get I J] is the first hop on that shortest path.  (To compute the
 135  * entire path, set K to be that first hop; the next hop is then [$path get K
 136  * J], and so on.)
 137  *
 138  * The adjacency matrix is given in VEC: negative entries indicate no path;
 139  * nonnegative entries are weights.  All entries must be integers.
 140  */
 141
 142 static int cmd_shortestpath(ClientData cd, Tcl_Interp *ti,
 143                             int objc, Tcl_Obj *const *objv)
 144 {
 145   vec *v, *lv = 0, *pv = 0;
 146   size_t n, i, j, k;
 147   unsigned long *a = 0, *p = 0;
 148   Tcl_Obj *o;
 149
 150   /* --- Read in the arguments --- */
 151
 152   if (objc != 2) {
 153     err(ti, "usage: graph-shortest-path VEC");
 154     goto fail;
 155   }
 156   if ((v = vec_find(ti, objv[1])) == 0 || import(ti, v, &a, &n) != TCL_OK)
 157     goto fail;
 158
 159   /* --- Set up the path table --- */
 160
 161   p = (void *)Tcl_Alloc(n * n * sizeof(*p));
 162   for (i = 0; i < n; i++) {
 163     for (j = 0; j < n; j++)
 164       p[i * n + j] = j;
 165     p[i * n + i] = INF;
 166   }
 167
 168   /* --- Do the main algorithm --- *
 169    *
 170    * Not so hard.  Just brute force and ignorance.
 171    */
 172
 173   for (k = 0; k < n; k++) {
 174     for (i = 0; i < n; i++) {
 175       for (j = 0; j < n; j++) {
 176         if (a[i * n + k] != INF && a[k * n + j] != INF &&
 177             a[i * n + k] + a[k * n + j] < a[i * n + j]) {
 178           a[i * n + j] = a[i * n + k] + a[k * n + j];
 179           p[i * n + j] = p[i * n + k];
 180         }
 181       }
 182     }
 183   }
 184
 185   /* --- Wrap up --- */
 186
 187   if ((lv = export(ti, a, n)) == 0 || (pv = export(ti, p, n)) == 0)
 188     goto fail;
 189   o = Tcl_NewListObj(0, 0);
 190   Tcl_ListObjAppendElement
 191     (ti, o, Tcl_NewStringObj(Tcl_GetCommandName(ti, lv->c), -1));
 192   Tcl_ListObjAppendElement
 193     (ti, o, Tcl_NewStringObj(Tcl_GetCommandName(ti, pv->c), -1));
 194   Tcl_SetObjResult(ti, o);
 195   Tcl_Free((void *)a);
 196   Tcl_Free((void *)p);
 197   return (TCL_OK);
 198
 199 fail:
 200   if (a) Tcl_Free((void *)a);
 201   if (p) Tcl_Free((void *)p);
 202   if (lv) vec_destroy(ti, lv);
 203   if (pv) vec_destroy(ti, pv);
 204   return (TCL_ERROR);
 205 }
 206
 207 /*----- Travelling Salesman Problem ---------------------------------------*/
 208
 209 /* --- @rrange@ --- *
 210  *
 211  * Arguments:   @size_t max@ = maximum number wanted
 212  *
 213  * Returns:     An integer uniformly distributed on %$[0, max)$%.
 214  */
 215
 216 static size_t rrange(size_t max)
 217 {
 218   size_t m, z, r;
 219
 220   z = RAND_MAX/max;
 221   m = z * max;
 222   do {
 223     r = rand();
 224   } while (r > m);
 225   r /= z;
 226   return (r);
 227 }
 228
 229 /* --- @graph-travelling-salesman [-OPTIONS] ADJ LIST@ --- *
 230  *
 231  * Solves the Travelling Salesman Problem approximately.  Returns a list
 232  * containing (firstly) the cost of the computed route, and secondly the
 233  * route itself.  Only the nodes in LIST are considered.  The OPTIONS affect
 234  * the algorithm in various ways.
 235  *
 236  * -cool FACTOR         Cooling factor.  Default is 1.001.  Must be greater
 237  *                      than 1 for the simulated annealing to work.
 238  *
 239  * -dead COUNT          Give up after COUNT cycles with no improvement.
 240  *                      Default is 200.
 241  *
 242  * -inner COUNT         Perform COUNT loops each cooling cycle.  Default is
 243  *                      10000.
 244  *
 245  * -temp TEMP           Set the initial temperature to TEMP.  Default is not
 246  *                      very helpful.  Initial setting should be well above
 247  *                      the maximum cost increase from a cycle.
 248  *
 249  * -cycle / -nocycle    If -cycle is set, solve the classical problem of
 250  *                      finding a minimal cyclic path.  If -nocycle is set,
 251  *                      then start at the first node in LIST, and minimize a
 252  *                      tour without caring where the end goes.  The default
 253  *                      is -cycle.
 254  */
 255
 256 static int cmd_tsp(ClientData cd, Tcl_Interp *ti,
 257                    int objc, Tcl_Obj *const *objv)
 258 {
 259   /* --- Initial algorithm parameters --- */
 260
 261   double cool = 1.001;
 262   double temp = 1024;
 263   long inner = 10000;
 264   long dead = 200;
 265   int cycle = 1;
 266
 267   /* --- Other variables --- */
 268
 269   vec *v;
 270   unsigned long *a = 0;
 271   size_t n;
 272   int nn;
 273   size_t *r = 0, *r_best = 0;
 274   unsigned long c_best = 0, c_curr, c;
 275   size_t i, j, t;
 276   long ii, d;
 277   int ok;
 278   int rc = TCL_ERROR;
 279   Tcl_Obj *o, *o2, **oo;
 280
 281   /* --- Parse the command line --- */
 282
 283   for (i = 1; i < objc; i++) {
 284     int len;
 285     char *p = Tcl_GetStringFromObj(objv[i], &len);
 286     if (strcmp(p, "-cool") == 0) {
 287       i++; if (i >= objc) goto args;
 288       if (Tcl_GetDoubleFromObj(ti, objv[i], &cool) != TCL_OK)
 289         goto done;
 290       if (cool <= 1) {
 291         err(ti, "cooling factor must be > 1");
 292         goto done;
 293       }
 294     } else if (strcmp(p, "-temp") == 0) {
 295       i++; if (i >= objc) goto args;
 296       if (Tcl_GetDoubleFromObj(ti, objv[i], &temp) != TCL_OK)
 297         goto done;
 298       if (temp <= 0) {
 299         err(ti, "initial temperature must be > 0");
 300         goto done;
 301       }
 302     } else if (strcmp(p, "-inner") == 0) {
 303       i++; if (i >= objc) goto args;
 304       if (Tcl_GetLongFromObj(ti, objv[i], &inner) != TCL_OK)
 305         goto done;
 306       if (inner <= 0) {
 307         err(ti, "inner loop count must be > 0");
 308         goto done;
 309       }
 310     } else if (strcmp(p, "-dead") == 0) {
 311       i++; if (i >= objc) goto args;
 312       if (Tcl_GetLongFromObj(ti, objv[i], &dead) != TCL_OK)
 313         goto done;
 314       if (dead <= 0) {
 315         err(ti, "dead cycles count must be > 0");
 316         goto done;
 317       }
 318     } else if (strcmp(p, "-cycle") == 0)
 319       cycle = 1;
 320     else if (strcmp(p, "-nocycle") == 0)
 321       cycle = 0;
 322     else if (strcmp(p, "--") == 0) {
 323       i++; break;
 324     } else if (*p != '-')
 325       break;
 326     else {
 327       err(ti, "bad option for graph-travelling-salesman");
 328       goto done;
 329     }
 330   }
 331
 332   /* --- Check the rest --- */
 333
 334   if (i + 2 != objc) {
 335     err(ti, "usage: graph-travelling-salesman [-OPTIONS] ADJ LIST");
 336     goto done;
 337   }
 338   if ((v = vec_find(ti, objv[i])) == 0 || import(ti, v, &a, &n) != TCL_OK)
 339     goto done;
 340   if (Tcl_ListObjGetElements(ti, objv[i + 1], &nn, &oo) != TCL_OK)
 341     goto done;
 342   if (!nn)
 343     goto wrap;
 344
 345   r = (void *)Tcl_Alloc(nn * sizeof(*r));
 346   r_best = (void *)Tcl_Alloc(nn * sizeof(*r_best));
 347   for (i = 0; i < nn; i++) {
 348     long l;
 349     if (Tcl_GetLongFromObj(ti, oo[i], &l) != TCL_OK)
 350       goto done;
 351     if (l < 0 || l >= n) {
 352       err(ti, "node index out of range");
 353       goto done;
 354     }
 355     r[i] = l;
 356   }
 357
 358   /* --- The one and two node problems are trivial --- *
 359    *
 360    * Avoiding these prevents us from having to mess with special cases later.
 361    */
 362
 363   if (nn <= 2) {
 364     memcpy(r_best, r, nn * sizeof(*r));
 365     if (nn == 1)
 366       c_best = a[r[0] * n + r[0]];
 367     else
 368       c_best = a[r[0] * n + r[1]];
 369     goto wrap;
 370   }
 371
 372   /* --- Randomize the initial vector --- *
 373    *
 374    * If we're not cycling, then nail the first item in place.
 375    */
 376
 377   for (i = cycle ? 0 : 1; i < nn; i++) {
 378     j = rrange(nn - i);
 379     t = r[i]; r[i] = r[i + j]; r[i + j] = t;
 380   }
 381
 382   /* --- Compute the initial cost --- *
 383    *
 384    * If we're not cycling, don't close off at the end.  The easiest way to do
 385    * that is to start at the end.  There are at least three elements.
 386    */
 387
 388   if (cycle) { j = 0; i = nn - 1; }
 389   else { j = nn - 1; i = j - 1; }
 390   c = 0;
 391   for (;;) {
 392     c += a[r[i] * n + r[j]];
 393     if (!i)
 394       break;
 395     j = i;
 396     i--;
 397   }
 398
 399 /*   printf("*** initial cost = %lu; n = %u; nn = %u\n", c, n, nn); */
 400   c_curr = c_best = c;
 401   memcpy(r_best, r, nn * sizeof(*r));
 402
 403   /* --- Embark on the main loop --- */
 404
 405   d = dead;
 406   while (d) {
 407     ok = 0;
 408     for (ii = inner; ii; ii--) {
 409       size_t i, j, ilo, ihi, jlo, jhi;
 410
 411       /* --- Decide on a change to make --- *
 412        *
 413        * We just swap two nodes around on the path.  This is simple and seems
 414        * to be effective.  Don't allow the first node to be moved if we're
 415        * not cycling.
 416        */
 417
 418       if (cycle) {
 419         i = rrange(nn);
 420         j = rrange(nn);
 421       } else {
 422         i = rrange(nn - 1) + 1;
 423         j = rrange(nn - 1) + 1;
 424       }
 425
 426       /* --- Compute the change in cost --- *
 427        *
 428        * Since we're only swapping two nodes, we can work out the change
 429        * without rescanning the entire path, by just looking at the local
 430        * effects.
 431        */
 432
 433       if (i == j)
 434         continue; /* No change */
 435       if (j < i) { t = i; i = j; j = t; }
 436       ilo = (i + nn - 1) % nn; ihi = (i + 1) % nn;
 437       jlo = (j + nn - 1) % nn; jhi = (j + 1) % nn;
 438
 439       c = c_curr;
 440       if (j == nn - 1) {
 441
 442         /* --- This is where the algorithms differ --- *
 443          *
 444          * If we're producing a cycle, then we need the cost function to wrap
 445          * around here.  Otherwise, it hits a barrier, and the last node only
 446          * has a partial effect.
 447          */
 448
 449         if (cycle) {
 450           if (i == 0) {
 451             c -= (a[r[jlo] * n + r[j]] +
 452                   a[r[j] * n + r[i]] +
 453                   a[r[i] * n + r[ihi]]);
 454             c += (a[r[jlo] * n + r[i]] +
 455                   a[r[i] * n + r[j]] +
 456                   a[r[j] * n + r[ihi]]);
 457           } else goto std;
 458         } else {
 459           if (i == j - 1) {
 460             c -= a[r[ilo] * n + r[i]] + a[r[i] * n + r[j]];
 461             c += a[r[ilo] * n + r[j]] + a[r[j] * n + r[i]];
 462           } else {
 463             c -= (a[r[ilo] * n + r[i]] +
 464                   a[r[i] * n + r[ihi]] +
 465                   a[r[jlo] * n + r[j]]);
 466             c += (a[r[ilo] * n + r[j]] +
 467                   a[r[j] * n + r[ihi]] +
 468                   a[r[jlo] * n + r[i]]);
 469           }
 470         }
 471       } else {
 472
 473         /* --- Usual case --- *
 474          *
 475          * This splits into two subcases, depending on whether the areas
 476          * overlap.
 477          */
 478
 479       std:
 480         if (i == j - 1) {
 481           c -= (a[r[ilo] * n + r[i]] +
 482                 a[r[i] * n + r[j]] +
 483                 a[r[j] * n + r[jhi]]);
 484           c += (a[r[ilo] * n + r[j]] +
 485                 a[r[j] * n + r[i]] +
 486                 a[r[i] * n + r[jhi]]);
 487         } else {
 488           c -= (a[r[ilo] * n + r[i]] +
 489                 a[r[i] * n + r[ihi]] +
 490                 a[r[jlo] * n + r[j]] +
 491                 a[r[j] * n + r[jhi]]);
 492           c += (a[r[ilo] * n + r[j]] +
 493                 a[r[j] * n + r[ihi]] +
 494                 a[r[jlo] * n + r[i]] +
 495                 a[r[i] * n + r[jhi]]);
 496         }
 497       }
 498
 499 #ifdef PARANOID_CHECKING /* Turn this on to check the shortcut */
 500       {
 501         unsigned long cc;
 502         size_t ii, jj;
 503         if (cycle) { jj = 0; ii = nn - 1; }
 504         else { jj = nn - 1; ii = jj - 1; }
 505         cc = 0;
 506         t = r[i]; r[i] = r[j]; r[j] = t;
 507         for (;;) {
 508           cc += a[r[ii] * n + r[jj]];
 509           if (!ii)
 510             break;
 511           jj = ii;
 512           ii--;
 513         }
 514         t = r[i]; r[i] = r[j]; r[j] = t;
 515         if (c != cc) {
 516           printf("i = %u; j = %u; c = %lu; cc = %lu\n", i, j, c, cc);
 517           abort();
 518         }
 519       }
 520 #endif
 521
 522       /* --- Decide what to do --- */
 523
 524       if (c > c_curr &&
 525           rrange(65536) >= (size_t)(exp(((double)c_curr -
 526                                          (double)c)/temp) * 65536))
 527         continue;
 528
 529       /* --- Accept the change --- */
 530
 531       if (c < c_curr)
 532         ok = 1;
 533       c_curr = c;
 534       t = r[i]; r[i] = r[j]; r[j] = t;
 535       if (c_curr < c_best) {
 536         c_best = c_curr;
 537 /*      printf("*** new best = %lu\n", c_best); */
 538         memcpy(r_best, r, nn * sizeof(*r));
 539       }
 540     }
 541     temp /= cool;
 542     if (ok)
 543       d = dead;
 544     else
 545       d--;
 546   }
 547
 548   /* --- Done --- */
 549
 550 wrap:
 551   o = Tcl_NewListObj(0, 0);
 552   o2 = Tcl_NewListObj(0, 0);
 553   Tcl_ListObjAppendElement(ti, o, Tcl_NewLongObj(c_best));
 554   for (i = 0; i < nn; i++)
 555     Tcl_ListObjAppendElement(ti, o2, Tcl_NewLongObj(r_best[i]));
 556   Tcl_ListObjAppendElement(ti, o, o2);
 557   Tcl_SetObjResult(ti, o);
 558   rc = TCL_OK;
 559
 560   /* --- Tidy up --- */
 561
 562 done:
 563   if (a) Tcl_Free((void *)a);
 564   if (r) Tcl_Free((void *)r);
 565   if (r_best) Tcl_Free((void *)r_best);
 566   return (rc);
 567
 568 args:
 569   err(ti, "missing argument for option");
 570   goto done;
 571 }
 572
 573 /*----- Initialization ----------------------------------------------------*/
 574
 575 int Graph_SafeInit(Tcl_Interp *ti)
 576 {
 577   static const struct cmd {
 578     /*const*/ char *name;
 579     Tcl_ObjCmdProc *proc;
 580   } cmds[] = {
 581     { "graph-shortest-path",    cmd_shortestpath },
 582     { "graph-travelling-salesman", cmd_tsp },
 583     { 0,                        0 }
 584   };
 585
 586   const struct cmd *c;
 587   if (Tcl_PkgRequire(ti, "vector", "1.0.0", 0) == 0)
 588     return (TCL_ERROR);
 589   for (c = cmds; c->name; c++)
 590     Tcl_CreateObjCommand(ti, c->name, c->proc, 0, 0);
 591   if (Tcl_PkgProvide(ti, "graph", "1.0.0"))
 592     return (TCL_ERROR);
 593   return (TCL_OK);
 594 }
 595
 596 int Graph_Init(Tcl_Interp *ti)
 597 {
 598   return (Graph_SafeInit(ti));
 599 }
 600
 601 /*----- That's all, folks -------------------------------------------------*/