/* An algorithm for drawing general undirected graphs
 * by Jonathan Ellis, based on an article of the same title by
 * Tomihisa Kamada and Satoru Kawai, Information Processing Letters, 4/12/89
 */

/* Tested on MS-DOS, SunOS systems */

/* usage: printgraph filename [closeEnough] */
/* Reads from specified file a graph and prints to stdout 
   a postscript layout for it.

   The input should have the form

   n
   v1 v2 w1
   v3 v4 w2
   v5 v6 w3
   etc

   where n is the number of vertices in the graph, v1, v2, ..., vn are
   vertices, and w1...wm are edges, where m is the number of edges in the
   graph.  (Graph is treated as unweighted if w1 = w2 = ... = wm = 1.)

   The results are undefined if you give it a flawed graph.  (An unconnected
   graph or a graph where you specify an incorrect number of vertices are
   examples of what the program considers flawed.  A graph with multiple
   edges between the same vertices is not flawed.)
*/

/* obvious problem--computes a local min.  Local mins aren't always
   the same or close to global mins though.  Example: a "square"
   4
   1 2 1
   1 3 1
   2 4 1
   3 4 1
   because of distributing points clockwise around a circle,
   tends to an "hourglass" instead of square.

   4
   1 2 1
   1 4 1
   2 3 1
   3 4 1
   produces the expected square, simply because of the initial positions.
*/

/* changing spring constant from ~ dist^-2 to ~ dist^-3 may help w/above. */

/* when moving points, move both ends */

/* give connections strength s.t. if graph *were* a circle, it would not
 * have to be changed from initial position */

/* need to do x reps on each point before recomputing maxDelta.  use prof. */

/* graph generator.  start with graph connected in circle, add random number
 * of connections to each vertex. */

/* numVertices ^2 memory used in computeDistances heap sucks.  don't think
 * there is better way that doesn't cost speed, but trading speed for
 * memory may be necessary at some point */

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

/* #define DEBUG 1 */

#define TOLERANCE 0.001
#define fzero(a) (fabs(a) < TOLERANCE)

/* #ifdef TURBO_C */
  #define PI 3.14159
/* #endif */
#define RADIUS 3.0
#define DRAW_WIDTH 400
#define SPRING_CONST 20

typedef struct {
  int distance, he, me; /* lower-distance nodes get popped first */
} heapNode;

typedef struct {
  int numNodes;
  int maxNodes;
  heapNode *data;
} heap;

/* distance = weight along edge from vertex he to vertex me */

typedef struct {
  int start, end; /* vertices that define the edge */
} answerNode;

typedef struct {
  int weighted, numVertices;
  int *weights, *drawn;
  float *x, *y, *distances, *springPower;
} graph;

int numTimes = 3; /* let user specify this, and find best value */

void newHeap(heap *h, int size);

void deleteHeap(heap h);

void upHeap(heap h, int k);

void downHeap(heap h, int k);

void insert(heap *h, heapNode v);

heapNode removeNode(heap *h);

void initGraph(graph *newGraph, int i);

void readGraph(FILE *fp, graph *g);

void destroyGraph(graph g);

void computeDistances(graph g);
/* determine shortest path length between each pair of vertices */

void computeDisplayDistances(graph g);
/* determine ideal length of springs connecting vertices */

void computeSpringPower(graph g);
/* call before computeDisplayDistances, which overwrites path lengths */

float dEdXm(graph g, int m); /* equation 7 */

float dEdYm(graph g, int m); /* t optimizable */ /* equation 8 */

float d2EdXm2(graph g, int m); /* equation 13 */

float d2EdXmdYm(graph g, int m); /* equation 14 */

float d2EdYmdXm(graph g, int m); /* equation 15 */

float d2EdYm2(graph g, int m); /* equation 16 */

void tuneGraph(graph g, float enough);
/* main loop */

void drawPair(graph g, int i, int j);
/* called by printGraph */

void printGraph(graph g);
/* print to stdout a postscript representation of graph */

