From 6f7ff7517fc2b5d67138ddd1055020385bf6baa7 Mon Sep 17 00:00:00 2001 From: mdw Date: Fri, 7 Mar 2003 00:45:13 +0000 Subject: [PATCH] Graph theory functions. --- graph.c | 588 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 588 insertions(+) create mode 100644 graph.c diff --git a/graph.c b/graph.c new file mode 100644 index 0000000..0f9fa84 --- /dev/null +++ b/graph.c @@ -0,0 +1,588 @@ +/* -*-c-*- + * + * $Id: graph.c,v 1.1 2003/03/07 00:45:13 mdw Exp $ + * + * Graph theory stuff + * + * (c) 2003 Mark Wooding + */ + +/*----- Licensing notice --------------------------------------------------* + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ + +/*----- Revision history --------------------------------------------------* + * + * $Log: graph.c,v $ + * Revision 1.1 2003/03/07 00:45:13 mdw + * Graph theory functions. + * + */ + +/*----- Header files ------------------------------------------------------*/ + +#include +#include +#include +#include +#include + +#include + +#include "vec.h" + +/*----- Static variables --------------------------------------------------*/ + +#define INF ((unsigned long)-1) + +/*----- Utility functions -------------------------------------------------*/ + +static int err(Tcl_Interp *ti, /*const*/ char *p) +{ + Tcl_SetResult(ti, p, TCL_STATIC); + return (TCL_ERROR); +} + +/* --- @import@ --- * + * + * Arguments: @Tcl_Interp *ti@ = interpreter to leave errors in + * @vec *v@ = pointer to input adjacency matrix + * @unsigned long *tt@ = pointer to output adjacency matrix + * @size_t *nn@ = where to put the table size + * + * Returns: Tcl return code. + * + * Use: Imports an adjacency matrix. + */ + +static int import(Tcl_Interp *ti, vec *v, unsigned long **tt, size_t *nn) +{ + size_t i; + unsigned long *t; + size_t n; + + /* --- Check the table is well-formed --- */ + + if (v->ndim != 2) + return (err(ti, "adjacency matrix must be two-dimensional")); + if (v->dim[0].lo != 0 || v->dim[1].lo || v->dim[0].hi != v->dim[1].hi) + return (err(ti, "adjacency matrix must be square and zero-origin")); + n = *nn = v->dim[0].hi; + + /* --- Copy the data over --- */ + + n *= n; + assert(n == v->n); + t = (void *)Tcl_Alloc(n * sizeof(*t)); + for (i = 0; i < n; i++) { + long l; + if (Tcl_GetLongFromObj(ti, v->v[i], &l) != TCL_OK) { + Tcl_Free((void *)t); + return (TCL_ERROR); + } + t[i] = l >= 0 ? l : INF; + } + *tt = t; + return (TCL_OK); +} + +/* --- @export@ --- * + * + * Arguments: @Tcl_Interp *ti@ = interpreter to create output vector + * @unsigned long *t@ = pointer to table + * @size_t n@ = size of the table + * + * Returns: A pointer to the vector, or null. + * + * Use: Exports an adjacency matrix. + */ + +static vec *export(Tcl_Interp *ti, unsigned long *t, size_t n) +{ + vec_bound b[2]; + vec *v; + size_t i; + Tcl_Obj *o; + + b[0].lo = b[1].lo = 0; + b[0].hi = b[1].hi = n; + if ((v = vec_create(ti, 2, b, 0)) == 0) + return (0); + o = Tcl_NewLongObj(-1); + Tcl_IncrRefCount(o); + for (i = 0; i < v->n; i++) { + v->v[i] = t[i] == INF ? o : Tcl_NewLongObj(t[i]); + Tcl_IncrRefCount(v->v[i]); + } + Tcl_DecrRefCount(o); + return (v); +} + +/*----- Floyd-Warshall all-points shortest path ---------------------------*/ + +/* --- @graph-shortest-path VEC@ --- * + * + * Returns a pair of vectors containing, respectively, the shortest path + * length and the successor element in the shortest path. If you say + * + * destructure {len path} [graph-shortest-path $v] + * + * then [$len get I J] is the shortest path length from node I to node J, and + * [$path get I J] is the first hop on that shortest path. (To compute the + * entire path, set K to be that first hop; the next hop is then [$path get K + * J], and so on.) + * + * The adjacency matrix is given in VEC: negative entries indicate no path; + * nonnegative entries are weights. All entries must be integers. + */ + +static int cmd_shortestpath(ClientData cd, Tcl_Interp *ti, + int objc, Tcl_Obj *const *objv) +{ + vec *v, *lv = 0, *pv = 0; + size_t n, i, j, k; + unsigned long *a = 0, *p = 0; + Tcl_Obj *o; + + /* --- Read in the arguments --- */ + + if (objc != 2) { + err(ti, "usage: graph-shortest-path VEC"); + goto fail; + } + if ((v = vec_find(ti, objv[1])) == 0 || import(ti, v, &a, &n) != TCL_OK) + goto fail; + + /* --- Set up the path table --- */ + + p = (void *)Tcl_Alloc(n * n * sizeof(*p)); + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) + p[i * n + j] = j; + p[i * n + i] = INF; + } + + /* --- Do the main algorithm --- * + * + * Not so hard. Just brute force and ignorance. + */ + + for (k = 0; k < n; k++) { + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + if (a[i * n + k] != INF && a[k * n + j] != INF && + a[i * n + k] + a[k * n + j] < a[i * n + j]) { + a[i * n + j] = a[i * n + k] + a[k * n + j]; + p[i * n + j] = p[i * n + k]; + } + } + } + } + + /* --- Wrap up --- */ + + if ((lv = export(ti, a, n)) == 0 || (pv = export(ti, p, n)) == 0) + goto fail; + o = Tcl_NewListObj(0, 0); + Tcl_ListObjAppendElement + (ti, o, Tcl_NewStringObj(Tcl_GetCommandName(ti, lv->c), -1)); + Tcl_ListObjAppendElement + (ti, o, Tcl_NewStringObj(Tcl_GetCommandName(ti, pv->c), -1)); + Tcl_SetObjResult(ti, o); + Tcl_Free((void *)a); + Tcl_Free((void *)p); + return (TCL_OK); + +fail: + if (a) Tcl_Free((void *)a); + if (p) Tcl_Free((void *)p); + if (lv) vec_destroy(ti, lv); + if (pv) vec_destroy(ti, pv); + return (TCL_ERROR); +} + +/*----- Travelling Salesman Problem ---------------------------------------*/ + +/* --- @rrange@ --- * + * + * Arguments: @size_t max@ = maximum number wanted + * + * Returns: An integer uniformly distributed on %$[0, max)$%. + */ + +static size_t rrange(size_t max) +{ + size_t m, z, r; + + z = RAND_MAX/max; + m = z * max; + do { + r = rand(); + } while (r > m); + r /= z; + return (r); +} + +/* --- @graph-travelling-salesman [-OPTIONS] ADJ LIST@ --- * + * + * Solves the Travelling Salesman Problem approximately. Returns a list + * containing (firstly) the cost of the computed route, and secondly the + * route itself. Only the nodes in LIST are considered. The OPTIONS affect + * the algorithm in various ways. + * + * -cool FACTOR Cooling factor. Default is 1.001. Must be greater + * than 1 for the simulated annealing to work. + * + * -dead COUNT Give up after COUNT cycles with no improvement. + * Default is 200. + * + * -inner COUNT Perform COUNT loops each cooling cycle. Default is + * 10000. + * + * -init TEMP Set the initial temperature to TEMP. Default is not + * very helpful. Initial setting should be well above + * the maximum cost increase from a cycle. + * + * -cycle / -nocycle If -cycle is set, solve the classical problem of + * finding a minimal cyclic path. If -nocycle is set, + * then start at the first node in LIST, and minimize a + * tour without caring where the end goes. The default + * is -cycle. + */ + +static int cmd_tsp(ClientData cd, Tcl_Interp *ti, + int objc, Tcl_Obj *const *objv) +{ + /* --- Initial algorithm parameters --- */ + + double cool = 1.001; + double temp = 1024; + long inner = 10000; + long dead = 200; + int cycle = 1; + + /* --- Other variables --- */ + + vec *v; + unsigned long *a = 0; + size_t n; + int nn; + size_t *r = 0, *r_best = 0; + unsigned long c_best = 0, c_curr, c; + size_t i, j, t; + long ii, d; + int ok; + int rc = TCL_ERROR; + Tcl_Obj *o, *o2, **oo; + + /* --- Parse the command line --- */ + + for (i = 1; i < objc; i++) { + int len; + char *p = Tcl_GetStringFromObj(objv[i], &len); + if (strcmp(p, "-cool") == 0) { + i++; if (i >= objc) goto args; + if (Tcl_GetDoubleFromObj(ti, objv[i], &cool) != TCL_OK) + goto done; + if (cool <= 1) { + err(ti, "cooling factor must be > 1"); + goto done; + } + } else if (strcmp(p, "-init") == 0) { + i++; if (i >= objc) goto args; + if (Tcl_GetDoubleFromObj(ti, objv[i], &temp) != TCL_OK) + goto done; + if (temp <= 0) { + err(ti, "initial temperature must be > 0"); + goto done; + } + } else if (strcmp(p, "-inner") == 0) { + i++; if (i >= objc) goto args; + if (Tcl_GetLongFromObj(ti, objv[i], &inner) != TCL_OK) + goto done; + if (inner <= 0) { + err(ti, "inner loop count must be > 0"); + goto done; + } + } else if (strcmp(p, "-dead") == 0) { + i++; if (i >= objc) goto args; + if (Tcl_GetLongFromObj(ti, objv[i], &dead) != TCL_OK) + goto done; + if (dead <= 0) { + err(ti, "dead cycles count must be > 0"); + goto done; + } + } else if (strcmp(p, "-cycle") == 0) + cycle = 1; + else if (strcmp(p, "-nocycle") == 0) + cycle = 0; + else if (strcmp(p, "--") == 0) { + i++; break; + } else if (*p != '-') + break; + else { + err(ti, "bad option for graph-travelling-salesman"); + goto done; + } + } + + /* --- Check the rest --- */ + + if (i + 2 != objc) { + err(ti, "usage: graph-travelling-salesman [-OPTIONS] ADJ LIST"); + goto done; + } + if ((v = vec_find(ti, objv[i])) == 0 || import(ti, v, &a, &n) != TCL_OK) + goto done; + if (Tcl_ListObjGetElements(ti, objv[i + 1], &nn, &oo) != TCL_OK) + goto done; + if (!nn) + goto wrap; + + r = (void *)Tcl_Alloc(nn * sizeof(*r)); + r_best = (void *)Tcl_Alloc(nn * sizeof(*r_best)); + for (i = 0; i < nn; i++) { + long l; + if (Tcl_GetLongFromObj(ti, oo[i], &l) != TCL_OK) + goto done; + if (l < 0 || l >= n) { + err(ti, "node index out of range"); + goto done; + } + r[i] = l; + } + + /* --- The one and two node problems are trivial --- * + * + * Avoiding these prevents us from having to mess with special cases later. + */ + + if (nn <= 2) { + memcpy(r_best, r, nn * sizeof(*r)); + if (n == 1) + c_best = a[r[0] * n + r[0]]; + else + c_best = a[r[0] * n + r[1]]; + goto wrap; + } + + /* --- Randomize the initial vector --- * + * + * If we're not cycling, then nail the first item in place. + */ + + for (i = cycle ? 0 : 1; i < nn; i++) { + j = rrange(nn - i); + t = r[i]; r[i] = r[i + j]; r[i + j] = t; + } + + /* --- Compute the initial cost --- * + * + * If we're not cycling, don't close off at the end. The easiest way to do + * that is to start at the end. There are at least three elements. + */ + + if (cycle) { j = 0; i = nn - 1; } + else { j = nn - 1; i = j - 1; } + c = 0; + for (;;) { + c += a[r[i] * n + r[j]]; + if (!i) + break; + j = i; + i--; + } + +/* printf("*** initial cost = %lu\n", c_best); */ + c_curr = c_best = c; + memcpy(r_best, r, nn * sizeof(*r)); + + /* --- Embark on the main loop --- */ + + d = dead; + while (d) { + ok = 0; + for (ii = inner; ii; ii--) { + size_t i, j, ilo, ihi, jlo, jhi; + + /* --- Decide on a change to make --- * + * + * We just swap two nodes around on the path. This is simple and seems + * to be effective. Don't allow the first node to be moved if we're + * not cycling. + */ + + if (cycle) { + i = rrange(nn); + j = rrange(nn); + } else { + i = rrange(nn - 1) + 1; + j = rrange(nn - 1) + 1; + } + + /* --- Compute the change in cost --- * + * + * Since we're only swapping two nodes, we can work out the change + * without rescanning the entire path, by just looking at the local + * effects. + */ + + if (i == j) + continue; /* No change */ + if (j < i) { t = i; i = j; j = t; } + ilo = (i - 1) % nn; ihi = (i + 1) % nn; + jlo = (j - 1) % nn; jhi = (j + 1) % nn; + + c = c_curr; + if (j == nn - 1) { + + /* --- This is where the algorithms differ --- * + * + * If we're producing a cycle, then we need the cost function to wrap + * around here. Otherwise, it hits a barrier, and the last node only + * has a partial effect. + */ + + if (cycle) { + if (i == 0) { + c -= (a[r[jlo] * n + r[j]] + + a[r[j] * n + r[i]] + + a[r[i] * n + r[ihi]]); + c += (a[r[jlo] * n + r[i]] + + a[r[i] * n + r[j]] + + a[r[j] * n + r[ihi]]); + } else goto std; + } else { + if (i == j - 1) { + c -= a[r[ilo] * n + r[i]] + a[r[i] * n + r[j]]; + c += a[r[ilo] * n + r[j]] + a[r[j] * n + r[i]]; + } else { + c -= (a[r[ilo] * n + r[i]] + + a[r[i] * n + r[ihi]] + + a[r[jlo] * n + r[j]]); + c += (a[r[ilo] * n + r[j]] + + a[r[j] * n + r[ihi]] + + a[r[jlo] * n + r[i]]); + } + } + } else { + + /* --- Usual case --- * + * + * This splits into two subcases, depending on whether the areas + * overlap. + */ + + std: + if (i == j - 1) { + c -= (a[r[ilo] * n + r[i]] + + a[r[i] * n + r[j]] + + a[r[j] * n + r[jhi]]); + c += (a[r[ilo] * n + r[j]] + + a[r[j] * n + r[i]] + + a[r[i] * n + r[jhi]]); + } else { + c -= (a[r[ilo] * n + r[i]] + + a[r[i] * n + r[ihi]] + + a[r[jlo] * n + r[j]] + + a[r[j] * n + r[jhi]]); + c += (a[r[ilo] * n + r[j]] + + a[r[j] * n + r[ihi]] + + a[r[jlo] * n + r[i]] + + a[r[i] * n + r[jhi]]); + } + } + + /* --- Decide what to do --- */ + + if (c > c_curr && + rrange(65536) >= (size_t)(exp(((double)c_curr - + (double)c)/temp) * 65536)) + continue; + + /* --- Accept the change --- */ + + if (c < c_curr) + ok = 1; + c_curr = c; + t = r[i]; r[i] = r[j]; r[j] = t; + if (c_curr < c_best) { + c_best = c_curr; +/* printf("*** new best = %lu\n", c_best); */ + memcpy(r_best, r, nn * sizeof(*r)); + } + } + temp /= cool; + if (ok) + d = dead; + else + d--; + } + + /* --- Done --- */ + +wrap: + o = Tcl_NewListObj(0, 0); + o2 = Tcl_NewListObj(0, 0); + Tcl_ListObjAppendElement(ti, o, Tcl_NewLongObj(c_best)); + for (i = 0; i < nn; i++) + Tcl_ListObjAppendElement(ti, o2, Tcl_NewLongObj(r_best[i])); + Tcl_ListObjAppendElement(ti, o, o2); + Tcl_SetObjResult(ti, o); + rc = TCL_OK; + + /* --- Tidy up --- */ + +done: + if (a) Tcl_Free((void *)a); + if (r) Tcl_Free((void *)r); + if (r_best) Tcl_Free((void *)r_best); + return (rc); + +args: + err(ti, "missing argument for option"); + goto done; +} + +/*----- Initialization ----------------------------------------------------*/ + +int Graph_SafeInit(Tcl_Interp *ti) +{ + static const struct cmd { + /*const*/ char *name; + Tcl_ObjCmdProc *proc; + } cmds[] = { + { "graph-shortest-path", cmd_shortestpath }, + { "graph-travelling-salesman", cmd_tsp }, + { 0, 0 } + }; + + const struct cmd *c; + if (Tcl_PkgRequire(ti, "vector", "1.0.0", 0) == 0) + return (TCL_ERROR); + for (c = cmds; c->name; c++) + Tcl_CreateObjCommand(ti, c->name, c->proc, 0, 0); + if (Tcl_PkgProvide(ti, "graph", "1.0.0")) + return (TCL_ERROR); + return (TCL_OK); +} + +int Graph_Init(Tcl_Interp *ti) +{ + return (Graph_SafeInit(ti)); +} + +/*----- That's all, folks -------------------------------------------------*/ -- 2.11.0