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)

46 files changed, 3131 insertions, 588 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 090b4c683..6358d35dd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -43,6 +43,8 @@ include libnr/Makefile_insert
 include libnrtype/Makefile_insert
 include libavoid/Makefile_insert
 include livarot/Makefile_insert
+include libvpsc/Makefile_insert
+include libcola/Makefile_insert
 include removeoverlap/Makefile_insert
 include graphlayout/Makefile_insert
 include svg/Makefile_insert
@@ -86,6 +88,8 @@ noinst_LIBRARIES =	\
 	libnr/libnr.a	\
 	libnrtype/libnrtype.a	\
 	libavoid/libavoid.a	\
+	libvpsc/libvpsc.a	\
+	libcola/libcola.a	\
 	livarot/libvarot.a	\
 	removeoverlap/libremoveoverlap.a	\
 	graphlayout/libgraphlayout.a	\
diff --git a/src/Makefile_insert b/src/Makefile_insert
index 4b797f233..b136a55e6 100644
--- a/src/Makefile_insert
+++ b/src/Makefile_insert
@@ -279,6 +279,8 @@ inkscape_private_libs =	\
 	ui/widget/libuiwidget.a	\
 	graphlayout/libgraphlayout.a	\
 	removeoverlap/libremoveoverlap.a		\
+	libcola/libcola.a	\
+	libvpsc/libvpsc.a	\
 	extension/libextension.a	\
 	extension/implementation/libimplementation.a	\
 	extension/internal/libinternal.a	\
diff --git a/src/graphlayout/graphlayout.cpp b/src/graphlayout/graphlayout.cpp
index 00829dac8..42b867a33 100644
--- a/src/graphlayout/graphlayout.cpp
+++ b/src/graphlayout/graphlayout.cpp
@@ -14,7 +14,6 @@
 #include <iostream>
 #include <config.h>
 
-#ifdef HAVE_BOOST_GRAPH_LIB
 #include "sp-path.h"
 #include "sp-item.h"
 #include "sp-item-transform.h"
@@ -22,24 +21,15 @@
 #include "conn-avoid-ref.h"
 #include "libavoid/connector.h"
 #include "libavoid/geomtypes.h"
-#include <boost/graph/kamada_kawai_spring_layout.hpp>
-#include <boost/graph/circle_layout.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/graphviz.hpp>
+#include "libcola/cola.h"
+#include "libvpsc/generate-constraints.h"
 #include <map>
 #include <vector>
 #include <algorithm>
 #include <float.h>
 
-using namespace boost;
-
-// create a typedef for the Graph type
-typedef adjacency_list<vecS, vecS, undirectedS, no_property, 
-	property<edge_weight_t, double> > Graph;
-typedef property_map<Graph, edge_weight_t>::type WeightMap;
-typedef graph_traits<Graph>::vertex_descriptor Vertex;
-typedef std::vector<Avoid::Point> PositionVec;
-typedef iterator_property_map<PositionVec::iterator, property_map<Graph, vertex_index_t>::type> PositionMap;
+using namespace std;
+using namespace cola;
 
 /**
  * Returns true if item is a connector
@@ -56,7 +46,7 @@ bool isConnector(SPItem const *const i) {
  * Scans the items list and places those items that are 
  * not connectors in filtered
  */
