From: erg <devnull@localhost>
Date: Fri, 28 Apr 2006 16:59:01 +0000 (+0000)
Subject: Add IPSEPCOLA additions from Tim Dwyer
X-Git-Tag: LAST_LIBGRAPH~32^2~6677
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=73c4c07aa0c16f2e8f267b7ae1b5c6b73181e609;p=graphviz

Add IPSEPCOLA additions from Tim Dwyer
---

diff --git a/lib/neatogen/Makefile.am b/lib/neatogen/Makefile.am
index a7b75f192..b71a97175 100644
--- a/lib/neatogen/Makefile.am
+++ b/lib/neatogen/Makefile.am
@@ -1,4 +1,3 @@
-# $Id$ $Revision$
 ## Process this file with automake to produce Makefile.in
 
 AM_CPPFLAGS = \
@@ -8,12 +7,14 @@ AM_CPPFLAGS = \
         -I$(top_srcdir)/lib/pack \
         -I$(top_srcdir)/lib/pathplan \
         -I$(top_srcdir)/lib/graph \
-        -I$(top_srcdir)/lib/cdt
+        -I$(top_srcdir)/lib/cdt \
+	-I$(top_srcdir)/lib/libvpsc
 
 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 quad_prog_solver.h digcola.h
+	matrix_ops.h pca.h stress.h quad_prog_solver.h digcola.h \
+	quad_prog_vpsc.h mosek_quad_solve.h
 noinst_LTLIBRARIES = libneatogen.la
 
 libneatogen_la_SOURCES = adjust.c circuit.c edges.c find_ints.c geometry.c \
@@ -22,6 +23,8 @@ libneatogen_la_SOURCES = adjust.c circuit.c edges.c find_ints.c geometry.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 quad_prog_solve.c \
 	smart_ini_x.c constrained_majorization.c opt_arrangement.c \
-	compute_hierarchy.c
+	compute_hierarchy.c \
+	quad_prog_vpsc.c constrained_majorization_ipsep.c \
+	mosek_quad_solve.c
 
 EXTRA_DIST = Makefile.old
diff --git a/lib/neatogen/constrained_majorization_ipsep.c b/lib/neatogen/constrained_majorization_ipsep.c
new file mode 100644
index 000000000..ad85c80db
--- /dev/null
+++ b/lib/neatogen/constrained_majorization_ipsep.c
@@ -0,0 +1,428 @@
+/**********************************************************
+ * Based on constrained_majorization.c
+ *
+ * Perform stress majorization subject
+ * to separation constraints, for background see the paper:
+ * "IPSep-CoLa: An Incremental Procedure for Separation Constraint Layout of Graphs"
+ * by Tim Dwyer, Yehuda Koren and Kim Marriott
+ *
+ * Available separation constraints so far are:
+ *  o Directed edge constraints
+ *  o Node non-overlap constraints
+ *  o Cluster containment constraints
+ *  o Cluster/node non-overlap constraints
+ *
+ * Tim Dwyer, 2006
+ **********************************************************/
+
+#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 <csolve_VPSC.h>
+#include "quad_prog_vpsc.h"
+#include "quad_prog_solver.h"
+#include "matrix_ops.h"
+
+#define localConstrMajorIterations 1000
+
+int 
+stress_majorization_cola(
+    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 model,          /* difference model */
+    int maxi,           /* max iterations */
+    ipsep_options* opt
+)
+{
+    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;
+	float * lap1 = NULL;
+	float * dist_accumulator = NULL;
+	float * tmp_coords = NULL;
+	float ** b = NULL;
+	double * degrees = NULL;
+	float * lap2=NULL;
+	int lap_length;
+	float * f_storage=NULL;
+	float ** coords=NULL;
+    int orig_n=n;
+
+	//double conj_tol=tolerance_cg;        /* tolerance of Conjugate Gradient */
+	CMajEnvVPSC *cMajEnvHor = NULL;
+	CMajEnvVPSC *cMajEnvVrt = NULL;
+	clock_t start_time;
+	double y_0;
+	int length;
+	DistType diameter;
+	float * Dij=NULL;
+	float constant_term;
+	int count;
+	double degree;
+	int step;
+	float val;
+	double old_stress, new_stress=0;
+	bool converged;
+	int len;
+    double nsizeScale=0;
+    float maxEdgeLen=0;
+
+    initLayout(graph, n, dim, d_coords);
+    if (n == 1) return 0;
+
+    for(i=0;i<n;i++) {
+        for(j=1;j<graph[i].nedges;j++) {
+            maxEdgeLen=MAX(graph[i].ewgts[j],maxEdgeLen);
+        }
+    }
+
+	/****************************************************
+	** 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\n ");
+        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];
+		}
+	}
+
+    /* 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;
+        }	
+    }
+
+	/**************************
+	** 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;
+	}
+
+	/**************************
+	** 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);
+    
+    if(opt->clusters->nclusters>0) {
+        int nn = n+opt->clusters->nclusters*2;
+        int clap_length = nn+nn*(nn-1)/2;
+        float *clap = N_GNEW(clap_length, float);
+        int c0,c1;
+        float v;
+        c0=c1=0;
+        for(i=0;i<nn;i++) {
+            for(j=0;j<nn-i;j++) {
+                if(i<n && j<n-i) {
+                    v=lap2[c0++];
+                } else {
+                    //v=j==1?i%2:0;
+                    if(j==1&&i%2==1) {
+                        v=maxEdgeLen;
+                        v*=v;
+                        if(v>0.01) {
+                            v=1.0/v;
+                        }
+                    }
+                    else v=0;
+                }
+                clap[c1++]=v;
+            }
+        }
+        free(lap2);
+        lap2=clap;
+        n=nn;
+        lap_length=clap_length;
+    }
+    /* 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];
+    }
+
+	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] = j<orig_n?(float)(d_coords[i][j]):0;
+		}
+	}
+
+	/* compute constant term in stress sum
+	 * which is \sum_{i<j} w_{ij}d_{ij}^2
+     */
+	constant_term=(float)(n*(n-1)/2);
+	
+	/*************************
+	** 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);
+	
+	old_stress=DBL_MAX; /* at least one iteration */
+
+	start_time = clock();
+
+	cMajEnvHor=initCMajVPSC(n,lap2,graph,opt,0);
+	cMajEnvVrt=initCMajVPSC(n,lap2,graph,opt,opt->diredges);
+
+	lap1 = N_GNEW(lap_length, float);
+
+	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);
+		sqrt_vecf(lap_length, lap2, lap1);
+		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		
+		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 */
+	    if (Verbose && (iterations % 1 == 0)) {
+            fprintf(stderr, "%.3f ", new_stress);
+            if (iterations % 10 == 0)
+                fprintf(stderr, "\n");
+	    }
+		converged = new_stress < old_stress && 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;
+
+        // in determining non-overlap constraints we gradually scale up the
+        // size of nodes to avoid local minima
+        if((iterations>=maxi-1||converged)&&opt->noverlap==1&&nsizeScale<0.999) {
+            nsizeScale+=0.1;
+            if(Verbose)
+                fprintf(stderr,"nsizescale=%f,iterations=%d\n",nsizeScale,iterations);
+            iterations=0;
+            converged = false;
+        }
+		
+
+        /* 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(opt->noverlap==1 && nsizeScale > 0.001) {
+            generateNonoverlapConstraints(cMajEnvHor,nsizeScale,coords,0,nsizeScale<0.5?false:true,opt);
+        }
+        if(cMajEnvHor->m > 0) {
+#ifdef MOSEK
+            if(opt->mosek) {
+                mosek_quad_solve_sep(cMajEnvHor->mosekEnv,n,b[0],coords[0]);
+            } else
+#endif // MOSEK
+            constrained_majorization_vpsc(cMajEnvHor,b[0],coords[0],localConstrMajorIterations);
+        } else {
+            // if there are no constraints then use conjugate gradient
+            // optimisation which should be considerably faster
+			conjugate_gradient_mkernel(lap2, coords[0], b[0], n, tolerance_cg, n);	
+        }
+        if(opt->noverlap==1 && nsizeScale > 0.001) {
+            generateNonoverlapConstraints(cMajEnvVrt,nsizeScale,coords,1,false,opt);
+        }
+        if(cMajEnvVrt->m > 0) {
+#ifdef MOSEK
+            if(opt->mosek) {
+                mosek_quad_solve_sep(cMajEnvVrt->mosekEnv,n,b[1],coords[1]);
+            } else
+#endif // MOSEK
+            constrained_majorization_vpsc(cMajEnvVrt,b[1],coords[1],localConstrMajorIterations);
+        } else {
+			conjugate_gradient_mkernel(lap2,coords[1],b[1],n,tolerance_cg,n);	
+        }
+	}
+    if (Verbose) {
+        fprintf(stderr, "\nfinal e = %f %d iterations %.2f sec\n",
+            new_stress, iterations, elapsed_sec());
+    }
+	deleteCMajEnvVPSC(cMajEnvHor);
+	deleteCMajEnvVPSC(cMajEnvVrt);
+
+    if (opt->noverlap==2) {
+        fprintf(stderr,"Removing overlaps as post-process...\n");
+        removeoverlaps(orig_n,coords,opt);
+    }
+	
+	if (coords!=NULL) {
+		for (i=0; i<dim; i++) {
+			for (j=0; j<orig_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);
+	free (lap1); 
+
+    return iterations;
+}
+#endif /* DIGCOLA */
+
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
diff --git a/lib/neatogen/defs.h b/lib/neatogen/defs.h
index 288f62224..d558452e3 100644
--- a/lib/neatogen/defs.h
+++ b/lib/neatogen/defs.h
@@ -33,6 +33,16 @@ extern "C" {
 	float *edists; /* directed dist reflecting the direction of the edge */
     };
 
