noinst_HEADERS = adjust.h edges.h geometry.h heap.h hedges.h info.h mem.h \
neato.h poly.h neatoprocs.h simple.h site.h voronoi.h \
bfs.h closest.h conjgrad.h defs.h dijkstra.h embed_graph.h kkutils.h \
- matrix_ops.h pca.h stress.h
+ matrix_ops.h pca.h stress.h quad_prog_solver.h
noinst_LTLIBRARIES = libneatogen.la
libneatogen_la_SOURCES = adjust.c circuit.c edges.c find_ints.c geometry.c \
heap.c hedges.c info.c neatoinit.c intersect.c legal.c lu.c matinv.c \
memory.c poly.c printvis.c site.c solve.c neatosplines.c stuff.c \
voronoi.c stress.c kkutils.c matrix_ops.c embed_graph.c dijkstra.c \
- conjgrad.c pca.c closest.c bfs.c constraint.c
+ conjgrad.c pca.c closest.c bfs.c constraint.c quad_prog_solve.c \
+ smart_ini_x.c constrained_majorization.c opt_arrangement.c \
+ compute_hierarchy.c
EXTRA_DIST = Makefile.old
neato.h neatoprocs.h poly.h simple.h site.h voronoi.h
NOBJS = stress.o kkutils.o pca.o matrix_ops.o embed_graph.o dijkstra.o \
- conjgrad.o closest.o bfs.o
+ conjgrad.o closest.o bfs.o compute_hierarchy.o opt_arrangement.o \
+ smart_ini_x.o constrained_majorization.o quad_prog_solve.o
OBJS = adjust.o circuit.o edges.o find_ints.o geometry.o heap.o hedges.o \
info.o neatoinit.o intersect.o legal.o lu.o matinv.o memory.o poly.o \
printvis.o site.o solve.o neatosplines.o stuff.o voronoi.o \
--- /dev/null
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#include <digcola.h>
+#ifdef DIGCOLA
+#include "kkutils.h"
+
+static int* given_levels = NULL;
+/*
+ * This function partitions the graph nodes into levels
+ * according to the minimizer of the hierarchy energy.
+ *
+ * To allow more flexibility we define a new level only
+ * when the hierarchy energy shows a significant jump
+ * (to compensate for noise).
+ * This is controlled by two parameters: 'abs_tol' and
+ * 'relative_tol'. The smaller these two are, the more
+ * levels we'll get.
+ * More speciffically:
+ * We never consider gaps smaller than 'abs_tol'
+ * Additionally, we never consider gaps smaller than 'abs_tol'*<avg_gap>
+ *
+ * The output is an ordering of the nodes according to
+ * their levels, as follows:
+ * First level:
+ * ordering[0],ordering[1],...ordering[levels[0]-1]
+ * Second level:
+ * ordering[levels[0]],ordering[levels[0]+1],...ordering[levels[1]-1]
+ * ...
+ * Last level:
+ * ordering[levels[num_levels-1]],ordering[levels[num_levels-1]+1],...ordering[n-1]
+ *
+ * Hence, the nodes were partitioned into 'num_levels'+1
+ * levels.
+ *
+ * Please note that 'ordering[levels[i]]' contains the first node at level i+1,
+ * and not the last node of level i.
+ */
+double
+compute_hierarchy(vtx_data* graph, int n, double abs_tol, double relative_tol,
+ double* given_coords, int** orderingp, int** levelsp, int* num_levelsp)
+{
+ double* y;
+ int i;
+ double spread;
+ int use_given_levels=FALSE;
+ int* ordering;
+ int* levels;
+ double tol; /* node 'i' precedes 'j' in hierachy iff y[i]-y[j]>tol */
+ double hierarchy_span;
+ int num_levels;
+
+ /* compute optimizer of hierarchy energy: 'y' */
+ if (given_coords) {
+ y = given_coords;
+ }
+ else {
+ y = N_GNEW(n,double);
+ compute_y_coords(graph, n, y, n);
+ }
+
+ /* sort nodes accoridng to their y-ordering */
+ *orderingp = ordering = N_NEW(n, int);
+ for (i=0; i<n; i++) {
+ ordering[i] = i;
+ }
+ quicksort_place(y, ordering, 0,n-1);
+
+ spread = y[ordering[n-1]]-y[ordering[0]];
+
+ /* after spread is computed, we may take the y-coords as the given levels */
+ if (given_levels) {
+ use_given_levels=TRUE;
+ for (i=0; i<n; i++) {
+ use_given_levels = use_given_levels && given_levels[i]>=0;
+ }
+ }
+ if (use_given_levels) {
+ for (i=0; i<n; i++) {
+ y[i] = given_levels[i];
+ ordering[i] = i;
+ }
+ quicksort_place(y, ordering, 0,n-1);
+ }
+
+ /* compute tolerance
+ * take the maximum between 'abs_tol' and a fraction of the average gap
+ * controlled by 'relative_tol'
+ */
+ hierarchy_span = y[ordering[n-1]]-y[ordering[0]];
+ tol = MAX(abs_tol, relative_tol*hierarchy_span/(n-1));
+ /* 'hierarchy_span/(n-1)' - average gap between consecutive nodes */
+
+
+ /* count how many levels the hierarchy contains (a SINGLE_LINK clustering */
+ /* alternatively we could use COMPLETE_LINK clustering) */
+ num_levels = 0;
+ for (i=1; i<n; i++) {
+ if (y[ordering[i]] - y[ordering[i-1]] > tol) {
+ num_levels++;
+ }
+ }
+ *num_levelsp = num_levels;
+ if (num_levels==0) {
+ *levelsp = levels = N_GNEW(1, int);
+ levels[0] = n;
+ }
+ else {
+ *levelsp = levels = N_GNEW(num_levels, int);
+ int count=0;
+ for (i=1; i<n; i++) {
+ if (y[ordering[i]] - y[ordering[i-1]] > tol) {
+ levels[count++] = i;
+ }
+ }
+ }
+ if (!given_coords) {
+ free(y);
+ }
+
+ return spread;
+}
+
+#endif /* DIGCOLA */
+
--- /dev/null
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#include "digcola.h"
+#ifdef DIGCOLA
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <stdio.h>
+#include <float.h>
+#include "stress.h"
+#include "dijkstra.h"
+#include "bfs.h"
+#include "matrix_ops.h"
+#include "kkutils.h"
+#include "conjgrad.h"
+#include "quad_prog_solver.h"
+#include "matrix_ops.h"
+
+#define localConstrMajorIterations 15
+#define levels_sep_tol 1e-1
+
+int
+stress_majorization_with_hierarchy(
+ vtx_data* graph, /* Input graph in sparse representation */
+ int n, /* Number of nodes */
+ int nedges_graph, /* Number of edges */
+ double** d_coords, /* Coordinates of nodes (output layout) */
+ int dim, /* Dimemsionality of layout */
+ int smart_ini, /* smart initialization */
+ int model, /* difference model */
+ int maxi, /* max iterations */
+ double levels_gap
+)
+{
+ int iterations = 0; /* Output: number of iteration of the process */
+
+ /*************************************************
+ ** Computation of full, dense, unrestricted k-D **
+ ** stress minimization by majorization **
+ ** This function imposes HIERARCHY CONSTRAINTS **
+ *************************************************/
+
+ int i,j,k;
+ bool directionalityExist = FALSE;
+ float * lap1 = NULL;
+ float * dist_accumulator = NULL;
+ float * tmp_coords = NULL;
+ float ** b = NULL;
+#ifdef NONCORE
+ FILE * fp=NULL;
+#endif
+ double * degrees = NULL;
+ float * lap2=NULL;
+ int lap_length;
+ float * f_storage=NULL;
+ float ** coords=NULL;
+
+ double conj_tol=tolerance_cg; /* tolerance of Conjugate Gradient */
+ float ** unpackedLap = NULL;
+ CMajEnv *cMajEnv = NULL;
+ clock_t start_time;
+ double y_0;
+ int length;
+ DistType diameter;
+ float * Dij=NULL;
+ /* to compensate noises, we never consider gaps smaller than 'abs_tol' */
+ double abs_tol=1e-2;
+ /* Additionally, we never consider gaps smaller than 'abs_tol'*<avg_gap> */
+ double relative_tol=levels_sep_tol;
+ int *ordering=NULL, *levels=NULL;
+ double hierarchy_spread;
+ float constant_term;
+ int count;
+ double degree;
+ int step;
+ float val;
+ double old_stress, new_stress;
+ bool converged;
+ int len;
+ int num_levels;
+ float *hierarchy_boundaries;
+
+ if (graph[0].edists!=NULL) {
+ for (i=0; i<n; i++) {
+ for (j=1; j<graph[i].nedges; j++) {
+ directionalityExist = directionalityExist || (graph[i].edists[j]!=0);
+ }
+ }
+ }
+ if (!directionalityExist) {
+ return stress_majorization_kD_mkernel(graph, n, nedges_graph, d_coords, dim, smart_ini, model, maxi);
+ }
+
+ /******************************************************************
+ ** First, partition nodes into layers: These are our constraints **
+ ******************************************************************/
+
+ if (smart_ini) {
+ double* x;
+ double* y;
+ if (dim>2) {
+ /* the dim==2 case is handled below */
+ stress_majorization_kD_mkernel(graph, n, nedges_graph, d_coords+1, dim-1, smart_ini, model, 15);
+ /* now copy the y-axis into the (dim-1)-axis */
+ for (i=0; i<n; i++) {
+ d_coords[dim-1][i] = d_coords[1][i];
+ }
+ }
+
+ x = d_coords[0]; y = d_coords[1];
+ compute_y_coords(graph, n, y, n);
+ hierarchy_spread = compute_hierarchy(graph, n, abs_tol, relative_tol, y, &ordering, &levels, &num_levels);
+ if (num_levels<=1) {
+ /* no hierarchy found, use faster algorithm */
+ return stress_majorization_kD_mkernel(graph, n, nedges_graph, d_coords, dim, smart_ini, model, maxi);
+ }
+
+ if (levels_gap>0) {
+ /* ensure that levels are separated in the initial layout */
+ double displacement = 0;
+ int stop;
+ for (i=0; i<num_levels; i++) {
+ displacement+=MAX((double)0,levels_gap-(y[ordering[levels[i]]]+displacement-y[ordering[levels[i]-1]]));
+ stop = i<num_levels-1 ? levels[i+1] : n;
+ for (j=levels[i]; j<stop; j++) {
+ y[ordering[j]] += displacement;
+ }
+ }
+ }
+ if (dim==2) {
+ IMDS_given_dim(graph, n, y, x, Epsilon);
+ }
+ }
+ else {
+ initLayout(graph, n, dim, d_coords);
+ hierarchy_spread = compute_hierarchy(graph, n, abs_tol, relative_tol, NULL, &ordering, &levels, &num_levels);
+ }
+
+ hierarchy_boundaries = N_GNEW(num_levels, float);
+
+ /****************************************************
+ ** Compute the all-pairs-shortest-distances matrix **
+ ****************************************************/
+
+ if (maxi==0)
+ return iterations;
+
+ if (model == MODEL_SUBSET) {
+ /* weight graph to separate high-degree nodes */
+ /* and perform slower Dijkstra-based computation */
+ if (Verbose)
+ fprintf(stderr, "Calculating subset model");
+ Dij = compute_apsp_artifical_weights_packed(graph, n);
+ } else if (model == MODEL_CIRCUIT) {
+ Dij = circuitModel(graph, n);
+ if (!Dij) {
+ agerr(AGWARN,
+ "graph is disconnected. Hence, the circuit model\n");
+ agerr(AGPREV,
+ "is undefined. Reverting to the shortest path model.\n");
+ }
+ }
+ if (!Dij) {
+ if (Verbose)
+ fprintf(stderr, "Calculating shortest paths");
+ Dij = compute_apsp_packed(graph, n);
+ }
+
+ diameter=-1;
+ length = n+n*(n-1)/2;
+ for (i=0; i<length; i++) {
+ if (Dij[i]>diameter) {
+ diameter = (int)Dij[i];
+ }
+ }
+
+ if (!smart_ini) {
+ /* for numerical stability, scale down layout */
+ /* No Jiggling, might conflict with constraints */
+ double max=1;
+ for (i=0; i<dim; i++) {
+ for (j=0; j<n; j++) {
+ if (fabs(d_coords[i][j])>max) {
+ max=fabs(d_coords[i][j]);
+ }
+ }
+ }
+ for (i=0; i<dim; i++) {
+ for (j=0; j<n; j++) {
+ d_coords[i][j]*=10/max;
+ }
+ }
+ }
+
+ if (levels_gap>0) {
+ int length = n+n*(n-1)/2;
+ double sum1, sum2, scale_ratio;
+ int count;
+ sum1=(float)(n*(n-1)/2);
+ sum2=0;
+ for (count=0, i=0; i<n-1; i++) {
+ count++; // skip self distance
+ for (j=i+1; j<n; j++,count++) {
+ sum2+=distance_kD(d_coords, dim, i, j)/Dij[count];
+ }
+ }
+ scale_ratio=sum2/sum1;
+ /* double scale_ratio=10; */
+ for (i=0; i<length; i++) {
+ Dij[i]*=(float)scale_ratio;
+ }
+ }
+
+ /**************************
+ ** Layout initialization **
+ **************************/
+
+ for (i=0; i<dim; i++) {
+ orthog1(n, d_coords[i]);
+ }
+
+ /* for the y-coords, don't center them, but translate them so y[0]=0 */
+ y_0 = d_coords[1][0];
+ for (i=0; i<n; i++) {
+ d_coords[1][i] -= y_0;
+ }
+
+ coords = N_GNEW(dim, float*);
+ f_storage = N_GNEW(dim*n, float);
+ for (i=0; i<dim; i++) {
+ coords[i] = f_storage+i*n;
+ for (j=0; j<n; j++) {
+ coords[i][j] = (float)(d_coords[i][j]);
+ }
+ }
+
+ /* compute constant term in stress sum
+ * which is \sum_{i<j} w_{ij}d_{ij}^2
+ */
+ constant_term=(float)(n*(n-1)/2);
+
+ /**************************
+ ** Laplacian computation **
+ **************************/
+
+ lap2 = Dij;
+ lap_length = n+n*(n-1)/2;
+ square_vec(lap_length, lap2);
+ /* compute off-diagonal entries */
+ invert_vec(lap_length, lap2);
+
+ /* compute diagonal entries */
+ count=0;
+ degrees = N_GNEW(n, double);
+ set_vector_val(n, 0, degrees);
+ for (i=0; i<n-1; i++) {
+ degree=0;
+ count++; // skip main diag entry
+ for (j=1; j<n-i; j++,count++) {
+ val = lap2[count];
+ degree+=val; degrees[i+j]-=val;
+ }
+ degrees[i]-=degree;
+ }
+ for (step=n,count=0,i=0; i<n; i++,count+=step,step--) {
+ lap2[count]=(float)degrees[i];
+ }
+
+#ifdef NONCORE
+ fpos_t pos;
+ if (n>max_nodes_in_mem) {
+ #define FILENAME "tmp_Dij$$$.bin"
+ fp = fopen(FILENAME, "wb");
+ fwrite(lap2, sizeof(float), lap_length, fp);
+ fclose(fp);
+ }
+#endif
+
+ /*************************
+ ** Layout optimization **
+ *************************/
+
+ b = N_GNEW (dim, float*);
+ b[0] = N_GNEW (dim*n, float);
+ for (k=1; k<dim; k++) {
+ b[k] = b[0]+k*n;
+ }
+
+ tmp_coords = N_GNEW(n, float);
+ dist_accumulator = N_GNEW(n, float);
+#ifdef NONCORE
+ if (n<=max_nodes_in_mem) {
+#endif
+ lap1 = N_GNEW(lap_length, float);
+#ifdef NONCORE
+ }
+ else {
+ lap1=lap2;
+ fp = fopen(FILENAME, "rb");
+ fgetpos(fp, &pos);
+ }
+#endif
+
+ old_stress=DBL_MAX; /* at least one iteration */
+
+ start_time = clock();
+
+ unpackedLap = unpackMatrix(lap2, n);
+ cMajEnv=initConstrainedMajorization(lap2, n, ordering, levels, num_levels);
+
+ for (converged=false,iterations=0; iterations<maxi && !converged; iterations++) {
+
+ /* First, construct Laplacian of 1/(d_ij*|p_i-p_j|) */
+ set_vector_val(n, 0, degrees);
+#ifdef NONCORE
+ if (n<=max_nodes_in_mem) {
+#endif
+ sqrt_vecf(lap_length, lap2, lap1);
+#ifdef NONCORE
+ }
+ else {
+ sqrt_vec(lap_length, lap1);
+ }
+#endif
+ for (count=0,i=0; i<n-1; i++) {
+ len = n-i-1;
+ /* init 'dist_accumulator' with zeros */
+ set_vector_valf(n, 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_valf(len, coords[k][i], tmp_coords);
+ vectors_mult_additionf(len, tmp_coords, -1, coords[k]+i+1);
+ square_vec(len, tmp_coords);
+ vectors_additionf(len, tmp_coords, dist_accumulator, dist_accumulator);
+ }
+
+ /* convert to 1/d_{ij} */
+ invert_sqrt_vec(len, dist_accumulator);
+ /* detect overflows */
+ for (j=0; j<len; j++) {
+ if (dist_accumulator[j]>=FLT_MAX || dist_accumulator[j]<0) {
+ dist_accumulator[j]=0;
+ }
+ }
+
+ count++; /* save place for the main diagonal entry */
+ degree=0;
+ for (j=0; j<len; j++,count++) {
+ val=lap1[count]*=dist_accumulator[j];
+ degree+=val; degrees[i+j+1]-=val;
+ }
+ degrees[i]-=degree;
+ }
+ for (step=n,count=0,i=0; i<n; i++,count+=step,step--) {
+ lap1[count]=(float)degrees[i];
+ }
+
+ /* Now compute b[] (L^(X(t))*X(t)) */
+ for (k=0; k<dim; k++) {
+ /* b[k] := lap1*coords[k] */
+ right_mult_with_vector_ff(lap1, n, coords[k], b[k]);
+ }
+
+ /* compute new stress
+ * remember that the Laplacians are negated, so we subtract
+ * instead of add and vice versa
+ */
+ new_stress=0;
+ for (k=0; k<dim; k++) {
+ new_stress+=vectors_inner_productf(n, coords[k], b[k]);
+ }
+ new_stress*=2;
+ new_stress+=constant_term; // only after mult by 2
+#ifdef NONCORE
+ if (n>max_nodes_in_mem) {
+ /* restore lap2 from disk */
+ fsetpos(fp, &pos);
+ fread(lap2, sizeof(float), lap_length, fp);
+ }
+#endif
+ for (k=0; k<dim; k++) {
+ right_mult_with_vector_ff(lap2, n, coords[k], tmp_coords);
+ new_stress-=vectors_inner_productf(n, coords[k], tmp_coords);
+ }
+
+#ifdef ALTERNATIVE_STRESS_CALC
+ {
+ double mat_stress=new_stress;
+ double compute_stress(float ** coords, float * lap, int dim, int n);
+ new_stress = compute_stress(coords, lap2, dim, n);
+ if (fabs(mat_stress-new_stress)/min(mat_stress,new_stress)>0.001) {
+ fprintf(stderr,"Diff stress vals: %lf %lf (iteration #%d)\n", mat_stress, new_stress,iterations);
+ }
+ }
+#endif
+ /* check for convergence */
+ converged = fabs(new_stress-old_stress)/fabs(old_stress+1e-10) < Epsilon;
+ converged = converged || (iterations>1 && new_stress>old_stress);
+ /* in first iteration we allowed stress increase, which
+ * might result ny imposing constraints
+ */
+ old_stress = new_stress;
+
+ for (k=0; k<dim; k++) {
+ /* now we find the optimizer of trace(X'LX)+X'B by solving 'dim'
+ * system of equations, thereby obtaining the new coordinates.
+ * If we use the constraints (given by the var's: 'ordering',
+ * 'levels' and 'num_levels'), we cannot optimize
+ * trace(X'LX)+X'B by simply solving equations, but we have
+ * to use a quadratic programming solver
+ * note: 'lap2' is a packed symmetric matrix, that is its
+ * upper-triangular part is arranged in a vector row-wise
+ * also note: 'lap2' is really the negated laplacian (the
+ * laplacian is -'lap2')
+ */
+
+ if (k==1) {
+ /* use quad solver in the y-dimension */
+ constrained_majorization_new_with_gaps(cMajEnv, b[k], coords, dim, k, localConstrMajorIterations, hierarchy_boundaries, levels_gap);
+
+ }
+ else {
+ /* use conjugate gradient for all dimensions except y */
+ conjugate_gradient_mkernel(lap2, coords[k], b[k], n, conj_tol, n);
+ }
+ }
+ }
+ free (hierarchy_boundaries);
+ deleteCMajEnv(cMajEnv);
+
+ if (coords!=NULL) {
+ for (i=0; i<dim; i++) {
+ for (j=0; j<n; j++) {
+ d_coords[i][j] = coords[i][j];
+ }
+ }
+ free (coords[0]); free (coords);
+ }
+
+ if (b) {
+ free (b[0]); free (b);
+ }
+ free (tmp_coords);
+ free (dist_accumulator);
+ free (degrees);
+ free (lap2);
+
+
+#ifdef NONCORE
+ if (n<=max_nodes_in_mem) {
+#endif
+ free (lap1);
+#ifdef NONCORE
+ }
+#endif
+
+ free (ordering);
+
+ free (levels);
+
+ if (unpackedLap) {
+ free (unpackedLap[0]); free (unpackedLap);
+ }
+ return iterations;
+}
+#endif /* DIGCOLA */
+
--- /dev/null
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#include "digcola.h"
+#ifdef DIGCOLA
+#include "matrix_ops.h"
+#include "conjgrad.h"
+
+static void
+construct_b(vtx_data* graph, int n, double *b)
+{
+ /* construct a vector - b s.t. -b[i]=\sum_j -w_{ij}*\delta_{ij}
+ * (the "balance vector")
+ * Note that we build -b and not b, since our matrix is not the
+ * real laplacian L, but its negation: -L.
+ * So instead of solving Lx=b, we will solve -Lx=-b
+ */
+ int i,j;
+
+ double b_i = 0;
+
+ for (i=0; i<n; i++) {
+ b_i = 0;
+ if (graph[0].edists==NULL) {
+ continue;
+ }
+ for (j=1; j<graph[i].nedges; j++) { /* skip the self loop */
+ b_i += graph[i].ewgts[j]*graph[i].edists[j];
+ }
+ b[i] = b_i;
+ }
+}
+
+#define hierarchy_cg_tol 1e-3
+
+void
+compute_y_coords(vtx_data* graph, int n, double* y_coords, int max_iterations)
+{
+ /* Find y coords of a directed graph by solving L*x = b */
+ int i,j;
+ double* b = N_NEW(n, double);
+ double tol = hierarchy_cg_tol;
+ int nedges = 0;
+ float* uniform_weights;
+ float* old_ewgts = graph[0].ewgts;
+
+ construct_b(graph, n, b);
+
+ init_vec_orth1(n, y_coords);
+
+ for (i=0; i<n; i++) {
+ nedges += graph[i].nedges;
+ }
+
+ /* replace original edge weights (which are lengths) with uniform weights */
+ /* for computing the optimal arrangement */
+ uniform_weights = N_GNEW(nedges,float);
+ for (i=0; i<n; i++) {
+ graph[i].ewgts = uniform_weights;
+ uniform_weights[0] = (float)-(graph[i].nedges-1);
+ for (j=1; j<graph[i].nedges; j++) {
+ uniform_weights[j] = 1;
+ }
+ uniform_weights += graph[i].nedges;
+ }
+
+ conjugate_gradient(graph, y_coords, b, n, tol, max_iterations);
+
+ /* restore original edge weights */
+ free (graph[0].ewgts);
+ for (i=0; i<n; i++) {
+ graph[i].ewgts = old_ewgts;
+ old_ewgts += graph[i].nedges;
+ }
+
+ free(b);
+}
+
+#endif /* DIGCOLA */
+
--- /dev/null
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#include "digcola.h"
+#ifdef DIGCOLA
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <stdio.h>
+#include <float.h>
+#include <assert.h>
+#include "matrix_ops.h"
+#include "kkutils.h"
+#include "quad_prog_solver.h"
+
+#define quad_prog_tol 1e-2
+
+float**
+unpackMatrix(float * packedMat, int n)
+{
+ float **mat;
+ int i,j,k;
+
+ mat = N_GNEW(n, float*);
+ mat[0] = N_GNEW(n*n, float);
+ set_vector_valf(n*n, 0, mat[0]);
+ for (i=1; i<n; i++) {
+ mat[i] = mat[0]+i*n;
+ }
+
+ for (i=0,k=0; i<n; i++) {
+ for (j=i; j<n; j++,k++) {
+ mat[j][i]=mat[i][j]=packedMat[k];
+ }
+ }
+ return mat;
+}
+
+static void
+ensureMonotonicOrdering(float * place, int n, int * ordering)
+{
+ /* ensure that 'ordering' is monotonically increasing by 'place', */
+ /* this also implies that levels are separated in the initial layout */
+ int i, node;
+ float lower_bound = place[ordering[0]];
+ for (i=1; i<n; i++) {
+ node = ordering[i];
+ if (place[node]<lower_bound) {
+ place[node]=lower_bound;
+ }
+ lower_bound = place[node];
+ }
+}
+
+static void
+ensureMonotonicOrderingWithGaps(float* place, int n, int* ordering,
+ int* levels, int num_levels, float levels_gap)
+{
+ /* ensure that levels are separated in the initial layout and that
+ * places are monotonic increasing within layer
+ */
+
+ int i;
+ int node, level, max_in_level;
+ float lower_bound;
+
+ level = -1; max_in_level=0;
+ for (i=0; i<n; i++) {
+ if (i>=max_in_level) {
+ /* we are entering a new level */
+ level++;
+ if (level==num_levels) {
+ /* last_level */
+ max_in_level=n;
+ }
+ else {
+ max_in_level = levels[level];
+ }
+ lower_bound = i>0 ? place[ordering[i-1]]+levels_gap : (float)-1e9;
+ quicksort_placef(place, ordering, i, max_in_level-1);
+ }
+
+ node = ordering[i];
+ if (place[node]<lower_bound) {
+ place[node]=lower_bound;
+ }
+ }
+}
+
+static void
+computeHierarchyBoundaries(float* place, int n, int* ordering, int* levels,
+ int num_levels, float* hierarchy_boundaries)
+{
+ int i;
+ for (i=0; i<num_levels; i++) {
+ hierarchy_boundaries[i] = place[ordering[levels[i]-1]];
+ }
+}
+
+
+int
+constrained_majorization_new(CMajEnv *e, float * b, float **coords,
+ int cur_axis, int dims, int max_iterations,
+ float* hierarchy_boundaries, float levels_gap)
+{
+ int n=e->n;
+ float *place=coords[cur_axis];
+ float **lap=e->A;
+ int * ordering=e->ordering;
+ int * levels=e->levels;
+ int num_levels=e->num_levels;
+ int i,j;
+ float new_place_i;
+ bool converged=false;
+ float upper_bound, lower_bound;
+ int node;
+ int left, right;
+ float cur_place;
+ float des_place_block;
+ float block_deg;
+ float toBlockConnectivity;
+ float * lap_node;
+ int block_len;
+ int first_next_level;
+ int level=-1, max_in_level=0;
+ float* desired_place;
+ float* prefix_desired_place;
+ float* suffix_desired_place;
+ int* block;
+ int* lev;
+ int counter;
+
+ if (max_iterations<=0) {
+ return 0;
+ }
+ if (levels_gap!=0) {
+ return constrained_majorization_new_with_gaps(e, b, coords,cur_axis,dims, max_iterations, hierarchy_boundaries, levels_gap);
+ }
+
+ /* ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels); */
+ ensureMonotonicOrdering(place, n, ordering);
+ /* it is important that in 'ordering' nodes are always sorted by layers,
+ * and within a layer by 'place'
+ */
+
+ /* the desired place of each individual node in the current block */
+ desired_place = e->fArray1;
+ /* the desired place of each prefix of current block */
+ prefix_desired_place = e->fArray2;
+ /* the desired place of each suffix of current block */
+ suffix_desired_place = e->fArray3;
+ /* current block (nodes connected by active constraints) */
+ block = e->iArray1;
+
+ lev = e->iArray2; /* level of each node */
+ for (i=0; i<n; i++) {
+ if (i>=max_in_level) {
+ /* we are entering a new level */
+ level++;
+ if (level==num_levels) {
+ /* last_level */
+ max_in_level=n;
+ }
+ else {
+ max_in_level = levels[level];
+ }
+ }
+ node = ordering[i];
+ lev[node]=level;
+ }
+
+ for (counter=0; counter<max_iterations&&!converged; counter++) {
+ converged=true;
+ lower_bound = -1e9; /* no lower bound for first level */
+ for (left=0; left<n; left=right) {
+ int best_i;
+ double max_movement;
+ double movement;
+ float prefix_des_place, suffix_des_place;
+ /* compute a block 'ordering[left]...ordering[right-1]' of
+ * nodes with the same coordinate:
+ */
+ cur_place = place[ordering[left]];
+ for (right=left+1; right<n; right++) {
+ if (place[ordering[right]]!=cur_place) {
+ break;
+ }
+ }
+
+ /* compute desired place of nodes in block: */
+ for (i=left; i<right; i++) {
+ node = ordering[i];
+ new_place_i = -b[node];
+ lap_node = lap[node];
+ for (j=0; j<n; j++) {
+ if (j==node) {
+ continue;
+ }
+ new_place_i += lap_node[j]*place[j];
+ }
+ desired_place[node] = new_place_i/(-lap_node[node]);
+ }
+
+ /* reorder block by levels, and within levels by "relaxed" desired position */
+ block_len = 0;
+ first_next_level=0;
+ for (i=left; i<right; i=first_next_level) {
+ level = lev[ordering[i]];
+ if (level==num_levels) {
+ /* last_level */
+ first_next_level=right;
+ }
+ else {
+ first_next_level = MIN(right,levels[level]);
+ }
+
+ /* First, collect all nodes with desired places smaller than current place */
+ for (j=i; j<first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] < cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Second, collect all nodes with desired places equal to current place */
+ for (j=i; j<first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] == cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Third, collect all nodes with desired places greater than current place */
+ for (j=i; j<first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] > cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ }
+
+ /* loop through block and compute desired places of its prefixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i=0; i<block_len; i++) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j=0; j<i; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity*=2;
+ /* update block stats */
+ des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
+ prefix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity-lap_node[node];
+ }
+
+ /* loop through block and compute desired places of its suffixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i=block_len-1; i>=0; i--) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j=i+1; j<block_len; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity*=2;
+ /* update block stats */
+ des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
+ suffix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity-lap_node[node];
+ }
+
+
+ /* now, find best place to split block */
+ best_i = -1;
+ max_movement = 0;
+ for (i=0; i<block_len; i++) {
+ suffix_des_place = suffix_desired_place[i];
+ prefix_des_place = i>0 ? prefix_desired_place[i-1] : suffix_des_place;
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place<prefix_des_place) {
+ if (suffix_des_place<cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
+ }
+ suffix_des_place = prefix_des_place;
+ }
+ else if (prefix_des_place>cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
+ movement = (block_len-i)*fabs(suffix_des_place-cur_place) + i*fabs(prefix_des_place-cur_place);
+ if (movement>max_movement) {
+ max_movement = movement;
+ best_i = i;
+ }
+ }
+ /* Actually move prefix and suffix */
+ if (best_i>=0) {
+ suffix_des_place = suffix_desired_place[best_i];
+ prefix_des_place = best_i>0 ? prefix_desired_place[best_i-1] : suffix_des_place;
+
+ /* compute right border of feasible move */
+ if (right>=n) {
+ /* no nodes after current block */
+ upper_bound = 1e9;
+ }
+ else {
+ upper_bound = place[ordering[right]];
+ }
+ suffix_des_place = MIN(suffix_des_place, upper_bound);
+ prefix_des_place = MAX(prefix_des_place, lower_bound);
+
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place<prefix_des_place) {
+ if (suffix_des_place<cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
+ }
+ suffix_des_place = prefix_des_place;
+ }
+ else if (prefix_des_place>cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
+
+ /* move prefix: */
+ for (i=0; i<best_i; i++) {
+ place[block[i]] = prefix_des_place;
+ }
+ /* move suffix: */
+ for (i=best_i; i<block_len; i++) {
+ place[block[i]] = suffix_des_place;
+ }
+
+ lower_bound = suffix_des_place; // lower bound for next block
+
+ /* reorder 'ordering' to reflect change of places
+ * Note that it is enough to reorder the level where
+ * the split was done
+ */
+#if 0
+ int max_in_level, min_in_level;
+
+ level = lev[best_i];
+ if (level==num_levels) {
+ // last_level
+ max_in_level = MIN(right,n);
+ }
+ else {
+ max_in_level = MIN(right,levels[level]);
+ }
+ if (level==0) {
+ // first level
+ min_in_level = MAX(left,0);
+ }
+ else {
+ min_in_level = MAX(left,levels[level-1]);
+ }
+#endif
+ for (i=left; i<right; i++) {
+ ordering[i] = block[i-left];
+ }
+ converged = converged && fabs(prefix_des_place-cur_place)<quad_prog_tol
+ && fabs(suffix_des_place-cur_place)<quad_prog_tol;
+
+ }
+ else {
+ /* no movement */
+ lower_bound = cur_place; // lower bound for next block
+ }
+
+ }
+ }
+
+ computeHierarchyBoundaries(place, n, ordering, levels, num_levels, hierarchy_boundaries);
+
+ return counter;
+}
+
+int
+constrained_majorization_new_with_gaps(CMajEnv *e, float* b, float** coords,
+ int ndims, int cur_axis, int max_iterations, float* hierarchy_boundaries,
+ float levels_gap)
+{
+ float *place=coords[cur_axis];
+ int i,j;
+ int n=e->n;
+ float** lap = e->A;
+ int* ordering = e->ordering;
+ int* levels = e->levels;
+ int num_levels = e->num_levels;
+ float new_place_i;
+ bool converged=false;
+ float upper_bound, lower_bound;
+ int node;
+ int left, right;
+ float cur_place;
+ float des_place_block;
+ float block_deg;
+ float toBlockConnectivity;
+ float * lap_node;
+ float * desired_place;
+ float * prefix_desired_place;
+ float * suffix_desired_place;
+ int* block;
+ int block_len;
+ int first_next_level;
+ int * lev;
+ int level=-1, max_in_level=0;
+ int counter;
+ float* gap;
+ float total_gap, target_place;
+
+ if (max_iterations<=0) {
+ return 0;
+ }
+
+ ensureMonotonicOrderingWithGaps(place, n, ordering, levels, num_levels, levels_gap);
+ /* it is important that in 'ordering' nodes are always sorted by layers,
+ * and within a layer by 'place'
+ */
+
+ /* the desired place of each individual node in the current block */
+ desired_place = e->fArray1;
+ /* the desired place of each prefix of current block */
+ prefix_desired_place = e->fArray2;
+ /* the desired place of each suffix of current block */
+ suffix_desired_place = e->fArray3;
+ /* current block (nodes connected by active constraints) */
+ block = e->iArray1;
+
+ lev = e->iArray2; /* level of each node */
+ for (i=0; i<n; i++) {
+ if (i>=max_in_level) {
+ /* we are entering a new level */
+ level++;
+ if (level==num_levels) {
+ /* last_level */
+ max_in_level=n;
+ }
+ else {
+ max_in_level = levels[level];
+ }
+ }
+ node = ordering[i];
+ lev[node]=level;
+ }
+
+ /* displacement of block's nodes from block's reference point */
+ gap = e->fArray4;
+
+ for (counter=0; counter<max_iterations&&!converged; counter++) {
+ converged=true;
+ lower_bound = -1e9; // no lower bound for first level
+ for (left=0; left<n; left=right) {
+ int best_i;
+ double max_movement;
+ double movement;
+ float prefix_des_place, suffix_des_place;
+ /* compute a block 'ordering[left]...ordering[right-1]' of
+ * nodes connected with active constraints
+ */
+ cur_place = place[ordering[left]];
+ total_gap = 0;
+ target_place = cur_place;
+ gap[ordering[left]] = 0;
+ for (right=left+1; right<n; right++) {
+ if (lev[right]>lev[right-1]) {
+ /* we are entering a new level */
+ target_place += levels_gap; /* note that 'levels_gap' may be negative */
+ total_gap += levels_gap;
+ }
+ node = ordering[right];
+#if 0
+ if (place[node]!=target_place)
+#endif
+ /* not sure if this is better than 'place[node]!=target_place' */
+ if (fabs(place[node]-target_place)>1e-9)
+ {
+ break;
+ }
+ gap[node] = place[node]-cur_place;
+ }
+
+ /* compute desired place of block's reference point according
+ * to each node in the block:
+ */
+ for (i=left; i<right; i++) {
+ node = ordering[i];
+ new_place_i = -b[node];
+ lap_node = lap[node];
+ for (j=0; j<n; j++) {
+ if (j==node) {
+ continue;
+ }
+ new_place_i += lap_node[j]*place[j];
+ }
+ desired_place[node] = new_place_i/(-lap_node[node]) - gap[node];
+ }
+
+ /* reorder block by levels, and within levels
+ * by "relaxed" desired position
+ */
+ block_len = 0;
+ first_next_level=0;
+ for (i=left; i<right; i=first_next_level) {
+ level = lev[ordering[i]];
+ if (level==num_levels) {
+ /* last_level */
+ first_next_level=right;
+ }
+ else {
+ first_next_level = MIN(right,levels[level]);
+ }
+
+ /* First, collect all nodes with desired places smaller
+ * than current place
+ */
+ for (j=i; j<first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] < cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Second, collect all nodes with desired places equal
+ * to current place
+ */
+ for (j=i; j<first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] == cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ /* Third, collect all nodes with desired places greater
+ * than current place
+ */
+ for (j=i; j<first_next_level; j++) {
+ node = ordering[j];
+ if (desired_place[node] > cur_place) {
+ block[block_len++] = node;
+ }
+ }
+ }
+
+ /* loop through block and compute desired places of its prefixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i=0; i<block_len; i++) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j=0; j<i; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity*=2;
+ /* update block stats */
+ des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
+ prefix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity-lap_node[node];
+ }
+
+ if (block_len==n) {
+ /* fix is needed since denominator was 0 in this case */
+ prefix_desired_place[n-1] = cur_place; /* a "neutral" value */
+ }
+
+ /* loop through block and compute desired places of its suffixes */
+ des_place_block = 0;
+ block_deg = 0;
+ for (i=block_len-1; i>=0; i--) {
+ node = block[i];
+ toBlockConnectivity = 0;
+ lap_node = lap[node];
+ for (j=i+1; j<block_len; j++) {
+ toBlockConnectivity -= lap_node[block[j]];
+ }
+ toBlockConnectivity*=2;
+ /* update block stats */
+ des_place_block = (block_deg*des_place_block+(-lap_node[node])*desired_place[node]+toBlockConnectivity*cur_place)/(block_deg-lap_node[node]+toBlockConnectivity);
+ suffix_desired_place[i] = des_place_block;
+ block_deg += toBlockConnectivity-lap_node[node];
+ }
+
+ if (block_len==n) {
+ /* fix is needed since denominator was 0 in this case */
+ suffix_desired_place[0] = cur_place; /* a "neutral" value? */
+ }
+
+
+ /* now, find best place to split block */
+ best_i = -1;
+ max_movement = 0;
+ for (i=0; i<block_len; i++) {
+ suffix_des_place = suffix_desired_place[i];
+ prefix_des_place = i>0 ? prefix_desired_place[i-1] : suffix_des_place;
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place<prefix_des_place) {
+ if (suffix_des_place<cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
+ }
+ suffix_des_place = prefix_des_place;
+ }
+ else if (prefix_des_place>cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
+ movement = (block_len-i)*fabs(suffix_des_place-cur_place) + i*fabs(prefix_des_place-cur_place);
+ if (movement>max_movement) {
+ max_movement = movement;
+ best_i = i;
+ }
+ }
+ /* Actually move prefix and suffix */
+ if (best_i>=0) {
+ suffix_des_place = suffix_desired_place[best_i];
+ prefix_des_place = best_i>0 ? prefix_desired_place[best_i-1] : suffix_des_place;
+
+ /* compute right border of feasible move */
+ if (right>=n) {
+ /* no nodes after current block */
+ upper_bound = 1e9;
+ }
+ else {
+ /* notice that we have to deduct 'gap[block[block_len-1]]'
+ * since all computations should be relative to
+ * the block's reference point
+ */
+ if (lev[ordering[right]]>lev[ordering[right-1]]) {
+ upper_bound = place[ordering[right]]-levels_gap-gap[block[block_len-1]];
+ }
+ else {
+ upper_bound = place[ordering[right]]-gap[block[block_len-1]];
+ }
+ }
+ suffix_des_place = MIN(suffix_des_place, upper_bound);
+ prefix_des_place = MAX(prefix_des_place, lower_bound);
+
+ /* limit moves to ensure that the prefix is placed before the suffix */
+ if (suffix_des_place<prefix_des_place) {
+ if (suffix_des_place<cur_place) {
+ if (prefix_des_place > cur_place) {
+ prefix_des_place = cur_place;
+ }
+ suffix_des_place = prefix_des_place;
+ }
+ else if (prefix_des_place>cur_place) {
+ prefix_des_place = suffix_des_place;
+ }
+ }
+
+ /* move prefix: */
+ for (i=0; i<best_i; i++) {
+ place[block[i]] = prefix_des_place+gap[block[i]];
+ }
+ /* move suffix: */
+ for (i=best_i; i<block_len; i++) {
+ place[block[i]] = suffix_des_place+gap[block[i]];
+ }
+
+
+ /* compute lower bound for next block */
+ if (right<n && lev[ordering[right]]>lev[ordering[right-1]]) {
+ lower_bound = place[block[block_len-1]]+levels_gap;
+ }
+ else {
+ lower_bound = place[block[block_len-1]];
+ }
+
+
+ /* reorder 'ordering' to reflect change of places
+ * Note that it is enough to reorder the level where
+ * the split was done
+ */
+#if 0
+ int max_in_level, min_in_level;
+
+ level = lev[best_i];
+ if (level==num_levels) {
+ /* last_level */
+ max_in_level = MIN(right,n);
+ }
+ else {
+ max_in_level = MIN(right,levels[level]);
+ }
+ if (level==0) {
+ /* first level */
+ min_in_level = MAX(left,0);
+ }
+ else {
+ min_in_level = MAX(left,levels[level-1]);
+ }
+#endif
+ for (i=left; i<right; i++) {
+ ordering[i] = block[i-left];
+ }
+ converged = converged && fabs(prefix_des_place-cur_place)<quad_prog_tol
+ && fabs(suffix_des_place-cur_place)<quad_prog_tol;
+
+
+ }
+ else {
+ /* no movement */
+ /* compute lower bound for next block */
+ if (right<n && lev[ordering[right]]>lev[ordering[right-1]]) {
+ lower_bound = place[block[block_len-1]]+levels_gap;
+ }
+ else {
+ lower_bound = place[block[block_len-1]];
+ }
+ }
+ }
+ orthog1f(n, place); /* for numerical stability, keep ||place|| small */
+ computeHierarchyBoundaries(place, n, ordering, levels, num_levels, hierarchy_boundaries);
+ }
+
+ return counter;
+}
+
+void
+deleteCMajEnv(CMajEnv *e)
+{
+ free (e->A[0]);
+ free (e->A);
+ free (e->lev);
+ free (e->fArray1);
+ free (e->fArray2);
+ free (e->fArray3);
+ free (e->fArray4);
+ free (e->iArray1);
+ free (e->iArray2);
+ free (e->iArray3);
+ free (e->iArray4);
+ free (e);
+}
+
+CMajEnv*
+initConstrainedMajorization(float* packedMat, int n, int* ordering,
+ int* levels, int num_levels)
+{
+ int i, level=-1, start_of_level_above=0;
+ CMajEnv *e = GNEW(CMajEnv);
+ e->A=NULL;
+ e->n=n;
+ e->ordering=ordering;
+ e->levels=levels;
+ e->num_levels=num_levels;
+ e->A = unpackMatrix(packedMat,n);
+ e->lev = N_GNEW(n, int);
+ for (i=0; i<e->n; i++) {
+ if (i>=start_of_level_above) {
+ level++;
+ start_of_level_above=(level==num_levels)?e->n:levels[level];
+ }
+ e->lev[ordering[i]]=level;
+ }
+ e->fArray1 = N_GNEW(n,float);
+ e->fArray2 = N_GNEW(n,float);
+ e->fArray3 = N_GNEW(n,float);
+ e->fArray4 = N_GNEW(n,float);
+ e->iArray1 = N_GNEW(n,int);
+ e->iArray2 = N_GNEW(n,int);
+ e->iArray3 = N_GNEW(n,int);
+ e->iArray4 = N_GNEW(n,int);
+ return e;
+}
+#endif /* DIGCOLA */
+
--- /dev/null
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#ifndef _CMAJ_H_
+#define _CMAJ_H_
+
+#ifdef DIGCOLA
+
+typedef struct {
+ float **A;
+ int n;
+ int *lev;
+ int *iArray1;
+ int *iArray2;
+ int *iArray3;
+ int *iArray4;
+ float *fArray1;
+ float *fArray2;
+ float *fArray3;
+ float *fArray4;
+ float *A_r;
+ int *ordering;
+ int *levels;
+ int num_levels;
+}CMajEnv;
+
+extern CMajEnv* initConstrainedMajorization(float *, int, int*, int*, int);
+
+extern int constrained_majorization_new(CMajEnv*, float*, float**,
+ int, int, int, float*, float);
+
+extern int constrained_majorization_new_with_gaps(CMajEnv*, float*, float**,
+ int, int, int, float*, float);
+extern void deleteCMajEnv(CMajEnv *e);
+
+extern float** unpackMatrix(float * packedMat, int n);
+
+#endif
+
+#endif /* _CMAJ_H_ */
+
+#ifdef __cplusplus
+}
+#endif
--- /dev/null
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#include "digcola.h"
+#ifdef DIGCOLA
+#include "kkutils.h"
+#include "matrix_ops.h"
+#include "conjgrad.h"
+
+static void
+standardize(double* orthog, int nvtxs)
+{
+ double len, avg = 0;
+ int i;
+ for (i=0; i<nvtxs; i++)
+ avg+=orthog[i];
+ avg/=nvtxs;
+
+ /* centralize: */
+ for (i=0; i<nvtxs; i++)
+ orthog[i]-=avg;
+
+ /* normalize: */
+ len = norm(orthog, 0, nvtxs-1);
+ vecscale(orthog, 0, nvtxs-1, 1.0 / len, orthog);
+}
+
+static void
+mat_mult_vec_orthog(float** mat, int dim1, int dim2, double* vec,
+ double* result, double* orthog)
+{
+ /* computes mat*vec, where mat is a dim1*dim2 matrix */
+ int i,j;
+ double sum;
+
+ for (i=0; i<dim1; i++) {
+ sum=0;
+ for (j=0; j<dim2; j++) {
+ sum += mat[i][j]*vec[j];
+ }
+ result[i]=sum;
+ }
+ if (orthog!=NULL) {
+ double alpha=-dot(result,0,dim1-1,orthog);
+ scadd(result, 0, dim1-1, alpha, orthog);
+ }
+}
+
+static void
+power_iteration_orthog(float** square_mat, int n, int neigs,
+ double** eigs, double* evals, double* orthog, double p_iteration_threshold)
+{
+ /*
+ * Power-Iteration with (I-orthog*orthog^T)*square_mat*(I-orthog*orthog^T)
+ */
+
+ int i,j;
+ double *tmp_vec = N_GNEW(n, double);
+ double *last_vec = N_GNEW(n, double);
+ double *curr_vector;
+ double len;
+ double angle;
+ double alpha;
+ int iteration;
+ int largest_index;
+ double largest_eval;
+
+ double tol=1-p_iteration_threshold;
+
+ if (neigs>=n) {
+ neigs=n;
+ }
+
+ for (i=0; i<neigs; i++) {
+ curr_vector = eigs[i];
+ /* guess the i-th eigen vector */
+choose:
+ for (j=0; j<n; j++) {
+ curr_vector[j] = rand()%100;
+ }
+
+ if (orthog!=NULL) {
+ alpha=-dot(orthog,0,n-1,curr_vector);
+ scadd(curr_vector, 0, n-1, alpha, orthog);
+ }
+ // orthogonalize against higher eigenvectors
+ for (j=0; j<i; j++) {
+ alpha = -dot(eigs[j], 0, n-1, curr_vector);
+ scadd(curr_vector, 0, n-1, alpha, eigs[j]);
+ }
+ len = norm(curr_vector, 0, n-1);
+ if (len<1e-10) {
+ /* We have chosen a vector colinear with prvious ones */
+ goto choose;
+ }
+ vecscale(curr_vector, 0, n-1, 1.0 / len, curr_vector);
+ iteration=0;
+ do {
+ iteration++;
+ cpvec(last_vec,0,n-1,curr_vector);
+
+ mat_mult_vec_orthog(square_mat,n,n,curr_vector,tmp_vec,orthog);
+ cpvec(curr_vector,0,n-1,tmp_vec);
+
+ /* orthogonalize against higher eigenvectors */
+ for (j=0; j<i; j++) {
+ alpha = -dot(eigs[j], 0, n-1, curr_vector);
+ scadd(curr_vector, 0, n-1, alpha, eigs[j]);
+ }
+ len = norm(curr_vector, 0, n-1);
+ if (len<1e-10) {
+ /* We have reached the null space (e.vec. associated
+ * with e.val. 0)
+ */
+ goto exit;
+ }
+
+ vecscale(curr_vector, 0, n-1, 1.0 / len, curr_vector);
+ angle = dot(curr_vector, 0, n-1, last_vec);
+ } while (fabs(angle)<tol);
+ /* the Rayleigh quotient (up to errors due to orthogonalization):
+ * u*(A*u)/||A*u||)*||A*u||, where u=last_vec, and ||u||=1
+ */
+ evals[i]=angle*len;
+ }
+exit:
+ for (; i<neigs; i++) {
+ /* compute the smallest eigenvector, which are
+ * probably associated with eigenvalue 0 and for
+ * which power-iteration is dangerous
+ */
+ curr_vector = eigs[i];
+ /* guess the i-th eigen vector */
+ for (j=0; j<n; j++)
+ curr_vector[j] = rand()%100;
+ /* orthogonalize against higher eigenvectors */
+ for (j=0; j<i; j++) {
+ alpha = -dot(eigs[j], 0, n-1, curr_vector);
+ scadd(curr_vector, 0, n-1, alpha, eigs[j]);
+ }
+ len = norm(curr_vector, 0, n-1);
+ vecscale(curr_vector, 0, n-1, 1.0 / len, curr_vector);
+ evals[i]=0;
+
+ }
+
+ /* sort vectors by their evals, for overcoming possible mis-convergence: */
+ for (i=0; i<neigs-1; i++) {
+ largest_index=i;
+ largest_eval=evals[largest_index];
+ for (j=i+1; j<neigs; j++) {
+ if (largest_eval<evals[j]) {
+ largest_index=j;
+ largest_eval=evals[largest_index];
+ }
+ }
+ if (largest_index!=i) { // exchange eigenvectors:
+ cpvec(tmp_vec,0,n-1,eigs[i]);
+ cpvec(eigs[i],0,n-1,eigs[largest_index]);
+ cpvec(eigs[largest_index],0,n-1,tmp_vec);
+
+ evals[largest_index]=evals[i];
+ evals[i]=largest_eval;
+ }
+ }
+
+ free (tmp_vec); free (last_vec);
+
+}
+
+static float*
+compute_avgs(DistType** Dij, int n, float* all_avg)
+{
+ float* row_avg = N_GNEW(n, float);
+ int i,j;
+ double sum=0, sum_row;
+
+ for (i=0; i<n; i++) {
+ sum_row=0;
+ for (j=0; j<n; j++) {
+ sum+=(double)Dij[i][j]*(double)Dij[i][j];
+ sum_row+=(double)Dij[i][j]*(double)Dij[i][j];
+ }
+ row_avg[i]=(float)sum_row/n;
+ }
+ *all_avg=(float)sum/(n*n);
+ return row_avg;
+}
+
+static float**
+compute_Bij(DistType** Dij, int n)
+{
+ int i,j;
+ float* storage = N_GNEW(n*n,float);
+ float** Bij = N_GNEW(n, float*);
+ float* row_avg;
+ float all_avg;
+
+ for (i=0; i<n; i++)
+ Bij[i] = storage+i*n;
+
+ row_avg = compute_avgs(Dij, n, &all_avg);
+ for (i=0; i<n; i++) {
+ for (j=0; j<=i; j++) {
+ Bij[i][j]=-(float)Dij[i][j]*Dij[i][j]+row_avg[i]+row_avg[j]-all_avg;
+ Bij[j][i]=Bij[i][j];
+ }
+ }
+ free (row_avg);
+ return Bij;
+}
+
+static void
+CMDS_orthog(vtx_data* graph, int n, int dim, double** eigs, double tol,
+ double* orthog, DistType** Dij)
+{
+ int i,j;
+ float** Bij = compute_Bij(Dij, n);
+ double* evals= N_GNEW(dim, double);
+
+ double * orthog_aux = NULL;
+ if (orthog!=NULL) {
+ orthog_aux = N_GNEW(n, double);
+ for (i=0; i<n; i++) {
+ orthog_aux[i]=orthog[i];
+ }
+ standardize(orthog_aux,n);
+ }
+ power_iteration_orthog(Bij, n, dim, eigs, evals, orthog_aux, tol);
+
+ for (i=0; i<dim; i++) {
+ for (j=0; j<n; j++) {
+ eigs[i][j]*=sqrt(fabs(evals[i]));
+ }
+ }
+ free (Bij[0]); free (Bij);
+ free (evals); free (orthog_aux);
+}
+
+#define SCALE_FACTOR 256
+
+void IMDS_given_dim(vtx_data* graph, int n, double* given_coords,
+ double* new_coords, double conj_tol)
+{
+ int iterations1,iterations2, mat_mult_count1;
+ int i,j;
+ DistType** Dij;
+ float* f_storage = NULL;
+ double* x = given_coords;
+ double uniLength;
+ double* orthog_aux = NULL;
+ double* y = new_coords;
+ float** lap = N_GNEW(n, float*);
+ float degree;
+ double pos_i;
+ double* balance = N_GNEW(n, double);
+ double b;
+ bool converged;
+
+ iterations1=mat_mult_count1=0; /* We don't compute the x-axis at all. */
+
+ Dij = compute_apsp(graph, n);
+
+ /* scaling up the distances to enable an 'sqrt' operation later
+ * (in case distances are integers)
+ */
+ for (i=0; i<n; i++)
+ for (j=0; j<n; j++)
+ Dij[i][j]*=SCALE_FACTOR;
+
+ if (x!=NULL) {
+ double sum1, sum2;
+ /* scale x (given axis) to minimize the stress */
+ orthog_aux = N_GNEW(n, double);
+ for (i=0; i<n; i++) {
+ orthog_aux[i]=x[i];
+ }
+ standardize(orthog_aux,n);
+
+ for (sum1=sum2=0,i=1; i<n; i++) {
+ for (j=0; j<i; j++) {
+ sum1+=1.0/(Dij[i][j])*fabs(x[i]-x[j]);
+ sum2+=1.0/(Dij[i][j]*Dij[i][j])*fabs(x[i]-x[j])*fabs(x[i]-x[j]);
+ }
+ }
+ uniLength=sum1/sum2;
+ for (i=0; i<n; i++)
+ x[i]*=uniLength;
+ }
+
+ /* smart ini: */
+ CMDS_orthog(graph, n, 1, &y, conj_tol, x, Dij);
+
+ /* Compute Laplacian: */
+ f_storage = N_GNEW(n*n, float);
+
+ for (i=0; i<n; i++) {
+ lap[i]=f_storage+i*n;
+ degree=0;
+ for (j=0; j<n; j++) {
+ if (j==i)
+ continue;
+ degree-=lap[i][j]=-1.0f/((float)Dij[i][j]*(float)Dij[i][j]); // w_{ij}
+
+ }
+ lap[i][i]=degree;
+ }
+
+
+ /* compute residual distances */
+ if (x!=NULL) {
+ double diff;
+ for (i=1; i<n; i++) {
+ pos_i=x[i];
+ for (j=0; j<i; j++) {
+ diff=(double)Dij[i][j]*(double)Dij[i][j]-(pos_i-x[j])*(pos_i-x[j]);
+ Dij[i][j]=Dij[j][i]=diff>0 ? (DistType)sqrt(diff) : 0;
+ }
+ }
+ }
+
+ /* Compute the balance vector: */
+ for (i=0; i<n; i++) {
+ pos_i=y[i];
+ balance[i]=0;
+ for (j=0; j<n; j++) {
+ if (j==i)
+ continue;
+ if (pos_i>=y[j]) {
+ balance[i]+=Dij[i][j]*(-lap[i][j]); // w_{ij}*delta_{ij}
+ }
+ else {
+ balance[i]-=Dij[i][j]*(-lap[i][j]); // w_{ij}*delta_{ij}
+ }
+ }
+ }
+
+ for (converged=false,iterations2=0; iterations2<200 && !converged; iterations2++) {
+ conjugate_gradient_f(lap, y, balance, n, conj_tol, n, true);
+ converged=true;
+ for (i=0; i<n; i++) {
+ pos_i=y[i];
+ b=0;
+ for (j=0; j<n; j++) {
+ if (j==i)
+ continue;
+ if (pos_i>=y[j]) {
+ b+=Dij[i][j]*(-lap[i][j]);
+
+ }
+ else {
+ b-=Dij[i][j]*(-lap[i][j]);
+
+ }
+ }
+ if ((b != balance[i]) && (fabs(1-b/balance[i])>1e-5)) {
+ converged=false;
+ balance[i]=b;
+ }
+ }
+ }
+
+ for (i=0; i<n; i++) {
+ x[i] /= uniLength;
+ y[i] /= uniLength;
+ }
+
+
+ free (Dij[0]); free (Dij);
+ free (lap[0]); free (lap);
+ free (orthog_aux); free (balance);
+}
+
+#endif /* DIGCOLA */
+
projects / graphviz / commitdiff