Release 1.1.6.

[rocl] / graph.c

/* -*-c-*-
 *
 * 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.
 */

/*----- 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.
 *
 * -temp 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, "-temp") == 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 (nn == 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 = %u; nn = %u\n", c, n, nn); */
  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 + nn - 1) % nn; ihi = (i + 1) % nn;
      jlo = (j + nn - 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]]);
	}
      }

#ifdef PARANOID_CHECKING /* Turn this on to check the shortcut */
      {
	unsigned long cc;
	size_t ii, jj;
	if (cycle) { jj = 0; ii = nn - 1; }
	else { jj = nn - 1; ii = jj - 1; }
	cc = 0;
	t = r[i]; r[i] = r[j]; r[j] = t;
	for (;;) {
	  cc += a[r[ii] * n + r[jj]];
	  if (!ii)
	    break;
	  jj = ii;
	  ii--;
	}
	t = r[i]; r[i] = r[j]; r[j] = t;
	if (c != cc) {
	  printf("i = %u; j = %u; c = %lu; cc = %lu\n", i, j, c, cc);
	  abort();
	}
      }
#endif

      /* --- 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 -------------------------------------------------*/