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

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