void main(int argc, char *argv[])
{
  int i, j;
  float closeEnough = (float)DRAW_WIDTH / 100.0;
  graph myGraph;
  FILE *fp;

  if (argc < 2 || argc > 3)
  {
    printf("usage: printgraph filename [closeEnough]\n");
    exit(1);
  }
  if (argc == 3)
    closeEnough = atof(argv[2]);
  fp = fopen(argv[1], "r");
  readGraph(fp, &myGraph);
  computeDistances(myGraph);
  computeSpringPower(myGraph);
  computeDisplayDistances(myGraph);
  tuneGraph(myGraph, closeEnough);
  printGraph(myGraph);
  destroyGraph(myGraph);
}

void initGraph(graph *newGraph, int i)
{
  int j;
  newGraph->numVertices = i;
  i *= i;
  newGraph->weights = malloc(i * sizeof(int));
  newGraph->drawn = malloc(newGraph->numVertices * sizeof(int));
  newGraph->distances = malloc(i * sizeof(float));
  newGraph->springPower = malloc(i * sizeof(float));
  newGraph->x = malloc(newGraph->numVertices * sizeof(float));
  newGraph->y = malloc(newGraph->numVertices * sizeof(float));
  assert(newGraph->weights);
  assert(newGraph->drawn);
  assert(newGraph->distances);
  assert(newGraph->springPower);
  assert(newGraph->x);
  assert(newGraph->y);
  newGraph->weighted = 0;

  for (i = 0; i < newGraph->numVertices; i++)
  {
    int dist;
    float angle;
    for (j = 0; j < newGraph->numVertices; j++)
    {
      newGraph->weights[i * newGraph->numVertices + j] = 0;
      newGraph->distances[i * newGraph->numVertices + j] = 0;
    }
    /* initial positions are spaced evenly around the perimeter of a circle */
    dist = DRAW_WIDTH / 2;
    angle = 2 * PI * (float)i / newGraph->numVertices;
    newGraph->x[i] = dist + dist * cos(angle);
    newGraph->y[i] = dist + dist * sin(angle);
  }
}

void destroyGraph(graph g)
{
  free(g.weights);
  free(g.drawn);
  free(g.distances);
  free(g.springPower);
  free(g.x);
  free(g.y);
}

void readGraph(FILE *fp, graph *newGraph)
{
  int i, j, l;
  fscanf(fp, "%d", &i);
  initGraph(newGraph, i);
#ifdef DEBUG
  printf("DEBUG: finished initGraph\n");
#endif
  while (fscanf(fp, "%d %d %d", &j, &i, &l) != EOF)
  {
    j--;  /* internal representation numbers vertices 0..n-1 */
    i--;
    if (l > 1)
      newGraph->weighted = 1;
    newGraph->weights[i * newGraph->numVertices + j] = l;
    newGraph->weights[j * newGraph->numVertices + i] = l; 
    /* symmetric, since undirected */
  }
}

void computeDistances(graph g)
/* could half computation by making use of symmetry */
{
  int numAdded, i, j;
  float *finished;
  heap distHeap;
  heapNode t, u;

  for (j = 0; j < g.numVertices; j++)
  {
    newHeap(&distHeap, g.numVertices * g.numVertices);
    numAdded = 1;
    /* set up for search: put vertex j in finished and its neighbors in queue */

    finished = g.distances + j * g.numVertices;
    finished[j] = 1;
    for (i = 0; i < g.numVertices; i++)
      if (g.weights[j * g.numVertices + i])
        {
        t.distance = g.weights[j * g.numVertices + i];
        t.he = j;
        t.me = i;
        insert(&distHeap, t);
      }

    /* main loop.  add the "nearest" vertex to the finished set, and add
       vertices adjacent to it to the queue. */
    while (numAdded < g.numVertices)
      {
      /* reject any duplicate vertices */
      for (u.me = j; finished[u.me]; u = removeNode(&distHeap)) ;

      finished[u.me] = u.distance;
      for (i = 0; i < g.numVertices; i++)
        if (g.weights[u.me * g.numVertices + i] && !finished[i])
          {
              t.distance = u.distance + g.weights[u.me * g.numVertices + i];
              t.he = u.me;
              t.me = i;
              insert(&distHeap, t);
	  }
      numAdded++;
    }
    finished[j] = 0; /* we don't want a spring from a node to itself */
    deleteHeap(distHeap);
  }
#ifdef DEBUG
  for (i = 0; i < g.numVertices; i++)
    {
        printf("\nDEBUG: distances: ");
        for (j = 0; j < g.numVertices; j++)
      printf("%f ", g.distances[i * g.numVertices + j]);
      }
#endif
}

