+/* -*-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 <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <tcl.h>
+
+#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 -------------------------------------------------*/
~mdw / rocl / commitdiff