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