void computeDisplayDistances(graph g)
{
  float edgeLength;
  float max = 0;
  int i, j;

  for (i = 0; i < g.numVertices; i++)
    for (j = 0; j < g.numVertices; j++)
      if (g.distances[i * g.numVertices + j] > max)
        max = g.distances[i * g.numVertices + j];

  edgeLength = (float)DRAW_WIDTH / max;

  for (i = 0; i < g.numVertices; i++)
    for (j = 0; j < g.numVertices; j++)
      g.distances[i * g.numVertices + j] *= edgeLength;
#ifdef DEBUG
  for (i = 0; i < g.numVertices; i++)
  {
    printf("\nDEBUG: ddistances: ");
    for (j = 0; j < g.numVertices; j++)
      printf("%f ", g.distances[i * g.numVertices + j]);
  }
  printf("\n"); /* flush linebuffering */
#endif
}

void computeSpringPower(graph g)
{
  int i, j;
  float k;
  for (i = 0; i < g.numVertices; i++)
    for (j = 0; j < g.numVertices; j++)
      if (i != j)
      {
        k = g.distances[i * g.numVertices + j];
        g.springPower[i * g.numVertices + j] = (float)SPRING_CONST / (k * k);
      }
}

float dEdXm(graph g, int m) /* equation 7 */
{
  float sum = 0;
  int i;
  int t = m; /* i * g.numVertices + m */
  for (i = 0; i < g.numVertices; i++)
  {
    if (i != m)
    {
      float a = g.x[m] - g.x[i];
      float b = g.y[m] - g.y[i];
      sum += g.springPower[t] * (a - g.distances[t] * a /
                                 sqrt(a * a + b * b));
    }
    t += g.numVertices;
  }
  return sum;
}

float dEdYm(graph g, int m) /* t optimizable */ /* equation 8 */
{
  float sum = 0;
  int i;
  int t = m; /* i * g.numVertices + m */
  for (i = 0; i < g.numVertices; i++)
  {
    if (i != m)
    {
      float a = g.x[m] - g.x[i];
      float b = g.y[m] - g.y[i];
      sum += g.springPower[t] * (b - g.distances[t] * b /
                                 sqrt(a * a + b * b));
    }
    t += g.numVertices;
  }
  return sum;
}

float d2EdXm2(graph g, int m) /* equation 13 */
{
  float sum = 0;
  int i;
  int t = m; /* t = i * numVertices + m */
  for (i = 0; i < g.numVertices; i++)
  {
    if (i != m)
    {
      float a = g.x[m] - g.x[i];
      float b = g.y[m] - g.y[i];
      sum += g.springPower[t] * (1 - g.distances[t] * b * b /
                                 pow(a * a + b * b, 1.5));
    }
    t += g.numVertices;
  }
  return sum;
}

float d2EdXmdYm(graph g, int m) /* equation 14 */
{
  float sum = 0;
  int i;
  int t = m; /* i * g.numVertices + m */
  for (i = 0; i < g.numVertices; i++)
  {
    if (i != m)
    {
      float a = g.x[m] - g.x[i];
      float b = g.y[m] - g.y[i];
      sum += g.springPower[t] * (g.distances[t] * a * b /
                                 pow(a * a + b * b, 1.5));
    }
    t += g.numVertices;
  }
  return sum;
}

float d2EdYmdXm(graph g, int m) /* equation 15 */
{
  return d2EdXmdYm(g, m);
}

float d2EdYm2(graph g, int m) /* equation 16 */
{
  float sum = 0;
  int i;
  int t = m;
  for (i = 0; i < g.numVertices; i++)
  {
    if (i != m)
    {
      float a = g.x[m] - g.x[i];
      float b = g.y[m] - g.y[i];
      sum += g.springPower[t] * (1 - g.distances[t] * a * a /
                                 pow(a * a + b * b, 1.5));
    }
    t += g.numVertices;
  }
  return sum;
}

#define DONE -1