+    struct cluster_data {
+	int nvars;         /* total count of vars in clusters */
+        int nclusters;     /* number of clusters */
+        int *clustersizes; /* number of vars in each cluster */
+        int **clusters;    /* list of var indices for constituents of each c */
+	int ntoplevel;     /* number of nodes not in any cluster */
+	int *toplevel;     /* array of nodes not in any cluster */
+	boxf *bb;	   /* bounding box of each cluster */
+    };
+
     typedef int DistType;	/* must be signed!! */
 
     inline double max(double x, double y) {
@@ -93,6 +103,17 @@ extern "C" {
 #endif
     } vtx_data;
 
+    typedef struct cluster_data {
+	int nvars;         /* total count of vars in clusters */
+        int nclusters;     /* number of clusters */
+        int *clustersizes; /* number of vars in each cluster */
+        int **clusters;    /* list of var indices for constituents of each c */
+	int ntoplevel;     /* number of nodes not in any cluster */
+	int *toplevel;     /* array of nodes not in any cluster */
+	boxf *bb;	   /* bounding box of each cluster */
+    } cluster_data;
+
+
     typedef int DistType;	/* must be signed!! */
 
 #ifdef UNUSED
diff --git a/lib/neatogen/digcola.h b/lib/neatogen/digcola.h
index a8e79bc46..db77d7923 100644
--- a/lib/neatogen/digcola.h
+++ b/lib/neatogen/digcola.h
@@ -26,6 +26,24 @@ extern double compute_hierarchy(vtx_data*, int, double, double,
 extern void IMDS_given_dim(vtx_data*, int, double*, double*, double);
 extern int stress_majorization_with_hierarchy(vtx_data*, int, int, double**, 
                                               int, int, int, int, double);
+typedef struct ipsep_options {
+    int diredges;       /* 1=generate directed edge constraints */
+                        /* 2=generate directed hierarchy level constraints (DiG-CoLa) */
+    double edge_gap;    /* amount to force vertical separation of */
+                        /* start/end nodes */
+    int noverlap;       /* 1=generate non-overlap constraints */
+                        /* 2=remove overlaps after layout */
+    pointf gap;         /* hor and vert gap to enforce when removing overlap*/
+    pointf* nsize;      /* node widths and heights */
+    cluster_data* clusters;
+                        /* list of node indices for each cluster */
+#ifdef MOSEK
+    int mosek;          /* use Mosek as constraint optimization engine */
+#endif //MOSEK
+} ipsep_options;
+
+// stress majorization, for Constraint Layout
+extern int stress_majorization_cola(vtx_data*, int, int, double**, int, int, int, ipsep_options*);
 #endif
 #endif
 
diff --git a/lib/neatogen/mosek_quad_solve.c b/lib/neatogen/mosek_quad_solve.c
new file mode 100644
index 000000000..06f196aff
--- /dev/null
+++ b/lib/neatogen/mosek_quad_solve.c
@@ -0,0 +1,485 @@
+/* 
+ * Interface to Mosek (www.mosek.com) quadratic programming solver for solving
+ * instances of the "Variable Placement with Separation Constraints" problem,
+ * and also DiG-CoLa style level constraints.
+ *
+ * Tim Dwyer, 2006
+ */
+#ifdef MOSEK
+#include <stdio.h>
+#include <assert.h>
+#include "defs.h"
+#include "mosek_quad_solve.h"
+#include "quad_prog_vpsc.h"
+
+//#define DUMP_CONSTRAINTS
+//#define EQUAL_WIDTH_LEVELS
+
+static FILE *logfile;
+static void MSKAPI printstr(void *handle, char str[])
+{
+	fprintf(logfile,"%s",str);
+} /* printstr */
+
+/**********************
+lap: the upper RHS of the symmetric graph laplacian matrix which will be transformed
+	to the hessian of the non-linear part of the optimisation function
+n: number of nodes (length of coords array)
+ordering: array containing sequences of nodes for each level,
+	ie, ordering[levels[i]] is first node of (i+1)th level
+level_indexes: array of starting node for each level in ordering
+	ie, levels[i] is index to first node of (i+1)th level
+	also, levels[0] is number of nodes in first level
+	and, levels[i]-levels[i-1] is number of nodes in ith level
+	and, n - levels[num_divisions-1] is number of nodes in last level
+num_divisions: number of divisions between levels, ie number of levels - 1
+separation: the minimum separation between nodes on different levels
+***********************/
+MosekEnv* mosek_init_hier(
+	float* lap, int n,int *ordering,int *level_indexes,int num_divisions, float separation)
+{
+	int i,j,num_levels=num_divisions+1;
+	MosekEnv *mskEnv = GNEW(MosekEnv);
+	// vars for nodes (except x0) + dummy nodes between levels
+	// x0 is fixed at 0, and therefore is not included in the opt problem
+	// add 2 more vars for top and bottom constraints
+	mskEnv->num_variables=n+num_divisions+1;
+
+	logfile = fopen("quad_solve_log","w");
+	DigColaLevel *levels=assign_digcola_levels(ordering,n,level_indexes,num_divisions);
+#ifdef DUMP_CONSTRAINTS
+	print_digcola_levels(logfile,levels,num_levels);
+#endif
+
+	/* nonlinear coefficients matrix of objective function */
+	//int lapsize=mskEnv->num_variables+(mskEnv->num_variables*(mskEnv->num_variables-1))/2;
+	int nonzero_lapsize=(n*(n-1))/2;
+	mskEnv->qval = N_GNEW(nonzero_lapsize,double);
+	mskEnv->qsubi = N_GNEW(nonzero_lapsize,int);
+	mskEnv->qsubj = N_GNEW(nonzero_lapsize,int);
+
+	/* solution vector */
+	mskEnv->xx = N_GNEW(mskEnv->num_variables,double);
+
+	/* constraint matrix */
+	separation/=2.0; // separation between each node and it's adjacent constraint
+	int num_constraints = get_num_digcola_constraints(levels, num_levels)+num_divisions+1;
+	// constraints of the form x_i - x_j >= sep so 2 non-zero entries per constraint in LHS matrix
+	// except x_0 (fixed at 0) constraints which have 1 nz val each.
+#ifdef EQUAL_WIDTH_LEVELS
+	num_constraints+=num_divisions;
+#endif
+	// pointer to beginning of nonzero sequence in a column
+
+	int count=0;
+	for(i=0;i<n-1;i++) {
+		for(j=i;j<n-1;j++) {
+			mskEnv->qval[count]=-2*lap[count+n];
+			assert(mskEnv->qval[count]!=0);
+			mskEnv->qsubi[count]=j;
+			mskEnv->qsubj[count]=i;
+			count++;
+		}
+	}
+#ifdef DUMP_CONSTRAINTS
+	fprintf(logfile,"Q=[");
+	int lapcntr=n;
+	for(i=0;i<mskEnv->num_variables;i++) {
+		if(i!=0) fprintf(logfile,";");
+		for(j=0;j<mskEnv->num_variables;j++) {
+			if(j<i||i>=n-1||j>=n-1) {
+				fprintf(logfile,"0 ");
+			} else {
+				fprintf(logfile,"%f ",-2*lap[lapcntr++]);
+			}
+		}
+	}
+	fprintf(logfile,"]\nQ=Q-diag(diag(Q))+Q'\n");
+#endif
+	fprintf(logfile,"\n");
+	/* Make the mosek environment. */
+	mskEnv->r = MSK_makeenv(&mskEnv->env,NULL,NULL,NULL,NULL);
+
+	/* Check whether the return code is ok. */
+	if ( mskEnv->r==MSK_RES_OK )
+	{
+		/* Directs the log stream to the user
+		specified procedure 'printstr'. */
+		MSK_linkfunctoenvstream(mskEnv->env,MSK_STREAM_LOG,NULL,printstr);
+	}
+
+	/* Initialize the environment. */   
+	mskEnv->r = MSK_initenv(mskEnv->env);
+	if (mskEnv->r==MSK_RES_OK)
+	{ 
+		/* Make the optimization task. */
+		mskEnv->r = MSK_maketask(mskEnv->env,num_constraints,mskEnv->num_variables,&mskEnv->task);
+
+		if ( mskEnv->r==MSK_RES_OK )
+		{
+			mskEnv->r = MSK_linkfunctotaskstream(mskEnv->task,MSK_STREAM_LOG,NULL,printstr);
+			/* Resize the task. */
+			if ( mskEnv->r==MSK_RES_OK )
+				mskEnv->r = MSK_resizetask(mskEnv->task,num_constraints,mskEnv->num_variables,
+				0, // no cones!!
+				// each constraint applies to 2 vars
+				2*num_constraints+num_divisions, 
+				nonzero_lapsize);
+
+			/* Append the constraints. */
+			if ( mskEnv->r==MSK_RES_OK )
+				mskEnv->r = MSK_append(mskEnv->task,1,num_constraints);
+
+			/* Append the variables. */
+			if ( mskEnv->r==MSK_RES_OK )
+				mskEnv->r = MSK_append(mskEnv->task,0,mskEnv->num_variables);
+			/* Put variable bounds. */
+			for(j=0; j<mskEnv->num_variables && mskEnv->r==MSK_RES_OK; ++j)
+				mskEnv->r = MSK_putbound(mskEnv->task,0,j,MSK_BK_RA,-MSK_INFINITY,MSK_INFINITY);
+			int c_ind=0;
+			int c_var=n-1;
+			for(j=0;j<levels[0].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
+				int node = levels[0].nodes[j]-1;
+				if(node>=0){
+					MSKidxt subi[]={c_var,node};
+					double vali[]={1.0,-1.0};
+					mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
+				} else { 
+					// constraint for y0 (fixed at 0)
+					mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,1.0);
+				}
+				mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
+				c_ind++;
+			}	
+			for(i=0;i<num_divisions && mskEnv->r==MSK_RES_OK;i++) {
+				c_var=n+i;
+				for(j=0;j<levels[i].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
+					// create separation constraint a>=b+separation
+					int node = levels[i].nodes[j]-1;
+					if(node>=0){ // no constraint for fixed node
+						MSKidxt subi[]={node,c_var};
+						double vali[]={1.0,-1.0};
+						mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
+					} else { 
+						// constraint for y0 (fixed at 0)
+						mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,-1.0);
+					}
+					mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
+					c_ind++;
+				}
+				for(j=0;j<levels[i+1].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
+					int node = levels[i+1].nodes[j]-1;
+					if(node>=0){
+						MSKidxt subi[]={c_var,node};
+						double vali[]={1.0,-1.0};
+						mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
+					} else { 
+						// constraint for y0 (fixed at 0)
+						mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,1.0);
+					}
+					mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
+					c_ind++;
+				}		
+			}
+			c_var=n+i;
+			for(j=0;j<levels[i].num_nodes && mskEnv->r==MSK_RES_OK;j++) {
+				// create separation constraint a>=b+separation
+				int node = levels[i].nodes[j]-1;
+				if(node>=0){ // no constraint for fixed node
+					MSKidxt subi[]={node,c_var};
+					double vali[]={1.0,-1.0};
+					mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
+				} else { 
+					// constraint for y0 (fixed at 0)
+					mskEnv->r=MSK_putaij(mskEnv->task,c_ind,c_var,-1.0);
+				}
+				mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,separation,MSK_INFINITY);
+				c_ind++;
+			}
+			// create constraints preserving the order of dummy vars
+			for(i=0;i<num_divisions+1 && mskEnv->r==MSK_RES_OK;i++) {
+				int c_var=n-1+i, c_var2=c_var+1;
+				MSKidxt subi[]={c_var,c_var2};
+				double vali[]={1.0,-1.0};
+				mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,2,subi,vali);
+				mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_LO,0,MSK_INFINITY);
+				c_ind++;
+			}
+#ifdef EQUAL_WIDTH_LEVELS
+			for(i=1;i<num_divisions+1 && mskEnv->r==MSK_RES_OK;i++) {
+				int c_var=n-1+i, c_var_lo=c_var-1, c_var_hi=c_var+1;
+				MSKidxt subi[]={c_var_lo,c_var,c_var_hi};
+				double vali[]={1.0,-2.0,1.0};
+				mskEnv->r=MSK_putavec(mskEnv->task,1,c_ind,3,subi,vali);
+				mskEnv->r=MSK_putbound(mskEnv->task,1,c_ind,MSK_BK_FX,0,0);
+				c_ind++;
+			}
+#endif
+			assert(c_ind==num_constraints);
+#ifdef DUMP_CONSTRAINTS
+			fprintf(logfile,"A=[");
+			for(i=0;i<num_constraints;i++) {
+				if(i!=0) fprintf(logfile,";");
+				for(j=0;j<mskEnv->num_variables;j++) {
+					double aij;
+					MSK_getaij(mskEnv->task,i,j,&aij);
+					fprintf(logfile,"%f ",aij);
+				}
+			}
+			fprintf(logfile,"]\n");
+			fprintf(logfile,"b=[");
+			for(i=0;i<num_constraints;i++) {
+				fprintf(logfile,"%f ",separation);
+			}
+			fprintf(logfile,"]\n");
+#endif
+			if ( mskEnv->r==MSK_RES_OK )
+			{
+				/*
+				* The lower triangular part of the Q
+				* matrix in the objective is specified.
+				*/
+				mskEnv->r = MSK_putqobj(mskEnv->task,nonzero_lapsize,mskEnv->qsubi,mskEnv->qsubj,mskEnv->qval);
+			}
+		}
+	}
+	delete_digcola_levels(levels,num_levels);
+	return mskEnv;
+}
+/*
+b: coefficients of linear part of optimisation function
+n: number of nodes
+coords: optimal y* vector, coord[i] is coordinate of node[i]
+hierarchy_boundaries: y coord of boundaries between levels 
+	(ie, solution values for the dummy variables used in constraints)
+*/
+void mosek_quad_solve_hier(MosekEnv *mskEnv, float *b,int n,float* coords, float *hierarchy_boundaries) {
+	int i,j;
+	for(i=1;i<n && mskEnv->r==MSK_RES_OK;i++) {
+		mskEnv->r = MSK_putcj(mskEnv->task,i-1,-2*b[i]);
+	}
+#ifdef DUMP_CONSTRAINTS
+	fprintf(logfile,"x0=[");
+	for(j=0;j<mskEnv->num_variables;j++) {
+		fprintf(logfile,"%f ",j<n?b[j]:0);
+	}
+	fprintf(logfile,"]\n");
+	fprintf(logfile,"f=[");
+	double *c=N_GNEW(mskEnv->num_variables,double);
+	MSK_getc(mskEnv->task,c);
+	for(j=0;j<mskEnv->num_variables;j++) {
+		fprintf(logfile,"%f ",c[j]);
+	}
+	free(c);
+	fprintf(logfile,"]\n");
+#endif
+	if ( mskEnv->r==MSK_RES_OK )
+		mskEnv->r = MSK_optimize(mskEnv->task);
+
+	if ( mskEnv->r==MSK_RES_OK )
+	{
+		MSK_getsolutionslice(mskEnv->task,
+			MSK_SOL_ITR,
+			MSK_SOL_ITEM_XX,
+			0,
+			mskEnv->num_variables,
+			mskEnv->xx);
+
+#ifdef DUMP_CONSTRAINTS
+		fprintf(logfile,"Primal solution\n");
+#endif
+		coords[0]=0;
+		for(j=0; j<mskEnv->num_variables; ++j) {
+#ifdef DUMP_CONSTRAINTS
+			fprintf(logfile,"x[%d]: %.2f\n",j,mskEnv->xx[j]);
+#endif
+			if(j<n-1) {
+				coords[j+1]=-mskEnv->xx[j];
+			} else if(j>=n&&j<mskEnv->num_variables-1) {
+				hierarchy_boundaries[j-n]=-mskEnv->xx[j];
+			}
+		}
+	}
+	fprintf(logfile,"Return code: %d\n",mskEnv->r);
+}
+/**********************
+lap: the upper RHS of the symmetric graph laplacian matrix which will be transformed
+	to the hessian of the non-linear part of the optimisation function
+	has dimensions num_variables, dummy vars do not have entries in lap
+cs: array of pointers to separation constraints
+***********************/
+MosekEnv* mosek_init_sep(
+	float* lap, int num_variables, int num_dummy_vars, Constraint **cs, int num_constraints)
+{
+	int i,j;
+	MosekEnv *mskEnv = GNEW(MosekEnv);
+    // fix var 0
+	mskEnv->num_variables=num_variables+num_dummy_vars-1;
+
+	fprintf(stderr,"MOSEK!\n");
+	logfile = fopen("quad_solve_log","w");
+
+	/* nonlinear coefficients matrix of objective function */
+	int nonzero_lapsize=num_variables*(num_variables-1)/2;
+	mskEnv->qval = N_GNEW(nonzero_lapsize,double);
+	mskEnv->qsubi = N_GNEW(nonzero_lapsize,int);
+	mskEnv->qsubj = N_GNEW(nonzero_lapsize,int);
+
+	/* solution vector */
+	mskEnv->xx = N_GNEW(mskEnv->num_variables,double);
+
+	// pointer to beginning of nonzero sequence in a column
+
+	int count=0;
+	for(i=0;i<num_variables-1;i++) {
+		for(j=i;j<num_variables-1;j++) {
+			mskEnv->qval[count]=-2*lap[count+num_variables];
+			//assert(mskEnv->qval[count]!=0);
+			mskEnv->qsubi[count]=j;
+			mskEnv->qsubj[count]=i;
+			count++;
+		}
+	}
+#ifdef DUMP_CONSTRAINTS
+	fprintf(logfile,"Q=[");
+	count=0;
+	for(i=0;i<num_variables-1;i++) {
+		if(i!=0) fprintf(logfile,";");
+		for(j=0;j<num_variables-1;j++) {
+			if(j<i) {
+				fprintf(logfile,"0 ");
+			} else {
+				fprintf(logfile,"%f ",-2*lap[num_variables+count++]);
+			}
+		}
+	}
+	fprintf(logfile,"]\nQ=Q-diag(diag(Q))+Q'\n");
+#endif
+	/* Make the mosek environment. */
+	mskEnv->r = MSK_makeenv(&mskEnv->env,NULL,NULL,NULL,NULL);
+
+	/* Check whether the return code is ok. */
+	if ( mskEnv->r==MSK_RES_OK )
+	{
+		/* Directs the log stream to the user
+		specified procedure 'printstr'. */
+		MSK_linkfunctoenvstream(mskEnv->env,MSK_STREAM_LOG,NULL,printstr);
+	}
+
+	/* Initialize the environment. */   
+	mskEnv->r = MSK_initenv(mskEnv->env);
+	if (mskEnv->r==MSK_RES_OK)
+	{ 
+		/* Make the optimization task. */
+		mskEnv->r = MSK_maketask(mskEnv->env,num_constraints,mskEnv->num_variables,&mskEnv->task);
+
+		if ( mskEnv->r==MSK_RES_OK )
+		{
+			mskEnv->r = MSK_linkfunctotaskstream(mskEnv->task,MSK_STREAM_LOG,NULL,printstr);
+			/* Resize the task. */
+			if ( mskEnv->r==MSK_RES_OK )
+				mskEnv->r = MSK_resizetask(mskEnv->task,num_constraints,mskEnv->num_variables,
+				0, // no cones!!
+				// number of non-zero constraint matrix entries:
+                //   each constraint applies to 2 vars
+				2*num_constraints, 
+				nonzero_lapsize);
+
+			/* Append the constraints. */
+			if ( mskEnv->r==MSK_RES_OK )
+				mskEnv->r = MSK_append(mskEnv->task,1,num_constraints);
+
+			/* Append the variables. */
+			if ( mskEnv->r==MSK_RES_OK )
+				mskEnv->r = MSK_append(mskEnv->task,0,mskEnv->num_variables);
+			/* Put variable bounds. */
+			for(j=0; j<mskEnv->num_variables && mskEnv->r==MSK_RES_OK; j++)
+				mskEnv->r = MSK_putbound(mskEnv->task,0,j,MSK_BK_RA,-MSK_INFINITY,MSK_INFINITY);
+            for(i=0;i<num_constraints;i++) {
+                int u = getLeftVarID(cs[i])-1;
+                int v = getRightVarID(cs[i])-1;
+                double separation = getSeparation(cs[i]);
+                if(u<0) {
+				    mskEnv->r = MSK_putbound(mskEnv->task,0,v,MSK_BK_RA,-MSK_INFINITY,-separation);
+			        assert ( mskEnv->r==MSK_RES_OK );
+                } else if(v<0) {
+				    mskEnv->r = MSK_putbound(mskEnv->task,0,u,MSK_BK_RA,separation,MSK_INFINITY);
+			        assert ( mskEnv->r==MSK_RES_OK );
+                } else {
+                    //fprintf(stderr,"u=%d,v=%d,sep=%f\n",u,v,separation);
+                    MSKidxt subi[]={u,v};
+                    double vali[]={1.0,-1.0};
+                    mskEnv->r=MSK_putavec(mskEnv->task,1,i,2,subi,vali);
+			        assert ( mskEnv->r==MSK_RES_OK );
+                    mskEnv->r=MSK_putbound(mskEnv->task,1,i,MSK_BK_LO,separation,MSK_INFINITY);
+			        assert ( mskEnv->r==MSK_RES_OK );
+                }
+            }
+			if ( mskEnv->r==MSK_RES_OK )
+			{
+				/*
+				* The lower triangular part of the Q
+				* matrix in the objective is specified.
+				*/
+				mskEnv->r = MSK_putqobj(mskEnv->task,nonzero_lapsize,mskEnv->qsubi,mskEnv->qsubj,mskEnv->qval);
+                assert ( mskEnv->r==MSK_RES_OK );
+			}
+		}
+	}
+	return mskEnv;
+}
+/*
+n: size of b and coords, may be smaller than mskEnv->num_variables if we
+have dummy vars
+b: coefficients of linear part of optimisation function
+coords: optimal y* vector, coord[i] is coordinate of node[i]
+*/
+void mosek_quad_solve_sep(MosekEnv *mskEnv,int n, float *b,float* coords) {
+	int i,j;
+	assert(n<=mskEnv->num_variables+1);
+	for(i=0;i<n-1 && mskEnv->r==MSK_RES_OK;i++) {
+		mskEnv->r = MSK_putcj(mskEnv->task,i,-2*b[i+1]);
+	}
+	if ( mskEnv->r==MSK_RES_OK )
+		mskEnv->r = MSK_optimize(mskEnv->task);
+
+	if ( mskEnv->r==MSK_RES_OK )
+	{
+		MSK_getsolutionslice(mskEnv->task,
+			MSK_SOL_ITR,
+			MSK_SOL_ITEM_XX,
+			0,
+			mskEnv->num_variables,
+			mskEnv->xx);
+
+#ifdef DUMP_CONSTRAINTS
+		fprintf(logfile,"Primal solution\n");
+#endif
+		coords[0]=0;
+		for(j=1; j<=n; j++) {
+#ifdef DUMP_CONSTRAINTS
+			fprintf(logfile,"x[%d]: %.2f\n",j,mskEnv->xx[j-1]);
+#endif
+			coords[j]=-mskEnv->xx[j-1];
+		}
+	}
+	fprintf(logfile,"Return code: %d\n",mskEnv->r);
+}
+/*
+please call to clean up
+*/
+void mosek_delete(MosekEnv *mskEnv) {
+	MSK_deletetask(&mskEnv->task);
+	MSK_deleteenv(&mskEnv->env);
+
+	if(logfile) {
+		fclose(logfile);
+		logfile=NULL;
+	}
+	free(mskEnv->qval);
+	free(mskEnv->qsubi);
+	free(mskEnv->qsubj);
+	free(mskEnv->xx);
+	free(mskEnv);
+}
+#endif /* MOSEK */
diff --git a/lib/neatogen/mosek_quad_solve.h b/lib/neatogen/mosek_quad_solve.h
new file mode 100644
index 000000000..571cda247
--- /dev/null
+++ b/lib/neatogen/mosek_quad_solve.h
@@ -0,0 +1,27 @@
+#ifdef MOSEK
+#ifndef _QSOLVE_H_
+#define _QSOLVE_H_
+
+#include <mosek.h> /* Include the MOSEK definition file. */
+#include "types.h"
+#include <csolve_VPSC.h>
+
+typedef struct {
+	int r;
+	MSKenv_t  env;
+	MSKtask_t task;
+	double *qval;
+	int *qsubi,*qsubj;
+	double *xx;
+	int num_variables;
+} MosekEnv;
+
+MosekEnv* mosek_init_hier(float* lap, int n,int *ordering,int *level_indexes,int num_divisions, float separation);
+void mosek_quad_solve_hier(MosekEnv*,float *b,int n,float* coords, float *hierarchy_boundaries);
+MosekEnv* mosek_init_sep(float* lap, int nv, int ndv, Constraint** cs, int m);
+void mosek_quad_solve_sep(MosekEnv*,int n,float *b,float* coords);
+void mosek_delete(MosekEnv*);
+
+
+#endif // _QSOLVE_H_
+#endif // MOSEK
diff --git a/lib/neatogen/neato.h b/lib/neatogen/neato.h
index 7671d95cb..fcd279678 100644
--- a/lib/neatogen/neato.h
+++ b/lib/neatogen/neato.h
@@ -29,6 +29,7 @@
 #define MODE_KK          0
 #define MODE_MAJOR       1
 #define MODE_HIER        2
