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

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