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