+#define MODE_IPSEP       3
 
 #define INIT_SELF        0
 #define INIT_REGULAR     1
diff --git a/lib/neatogen/neatoinit.c b/lib/neatogen/neatoinit.c
index 065f4eb79..e46cb8c2e 100644
--- a/lib/neatogen/neatoinit.c
+++ b/lib/neatogen/neatoinit.c
@@ -225,6 +225,78 @@ static void set_elabel(edge_t * e, textlabel_t * l, char *name)
     }
 }
 
+static cluster_data* cluster_map(graph_t *mastergraph, graph_t *g) {
+    /* search meta-graph to find clusters */
+    graph_t *mg, *subg;
+    node_t *mm, *mn;
+    edge_t *me;
+    // array of arrays of node indices in each cluster
+    int **cs,*cn;
+    int i,j,nclusters=0;
+    bool assigned[agnnodes(g)];
+    for(i=0;i<agnnodes(g);i++) {
+        assigned[i]=false;
+    }
+    cluster_data *cdata = GNEW(cluster_data);
+    cdata->ntoplevel=agnnodes(g);
+    mm = mastergraph->meta_node;
+    mg = mm->graph;
+    for (me = agfstout(mg, mm); me; me = agnxtout(mg, me)) {
+        mn = me->head;
+        subg = agusergraph(mn);
+        if (!strncmp(subg->name, "cluster", 7)) {
+            nclusters++;
+        }
+    }
+    cdata->nvars=0;
+    cdata->nclusters=nclusters;
+    cs=cdata->clusters=N_GNEW(nclusters,int*);
+    cn=cdata->clustersizes=N_GNEW(nclusters,int);
+    //fprintf(stderr,"search %d clusters...\n",nclusters);
+    for (me = agfstout(mg, mm); me; me = agnxtout(mg, me)) {
+        mn = me->head;
+        subg = agusergraph(mn);
+        /* clusters are processed by separate calls to ordered_edges */
+        if (!strncmp(subg->name, "cluster", 7)) {
+            *cn=agnnodes(subg);
+            cdata->nvars+=*cn;
+            int *c=*cs++=N_GNEW(*cn++,int);
+            node_t *n;
+            //fprintf(stderr,"Cluster with %d nodes...\n",agnnodes(subg));
+            for (n = agfstnode(subg); n; n = agnxtnode(subg, n)) {
+                node_t *gn;
+                int ind=0;
+                for (gn = agfstnode(g); gn; gn = agnxtnode(g, gn)) {
+                    if(gn->id==n->id) break;
+                    ind++;
+                }
+                //fprintf(stderr,"  node=%s, id=%d, ind=%d\n",n->name,n->id,ind);
+                *c++=ind;
+                assigned[ind]=true;
+                cdata->ntoplevel--;
+            }
+        }
+    }
+    cdata->bb=N_GNEW(cdata->nclusters,boxf);
+    cdata->toplevel=N_GNEW(cdata->ntoplevel,int);
+    for(i=j=0;i<agnnodes(g);i++) {
+        if(!assigned[i]) {
+            cdata->toplevel[j++]=i;
+        }
+    }
+    assert(cdata->ntoplevel==agnnodes(g)-cdata->nvars);
+    return cdata;
+}
+static void freeClusterData(cluster_data *c) {
+    if(c->nclusters>0) {
+        free(c->clusters[0]);
+        free(c->clusters);
+        free(c->clustersizes);
+        free(c->toplevel);
+        free(c->bb);
+    }
+    free(c);
+}
 /* user_spline:
  * Attempt to use already existing pos info for spline
  * Return 1 if successful, 0 otherwise.
@@ -638,6 +710,8 @@ static int neatoMode(graph_t * g)
 #ifdef DIGCOLA
 	else if (streq(str, "hier"))
 	    mode = MODE_HIER;
+        else if (streq(str, "ipsep"))
+            mode = MODE_IPSEP;
 #endif
 	else
 	    agerr(AGWARN,
@@ -670,7 +744,7 @@ static int checkEdge(PointMap * pm, edge_t * ep, int idx)
  * dfs for breaking cycles in vtxdata
  */
 static void
