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