chiark / gitweb /
Print seed -- it's useful.
[rocl] / graph.c
CommitLineData
6f7ff751 1/* -*-c-*-
2 *
a9dd7681 3 * $Id: graph.c,v 1.3 2003/03/10 23:37:21 mdw Exp $
6f7ff751 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 $
a9dd7681 30 * Revision 1.3 2003/03/10 23:37:21 mdw
31 * Fix single point TSP.
32 *
b758c343 33 * Revision 1.2 2003/03/08 00:40:32 mdw
34 * Fix unsigned crapness in travelling-salesman solver.
35 *
6f7ff751 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
59static 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
77static 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
119static 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
158static 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
215fail:
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
232static 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 *
b758c343 261 * -temp TEMP Set the initial temperature to TEMP. Default is not
6f7ff751 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
272static 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 }
b758c343 310 } else if (strcmp(p, "-temp") == 0) {
6f7ff751 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));
a9dd7681 381 if (nn == 1)
6f7ff751 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
b758c343 415/* printf("*** initial cost = %lu; n = %u; nn = %u\n", c, n, nn); */
6f7ff751 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; }
b758c343 452 ilo = (i + nn - 1) % nn; ihi = (i + 1) % nn;
453 jlo = (j + nn - 1) % nn; jhi = (j + 1) % nn;
6f7ff751 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
b758c343 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
6f7ff751 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
566wrap:
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
578done:
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
584args:
585 err(ti, "missing argument for option");
586 goto done;
587}
588
589/*----- Initialization ----------------------------------------------------*/
590
591int 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
612int Graph_Init(Tcl_Interp *ti)
613{
614 return (Graph_SafeInit(ti));
615}
616
617/*----- That's all, folks -------------------------------------------------*/