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