freeHeap(&H);
free(index);
}
+
+// single source shortest paths that also builds terms as it goes
+// mostly copied from dijkstra_f above
+// returns the number of terms built
+int dijkstra_sgd(graph_sgd *graph, int source, term_sgd *terms) {
+ heap h;
+#ifdef OBSOLETE
+ mkHeap(&h, graph->n);
+#endif
+ int *indices = N_GNEW(graph->n, int);
+ float *dists = N_GNEW(graph->n, float);
+ int i;
+ for (i=0; i<graph->n; i++) {
+ dists[i] = MAXFLOAT;
+ }
+ dists[source] = 0;
+ for (i=graph->sources[source]; i<graph->sources[source+1]; i++) {
+ int target = graph->targets[i];
+ dists[target] = graph->weights[target];
+ }
+ initHeap_f(&h, source, indices, dists, graph->n);
+
+ int closest = 0, offset = 0;
+ while (extractMax_f(&h, &closest, indices, dists)) {
+ float d = dists[closest];
+ if (d == MAXFLOAT) {
+ break;
+ }
+ // 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 (graph->pinneds[closest] || closest<source) {
+ terms[offset].i = source;
+ terms[offset].j = closest;
+ terms[offset].d = d;
+ terms[offset].w = 1 / (d*d);
+ offset++;
+ }
+ for (i=graph->sources[closest]; i<graph->sources[closest+1]; i++) {
+ int target = graph->targets[i];
+ int weight = graph->weights[i];
+ increaseKey_f(&h, target, d+weight, indices, dists);
+ }
+ }
+ freeHeap(&h);
+ free(indices);
+ free(dists);
+ return offset;
+}
+
#include "sgd.h"
+#include "dijkstra.h"
+#include "randomkit.h"
#include "neatoprocs.h"
#include <math.h>
#include <stdlib.h>
-#include "randomkit.h"
-
-typedef struct term {
- node_t *i, *j;
- float d, w;
-} term;
-// TODO: convert to just indices rather than node_t, because working with the auxiliary graph_t is slow
-float calculate_stress(term *terms, int n_terms) {
+float calculate_stress(float *pos, term_sgd *terms, int n_terms) {
float stress = 0;
int 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 dx = pos[2*terms[ij].i] - pos[2*terms[ij].j];
+ float dy = pos[2*terms[ij].i+1] - pos[2*terms[ij].j+1];
float r = sqrt(dx*dx + dy*dy) - terms[ij].d;
stress += terms[ij].w * (r * r);
}
}
// it is much faster to shuffle term rather than pointers to term, even though the swap is more expensive
static rk_state rstate;
-static void fisheryates_shuffle(term *terms, int n_terms) {
+static void fisheryates_shuffle(term_sgd *terms, int n_terms) {
int i;
for (i=n_terms-1; i>=1; i--) {
// srand48() is called in neatoinit.c, so no need to seed here
//int j = (int)(drand48() * (i+1));
int j = rk_interval(i, &rstate);
- term temp = terms[i];
+ term_sgd 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
-// returns the number of terms built
-int dijkstra_single_source(graph_t *G, node_t *source, term *terms) {
- 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;
- neato_enqueue(G, source);
- edge_t *e;
- int offset = 0;
- 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++;
+// to make dijkstras faster
+graph_sgd * extract_adjacency(graph_t *G) {
+ node_t *np;
+ edge_t *ep;
+ int n_nodes = 0, n_edges = 0;
+ for (np = agfstnode(G); np; np = agnxtnode(G,np)) {
+ assert(ND_id(np) == n_nodes);
+ n_nodes++;
+ for (ep = agfstedge(G, np); ep; ep = agnxtedge(G, ep, np)) {
+ n_edges++;
}
- 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 {
- neato_enqueue(G, u);
- }
- }
+ }
+ graph_sgd *graph = N_NEW(1, graph_sgd);
+ graph->sources = N_NEW(n_nodes+1, int);
+ graph->pinneds = N_NEW(n_nodes, bool);
+ graph->targets = N_NEW(n_edges, int);
+ graph->weights = N_NEW(n_edges, float);
+
+ graph->n = n_nodes;
+ graph->sources[graph->n] = n_edges; // to make looping nice
+
+ n_nodes = 0, n_edges = 0;
+ for (np = agfstnode(G); np; np = agnxtnode(G,np)) {
+ graph->sources[n_nodes] = n_edges;
+ graph->pinneds[n_nodes] = isFixed(np);
+ for (ep = agfstedge(G, np); ep; ep = agnxtedge(G, ep, np)) {
+ node_t *target = (agtail(ep) == np) ? aghead(ep) : agtail(ep); // in case edge is reversed
+ graph->targets[n_edges] = ND_id(target);
+ graph->weights[n_edges] = ED_dist(ep);
+ /*if (model == MODEL_MDS) {
+ graph->weights[n_edges] = ED_dist(ep);
+ } else if (model == MODEL_SHORTPATH) {
+ graph->weights[n_edges] = 1;
+ } else {
+ assert(false); // not supported
+ }*/
+ assert(ED_dist(ep) > 0);
+ n_edges++;
}
+ n_nodes++;
}
- return offset;
+ assert(n_nodes == graph->n);
+ assert(n_edges == graph->sources[graph->n]);
+
+ graph->sources[n_nodes] = n_edges;
+ return graph;
+}
+void free_adjacency(graph_sgd *graph) {
+ free(graph->sources);
+ free(graph->pinneds);
+ free(graph->targets);
+ free(graph->weights);
+ free(graph);
}
+
void sgd(graph_t *G, /* input graph */
- int n, /* number of nodes */
int model /* distance model */)
{
- 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) {
+ 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");
}
+ int n = agnnodes(G);
+ if (Verbose) {
+ fprintf(stderr, "calculating shortest paths:");
+ start_timer();
+ }
// 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++) {
n_terms += n-n_fixed;
}
}
- term *terms = N_NEW(n_terms, term);
-
+ term_sgd *terms = N_NEW(n_terms, term_sgd);
// calculate term values through shortest paths
int offset = 0;
- GD_heap(G) = N_NEW(n, node_t *);
- GD_heapsize(G) = 0;
- if (Verbose) {
- fprintf(stderr, "calculating shortest paths:");
- start_timer();
- }
+ graph_sgd *graph = extract_adjacency(G);
for (i=0; i<n; i++) {
- node_t *source = GD_neato_nlist(G)[i];
- if (!isFixed(source)) {
- offset += dijkstra_single_source(G, source, terms+offset);
+ if (!isFixed(GD_neato_nlist(G)[i])) {
+ offset += dijkstra_sgd(graph, i, terms+offset);
}
}
- free(GD_heap(G));
assert(offset == n_terms);
+ free(graph);
if (Verbose) {
fprintf(stderr, " %.2f\n", elapsed_sec());
}
// initialise starting positions (from neatoprocs)
initial_positions(G, n);
+ // copy initial positions and state into temporary space for speed
+ float *pos = N_NEW(2*n, float);
+ bool *unfixed = N_NEW(n, bool);
+ for (i=0; i<n; i++) {
+ node_t *node = GD_neato_nlist(G)[i];
+ pos[2*i] = ND_pos(node)[0];
+ pos[2*i+1] = ND_pos(node)[1];
+ unfixed[i] = !isFixed(node);
+ }
// perform optimisation
if (Verbose) {
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 dx = pos[2*terms[ij].i] - pos[2*terms[ij].j];
+ float dy = pos[2*terms[ij].i+1] - pos[2*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 (!isFixed(terms[ij].i)) {
- ND_pos(terms[ij].i)[0] -= r_x;
- ND_pos(terms[ij].i)[1] -= r_y;
+ if (unfixed[terms[ij].i]) {
+ pos[2*terms[ij].i] -= r_x;
+ pos[2*terms[ij].i+1] -= r_y;
}
- if (!isFixed(terms[ij].j)) {
- ND_pos(terms[ij].j)[0] += r_x;
- ND_pos(terms[ij].j)[1] += r_y;
+ if (unfixed[terms[ij].j]) {
+ pos[2*terms[ij].j] += r_x;
+ pos[2*terms[ij].j+1] += r_y;
}
}
if (Verbose) {
- fprintf(stderr, " %.3f", calculate_stress(terms, n_terms));
+ fprintf(stderr, " %.3f", calculate_stress(pos, terms, n_terms));
}
}
if (Verbose) {
fprintf(stderr, "\nfinished in %.2f sec\n", elapsed_sec());
}
free(terms);
+
+ // copy temporary positions back into graph_t
+ for (i=0; i<n; i++) {
+ node_t *node = GD_neato_nlist(G)[i];
+ ND_pos(node)[0] = pos[2*i];
+ ND_pos(node)[1] = pos[2*i+1];
+ }
+ free(pos);
+ free(unfixed);
}
projects / graphviz / commitdiff