Add hierarchical constraints. This is enabled only if the compile-time

variable DIGCOLA is defined. At present, by default this is undefined.

diff --git a/lib/neatogen/matrix_ops.c b/lib/neatogen/matrix_ops.c
index 024cae5f9f5d7f6b6ac42d8a2e7a48cd2d195989..5c85f56ba3355f88880a12c9144404c44de9fb33 100644 (file)
--- a/lib/neatogen/matrix_ops.c
+++ b/lib/neatogen/matrix_ops.c
@@ -676,13 +676,20 @@ double vectors_inner_productf(int n, float *vector1, float *vector2)
 }
 
 /* inline */
-void set_vector_val(int n, float val, float *result)
+void set_vector_val(int n, double val, double *result)
 {
     int i;
     for (i = 0; i < n; i++)
        result[i] = val;
 }
 
+/* inline */
+void set_vector_valf(int n, float val, float * result) {
+    int i;
+    for (i=0; i<n; i++)
+        result[i]=val;
+}
+
 /* inline */
 double max_absf(int n, float *vector)
 {

diff --git a/lib/neatogen/matrix_ops.h b/lib/neatogen/matrix_ops.h
index cfe21a1f7060770c2096e5efcbce8655640d29b7..f851036d45b1858b50225de37da31b2a2b253133 100644 (file)
--- a/lib/neatogen/matrix_ops.h
+++ b/lib/neatogen/matrix_ops.h
@@ -91,7 +91,8 @@ extern "C" {
     extern void copy_vectorf(int n, float *source, float *dest);
     extern double vectors_inner_productf(int n, float *vector1,
                                         float *vector2);
-    extern void set_vector_val(int n, float val, float *result);
+    extern void set_vector_val(int n, double val, double *result);
+    extern void set_vector_valf(int n, float val, float * result);
     extern double max_absf(int n, float *vector);
     extern void square_vec(int n, float *vec);
     extern void invert_vec(int n, float *vec);

diff --git a/lib/neatogen/neato.h b/lib/neatogen/neato.h
index 81904b6995fe822ceaa797ef6db4a3057f3e80b2..7671d95cb4efa4a9bf0bcae20ad5829147a41568 100644 (file)
--- a/lib/neatogen/neato.h
+++ b/lib/neatogen/neato.h
@@ -28,6 +28,7 @@
 
 #define MODE_KK          0
 #define MODE_MAJOR       1
+#define MODE_HIER        2
 
 #define INIT_SELF        0
 #define INIT_REGULAR     1

diff --git a/lib/neatogen/neatoinit.c b/lib/neatogen/neatoinit.c
index c3ce0df6a7162afc1d4d084a9badeafe7d58233e..afe0cbcc64ab3335ae8cc50cb836fcbc2b669f3b 100644 (file)
--- a/lib/neatogen/neatoinit.c
+++ b/lib/neatogen/neatoinit.c
@@ -26,6 +26,9 @@
 #include "neato.h"
 #include "pack.h"
 #include "stress.h"
+#ifdef DIGCOLA
+#include "digcola.h"
+#endif
 #include "kkutils.h"
 #include "pointset.h"
 #include <ctype.h>
@@ -628,6 +631,10 @@ static int neatoMode(graph_t * g)
            mode = MODE_KK;
        else if (streq(str, "major"))
            mode = MODE_MAJOR;
+#ifdef DIGCOLA
+       else if (streq(str, "hier"))
+           mode = MODE_HIER;
+#endif
        else
            agerr(AGWARN,
                  "Illegal value %s for attribute \"mode\" in graph %s - ignored\n",
@@ -654,6 +661,57 @@ static int checkEdge(PointMap * pm, edge_t * ep, int idx)
     return insertPM(pm, i, j, idx);
 }
 
+#ifdef DIGCOLA
+/* dfsCycle:
+ * dfs for breaking cycles in vtxdata
+ */
+static void
+dfsCycle (vtx_data* graph, int i)
+{
+    node_t* np = graph[i].np;
+    node_t* hp;
+    ND_mark(np) = TRUE;
+    ND_onstack(np) = TRUE;
+    int j, e, f;
+
+    for (e = 1; e < graph[i].nedges; e++) {
+       if (graph[i].edists[e] == 1.0) continue;  /* in edge */
+       j = graph[i].edges[e];
+       hp = graph[j].np;
+       if (ND_onstack(hp)) {  /* back edge: reverse it */
+           graph[i].edists[e] = 1.0;
+           for (f = 1; (f < graph[j].nedges) &&(graph[j].edges[f] != i); f++) ;
+           assert (f < graph[j].nedges);
+           graph[j].edists[f] = -1.0;
+        }
+       else if (ND_mark(hp) == FALSE) dfsCycle(graph, j);
+
+    }
+    ND_onstack(np) = FALSE;
+}
+
+/* acyclic:
+ * Do a dfs of the vtx_data, looking for cycles, reversing edges.
+ */
+static void
+acyclic (vtx_data* graph, int nv)
+{
+    int i;
+    node_t* np;
+
+    for (i = 0; i < nv; i++) {
+       np = graph[i].np;
+       ND_mark(np) = FALSE;
+       ND_onstack(np) = FALSE;
+    }
+    for (i = 0; i < nv; i++) {
+       if (ND_mark(graph[i].np)) continue;
+       dfsCycle (graph, i);    
+    }
+
+}
+#endif
+
 /* makeGraphData:
  * Create sparse graph representation via arrays.
  * Each node is represented by a vtx_data.
@@ -663,10 +721,16 @@ static int checkEdge(PointMap * pm, edge_t * ep, int idx)
  * By convention, graph[i].edges[0] == i.
  * The values graph[i].ewgts[0] and graph[i].eweights[0] are left undefined.
  *
- * In construct graph from g, we neglect loops. We track multiedges (ignoring
+ * In constructing graph from g, we neglect loops. We track multiedges (ignoring
  * direction). Edge weights are additive; the final edge length is the max.
+ *
+ * If direction is used, we set the edists field, -1 for tail, +1 for head. 
+ * graph[i].edists[0] is left undefined. If multiedges exist, the direction
+ * of the first one encountered is used. Finally, a pass is made to guarantee
+ * the graph is acyclic.
+ *
  */
-static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int model)
+static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int mode, int model)
 {
     vtx_data *graph;
     int ne = agnedges(g);      /* upper bound */
@@ -675,8 +739,12 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int model)
     node_t *np;
     edge_t *ep;
     float *eweights = NULL;