-dfsCycle (vtx_data* graph, int i)
+dfsCycle (vtx_data* graph, int i,int mode)
 {
     node_t *np, *hp;
     int j, e, f;
@@ -683,12 +757,14 @@ dfsCycle (vtx_data* graph, int i)
 	j = graph[i].edges[e];
 	hp = graph[j].np;
 	if (ND_onstack(hp)) {  /* back edge: reverse it */
-	    graph[i].edists[e] = 1.0;
-	    for (f = 1; (f < graph[j].nedges) &&(graph[j].edges[f] != i); f++) ;
-	    assert (f < graph[j].nedges);
-	    graph[j].edists[f] = -1.0;
+            // if mode is IPSEP make it an in-edge 
+            // at both ends, so that an edge constraint won't be generated!
+            graph[i].edists[e] = mode==MODE_IPSEP?-1.0:1.0;
+            for (f = 1; (f < graph[j].nedges) &&(graph[j].edges[f] != i); f++) ;
+            assert (f < graph[j].nedges);
+            graph[j].edists[f] = -1.0;
         }
-	else if (ND_mark(hp) == FALSE) dfsCycle(graph, j);
+	else if (ND_mark(hp) == FALSE) dfsCycle(graph, j, mode);
 
     }
     ND_onstack(np) = FALSE;
@@ -698,7 +774,7 @@ dfsCycle (vtx_data* graph, int i)
  * Do a dfs of the vtx_data, looking for cycles, reversing edges.
  */
 static void
-acyclic (vtx_data* graph, int nv)
+acyclic (vtx_data* graph, int nv, int mode)
 {
     int i;
     node_t* np;
@@ -710,7 +786,7 @@ acyclic (vtx_data* graph, int nv)
     }
     for (i = 0; i < nv; i++) {
 	if (ND_mark(graph[i].np)) continue;
-	dfsCycle (graph, i);	
+	dfsCycle (graph, i, mode);	
     }
 
 }
@@ -760,7 +836,7 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int mode, int m
 	haveLen = (agindex(g->root->proto->e, "len") >= 0);
 	haveWt = (E_weight != 0);
     }
-    if (mode == MODE_HIER)
+    if (mode == MODE_HIER || mode == MODE_IPSEP)
 	haveDir = TRUE;
     else
 	haveDir = FALSE;
@@ -825,9 +901,14 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int mode, int m
 		else if (haveDir)
 		    *ewgts++ = 1.0;
 #ifdef DIGCOLA
-		if (haveDir) {
-			*edists++ = (np == ep->head ? 1.0 : -1.0);
-		}
+                if (haveDir) {
+                    char *s=agget(ep,"dir");
+                    if(s&&!strncmp(s,"none",4)) {
+                        *edists++ = 0;
+                    } else {
+                        *edists++ = (np == ep->head ? 1.0 : -1.0);
+                    }
+                }
 #endif
 		i_nedges++;
 	    }
