3 * $Id: graph.c,v 1.2 2003/03/08 00:40:32 mdw Exp $
7 * (c) 2003 Mark Wooding
10 /*----- Licensing notice --------------------------------------------------*
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.
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.
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.
27 /*----- Revision history --------------------------------------------------*
30 * Revision 1.2 2003/03/08 00:40:32 mdw
31 * Fix unsigned crapness in travelling-salesman solver.
33 * Revision 1.1 2003/03/07 00:45:13 mdw
34 * Graph theory functions.
38 /*----- Header files ------------------------------------------------------*/
50 /*----- Static variables --------------------------------------------------*/
52 #define INF ((unsigned long)-1)
54 /*----- Utility functions -------------------------------------------------*/
56 static int err(Tcl_Interp
*ti
, /*const*/ char *p
)
58 Tcl_SetResult(ti
, p
, TCL_STATIC
);
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
69 * Returns: Tcl return code.
71 * Use: Imports an adjacency matrix.
74 static int import(Tcl_Interp
*ti
, vec
*v
, unsigned long **tt
, size_t *nn
)
80 /* --- Check the table is well-formed --- */
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
;
88 /* --- Copy the data over --- */
92 t
= (void *)Tcl_Alloc(n
* sizeof(*t
));
93 for (i
= 0; i
< n
; i
++) {
95 if (Tcl_GetLongFromObj(ti
, v
->v
[i
], &l
) != TCL_OK
) {
99 t
[i
] = l
>= 0 ? l
: INF
;
105 /* --- @export@ --- *
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
111 * Returns: A pointer to the vector, or null.
113 * Use: Exports an adjacency matrix.
116 static vec
*export(Tcl_Interp
*ti
, unsigned long *t
, size_t n
)
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)
127 o
= Tcl_NewLongObj(-1);
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
]);
137 /*----- Floyd-Warshall all-points shortest path ---------------------------*/
139 /* --- @graph-shortest-path VEC@ --- *
141 * Returns a pair of vectors containing, respectively, the shortest path
142 * length and the successor element in the shortest path. If you say
144 * destructure {len path} [graph-shortest-path $v]
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
151 * The adjacency matrix is given in VEC: negative entries indicate no path;
152 * nonnegative entries are weights. All entries must be integers.
155 static int cmd_shortestpath(ClientData cd
, Tcl_Interp
*ti
,
156 int objc
, Tcl_Obj
*const *objv
)
158 vec
*v
, *lv
= 0, *pv
= 0;
160 unsigned long *a
= 0, *p
= 0;
163 /* --- Read in the arguments --- */
166 err(ti
, "usage: graph-shortest-path VEC");
169 if ((v
= vec_find(ti
, objv
[1])) == 0 || import(ti
, v
, &a
, &n
) != TCL_OK
)
172 /* --- Set up the path table --- */
174 p
= (void *)Tcl_Alloc(n
* n
* sizeof(*p
));
175 for (i
= 0; i
< n
; i
++) {
176 for (j
= 0; j
< n
; j
++)
181 /* --- Do the main algorithm --- *
183 * Not so hard. Just brute force and ignorance.
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
];
198 /* --- Wrap up --- */
200 if ((lv
= export(ti
, a
, n
)) == 0 || (pv
= export(ti
, p
, n
)) == 0)
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
);
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
);
220 /*----- Travelling Salesman Problem ---------------------------------------*/
222 /* --- @rrange@ --- *
224 * Arguments: @size_t max@ = maximum number wanted
226 * Returns: An integer uniformly distributed on %$[0, max)$%.
229 static size_t rrange(size_t max
)
242 /* --- @graph-travelling-salesman [-OPTIONS] ADJ LIST@ --- *
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.
249 * -cool FACTOR Cooling factor. Default is 1.001. Must be greater
250 * than 1 for the simulated annealing to work.
252 * -dead COUNT Give up after COUNT cycles with no improvement.
255 * -inner COUNT Perform COUNT loops each cooling cycle. Default is
258 * -temp TEMP Set the initial temperature to TEMP. Default is not
259 * very helpful. Initial setting should be well above
260 * the maximum cost increase from a cycle.
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
269 static int cmd_tsp(ClientData cd
, Tcl_Interp
*ti
,
270 int objc
, Tcl_Obj
*const *objv
)
272 /* --- Initial algorithm parameters --- */
280 /* --- Other variables --- */
283 unsigned long *a
= 0;
286 size_t *r
= 0, *r_best
= 0;
287 unsigned long c_best
= 0, c_curr
, c
;
292 Tcl_Obj
*o
, *o2
, **oo
;
294 /* --- Parse the command line --- */
296 for (i
= 1; i
< objc
; i
++) {
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
)
304 err(ti
, "cooling factor must be > 1");
307 } else if (strcmp(p
, "-temp") == 0) {
308 i
++; if (i
>= objc
) goto args
;
309 if (Tcl_GetDoubleFromObj(ti
, objv
[i
], &temp
) != TCL_OK
)
312 err(ti
, "initial temperature must be > 0");
315 } else if (strcmp(p
, "-inner") == 0) {
316 i
++; if (i
>= objc
) goto args
;
317 if (Tcl_GetLongFromObj(ti
, objv
[i
], &inner
) != TCL_OK
)
320 err(ti
, "inner loop count must be > 0");
323 } else if (strcmp(p
, "-dead") == 0) {
324 i
++; if (i
>= objc
) goto args
;
325 if (Tcl_GetLongFromObj(ti
, objv
[i
], &dead
) != TCL_OK
)
328 err(ti
, "dead cycles count must be > 0");
331 } else if (strcmp(p
, "-cycle") == 0)
333 else if (strcmp(p
, "-nocycle") == 0)
335 else if (strcmp(p
, "--") == 0) {
337 } else if (*p
!= '-')
340 err(ti
, "bad option for graph-travelling-salesman");
345 /* --- Check the rest --- */
348 err(ti
, "usage: graph-travelling-salesman [-OPTIONS] ADJ LIST");
351 if ((v
= vec_find(ti
, objv
[i
])) == 0 || import(ti
, v
, &a
, &n
) != TCL_OK
)
353 if (Tcl_ListObjGetElements(ti
, objv
[i
+ 1], &nn
, &oo
) != TCL_OK
)
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
++) {
362 if (Tcl_GetLongFromObj(ti
, oo
[i
], &l
) != TCL_OK
)
364 if (l
< 0 || l
>= n
) {
365 err(ti
, "node index out of range");
371 /* --- The one and two node problems are trivial --- *
373 * Avoiding these prevents us from having to mess with special cases later.
377 memcpy(r_best
, r
, nn
* sizeof(*r
));
379 c_best
= a
[r
[0] * n
+ r
[0]];
381 c_best
= a
[r
[0] * n
+ r
[1]];
385 /* --- Randomize the initial vector --- *
387 * If we're not cycling, then nail the first item in place.
390 for (i
= cycle ?
0 : 1; i
< nn
; i
++) {
392 t
= r
[i
]; r
[i
] = r
[i
+ j
]; r
[i
+ j
] = t
;
395 /* --- Compute the initial cost --- *
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.
401 if (cycle
) { j
= 0; i
= nn
- 1; }
402 else { j
= nn
- 1; i
= j
- 1; }
405 c
+= a
[r
[i
] * n
+ r
[j
]];
412 /* printf("*** initial cost = %lu; n = %u; nn = %u\n", c, n, nn); */
414 memcpy(r_best
, r
, nn
* sizeof(*r
));
416 /* --- Embark on the main loop --- */
421 for (ii
= inner
; ii
; ii
--) {
422 size_t i
, j
, ilo
, ihi
, jlo
, jhi
;
424 /* --- Decide on a change to make --- *
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
435 i
= rrange(nn
- 1) + 1;
436 j
= rrange(nn
- 1) + 1;
439 /* --- Compute the change in cost --- *
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
447 continue; /* No change */
448 if (j
< i
) { t
= i
; i
= j
; j
= t
; }
449 ilo
= (i
+ nn
- 1) % nn
; ihi
= (i
+ 1) % nn
;
450 jlo
= (j
+ nn
- 1) % nn
; jhi
= (j
+ 1) % nn
;
455 /* --- This is where the algorithms differ --- *
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.
464 c
-= (a
[r
[jlo
] * n
+ r
[j
]] +
466 a
[r
[i
] * n
+ r
[ihi
]]);
467 c
+= (a
[r
[jlo
] * n
+ r
[i
]] +
469 a
[r
[j
] * n
+ r
[ihi
]]);
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
]];
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
]]);
486 /* --- Usual case --- *
488 * This splits into two subcases, depending on whether the areas
494 c
-= (a
[r
[ilo
] * n
+ r
[i
]] +
496 a
[r
[j
] * n
+ r
[jhi
]]);
497 c
+= (a
[r
[ilo
] * n
+ r
[j
]] +
499 a
[r
[i
] * n
+ r
[jhi
]]);
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
]]);
512 #ifdef PARANOID_CHECKING /* Turn this on to check the shortcut */
516 if (cycle
) { jj
= 0; ii
= nn
- 1; }
517 else { jj
= nn
- 1; ii
= jj
- 1; }
519 t
= r
[i
]; r
[i
] = r
[j
]; r
[j
] = t
;
521 cc
+= a
[r
[ii
] * n
+ r
[jj
]];
527 t
= r
[i
]; r
[i
] = r
[j
]; r
[j
] = t
;
529 printf("i = %u; j = %u; c = %lu; cc = %lu\n", i
, j
, c
, cc
);
535 /* --- Decide what to do --- */
538 rrange(65536) >= (size_t)(exp(((double)c_curr
-
539 (double)c
)/temp
) * 65536))
542 /* --- Accept the change --- */
547 t
= r
[i
]; r
[i
] = r
[j
]; r
[j
] = t
;
548 if (c_curr
< c_best
) {
550 /* printf("*** new best = %lu\n", c_best); */
551 memcpy(r_best
, r
, nn
* sizeof(*r
));
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
);
573 /* --- Tidy up --- */
576 if (a
) Tcl_Free((void *)a
);
577 if (r
) Tcl_Free((void *)r
);
578 if (r_best
) Tcl_Free((void *)r_best
);
582 err(ti
, "missing argument for option");
586 /*----- Initialization ----------------------------------------------------*/
588 int Graph_SafeInit(Tcl_Interp
*ti
)
590 static const struct cmd
{
591 /*const*/ char *name
;
592 Tcl_ObjCmdProc
*proc
;
594 { "graph-shortest-path", cmd_shortestpath
},
595 { "graph-travelling-salesman", cmd_tsp
},
600 if (Tcl_PkgRequire(ti
, "vector", "1.0.0", 0) == 0)
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"))
609 int Graph_Init(Tcl_Interp
*ti
)
611 return (Graph_SafeInit(ti
));
614 /*----- That's all, folks -------------------------------------------------*/