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)

1 files changed, 242 insertions, 0 deletions

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