extern double fpow32(double);
extern Ppolyline_t getPath(edge_t *, vconfig_t *, int, Ppoly_t **,
int);
- extern void heapdown(Agnode_t *);
- extern void heapup(Agnode_t *);
+ extern void heapdown(graph_t *, Agnode_t *);
+ extern void heapup(graph_t *, Agnode_t *);
extern void initial_positions(graph_t *, int);
extern int init_port(Agnode_t *, Agedge_t *, char *, boolean);
extern void jitter3d(Agnode_t *, int);
extern void move_node(graph_t *, int, Agnode_t *);
extern int init_nop(graph_t * g, int);
extern void neato_cleanup(graph_t * g);
- extern node_t *neato_dequeue(void);
- extern void neato_enqueue(node_t *);
+ extern node_t *neato_dequeue(graph_t *);
+ extern void neato_enqueue(graph_t *, node_t *);
extern void neato_init_node(node_t * n);
extern void neato_layout(Agraph_t * g);
extern int Plegal_arrangement(Ppoly_t ** polys, int n_polys);
#include "sgd.h"
#include "neatoprocs.h"
#include <math.h>
+#include <stdlib.h>
typedef struct term {
node_t *i, *j;
float d, w;
} term;
-void fisheryates_shuffle(term *terms, int nC2) {
- int i;
- for (i=nC2-1; i>=1; i--) {
- // unsigned j = rk_interval(i, &rstate);
- int j = (int)(drand48() * (i+1)); // TODO: better rng
-
- term temp = terms[i];
- terms[i] = terms[j];
- terms[j] = temp;
- }
-}
-float calculate_stress(term *terms, int nC2) {
+float calculate_stress(term *terms, int n_terms) {
float stress = 0;
int ij;
- for (ij=0; ij<nC2; ij++) {
+ for (ij=0; ij<n_terms; ij++) {
float dx = ND_pos(terms[ij].i)[0] - ND_pos(terms[ij].j)[0];
float dy = ND_pos(terms[ij].i)[1] - ND_pos(terms[ij].j)[1];
float r = sqrt(dx*dx + dy*dy) - terms[ij].d;
}
return stress;
}
+void fisheryates_shuffle(term *terms, int n_terms) {
+ int i;
+ for (i=n_terms-1; i>=1; i--) {
+ int j = (int)(drand48() * (i+1)); // TODO: better rng?
+
+ term temp = terms[i];
+ terms[i] = terms[j];
+ terms[j] = temp;
+ }
+}
+// single source shortest paths that also builds terms as it goes
+// mostly copied from from stuff.c
+void dijkstra_single_source(graph_t *G, node_t *source, term *terms, int *offset) {
+ int t;
+ node_t *v, *u;
+ for (t = 0; (v = GD_neato_nlist(G)[t]); t++) {
+ ND_dist(v) = Initial_dist;
+ ND_heapindex(v) = -1;
+ }
+
+ ND_dist(source) = 0;
+ ND_hops(source) = 0;
+ neato_enqueue(G, source);
+
+ edge_t *e;
+ while ((v = neato_dequeue(G))) {
+ // if the target is fixed then always create a term as shortest paths are not calculated from there
+ // if not fixed then only create a term if the target index is lower
+ if (isFixed(v) || ND_id(v)<ND_id(source)) {
+ terms[*offset].i = source;
+ terms[*offset].j = v;
+ terms[*offset].d = ND_dist(v);
+ terms[*offset].w = 1 / (terms[*offset].d*terms[*offset].d);
+ (*offset)++;
+ }
+ for (e = agfstedge(G, v); e; e = agnxtedge(G, e, v)) {
+ if ((u = agtail(e)) == v)
+ u = aghead(e);
+
+ double f = ND_dist(v) + ED_dist(e);
+ if (ND_dist(u) > f) {
+ ND_dist(u) = f;
+ if (ND_heapindex(u) >= 0) {
+ heapup(G, u);
+ } else {
+ ND_hops(u) = ND_hops(v) + 1;
+ neato_enqueue(G, u);
+ }
+ }
+ }
+ }
+}
void sgd(graph_t *G, /* input graph */
- int n, /* number of nodes */
- int model /* distance model */)
+ int n, /* number of nodes */
+ int model /* distance model */)
{
- // calculate shortest paths (from neatoprocs)
- shortest_path(G, n);
+ if (model == MODEL_SHORTPATH) {
+ agerr(AGWARN, "shortest path model not yet supported in Gmode=sgd, reverting to MDS model\n");
+ } else if (model == MODEL_SUBSET) {
+ agerr(AGWARN, "subset model not yet supported in Gmode=sgd, reverting to MDS model\n");
+ } else if (model == MODEL_CIRCUIT) {
+ agerr(AGWARN, "circuit model not yet supported in Gmode=sgd, reverting to MDS model\n");
+ }
- // initialise array of terms
- int nC2 = (n*(n-1))/2;
- term *terms = N_NEW(nC2, term);
- int i, j;
- int ij = 0;
+ // calculate how many terms will be needed as fixed nodes can be ignored
+ int i, n_fixed = 0, n_terms = 0;
for (i=0; i<n; i++) {
- for (j=0; j<i; j++) {
- terms[ij].i = GD_neato_nlist(G)[i];
- terms[ij].j = GD_neato_nlist(G)[j];
- terms[ij].d = GD_dist(G)[ND_id(terms[ij].i)][ND_id(terms[ij].j)];
- terms[ij].w = 1 / (terms[ij].d*terms[ij].d);
- ij++;
+ if (!isFixed(GD_neato_nlist(G)[i])) {
+ n_fixed++;
+ n_terms += n-n_fixed;
}
}
+ term *terms = N_NEW(n_terms, term);
+
+ // calculate term values through shortest paths
+ int offset = 0;
+ GD_heap(G) = N_NEW(n + 1, node_t *);
+ GD_heapsize(G) = 0;
+ if (Verbose) {
+ fprintf(stderr, "calculating shortest paths:");
+ start_timer();
+ }
+ for (i=0; i<n; i++) {
+ node_t *source = GD_neato_nlist(G)[i];
+ if (!isFixed(source))
+ dijkstra_single_source(G, source, terms, &offset);
+ }
+ free(GD_heap(G));
+ assert(offset == n_terms);
+ if (Verbose) {
+ fprintf(stderr, " %.2f\n", elapsed_sec());
+ }
// initialise annealing schedule
float w_min = terms[0].w, w_max = terms[0].w;
- for (ij=1; ij<nC2; ij++) {
+ int ij;
+ for (ij=1; ij<n_terms; ij++) {
if (terms[ij].w < w_min)
w_min = terms[ij].w;
if (terms[ij].w > w_max)
w_max = terms[ij].w;
}
+ // note: Epsilon is different from MODE_KK and MODE_MAJOR as it is a minimum step size rather than energy threshold
+ // MaxIter is also different as it is a fixed number of iterations rather than a maximum
float eta_max = 1 / w_min;
- float eta_min = .01 / w_max; // TODO: read in epsilon
- float lambda = log(eta_max/eta_min) / (30-1); // TODO: read in t_max
+ float eta_min = Epsilon / w_max;
+ float lambda = log(eta_max/eta_min) / (MaxIter-1);
// initialise starting positions (from neatoprocs)
initial_positions(G, n);
// perform optimisation
int t;
if (Verbose) {
- fprintf(stderr, "Solving model:");
+ fprintf(stderr, "solving model:");
start_timer();
}
- for (t=0; t<30; t++) {
- fisheryates_shuffle(terms, nC2);
+ for (t=0; t<MaxIter; t++) {
+ fisheryates_shuffle(terms, n_terms);
float eta = eta_max * exp(-lambda * t);
- for (ij=0; ij<nC2; ij++) {
+ for (ij=0; ij<n_terms; ij++) {
// cap step size
float mu = eta * terms[ij].w;
if (mu > 1)
mu = 1;
-
+
float dx = ND_pos(terms[ij].i)[0] - ND_pos(terms[ij].j)[0];
float dy = ND_pos(terms[ij].i)[1] - ND_pos(terms[ij].j)[1];
float mag = sqrt(dx*dx + dy*dy);
-
+
float r = (mu * (mag-terms[ij].d)) / (2*mag);
float r_x = r * dx;
float r_y = r * dy;
- if (!ND_pinned(terms[ij].i)) {
+ if (!isFixed(terms[ij].i)) {
ND_pos(terms[ij].i)[0] -= r_x;
ND_pos(terms[ij].i)[1] -= r_y;
}
- if (!ND_pinned(terms[ij].j)) {
+ if (!isFixed(terms[ij].j)) {
ND_pos(terms[ij].j)[0] += r_x;
ND_pos(terms[ij].j)[1] += r_y;
}
}
if (Verbose) {
- fprintf(stderr, " %.3f", calculate_stress(terms, nC2));
+ fprintf(stderr, " %.3f", calculate_stress(terms, n_terms));
}
}
if (Verbose) {
- fprintf(stderr, "\n %.2f sec\n", elapsed_sec());
+ fprintf(stderr, "\nfinished in %.2f sec\n", elapsed_sec());
}
free(terms);
}
nE = agnedges(G);
lenx = agattr(G, AGEDGE, "len", 0);
- if (mode == MODE_KK || mode == MODE_SGD) {
+ if (mode == MODE_KK) {
Epsilon = .0001 * nV;
getdouble(G, "epsilon", &Epsilon);
if ((str = agget(G->root, "Damping")))
ND_heapindex(np) = -1;
total_len += setEdgeLen(G, np, lenx, dfltlen);
}
+ } else if (mode == MODE_SGD) {
+ Epsilon = .01;
+ getdouble(G, "epsilon", &Epsilon);
+ GD_neato_nlist(G) = N_NEW(nV + 1, node_t *);
+ for (i = 0, np = agfstnode(G); np; np = agnxtnode(G, np)) {
+ GD_neato_nlist(G)[i] = np;
+ ND_id(np) = i++;
+ total_len += setEdgeLen(G, np, lenx, dfltlen);
+ }
} else {
Epsilon = DFLT_TOLERANCE;
getdouble(G, "epsilon", &Epsilon);
else
Initial_dist = total_len / (nE > 0 ? nE : 1) * sqrt(nV) + 1;
- if (!Nop && (mode == MODE_KK || mode == MODE_SGD)) {
+ if (!Nop && (mode == MODE_KK)) {
GD_dist(G) = new_array(nV, nV, Initial_dist);
GD_spring(G) = new_array(nV, nV, 1.0);
GD_sum_t(G) = new_array(nV, Ndim, 1.0);
}
}
-static node_t **Heap;
-static int Heapsize;
-static node_t *Src;
-
-void heapup(node_t * v)
+void heapup(graph_t * G, node_t * v)
{
int i, par;
node_t *u;
for (i = ND_heapindex(v); i > 0; i = par) {
par = (i - 1) / 2;
- u = Heap[par];
+ u = GD_heap(G)[par];
if (ND_dist(u) <= ND_dist(v))
break;
- Heap[par] = v;
+ GD_heap(G)[par] = v;
ND_heapindex(v) = par;
- Heap[i] = u;
+ GD_heap(G)[i] = u;
ND_heapindex(u) = i;
}
}
-void heapdown(node_t * v)
+void heapdown(graph_t * G, node_t * v)
{
int i, left, right, c;
node_t *u;
i = ND_heapindex(v);
- while ((left = 2 * i + 1) < Heapsize) {
+ while ((left = 2 * i + 1) < GD_heapsize(G)) {
right = left + 1;
- if ((right < Heapsize)
- && (ND_dist(Heap[right]) < ND_dist(Heap[left])))
+ if ((right < GD_heapsize(G))
+ && (ND_dist(GD_heap(G)[right]) < ND_dist(GD_heap(G)[left])))
c = right;
else
c = left;
- u = Heap[c];
+ u = GD_heap(G)[c];
if (ND_dist(v) <= ND_dist(u))
break;
- Heap[c] = v;
+ GD_heap(G)[c] = v;
ND_heapindex(v) = c;
- Heap[i] = u;
+ GD_heap(G)[i] = u;
ND_heapindex(u) = i;
i = c;
}
}
-void neato_enqueue(node_t * v)
+void neato_enqueue(graph_t * G, node_t * v)
{
int i;
assert(ND_heapindex(v) < 0);
- i = Heapsize++;
+ i = GD_heapsize(G)++;
ND_heapindex(v) = i;
- Heap[i] = v;
+ GD_heap(G)[i] = v;
if (i > 0)
- heapup(v);
+ heapup(G, v);
}
-node_t *neato_dequeue(void)
+node_t *neato_dequeue(graph_t * G)
{
int i;
node_t *rv, *v;
- if (Heapsize == 0)
+ if (GD_heapsize(G) == 0)
return NULL;
- rv = Heap[0];
- i = --Heapsize;
- v = Heap[i];
- Heap[0] = v;
+ rv = GD_heap(G)[0];
+ i = --GD_heapsize(G);
+ v = GD_heap(G)[i];
+ GD_heap(G)[0] = v;
ND_heapindex(v) = 0;
if (i > 1)
- heapdown(v);
+ heapdown(G, v);
ND_heapindex(rv) = -1;
return rv;
}
{
node_t *v;
- Heap = N_NEW(nG + 1, node_t *);
+ GD_heap(G) = N_NEW(nG + 1, node_t *);
+ GD_heapsize(G) = 0;
if (Verbose) {
fprintf(stderr, "Calculating shortest paths: ");
start_timer();
if (Verbose) {
fprintf(stderr, "%.2f sec\n", elapsed_sec());
}
- free(Heap);
+ free(GD_heap(G));
}
void s1(graph_t * G, node_t * node)
for (t = 0; (v = GD_neato_nlist(G)[t]); t++)
ND_dist(v) = Initial_dist;
- Src = node;
+ node_t * Src = node;
ND_dist(Src) = 0;
ND_hops(Src) = 0;
- neato_enqueue(Src);
+ neato_enqueue(G, Src);
- while ((v = neato_dequeue())) {
+ while ((v = neato_dequeue(G))) {
if (v != Src)
make_spring(G, Src, v, ND_dist(v));
for (e = agfstedge(G, v); e; e = agnxtedge(G, e, v)) {
if (ND_dist(u) > f) {
ND_dist(u) = f;
if (ND_heapindex(u) >= 0)
- heapup(u);
+ heapup(G, u);
else {
ND_hops(u) = ND_hops(v) + 1;
- neato_enqueue(u);
+ neato_enqueue(G, u);
}
}
}
projects / graphviz / commitdiff