From: Jonathan Zheng <jonathanzheng@hotmail.co.uk>
Date: Thu, 16 Jan 2020 21:01:45 +0000 (+0000)
Subject: dijkstra is now fast!
X-Git-Tag: stable_release_2.44.0~12^2~12
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=acab13dfc346a81d12ffe88791216bbc53132e16;p=graphviz

dijkstra is now fast!
---

diff --git a/lib/neatogen/dijkstra.c b/lib/neatogen/dijkstra.c
index 4fd8605e5..198523556 100644
--- a/lib/neatogen/dijkstra.c
+++ b/lib/neatogen/dijkstra.c
@@ -394,3 +394,52 @@ void dijkstra_f(int vertex, vtx_data * graph, int n, float *dist)
     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;
+}
+
diff --git a/lib/neatogen/dijkstra.h b/lib/neatogen/dijkstra.h
index 2eb3894e8..c3aaeeea9 100644
--- a/lib/neatogen/dijkstra.h
+++ b/lib/neatogen/dijkstra.h
@@ -21,6 +21,7 @@ extern "C" {
 #define _DIJKSTRA_H_
 
 #include "defs.h"
+#include "sgd.h"
 
 #ifdef __cplusplus
     void dijkstra(int vertex, vtx_data * graph, int n, DistType * dist);
@@ -38,6 +39,7 @@ extern "C" {
     extern int dijkstra_bounded(int, vtx_data *, int, DistType *, int,
 				int *);
 #endif
+    extern int dijkstra_sgd(graph_sgd *, int, term_sgd *);
 
 #endif
 
diff --git a/lib/neatogen/neatoinit.c b/lib/neatogen/neatoinit.c
index 2f63bb90c..7d9e3c464 100644
--- a/lib/neatogen/neatoinit.c
+++ b/lib/neatogen/neatoinit.c
@@ -1362,7 +1362,7 @@ neatoLayout(Agraph_t * mg, Agraph_t * g, int layoutMode, int layoutModel,
     if (layoutMode == MODE_KK)
 	kkNeato(g, nG, layoutModel);
     else if (layoutMode == MODE_SGD)
-	sgd(g, nG, layoutModel);
+	sgd(g, layoutModel);
     else
 	majorization(mg, g, nG, layoutMode, layoutModel, Ndim, MaxIter, am);
 }
diff --git a/lib/neatogen/sgd.c b/lib/neatogen/sgd.c
index 21269b15f..5a470b928 100644
--- a/lib/neatogen/sgd.c
+++ b/lib/neatogen/sgd.c
@@ -1,22 +1,17 @@
 #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);
     }
@@ -24,74 +19,89 @@ float calculate_stress(term *terms, int n_terms) {
 }
 // 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++) {
@@ -100,24 +110,17 @@ void sgd(graph_t *G, /* input graph */
             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());
     }
@@ -139,6 +142,15 @@ void sgd(graph_t *G, /* input graph */
 
     // 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) {
@@ -156,29 +168,38 @@ void sgd(graph_t *G, /* input graph */
             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);
 }
diff --git a/lib/neatogen/sgd.h b/lib/neatogen/sgd.h
index b99263ece..9b2b12f07 100644
--- a/lib/neatogen/sgd.h
+++ b/lib/neatogen/sgd.h
@@ -3,6 +3,20 @@
 
 #include "neato.h"
 
-extern void sgd(graph_t *, int, int);
+typedef struct term_sgd {
+    int i, j;
+    float d, w;
+} term_sgd;
+
+typedef struct graph_sgd {
+    int n; // number of nodes
+    int *sources; // index of first edge in *targets for each node (length n+1)
+    bool *pinneds; // whether a node is fixed or not
+
+    int *targets; // index of targets for each node (length sources[n])
+    float *weights; // weights of edges (length sources[n])
+} graph_sgd;
+
+extern void sgd(graph_t *, int);
 
 #endif /* SGD_H */