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

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