+#ifdef DIGCOLA
+    float *edists = NULL;
+#endif
     int haveLen;
     int haveWt;
+    int haveDir;
     PointMap *ps = newPM();
     int i, i_nedges, idx;
 
@@ -688,13 +756,21 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int model)
        haveLen = (agindex(g->root->proto->e, "len") >= 0);
        haveWt = (E_weight != 0);
     }
+    if (mode == MODE_HIER)
+       haveDir = TRUE;
+    else
+       haveDir = FALSE;
 
     graph = N_GNEW(nv, vtx_data);
     edges = N_GNEW(2 * ne + nv, int);  /* reserve space for self loops */
-    if (haveLen)
+    if (haveLen || haveDir)
        ewgts = N_GNEW(2 * ne + nv, float);
     if (haveWt)
        eweights = N_GNEW(2 * ne + nv, float);
+#ifdef DIGCOLA
+    if (haveDir)
+       edists = N_GNEW(2*ne+nv,float);
+#endif
 
     i = 0;
     ne = 0;
@@ -704,7 +780,7 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int model)
        assert(ND_id(np) == i);
        graph[i].np = np;
        graph[i].edges = edges++;       /* reserve space for the self loop */
-       if (haveLen)
+       if (haveLen || haveDir)
            graph[i].ewgts = ewgts++;
        else
            graph[i].ewgts = NULL;
@@ -712,6 +788,13 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int model)
            graph[i].eweights = eweights++;
        else
            graph[i].eweights = NULL;
+#ifdef DIGCOLA
+       if (haveDir) {
+           graph[i].edists = edists++;
+       }
+       else
+           graph[i].edists = NULL;
+#endif
        i_nedges = 1;           /* one for the self */
 
        for (ep = agfstedge(g, np); ep; ep = agnxtedge(g, ep, np)) {
@@ -735,6 +818,13 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int model)
                    *eweights++ = ED_factor(ep);
                if (haveLen)
                    *ewgts++ = ED_dist(ep);
+               else if (haveDir)
+                   *ewgts++ = 1.0;
+#ifdef DIGCOLA
+               if (haveDir) {
+                   *edists++ = (np == ep->head ? 1.0 : -1.0);
+               }
+#endif
                i_nedges++;
            }
        }
@@ -746,6 +836,12 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int model)
 #endif
        i++;
     }
+#ifdef DIGCOLA
+    if (haveDir) {
+    /* Make graph acyclic */
+       acyclic (graph, nv);
+    }
+#endif
 
     ne /= 2;                   /* every edge is counted twice */
 
@@ -919,7 +1015,7 @@ majorization(graph_t * g, int nv, int mode, int model, int dim, int steps)
     vtx_data *gp;
     int init;
 
-    init = checkStart(g, nv, INIT_RANDOM);
+    init = checkStart(g, nv, (mode == MODE_HIER ? INIT_SELF : INIT_RANDOM));
     coords = N_GNEW(dim, double *);
     coords[0] = N_GNEW(nv * dim, double);
     for (i = 1; i < Ndim; i++) {
@@ -931,12 +1027,20 @@ majorization(graph_t * g, int nv, int mode, int model, int dim, int steps)
        fprintf(stderr, "convert graph: ");
        start_timer();
     }
-    gp = makeGraphData(g, nv, &ne, model);
+    gp = makeGraphData(g, nv, &ne, mode, model);
     if (Verbose) {
        fprintf(stderr, "%d nodes %.2f sec\n", nv, elapsed_sec());
     }
 
