From 12b21e1d27f43deaa748419919b40b80cedd0ddd Mon Sep 17 00:00:00 2001
From: Tim Dwyer <tgdwyer@gmail.com>
Date: Wed, 12 Jul 2006 00:55:58 +0000
Subject: 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)
---
 src/libcola/gradient_projection.cpp | 234 ++++++++++++++++++++++++++++++++++++
 1 file changed, 234 insertions(+)
 create mode 100644 src/libcola/gradient_projection.cpp

(limited to 'src/libcola/gradient_projection.cpp')

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 :
-- 
cgit v1.2.3