chiark / gitweb /
Fix output formatting a little.
[rocl] / graph.c
CommitLineData
6f7ff751 1/* -*-c-*-
2 *
c53e2dd3 3 * $Id$
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
6f7ff751 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
45static 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
63static 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
105static 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
144static 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
201fail:
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
218static 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 *
b758c343 247 * -temp TEMP Set the initial temperature to TEMP. Default is not
6f7ff751 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
258static 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 }
b758c343 296 } else if (strcmp(p, "-temp") == 0) {
6f7ff751 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));
a9dd7681 367 if (nn == 1)
6f7ff751 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
b758c343 401/* printf("*** initial cost = %lu; n = %u; nn = %u\n", c, n, nn); */
6f7ff751 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; }
b758c343 438 ilo = (i + nn - 1) % nn; ihi = (i + 1) % nn;
439 jlo = (j + nn - 1) % nn; jhi = (j + 1) % nn;
6f7ff751 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
b758c343 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
6f7ff751 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
552wrap:
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
564done:
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
570args:
571 err(ti, "missing argument for option");
572 goto done;
573}
574
575/*----- Initialization ----------------------------------------------------*/
576
577int 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
598int Graph_Init(Tcl_Interp *ti)
599{
600 return (Graph_SafeInit(ti));
601}
602
603/*----- That's all, folks -------------------------------------------------*/