-    stress_majorization_kD_mkernel(gp, nv, ne, coords, Ndim,
+#ifdef DIGCOLA
+    if (mode == MODE_HIER) {
+       double lgap = late_double(g, agfindattr(g, "levelsgap"), 0.0, -MAXDOUBLE);
+       stress_majorization_with_hierarchy(gp, nv, ne, coords, Ndim,
+                                  (init == INIT_SELF), model, MaxIter, lgap);
+    }
+    else
+#endif
+       stress_majorization_kD_mkernel(gp, nv, ne, coords, Ndim,
                                   (init == INIT_SELF), model, MaxIter);
 
     /* store positions back in nodes */
@@ -958,7 +1062,7 @@ static void subset_model(Agraph_t * G, int nG)
     DistType **Dij;
     vtx_data *gp;
 
-    gp = makeGraphData(G, nG, &ne, MODEL_SUBSET);
+    gp = makeGraphData(G, nG, &ne, MODE_KK, MODEL_SUBSET);
     Dij = compute_apsp_artifical_weights(gp, nG);
     for (i = 0; i < nG; i++) {
        for (j = 0; j < nG; j++) {

diff --git a/lib/neatogen/stress.c b/lib/neatogen/stress.c
index fa8a8c13528be09dddca2950409cf7e64d6db7c6..a44c38a5f770713a6c7d40b5f8198da4da638cd9 100644 (file)
--- a/lib/neatogen/stress.c
+++ b/lib/neatogen/stress.c
@@ -86,9 +86,6 @@ extern double drand48(void);
 #define max_nodes_in_mem 18000
 #endif
 
- /* conjugate gradient related: */
-#define tolerance_cg 1e-3
-
  /* relevant when using sparse distance matrix not within subspace */
 #define smooth_pivots true
 
@@ -206,7 +203,8 @@ compute_stress1(double **coords, dist_data * distances, int dim, int n)
  * a position, use it.
  * Return true if some node is fixed.
  */
-static int initLayout(vtx_data * graph, int n, int dim, double **coords)
+int
+initLayout(vtx_data * graph, int n, int dim, double **coords)
 {
     node_t *np;
     double *xp;
@@ -245,7 +243,8 @@ static int initLayout(vtx_data * graph, int n, int dim, double **coords)
     return pinned;
 }
 
-static float *circuitModel(vtx_data * graph, int nG)
+float*
+circuitModel(vtx_data * graph, int nG)
 {
     int i, j, e, rv, count;
     float *Dij = N_NEW(nG * (nG + 1) / 2, float);
@@ -1362,7 +1361,7 @@ static float *compute_weighted_apsp_packed(vtx_data * graph, int n)
 /* compute_apsp_packed:
  * Assumes integral weights > 0.
  */
-static float *compute_apsp_packed(vtx_data * graph, int n)
+float *compute_apsp_packed(vtx_data * graph, int n)
 {
     int i, j, count;
     float *Dij = N_NEW(n * (n + 1) / 2, float);
@@ -1386,8 +1385,7 @@ static float *compute_apsp_packed(vtx_data * graph, int n)
 
 #define max(x,y) ((x)>(y)?(x):(y))
 
-static float *compute_apsp_artifical_weights_packed(vtx_data * graph,
-                                                   int n)
+float *compute_apsp_artifical_weights_packed(vtx_data * graph, int n)
 {
     /* compute all-pairs-shortest-path-length while weighting the graph */
     /* so high-degree nodes are distantly located */
@@ -1707,11 +1705,11 @@ int stress_majorization_kD_mkernel(vtx_data * graph,    /* Input graph in sparse re
        for (count = 0, i = 0; i < n - 1; i++) {
            len = n - i - 1;
            /* init 'dist_accumulator' with zeros */
-           set_vector_val(len, 0, dist_accumulator);
+           set_vector_valf(len, 0, dist_accumulator);
 
            /* put into 'dist_accumulator' all squared distances between 'i' and 'i'+1,...,'n'-1 */
            for (k = 0; k < dim; k++) {
-               set_vector_val(len, coords[k][i], tmp_coords);
+               set_vector_valf(len, coords[k][i], tmp_coords);
                vectors_mult_additionf(len, tmp_coords, -1,
                                       coords[k] + i + 1);
                square_vec(len, tmp_coords);

diff --git a/lib/neatogen/stress.h b/lib/neatogen/stress.h
index 1bd0c1b534e9ff4ca85acf6dea3e82efb2e6432e..26c50d7b6c482a184ab06a0adbbb5c77c1a2ae91 100644 (file)
--- a/lib/neatogen/stress.h
+++ b/lib/neatogen/stress.h
@@ -24,6 +24,8 @@ extern "C" {
 
 #include "defs.h"
 
+#define tolerance_cg 1e-3
+
 #define DFLT_ITERATIONS 200
 
 #define DFLT_TOLERANCE 1e-4
@@ -54,6 +56,10 @@ extern "C" {
                                              int maxi  /* max iterations */
        );
 
+extern float *compute_apsp_packed(vtx_data * graph, int n);
+extern float *compute_apsp_artifical_weights_packed(vtx_data * graph, int n);
+extern float* circuitModel(vtx_data * graph, int nG);
+extern int initLayout(vtx_data * graph, int n, int dim, double **coords);
 
 #endif