-void filterConnectors(GSList const *const items, std::list<SPItem *> &filtered) {
+void filterConnectors(GSList const *const items, list<SPItem *> &filtered) {
 	for(GSList *i=(GSList *)items; i!=NULL; i=i->next) {
 		SPItem *item=SP_ITEM(i->data);
 		if(!isConnector(item)) {
@@ -75,90 +65,79 @@ void graphlayout(GSList const *const items) {
 	}
 
 	using Inkscape::Util::GSListConstIterator;
-	std::list<SPItem *> selected;
+	list<SPItem *> selected;
 	filterConnectors(items,selected);
 	if (selected.empty()) return;
 
-	int n=selected.size();
-	std::cout<<"|V|="<<n<<std::endl;
-
-	Graph g;
+	const unsigned n=selected.size();
+	cout<<"|V|="<<n<<endl;
+	//Check 2 or more selected objects
+	if (n < 2) return;
 
 	double minX=DBL_MAX, minY=DBL_MAX, maxX=-DBL_MAX, maxY=-DBL_MAX;
 
-	std::map<std::string,Vertex> nodelookup;
-	for (std::list<SPItem *>::iterator i(selected.begin());
+	map<string,unsigned> nodelookup;
+	vector<Rectangle*> rs;
+	vector<Edge> es;
+	for (list<SPItem *>::iterator i(selected.begin());
 		i != selected.end();
 		++i)
 	{
 		SPItem *u=*i;
-		std::cout<<"Creating node for id: "<<u->id<<std::endl;
-		nodelookup[u->id]=add_vertex(g);
+		NR::Rect const item_box(sp_item_bbox_desktop(u));
+		NR::Point ll(item_box.min());
+		NR::Point ur(item_box.max());
+		minX=min(ll[0],minX); minY=min(ll[1],minY);
+	        maxX=max(ur[0],maxX); maxY=max(ur[1],maxY);
+		cout<<"Creating node for id: "<<u->id<<endl;
+		nodelookup[u->id]=rs.size();
+		rs.push_back(new Rectangle(ll[0],ur[0],ll[1],ur[1]));
 	}
 
-	//Check 2 or more selected objects
-	if (n < 2) return;
 
-	WeightMap weightmap=get(edge_weight, g);
-	for (std::list<SPItem *>::iterator i(selected.begin());
+	for (list<SPItem *>::iterator i(selected.begin());
 		i != selected.end();
 		++i)
 	{
-		using NR::X; using NR::Y;
 		SPItem *iu=*i;
-		std::cout<<"Getting neighbours for id: "<<iu->id<<std::endl;
-		Vertex u=nodelookup[iu->id];
+		cout<<"Getting neighbours for id: "<<iu->id<<endl;
+		unsigned u=nodelookup[iu->id];
 		GSList *nlist=iu->avoidRef->getAttachedShapes(Avoid::runningFrom);
-		std::list<SPItem *> neighbours;
+		list<SPItem *> neighbours;
 		neighbours.insert<GSListConstIterator<SPItem *> >(neighbours.end(),nlist,NULL);
-		for (std::list<SPItem *>::iterator j(neighbours.begin());
+		for (list<SPItem *>::iterator j(neighbours.begin());
 				j != neighbours.end();
 				++j) {
 			
 			SPItem *iv=*j;
-			Vertex v=nodelookup[iv->id];
-			Graph::edge_descriptor e; bool inserted;
-			tie(e, inserted)=add_edge(u,v,g);
-			weightmap[e]=1.0;
+			// What do we do if iv not in nodelookup?!?!
+			unsigned v=nodelookup[iv->id];
+			es.push_back(make_pair(u,v));
 		}
 		if(nlist) {
 			g_slist_free(nlist);
 		}
-		NR::Rect const item_box(sp_item_bbox_desktop(*i));
-			
-		NR::Point ll(item_box.min());
-		minX=std::min(ll[0],minX);
-		minY=std::min(ll[1],minY);
-		NR::Point ur(item_box.max());
-		maxX=std::max(ur[0],maxX);
-		maxY=std::max(ur[1],maxY);
 	}
 	double width=maxX-minX;
 	double height=maxY-minY;
-	std::cout<<"Graph has |V|="<<num_vertices(g)<<" Width="<<width<<" Height="<<height<<std::endl;
-  	PositionVec position_vec(num_vertices(g));
-  	PositionMap position(position_vec.begin(), get(vertex_index, g));
-  	write_graphviz(std::cout, g);
-	circle_graph_layout(g, position, width/2.0);
-	kamada_kawai_spring_layout(g, position, weightmap, side_length(width));
+	const unsigned E = es.size();
+	double eweights[E];
+	fill(eweights,eweights+E,1);
 
-	graph_traits<Graph>::vertex_iterator vi, vi_end;
-	for (std::list<SPItem *>::iterator it(selected.begin());
+	ConstrainedMajorizationLayout alg(rs,es,eweights,width/n);
+	alg.run();
+	
+	for (list<SPItem *>::iterator it(selected.begin());
 		it != selected.end();
 		++it)
 	{
 		SPItem *u=*it;
 		if(!isConnector(u)) {
+			Rectangle* r=rs[nodelookup[u->id]];
 			NR::Rect const item_box(sp_item_bbox_desktop(u));
 			NR::Point const curr(item_box.midpoint());
-			NR::Point const dest(minX+width/2.0+position[nodelookup[u->id]].x,
-				     	minY+height/2.0+position[nodelookup[u->id]].y);
+			NR::Point const dest(r->getCentreX(),r->getCentreY());
 			sp_item_move_rel(u, NR::translate(dest - curr));
 		}
 	}
 }
-#else
-void graphlayout(GSList const *const items) {
-	std::cout<<"Connector network layout not available!  Install boost graph library and recompile to enable."<<std::endl;
-}
-#endif // HAVE_BOOST_GRAPH_LIB
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
diff --git a/src/libvpsc/COPYING b/src/libvpsc/COPYING
new file mode 100644
index 000000000..6ff1d7b6f
--- /dev/null
+++ b/src/libvpsc/COPYING
@@ -0,0 +1,505 @@
+		  GNU LESSER GENERAL PUBLIC LICENSE
+		       Version 2.1, February 1999
+
+ Copyright (C) 1991, 1999 Free Software Foundation, Inc.
+     51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+[This is the first released version of the Lesser GPL.  It also counts
+ as the successor of the GNU Library Public License, version 2, hence
+ the version number 2.1.]
+
+			    Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+Licenses are intended to guarantee your freedom to share and change
+free software--to make sure the software is free for all its users.
+
+  This license, the Lesser General Public License, applies to some
+specially designated software packages--typically libraries--of the
+Free Software Foundation and other authors who decide to use it.  You
+can use it too, but we suggest you first think carefully about whether
+this license or the ordinary General Public License is the better
+strategy to use in any particular case, based on the explanations below.
+
+  When we speak of free software, we are referring to freedom of use,
+not price.  Our General Public Licenses are designed to make sure that
+you have the freedom to distribute copies of free software (and charge
+for this service if you wish); that you receive source code or can get
+it if you want it; that you can change the software and use pieces of
+it in new free programs; and that you are informed that you can do
+these things.
+
+  To protect your rights, we need to make restrictions that forbid
+distributors to deny you these rights or to ask you to surrender these
+rights.  These restrictions translate to certain responsibilities for
+you if you distribute copies of the library or if you modify it.
+
+  For example, if you distribute copies of the library, whether gratis
+or for a fee, you must give the recipients all the rights that we gave
+you.  You must make sure that they, too, receive or can get the source
+code.  If you link other code with the library, you must provide
+complete object files to the recipients, so that they can relink them
+with the library after making changes to the library and recompiling
+it.  And you must show them these terms so they know their rights.
+
+  We protect your rights with a two-step method: (1) we copyright the
+library, and (2) we offer you this license, which gives you legal
+permission to copy, distribute and/or modify the library.
+
+  To protect each distributor, we want to make it very clear that
+there is no warranty for the free library.  Also, if the library is
+modified by someone else and passed on, the recipients should know
+that what they have is not the original version, so that the original
+author's reputation will not be affected by problems that might be
+introduced by others.
+
+  Finally, software patents pose a constant threat to the existence of
+any free program.  We wish to make sure that a company cannot
+effectively restrict the users of a free program by obtaining a
+restrictive license from a patent holder.  Therefore, we insist that
+any patent license obtained for a version of the library must be
+consistent with the full freedom of use specified in this license.
+
+  Most GNU software, including some libraries, is covered by the
+ordinary GNU General Public License.  This license, the GNU Lesser
+General Public License, applies to certain designated libraries, and
+is quite different from the ordinary General Public License.  We use
+this license for certain libraries in order to permit linking those
+libraries into non-free programs.
+
+  When a program is linked with a library, whether statically or using
+a shared library, the combination of the two is legally speaking a
+combined work, a derivative of the original library.  The ordinary
+General Public License therefore permits such linking only if the
+entire combination fits its criteria of freedom.  The Lesser General
+Public License permits more lax criteria for linking other code with
+the library.
+
+  We call this license the "Lesser" General Public License because it
+does Less to protect the user's freedom than the ordinary General
+Public License.  It also provides other free software developers Less
+of an advantage over competing non-free programs.  These disadvantages
+are the reason we use the ordinary General Public License for many
+libraries.  However, the Lesser license provides advantages in certain
+special circumstances.
+
+  For example, on rare occasions, there may be a special need to
+encourage the widest possible use of a certain library, so that it becomes
+a de-facto standard.  To achieve this, non-free programs must be
+allowed to use the library.  A more frequent case is that a free
+library does the same job as widely used non-free libraries.  In this
+case, there is little to gain by limiting the free library to free
+software only, so we use the Lesser General Public License.
+
+  In other cases, permission to use a particular library in non-free
+programs enables a greater number of people to use a large body of
+free software.  For example, permission to use the GNU C Library in
+non-free programs enables many more people to use the whole GNU
+operating system, as well as its variant, the GNU/Linux operating
+system.
+
+  Although the Lesser General Public License is Less protective of the
+users' freedom, it does ensure that the user of a program that is
+linked with the Library has the freedom and the wherewithal to run
+that program using a modified version of the Library.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.  Pay close attention to the difference between a
+"work based on the library" and a "work that uses the library".  The
+former contains code derived from the library, whereas the latter must
+be combined with the library in order to run.
+
+		  GNU LESSER GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License Agreement applies to any software library or other
+program which contains a notice placed by the copyright holder or
+other authorized party saying it may be distributed under the terms of
+this Lesser General Public License (also called "this License").
+Each licensee is addressed as "you".
+
+  A "library" means a collection of software functions and/or data
+prepared so as to be conveniently linked with application programs
+(which use some of those functions and data) to form executables.
+
+  The "Library", below, refers to any such software library or work
+which has been distributed under these terms.  A "work based on the
+Library" means either the Library or any derivative work under
+copyright law: that is to say, a work containing the Library or a
+portion of it, either verbatim or with modifications and/or translated
+straightforwardly into another language.  (Hereinafter, translation is
+included without limitation in the term "modification".)
+
+  "Source code" for a work means the preferred form of the work for
+making modifications to it.  For a library, complete source code means
+all the source code for all modules it contains, plus any associated
+interface definition files, plus the scripts used to control compilation
+and installation of the library.
+
+  Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running a program using the Library is not restricted, and output from
+such a program is covered only if its contents constitute a work based
+on the Library (independent of the use of the Library in a tool for
+writing it).  Whether that is true depends on what the Library does
+and what the program that uses the Library does.
+  
+  1. You may copy and distribute verbatim copies of the Library's
+complete source code as you receive it, in any medium, provided that
+you conspicuously and appropriately publish on each copy an
+appropriate copyright notice and disclaimer of warranty; keep intact
+all the notices that refer to this License and to the absence of any
+warranty; and distribute a copy of this License along with the
+Library.
+
+  You may charge a fee for the physical act of transferring a copy,
+and you may at your option offer warranty protection in exchange for a
+fee.
+
+  2. You may modify your copy or copies of the Library or any portion
+of it, thus forming a work based on the Library, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) The modified work must itself be a software library.
+
+    b) You must cause the files modified to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    c) You must cause the whole of the work to be licensed at no
+    charge to all third parties under the terms of this License.
+
+    d) If a facility in the modified Library refers to a function or a
+    table of data to be supplied by an application program that uses
+    the facility, other than as an argument passed when the facility
+    is invoked, then you must make a good faith effort to ensure that,
+    in the event an application does not supply such function or
+    table, the facility still operates, and performs whatever part of
+    its purpose remains meaningful.
+
+    (For example, a function in a library to compute square roots has
+    a purpose that is entirely well-defined independent of the
+    application.  Therefore, Subsection 2d requires that any
+    application-supplied function or table used by this function must
+    be optional: if the application does not supply it, the square
+    root function must still compute square roots.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Library,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Library, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote
+it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Library.
+
+In addition, mere aggregation of another work not based on the Library
+with the Library (or with a work based on the Library) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may opt to apply the terms of the ordinary GNU General Public
+License instead of this License to a given copy of the Library.  To do
+this, you must alter all the notices that refer to this License, so
+that they refer to the ordinary GNU General Public License, version 2,
+instead of to this License.  (If a newer version than version 2 of the
+ordinary GNU General Public License has appeared, then you can specify
+that version instead if you wish.)  Do not make any other change in
+these notices.
+
+  Once this change is made in a given copy, it is irreversible for
+that copy, so the ordinary GNU General Public License applies to all
+subsequent copies and derivative works made from that copy.
+
+  This option is useful when you wish to copy part of the code of
+the Library into a program that is not a library.
+
+  4. You may copy and distribute the Library (or a portion or
+derivative of it, under Section 2) in object code or executable form
+under the terms of Sections 1 and 2 above provided that you accompany
+it with the complete corresponding machine-readable source code, which
+must be distributed under the terms of Sections 1 and 2 above on a
+medium customarily used for software interchange.
+
+  If distribution of object code is made by offering access to copy
+from a designated place, then offering equivalent access to copy the
+source code from the same place satisfies the requirement to
+distribute the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  5. A program that contains no derivative of any portion of the
+Library, but is designed to work with the Library by being compiled or
+linked with it, is called a "work that uses the Library".  Such a
+work, in isolation, is not a derivative work of the Library, and
+therefore falls outside the scope of this License.
+
+  However, linking a "work that uses the Library" with the Library
+creates an executable that is a derivative of the Library (because it
+contains portions of the Library), rather than a "work that uses the
+library".  The executable is therefore covered by this License.
+Section 6 states terms for distribution of such executables.
+
+  When a "work that uses the Library" uses material from a header file
+that is part of the Library, the object code for the work may be a
+derivative work of the Library even though the source code is not.
+Whether this is true is especially significant if the work can be
+linked without the Library, or if the work is itself a library.  The
+threshold for this to be true is not precisely defined by law.
+
+  If such an object file uses only numerical parameters, data
+structure layouts and accessors, and small macros and small inline
+functions (ten lines or less in length), then the use of the object
+file is unrestricted, regardless of whether it is legally a derivative
+work.  (Executables containing this object code plus portions of the
+Library will still fall under Section 6.)
+
+  Otherwise, if the work is a derivative of the Library, you may
+distribute the object code for the work under the terms of Section 6.
+Any executables containing that work also fall under Section 6,
+whether or not they are linked directly with the Library itself.
+
+  6. As an exception to the Sections above, you may also combine or
+link a "work that uses the Library" with the Library to produce a
+work containing portions of the Library, and distribute that work
+under terms of your choice, provided that the terms permit
+modification of the work for the customer's own use and reverse
+engineering for debugging such modifications.
+
+  You must give prominent notice with each copy of the work that the
+Library is used in it and that the Library and its use are covered by
+this License.  You must supply a copy of this License.  If the work
+during execution displays copyright notices, you must include the
+copyright notice for the Library among them, as well as a reference
+directing the user to the copy of this License.  Also, you must do one
+of these things:
+
+    a) Accompany the work with the complete corresponding
+    machine-readable source code for the Library including whatever
+    changes were used in the work (which must be distributed under
+    Sections 1 and 2 above); and, if the work is an executable linked
+    with the Library, with the complete machine-readable "work that
+    uses the Library", as object code and/or source code, so that the
+    user can modify the Library and then relink to produce a modified
+    executable containing the modified Library.  (It is understood
+    that the user who changes the contents of definitions files in the
+    Library will not necessarily be able to recompile the application
+    to use the modified definitions.)
+
+    b) Use a suitable shared library mechanism for linking with the
+    Library.  A suitable mechanism is one that (1) uses at run time a
+    copy of the library already present on the user's computer system,
+    rather than copying library functions into the executable, and (2)
+    will operate properly with a modified version of the library, if
+    the user installs one, as long as the modified version is
+    interface-compatible with the version that the work was made with.
+
+    c) Accompany the work with a written offer, valid for at
+    least three years, to give the same user the materials
+    specified in Subsection 6a, above, for a charge no more
+    than the cost of performing this distribution.
+
+    d) If distribution of the work is made by offering access to copy
+    from a designated place, offer equivalent access to copy the above
+    specified materials from the same place.
+
+    e) Verify that the user has already received a copy of these
+    materials or that you have already sent this user a copy.
+
+  For an executable, the required form of the "work that uses the
+Library" must include any data and utility programs needed for
+reproducing the executable from it.  However, as a special exception,
+the materials to be distributed need not include anything that is
+normally distributed (in either source or binary form) with the major
+components (compiler, kernel, and so on) of the operating system on
+which the executable runs, unless that component itself accompanies
+the executable.
+
+  It may happen that this requirement contradicts the license
+restrictions of other proprietary libraries that do not normally
+accompany the operating system.  Such a contradiction means you cannot
+use both them and the Library together in an executable that you
+distribute.
+
+  7. You may place library facilities that are a work based on the
+Library side-by-side in a single library together with other library
+facilities not covered by this License, and distribute such a combined
+library, provided that the separate distribution of the work based on
+the Library and of the other library facilities is otherwise
+permitted, and provided that you do these two things:
+
+    a) Accompany the combined library with a copy of the same work
+    based on the Library, uncombined with any other library
+    facilities.  This must be distributed under the terms of the
+    Sections above.
+
+    b) Give prominent notice with the combined library of the fact
+    that part of it is a work based on the Library, and explaining
+    where to find the accompanying uncombined form of the same work.
+
+  8. You may not copy, modify, sublicense, link with, or distribute
+the Library except as expressly provided under this License.  Any
+attempt otherwise to copy, modify, sublicense, link with, or
+distribute the Library is void, and will automatically terminate your
+rights under this License.  However, parties who have received copies,
+or rights, from you under this License will not have their licenses
+terminated so long as such parties remain in full compliance.
+
+  9. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Library or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Library (or any work based on the
+Library), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Library or works based on it.
+
+  10. Each time you redistribute the Library (or any work based on the
+Library), the recipient automatically receives a license from the
+original licensor to copy, distribute, link with or modify the Library
+subject to these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties with
+this License.
+
+  11. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Library at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Library by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Library.
+
+If any portion of this section is held invalid or unenforceable under any
+particular circumstance, the balance of the section is intended to apply,
+and the section as a whole is intended to apply in other circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  12. If the distribution and/or use of the Library is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Library under this License may add
+an explicit geographical distribution limitation excluding those countries,
+so that distribution is permitted only in or among countries not thus
+excluded.  In such case, this License incorporates the limitation as if
+written in the body of this License.
+
+  13. The Free Software Foundation may publish revised and/or new
+versions of the Lesser General Public License from time to time.
+Such new versions will be similar in spirit to the present version,
+but may differ in detail to address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Library
+specifies a version number of this License which applies to it and
+"any later version", you have the option of following the terms and
+conditions either of that version or of any later version published by
+the Free Software Foundation.  If the Library does not specify a
+license version number, you may choose any version ever published by
+the Free Software Foundation.
+
+  14. If you wish to incorporate parts of the Library into other free
+programs whose distribution conditions are incompatible with these,
+write to the author to ask for permission.  For software which is
+copyrighted by the Free Software Foundation, write to the Free
+Software Foundation; we sometimes make exceptions for this.  Our
+decision will be guided by the two goals of preserving the free status
+of all derivatives of our free software and of promoting the sharing
+and reuse of software generally.
+
+			    NO WARRANTY
+
+  15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO
+WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW.
+EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR
+OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY
+KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE
+LIBRARY IS WITH YOU.  SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME
+THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN
+WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY
+AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU
+FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR
+CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE
+LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING
+RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A
+FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF
+SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH
+DAMAGES.
+
+		     END OF TERMS AND CONDITIONS
+
+           How to Apply These Terms to Your New Libraries
+
+  If you develop a new library, and you want it to be of the greatest
+possible use to the public, we recommend making it free software that
+everyone can redistribute and change.  You can do so by permitting
+redistribution under these terms (or, alternatively, under the terms of the
+ordinary General Public License).
+
+  To apply these terms, attach the following notices to the library.  It is
+safest to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least the
+"copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the library's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This library is free software; you can redistribute it and/or
+    modify it under the terms of the GNU Lesser General Public
+    License as published by the Free Software Foundation; either
+    version 2.1 of the License, or (at your option) any later version.
+
+    This library is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public
+    License along with this library; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+Also add information on how to contact you by electronic and paper mail.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the library, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the
+  library `Frob' (a library for tweaking knobs) written by James Random Hacker.
+
+  <signature of Ty Coon>, 1 April 1990
+  Ty Coon, President of Vice
+
+That's all there is to it!
+
+
+
diff --git a/src/libvpsc/Makefile_insert b/src/libvpsc/Makefile_insert
new file mode 100644
index 000000000..78b825b13
--- /dev/null
+++ b/src/libvpsc/Makefile_insert
@@ -0,0 +1,26 @@
+## Makefile.am fragment sourced by src/Makefile.am.
+libvpsc/all: libvpsc/libvpsc.a
+
+libvpsc/clean:
+	rm -f libvpsc/libvpsc.a $(libvpsc_libvpsc_a_OBJECTS)
+
+libvpsc_libvpsc_a_SOURCES = libvpsc/block.cpp\
+	libvpsc/blocks.cpp\
+	libvpsc/constraint.cpp\
+	libvpsc/generate-constraints.cpp\
+	libvpsc/pairingheap/PairingHeap.cpp\
+	libvpsc/remove_rectangle_overlap.cpp\
+	libvpsc/solve_VPSC.cpp\
+	libvpsc/csolve_VPSC.cpp\
+	libvpsc/variable.cpp\
+	libvpsc/isnan.h\
+	libvpsc/block.h\
+	libvpsc/blocks.h\
+	libvpsc/constraint.h\
+	libvpsc/generate-constraints.h\
+	libvpsc/pairingheap/PairingHeap.h\
+	libvpsc/pairingheap/dsexceptions.h\
+	libvpsc/remove_rectangle_overlap.h\
+	libvpsc/solve_VPSC.h\
+	libvpsc/csolve_VPSC.h\
+	libvpsc/variable.h
diff --git a/src/removeoverlap/block.cpp b/src/libvpsc/block.cpp
index 042d9fc3c..69a439cd7 100644
--- a/src/removeoverlap/block.cpp
+++ b/src/libvpsc/block.cpp
@@ -23,16 +23,14 @@ using std::endl;
 #endif
 using std::vector;
 
-typedef vector<Constraint*>::iterator Cit;
-
-void Block::addVariable(Variable *v) {
+void Block::addVariable(Variable* const v) {
 	v->block=this;
 	vars->push_back(v);
 	weight+=v->weight;
 	wposn += v->weight * (v->desiredPosition - v->offset);
 	posn=wposn/weight;
 }
-Block::Block(Variable *v) {
+Block::Block(Variable* const v) {
 	timeStamp=0;
 	posn=weight=wposn=0;
 	in=NULL;
@@ -47,7 +45,7 @@ Block::Block(Variable *v) {
 
 double Block::desiredWeightedPosition() {
 	double wp = 0;
-	for (vector<Variable*>::iterator v=vars->begin();v!=vars->end();++v) {
+	for (Vit v=vars->begin();v!=vars->end();++v) {
 		wp += ((*v)->desiredPosition - (*v)->offset) * (*v)->weight;
 	}
 	return wp;
@@ -67,7 +65,7 @@ void Block::setUpOutConstraints() {
 void Block::setUpConstraintHeap(PairingHeap<Constraint*>* &h,bool in) {
 	delete h;
 	h = new PairingHeap<Constraint*>(&compareConstraints);
-	for (vector<Variable*>::iterator i=vars->begin();i!=vars->end();++i) {
+	for (Vit i=vars->begin();i!=vars->end();++i) {
 		Variable *v=*i;
 		vector<Constraint*> *cs=in?&(v->in):&(v->out);
 		for (Cit j=cs->begin();j!=cs->end();++j) {
@@ -112,7 +110,7 @@ void Block::merge(Block *b, Constraint *c, double dist) {
 	wposn+=b->wposn-dist*b->weight;
 	weight+=b->weight;
 	posn=wposn/weight;
-	for(vector<Variable*>::iterator i=b->vars->begin();i!=b->vars->end();++i) {
+	for(Vit i=b->vars->begin();i!=b->vars->end();++i) {
 		Variable *v=*i;
 		v->block=this;
 		v->offset+=dist;
@@ -209,10 +207,10 @@ void Block::deleteMinInConstraint() {
 void Block::deleteMinOutConstraint() {
 	out->deleteMin();
 }
-inline bool Block::canFollowLeft(Constraint *c, Variable *last) {
+inline bool Block::canFollowLeft(Constraint *c, const Variable* const last) {
 	return c->left->block==this && c->active && last!=c->left;
 }
-inline bool Block::canFollowRight(Constraint *c, Variable *last) {
+inline bool Block::canFollowRight(Constraint *c, const Variable* const last) {
 	return c->right->block==this && c->active && last!=c->right;
 }
 
@@ -221,7 +219,8 @@ inline bool Block::canFollowRight(Constraint *c, Variable *last) {
 // Does not backtrack over u.
 // also records the constraint with minimum lagrange multiplier
 // in min_lm
-double Block::compute_dfdv(Variable *v, Variable *u, Constraint *&min_lm) {
+double Block::compute_dfdv(Variable* const v, Variable* const u,
+	       	Constraint *&min_lm) {
 	double dfdv=v->weight*(v->position() - v->desiredPosition);
 	for(Cit it=v->out.begin();it!=v->out.end();++it) {
 		Constraint *c=*it;
@@ -254,8 +253,9 @@ double Block::compute_dfdv(Variable *v, Variable *u, Constraint *&min_lm) {
 // the recursion (checking only constraints traversed left-to-right 
 // in order to avoid creating any new violations).
 // We also do not consider equality constraints as potential split points
-Block::Pair Block::compute_dfdv_between(Variable* r, Variable* v, Variable* u, 
-		Direction dir = NONE, bool changedDirection = false) {
+Block::Pair Block::compute_dfdv_between(
+		Variable* r, Variable* const v, Variable* const u, 
+		const Direction dir = NONE, bool changedDirection = false) {
 	double dfdv=v->weight*(v->position() - v->desiredPosition);
 	Constraint *m=NULL;
 	for(Cit it(v->in.begin());it!=v->in.end();++it) {
@@ -300,7 +300,7 @@ Block::Pair Block::compute_dfdv_between(Variable* r, Variable* v, Variable* u,
 // resets LMs for all active constraints to 0 by
 // traversing active constraint tree starting from v,
 // not back tracking over u
-void Block::reset_active_lm(Variable *v, Variable *u) {
+void Block::reset_active_lm(Variable* const v, Variable* const u) {
 	for(Cit it=v->out.begin();it!=v->out.end();++it) {
 		Constraint *c=*it;
 		if(canFollowRight(c,u)) {
@@ -326,7 +326,7 @@ Constraint *Block::findMinLM() {
 	compute_dfdv(vars->front(),NULL,min_lm);
 	return min_lm;
 }
-Constraint *Block::findMinLMBetween(Variable* lv, Variable* rv) {
+Constraint *Block::findMinLMBetween(Variable* const lv, Variable* const rv) {
 	Constraint *min_lm=NULL;
 	reset_active_lm(vars->front(),NULL);
 	min_lm=compute_dfdv_between(rv,lv,NULL).second;
@@ -335,7 +335,7 @@ Constraint *Block::findMinLMBetween(Variable* lv, Variable* rv) {
 
 // populates block b by traversing the active constraint tree adding variables as they're 
 // visited.  Starts from variable v and does not backtrack over variable u.
-void Block::populateSplitBlock(Block *b, Variable *v, Variable *u) {
+void Block::populateSplitBlock(Block *b, Variable* const v, Variable* const u) {
 	b->addVariable(v);
 	for (Cit c=v->in.begin();c!=v->in.end();++c) {
 		if (canFollowLeft(*c,u))
@@ -352,7 +352,8 @@ void Block::populateSplitBlock(Block *b, Variable *v, Variable *u) {
  * with min lagrangrian multiplier and split at that point.
  * Returns the split constraint
  */
-Constraint* Block::splitBetween(Variable* vl, Variable* vr, Block* &lb, Block* &rb) {
+Constraint* Block::splitBetween(Variable* const vl, Variable* const vr,
+	       	Block* &lb, Block* &rb) {
 #ifdef RECTANGLE_OVERLAP_LOGGING
 	ofstream f(LOGFILE,ios::app);
 	f<<"  need to split between: "<<*vl<<" and "<<*vr<<endl;
@@ -384,16 +385,16 @@ void Block::split(Block* &l, Block* &r, Constraint* c) {
  */
 double Block::cost() {
 	double c = 0;
-	for (vector<Variable*>::iterator v=vars->begin();v!=vars->end();++v) {
+	for (Vit v=vars->begin();v!=vars->end();++v) {
 		double diff = (*v)->position() - (*v)->desiredPosition;
 		c += (*v)->weight * diff * diff;
 	}
 	return c;
 }
-ostream& operator <<(ostream &os, const Block &b)
+ostream& operator <<(ostream &os, const Block& b)
 {
 	os<<"Block:";
-	for(vector<Variable*>::iterator v=b.vars->begin();v!=b.vars->end();++v) {
+	for(Block::Vit v=b.vars->begin();v!=b.vars->end();++v) {
 		os<<" "<<**v;
 	}
 	if(b.deleted) {
diff --git a/src/removeoverlap/block.h b/src/libvpsc/block.h
index 997009a55..81e6c7637 100644
--- a/src/removeoverlap/block.h
+++ b/src/libvpsc/block.h
@@ -23,16 +23,20 @@ class StupidPriorityQueue;
 
 class Block
 {
+	typedef std::vector<Variable*> Variables;
+	typedef std::vector<Constraint*>::iterator Cit;
+	typedef std::vector<Variable*>::iterator Vit;
+
 	friend std::ostream& operator <<(std::ostream &os,const Block &b);
 public:
-	std::vector<Variable*> *vars;
+	Variables *vars;
 	double posn;
 	double weight;
 	double wposn;
-	Block(Variable *v=NULL);
+	Block(Variable* const v=NULL);
 	~Block(void);
 	Constraint* findMinLM();
-	Constraint* findMinLMBetween(Variable* lv, Variable* rv);
+	Constraint* findMinLMBetween(Variable* const lv, Variable* const rv);
 	Constraint* findMinInConstraint();
 	Constraint* findMinOutConstraint();
 	void deleteMinInConstraint();
@@ -54,14 +58,16 @@ public:
 private:
 	typedef enum {NONE, LEFT, RIGHT} Direction;
 	typedef std::pair<double, Constraint*> Pair;
-	void reset_active_lm(Variable *v, Variable *u);
-	double compute_dfdv(Variable *v, Variable *u, Constraint *&min_lm);
+	void reset_active_lm(Variable* const v, Variable* const u);
+	double compute_dfdv(Variable* const v, Variable* const u,
+		       	Constraint *&min_lm);
 	Pair compute_dfdv_between(
-			Variable*, Variable*, Variable*, Direction, bool);
-	bool canFollowLeft(Constraint *c, Variable *last);
-	bool canFollowRight(Constraint *c, Variable *last);
-	void populateSplitBlock(Block *b, Variable *v, Variable *u);
-	void addVariable(Variable *v);
+			Variable*, Variable* const, Variable* const,
+		       	const Direction, bool);
+	bool canFollowLeft(Constraint *c, const Variable* const last);
+	bool canFollowRight(Constraint *c, const Variable* const last);
+	void populateSplitBlock(Block *b, Variable* const v, Variable* const u);
+	void addVariable(Variable* const v);
 	void setUpConstraintHeap(PairingHeap<Constraint*>* &h,bool in);
 };
 
diff --git a/src/removeoverlap/blocks.cpp b/src/libvpsc/blocks.cpp
index 62b7e99f5..48f0237bf 100644
--- a/src/removeoverlap/blocks.cpp
+++ b/src/libvpsc/blocks.cpp
@@ -30,7 +30,7 @@ using std::copy;
 
 long blockTimeCtr;
 
-Blocks::Blocks(const int n, Variable *vs[]) : vs(vs),nvs(n) {
+Blocks::Blocks(const int n, Variable* const vs[]) : vs(vs),nvs(n) {
 	blockTimeCtr=0;
 	for(int i=0;i<nvs;i++) {
 		insert(new Block(vs[i]));
diff --git a/src/removeoverlap/blocks.h b/src/libvpsc/blocks.h
index 1a603eb41..0be1d7636 100644
--- a/src/removeoverlap/blocks.h
+++ b/src/libvpsc/blocks.h
@@ -34,7 +34,7 @@ class Constraint;
 class Blocks : public std::set<Block*>
 {
 public:
-	Blocks(const int n, Variable *vs[]);
+	Blocks(const int n, Variable* const vs[]);
 	~Blocks(void);
 	void mergeLeft(Block *r);
 	void mergeRight(Block *l);
@@ -45,7 +45,7 @@ public:
 private:
 	void dfsVisit(Variable *v, std::list<Variable*> *order);
 	void removeBlock(Block *doomed);
-	Variable **vs;
+	Variable* const *vs;
 	int nvs;
 };
 
diff --git a/src/removeoverlap/constraint.cpp b/src/libvpsc/constraint.cpp
index 7c200878b..7c200878b 100644
--- a/src/removeoverlap/constraint.cpp
+++ b/src/libvpsc/constraint.cpp
diff --git a/src/removeoverlap/constraint.h b/src/libvpsc/constraint.h
index 3da7449cd..3da7449cd 100644
--- a/src/removeoverlap/constraint.h
+++ b/src/libvpsc/constraint.h
diff --git a/src/libvpsc/csolve_VPSC.cpp b/src/libvpsc/csolve_VPSC.cpp
new file mode 100644
index 000000000..b78b01054
--- /dev/null
+++ b/src/libvpsc/csolve_VPSC.cpp
@@ -0,0 +1,124 @@
+/**
+ * \brief Bridge for C programs to access solve_VPSC (which is in C++)
+ *
+ * Authors:
+ *   Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * Released under GNU LGPL.  Read the file 'COPYING' for more information.
+ */
+#include <iostream>
+#include <cassert>
+#include "variable.h"
+#include "constraint.h"
+#include "generate-constraints.h"
+#include "solve_VPSC.h"
+#include "csolve_VPSC.h"
+extern "C" {
+Variable* newVariable(int id, double desiredPos, double weight) {
+	return new Variable(id,desiredPos,weight);
+}
+Constraint* newConstraint(Variable* left, Variable* right, double gap) {
+	return new Constraint(left,right,gap);
+}
+VPSC* newVPSC(int n, Variable* vs[], int m, Constraint* cs[]) {
+	return new VPSC(n,vs,m,cs);
+}
+VPSC* newIncVPSC(int n, Variable* vs[], int m, Constraint* cs[]) {
+	return (VPSC*)new IncVPSC(n,vs,m,cs);
+}
+
+int genXConstraints(int n, boxf* bb, Variable** vs, Constraint*** cs,int transitiveClosure) {
+	Rectangle* rs[n];
+	for(int i=0;i<n;i++) {
+		rs[i]=new Rectangle(bb[i].LL.x,bb[i].UR.x,bb[i].LL.y,bb[i].UR.y);
+	}
+	int m = generateXConstraints(n,rs,vs,*cs,transitiveClosure);
+	for(int i=0;i<n;i++) {
+		delete rs[i];
+	}
+	return m;
+}
+int genYConstraints(int n, boxf* bb, Variable** vs, Constraint*** cs) {
+	Rectangle* rs[n];
+	for(int i=0;i<n;i++) {
+		rs[i]=new Rectangle(bb[i].LL.x,bb[i].UR.x,bb[i].LL.y,bb[i].UR.y);
+	}
+	int m = generateYConstraints(n,rs,vs,*cs);
+	for(int i=0;i<n;i++) {
+		delete rs[i];
+	}
+	return m;
+}
+
+Constraint** newConstraints(int m) {
+	return new Constraint*[m];
+}
+void deleteConstraints(int m, Constraint **cs) {
+	for(int i=0;i<m;i++) {
+		delete cs[i];
+	}
+	delete [] cs;
+}
+void deleteConstraint(Constraint* c) {
+	delete c;
+}
+void deleteVariable(Variable* v) {
+	delete v;
+}
+void satisfyVPSC(VPSC* vpsc) {
+	try {
+		vpsc->satisfy();
+	} catch(const char *e) {
+		std::cerr << e << std::endl;
+		exit(1);
+	}
+}
+int getSplitCnt(IncVPSC *vpsc) {
+	return vpsc->splitCnt;
+}
+void deleteVPSC(VPSC *vpsc) {
+	assert(vpsc!=NULL);
+	delete vpsc;
+}
+void solveVPSC(VPSC* vpsc) {
+	vpsc->solve();
+}
+void splitIncVPSC(IncVPSC* vpsc) {
+	vpsc->splitBlocks();
+}
+void setVariableDesiredPos(Variable *v, double desiredPos) {
+	v->desiredPosition = desiredPos;
+}
+double getVariablePos(Variable *v) {
+	return v->position();
+}
+void remapInConstraints(Variable *u, Variable *v, double dgap) {
+	for(Constraints::iterator i=u->in.begin();i!=u->in.end();i++) {
+		Constraint* c=*i;	
+		c->right=v;
+		c->gap+=dgap;
+		v->in.push_back(c);
+	}
+	u->in.clear();
+}
+void remapOutConstraints(Variable *u, Variable *v, double dgap) {
+	for(Constraints::iterator i=u->out.begin();i!=u->out.end();i++) {
+		Constraint* c=*i;	
+		c->left=v;
+		c->gap+=dgap;
+		v->out.push_back(c);
+	}
+	u->out.clear();
+}
+int getLeftVarID(Constraint *c) {
+	return c->left->id;
+}
+int getRightVarID(Constraint *c){
+	return c->right->id;
+}
+double getSeparation(Constraint *c){
+	return c->gap;
+}
+}
diff --git a/src/libvpsc/csolve_VPSC.h b/src/libvpsc/csolve_VPSC.h
new file mode 100644
index 000000000..cd879effe
--- /dev/null
+++ b/src/libvpsc/csolve_VPSC.h
@@ -0,0 +1,54 @@
+/**
+ * \brief Bridge for C programs to access solve_VPSC (which is in C++)
+ *
+ * Authors:
+ *   Tim Dwyer <tgdwyer@gmail.com>
+ *
+ * Copyright (C) 2005 Authors
+ *
+ * Released under GNU LGPL.  Read the file 'COPYING' for more information.
+ */
+#ifndef _CSOLVE_VPSC_H_
+#define _CSOLVE_VPSC_H_
+#ifdef __cplusplus
+extern "C" {
+#endif
+typedef struct Variable Variable;
+Variable* newVariable(int id, double desiredPos, double weight);
+void setVariableDesiredPos(Variable *, double desiredPos);
+double getVariablePos(Variable*);
+
+typedef struct Constraint Constraint;
+Constraint* newConstraint(Variable* left, Variable* right, double gap);
+
+typedef struct VPSC VPSC;
+VPSC* newVPSC(int n, Variable* vs[], int m, Constraint* cs[]);
+void deleteVPSC(VPSC*);
+void deleteConstraint(Constraint*);
+void deleteVariable(Variable*);
+Constraint** newConstraints(int m);
+void deleteConstraints(int m,Constraint**);
+void remapInConstraints(Variable *u, Variable *v, double dgap);
+void remapOutConstraints(Variable *u, Variable *v, double dgap);
+int getLeftVarID(Constraint *c);
+int getRightVarID(Constraint *c);
+double getSeparation(Constraint *c);
+
+#ifndef HAVE_POINTF_S
+typedef struct pointf_s { double x, y; } pointf;
+typedef struct { pointf LL, UR; } boxf;
+#endif
+int genXConstraints(int n, boxf[], Variable** vs, Constraint*** cs,
+		int transitiveClosure);
+int genYConstraints(int n, boxf[], Variable** vs, Constraint*** cs);
+
+void satisfyVPSC(VPSC*);
+void solveVPSC(VPSC*);
+typedef struct IncVPSC IncVPSC;
+VPSC* newIncVPSC(int n, Variable* vs[], int m, Constraint* cs[]);
+void splitIncVPSC(IncVPSC*);
+int getSplitCnt(IncVPSC *vpsc);
+#ifdef __cplusplus
+}
+#endif
+#endif /* _CSOLVE_VPSC_H_ */
diff --git a/src/removeoverlap/generate-constraints.cpp b/src/libvpsc/generate-constraints.cpp
index 312ad960b..312ad960b 100644
--- a/src/removeoverlap/generate-constraints.cpp
+++ b/src/libvpsc/generate-constraints.cpp
diff --git a/src/removeoverlap/generate-constraints.h b/src/libvpsc/generate-constraints.h
index 56ee9536a..56ee9536a 100644
--- a/src/removeoverlap/generate-constraints.h
+++ b/src/libvpsc/generate-constraints.h
diff --git a/src/libvpsc/isnan.h b/src/libvpsc/isnan.h
new file mode 100644
index 000000000..388f9144b
--- /dev/null
+++ b/src/libvpsc/isnan.h
@@ -0,0 +1,57 @@
+#ifndef __ISNAN_H__
+#define __ISNAN_H__
+
+/*
+ * Temporary fix for various misdefinitions of isnan().
+ * isnan() is becoming undef'd in some .h files. 
+ * #include this last in your .cpp file to get it right.
+ *
+ * The problem is that isnan and isfinite are part of C99 but aren't part of
+ * the C++ standard (which predates C99).
+ *
+ * Authors:
+ *   Inkscape groupies and obsessive-compulsives
+ *
+ * Copyright (C) 2004 authors
+ *
+ * Released under GNU LGPL, read the file 'COPYING' for more information
+ *
+ * 2005 modification hereby placed in public domain.  Probably supercedes the 2004 copyright
+ * for the code itself.
+ */
+
+#include <math.h>
+/* You might try changing the above to <cmath> if you have problems.
+ * Whether you use math.h or cmath, you may need to edit the .cpp file
+ * and/or other .h files to use the same header file.
+ */
+
+#if defined(__isnan)
+# define isNaN(_a) (__isnan(_a))	/* MacOSX/Darwin definition < 10.4 */
+#elif defined(WIN32) || defined(_isnan)
+# define isNaN(_a) (_isnan(_a)) 	/* Win32 definition */
+#elif defined(isnan) || defined(__FreeBSD__)
+# define isNaN(_a) (isnan(_a))		/* GNU definition */
+#else
+# define isNaN(_a) (std::isnan(_a))
+#endif
+/* If the above doesn't work, then try (a != a).
+ * Also, please report a bug as per http://www.inkscape.org/report_bugs.php,
+ * giving information about what platform and compiler version you're using.
+ */
+
+
+#if defined(__isfinite)
+# define isFinite(_a) (__isfinite(_a))	/* MacOSX/Darwin definition < 10.4 */
+#elif defined(isfinite)
+# define isFinite(_a) (isfinite(_a))
+#else
+# define isFinite(_a) (std::isfinite(_a))
+#endif
+/* If the above doesn't work, then try (finite(_a) && !isNaN(_a)) or (!isNaN((_a) - (_a))).
+ * Also, please report a bug as per http://www.inkscape.org/report_bugs.php,
+ * giving information about what platform and compiler version you're using.
+ */
+
+
+#endif /* __ISNAN_H__ */
diff --git a/src/libvpsc/pairingheap/.dirstamp b/src/libvpsc/pairingheap/.dirstamp
new file mode 100644
index 000000000..e69de29bb
--- /dev/null
+++ b/src/libvpsc/pairingheap/.dirstamp
diff --git a/src/removeoverlap/pairingheap/PairingHeap.cpp b/src/libvpsc/pairingheap/PairingHeap.cpp
index 42d009c6a..202980b70 100644
--- a/src/removeoverlap/pairingheap/PairingHeap.cpp
+++ b/src/libvpsc/pairingheap/PairingHeap.cpp
@@ -2,7 +2,7 @@
  * \brief Pairing heap datastructure implementation
  *
  * Based on example code in "Data structures and Algorithm Analysis in C++"
- * by Mark Allen Weiss, used and released under the GPL by permission
+ * by Mark Allen Weiss, used and released under the LGPL by permission
  * of the author.
  *
  * No promises about correctness.  Use at your own risk!
@@ -13,7 +13,7 @@
  *
  * Copyright (C) 2005 Authors
  *
- * Released under GNU GPL.  Read the file 'COPYING' for more information.
+ * Released under GNU LGPL.  Read the file 'COPYING' for more information.
  */
 
 #include <vector>
@@ -180,7 +180,7 @@ template <class T>
 void PairingHeap<T>::decreaseKey( PairNode<T> *p,
 										  const T & newVal )
 {
-	if( p->element < newVal )
+	if( lessThan(p->element,newVal) )
 		return;    // newVal cannot be bigger
 	p->element = newVal;
 	if( p != root )
diff --git a/src/removeoverlap/pairingheap/PairingHeap.h b/src/libvpsc/pairingheap/PairingHeap.h
index 52941873b..63c9badcf 100644
--- a/src/removeoverlap/pairingheap/PairingHeap.h
+++ b/src/libvpsc/pairingheap/PairingHeap.h
@@ -2,7 +2,7 @@
  * \brief Pairing heap datastructure implementation
  *
  * Based on example code in "Data structures and Algorithm Analysis in C++"
- * by Mark Allen Weiss, used and released under the GPL by permission
+ * by Mark Allen Weiss, used and released under the LGPL by permission
  * of the author.
  *
  * No promises about correctness.  Use at your own risk!
@@ -13,7 +13,7 @@
  *
  * Copyright (C) 2005 Authors
  *
- * Released under GNU GPL.  Read the file 'COPYING' for more information.
+ * Released under GNU LGPL.  Read the file 'COPYING' for more information.
  */
 #ifndef PAIRING_HEAP_H_
 #define PAIRING_HEAP_H_
@@ -83,6 +83,11 @@ public:
 	PairNode<T> *insert( const T & x );
 	const T & findMin( ) const;
 	void deleteMin( );
+	const T extractMin( ) {
+		T v = findMin();
+		deleteMin();
+		return v;
+	}
 	void makeEmpty( );
 	void decreaseKey( PairNode<T> *p, const T & newVal );
 	void merge( PairingHeap<T> *rhs )
diff --git a/src/removeoverlap/pairingheap/dsexceptions.h b/src/libvpsc/pairingheap/dsexceptions.h
index bef2c78c5..bef2c78c5 100644
--- a/src/removeoverlap/pairingheap/dsexceptions.h
+++ b/src/libvpsc/pairingheap/dsexceptions.h
diff --git a/src/removeoverlap/placement_SolveVPSC.h b/src/libvpsc/placement_SolveVPSC.h
index 9f1c10cda..9f1c10cda 100755..100644
--- a/src/removeoverlap/placement_SolveVPSC.h
+++ b/src/libvpsc/placement_SolveVPSC.h
diff --git a/src/removeoverlap/remove_rectangle_overlap.cpp b/src/libvpsc/remove_rectangle_overlap.cpp
index 9fbef647b..6f6ace03a 100755..100644
--- a/src/removeoverlap/remove_rectangle_overlap.cpp
+++ b/src/libvpsc/remove_rectangle_overlap.cpp
@@ -40,18 +40,19 @@ double Rectangle::yBorder=0;
  *    too much in the first pass.
  */
 void removeRectangleOverlap(unsigned n, Rectangle *rs[], double xBorder, double yBorder) {
+	assert(0 <= n);
 	try {
 	// The extra gap avoids numerical imprecision problems
 	Rectangle::setXBorder(xBorder+EXTRA_GAP);
 	Rectangle::setYBorder(yBorder+EXTRA_GAP);
 	Variable **vs=new Variable*[n];
-	for(unsigned int i=0;i<n;i++) {
+	for(int i=0;i<n;i++) {
 		vs[i]=new Variable(i,0,1);
 	}
 	Constraint **cs;
 	double *oldX = new double[n];
 	int m=generateXConstraints(n,rs,vs,cs,true);
-	for(unsigned int i=0;i<n;i++) {
+	for(int i=0;i<n;i++) {
 		oldX[i]=vs[i]->desiredPosition;
 	}
 	VPSC vpsc_x(n,vs,m,cs);
@@ -61,7 +62,7 @@ void removeRectangleOverlap(unsigned n, Rectangle *rs[], double xBorder, double
 	f.close();
 #endif
 	vpsc_x.solve();
-	for(unsigned int i=0;i<n;i++) {
+	for(int i=0;i<n;i++) {
 		rs[i]->moveCentreX(vs[i]->position());
 	}
 	for(int i = 0; i < m; ++i) {
@@ -79,7 +80,7 @@ void removeRectangleOverlap(unsigned n, Rectangle *rs[], double xBorder, double
 	f.close();
 #endif
 	vpsc_y.solve();
-	for(unsigned int i=0;i<n;i++) {
+	for(int i=0;i<n;i++) {
 		rs[i]->moveCentreY(vs[i]->position());
 		rs[i]->moveCentreX(oldX[i]);
 	}
@@ -101,14 +102,14 @@ void removeRectangleOverlap(unsigned n, Rectangle *rs[], double xBorder, double
 		delete cs[i];
 	}
 	delete [] cs;
-	for(unsigned int i=0;i<n;i++) {
+	for(int i=0;i<n;i++) {
 		rs[i]->moveCentreX(vs[i]->position());
 		delete vs[i];
 	}
 	delete [] vs;
 	} catch (char const *str) {
 		std::cerr<<str<<std::endl;
-		for(unsigned int i=0;i<n;i++) {
+		for(int i=0;i<n;i++) {
 			std::cerr << *rs[i]<<std::endl;
 		}
 	}
diff --git a/src/removeoverlap/remove_rectangle_overlap.h b/src/libvpsc/remove_rectangle_overlap.h
index 08b035e31..08b035e31 100755..100644
--- a/src/removeoverlap/remove_rectangle_overlap.h
+++ b/src/libvpsc/remove_rectangle_overlap.h
diff --git a/src/removeoverlap/solve_VPSC.cpp b/src/libvpsc/solve_VPSC.cpp
index 77279c8a8..44918951d 100644
--- a/src/removeoverlap/solve_VPSC.cpp
+++ b/src/libvpsc/solve_VPSC.cpp
@@ -28,14 +28,16 @@ using std::ostringstream;
 using std::list;
 using std::set;
 
-IncVPSC::IncVPSC(const unsigned n, Variable *vs[], const unsigned m, Constraint *cs[]) 
+static const double ZERO_UPPERBOUND=-0.0000001;
+
+IncVPSC::IncVPSC(const unsigned n, Variable* const vs[], const unsigned m, Constraint *cs[]) 
 	: VPSC(n,vs,m,cs) {
 	inactive.assign(cs,cs+m);
 	for(ConstraintList::iterator i=inactive.begin();i!=inactive.end();++i) {
 		(*i)->active=false;
 	}
 }
-VPSC::VPSC(const unsigned n, Variable *vs[], const unsigned m, Constraint *cs[]) : cs(cs), m(m), vs(vs) {
+VPSC::VPSC(const unsigned n, Variable* const vs[], const unsigned m, Constraint *cs[]) : m(m), cs(cs), n(n), vs(vs) {
 	bs=new Blocks(n, vs);
 #ifdef RECTANGLE_OVERLAP_LOGGING
 	printBlocks();
@@ -79,7 +81,7 @@ void VPSC::satisfy() {
 	}
 	bs->cleanup();
 	for(unsigned i=0;i<m;i++) {
-		if(cs[i]->slack()<-0.0000001) {
+		if(cs[i]->slack() < ZERO_UPPERBOUND) {
 #ifdef RECTANGLE_OVERLAP_LOGGING
 			ofstream f(LOGFILE,ios::app);
 			f<<"Error: Unsatisfied constraint: "<<*cs[i]<<endl;
@@ -96,7 +98,7 @@ void VPSC::refine() {
 	// Solve shouldn't loop indefinately
 	// ... but just to make sure we limit the number of iterations
 	unsigned maxtries=100;
-	while(!solved&&maxtries>0) {
+	while(!solved&&maxtries>=0) {
 		solved=true;
 		maxtries--;
 		for(set<Block*>::const_iterator i=bs->begin();i!=bs->end();++i) {
@@ -123,8 +125,8 @@ void VPSC::refine() {
 		}
 	}
 	for(unsigned i=0;i<m;i++) {
-		if(cs[i]->slack()<-0.0000001) {
-			assert(cs[i]->slack()>-0.0000001);
+		if(cs[i]->slack() < ZERO_UPPERBOUND) {
+			assert(cs[i]->slack()>ZERO_UPPERBOUND);
 			throw "Unsatisfied constraint";
 		}
 	}
@@ -177,7 +179,7 @@ void IncVPSC::satisfy() {
 	splitBlocks();
 	long splitCtr = 0;
 	Constraint* v = NULL;
-	while((v=mostViolated(inactive))&&(v->equality || v->slack()<-0.000001)) {
+	while((v=mostViolated(inactive))&&(v->equality || v->slack() < ZERO_UPPERBOUND)) {
 		assert(!v->active);
 		Block *lb = v->left->block, *rb = v->right->block;
 		if(lb != rb) {
@@ -198,10 +200,13 @@ void IncVPSC::satisfy() {
 	bs->cleanup();
 	for(unsigned i=0;i<m;i++) {
 		v=cs[i];
-		if(v->slack()<-0.0000001) {
-			//assert(cs[i]->slack()>-0.0000001);
+		if(v->slack() < ZERO_UPPERBOUND) {
 			ostringstream s;
 			s<<"Unsatisfied constraint: "<<*v;
+#ifdef RECTANGLE_OVERLAP_LOGGING
+			ofstream f(LOGFILE,ios::app);
+			f<<s.str()<<endl;
+#endif
 			throw s.str().c_str();
 		}
 	}
@@ -234,7 +239,7 @@ void IncVPSC::splitBlocks() {
 	for(set<Block*>::const_iterator i(bs->begin());i!=bs->end();++i) {
 		Block* b = *i;
 		Constraint* v=b->findMinLM();
-		if(v!=NULL && v->lm < -0.0000001) {
+		if(v!=NULL && v->lm < ZERO_UPPERBOUND) {
 			assert(!v->equality);
 #ifdef RECTANGLE_OVERLAP_LOGGING
 			f<<"    found split point: "<<*v<<" lm="<<v->lm<<endl;
@@ -289,7 +294,7 @@ Constraint* IncVPSC::mostViolated(ConstraintList &l) {
 	// move the last element over the deletePoint and resize
 	// downwards.  There is always at least 1 element in the
 	// vector because of search.
-	if(deletePoint != end && (minSlack<-0.0000001||v->equality)) {
+	if(deletePoint != end && (minSlack<ZERO_UPPERBOUND||v->equality)) {
 		*deletePoint = l[l.size()-1];
 		l.resize(l.size()-1);
 	}
diff --git a/src/removeoverlap/solve_VPSC.h b/src/libvpsc/solve_VPSC.h
index 9f6244a5a..4cd5559d6 100644
--- a/src/removeoverlap/solve_VPSC.h
+++ b/src/libvpsc/solve_VPSC.h
@@ -9,6 +9,14 @@
  *
  * Released under GNU LGPL.  Read the file 'COPYING' for more information.
  */
+
+//
+// TODO: Really, we should have three classes: VPSC, IncrementalVPSC and
+// StaticVPSC, where the latter two inherit from VPSC.  StaticVPSC would be
+// the equivalent of what is currently VPSC.
+// Also, a lot of the code specific to one or other of these concrete
+// implementations should be moved from Block and Blocks: e.g. mergeLeft etc.
+//
 #ifndef SEEN_REMOVEOVERLAP_SOLVE_VPSC_H
 #define SEEN_REMOVEOVERLAP_SOLVE_VPSC_H
 
@@ -25,15 +33,16 @@ public:
 	virtual void satisfy();
 	virtual void solve();
 
-	VPSC(const unsigned n, Variable *vs[], const unsigned m, Constraint *cs[]);
+	VPSC(const unsigned n, Variable* const vs[], const unsigned m, Constraint *cs[]);
 	virtual ~VPSC();
-	Constraint** getConstraints() { return cs; }
-	Variable** getVariables() { return vs; }
+	Constraint** getConstraints(unsigned &m) { m=this->m; return cs; }
+	const Variable* const * getVariables(unsigned &n) { n=this->n; return vs; }
 protected:
 	Blocks *bs;
-	Constraint **cs;
 	unsigned m;
-	Variable **vs;
+	Constraint **cs;
+	unsigned n;
+	const Variable* const *vs;
 	void printBlocks();
 private:
 	void refine();
@@ -48,7 +57,7 @@ public:
 	void solve();
 	void moveBlocks();
 	void splitBlocks();
-	IncVPSC(const unsigned n, Variable *vs[], const unsigned m, Constraint *cs[]);
+	IncVPSC(const unsigned n, Variable* const vs[], const unsigned m, Constraint *cs[]);
 private:
 	typedef std::vector<Constraint*> ConstraintList;
 	ConstraintList inactive;
diff --git a/src/removeoverlap/variable.cpp b/src/libvpsc/variable.cpp
index 1890f788e..1890f788e 100644
--- a/src/removeoverlap/variable.cpp
+++ b/src/libvpsc/variable.cpp
diff --git a/src/removeoverlap/variable.h b/src/libvpsc/variable.h
index b1ab66c95..b1ab66c95 100644
--- a/src/removeoverlap/variable.h
+++ b/src/libvpsc/variable.h
diff --git a/src/removeoverlap/Makefile_insert b/src/removeoverlap/Makefile_insert
index e28304431..9df2ca2e3 100644
--- a/src/removeoverlap/Makefile_insert
+++ b/src/removeoverlap/Makefile_insert
@@ -6,26 +6,5 @@ removeoverlap/clean:
 	rm -f removeoverlap/libremoveoverlap.a $(removeoverlap_libremoveoverlap_a_OBJECTS)
 
 removeoverlap_libremoveoverlap_a_SOURCES =	\
-	removeoverlap/block.cpp	\
-	removeoverlap/block.h	\
-	removeoverlap/blocks.cpp	\
-	removeoverlap/blocks.h	\
-	removeoverlap/constraint.cpp	\
-	removeoverlap/constraint.h	\
-	removeoverlap/generate-constraints.cpp	\
-	removeoverlap/generate-constraints.h	\
-	removeoverlap/remove_rectangle_overlap.cpp	\
-	removeoverlap/remove_rectangle_overlap.h	\
 	removeoverlap/removeoverlap.cpp	\
-	removeoverlap/removeoverlap.h	\
-	removeoverlap/solve_VPSC.cpp	\
-	removeoverlap/solve_VPSC.h	\
-	removeoverlap/variable.cpp	\
-	removeoverlap/variable.h	\
-	removeoverlap/pairingheap/dsexceptions.h	\
-	removeoverlap/pairingheap/PairingHeap.cpp	\
-	removeoverlap/pairingheap/PairingHeap.h
-
-removeoverlap_remove_rectangle_overlap_test_SOURCES =	\
-	removeoverlap/remove_rectangle_overlap-test.cpp
-removeoverlap_remove_rectangle_overlap_test_LDADD = removeoverlap/libremoveoverlap.a -lglib-2.0
+	removeoverlap/removeoverlap.h
diff --git a/src/removeoverlap/pairingheap/.cvsignore b/src/removeoverlap/pairingheap/.cvsignore
deleted file mode 100644
index e8014d011..000000000
--- a/src/removeoverlap/pairingheap/.cvsignore
+++ /dev/null
@@ -1,5 +0,0 @@
-Makefile
-Makefile.in
-.deps
-makefile
-.dirstamp
diff --git a/src/removeoverlap/placement_SolveVPSC.cpp b/src/removeoverlap/placement_SolveVPSC.cpp
deleted file mode 100755
index a9f4344c8..000000000
--- a/src/removeoverlap/placement_SolveVPSC.cpp
+++ /dev/null
@@ -1,130 +0,0 @@
-#include <jni.h>
-#include "placement_SolveVPSC.h"
-#include <stdio.h>
-#include "solve_VPSC.h"
-#include "variable.h"
-#include "constraint.h"
-#include "remove_rectangle_overlap.h"
-#include "generate-constraints.h"
-#include <assert.h>
-#include <map>
-#define MaxSize 500
-
-JNIEXPORT jdouble JNICALL Java_placement_SolveVPSC_solve
-  (JNIEnv *env, jobject obj, jobjectArray vName, jdoubleArray vWeight, jdoubleArray vDesPos, jintArray cLeft, jintArray cRight, jdoubleArray cGap, jdoubleArray vResult, jint mode)
-{
-	jsize n = env->GetArrayLength(vWeight);
-	jsize m = env->GetArrayLength(cLeft);
-	int i;
-	double *lvWeight = env->GetDoubleArrayElements(vWeight, 0);
-	double *lvDesPos = env->GetDoubleArrayElements(vDesPos, 0);
-	long *lcLeft = env->GetIntArrayElements(cLeft, 0);
-	long *lcRight = env->GetIntArrayElements(cRight, 0);
-	double *lcGap = env->GetDoubleArrayElements(cGap, 0);
-	Variable **vs=new Variable*[n];
-	Constraint **cs=new Constraint*[m];
-	for (i=0; i<n; i++) {
-		jstring lvName = (jstring)env->GetObjectArrayElement(vName, i);
-		const char *name = env->GetStringUTFChars(lvName, NULL); 
-		// once upon a time variables had real names, now you'll have to 
-		// track them by number.
-		vs[i]=new Variable(i,lvDesPos[i],lvWeight[i]);
-	}
-	for (i=0; i<m; i++) {
-		cs[i]=new Constraint(vs[lcLeft[i]],vs[lcRight[i]],lcGap[i]);
-	}
-	double cost=0;
-	VPSC vpsc(vs,n,cs,m);
-	if(mode==0) {
-		vpsc.satisfy();
-	} else {
-		vpsc.solve();
-	}
-	for (i=0; i<n; i++) {
-		double p=vs[i]->position();
-		env->SetDoubleArrayRegion(vResult, i,1,&p);
-	}
-	for (i=0; i<m; i++) {
-		delete cs[i];
-	}
-	delete [] cs;
-	for (i=0; i<n; i++) {
-		delete vs[i];
-	}
-	env->ReleaseIntArrayElements(cLeft, lcLeft, 0);
-	env->ReleaseIntArrayElements(cRight, lcRight, 0);
-	env->ReleaseDoubleArrayElements(cGap, lcGap, 0);
-	env->ReleaseDoubleArrayElements(vWeight, lvWeight, 0);
-	env->ReleaseDoubleArrayElements(vDesPos, lvDesPos, 0);
-	delete [] vs;
-	return cost;
-}
-
-static Variable **vs;
-static Constraint **cs;
-static int m,n;
-JNIEXPORT jint JNICALL Java_placement_SolveVPSC_generateXConstraints
-(JNIEnv *env, jobject obj, jdoubleArray rMinX, jdoubleArray rMaxX, jdoubleArray rMinY, jdoubleArray rMaxY, jdoubleArray rWeight) {
-	n = (int)env->GetArrayLength(rWeight);
-	Rectangle **rs=new Rectangle*[n];
-	double *ws = env->GetDoubleArrayElements(rWeight, 0);
-	double *minX = env->GetDoubleArrayElements(rMinX, 0);
-	double *maxX = env->GetDoubleArrayElements(rMaxX, 0);
-	double *minY = env->GetDoubleArrayElements(rMinY, 0);
-	double *maxY = env->GetDoubleArrayElements(rMaxY, 0);
-	for(int i=0;i<n;i++) rs[i]=new Rectangle(minX[i],maxX[i],minY[i],maxY[i]);
-	m = generateXConstraints(rs, ws, n, vs, cs, true);
-	return m;
-}
-
-JNIEXPORT jint JNICALL Java_placement_SolveVPSC_generateYConstraints
-(JNIEnv *env, jobject obj, jdoubleArray rMinX, jdoubleArray rMaxX, jdoubleArray rMinY, jdoubleArray rMaxY, jdoubleArray rWeight) {
-	n = (int)env->GetArrayLength(rWeight);
-	Rectangle **rs=new Rectangle*[n];
-	double *ws = env->GetDoubleArrayElements(rWeight, 0);
-	double *minX = env->GetDoubleArrayElements(rMinX, 0);
-	double *maxX = env->GetDoubleArrayElements(rMaxX, 0);
-	double *minY = env->GetDoubleArrayElements(rMinY, 0);
-	double *maxY = env->GetDoubleArrayElements(rMaxY, 0);
-	for(int i=0;i<n;i++) rs[i]=new Rectangle(minX[i],maxX[i],minY[i],maxY[i]);
-	m = generateYConstraints(rs, ws, n, vs, cs);
-	return m;
-}
-using namespace std;
-JNIEXPORT void JNICALL Java_placement_SolveVPSC_getConstraints
-(JNIEnv *env, jobject obj, jintArray cLeft, jintArray cRight, jdoubleArray cGap) {
-	map<Variable*,int> vmap;
-	for(int i=0;i<n;i++) {
-		vmap[vs[i]]=i;
-	}
-	
-	for(int i=0;i<m;i++) {
-		jint l=vmap[cs[i]->left];
-		jint r=vmap[cs[i]->right];
-		double g=cs[i]->gap;
-		env->SetIntArrayRegion(cLeft, i,1,&l);
-		env->SetIntArrayRegion(cRight, i,1,&r);
-		env->SetDoubleArrayRegion(cGap, i,1,&g);
-	}
-}
-JNIEXPORT void JNICALL Java_placement_SolveVPSC_removeOverlaps
-(JNIEnv *env, jobject obj, jdoubleArray rMinX, jdoubleArray rMaxX, jdoubleArray rMinY, jdoubleArray rMaxY) {
-	//assert(1==2); //break for debugging
-	n = (int)env->GetArrayLength(rMinX);
-	Rectangle **rs=new Rectangle*[n];
-	double *minX = env->GetDoubleArrayElements(rMinX, 0);
-	double *maxX = env->GetDoubleArrayElements(rMaxX, 0);
-	double *minY = env->GetDoubleArrayElements(rMinY, 0);
-	double *maxY = env->GetDoubleArrayElements(rMaxY, 0);
-	for(int i=0;i<n;i++) rs[i]=new Rectangle(minX[i],maxX[i],minY[i],maxY[i]);
-	removeRectangleOverlap(rs,n,0,0);
-	for (i=0; i<n; i++) {
-		double x=rs[i]->getMinX();
-		double y=rs[i]->getMinY();
-		env->SetDoubleArrayRegion(rMinX, i,1,&x);
-		env->SetDoubleArrayRegion(rMinY, i,1,&y);
-	}
-	delete [] rs;
-	env->ReleaseDoubleArrayElements(rMaxX, maxX, 0);
-	env->ReleaseDoubleArrayElements(rMaxY, maxY, 0);
-}
\ No newline at end of file
diff --git a/src/removeoverlap/remove_rectangle_overlap-test.cpp b/src/removeoverlap/remove_rectangle_overlap-test.cpp
deleted file mode 100644
index 20f5489cc..000000000
--- a/src/removeoverlap/remove_rectangle_overlap-test.cpp
+++ /dev/null
@@ -1,308 +0,0 @@
-#include "removeoverlap/remove_rectangle_overlap.h"
-#include <unistd.h>  // for alarm()
-#include <time.h>  // for srand seed and clock().
-#include <glib/gmacros.h>
-#include <glib/gmem.h>
-#include <cstdlib>
-#include <cmath>
-#include "removeoverlap/generate-constraints.h"
-#include "utest/utest.h"
-using std::abs;
-using std::rand;
-
-static bool
-possibly_eq(double const a, double const b)
-{
-    return abs(a - b) < 1e-13;
-}
-
-static bool
-possibly_le(double const a, double const b)
-{
-    return a - b < 1e-13;
-}
-
-static void
-show_rects(unsigned const n_rects, double const rect2coords[][4])
-{
-    for (unsigned i = 0; i < n_rects; ++i) {
-        printf("{%g, %g, %g, %g},\n",
-               rect2coords[i][0],
-               rect2coords[i][1],
-               rect2coords[i][2],
-               rect2coords[i][3]);
-    }
-}
-
-/**
- * Returns the signum of x, but erring towards returning 0 if x is "not too far" from 0.  ("Not too
- * far from 0" means [-0.9, 0.9] in current version.)
- */
-static int
-sgn0(double const x)
-{
-    if (x <= -0.9) {
-        return -1;
-    } else if (0.9 <= x) {
-        return 1;
-    } else {
-        return 0;
-    }
-}
-
-static void
-test_case(unsigned const n_rects, double const rect2coords[][4])
-{
-    Rectangle **rs = (Rectangle **) g_malloc(sizeof(Rectangle*) * n_rects);
-    for (unsigned i = 0; i < n_rects; ++i) {
-        rs[i] = new Rectangle(rect2coords[i][0],
-                              rect2coords[i][1],
-                              rect2coords[i][2],
-                              rect2coords[i][3]);
-    }
-    removeRectangleOverlap(n_rects,rs,0.0, 0.0);
-    for (unsigned i = 0; i < n_rects; ++i) {
-        UTEST_ASSERT(possibly_eq(rs[i]->width(), (rect2coords[i][1] -
-                                                  rect2coords[i][0]  )));
-        UTEST_ASSERT(possibly_eq(rs[i]->height(), (rect2coords[i][3] -
-                                                   rect2coords[i][2]  )));
-        for (unsigned j = 0; j < i; ++j) {
-            if (!( possibly_le(rs[i]->getMaxX(), rs[j]->getMinX()) ||
-                   possibly_le(rs[j]->getMaxX(), rs[i]->getMinX()) ||
-                   possibly_le(rs[i]->getMaxY(), rs[j]->getMinY()) ||
-                   possibly_le(rs[j]->getMaxY(), rs[i]->getMinY())   )) {
-                show_rects(n_rects, rect2coords);
-                char buf[32];
-                sprintf(buf, "[%u],[%u] of %u", j, i, n_rects);
-                utest__fail("Found overlap among ", buf, " rectangles");
-            }
-        }
-
-        /* Optimality test. */
-        {
-            bool found_block[2] = {false, false};
-            int const desired_movement[2] = {sgn0(rect2coords[i][0] - rs[i]->getMinX()),
-                                             sgn0(rect2coords[i][2] - rs[i]->getMinY())};
-            for (unsigned j = 0; j < n_rects; ++j) {
-                if (j == i)
-                    continue;
-                for (unsigned d = 0; d < 2; ++d) {
-                    if ( ( desired_movement[d] < 0
-                           ? abs(rs[j]->getMaxD(d) - rs[i]->getMinD(d))
-                           : abs(rs[i]->getMaxD(d) - rs[j]->getMinD(d)) )
-                         < .002 ) {
-                        found_block[d] = true;
-                    }
-                }
-            }
-
-            for (unsigned d = 0; d < 2; ++d) {
-                if ( !found_block[d]
-                     && desired_movement[d] != 0 ) {
-                    show_rects(n_rects, rect2coords);
-                    char buf[32];
-                    sprintf(buf, "%c in rectangle [%u] of %u", "XY"[d], i, n_rects);
-                    utest__fail("Found clear non-optimality in ", buf, " rectangles");
-                }
-            }
-        }
-    }
-    for (unsigned i = 0; i < n_rects; ++i) {
-        delete rs[i];
-    }
-    g_free(rs);
-}
-
-int main()
-{
-    srand(time(NULL));
-
-    /* Ensure that the program doesn't run for more than 60 seconds. */
-    alarm(60);
-
-    utest_start("removeRectangleOverlap(zero gaps)");
-
-    /* Derived from Bulia's initial test case.  This used to crash. */
-    UTEST_TEST("eg0") {
-        double case0[][4] = {
-            {-180.5, 69.072, 368.071, 629.071},
-            {99.5, 297.644, 319.5, 449.071},
-            {199.5, 483.358, 450.929, 571.929},
-            {168.071, 277.644, 462.357, 623.357},
-            {99.5, 99.751, 479.5, 674.786},
-            {-111.929, 103.358, 453.786, 611.929},
-            {-29.0714, 143.358, 273.786, 557.643},
-            {122.357, 269.072, 322.357, 531.929},
-            {256.643, 357.644, 396.643, 520.5}
-        };
-        test_case(G_N_ELEMENTS(case0), case0);
-    }
-
-#if 0 /* This involves a zero-height rect, so we'll ignore for the moment. */
-    UTEST_TEST("eg1") {
-        double case1[][4] = {
-            {5, 14, 9, 14},
-            {6, 13, 6, 8},
-            {11, 12, 5, 5},
-            {5, 8, 5, 7},
-            {12, 14, 14, 15},
-            {12, 14, 1, 14},
-            {1, 15, 14, 15},
-            {5, 6, 13, 13}
-        };
-        test_case(G_N_ELEMENTS(case1), case1);
-    }
-#endif
-
-    /* The next few examples used to result in overlaps. */
-    UTEST_TEST("eg2") {
-        double case2[][4] = {
-            {3, 4, 6, 13},
-            {0, 1, 0, 5},
-            {0, 4, 1, 6},
-            {2, 5, 0, 6},
-            {0, 10, 9, 13},
-            {5, 11, 1, 13},
-            {1, 2, 3, 8}
-        };
-        test_case(G_N_ELEMENTS(case2), case2);
-    }
-
-    UTEST_TEST("eg3") {
-        double case3[][4] = {
-            {0, 5, 0, 3},
-            {1, 2, 1, 3},
-            {3, 7, 4, 7},
-            {0, 9, 4, 5},
-            {3, 7, 0, 3}
-        };
-        test_case(G_N_ELEMENTS(case3), case3);
-    }
-
-    UTEST_TEST("eg4") {
-        double case4[][4] = {
-            {0, 1, 2, 3},
-            {0, 4, 0, 4},
-            {1, 6, 0, 4},
-            {2, 3, 4, 5},
-            {0, 5, 4, 6}
-        };
-        test_case(G_N_ELEMENTS(case4), case4);
-    }
-
-    UTEST_TEST("eg5") {
-        double case5[][4] = {
-            {1, 5, 1, 2},
-            {1, 6, 5, 7},
-            {6, 8, 1, 2},
-            {2, 3, 1, 4},
-            {5, 8, 2, 6}
-        };
-        test_case(G_N_ELEMENTS(case5), case5);
-    }
-
-    /* This one causes overlap in 2005-12-19 04:00 UTC version. */
-    UTEST_TEST("olap6") {
-        double case6[][4] = {
-            {7, 22, 39, 54},
-            {7, 33, 0, 59},
-            {3, 26, 16, 56},
-            {7, 17, 18, 20},
-            {1, 59, 11, 26},
-            {19, 20, 13, 49},
-            {1, 10, 0, 4},
-            {47, 52, 1, 3}
-        };
-        test_case(G_N_ELEMENTS(case6), case6);
-    }
-
-    /* The next two examples caused loops in the version at 2005-12-07 04:00 UTC. */
-    UTEST_TEST("loop0") {
-        double loop0[][4] = {
-            {13, 16, 6, 27},
-            {0, 6, 0, 12},
-            {11, 14, 1, 10},
-            {12, 39, 5, 24},
-            {14, 34, 4, 7},
-            {1, 30, 20, 27},
-            {1, 6, 1, 2},
-            {19, 28, 10, 24},
-            {4, 34, 15, 21},
-            {7, 13, 13, 34}
-        };
-        test_case(G_N_ELEMENTS(loop0), loop0);
-    }
-
-    UTEST_TEST("loop1") {
-        double loop1[][4] = {
-            {6, 18, 9, 16},
-            {8, 26, 10, 13},
-            {3, 10, 0, 14},
-            {0, 5, 16, 22},
-            {1, 8, 11, 21},
-            {1, 5, 0, 13},
-            {24, 25, 0, 2}
-        };
-        test_case(G_N_ELEMENTS(loop1), loop1);
-    }
-
-    UTEST_TEST("loop2") {
-        double loop2[][4] = {
-            {16, 22, 9, 16},
-            {8, 9, 14, 19},
-            {17, 25, 8, 13},
-            {10, 26, 26, 29},
-            {14, 19, 9, 19},
-            {0, 18, 3, 12},
-            {7, 8, 14, 22},
-            {14, 20, 25, 29}
-        };
-        test_case(G_N_ELEMENTS(loop2), loop2);
-    }
-
-    /* Random cases of up to 10 rectangles, with small non-neg int coords. */
-    for (unsigned n = 0; n <= 10; ++n) {
-        char buf[64];
-        sprintf(buf, "random ints with %u rectangles", n);
-        UTEST_TEST(buf) {
-            unsigned const fld_size = 8 * n;
-            double (*coords)[4] = (double (*)[4]) g_malloc(n * 4 * sizeof(double));
-            clock_t const clock_stop = clock() + CLOCKS_PER_SEC;
-            for (unsigned repeat = (n == 0 ? 1
-                                    : n == 1 ? 36
-                                    : (1 << 16)  ); repeat--;) {
-                for (unsigned i = 0; i < n; ++i) {
-                    for (unsigned d = 0; d < 2; ++d) {
-                        //unsigned const start = rand() % fld_size;
-                        //unsigned const end = start + rand() % (fld_size - start);
-                        unsigned const end = 1 + (rand() % (fld_size - 1));
-                        unsigned const start = rand() % end;
-                        coords[i][2 * d] = start;
-                        coords[i][2 * d + 1] = end;
-                    }
-                }
-                test_case(n, coords);
-                if (clock() >= clock_stop) {
-                    break;
-                }
-            }
-            g_free(coords);
-        }
-    }
-
-    return ( utest_end()
-             ? EXIT_SUCCESS
-             : EXIT_FAILURE );
-}
-
-
-/*
-  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=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/removeoverlap/removeoverlap.cpp b/src/removeoverlap/removeoverlap.cpp
index d79fa9eab..3a8481db2 100644
--- a/src/removeoverlap/removeoverlap.cpp
+++ b/src/removeoverlap/removeoverlap.cpp
@@ -12,8 +12,8 @@
 #include "util/glib-list-iterators.h"
 #include "sp-item.h"
 #include "sp-item-transform.h"
-#include "removeoverlap/generate-constraints.h"
-#include "removeoverlap/remove_rectangle_overlap.h"
+#include "libvpsc/generate-constraints.h"
+#include "libvpsc/remove_rectangle_overlap.h"
 
 /**
 * Takes a list of inkscape items and moves them as little as possible
@@ -29,7 +29,7 @@ void removeoverlap(GSList const *const items, double const xGap, double const yG
 	std::list<SPItem *> selected;
 	selected.insert<GSListConstIterator<SPItem *> >(selected.end(), items, NULL);
 	if (selected.empty()) return;
-	unsigned n=selected.size();
+	int n=selected.size();
 
 	//Check 2 or more selected objects
 	if (n < 2) return;