chiark / gitweb /
And then screwed the code over...
[rocl] / 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 -------------------------------------------------*/