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

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