int maxDelta(graph g, float enough) 
/* find which particle has largest delta; eq'n 9 */
{
  float delta = 0;
  int which = 0, i;
  for (i = 0; i < g.numVertices; i++)
    {
    float a, b, d;
    a = dEdXm(g, i);
    b = dEdYm(g, i);
    d = sqrt(a * a + b * b);
    if (d > delta)
      {
      delta = d;
      which = i;
    }
  }
  if (delta > enough)
    return which;
  else
      return DONE;
}

void tuneGraph(graph g, float enough)
/* paper picks one point, and stays with it until good enough, then next pt.
   we pick one point, adjust, pick the new worst, adjust, until done. */
{
  int m, i;
  float sigmaY, sigmaX;
  while ((m = maxDelta(g, enough)) != DONE)
    for (i = 0; i < numTimes; i++)
    {
    float a, b, c, d, e, f, h; /* coefficients in eq'ns 11, 12 */
    a = d2EdXm2(g, m);
    b = d2EdXmdYm(g, m);
    c = d2EdYmdXm(g, m);
    d = d2EdYm2(g, m);
    e = -dEdXm(g, m);
    f = -dEdYm(g, m);
    if (fzero(a)) /* first pivot close to zero => exchange rows */
      {
      h = a;
      a = c;
      c = h;
      h = b;
      b = d;
      d = h;
      h = e;
      e = f;
      f = h;
    }
    /* gaussian elimination */
    h = c / a;
    d -= b * h;
    f -= e * h;
    sigmaY = f / d;
    sigmaX = (e - b * sigmaY) / a;
    g.x[m] += sigmaX;
    g.y[m] += sigmaY;
  }
}

void newHeap(heap *h, int size)
{
  h->maxNodes = size;
  h->numNodes = 0;
  h->data = malloc((size + 1) * sizeof(heapNode)); /* leave room for sentinel */

  assert(h->data);
}

void deleteHeap(heap h)
{
  free(h.data);
}

void upHeap(heap h, int k)
{
  int v;
  heapNode t;
  v = h.data[k].distance;
  t = h.data[k];
  h.data[0].distance = -1;
  while (h.data[k / 2].distance >= v) {
    h.data[k] = h.data[k / 2];
    k /= 2;
  }
  h.data[k] = t;
}

void downHeap(heap h, int k)
{
  int j;
  heapNode t;
  t = h.data[k];
  while (k <= h.numNodes / 2) {
    j = k + k;
    if (j < h.numNodes && h.data[j].distance > h.data[j + 1].distance) j++;
    if (t.distance <= h.data[j].distance) break;
    h.data[k] = h.data[j];
    k = j;
  }
  h.data[k] = t;
}

void insert(heap *h, heapNode t)
{
  h->numNodes++;
  assert(h->numNodes <= h->maxNodes);
  h->data[h->numNodes] = t;
  upHeap(*h, h->numNodes);
}

heapNode removeNode(heap *h)
{
  heapNode t;
  t = h->data[1];
  h->data[1] = h->data[h->numNodes];
  h->numNodes--;
  assert(h->numNodes >= 0);
  downHeap(*h, 1);
  return t;
}

#define MARGIN 70
#define HEIGHT (792 - MARGIN)

void drawPair(graph g, int i, int j)
{
  float px, py, cx, cy;
  px = MARGIN + g.x[i];
  py = HEIGHT - g.y[i];
  cx = MARGIN + g.x[j];
  cy = HEIGHT - g.y[j];
  printf("%f %f moveto %f %f lineto stroke\n", px, py, cx, cy);
  if (!g.drawn[i])
  {
    printf("%f %f %f 0 360 arc fill\n",px,py, RADIUS);
    g.drawn[i] = 1;
  }
  if (!g.drawn[j])
  {
    printf("%f %f %f 0 360 arc fill\n",cx,cy, RADIUS);
    g.drawn[j] = 1;
  }
}

void printGraph(graph g)
{
  int i, j;
  for (i = 0; i < g.numVertices; i++)
    g.drawn[i] = 0;

  printf("%%!\n");
  printf("0 setgray\n");

  for (i = 0; i < g.numVertices; i++)
    for (j = i + 1; j < g.numVertices; j++)
      if (g.weights[i * g.numVertices + j])
      {
#ifdef DEBUG
        printf("%d %d ", i, j);
#endif
        drawPair(g, i, j);
      }
  printf("showpage\n");
}