Previously graph layout was done using the Kamada-Kawai layout algorithm

implemented in Boost.  I am replacing this with a custom implementation of
a constrained stress-majorization algorithm.

The stress-majorization algorithm is more robust and has better convergence
characteristics than Kamada-Kawai, and also simple constraints can be placed
on node position (for example, to enforce downward-pointing edges, non-overlap constraints, or cluster constraints).

Another big advantage is that we no longer need Boost.

I've tested the basic functionality, but I have yet to properly handle
disconnected graphs or to properly scale the resulting layout.

This commit also includes significant refactoring... the quadratic program solver - libvpsc (Variable Placement with Separation Constraints) has been moved to src/libvpsc and the actual graph layout algorithm is in libcola.

(bzr r1394)

14 files changed, 2229 insertions, 0 deletions

diff --git a/src/libcola/Makefile_insert b/src/libcola/Makefile_insert
new file mode 100644
index 000000000..f8f9de20c
--- /dev/null
+++ b/src/libcola/Makefile_insert
@@ -0,0 +1,17 @@
+## Makefile.am fragment sourced by src/Makefile.am.
+
+libcola/all: libcola.a
+
+libcola/clean:
+	rm -f libcola/libcola.a $(libcola_libcola_a_OBJECTS)
+
+libcola_libcola_a_SOURCES = libcola/cola.h\
+	libcola/cola.cpp\
+	libcola/conjugate_gradient.cpp\
+	libcola/conjugate_gradient.h\
+	libcola/gradient_projection.cpp\
+	libcola/gradient_projection.h\
+	libcola/shortest_paths.cpp\
+	libcola/shortest_paths.h\
+	libcola/straightener.h\
+	libcola/straightener.cpp
diff --git a/src/libcola/cola.cpp b/src/libcola/cola.cpp
new file mode 100644
index 000000000..74663f501
--- /dev/null
+++ b/src/libcola/cola.cpp
@@ -0,0 +1,299 @@
+#include "cola.h"
+#include "conjugate_gradient.h"
+#include "straightener.h"
+
+namespace cola {
+
+/**
+ * Find the euclidean distance between a pair of dummy variables
+ */
+inline double dummy_var_euclidean_dist(GradientProjection* gpx, GradientProjection* gpy, unsigned i) {
+    double dx = gpx->dummy_vars[i]->place_r - gpx->dummy_vars[i]->place_l,
+        dy = gpy->dummy_vars[i]->place_r - gpy->dummy_vars[i]->place_l;
+    return sqrt(dx*dx + dy*dy);
+}
+
+void 
+ConstrainedMajorizationLayout
+::setupDummyVars() {
+    if(clusters==NULL) return;
+    double* coords[2]={X,Y};
+    GradientProjection* gp[2]={gpX,gpY};
+    for(unsigned k=0;k<2;k++) {
+        gp[k]->clearDummyVars();
+        if(clusters) {
+            for(Clusters::iterator cit=clusters->begin();
+                    cit!=clusters->end(); ++cit) {
+                Cluster *c = *cit;
+                DummyVarPair* p = new DummyVarPair(edge_length);
+                gp[k]->dummy_vars.push_back(p);
+                double minPos=DBL_MAX, maxPos=-DBL_MAX;
+                for(Cluster::iterator vit=c->begin();
+                        vit!=c->end(); ++vit) {
+                    double pos = coords[k][*vit];
+                    minPos=min(pos,minPos);
+                    maxPos=max(pos,maxPos);
+                    p->leftof.push_back(make_pair(*vit,0)); 
+                    p->rightof.push_back(make_pair(*vit,0)); 
+                }
+                p->place_l = minPos;
+                p->place_r = maxPos;
+            }
+        }
+    }
+    for(unsigned k=0;k<2;k++) {
+        unsigned n_d = gp[k]->dummy_vars.size();
+        if(n_d > 0) {
+            for(unsigned i=0; i<n_d; i++) {
+                gp[k]->dummy_vars[i]->computeLinearTerm(dummy_var_euclidean_dist(gpX, gpY, i));
+            }
+        }
+    }
+}
+void ConstrainedMajorizationLayout::majlayout(
+        double** Dij, GradientProjection* gp, double* coords) 
+{
+    double b[n];
+    fill(b,b+n,0);
+    majlayout(Dij,gp,coords,b);
+}
+void ConstrainedMajorizationLayout::majlayout(
+        double** Dij, GradientProjection* gp, double* coords, double* b) 
+{
+    double L_ij,dist_ij,degree;
+    /* compute the vector b */
+    /* multiply on-the-fly with distance-based laplacian */
+    for (unsigned i = 0; i < n; i++) {
+        degree = 0;
+        if(i<lapSize) {
+            for (unsigned j = 0; j < lapSize; j++) {
+                if (j == i) continue;
+                dist_ij = euclidean_distance(i, j);
+                if (dist_ij > 1e-30 && Dij[i][j] > 1e-30) {	/* skip zero distances */
+                    /* calculate L_ij := w_{ij}*d_{ij}/dist_{ij} */
+                    L_ij = 1.0 / (dist_ij * Dij[i][j]);
+                    degree -= L_ij;
+                    b[i] += L_ij * coords[j];
+                }
+            }
+            b[i] += degree * coords[i];
+        }
+        assert(!isnan(b[i]));
+    }
+    if(constrainedLayout) {
+        setupDummyVars();
+        gp->solve(b);
+    } else {
+        conjugate_gradient(lap2, coords, b, n, tol, n);
+    }
+    moveBoundingBoxes();
+}
+inline double ConstrainedMajorizationLayout
+::compute_stress(double **Dij) {
+    double sum = 0, d, diff;
+    for (unsigned i = 1; i < lapSize; i++) {
+        for (unsigned j = 0; j < i; j++) {
+            d = Dij[i][j];
+            diff = d - euclidean_distance(i,j);
+            sum += diff*diff / (d*d);
+        }
+    }
+    if(clusters!=NULL) {
+        for(unsigned i=0; i<gpX->dummy_vars.size(); i++) {
+            sum += gpX->dummy_vars[i]->stress(dummy_var_euclidean_dist(gpX, gpY, i));
+        }
+    }
+    return sum;
+}
+/*
+void ConstrainedMajorizationLayout
+::addLinearConstraints(LinearConstraints* linearConstraints) {
+    n=lapSize+linearConstraints->size();
+    Q=new double*[n];
+    X=new double[n];
+    Y=new double[n];
+    for(unsigned i = 0; i<n; i++) {
+        X[i]=rs[i]->getCentreX();
+        Y[i]=rs[i]->getCentreY();
+        Q[i]=new double[n];
+        for(unsigned j=0; j<n; j++) {
+            if(i<lapSize&&j<lapSize) {
+                Q[i][j]=lap2[i][j];
+            } else {
+                Q[i][j]=0;
+            }
+        }
+    }
+    for(LinearConstraints::iterator i=linearConstraints->begin();
+           i!= linearConstraints->end();i++) {
+        LinearConstraint* c=*i;
+        Q[c->u][c->u]+=c->w*c->duu;
+        Q[c->u][c->v]+=c->w*c->duv;
+        Q[c->u][c->b]+=c->w*c->dub;
+        Q[c->v][c->u]+=c->w*c->duv;
+        Q[c->v][c->v]+=c->w*c->dvv;
+        Q[c->v][c->b]+=c->w*c->dvb;
+        Q[c->b][c->b]+=c->w*c->dbb;
+        Q[c->b][c->u]+=c->w*c->dub;
+        Q[c->b][c->v]+=c->w*c->dvb;
+    }
+}
+*/
+
+bool ConstrainedMajorizationLayout::run() {
+    /*
+    for(unsigned i=0;i<n;i++) {
+        for(unsigned j=0;j<n;j++) {
+            cout << lap2[i][j] << " ";
+        }
+        cout << endl;
+    }
+    */
+    do {
+        /* Axis-by-axis optimization: */
+        if(straightenEdges) {
+            straighten(*straightenEdges,HORIZONTAL);
+            straighten(*straightenEdges,VERTICAL);
+        } else {
+            majlayout(Dij,gpX,X);
+            majlayout(Dij,gpY,Y);
+        }
+    } while(!done(compute_stress(Dij),X,Y));
+    return true;
+}
+static bool straightenToProjection=true;
+void ConstrainedMajorizationLayout::straighten(vector<straightener::Edge*>& sedges, Dim dim) {
+	vector<straightener::Node*> snodes;
+	for (unsigned i=0;i<lapSize;i++) {
+		snodes.push_back(new straightener::Node(i,boundingBoxes[i]));
+	}
+    SimpleConstraints cs;
+    straightener::generateConstraints(snodes,sedges,cs,dim);
+    n=snodes.size();
+    Q=new double*[n];
+    delete [] X;
+    delete [] Y;
+    X=new double[n];
+    Y=new double[n];
+    for(unsigned i = 0; i<n; i++) {
+        X[i]=snodes[i]->x;
+        Y[i]=snodes[i]->y;
+        Q[i]=new double[n];
+        for(unsigned j=0; j<n; j++) {
+            if(i<lapSize&&j<lapSize) {
+                Q[i][j]=lap2[i][j];
+            } else {
+                Q[i][j]=0;
+            }
+        }
+    }
+    LinearConstraints linearConstraints;
+    for(unsigned i=0;i<sedges.size();i++) {
+        sedges[i]->nodePath(snodes);
+        vector<unsigned>& path=sedges[i]->path;
+        // take u and v as the ends of the line
+        //unsigned u=path[0];
+        //unsigned v=path[path.size()-1];
+        double total_length=0;
+        for(unsigned j=1;j<path.size();j++) {
+            unsigned u=path[j-1], v=path[j];
+            total_length+=euclidean_distance(u,v);
+        }
+        for(unsigned j=1;j<path.size()-1;j++) {
+            // find new u and v for each line segment
+            unsigned u=path[j-1];
+            unsigned b=path[j];
+            unsigned v=path[j+1];
+            double weight=-0.01;
+            double wub=euclidean_distance(u,b)/total_length;
+            double wbv=euclidean_distance(b,v)/total_length;
+            linearConstraints.push_back(new cola::LinearConstraint(u,v,b,weight,wub,wbv,X,Y));
+        }
+    }
+    //cout << "Generated "<<linearConstraints.size()<< " linear constraints"<<endl;
+    assert(snodes.size()==lapSize+linearConstraints.size());
+    double b[n],*coords=dim==HORIZONTAL?X:Y,dist_ub,dist_bv;
+    fill(b,b+n,0);
+    for(LinearConstraints::iterator i=linearConstraints.begin();
+           i!= linearConstraints.end();i++) {
+        LinearConstraint* c=*i;
+        if(straightenToProjection) {
+            Q[c->u][c->u]+=c->w*c->duu;
+            Q[c->u][c->v]+=c->w*c->duv;
+            Q[c->u][c->b]+=c->w*c->dub;
+            Q[c->v][c->u]+=c->w*c->duv;
+            Q[c->v][c->v]+=c->w*c->dvv;
+            Q[c->v][c->b]+=c->w*c->dvb;
+            Q[c->b][c->b]+=c->w*c->dbb;
+            Q[c->b][c->u]+=c->w*c->dub;
+            Q[c->b][c->v]+=c->w*c->dvb;
+        } else {
+            double wub=edge_length*c->frac_ub;
+            double wbv=edge_length*c->frac_bv;
+            dist_ub=euclidean_distance(c->u,c->b)*wub;
+            dist_bv=euclidean_distance(c->b,c->v)*wbv;
+            wub=max(wub,0.00001);
+            wbv=max(wbv,0.00001);
+            dist_ub=max(dist_ub,0.00001);
+            dist_bv=max(dist_bv,0.00001);
+            wub=1/(wub*wub);
+            wbv=1/(wbv*wbv);
+            Q[c->u][c->u]-=wub;
+            Q[c->u][c->b]+=wub;
+            Q[c->v][c->v]-=wbv;
+            Q[c->v][c->b]+=wbv;
+            Q[c->b][c->b]-=wbv+wub;
+            Q[c->b][c->u]+=wub;
+            Q[c->b][c->v]+=wbv;
+
+            b[c->u]+=(coords[c->b]-coords[c->u]) / dist_ub;
+            b[c->v]+=(coords[c->b]-coords[c->v]) / dist_bv;
+            b[c->b]+=coords[c->u] / dist_ub + coords[c->v] / dist_bv
+                   - coords[c->b] / dist_ub - coords[c->b] / dist_bv;
+        }
+    }
+	GradientProjection gp(dim,n,Q,coords,tol,100,
+            (AlignmentConstraints*)NULL,false,(Rectangle**)NULL,(PageBoundaryConstraints*)NULL,&cs);
+    constrainedLayout = true;
+    majlayout(Dij,&gp,coords,b);
+    for(unsigned i=0;i<sedges.size();i++) {
+        sedges[i]->createRouteFromPath(X,Y);
+        sedges[i]->dummyNodes.clear();
+        sedges[i]->path.clear();
+    }
+    for(unsigned i=0;i<cs.size();i++) {
+        delete cs[i];
+    }
+    for(unsigned i=0;i<linearConstraints.size();i++) {
+        delete linearConstraints[i];
+    }
+    for(unsigned i=0;i<snodes.size();i++) {
+        delete snodes[i];
+    }
+    for(unsigned i = 0; i<n; i++) {
+        delete [] Q[i];
+    }
+    delete [] Q;
+    snodes.resize(lapSize);
+}
+
+void ConstrainedMajorizationLayout::setupConstraints(
+        AlignmentConstraints* acsx, AlignmentConstraints* acsy,
+        bool avoidOverlaps, 
+        PageBoundaryConstraints* pbcx, PageBoundaryConstraints* pbcy,
+        SimpleConstraints* scx, SimpleConstraints* scy,
+        Clusters* cs,
+        vector<straightener::Edge*>* straightenEdges) {
+    constrainedLayout = true;
+    this->avoidOverlaps = avoidOverlaps;
+    if(cs) {
+        clusters=cs;
+    }
+	gpX=new GradientProjection(
+            HORIZONTAL,n,Q,X,tol,100,acsx,avoidOverlaps,boundingBoxes,pbcx,scx);
+	gpY=new GradientProjection(
+            VERTICAL,n,Q,Y,tol,100,acsy,avoidOverlaps,boundingBoxes,pbcy,scy);
+    this->straightenEdges = straightenEdges;
+}
+} // namespace cola
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4
diff --git a/src/libcola/cola.h b/src/libcola/cola.h
new file mode 100644
index 000000000..d4b0d1421
--- /dev/null
+++ b/src/libcola/cola.h
@@ -0,0 +1,242 @@
+#ifndef COLA_H
+#define COLA_H
+
+#include <utility>
+#include <iterator>
+#include <vector>
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <cassert>
+#include "shortest_paths.h"
+#include "gradient_projection.h"
+#include <libvpsc/generate-constraints.h>
+#include "straightener.h"
+
+
+typedef vector<unsigned> Cluster;
+typedef vector<Cluster*> Clusters;
+
+namespace cola {
+    typedef pair<unsigned, unsigned> Edge;
+
+    // defines references to three variables for which the goal function
+    // will be altered to prefer points u-b-v are in a linear arrangement
+    // such that b is placed at u+t(v-u).
+    struct LinearConstraint {
+        LinearConstraint(unsigned u, unsigned v, unsigned b, double w, 
+                double frac_ub, double frac_bv,
+                double* X, double* Y) 
+            : u(u),v(v),b(b),w(w),frac_ub(frac_ub),frac_bv(frac_bv),
+              tAtProjection(true) 
+        {
+            assert(frac_ub<=1.0);
+            assert(frac_bv<=1.0);
+            assert(frac_ub>=0);
+            assert(frac_bv>=0);
+            if(tAtProjection) {
+                double uvx = X[v] - X[u],
+                       uvy = Y[v] - Y[u],
+                       vbx = X[b] - X[u],
+                       vby = Y[b] - Y[u];
+                t = uvx * vbx + uvy * vby;
+                t/= uvx * uvx + uvy * uvy;
+                // p is the projection point of b on line uv
+                //double px = scalarProj * uvx + X[u];
+                //double py = scalarProj * uvy + Y[u];
+                // take t=|up|/|uv|
+            } else {
+                double numerator=X[b]-X[u];
+                double denominator=X[v]-X[u];
+                if(fabs(denominator)<0.001) {
+                    // if line is close to vertical then use Y coords to compute T
+                    numerator=Y[b]-Y[u];
+                    denominator=Y[v]-Y[u];
+                }
+                if(fabs(denominator)<0.0001) {
+                    denominator=1;
+                }
+                t=numerator/denominator;
+            }
+            duu=(1-t)*(1-t);
+            duv=t*(1-t);
+            dub=t-1;
+            dvv=t*t;
+            dvb=-t;
+            dbb=1;
+             //printf("New LC: t=%f\n",t); 
+        }
+        unsigned u;
+        unsigned v;
+        unsigned b;
+        double w; // weight
+        double t;
+        // 2nd partial derivatives of the goal function
+        //   (X[b] - (1-t) X[u] - t X[v])^2
+        double duu;
+        double duv;
+        double dub;
+        double dvv;
+        double dvb;
+        double dbb;
+        // Length of each segment as a fraction of the total edge length
+        double frac_ub;
+        double frac_bv;
+        bool tAtProjection;
+    };
+    typedef vector<LinearConstraint*> LinearConstraints;
+	
+	class TestConvergence {
+    public:
+        double old_stress;
+		TestConvergence(const double& tolerance = 0.001, const unsigned maxiterations = 1000)
+			: old_stress(DBL_MAX),
+              tolerance(tolerance),
+              maxiterations(maxiterations),
+              iterations(0) { }
+        virtual ~TestConvergence() {}
+
+		virtual bool operator()(double new_stress, double* X, double* Y) {
+            //std::cout<<"iteration="<<iterations<<", new_stress="<<new_stress<<std::endl;
+			if (old_stress == DBL_MAX) {
+				old_stress = new_stress;
+                if(++iterations>=maxiterations) {;
+                    return true;
+                } else {
+          			return false;
+                }
+			}
+            bool converged = 
+                fabs(new_stress - old_stress) / (new_stress + 1e-10) < tolerance
+                || ++iterations > maxiterations;
+            old_stress = new_stress;
+			return converged;
+		}
+  	private:
+        double tolerance;
+        unsigned maxiterations;
+        unsigned iterations;
+  	};
+    static TestConvergence defaultTest(0.0001,100);
+	class ConstrainedMajorizationLayout {
+    public:
+		ConstrainedMajorizationLayout(
+                vector<Rectangle*>& rs,
+                vector<Edge>& es,
+				double* eweights,
+                double idealLength,
+				TestConvergence& done=defaultTest)
+			: constrainedLayout(false),
+              n(rs.size()),
+              lapSize(n), lap2(new double*[lapSize]), 
+              Q(lap2), Dij(new double*[lapSize]),
+              tol(0.0001),
+              done(done),
+              X(new double[n]),
+              Y(new double[n]),
+              clusters(NULL),
+              linearConstraints(NULL),
+              gpX(NULL),
+              gpY(NULL),
+              straightenEdges(NULL)
+        {
+            assert(rs.size()==n);
+            boundingBoxes = new Rectangle*[rs.size()];
+            copy(rs.begin(),rs.end(),boundingBoxes);
+
+            double** D=new double*[n];
+            for(unsigned i=0;i<n;i++) {
+                D[i]=new double[n];
+            }
+            shortest_paths::johnsons(n,D,es,eweights);
+            edge_length = idealLength;
+            // Lij_{i!=j}=1/(Dij^2)
+            //
+            for(unsigned i = 0; i<n; i++) {
+                X[i]=rs[i]->getCentreX();
+                Y[i]=rs[i]->getCentreY();
+                double degree = 0;
+                lap2[i]=new double[n];
+                Dij[i]=new double[n];
+                for(unsigned j=0;j<n;j++) {
+                    double w = edge_length * D[i][j];
+                    Dij[i][j]=w;
+                    if(i==j) continue;
+                    degree+=lap2[i][j]=w>1e-30?1.f/(w*w):0;
+                }
+                lap2[i][i]=-degree;
+                delete [] D[i];
+            }
+            delete [] D;
+        }
+
+        void moveBoundingBoxes() {
+            for(unsigned i=0;i<lapSize;i++) {
+                boundingBoxes[i]->moveCentreX(X[i]);
+                boundingBoxes[i]->moveCentreY(Y[i]);
+            }
+        }
+
+        void setupConstraints(
+                AlignmentConstraints* acsx, AlignmentConstraints* acsy,
+                bool avoidOverlaps, 
+                PageBoundaryConstraints* pbcx = NULL,
+                PageBoundaryConstraints* pbcy = NULL,
+                SimpleConstraints* scx = NULL,
+                SimpleConstraints* scy = NULL,
+                Clusters* cs = NULL,
+                vector<straightener::Edge*>* straightenEdges = NULL);
+
+        void addLinearConstraints(LinearConstraints* linearConstraints);
+
+        void setupDummyVars();
+
+        ~ConstrainedMajorizationLayout() {
+            if(boundingBoxes) {
+                delete [] boundingBoxes;
+            }
+            if(constrainedLayout) {
+                delete gpX;
+                delete gpY;
+            }
+            for(unsigned i=0;i<lapSize;i++) {
+                delete [] lap2[i];
+                delete [] Dij[i];
+            }
+            delete [] lap2;
+            delete [] Dij;
+            delete [] X;
+            delete [] Y;
+        }
+		bool run();
+        void straighten(vector<straightener::Edge*>&, Dim);
+        bool avoidOverlaps;
+        bool constrainedLayout;
+    private:
+        double euclidean_distance(unsigned i, unsigned j) {
+            return sqrt(
+                (X[i] - X[j]) * (X[i] - X[j]) +
+                (Y[i] - Y[j]) * (Y[i] - Y[j]));
+        }
+        double compute_stress(double **Dij);
+        void majlayout(double** Dij,GradientProjection* gp, double* coords);
+        void majlayout(double** Dij,GradientProjection* gp, double* coords, 
+                double* b);
+        unsigned n; // is lapSize + dummyVars
+        unsigned lapSize; // lapSize is the number of variables for actual nodes
+        double** lap2; // graph laplacian
+        double** Q; // quadratic terms matrix used in computations
+        double** Dij;
+        double tol;
+		TestConvergence& done;
+        Rectangle** boundingBoxes;
+        double *X, *Y;
+        Clusters* clusters;
+        double edge_length;
+        LinearConstraints *linearConstraints;
+        GradientProjection *gpX, *gpY;
+        vector<straightener::Edge*>* straightenEdges;
+	};
+}
+#endif				// COLA_H
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4
diff --git a/src/libcola/conjugate_gradient.cpp b/src/libcola/conjugate_gradient.cpp
new file mode 100644
index 000000000..5dfb4363d
--- /dev/null
+++ b/src/libcola/conjugate_gradient.cpp
@@ -0,0 +1,113 @@
+#include <math.h>
+#include <stdlib.h>
+#include <valarray>
+#include <cassert>
+#include "conjugate_gradient.h"
+
+/*
+* Authors:
+*   Nathan Hurst <njh@njhurst.com>
+*   Tim Dwyer <tgdwyer@gmail.com>
+*
+* Copyright (C) 2006 Authors
+*
+* Released under GNU LGPL.
+*/
+
+/* lifted wholely from wikipedia.  Well, apart from the bug in the wikipedia version. */
+
+using std::valarray;
+
+static void 
+matrix_times_vector(valarray<double> const &matrix, /* m * n */
+		    valarray<double> const &vec,  /* n */
+		    valarray<double> &result) /* m */
+{
+    unsigned n = vec.size();
+    unsigned m = result.size();
+    assert(m*n == matrix.size());
+    const double* mp = &matrix[0];
+    for (unsigned i = 0; i < m; i++) {
+        double res = 0;
+        for (unsigned j = 0; j < n; j++)
+            res += *mp++ * vec[j];
+        result[i] = res;
+    }
+}
+
+static double Linfty(valarray<double> const &vec) {
+    return std::max(vec.max(), -vec.min());
+}
+
+double
+inner(valarray<double> const &x, 
+      valarray<double> const &y) {
+    double total = 0;
+    for(unsigned i = 0; i < x.size(); i++)
+        total += x[i]*y[i];
+    return total;// (x*y).sum(); <- this is more concise, but ineff
+}
+
+void 
+conjugate_gradient(double **A, 
+                   double *x, 
+                   double *b, 
+                   unsigned n, 
+                   double tol,
+                   unsigned max_iterations) {
+    valarray<double> vA(n*n);
+    valarray<double> vx(n);
+    valarray<double> vb(n);
+    for(unsigned i=0;i<n;i++) {
+        vx[i]=x[i];
+        vb[i]=b[i];
+        for(unsigned j=0;j<n;j++) {
+            vA[i*n+j]=A[i][j];
+        }
+    }
+    conjugate_gradient(vA,vx,vb,n,tol,max_iterations);
+    for(unsigned i=0;i<n;i++) {
+        x[i]=vx[i];
+    }
+}
+void 
+conjugate_gradient(valarray<double> const &A, 
+		   valarray<double> &x, 
+		   valarray<double> const &b, 
+		   unsigned n, double tol,
+		   unsigned max_iterations) {
+    valarray<double> Ap(n), p(n), r(n);
+    matrix_times_vector(A,x,Ap);
+    r=b-Ap; 
+    double r_r = inner(r,r);
+    unsigned k = 0;
+    tol *= tol;
+    while(k < max_iterations && r_r > tol) {
+        k++;
+        double r_r_new = r_r;
+        if(k == 1)
+            p = r;
+        else {
+            r_r_new = inner(r,r);
+            p = r + (r_r_new/r_r)*p;
+        }
+        matrix_times_vector(A, p, Ap);
+        double alpha_k = r_r_new / inner(p, Ap);
+        x += alpha_k*p;
+        r -= alpha_k*Ap;
+        r_r = r_r_new;
+    }
+    printf("njh: %d iters, Linfty = %g L2 = %g\n", k, 
+    std::max(-r.min(), r.max()), sqrt(r_r));
+    // x is solution
+}
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4
diff --git a/src/libcola/conjugate_gradient.h b/src/libcola/conjugate_gradient.h
new file mode 100644
index 000000000..17e59c9af
--- /dev/null
+++ b/src/libcola/conjugate_gradient.h
@@ -0,0 +1,23 @@
+#ifndef _CONJUGATE_GRADIENT_H
+#define _CONJUGATE_GRADIENT_H
+
+#include <valarray>
+
+double
+inner(std::valarray<double> const &x, 
+      std::valarray<double> const &y);
+
+void 
+conjugate_gradient(double **A, 
+		   double *x, 
+		   double *b, 
+		   unsigned n, 
+		   double tol,
+           unsigned max_iterations);
+void 
+conjugate_gradient(std::valarray<double> const &A, 
+		   std::valarray<double> &x, 
+		   std::valarray<double> const &b, 
+		   unsigned n, double tol,
+		   unsigned max_iterations);
+#endif // _CONJUGATE_GRADIENT_H
diff --git a/src/libcola/cycle_detector.cpp b/src/libcola/cycle_detector.cpp
new file mode 100644
index 000000000..ddc056e4d
--- /dev/null
+++ b/src/libcola/cycle_detector.cpp
@@ -0,0 +1,228 @@
+/* Cycle detector that returns a list of 
+ * edges involved in cycles in a digraph.
+ *
+ * Kieran Simpson 2006
+*/
+#include <iostream>
+#include <stack>
+#include <vector>
+#include <cassert>
+#include <cycle_detector.h>
+
+#define RUN_DEBUG
+
+using namespace std;
+using namespace cola;
+
+// a global var representing time
+TimeStamp Time;
+
+CycleDetector::CycleDetector(unsigned numVertices, Edges *edges)  {
+  this->V = numVertices; 
+  this->edges = edges;
+
+  // make the adjacency matrix
+  this->make_matrix();
+  assert(nodes.size() == this->V);
+}
+
+CycleDetector::~CycleDetector()  {
+  if (!nodes.empty())  { for (unsigned i = 0; i < nodes.size(); i++)  { delete nodes[i]; } }
+}
+
+void CycleDetector::make_matrix()  {
+  Edges::iterator ei;
+  Edge anEdge;
+  Node *newNode;
+
+  if (!nodes.empty())  { for (map<unsigned, Node *>::iterator ni = nodes.begin(); ni != nodes.end(); ni++)  { delete nodes[ni->first]; } }
+  nodes.clear();
+  traverse.clear();
+
+  // we should have no nodes in the list
+  assert(nodes.empty());
+  assert(traverse.empty());
+
+  // from the edges passed, fill the adjacency matrix
+  for (ei = edges->begin(); ei != edges->end(); ei++)  {
+    anEdge = *ei;
+    // the matrix is indexed by the first vertex of the edge
+    // the second vertex of the edge is pushed onto another
+    // vector of all vertices connected to the first vertex
+    // with a directed edge
+    #ifdef ADJMAKE_DEBUG
+      cout << "vertex1: " << anEdge.first << ", vertex2: " << anEdge.second << endl;
+    #endif
+    if (nodes.find(anEdge.first) == nodes.end())  {
+      #ifdef ADJMAKE_DEBUG
+        cout << "Making a new vector indexed at: " << anEdge.first << endl;
+      #endif
+
+      newNode = new Node(anEdge.first);
+      newNode->dests.push_back(anEdge.second);
+      nodes[anEdge.first] = newNode;
+    }
+    else  {
+      nodes[anEdge.first]->dests.push_back(anEdge.second);
+    }
+
+    // check if the destination vertex exists in the nodes map
+    if (nodes.find(anEdge.second) == nodes.end())  {
+      #ifdef ADJMAKE_DEBUG
+        cerr << "Making a new vector indexed at: " << anEdge.second << endl;
+      #endif
+
+      newNode = new Node(anEdge.second);
+      nodes[anEdge.second] = newNode;
+    }
+  }
+
+  assert(!nodes.empty());
+
+  // the following block is code to print out
+  // the adjacency matrix.
+  #ifdef ADJMAKE_DEBUG
+    for (map<unsigned, Node *>::iterator ni = nodes.begin(); ni != nodes.end(); ni++)  {
+      Node *node = ni->second;
+      cout << "nodes[" << node->id << "]: ";
+      
+      if (isSink(node))  { cout << "SINK"; }
+      else  {
+        for (unsigned j = 0; j < node->dests.size(); j++)  { cout << node->dests[j] << " "; }
+      }
+      cout << endl;
+    }
+  #endif
+}
+
+vector<bool> *CycleDetector::detect_cycles()  {
+  vector<bool> *cyclicEdgesMapping = NULL;
+
+  assert(!nodes.empty());
+  assert(!edges->empty());
+
+  // make a copy of the graph to ensure that we have visited all
+  // vertices
+  traverse.clear(); assert(traverse.empty());
+  for (unsigned i = 0; i < V; i++)  { traverse[i] = false; }
+  #ifdef SETUP_DEBUG
+    for (map<unsigned, bool>::iterator ivi = traverse.begin(); ivi != traverse.end(); ivi++)  {
+      cout << "traverse{" << ivi->first << "}: " << ivi->second << endl;
+    }
+  #endif
+
+  // set up the mapping between the edges and their cyclic truth
+  for(unsigned i = 0; i < edges->size(); i++)  { cyclicEdges[(*edges)[i]] = false; }
+
+  // find the cycles
+  assert(nodes.size() > 1);
+
+  // while we still have vertices to visit, visit.
+  while (!traverse.empty())  {
+    // mark any vertices seen in a previous run as closed
+    while (!seenInRun.empty())  {
+      unsigned v = seenInRun.top();
+      seenInRun.pop();
+      #ifdef RUN_DEBUG
+        cout << "Marking vertex(" << v << ") as CLOSED" << endl;
+      #endif
+      nodes[v]->status = Node::DoneVisiting;
+    }
+
+    assert(seenInRun.empty());
+
+    #ifdef VISIT_DEBUG
+      cout << "begining search at vertex(" << traverse.begin()->first << ")" << endl;
+    #endif
+
+    Time = 0;
+
+    // go go go
+    visit(traverse.begin()->first);
+  }
+
+  // clean up
+  while (!seenInRun.empty())  { seenInRun.pop(); }
+  assert(seenInRun.empty());
+  assert(traverse.empty());
+
+  cyclicEdgesMapping = new vector<bool>(edges->size(), false);
+
+  for (unsigned i = 0; i < edges->size(); i++)  {
+    if (cyclicEdges[(*edges)[i]])  { 
+      (*cyclicEdgesMapping)[i] = true;
+      #ifdef OUTPUT_DEBUG
+        cout << "Setting cyclicEdgesMapping[" << i << "] to true" << endl;
+      #endif
+    }
+  }
+
+  return cyclicEdgesMapping;
+}
+
+void CycleDetector::mod_graph(unsigned numVertices, Edges *edges)  {
+  this->V = numVertices;
+  this->edges = edges;
+  // remake the adjaceny matrix
+  this->make_matrix();
+  assert(nodes.size() == this->V);
+}
+
+void CycleDetector::visit(unsigned k)  {
+  unsigned cycleOpen;
+
+  // state that we have seen this vertex
+  if (traverse.find(k) != traverse.end())  {
+    #ifdef VISIT_DEBUG
+      cout << "Visiting vertex(" << k << ") for the first time" << endl;
+    #endif
+    traverse.erase(k);
+  }
+
+  seenInRun.push(k);
+
+  // set up this node as being visited.
+  nodes[k]->stamp = ++Time;
+  nodes[k]->status = Node::BeingVisited;
+
+  // traverse to all the vertices adjacent to this vertex.
+  for (unsigned n = 0; n < nodes[k]->dests.size(); n++)  {
+    unsigned v = nodes[k]->dests[n];
+
+    if (nodes[v]->status != Node::DoneVisiting)  {
+      if (nodes[v]->status == Node::NotVisited)  {  
+        // visit this node
+	#ifdef VISIT_DEBUG
+          cout << "traversing from vertex(" << k << ") to vertex(" << v << ")" << endl;
+	#endif
+        visit(v);
+      }
+
+      // if we are part of a cycle get the timestamp of the ancestor
+      if (nodes[v]->cyclicAncestor != NULL)  { cycleOpen = nodes[v]->cyclicAncestor->stamp; }
+      // else just get the timestamp of the node we just visited
+      else  { cycleOpen = nodes[v]->stamp; }
+
+      // compare the stamp of the traversal with our stamp
+      if (cycleOpen <= nodes[k]->stamp)  {
+	if (nodes[v]->cyclicAncestor == NULL)  { nodes[v]->cyclicAncestor = nodes[v]; }
+	// store the cycle
+	cyclicEdges[Edge(k,v)] = true;
+        // this node is part of a cycle
+        if (nodes[k]->cyclicAncestor == NULL)  { nodes[k]->cyclicAncestor = nodes[v]->cyclicAncestor; }
+
+	// see if we are part of a cycle with a cyclicAncestor that possesses a lower timestamp
+	if (nodes[v]->cyclicAncestor->stamp < nodes[k]->cyclicAncestor->stamp)  { nodes[k]->cyclicAncestor = nodes[v]->cyclicAncestor; }
+      }
+    }
+  }
+}
+
+
+// determines whether or not a vertex is a sink
+bool CycleDetector::isSink(Node *node)  {
+  // a vertex is a sink if it has no outgoing edges,
+  // or that the adj entry is empty
+  if (node->dests.empty())  { return true; }
+  else  { return false; }
+}
diff --git a/src/libcola/cycle_detector.h b/src/libcola/cycle_detector.h
new file mode 100644
index 000000000..5cd52e324
--- /dev/null
+++ b/src/libcola/cycle_detector.h
@@ -0,0 +1,54 @@
+#ifndef CYCLE_DETECTOR_H
+#define CYCLE_DETECTOR_H
+
+#include <map>
+#include <vector>
+#include <stack>
+#include "cola.h"
+
+typedef unsigned TimeStamp;
+typedef std::vector<cola::Edge> Edges;
+typedef std::vector<bool> CyclicEdges;
+class Node;
+
+class CycleDetector  {
+  public:
+    CycleDetector(unsigned numVertices, Edges *edges);
+    ~CycleDetector();
+    std::vector<bool> *detect_cycles();
+    void mod_graph(unsigned numVertices, Edges *edges);
+    unsigned getV()  { return this->V; }
+    Edges *getEdges()  { return this->edges; }
+
+  private:
+    // attributes
+    unsigned V;
+    Edges *edges;
+
+    // internally used variables.
+    std::map<unsigned, Node *> nodes; // the nodes in the graph
+    std::map<cola::Edge, bool> cyclicEdges; // the cyclic edges in the graph.
+    std::map<unsigned, bool> traverse; // nodes still left to visit in the graph
+    std::stack<unsigned> seenInRun; // nodes visited in a single pass.
+
+    // internally used methods
+    void make_matrix();
+    void visit(unsigned k);
+    bool isSink(Node *node);
+};
+
+class Node  {
+  public:
+    enum StatusType { NotVisited, BeingVisited, DoneVisiting };
+
+    unsigned id;
+    TimeStamp stamp;
+    Node *cyclicAncestor;
+    vector<unsigned> dests;
+    StatusType status;
+
+    Node(unsigned id)  { this->id = id; cyclicAncestor = NULL; status = NotVisited; }
+    ~Node()  {}
+};
+
+#endif
diff --git a/src/libcola/defs.h b/src/libcola/defs.h
new file mode 100644
index 000000000..cd8084c2c
--- /dev/null
+++ b/src/libcola/defs.h
@@ -0,0 +1,132 @@
+/* $Id: defs.h,v 1.5 2005/10/18 18:42:59 ellson Exp $ $Revision: 1.5 $ */
+/* vim:set shiftwidth=4 ts=8: */
+
+/**********************************************************
+*      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 _DEFS_H_
+#define _DEFS_H_
+
+#include "neato.h"
+
+#ifdef __cplusplus
+    enum Style { regular, invisible };
+    struct vtx_data {
+	int nedges;
+	int *edges;
+	float *ewgts;
+	Style *styles;
+	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) {
+	if (x >= y)
+	    return x;
+	else
+	    return y;
+    } inline double min(double x, double y) {
+	if (x <= y)
+	    return x;
+	else
+	    return y;
+    }
+
+    inline int max(int x, int y) {
+	if (x >= y)
+	    return x;
+	else
+	    return y;
+    }
+
+    inline int min(int x, int y) {
+	if (x <= y)
+	    return x;
+	else
+	    return y;
+    }
+
+    struct Point {
+	double x;
+	double y;
+	int operator==(Point other) {
+	    return x == other.x && y == other.y;
+    }};
+#else
+#undef inline
+#define inline
+#define NOTUSED(var)      (void) var
+
+#include <macros.h>
+    extern void *gmalloc(size_t);
+#define DIGCOLA 1
+
+#ifdef USE_STYLES
+    typedef enum { regular, invisible } Style;
+#endif
+    typedef struct {
+	int nedges;		/* no. of neighbors, including self */
+	int *edges;		/* edges[0..(nedges-1)] are neighbors; edges[0] is self */
+	float *ewgts;		/* preferred edge lengths */
+	float *eweights;	/* edge weights */
+	node_t *np;		/* original node */
+#ifdef USE_STYLES
+	Style *styles;
+#endif
+#ifdef DIGCOLA
+	float *edists; /* directed dist reflecting the direction of the edge */
+#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
+    typedef struct {
+	double x;
+	double y;
+    } Point;
+#endif
+
+#endif
+
+#endif
+
+#ifdef __cplusplus
+}
+#endif
diff --git a/src/libcola/gradient_projection.cpp b/src/libcola/gradient_projection.cpp
new file mode 100644
index 000000000..061ba0f1a
--- /dev/null
+++ b/src/libcola/gradient_projection.cpp
@@ -0,0 +1,234 @@
+/**********************************************************
+ *
+ * Solve a quadratic function f(X) = X' A X + b X
+ * subject to a set of separation constraints cs
+ *
+ * Tim Dwyer, 2006
+ **********************************************************/
+
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <stdio.h>
+#include <float.h>
+#include <cassert>
+#include <libvpsc/solve_VPSC.h>
+#include <libvpsc/variable.h>
+#include <libvpsc/constraint.h>
+#include "gradient_projection.h"
+#include <iostream>
+
+using namespace std;
+//#define CONMAJ_LOGGING 1
+
+static void dumpVPSCException(char const *str, IncVPSC* vpsc) {
+    cerr<<str<<endl;
+    unsigned m;
+    Constraint** cs = vpsc->getConstraints(m);
+    for(unsigned i=0;i<m;i++) {
+        cerr << *cs[i] << endl;
+    }
+}
+/*
+ * Use gradient-projection to solve an instance of
+ * the Variable Placement with Separation Constraints problem.
+ * Uses sparse matrix techniques to handle pairs of dummy
+ * vars.
+ */
+unsigned GradientProjection::solve(double * b) {
+	unsigned i,j,counter;
+	if(max_iterations==0) return 0;
+
+	bool converged=false;
+
+    IncVPSC* vpsc=NULL;
+
+    vpsc = setupVPSC();
+    //cerr << "in gradient projection: n=" << n << endl;
+    for (i=0;i<n;i++) {
+        assert(!isnan(place[i]));
+        assert(!isinf(place[i]));
+        vars[i]->desiredPosition=place[i];
+    }
+    try {
+        vpsc->satisfy();
+    } catch (char const *str) {
+        dumpVPSCException(str,vpsc);
+	}
+
+    for (i=0;i<n;i++) {
+        place[i]=vars[i]->position();
+    }
+    for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
+        (*it)->updatePosition();
+    }
+    	
+	for (counter=0; counter<max_iterations&&!converged; counter++) {
+		converged=true;		
+		// find steepest descent direction
+        //  g = 2 ( b - Ax )
+		for (i=0; i<n; i++) {
+			old_place[i]=place[i];
+			g[i] = b[i];
+			for (j=0; j<n; j++) {
+				g[i] -= A[i][j]*place[j];
+			}
+            g[i] *= 2.0;
+		}		
+        for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
+            (*it)->computeDescentVector();
+        }
+        // compute step size: alpha = ( g' g ) / ( 2 g' A g )
+        //   g terms for dummy vars cancel out so don't consider
+		double 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 += A[i][j]*g[j];
+			}
+			denominator -= 2.0 * r*g[i];
+		}
+		double alpha = numerator/denominator;
+
+        // move to new unconstrained position
+		for (i=0; i<n; i++) {
+			place[i]-=alpha*g[i];
+            assert(!isnan(place[i]));
+            assert(!isinf(place[i]));
+            vars[i]->desiredPosition=place[i];
+		}
+        for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
+            (*it)->steepestDescent(alpha);
+        }
+
+        //project to constraint boundary
+        try {
+            vpsc->satisfy();
+        } catch (char const *str) {
+            dumpVPSCException(str,vpsc);
+        }
+        for (i=0;i<n;i++) {
+            place[i]=vars[i]->position();
+        }
+        for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
+            (*it)->updatePosition();
+        }
+        // compute d, the vector from last pnt to projection pnt
+		for (i=0; i<n; i++) {
+			d[i]=place[i]-old_place[i];
+		}	
+		// now compute beta, optimal step size from last pnt to projection pnt
+        //   beta = ( g' d ) / ( 2 d' A d )
+		numerator = 0, denominator = 0;
+		for (i=0; i<n; i++) {
+			numerator += g[i] * d[i];
+			r=0;
+			for (j=0; j<n; j++) {
+				r += A[i][j] * d[j];
+			}
+			denominator += 2.0 * r * d[i];
+		}
+        for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
+            (*it)->betaCalc(numerator,denominator);
+        }
+		double beta = numerator/denominator;
+
+        // 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) {
+            for (i=0; i<n; i++) {
+                place[i]=old_place[i]+beta*d[i];
+            }
+            for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
+                (*it)->feasibleDescent(beta);
+            }
+        }
+		double test=0;
+		for (i=0; i<n; i++) {
+			test += fabs(place[i]-old_place[i]);
+		}
+        for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
+            test += (*it)->absoluteDisplacement();
+        }
+		if(test>tolerance) {
+			converged=false;
+		}
+	}
+    destroyVPSC(vpsc);
+	return counter;
+}
+// Setup an instance of the Variable Placement with Separation Constraints
+// for one iteration.
+// Generate transient local constraints --- such as non-overlap constraints 
+// --- that are only relevant to one iteration, and merge these with the
+// global constraint list (including alignment constraints,
+// dir-edge constraints, containment constraints, etc).
+IncVPSC* GradientProjection::setupVPSC() {
+    Constraint **cs;
+    //assert(lcs.size()==0);
+    
+    for(DummyVars::iterator dit=dummy_vars.begin();
+            dit!=dummy_vars.end(); ++dit) {
+        (*dit)->setupVPSC(vars,lcs);
+    }
+    Variable** vs = new Variable*[vars.size()];
+    for(unsigned i=0;i<vars.size();i++) {
+        vs[i]=vars[i];
+    }
+    if(nonOverlapConstraints) {
+        Constraint** tmp_cs=NULL;
+        unsigned m=0;
+        if(k==HORIZONTAL) {
+            Rectangle::setXBorder(0.0001);
+            m=generateXConstraints(n,rs,vs,tmp_cs,true); 
+            Rectangle::setXBorder(0);
+        } else {
+            m=generateYConstraints(n,rs,vs,tmp_cs); 
+        }
+        for(unsigned i=0;i<m;i++) {
+            lcs.push_back(tmp_cs[i]);
+        }
+    }
+    cs = new Constraint*[lcs.size() + gcs.size()];
+    unsigned m = 0 ;
+    for(Constraints::iterator ci = lcs.begin();ci!=lcs.end();++ci) {
+        cs[m++] = *ci;
+    }
+    for(Constraints::iterator ci = gcs.begin();ci!=gcs.end();++ci) {
+        cs[m++] = *ci;
+    }
+    return new IncVPSC(vars.size(),vs,m,cs);
+}
+void GradientProjection::clearDummyVars() {
+    for(DummyVars::iterator i=dummy_vars.begin();i!=dummy_vars.end();++i) {
+        delete *i;
+    }
+    dummy_vars.clear();
+}
+void GradientProjection::destroyVPSC(IncVPSC *vpsc) {
+    if(acs) {
+        for(AlignmentConstraints::iterator ac=acs->begin(); ac!=acs->end();++ac) {
+            (*ac)->updatePosition();
+        }
+    }
+    unsigned m,n;
+    Constraint** cs = vpsc->getConstraints(m);
+    const Variable* const* vs = vpsc->getVariables(n);
+    delete vpsc;
+    delete [] cs;
+    delete [] vs;
+    for(Constraints::iterator i=lcs.begin();i!=lcs.end();i++) {
+            delete *i;
+    }
+    lcs.clear();
+    //cout << " Vars count = " << vars.size() << " Dummy vars cnt=" << dummy_vars.size() << endl;
+    vars.resize(vars.size()-dummy_vars.size()*2);
+    for(DummyVars::iterator i=dummy_vars.begin();i!=dummy_vars.end();++i) {
+        DummyVarPair* p = *i;
+        delete p->left;
+        delete p->right;
+    }
+}
+
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
diff --git a/src/libcola/gradient_projection.h b/src/libcola/gradient_projection.h
new file mode 100644
index 000000000..e8b72180b
--- /dev/null
+++ b/src/libcola/gradient_projection.h
@@ -0,0 +1,266 @@
+#ifndef _GRADIENT_PROJECTION_H
+#define _GRADIENT_PROJECTION_H
+
+#include <libvpsc/solve_VPSC.h>
+#include <libvpsc/variable.h>
+#include <libvpsc/constraint.h>
+#include <libvpsc/generate-constraints.h>
+#include <vector>
+#include <iostream>
+#include <math.h>
+
+using namespace std;
+
+typedef vector<Constraint*> Constraints;
+typedef vector<Variable*> Variables;
+typedef vector<pair<unsigned,double> > OffsetList;
+
+class SimpleConstraint {
+public:
+    SimpleConstraint(unsigned l, unsigned r, double g) 
+        : left(l), right(r), gap(g)  {}
+    unsigned left;
+    unsigned right;
+    double gap;
+};
+typedef vector<SimpleConstraint*> SimpleConstraints;
+class AlignmentConstraint {
+friend class GradientProjection;
+public:
+    AlignmentConstraint(double pos) : position(pos), variable(NULL) {}
+    void updatePosition() {
+        position = variable->position();
+    }
+    OffsetList offsets;
+    void* guide;
+    double position;
+private:
+    Variable* variable;
+};
+typedef vector<AlignmentConstraint*> AlignmentConstraints;
+
+class PageBoundaryConstraints {
+public:
+    PageBoundaryConstraints(double lm, double rm, double w)
+        : leftMargin(lm), rightMargin(rm), weight(w) { }
+    void createVarsAndConstraints(Variables &vs, Constraints &cs) {
+        Variable* vl, * vr;
+        // create 2 dummy vars, based on the dimension we are in
+        vs.push_back(vl=new Variable(vs.size(), leftMargin, weight));
+        vs.push_back(vr=new Variable(vs.size(), rightMargin, weight));
+
+        // for each of the "real" variables, create a constraint that puts that var
+        // between our two new dummy vars, depending on the dimension.
+        for(OffsetList::iterator o=offsets.begin(); o!=offsets.end(); ++o)  {
+            cs.push_back(new Constraint(vl, vs[o->first], o->second));
+            cs.push_back(new Constraint(vs[o->first], vr, o->second));
+        }
+    }
+    OffsetList offsets;
+private:
+    double leftMargin;
+    double rightMargin;
+    double weight;
+};
+
+typedef vector<pair<unsigned,double> > CList;
+/**
+ * A DummyVarPair is a pair of variables with an ideal distance between them and which have no
+ * other interaction with other variables apart from through constraints.  This means that
+ * the entries in the laplacian matrix for dummy vars and other vars would be 0 - thus, sparse
+ * matrix techniques can be used in laplacian operations.
+ * The constraints are specified by a two lists of pairs of variable indexes and required separation.
+ * The two lists are:
+ *   leftof: variables to which left must be to the left of, 
+ *   rightof: variables to which right must be to the right of. 
+ */
+class DummyVarPair {
+public:
+    DummyVarPair(double desiredDist) : dist(desiredDist), lap2(1.0/(desiredDist*desiredDist)) { }
+    CList leftof; // variables to which left dummy var must be to the left of
+    CList rightof; // variables to which right dummy var must be to the right of
+    double place_l;
+    double place_r;
+    void computeLinearTerm(double euclideanDistance) {   
+        if(euclideanDistance > 1e-30) {
+            b = place_r - place_l;
+            b /= euclideanDistance * dist;
+        } else { b=0; }
+    }
+    double stress(double euclideanDistance) {
+        double diff = dist - euclideanDistance;
+        return diff*diff / (dist*dist);
+    }
+private:
+friend class GradientProjection; 
+    /**
+     * Setup vars and constraints for an instance of the VPSC problem.
+     * Adds generated vars and constraints to the argument vectors.
+     */
+    void setupVPSC(Variables &vars, Constraints &cs) {
+        double weight=1;
+        left = new Variable(vars.size(),place_l,weight);
+        vars.push_back(left);
+        right = new Variable(vars.size(),place_r,weight);
+        vars.push_back(right);
+        for(CList::iterator cit=leftof.begin();
+                cit!=leftof.end(); ++cit) {
+            Variable* v = vars[(*cit).first];
+            cs.push_back(new Constraint(left,v,(*cit).second)); 
+        }
+        for(CList::iterator cit=rightof.begin();
+                cit!=rightof.end(); ++cit) {
+            Variable* v = vars[(*cit).first];
+            cs.push_back(new Constraint(v,right,(*cit).second)); 
+        }
+    }
+    /**
+     * Extract the result of a VPSC solution to the variable positions
+     */
+    void updatePosition() {
+        place_l=left->position();
+        place_r=right->position();
+    }
+    /**
+     * Compute the descent vector, also copying the current position to old_place for
+     * future reference
+     */
+    void computeDescentVector() {
+        old_place_l=place_l;
+        old_place_r=place_r;
+        g = 2.0 * ( b + lap2 * ( place_l - place_r ) );
+    }
+    /** 
+     * move in the direction of steepest descent (based on g computed by computeGradient)
+     * alpha is the step size.
+     */
+    void steepestDescent(double alpha) {
+        place_l -= alpha*g;
+        place_r += alpha*g;
+        left->desiredPosition=place_l;
+        right->desiredPosition=place_r;
+    }
+    /**
+     * add dummy vars' contribution to numerator and denominator for 
+     * beta step size calculation
+     */
+    void betaCalc(double &numerator, double &denominator) {
+        double dl = place_l-old_place_l,
+               dr = place_r-old_place_r,
+               r = 2.0 * lap2 * ( dr - dl );
+        numerator += g * ( dl - dr );
+        denominator += r*dl - r * dr;
+    }
+    /**
+     * move by stepsize beta from old_place to place
+     */
+    void feasibleDescent(double beta) {
+        left->desiredPosition = place_l = old_place_l + beta*(place_l - old_place_l);
+        right->desiredPosition = place_r = old_place_r + beta*(place_r - old_place_r);
+    }
+    double absoluteDisplacement() {
+        return fabs(place_l - old_place_l) + fabs(place_r - old_place_r);
+    }
+    double dist; // ideal distance between vars
+    double b; // linear coefficient in quad form for left (b_right = -b)
+    Variable* left; // Variables used in constraints
+    Variable* right;
+    double lap2; // laplacian entry
+    double g; // descent vec for quad form for left (g_right = -g)
+    double old_place_l; // old_place is where the descent vec g was computed
+    double old_place_r;
+};
+typedef vector<DummyVarPair*> DummyVars;
+
+enum Dim { HORIZONTAL, VERTICAL };
+
+class GradientProjection {
+public:
+	GradientProjection(
+        const Dim k,
+		unsigned n, 
+		double** A,
+		double* x,
+		double tol,
+		unsigned max_iterations,
+        AlignmentConstraints* acs=NULL,
+        bool nonOverlapConstraints=false,
+        Rectangle** rs=NULL,
+        PageBoundaryConstraints *pbc = NULL,
+        SimpleConstraints *sc = NULL)
+            : k(k), n(n), A(A), place(x), rs(rs),
+              nonOverlapConstraints(nonOverlapConstraints),
+              tolerance(tol), acs(acs), max_iterations(max_iterations),
+              g(new double[n]), d(new double[n]), old_place(new double[n]),
+              constrained(false)
+    {
+        for(unsigned i=0;i<n;i++) {
+            vars.push_back(new Variable(i,1,1));
+        }
+        if(acs) {
+            for(AlignmentConstraints::iterator iac=acs->begin();
+                    iac!=acs->end();++iac) {
+                AlignmentConstraint* ac=*iac;
+                Variable *v=ac->variable=new Variable(vars.size(),ac->position,0.0001);
+                vars.push_back(v);
+                for(OffsetList::iterator o=ac->offsets.begin();
+                        o!=ac->offsets.end();
+                        o++) {
+                    gcs.push_back(new Constraint(v,vars[o->first],o->second,true));
+                }
+            }
+        }
+        if (pbc)  {          
+            pbc->createVarsAndConstraints(vars,gcs);
+        }
+        if (sc) {
+            for(SimpleConstraints::iterator c=sc->begin(); c!=sc->end();++c) {
+                gcs.push_back(new Constraint(
+                        vars[(*c)->left],vars[(*c)->right],(*c)->gap));
+            }
+        }
+        if(!gcs.empty() || nonOverlapConstraints) {
+            constrained=true;
+        }
+	}
+    ~GradientProjection() {
+        delete [] g;
+        delete [] d;
+        delete [] old_place;
+        for(Constraints::iterator i(gcs.begin()); i!=gcs.end(); i++) {
+            delete *i;
+        }
+        gcs.clear();
+        for(unsigned i=0;i<vars.size();i++) {
+            delete vars[i];
+        }
+    }
+    void clearDummyVars();
+	unsigned solve(double* b);
+    DummyVars dummy_vars; // special vars that must be considered in Lapl.
+private:
+    IncVPSC* setupVPSC();
+    void destroyVPSC(IncVPSC *vpsc);
+    Dim k;
+	unsigned n; // number of actual vars
+	double** A; // Graph laplacian matrix
+    double* place;
+	Variables vars; // all variables
+                          // computations
+    Constraints gcs; /* global constraints - persist throughout all
+                                iterations */
+    Constraints lcs; /* local constraints - only for current iteration */
+    Rectangle** rs;
+    bool nonOverlapConstraints;
+    double tolerance;
+    AlignmentConstraints* acs;
+    unsigned max_iterations;
+	double* g; /* gradient */
+	double* d;
+	double* old_place;
+    bool constrained;
+};
+
+#endif /* _GRADIENT_PROJECTION_H */
+
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
diff --git a/src/libcola/shortest_paths.cpp b/src/libcola/shortest_paths.cpp
new file mode 100644
index 000000000..4f4183b07
--- /dev/null
+++ b/src/libcola/shortest_paths.cpp
@@ -0,0 +1,100 @@
+// vim: set cindent 
+// vim: ts=4 sw=4 et tw=0 wm=0
+#include "shortest_paths.h"
+#include <float.h>
+#include <cassert>
+#include <iostream>
+#include <libvpsc/pairingheap/PairingHeap.h>
+using namespace std;
+namespace shortest_paths {
+// O(n^3) time.  Slow, but fool proof.  Use for testing.
+void floyd_warshall(
+        unsigned n,
+        double** D, 
+        vector<Edge>& es,
+        double* eweights) 
+{
+    for(unsigned i=0;i<n;i++) {
+        for(unsigned j=0;j<n;j++) {
+            if(i==j) D[i][j]=0;
+            else D[i][j]=DBL_MAX;
+        }
+    }
+    for(unsigned i=0;i<es.size();i++) {
+        unsigned u=es[i].first, v=es[i].second;
+        assert(u<n&&v<n);
+        D[u][v]=D[v][u]=eweights[i];
+    }
+    for(unsigned k=0; k<n; k++) {
+        for(unsigned i=0; i<n; i++) {
+            for(unsigned j=0; j<n; j++) {
+                D[i][j]=min(D[i][j],D[i][k]+D[k][j]);
+            }
+        }
+    }
+}
+void dijkstra_init(Node* vs, vector<Edge>& es, double* eweights) {
+    for(unsigned i=0;i<es.size();i++) {
+        unsigned u=es[i].first, v=es[i].second;
+        vs[u].neighbours.push_back(&vs[v]);
+        vs[u].nweights.push_back(eweights[i]);
+        vs[v].neighbours.push_back(&vs[u]);
+        vs[v].nweights.push_back(eweights[i]);
+    }
+}
+void dijkstra(
+        unsigned s,
+        unsigned n,
+        Node* vs,
+        double* d)
+{
+    assert(s<n);
+    for(unsigned i=0;i<n;i++) {
+        vs[i].id=i;
+        vs[i].d=DBL_MAX;
+        vs[i].p=NULL;
+    }
+    vs[s].d=0;
+    PairingHeap<Node*> Q(&compareNodes);
+    for(unsigned i=0;i<n;i++) {
+        vs[i].qnode = Q.insert(&vs[i]);
+    }
+    while(!Q.isEmpty()) {
+        Node *u=Q.extractMin();
+        d[u->id]=u->d;
+        for(unsigned i=0;i<u->neighbours.size();i++) {
+            Node *v=u->neighbours[i];
+            double w=u->nweights[i];
+            if(v->d>u->d+w) {
+                v->p=u;
+                v->d=u->d+w;
+                Q.decreaseKey(v->qnode,v);
+            }
+        }
+    }
+}
+void dijkstra(
+        unsigned s,
+        unsigned n,
+        double* d,
+        vector<Edge>& es,
+        double* eweights)
+{
+    assert(s<n);
+    Node vs[n];
+    dijkstra_init(vs,es,eweights);
+    dijkstra(s,n,vs,d);
+}
+void johnsons(
+        unsigned n,
+        double** D, 
+        vector<Edge>& es,
+        double* eweights) 
+{
+    Node vs[n];
+    dijkstra_init(vs,es,eweights);
+    for(unsigned k=0;k<n;k++) {
+        dijkstra(k,n,vs,D[k]);
+    }
+}
+}
diff --git a/src/libcola/shortest_paths.h b/src/libcola/shortest_paths.h
new file mode 100644
index 000000000..20107caf0
--- /dev/null
+++ b/src/libcola/shortest_paths.h
@@ -0,0 +1,28 @@
+// vim: set cindent 
+// vim: ts=4 sw=4 et tw=0 wm=0
+#include <vector>
+using namespace std;
+template <class T>
+class PairNode;
+namespace shortest_paths {
+
+struct Node {
+    unsigned id;
+    double d;
+    Node* p; // predecessor    
+    vector<Node*> neighbours;
+    vector<double> nweights;
+    PairNode<Node*>* qnode;
+};
+inline bool compareNodes(Node *const &u, Node *const &v) {
+	return u->d < v->d;
+}
+
+typedef pair<unsigned,unsigned> Edge;
+void floyd_warshall(unsigned n, double** D,
+        vector<Edge>& es,double* eweights); 
+void johnsons(unsigned n, double** D,
+        vector<Edge>& es, double* eweights);
+void dijkstra(unsigned s, unsigned n, double* d, 
+        vector<Edge>& es, double* eweights);
+}
diff --git a/src/libcola/straightener.cpp b/src/libcola/straightener.cpp
new file mode 100644
index 000000000..6b062eb32
--- /dev/null
+++ b/src/libcola/straightener.cpp
@@ -0,0 +1,360 @@
+/*
+** vim: set cindent 
+** vim: ts=4 sw=4 et tw=0 wm=0
+*/
+/**
+ * \brief Functions to automatically generate constraints for the
+ * rectangular node overlap removal problem.
+ *
+ * Authors:
+ *   Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * Released under GNU LGPL.  Read the file 'COPYING' for more information.
+ */
+
+#include <set>
+#include <list>
+#include <cassert>
+#include "straightener.h"
+#include <iostream>
+#include <cmath>
+
+using std::set;
+using std::vector;
+using std::list;
+
+namespace straightener {
+
+    // is point p on line a-b?
+    static bool pointOnLine(double px,double py, double ax, double ay, double bx, double by, double& tx) {
+        double dx=bx-ax;
+        double dy=by-ay;
+        double ty=0;
+        if(fabs(dx)<0.0001&&fabs(dy)<0.0001) {
+            // runty line!
+            tx=px-ax;
+            ty=py-ay;
+        } else {
+            if(fabs(dx)<0.0001) {
+                //vertical line
+                if(fabs(px-ax)<0.01) {
+                   tx=(py-ay)/dy;
+                } 
+            } else {
+                tx=(px-ax)/dx;
+            } 
+            if(fabs(dy)<0.0001) {
+                //horizontal line
+                if(fabs(py-ay)<0.01) {
+                   ty=tx;
+                } 
+            } else {
+                ty=(py-ay)/dy;
+            }
+        }
+        //printf("      tx=%f,ty=%f\n",tx,ty);
+        if(fabs(tx-ty)<0.001 && tx>0 && tx<=1) {
+            return true;
+        }
+        return false;
+    }
+    void Edge::nodePath(vector<Node*>& nodes) {
+        list<unsigned> ds(dummyNodes.size());
+        copy(dummyNodes.begin(),dummyNodes.end(),ds.begin());
+        //printf("Edge::nodePath: (%d,%d) dummyNodes:%d\n",startNode,endNode,ds.size());
+        path.clear();
+        path.push_back(startNode);
+        for(unsigned i=1;i<route->n;i++) {
+            //printf("  checking segment %d-%d\n",i-1,i);
+            set<pair<double,unsigned> > pntsOnLineSegment;
+            for(list<unsigned>::iterator j=ds.begin();j!=ds.end();) {
+                double px=nodes[*j]->x;
+                double py=nodes[*j]->y;
+                double ax=route->xs[i-1];
+                double ay=route->ys[i-1];
+                double bx=route->xs[i];
+                double by=route->ys[i];
+                double t=0;
+                list<unsigned>::iterator copyit=j++;
+                //printf("     px=%f, py=%f, ax=%f, ay=%f, bx=%f, by=%f\n",px,py,ax,ay,bx,by);
+                if(pointOnLine(px,py,ax,ay,bx,by,t)) {
+                    //printf(" got node %d\n",*copyit);
+                    pntsOnLineSegment.insert(make_pair(t,*copyit));
+                    ds.erase(copyit);
+                }
+            }
+            for(set<pair<double,unsigned> >::iterator j=pntsOnLineSegment.begin();j!=pntsOnLineSegment.end();j++) {
+                path.push_back(j->second);
+            }
+            //printf("\n");
+        }
+        path.push_back(endNode);
+        assert(ds.empty());
+    }
+
+    typedef enum {Open, Close} EventType;
+    struct Event {
+        EventType type;
+        Node *v;
+        Edge *e;
+        double pos;
+        Event(EventType t, Node *v, double p) : type(t),v(v),e(NULL),pos(p) {};
+        Event(EventType t, Edge *e, double p) : type(t),v(NULL),e(e),pos(p) {};
+    };
+    Event **events;
+    int compare_events(const void *a, const void *b) {
+        Event *ea=*(Event**)a;
+        Event *eb=*(Event**)b;
+        if(ea->v!=NULL&&ea->v==eb->v||ea->e!=NULL&&ea->e==eb->e) {
+            // when comparing opening and closing from object
+            // open must come first
+            if(ea->type==Open) return -1;
+            return 1;
+        } else if(ea->pos > eb->pos) {
+            return 1;
+        } else if(ea->pos < eb->pos) {
+            return -1;
+        }	
+        return 0;
+    }
+
+    void sortNeighbours(Node* v, Node* l, Node* r, 
+            double conjpos, vector<Edge*>& openEdges, 
+            vector<Node*>& L,vector<Node*>& nodes, Dim dim) {
+        double minpos=-DBL_MAX, maxpos=DBL_MAX;
+        if(l!=NULL) {
+            L.push_back(l);
+            minpos=l->scanpos;
+        }
+        typedef pair<double,Edge*> PosEdgePair;
+        set<PosEdgePair> sortedEdges;
+        for(unsigned i=0;i<openEdges.size();i++) {
+            Edge *e=openEdges[i];
+            vector<double> bs;
+            if(dim==HORIZONTAL) {
+                e->xpos(conjpos,bs);
+            } else {
+                e->ypos(conjpos,bs);
+            }
+            //cerr << "edge(intersections="<<bs.size()<<":("<<e->startNode<<","<<e->endNode<<"))"<<endl;
+            for(vector<double>::iterator it=bs.begin();it!=bs.end();it++) {
+                sortedEdges.insert(make_pair(*it,e));
+            }
+        }
+        for(set<PosEdgePair>::iterator i=sortedEdges.begin();i!=sortedEdges.end();i++) {
+            double pos=i->first;
+            if(pos < minpos) continue;
+            if(pos > v->scanpos) break;
+            // if edge is connected (start or end) to v then skip
+            // need to record start and end positions of edge segment!
+            Edge* e=i->second; 
+            if(e->startNode==v->id||e->endNode==v->id) continue;
+            //if(l!=NULL&&(e->startNode==l->id||e->endNode==l->id)) continue;
+            //cerr << "edge("<<e->startNode<<","<<e->endNode<<",pts="<<e->pts<<")"<<endl;
+            Node* d=dim==HORIZONTAL?
+                new Node(nodes.size(),pos,conjpos,e):
+                new Node(nodes.size(),conjpos,pos,e);
+            L.push_back(d);
+            nodes.push_back(d);
+        }
+        L.push_back(v);
+
+        if(r!=NULL) {
+            maxpos=r->scanpos;
+        }
+        for(set<PosEdgePair>::iterator i=sortedEdges.begin();i!=sortedEdges.end();i++) {
+            if(i->first < v->scanpos) continue;
+            if(i->first > maxpos) break;
+            double pos=i->first;
+            // if edge is connected (start or end) to v then skip
+            // need to record start and end positions of edge segment!
+            Edge* e=i->second; 
+            if(e->startNode==v->id||e->endNode==v->id) continue;
+            //if(r!=NULL&&(e->startNode==r->id||e->endNode==r->id)) continue;
+            //cerr << "edge("<<e->startNode<<","<<e->endNode<<",pts="<<e->pts<<")"<<endl;
+            Node* d=dim==HORIZONTAL?
+                new Node(nodes.size(),pos,conjpos,e):
+                new Node(nodes.size(),conjpos,pos,e);
+            L.push_back(d);
+            nodes.push_back(d);
+        }
+        if(r!=NULL) {
+            L.push_back(r);
+        }
+    }
+    static SimpleConstraint* createConstraint(Node* u, Node* v, Dim dim) {
+        double g=dim==HORIZONTAL?(u->width+v->width):(u->height+v->height);
+        g/=2;
+        //cerr << "Constraint: "<< u->id << "+"<<g<<"<="<<v->id<<endl;
+        return new SimpleConstraint(u->id,v->id,g);
+    }
+
+    void generateConstraints(vector<Node*>& nodes, vector<Edge*>& edges,vector<SimpleConstraint*>& cs,Dim dim) {
+        unsigned nevents=2*nodes.size()+2*edges.size();
+        events=new Event*[nevents];
+        unsigned ctr=0;
+        if(dim==HORIZONTAL) {
+            //cout << "Scanning top to bottom..." << endl;
+            for(unsigned i=0;i<nodes.size();i++) {
+                Node *v=nodes[i];
+                v->scanpos=v->x;
+                events[ctr++]=new Event(Open,v,v->ymin+0.01);
+                events[ctr++]=new Event(Close,v,v->ymax-0.01);
+            }
+            for(unsigned i=0;i<edges.size();i++) {
+                Edge *e=edges[i];
+                events[ctr++]=new Event(Open,e,e->ymin-1);
+                events[ctr++]=new Event(Close,e,e->ymax+1);
+            }
+        } else {
+            //cout << "Scanning left to right..." << endl;
+            for(unsigned i=0;i<nodes.size();i++) {
+                Node *v=nodes[i];
+                v->scanpos=v->y;
+                events[ctr++]=new Event(Open,v,v->xmin+0.01);
+                events[ctr++]=new Event(Close,v,v->xmax-0.01);
+            }
+            for(unsigned i=0;i<edges.size();i++) {
+                Edge *e=edges[i];
+                events[ctr++]=new Event(Open,e,e->xmin-1);
+                events[ctr++]=new Event(Close,e,e->xmax+1);
+            }
+        }
+        qsort((Event*)events, (size_t)nevents, sizeof(Event*), compare_events );
+
+        NodeSet openNodes;
+        vector<Edge*> openEdges;
+        for(unsigned i=0;i<nevents;i++) {
+            Event *e=events[i];
+            Node *v=e->v;
+            if(v!=NULL) {
+                v->open = true;
+                //printf("NEvent@%f,nid=%d,(%f,%f),w=%f,h=%f,openn=%d,opene=%d\n",e->pos,v->id,v->x,v->y,v->width,v->height,(int)openNodes.size(),(int)openEdges.size());
+                Node *l=NULL, *r=NULL;
+                if(!openNodes.empty()) {
+                    // it points to the first node to the right of v
+                    NodeSet::iterator it=openNodes.lower_bound(v);
+                    // step left to find the first node to the left of v
+                    if(it--!=openNodes.begin()) {
+                        l=*it;
+                        //printf("l=%d\n",l->id);
+                    }
+                    it=openNodes.upper_bound(v);
+                    if(it!=openNodes.end()) {
+                        r=*it;
+                    }
+                }
+                vector<Node*> L;
+                sortNeighbours(v,l,r,e->pos,openEdges,L,nodes,dim);
+                //printf("L=[");
+                for(unsigned i=0;i<L.size();i++) {
+                    //printf("%d ",L[i]->id);
+                }
+                //printf("]\n");
+                
+                // Case A: create constraints between adjacent edges skipping edges joined
+                // to l,v or r.
+                Node* lastNode=NULL;
+                for(vector<Node*>::iterator i=L.begin();i!=L.end();i++) {
+                    if((*i)->dummy) {
+                        // node is on an edge
+                        Edge *edge=(*i)->edge;
+                        if(!edge->isEnd(v->id)
+                                &&(l!=NULL&&!edge->isEnd(l->id)||l==NULL)
+                                &&(r!=NULL&&!edge->isEnd(r->id)||r==NULL)) {
+                            if(lastNode!=NULL) {
+                                //printf("  Rule A: Constraint: v%d +g <= v%d\n",lastNode->id,(*i)->id);
+                                cs.push_back(createConstraint(lastNode,*i,dim));
+                            }
+                            lastNode=*i;
+                        }
+                    } else {
+                        // is an actual node
+                        lastNode=NULL;
+                    }
+                }
+                // Case B: create constraints for all the edges connected to the right of
+                // their own end, also in the scan line
+                vector<Node*> skipList;
+                lastNode=NULL;
+                for(vector<Node*>::iterator i=L.begin();i!=L.end();i++) {
+                    if((*i)->dummy) {
+                        // node is on an edge
+                        if(lastNode!=NULL) {
+                            if((*i)->edge->isEnd(lastNode->id)) {
+                                skipList.push_back(*i);
+                            } else {
+                                for(vector<Node*>::iterator j=skipList.begin();
+                                        j!=skipList.end();j++) {
+                                    //printf("  Rule B: Constraint: v%d +g <= v%d\n",(*j)->id,(*i)->id);
+                                    cs.push_back(createConstraint(*j,*i,dim));
+                                }
+                                skipList.clear();
+                            }
+                        }
+                    } else {
+                        // is an actual node
+                        skipList.clear();
+                        skipList.push_back(*i);
+                        lastNode=*i;
+                    }
+                }
+                skipList.clear();
+                // Case C: reverse of B
+                lastNode=NULL;
+                for(vector<Node*>::reverse_iterator i=L.rbegin();i!=L.rend();i++) {
+                    if((*i)->dummy) {
+                        // node is on an edge
+                        if(lastNode!=NULL) {
+                            if((*i)->edge->isEnd(lastNode->id)) {
+                                skipList.push_back(*i);
+                            } else {
+                                for(vector<Node*>::iterator j=skipList.begin();
+                                        j!=skipList.end();j++) {
+                                    //printf("  Rule C: Constraint: v%d +g <= v%d\n",(*i)->id,(*j)->id);
+                                    cs.push_back(createConstraint(*i,*j,dim));
+                                }
+                                skipList.clear();
+                            }
+                        }
+                    } else {
+                        // is an actual node
+                        skipList.clear();
+                        skipList.push_back(*i);
+                        lastNode=*i;
+                    }
+                }
+                if(e->type==Close) {
+                    if(l!=NULL) cs.push_back(createConstraint(l,v,dim));
+                    if(r!=NULL) cs.push_back(createConstraint(v,r,dim));
+                }
+            }
+            if(e->type==Open) {
+                if(v!=NULL) {
+                    openNodes.insert(v);
+                } else {
+                    //printf("EdgeOpen@%f,eid=%d,(u,v)=(%d,%d)\n", e->pos,e->e->id,e->e->startNode,e->e->endNode);
+                    e->e->openInd=openEdges.size();
+                    openEdges.push_back(e->e);
+                }
+            } else {
+                // Close
+                if(v!=NULL) {
+                    openNodes.erase(v);
+                    v->open=false;
+                } else {
+                    //printf("EdgeClose@%f,eid=%d,(u,v)=(%d,%d)\n", e->pos,e->e->id,e->e->startNode,e->e->endNode);
+                    unsigned i=e->e->openInd;
+                    openEdges[i]=openEdges[openEdges.size()-1];
+                    openEdges[i]->openInd=i;
+                    openEdges.resize(openEdges.size()-1);
+                }
+            }
+            delete e;
+        }
+        delete [] events;
+    }
+}
+
diff --git a/src/libcola/straightener.h b/src/libcola/straightener.h
new file mode 100644
index 000000000..33af0c697
--- /dev/null
+++ b/src/libcola/straightener.h
@@ -0,0 +1,133 @@
+/*
+** vim: set cindent 
+** vim: ts=4 sw=4 et tw=0 wm=0
+*/
+#ifndef STRAIGHTENER_H
+#define STRAIGHTENER_H
+#include <set>
+#include <libvpsc/generate-constraints.h>
+#include "gradient_projection.h"
+namespace straightener {
+    struct Route {
+        Route(unsigned n) : n(n), xs(new double[n]), ys(new double[n]) {}
+        ~Route() {
+            delete [] xs;
+            delete [] ys;
+        }
+        void boundingBox(double &xmin,double &ymin,double &xmax,double &ymax) {
+            xmin=ymin=DBL_MAX;
+            xmax=ymax=-DBL_MAX;
+            for(unsigned i=0;i<n;i++) {
+                xmin=min(xmin,xs[i]);
+                xmax=max(xmax,xs[i]);
+                ymin=min(ymin,ys[i]);
+                ymax=max(ymax,ys[i]);
+            } 
+        }
+        unsigned n;
+        double *xs;
+        double *ys;
+    };
+    class Node;
+    struct Edge {
+        unsigned id;
+        unsigned openInd; // position in openEdges
+        unsigned startNode, endNode;
+        Route* route;
+        double xmin, xmax, ymin, ymax;
+        vector<unsigned> dummyNodes;
+        vector<unsigned> path;
+        Edge(unsigned id, unsigned start, unsigned end, Route* route)
+        : id(id), startNode(start), endNode(end), route(route)
+        {
+            route->boundingBox(xmin,ymin,xmax,ymax);
+        }
+        ~Edge() {
+            delete route;
+        }
+        void setRoute(Route* r) {
+            delete route;
+            route=r;
+            route->boundingBox(xmin,ymin,xmax,ymax);
+        }
+        bool isEnd(unsigned n) {
+            if(startNode==n||endNode==n) return true;
+            return false;
+        }
+        void nodePath(vector<Node*>& nodes);
+        void createRouteFromPath(double* X, double* Y) {
+            Route* r=new Route(path.size());
+            for(unsigned i=0;i<path.size();i++) {
+                r->xs[i]=X[path[i]];
+                r->ys[i]=Y[path[i]];
+            }
+            setRoute(r);
+        }
+        void xpos(double y, vector<double>& xs) {
+            // search line segments for intersection points with y pos
+            for(unsigned i=1;i<route->n;i++) {
+                double ax=route->xs[i-1], bx=route->xs[i], ay=route->ys[i-1], by=route->ys[i];
+                double r=(y-ay)/(by-ay);
+                // as long as y is between ay and by then r>0
+                if(r>0&&r<=1) {
+                    xs.push_back(ax+(bx-ax)*r);
+                }
+            }
+        }
+        void ypos(double x, vector<double>& ys) {
+            // search line segments for intersection points with x pos
+            for(unsigned i=1;i<route->n;i++) {
+                double ax=route->xs[i-1], bx=route->xs[i], ay=route->ys[i-1], by=route->ys[i];
+                double r=(x-ax)/(bx-ax);
+                // as long as y is between ax and bx then r>0
+                if(r>0&&r<=1) {
+                    ys.push_back(ay+(by-ay)*r);
+                }
+            }
+        }
+    };
+    class Node {
+    public:
+        unsigned id;
+        double x,y;
+        double scanpos;
+        double width, height;
+        double xmin, xmax, ymin, ymax;
+        Edge *edge;
+        bool dummy;
+        double weight;
+        bool open;
+        Node(unsigned id, Rectangle* r) :
+            id(id),x(r->getCentreX()),y(r->getCentreY()), width(r->width()), height(r->height()),
+            xmin(x-width/2),xmax(x+width/2),
+            ymin(y-height/2),ymax(y+height/2),
+            edge(NULL),dummy(false),weight(-0.1),open(false) { }
+    private:
+        friend void sortNeighbours(Node* v, Node* l, Node* r, 
+            double conjpos, vector<Edge*>& openEdges, 
+            vector<Node*>& L,vector<Node*>& nodes, Dim dim);
+        Node(unsigned id, double x, double y, Edge* e) : 
+            id(id),x(x),y(y), width(4), height(width),
+            xmin(x-width/2),xmax(x+width/2),
+            ymin(y-height/2),ymax(y+height/2),
+            edge(e),dummy(true),weight(-0.1)  {
+                e->dummyNodes.push_back(id);
+            }
+    };
+    struct CmpNodePos {
+        bool operator() (const Node* u, const Node* v) const {
+            if (u->scanpos < v->scanpos) {
+                return true;
+            }
+            if (v->scanpos < u->scanpos) {
+                return false;
+            }
+            return u < v;
+        }
+    };
+    typedef std::set<Node*,CmpNodePos> NodeSet;
+    void generateConstraints(vector<Node*>& nodes, vector<Edge*>& edges,vector<SimpleConstraint*>& cs, Dim dim);
+    void nodePath(Edge& e,vector<Node*>& nodes, vector<unsigned>& path);
+}
+
+#endif