@@ -843,7 +924,7 @@ static vtx_data *makeGraphData(graph_t * g, int nv, int *nedges, int mode, int m
 #ifdef DIGCOLA
     if (haveDir) {
     /* Make graph acyclic */
-	acyclic (graph, nv);
+	acyclic (graph, nv, mode);
     }
 #endif
 
@@ -1010,13 +1091,15 @@ void dumpData(graph_t * g, vtx_data * gp, int nv, int ne)
  *  weight?
  */
 static void
-majorization(graph_t * g, int nv, int mode, int model, int dim, int steps)
+majorization(graph_t *mg, graph_t * g, int nv, int mode, int model, int dim, int steps)
 {
     double **coords;
     int ne;
     int i;
     node_t *v;
     vtx_data *gp;
+    cluster_data *cs;
+
     int init;
 
     init = checkStart(g, nv, (mode == MODE_HIER ? INIT_SELF : INIT_RANDOM));
@@ -1032,15 +1115,84 @@ majorization(graph_t * g, int nv, int mode, int model, int dim, int steps)
 	start_timer();
     }
     gp = makeGraphData(g, nv, &ne, mode, model);
+
+    cs=(cluster_data*)cluster_map(mg,g);
+
     if (Verbose) {
 	fprintf(stderr, "%d nodes %.2f sec\n", nv, elapsed_sec());
     }
 
 #ifdef DIGCOLA
-    if (mode == MODE_HIER) {
-	double lgap = late_double(g, agfindattr(g, "levelsgap"), 0.0, -MAXDOUBLE);
-	stress_majorization_with_hierarchy(gp, nv, ne, coords, Ndim,
-				   (init == INIT_SELF), model, MaxIter, lgap);
+    if (mode == MODE_HIER || mode == MODE_IPSEP) {
+        double lgap = late_double(g, agfindattr(g, "levelsgap"), 0.0, -MAXDOUBLE);
+        if (mode == MODE_HIER) {        
+            stress_majorization_with_hierarchy(gp, nv, ne, coords, Ndim,
+                       (init == INIT_SELF), model, MaxIter, lgap);
+        } else {
+            char* str;
+            ipsep_options opt;
+            pointf nsize[nv];
+            opt.edge_gap=lgap;
+#ifdef MOSEK
+            opt.mosek=0;
+#endif // MOSEK
+            opt.noverlap=opt.diredges=0;
+            opt.gap.x=opt.gap.y=0;
+            opt.nsize=nsize;
+            opt.clusters=cs;
+            str = agget(g, "diredgeconstraints");
+            if(str && !strncmp(str,"true",4)) {
+                opt.diredges = 1;
+                if(Verbose)
+                    fprintf(stderr,"Generating Edge Constraints...\n");
+            } else if(str && !strncmp(str,"hier",4)) {
+                opt.diredges = 2;
+                if(Verbose)
+                    fprintf(stderr,"Generating DiG-CoLa Edge Constraints...\n");
+            }
+            str = agget(g, "overlapconstraints");
+            if(str && !strncmp(str,"true",4)) {
+                opt.noverlap = 1;
+                if(Verbose)
+                    fprintf(stderr,"Generating Non-overlap Constraints...\n");
+            } else if(str && !strncmp(str,"post",4)) {
+                opt.noverlap = 2;
+                if(Verbose)
+                    fprintf(stderr,"Removing overlaps as postprocess...\n");
+            }  
+#ifdef MOSEK
+            str = agget(g, "mosek");
+            if(str && !strncmp(str,"true",4)) {
+                opt.mosek = 1;
+                if(Verbose)
+                    fprintf(stderr,"Using Mosek for constraint optimization...\n");
+            }
+#endif // MOSEK
+            if ((str = agget(g, "sep"))) {
+                char* c;
+                if ((c=strchr(str,','))) {
+                    i=c-str;
+                    if(Verbose)
+                        fprintf(stderr,"gap=%s,%d\n",str,i);
+                    char s[strlen(str)];
+                    strncpy(s,str,i);
+                    s[i]=0;
+                    opt.gap.x=atof(s);
+                    strcpy(s,c+1);
+                    opt.gap.y=atof(s);
+                    if(Verbose)
+                        fprintf(stderr,"gap=%f,%f\n",opt.gap.x,opt.gap.y);
+                } else {
+                    opt.gap.x =opt.gap.y = atof(str);
+                }
+            }
+            for (i=0, v = agfstnode(g); v; v = agnxtnode(g, v),i++) {
+                nsize[i].x=ND_width(v);
+                nsize[i].y=ND_height(v);
+            }
+
+            stress_majorization_cola(gp, nv, ne, coords, Ndim, model, MaxIter, &opt);
+        }
     }
     else
 #endif
@@ -1056,6 +1208,7 @@ majorization(graph_t * g, int nv, int mode, int model, int dim, int steps)
 	}
     }
     freeGraphData(gp);
+    freeClusterData(cs);
     free(coords[0]);
     free(coords);
 }
@@ -1112,7 +1265,7 @@ static void kkNeato(Agraph_t * g, int nG, int model)
 /* neatoLayout:
  * Use stress optimization to layout a single component
  */
-void neatoLayout(Agraph_t * g, int layoutMode, int layoutModel)
+void neatoLayout(Agraph_t * mg, Agraph_t * g, int layoutMode, int layoutModel)
 {
     int nG;
     char *str;
@@ -1128,7 +1281,7 @@ void neatoLayout(Agraph_t * g, int layoutMode, int layoutModel)
     if (!nG)
 	return;
     if (layoutMode)
-	majorization(g, nG, layoutMode, layoutModel, Ndim, MaxIter);
+	majorization(mg, g, nG, layoutMode, layoutModel, Ndim, MaxIter);
     else
 	kkNeato(g, nG, layoutModel);
 }
@@ -1200,7 +1353,7 @@ void neato_layout(Agraph_t * g)
 	    for (i = 0; i < n_cc; i++) {
 		gc = cc[i];
 		nodeInduce(gc);
-		neatoLayout(gc, layoutMode, model);
+		neatoLayout(g, gc, layoutMode, model);
 		adjustNodes(gc);
 	    }
 	    if (n_cc > 1) {
@@ -1228,8 +1381,23 @@ void neato_layout(Agraph_t * g)
 		free_scan_graph(gc);
 		agdelete(g, gc);
 	    }
