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 | |
45 | static 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 | |
63 | static 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 | |
105 | static 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 | |
144 | static 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 | |
201 | fail: |
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 | |
218 | static 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 | |
258 | static 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 | |
552 | wrap: |
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 | |
564 | done: |
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 | |
570 | args: |
571 | err(ti, "missing argument for option"); |
572 | goto done; |
573 | } |
574 | |
575 | /*----- Initialization ----------------------------------------------------*/ |
576 | |
577 | int 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 | |
598 | int Graph_Init(Tcl_Interp *ti) |
599 | { |
600 | return (Graph_SafeInit(ti)); |
601 | } |
602 | |
603 | /*----- That's all, folks -------------------------------------------------*/ |