From: erg <devnull@localhost>
Date: Wed, 13 Jul 2005 20:35:21 +0000 (+0000)
Subject: Commit new files for hierarchical neato, plus changed makefiles.
X-Git-Tag: LAST_LIBGRAPH~32^2~7459
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=58c611ed54cec49829f4fda2551b85c6f9063a27;p=graphviz

Commit new files for hierarchical neato, plus changed makefiles.
Note that unless DIGCOLA is defined, the .c files below will be empty
compilation units.
---

diff --git a/lib/neatogen/Makefile.am b/lib/neatogen/Makefile.am
index be754761b..0353cf551 100644
--- a/lib/neatogen/Makefile.am
+++ b/lib/neatogen/Makefile.am
@@ -14,13 +14,15 @@ AM_CPPFLAGS = \
 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
diff --git a/lib/neatogen/Makefile.old b/lib/neatogen/Makefile.old
index 20c58e04e..e054437d0 100644
--- a/lib/neatogen/Makefile.old
+++ b/lib/neatogen/Makefile.old
@@ -18,7 +18,8 @@ HDRS = adjust.h edges.h geometry.h heap.h hedges.h info.h mem.h \
        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 \
diff --git a/lib/neatogen/compute_hierarchy.c b/lib/neatogen/compute_hierarchy.c
new file mode 100644
index 000000000..683b3b286
--- /dev/null
+++ b/lib/neatogen/compute_hierarchy.c
@@ -0,0 +1,135 @@
+/**********************************************************
+*      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 */
+
diff --git a/lib/neatogen/constrained_majorization.c b/lib/neatogen/constrained_majorization.c
new file mode 100644
index 000000000..9d948aaa8
--- /dev/null
+++ b/lib/neatogen/constrained_majorization.c
@@ -0,0 +1,481 @@
+/**********************************************************
+*      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 */
+
diff --git a/lib/neatogen/opt_arrangement.c b/lib/neatogen/opt_arrangement.c
new file mode 100644
index 000000000..79349c920
--- /dev/null
+++ b/lib/neatogen/opt_arrangement.c
@@ -0,0 +1,90 @@
+/**********************************************************
+*      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 */
+
diff --git a/lib/neatogen/quad_prog_solve.c b/lib/neatogen/quad_prog_solve.c
new file mode 100644
index 000000000..903bd03b9
--- /dev/null
+++ b/lib/neatogen/quad_prog_solve.c
@@ -0,0 +1,780 @@
+/**********************************************************
+*      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 */
+
diff --git a/lib/neatogen/quad_prog_solver.h b/lib/neatogen/quad_prog_solver.h
new file mode 100644
index 000000000..e0084b2fd
--- /dev/null
+++ b/lib/neatogen/quad_prog_solver.h
@@ -0,0 +1,58 @@
+/**********************************************************
+*      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
diff --git a/lib/neatogen/smart_ini_x.c b/lib/neatogen/smart_ini_x.c
new file mode 100644
index 000000000..5500aeda7
--- /dev/null
+++ b/lib/neatogen/smart_ini_x.c
@@ -0,0 +1,385 @@
+/**********************************************************
+*      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 */
+