+            {
+                graph_t *mg, *subg;
+                node_t *mm, *mn;
+                edge_t *me;
+                mm = g->meta_node;
+                mg = mm->graph;
+                for (me = agfstout(mg, mm); me; me = agnxtout(mg, me)) {
+                    mn = me->head;
+                    subg = agusergraph(mn);
+                    if (!strncmp(subg->name, "cluster", 7)) {
+                        add_cluster(g,subg);
+                        compute_bb(subg);
+                    }
+                }
+            }
 	} else {
-	    neatoLayout(g, layoutMode, model);
+	    neatoLayout(g, g, layoutMode, model);
 	    adjustNodes(g);
 	    addZ (g);
 	    spline_edges(g);
diff --git a/lib/neatogen/neatoprocs.h b/lib/neatogen/neatoprocs.h
index 93c4247d9..0b11bcb75 100644
--- a/lib/neatogen/neatoprocs.h
+++ b/lib/neatogen/neatoprocs.h
@@ -59,7 +59,7 @@ extern "C" {
     extern void neato_init_edge(Agedge_t *);
     extern void neato_init_node(Agnode_t *);
     extern void neato_layout(Agraph_t * g);
-    extern void neatoLayout(Agraph_t * g, int layoutMode, int layoutModel);
+    extern void neatoLayout(Agraph_t *mg, Agraph_t * g, int layoutMode, int layoutModel);
     extern void neato_init_graph(graph_t * g);
     extern void neato_init_graphn(graph_t * g, int);
     extern int Plegal_arrangement(Ppoly_t ** polys, int n_polys);
diff --git a/lib/neatogen/quad_prog_solve.c b/lib/neatogen/quad_prog_solve.c
index e88f8120a..0301d420c 100644
--- a/lib/neatogen/quad_prog_solve.c
+++ b/lib/neatogen/quad_prog_solve.c
@@ -389,6 +389,176 @@ constrained_majorization_new(CMajEnv *e, float * b, float **coords,
 	return counter;
 }
 
+static float *place;
+static int compare_incr(const void *a, const void *b) {
+	if(place[*(int*)a]>place[*(int*)b]) {
+		return 1;
+	} else if(place[*(int*)a]<place[*(int*)b]) {
+		return -1;
+	}
+	return 0;
+}
+/*
+While not converged: move everything towards the optimum, then satisfy constraints with as little displacement as possible.
+Returns number of iterations before convergence.
+*/
+int constrained_majorization_gradient_projection(CMajEnv *e,
+		float * b, float ** coords, int ndims, int cur_axis, int max_iterations,
+		float * hierarchy_boundaries,float levels_gap) {
+	
+	int i,j,counter;
+	if(max_iterations==0) return 0;
+
+	int * ordering=e->ordering;
+	int * levels=e->levels;
+	int num_levels=e->num_levels; 
+	bool converged=false;
+	place = coords[cur_axis];
+	float *g = e->fArray1;
+	float *old_place = e->fArray2;
+	float *d = e->fArray4;
+
+#ifdef CONMAJ_LOGGING
+	double prev_stress=0;
+	static int call_no=0;
+	for (i=0; i<e->n; i++) {
+		prev_stress += 2*b[i]*place[i];
+		for (j=0; j<e->n; j++) {
+			prev_stress -= e->A[i][j]*place[j]*place[i];
+		}
+	}
+	FILE *logfile = fopen("constrained_majorization_log","a");
+
+	fprintf(logfile,"grad proj %d: stress=%f\n",call_no,prev_stress);
+#endif
+	for (counter=0; counter<max_iterations&&!converged; counter++) {
+		converged=true;		
+		// find steepest descent direction
+		for (i=0; i<e->n; i++) {
+			old_place[i]=place[i];
+			g[i] = 2*b[i];
+			for (j=0; j<e->n; j++) {
+				g[i] -= 2*e->A[i][j]*place[j];
+			}
+		}		
+		float numerator = 0, denominator = 0, r;
+		for (i=0; i<e->n; i++) {
+			numerator += g[i]*g[i];
+			r=0;
+			for (j=0; j<e->n; j++) {
+				r += 2*e->A[i][j]*g[j];
+			}
+			denominator -= r*g[i];
+		}
+		float alpha = numerator/denominator;
+		for (i=0; i<e->n; i++) {
+			if(alpha > 0 && alpha < 1000) {
+				place[i]-=alpha*g[i];
+			}
+		}
+		if(num_levels) 
+			qsort((int *)ordering, (size_t)levels[0], sizeof(int), compare_incr );
+		//project to constraint boundary
+		for (i=0; i<num_levels; i++) {
+			int endOfLevel=i==num_levels-1?e->n:levels[i+1];
+			// ensure monotic increase in position within levels
+			qsort((int *)ordering+levels[i], (size_t)endOfLevel-levels[i], sizeof(int), compare_incr );
+			// If there are overlapping levels find offending nodes and place at average position
+			int ui=levels[i],li=ui-1;
+			int l = ordering[li--],u = ordering[ui++];
+			if(place[l]+levels_gap>place[u]) {
+				float sum=place[l]+place[u]-levels_gap*(e->lev[l]+e->lev[u]),w=2;
+				float avgPos=sum/w;
+				float pos;
+				bool finished;
+				do {
+					finished=true;
+					if(ui<endOfLevel) {
+						u=ordering[ui];
+						pos=place[u]-levels_gap*e->lev[u];
+						if(pos<avgPos) {
+							ui++; w++;
+							sum+=pos;
+							avgPos=sum/w;
+							finished=false;
+						}
+					}
+					
+					if(li>=0) {
+						l=ordering[li];
+						pos=place[l]-levels_gap*e->lev[l];
+						if(pos>avgPos) {
+							li--; w++;
+							sum+=pos;
+							avgPos=sum/w;
+							finished=false;
+						}
+					}
+				} while (!finished);
+				for(j=li+1;j<ui;j++) {
+					place[ordering[j]]=avgPos+levels_gap*e->lev[ordering[j]];
+				}
+			}
+		}
+		// set place to the intersection of old_place-g and boundary and compute d, the vector from intersection pnt to projection pnt
+		for (i=0; i<e->n; i++) {
+			d[i]=place[i]-old_place[i];
+		}	
+		// now compute beta
+		numerator = 0, denominator = 0;
+		for (i=0; i<e->n; i++) {
+			numerator += g[i]*d[i];
+			r=0;
+			for (j=0; j<e->n; j++) {
+				r += 2*e->A[i][j]*d[j];
+			}
+			denominator += r*d[i];
+		}
+		float beta = numerator/denominator;
+
+		float test=0,tmptest=0;
+		for (i=0; i<e->n; i++) {
+			if(beta>0&&beta<1.0) {
+				place[i]=old_place[i]+beta*d[i];
+			}
+			tmptest = fabs(place[i]-old_place[i]);
+			if(test<tmptest) test=tmptest;
+		}
+		computeHierarchyBoundaries(place, e->n, ordering, levels, num_levels, hierarchy_boundaries);
+		/*
+		if(num_levels) 
+			qsort((int *)ordering, (size_t)levels[0], sizeof(int), compare_incr );
+		for (i=0; i<num_levels; i++) {
+			int endOfLevel=i==num_levels-1?e->n:levels[i+1];
+			// ensure monotic increase in position within levels
+			qsort((int *)ordering+levels[i], (size_t)endOfLevel-levels[i], sizeof(int), compare_incr );
+			// If there are overlapping levels find offending nodes and place at average position
+			int l = ordering[levels[i]-1],u = ordering[levels[i]];
+			//assert(place[l]+levels_gap<=place[u]+0.00001);
+		}
+		*/
+#ifdef CONMAJ_LOGGING
+		double stress=0;
+		for (i=0; i<e->n; i++) {
+			stress += 2*b[i]*place[i];
+			for (j=0; j<e->n; j++) {
+				stress -= e->A[i][j]*place[j]*place[i];
+			}
+		}
+		fprintf(logfile,"%d: stress=%f, test=%f, %s\n",call_no,stress,test,(stress>=prev_stress)?"No Improvement":"");
+		prev_stress=stress;
+#endif
+		if(test>quad_prog_tol) {
+			converged=false;
+		}
+	}
+#ifdef CONMAJ_LOGGING
+	call_no++;
+	fclose(logfile);
+#endif
+	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,
diff --git a/lib/neatogen/quad_prog_solver.h b/lib/neatogen/quad_prog_solver.h
index e0084b2fd..fb4d894eb 100644
--- a/lib/neatogen/quad_prog_solver.h
+++ b/lib/neatogen/quad_prog_solver.h
@@ -43,6 +43,12 @@ 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 int constrained_majorization_gradient_projection(CMajEnv *e,
+	float * b, float ** coords, int ndims, int cur_axis, int max_iterations,
+	float * hierarchy_boundaries,float levels_gap);
+
 extern int constrained_majorization_new_with_gaps(CMajEnv*, float*, float**, 
                                                   int, int, int,  float*, float);
 extern void deleteCMajEnv(CMajEnv *e);
diff --git a/lib/neatogen/quad_prog_vpsc.c b/lib/neatogen/quad_prog_vpsc.c
new file mode 100644
index 000000000..704e57f75
--- /dev/null
+++ b/lib/neatogen/quad_prog_vpsc.c
@@ -0,0 +1,594 @@
+/**********************************************************
+ *
+ * Solve a quadratic function f(X) = X' e->A X + b X
+ * subject to a set of separation constraints e->cs
+ *
+ * Tim Dwyer, 2006
+ **********************************************************/
+
+#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 <csolve_VPSC.h>
+#include "quad_prog_vpsc.h"
+#include "quad_prog_solver.h"
+
+//#define CONMAJ_LOGGING 1
+#define quad_prog_tol 1e-4
+
+/*
+ * Use gradient-projection to solve Variable Placement with Separation Constraints problem.
+ */
+int 
+constrained_majorization_vpsc(CMajEnvVPSC *e, float * b, float *place, 
+                             int max_iterations)
+{
+	int i,j,counter;
+    int n=e->nv+e->nldv; // for laplacian computation need number of real vars and those
+                         // dummy vars included in lap
+	if(max_iterations==0) return 0;
+
+	bool converged=false;
+#ifdef CONMAJ_LOGGING
+	static int call_no=0;
+#endif // CONMAJ_LOGGING
+	float *g = e->fArray1;
+	float *old_place = e->fArray2;
+	float *d = e->fArray3;
+    //fprintf(stderr,"Entered: constrained_majorization_vpsc, #constraints=%d\n",e->m);
+    if(e->m>0) {
+	    for (i=0;i<n;i++) {
+		    setVariableDesiredPos(e->vs[i],place[i]);
+	    }
+        //fprintf(stderr,"  calling satisfyVPSC...\n");
+        satisfyVPSC(e->vpsc);	
+        for (i=0;i<n;i++) {
+            place[i]=getVariablePos(e->vs[i]);
+            //fprintf(stderr,"vs[%d]=%f\n",i,place[i]);
+        }
+        //fprintf(stderr,"    done.\n");
+    }
+#ifdef CONMAJ_LOGGING
+	float prev_stress=0;
+	for (i=0; i<n; i++) {
+		prev_stress += 2*b[i]*place[i];
+		for (j=0; j<n; j++) {
+			prev_stress -= e->A[i][j]*place[j]*place[i];
+		}
+	}
+	FILE *logfile = fopen("constrained_majorization_log","a");
+
+	//fprintf(logfile,"grad proj %d: stress=%f\n",call_no,prev_stress);
+#endif
+	
+	for (counter=0; counter<max_iterations&&!converged; counter++) {
+        //fprintf(stderr,".");
+		converged=true;		
+		// find steepest descent direction
+		for (i=0; i<n; i++) {
+			old_place[i]=place[i];
+			g[i] = 2*b[i];
+			for (j=0; j<n; j++) {
+				g[i] -= 2*e->A[i][j]*place[j];
+			}
+		}		
+		float numerator = 0, denominator = 0, r;
+		for (i=0; i<n; i++) {
+			numerator += g[i]*g[i];
+			r=0;
+			for (j=0; j<n; j++) {
+				r += 2*e->A[i][j]*g[j];
+			}
+			denominator -= r*g[i];
+		}
+		float alpha = numerator/denominator;
+		for (i=0; i<n; i++) {
+			place[i]-=alpha*g[i];
+		}
+        if(e->m>0) {
+            //project to constraint boundary
+            for (i=0;i<n;i++) {
+                setVariableDesiredPos(e->vs[i],place[i]);
+            }
+            satisfyVPSC(e->vpsc);
+            for (i=0;i<n;i++) {
+                place[i]=getVariablePos(e->vs[i]);
+            }
+        }
+		// set place to the intersection of old_place-g and boundary and compute d, the vector from intersection pnt to projection pnt
+		for (i=0; i<n; i++) {
+			d[i]=place[i]-old_place[i];
+		}	
+		// now compute beta
+		numerator = 0, denominator = 0;
+		for (i=0; i<n; i++) {
+			numerator += g[i]*d[i];
+			r=0;
+			for (j=0; j<n; j++) {
+				r += 2*e->A[i][j]*d[j];
+			}
+			denominator += r*d[i];
+		}
+		float beta = numerator/denominator;
+
+		float test=0;
+		for (i=0; i<n; i++) {
+			// beta > 1.0 takes us back outside the feasible region
+			// beta < 0 clearly not useful and may happen due to numerical imp.
+			if(beta>0&&beta<1.0) {
+				place[i]=old_place[i]+beta*d[i];
+			}
+			test+= fabs(place[i]-old_place[i]);
+		}
+#ifdef CONMAJ_LOGGING
+		float stress=0;
+		for (i=0; i<n; i++) {
+			stress += 2*b[i]*place[i];
+			for (j=0; j<n; j++) {
+				stress -= e->A[i][j]*place[j]*place[i];
+			}
+		}
+		fprintf(logfile,"%d: stress=%f, test=%f, %s\n",call_no,stress,test,(stress>=prev_stress)?"No Improvement":"");
+		prev_stress=stress;
+#endif
+		if(test>quad_prog_tol) {
+			converged=false;
+		}
+	}
+#ifdef CONMAJ_LOGGING
+	call_no++;
+	fclose(logfile);
+#endif
+	return counter;
+}
+
+/*
+ * Set up environment and global constraints (dir-edge constraints, containment constraints
+ * etc).
+ *
+ * diredges: 0=no dir edge constraints
+ *           1=one separation constraint for each edge (in acyclic subgraph)
+ *           2=DiG-CoLa level constraints
+ */
+CMajEnvVPSC*
+initCMajVPSC(int n, float* packedMat, vtx_data* graph,
+        ipsep_options *opt, int diredges)
+{
+    int i,j;
+    // nv is the number of real nodes
+    int nConCs;
+    //fprintf(stderr,"Entered initCMajVPSC\n");
+	CMajEnvVPSC *e = GNEW(CMajEnvVPSC);
+	e->A=NULL;
+    e->packedMat=packedMat;
+    /* if we have clusters then we'll need two constraints for each var in
+     * a cluster */
+    e->nldv = 2*opt->clusters->nclusters;
+    e->nv = n - e->nldv;
+    e->ndv = 0;
+
+    e->gcs=NULL;
+    e->vs=N_GNEW(n, Variable*);
+    for(i=0;i<n;i++) {
+        e->vs[i]=newVariable(i,1.0,1.0);
+    }
+    e->gm=0;
+    if(diredges==1) {
+        if(Verbose)
+            fprintf(stderr,"  generate edge constraints...\n");
+        for(i=0;i<e->nv;i++) {
+            for(j=1;j<graph[i].nedges;j++) {
+                //fprintf(stderr,"edist=%f\n",graph[i].edists[j]);
+                if(graph[i].edists[j]>0.01) {
+                    e->gm++;
+                }
+            }
+        }
+        e->gcs=newConstraints(e->gm);
+        e->gm=0;
+        for(i=0;i<e->nv;i++) {
+            for(j=1;j<graph[i].nedges;j++) {
+                int u=i,v=graph[i].edges[j];
+                if(graph[i].edists[j]>0) {
+                    e->gcs[e->gm++]=newConstraint(e->vs[u],e->vs[v],opt->edge_gap);
+                }
+            }
+        }
+    } else if(diredges==2) {
+	    int *ordering=NULL, *ls=NULL, cvar;
+        DigColaLevel* levels;
+        Variable** vs = e->vs;
+        // e->ndv is the number of dummy variables required, one for each boundary
+		compute_hierarchy(graph, n, 1e-2, 1e-1, NULL, &ordering, &ls, &e->ndv);		
+        levels=assign_digcola_levels(ordering,n,ls,e->ndv);
+        if(Verbose)
+            fprintf(stderr,"Found %d DiG-CoLa boundaries\n",e->ndv);
+        e->gm=get_num_digcola_constraints(levels,e->ndv+1)+e->ndv-1;
+        e->gcs=newConstraints(e->gm);
+        e->gm=0;
+        e->vs = N_GNEW(n+e->ndv,Variable*);
+        for(i=0;i<n;i++) {
+            e->vs[i]=vs[i];
+        }
+        free(vs);
+        // create dummy vars
+        for(i=0;i<e->ndv;i++) {
+            //  dummy vars should have 0 weight
+            cvar=n+i;
+            e->vs[cvar]=newVariable(cvar,1.0,0.000001);
+        }
+        double halfgap=opt->edge_gap;
+        for(i=0;i<e->ndv;i++) {
+            cvar=n+i;
+            // outgoing constraints for each var in level below boundary
+            for(j=0;j<levels[i].num_nodes;j++) {
+                e->gcs[e->gm++]=newConstraint(e->vs[levels[i].nodes[j]],e->vs[cvar],halfgap);
+            }
+            // incoming constraints for each var in level above boundary
+            for(j=0;j<levels[i+1].num_nodes;j++) {
+                e->gcs[e->gm++]=newConstraint(e->vs[cvar],e->vs[levels[i+1].nodes[j]],halfgap);
+            }
+        }
+        // constraints between adjacent boundary dummy vars
+        for(i=0;i<e->ndv-1;i++) {
+            e->gcs[e->gm++]=newConstraint(e->vs[n+i],e->vs[n+i+1],0);
+        }
+    }
+    //fprintf(stderr,"  generate edge constraints... done: n=%d,m=%d\n",e->n,e->gm);
+    if(opt->clusters->nclusters>0) {
+        //fprintf(stderr,"  generate cluster containment constraints...\n");
+        Constraint** ecs=e->gcs;
+        nConCs=2*opt->clusters->nvars;
+        e->gcs=newConstraints(e->gm+nConCs);
+        for(i=0;i<e->gm;i++) {
+            e->gcs[i]=ecs[i];
+        }
+        if(ecs!=NULL) deleteConstraints(0,ecs);
+        for(i=0;i<opt->clusters->nclusters;i++) {
+            for(j=0;j<opt->clusters->clustersizes[i];j++) {
+                Variable* v=e->vs[opt->clusters->clusters[i][j]];
+                Variable* cl=e->vs[e->nv+2*i];
+                Variable* cr=e->vs[e->nv+2*i+1];
+                e->gcs[e->gm++]=newConstraint(cl,v,0);
+                e->gcs[e->gm++]=newConstraint(v,cr,0);
+            }
+        }
+        //fprintf(stderr,"  containment constraints... done: \n");
+    }
+
+    e->m=0;
+    e->cs=NULL;
+    if(e->gm>0) {
+        e->vpsc = newIncVPSC(n+e->ndv,e->vs,e->gm,e->gcs);
+        e->m=e->gm;
+        e->cs=e->gcs;
+    }
+    if(packedMat!=NULL) {
+	    e->A = unpackMatrix(packedMat,n);
+    }
+#ifdef MOSEK
+    e->mosekEnv=NULL;
+    if(opt->mosek) {
+        e->mosekEnv = mosek_init_sep(e->packedMat,n,e->ndv,e->gcs,e->gm);
+    }
+#endif
+
+	e->fArray1 = N_GNEW(n,float);
+	e->fArray2 = N_GNEW(n,float);
+	e->fArray3 = N_GNEW(n,float);
+    if(Verbose)
+        fprintf(stderr,"  initCMajVPSC done: %d global constraints generated.\n",e->m);
+	return e;
+}
+void
+deleteCMajEnvVPSC(CMajEnvVPSC *e) 
+{
+    int i;
+    if(e->A!=NULL) {
+        free (e->A[0]); 
+        free (e->A);
+    }
+    if(e->m>0) {
+        deleteVPSC(e->vpsc);
+        if(e->cs!=e->gcs && e->gcs!=NULL) deleteConstraints(0,e->gcs);
+        deleteConstraints(e->m,e->cs);
+        for(i=0;i<e->nv+e->nldv+e->ndv;i++) {
+            deleteVariable(e->vs[i]);
+        }
+        free (e->vs);
+    }
+	free (e->fArray1);
+	free (e->fArray2);
+	free (e->fArray3);
+#ifdef MOSEK
+    if(e->mosekEnv) {
+        mosek_delete(e->mosekEnv);
+    }
+#endif //MOSEK
+	free (e);
+}
+
+// generate non-overlap constraints inside each cluster, including dummy
+// nodes at bounds of cluster
+// generate constraints again for top level nodes and clusters treating
+// clusters as rectangles of dim (l,r,b,t)
+// for each cluster map in-constraints to l out-constraints to r 
+//
+// For now, we'll keep the global containment constraints already
+// generated for each cluster, and simply generate non-overlap constraints
+// for all nodes and then an additional set of non-overlap constraints for
+// clusters that we'll map back to the dummy vars as above.
+void generateNonoverlapConstraints(
+        CMajEnvVPSC* e,
+        float nsizeScale,
+        float** coords,
+        int k,
+        bool transitiveClosure,
+	    ipsep_options* opt
+) {
+    Constraint **csol, **csolptr;
+    int i,j,mol=0;
+    int n=e->nv+e->nldv;
+    boxf bb[n];
+    bool genclusters=opt->clusters->nclusters>0;
+    if(genclusters) {
+        // n is the number of real variables, not dummy cluster vars
+        n-=2*opt->clusters->nclusters;
+    }
+    if(k==0) {
+        // grow a bit in the x dimension, so that if overlap is resolved
+        // horizontally then it won't be considered overlapping vertically
+        nsizeScale*=1.0001;
+    }
+    for(i=0;i<n;i++) {
+        bb[i].LL.x=coords[0][i]-nsizeScale*opt->nsize[i].x/2.0-opt->gap.x/2.0; 
+        bb[i].UR.x=coords[0][i]+nsizeScale*opt->nsize[i].x/2.0+opt->gap.x/2.0; 
+        bb[i].LL.y=coords[1][i]-nsizeScale*opt->nsize[i].y/2.0-opt->gap.y/2.0; 
+        bb[i].UR.y=coords[1][i]+nsizeScale*opt->nsize[i].y/2.0+opt->gap.y/2.0; 
+    }
+    if(genclusters) {
+        Constraint **cscl[opt->clusters->nclusters+1];
+        int cm[opt->clusters->nclusters+1];
+        for(i=0;i<opt->clusters->nclusters;i++) {
+            int cn=opt->clusters->clustersizes[i];
+            Variable* cvs[cn+2];
+            boxf cbb[cn+2];
+            // compute cluster bounding bb
+            boxf container;
+            container.LL.x=container.LL.y=DBL_MAX;
+            container.UR.x=container.UR.y=-DBL_MAX;
+            for(j=0;j<cn;j++) {
+                int iv=opt->clusters->clusters[i][j];
+                cvs[j]=e->vs[iv];
+                B2BF(bb[iv],cbb[j]);
+                EXPANDBB(container,bb[iv]);
+            }
+            B2BF(container,opt->clusters->bb[i]);
+            cvs[cn]=e->vs[n+2*i];
+            cvs[cn+1]=e->vs[n+2*i+1];
+            B2BF(container,cbb[cn]);
+            B2BF(container,cbb[cn+1]);
+            if(k==0) {
+                cbb[cn].UR.x=container.LL.x+0.0001;
+                cbb[cn+1].LL.x=container.UR.x-0.0001;
+                cm[i] = genXConstraints(cn+2,cbb,cvs,&cscl[i],transitiveClosure);
+            } else {
+                cbb[cn].UR.y=container.LL.y+0.0001;
+                cbb[cn+1].LL.y=container.UR.y-0.0001;
+                cm[i] = genYConstraints(cn+2,cbb,cvs,&cscl[i]);
+            }
+            mol+=cm[i];
+        }
+        // generate top level constraints
+        {
+            int cn=opt->clusters->ntoplevel+opt->clusters->nclusters;
+            Variable* cvs[cn];
+            boxf cbb[cn];
+            for(i=0;i<opt->clusters->ntoplevel;i++) {
+                int iv=opt->clusters->toplevel[i];
+                cvs[i]=e->vs[iv];
+                B2BF(bb[iv],cbb[i]);
+            }
+            // make dummy variables for clusters
+            for(i=opt->clusters->ntoplevel;i<cn;i++) {
+                cvs[i]=newVariable(123+i,1,1);
+                j=i-opt->clusters->ntoplevel;
+                B2BF(opt->clusters->bb[j],cbb[i]);
+            }
+            i=opt->clusters->nclusters;
+            if(k==0) {
+                cm[i] = genXConstraints(cn,cbb,cvs,&cscl[i],transitiveClosure);
+            } else {
+                cm[i] = genYConstraints(cn,cbb,cvs,&cscl[i]);
+            }
+            // remap constraints from tmp dummy vars to cluster l and r vars
+            for(i=opt->clusters->ntoplevel;i<cn;i++) {
+                j=i-opt->clusters->ntoplevel;
+                // dgap is the change in required constraint gap.
+                // since we are going from a source rectangle the size
+                // of the cluster bounding box to a zero width (in x dim,
+                // zero height in y dim) rectangle, the change will be
+                // half the bb width.
+                double dgap;
+                if(k==0) {
+                    dgap=-(cbb[i].UR.x-cbb[i].LL.x)/2.0;
+                } else {
+                    dgap=-(cbb[i].UR.y-cbb[i].LL.y)/2.0;
+                }
+                remapInConstraints(cvs[i],e->vs[n+2*j],dgap);
+                remapOutConstraints(cvs[i],e->vs[n+2*j+1],dgap);
+                // there may be problems with cycles between
+                // cluster non-overlap and diredge constraints,
+                // to resolve:
+                // 
+                // for each constraint c:v->cvs[i]:
+                //   if exists diredge constraint u->v where u in c:
+                //     remap v->cl to cr->v (gap = height(v)/2)
+                //
+                // in=getInConstraints(cvs[i])
+                // for(c : in) {
+                //   assert(c.right==cvs[i]);
+                //   vin=getOutConstraints(v=c.left)
+                //   for(d : vin) {
+                //     if(d.left.cluster==i):
+                //       tmp=d.left
+                //       d.left=d.right
+                //       d.right=tmp
+                //       d.gap==height(d.right)/2
+                //   }
+                // }
+                //       
+                deleteVariable(cvs[i]);
+            }
+            mol+=cm[opt->clusters->nclusters];
+        }
+        csolptr=csol=newConstraints(mol);
+        for(i=0;i<opt->clusters->nclusters+1;i++) {
+            // copy constraints into csol
+            for(j=0;j<cm[i];j++) {
+                *csolptr++=cscl[i][j];
+            }
+            deleteConstraints(0,cscl[i]);
+        }
+    } else {
+        if(k==0) {
+            mol = genXConstraints(n,bb,e->vs,&csol,transitiveClosure);
+        } else {
+            mol = genYConstraints(n,bb,e->vs,&csol);
+        }
+    }
+    // remove constraints from previous iteration
+    if(e->m>0) {
+        // can't reuse instance of VPSC when constraints change!
+        deleteVPSC(e->vpsc);
+        for(i=e->gm==0?0:e->gm;i<e->m;i++) {
+            // delete previous overlap constraints
+            deleteConstraint(e->cs[i]);
+        }
+        // just delete the array, not the elements
+        if(e->cs!=e->gcs) deleteConstraints(0,e->cs);
+    }
+    // if we have no global constraints then the overlap constraints 
+    // are all we have to worry about.
+    // Otherwise, we have to copy the global and overlap constraints 
+    // into the one array
+    if(e->gm==0) {
+        e->m = mol;
+        e->cs = csol;
+    } else {
+        e->m = mol + e->gm;
+        e->cs = newConstraints(e->m);
+        for(i=0;i<e->m;i++) {
+            if(i<e->gm) {
+                e->cs[i]=e->gcs[i];
+            } else {
+                e->cs[i]=csol[i-e->gm];
+            }
+        }
+        // just delete the array, not the elements
+        deleteConstraints(0,csol);
+    }
+    if(Verbose)
+        fprintf(stderr,"  generated %d constraints\n",e->m);
+    e->vpsc = newIncVPSC(e->nv+e->nldv+e->ndv,e->vs,e->m,e->cs);
+#ifdef MOSEK
+    if(opt->mosek) {
+        if(e->mosekEnv!=NULL) {
+            mosek_delete(e->mosekEnv);
+        }
+        e->mosekEnv = mosek_init_sep(e->packedMat,e->nv+e->nldv,e->ndv,e->cs,e->m);
+    }
+#endif
+}
+
+/*
+ * Statically remove overlaps, that is remove all overlaps by moving each node as
+ * little as possible.
+ */
+void removeoverlaps(int n,float** coords, ipsep_options* opt) {
+    int i;
+	CMajEnvVPSC *e = initCMajVPSC(n,NULL,NULL,opt,0);
+    generateNonoverlapConstraints(e,1.0,coords,0,true,opt);
+    solveVPSC(e->vpsc);
+    for(i=0;i<n;i++) {
+        coords[0][i]=getVariablePos(e->vs[i]);
+    }
+    generateNonoverlapConstraints(e,1.0,coords,1,false,opt);
+    solveVPSC(e->vpsc);
+    for(i=0;i<n;i++) {
+        coords[1][i]=getVariablePos(e->vs[i]);
+    }
+    deleteCMajEnvVPSC(e);
+}
+
+/*
+ unpack the "ordering" array into an array of DigColaLevel
+*/
+DigColaLevel* assign_digcola_levels(int *ordering, int n, int *level_inds, int num_divisions) {
+	int i,j;
+	DigColaLevel *l=N_GNEW(num_divisions+1,DigColaLevel);
+	// first level
+	l[0].num_nodes=level_inds[0];
+	l[0].nodes=N_GNEW(l[0].num_nodes,int);
+	for(i=0;i<l[0].num_nodes;i++) {
+		l[0].nodes[i]=ordering[i];
+	}
+	// second through second last level
+	for(i=1;i<num_divisions;i++) {
+		l[i].num_nodes=level_inds[i]-level_inds[i-1];
+		l[i].nodes=N_GNEW(l[i].num_nodes,int);
+		for(j=0;j<l[i].num_nodes;j++) {
+			l[i].nodes[j]=ordering[level_inds[i-1]+j];
+		}
+	}
+	// last level
+	if (num_divisions>0) {
+		l[num_divisions].num_nodes=n-level_inds[num_divisions-1];
+		l[num_divisions].nodes=N_GNEW(l[num_divisions].num_nodes,int);
+		for(i=0;i<l[num_divisions].num_nodes;i++) {
+			l[num_divisions].nodes[i]=ordering[level_inds[num_divisions-1]+i];
+		}
+	}
+	return l;
+}
+void delete_digcola_levels(DigColaLevel *l, int num_levels) {
+	int i;
+	for(i=0;i<num_levels;i++) {
+		free(l[i].nodes);
+	}
+	free(l);
+}
+void print_digcola_levels(FILE *logfile, DigColaLevel *levels, int num_levels) {
+	fprintf(logfile,"levels:\n");
+	int i,j;
+	for(i=0;i<num_levels;i++) {
+		fprintf(logfile,"  l[%d]:",i);
+		for(j=0;j<levels[i].num_nodes;j++) {
+			fprintf(logfile,"%d ",levels[i].nodes[j]);
+		}
+		fprintf(logfile,"\n");
+	}
+}
+/*********************
+get number of separation constraints based on the number of nodes in each level
+ie, num_sep_constraints = sum_i^{num_levels-1} (|L[i]|+|L[i+1]|)
+**********************/
+int get_num_digcola_constraints(DigColaLevel *levels, int num_levels) {
+	int i,nc=0;
+	for(i=1;i<num_levels;i++) {
+		nc+=levels[i].num_nodes+levels[i-1].num_nodes;
+	}
+	nc+=levels[0].num_nodes+levels[num_levels-1].num_nodes;
+	return nc;
+}
+
+#endif /* DIGCOLA */
+
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
diff --git a/lib/neatogen/quad_prog_vpsc.h b/lib/neatogen/quad_prog_vpsc.h
new file mode 100644
index 000000000..fc530ffe6
--- /dev/null
+++ b/lib/neatogen/quad_prog_vpsc.h
@@ -0,0 +1,77 @@
+/**********************************************************
+*      Written by Tim Dwyer for the graphviz package      *
+*                http://www.graphviz.org/                 *
+*                                                         *
+**********************************************************/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#ifndef _QUAD_PROG_VPSC_H_
+#define _QUAD_PROG_VPSC_H_
+
+#ifdef DIGCOLA
+
+#include "defs.h"
+#include "digcola.h"
+#ifdef MOSEK
+#include "mosek_quad_solve.h"
+#endif //MOSEK
+
+typedef struct CMajEnvVPSC {
+	float **A;
+	float *packedMat;
+	int nv;   // number of actual vars
+	int nldv; // number of dummy nodes included in lap matrix
+	int ndv;  // number of dummy nodes not included in lap matrix
+	Variable **vs;
+	int m; /* total number of constraints for next iteration */
+	int gm; /* number of global constraints */
+	Constraint **cs;
+	/* global constraints are persistant throughout optimisation process */
+	Constraint **gcs;
+	VPSC *vpsc;
+	float *fArray1; /* utility arrays - reusable memory */
+	float *fArray2;
+	float *fArray3;
+#ifdef MOSEK
+	MosekEnv *mosekEnv;
+#endif // MOSEK
+} CMajEnvVPSC;
+
+extern CMajEnvVPSC* initCMajVPSC(int n, float *packedMat, vtx_data* graph, ipsep_options *opt, int diredges);
+
+extern int constrained_majorization_vpsc(CMajEnvVPSC*, float*, float*, int);
+
+extern void deleteCMajEnvVPSC(CMajEnvVPSC *e);
+extern void generateNonoverlapConstraints(
+        CMajEnvVPSC* e,
+        float nsizeScale,
+        float** coords,
+        int k,
+	bool transitiveClosure,
+	ipsep_options* opt
+);
+
+extern void removeoverlaps(int,float**,ipsep_options*);
+
+typedef struct {
+	int *nodes;
+	int num_nodes;
+} DigColaLevel;
+
+/*
+ unpack the "ordering" array into an array of DigColaLevel (as defined above)
+*/
+extern DigColaLevel* assign_digcola_levels(int *ordering, int n, int *level_inds, int num_divisions);
+extern void delete_digcola_levels(DigColaLevel *l, int num_levels);
+extern void print_digcola_levels(FILE* logfile, DigColaLevel *levels, int num_levels);
+int get_num_digcola_constraints(DigColaLevel *levels, int num_levels);
+#endif 
+
+#endif /* _QUAD_PROG_VPSC_H_ */
+
+#ifdef __cplusplus
+}
+#endif