6f7ff751 |
1 | /* -*-c-*- |
2 | * |
6f7ff751 |
3 | * Graph theory stuff |
4 | * |
5 | * (c) 2003 Mark Wooding |
6 | */ |
7 | |
378b623c |
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. |
378b623c |
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. |
378b623c |
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 | |
43 | static 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 | |
61 | static 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 | |
103 | static 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 | |
142 | static 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 | |
199 | fail: |
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 | |
216 | static 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 | |
256 | static 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; |
378b623c |
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 | |
550 | wrap: |
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 | |
562 | done: |
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 | |
568 | args: |
569 | err(ti, "missing argument for option"); |
570 | goto done; |
571 | } |
572 | |
573 | /*----- Initialization ----------------------------------------------------*/ |
574 | |
575 | int 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 | |
596 | int Graph_Init(Tcl_Interp *ti) |
597 | { |
598 | return (Graph_SafeInit(ti)); |
599 | } |
600 | |
601 | /*----- That's all, folks -------------------------------------------------*/ |