Merge https://gitlab.com/inkscape/inkscape into selectable-knots

1 files changed, 430 insertions, 184 deletions

diff --git a/src/libcola/gradient_projection.cpp b/src/libcola/gradient_projection.cpp
index c079c3149..50b8d27d7 100644
--- a/src/libcola/gradient_projection.cpp
+++ b/src/libcola/gradient_projection.cpp
@@ -1,165 +1,372 @@
+/*
+ * vim: ts=4 sw=4 et tw=0 wm=0
+ *
+ * libcola - A library providing force-directed network layout using the 
+ *           stress-majorization method subject to separation constraints.
+ *
+ * Copyright (C) 2006-2008  Monash University
+ *
+ * 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.
+ * See the file LICENSE.LGPL distributed with the library.
+ *
+ * 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.
+ *
+*/
+
 /**********************************************************
  *
  * 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 <iostream>
+#include <cmath>
+#include <ctime>
+
 #include <libvpsc/solve_VPSC.h>
 #include <libvpsc/variable.h>
+#include <libvpsc/isnan.h>
 #include <libvpsc/constraint.h>
-#include "gradient_projection.h"
-#include <iostream>
-#include <2geom/math-utils.h>
-#include <math.h>
+#include <libvpsc/assertions.h>
+
+#include "libcola/cluster.h"
+#include "libcola/gradient_projection.h"
+#include "libcola/straightener.h"
+#include "libcola/commondefs.h"
+//#define CHECK_CONVERGENCE_BY_COST 1
+
 
 using namespace std;
 using namespace vpsc;
-//#define CONMAJ_LOGGING 1
+namespace cola {
+GradientProjection::GradientProjection(
+    const Dim k,
+    std::valarray<double> *denseQ,
+    const double tol,
+    const unsigned max_iterations,
+    CompoundConstraints const *ccs,
+    UnsatisfiableConstraintInfos *unsatisfiableConstraints,
+    NonOverlapConstraintsMode nonOverlapConstraints,
+    RootCluster* clusterHierarchy,
+    vpsc::Rectangles* rs,
+    const bool scaling,
+    SolveWithMosek solveWithMosek) 
+        : k(k), 
+          denseSize(static_cast<unsigned>((floor(sqrt(static_cast<double>(denseQ->size())))))),
+          denseQ(denseQ), 
+          rs(rs),
+          ccs(ccs),
+          unsatisfiableConstraints(unsatisfiableConstraints),
+          nonOverlapConstraints(nonOverlapConstraints),
+          clusterHierarchy(clusterHierarchy),
+          tolerance(tol), 
+          max_iterations(max_iterations),
+          sparseQ(NULL),
+          solveWithMosek(solveWithMosek),
+          scaling(scaling)
+{
+    printf("GP Instance: scaling=%d, mosek=%d\n",scaling,solveWithMosek);
+    for(unsigned i=0;i<denseSize;i++) {
+        vars.push_back(new vpsc::Variable(i,1,1));
+    }
+    if(scaling) {
+        scaledDenseQ.resize(denseSize*denseSize);
+        for(unsigned i=0;i<denseSize;i++) {
+            vars[i]->scale=1./sqrt(fabs((*denseQ)[i*denseSize+i]));
+            // XXX: Scale can sometimes be set to infinity here when 
+            //      there are nodes not connected to any other node.
+            //      Thus we just set the scale for such a variable to 1.
+            if (!isFinite(vars[i]->scale))
+            {
+                vars[i]->scale = 1;
+            }
+        }
+        // the following computes S'QS for Q=denseQ
+        // and S is diagonal matrix of scale factors
+        for(unsigned i=0;i<denseSize;i++) {
+            for(unsigned j=0;j<denseSize;j++) {
+                scaledDenseQ[i*denseSize+j]=(*denseQ)[i*denseSize+j]*vars[i]->scale
+                    *vars[j]->scale;
+            }
+        }
+        this->denseQ = &scaledDenseQ;
+    }
+    //dumpSquareMatrix(*this->denseQ);
+    //dumpSquareMatrix(scaledDenseQ);
+
+    if(ccs) {
+        for(CompoundConstraints::const_iterator c=ccs->begin();
+                c!=ccs->end();++c) {
+            (*c)->generateVariables(k, vars);
+            OrthogonalEdgeConstraint* e=dynamic_cast<OrthogonalEdgeConstraint*>(*c);
+            if(e) {
+                orthogonalEdges.push_back(e);
+            }
+        }
+        for(CompoundConstraints::const_iterator c=ccs->begin();
+                c!=ccs->end();++c) {
+            (*c)->generateSeparationConstraints(k, vars, gcs, *rs);
+        }
+    }
+    /*
+    if(clusterHierarchy) {
+        clusterHierarchy->createVars(k,*rs,vars);
+    }
+    */
+    numStaticVars=vars.size();
+    //solver=setupVPSC();
+}
+static inline double dotProd(valarray<double> const & a, valarray<double> const & b) {
+    double p = 0;
+    for (unsigned i=0; i<a.size(); i++) {
+        p += a[i]*b[i];
+    }
+    return p;
+}
+double GradientProjection::computeCost(
+        valarray<double> const &b,
+        valarray<double> const &x) const {
+    // computes cost = 2 b x - x A x
+    double cost = 2. * dotProd(b,x);
+    valarray<double> Ax(x.size());
+    for (unsigned i=0; i<denseSize; i++) {
+        Ax[i] = 0;
+        for (unsigned j=0; j<denseSize; j++) {
+            Ax[i] += (*denseQ)[i*denseSize+j]*x[j];
+        }
+    }
+    if(sparseQ) {
+        valarray<double> r(x.size());
+        sparseQ->rightMultiply(x,r);
+        Ax+=r;
+    }
+    return cost - dotProd(x,Ax);
+}
 
-static void dumpVPSCException(char const *str, IncSolver* solver) {
-    cerr<<str<<endl;
-    unsigned m;
-    Constraint** cs = solver->getConstraints(m);
-    for(unsigned i=0;i<m;i++) {
-        cerr << *cs[i] << endl;
+double GradientProjection::computeSteepestDescentVector(
+        valarray<double> const &b,
+        valarray<double> const &x,
+        valarray<double> &g) const {
+    // find steepest descent direction
+    //  g = 2 ( b - A x )
+    //    where: A = denseQ + sparseQ
+    //  g = 2 ( b - denseQ x) - 2 sparseQ x
+    //
+    //  except the 2s don't matter because we compute 
+    //  the optimal stepsize anyway
+    COLA_ASSERT(x.size()==b.size() && b.size()==g.size());
+    g = b;
+    for (unsigned i=0; i<denseSize; i++) {
+        for (unsigned j=0; j<denseSize; j++) {
+            g[i] -= (*denseQ)[i*denseSize+j]*x[j];
+        }
     }
+    // sparse part:
+    if(sparseQ) {
+        valarray<double> r(x.size());
+        sparseQ->rightMultiply(x,r);
+        g-=r;
+    }
+    return computeStepSize(g,g);
 }
+// compute optimal step size along descent vector d relative to
+// a gradient related vector g 
+//    stepsize = ( g' d ) / ( d' A d )
+double GradientProjection::computeStepSize(
+        valarray<double> const & g, valarray<double> const & d) const {
+    COLA_ASSERT(g.size()==d.size());
+    valarray<double> Ad;
+    if(sparseQ) {
+        Ad.resize(g.size());
+        sparseQ->rightMultiply(d,Ad);
+    }
+    double const numerator = dotProd(g, d);
+    double denominator = 0;
+    for (unsigned i=0; i<g.size(); i++) {
+        double r = sparseQ ? Ad[i] : 0;
+        if(i<denseSize) { for (unsigned j=0; j<denseSize; j++) {
+            r += (*denseQ)[i*denseSize+j] * d[j];
+        } }
+        denominator += r * d[i];
+    }
+    if(denominator==0) {
+        return 0;
+    }
+    return numerator/(2.*denominator);
+}
+
+bool GradientProjection::runSolver(valarray<double> & result) {
+    bool activeConstraints = false;
+    switch(solveWithMosek) {
+        case Off:
+            activeConstraints = solver->satisfy();
+            for (unsigned i=0;i<vars.size();i++) {
+                result[i]=vars[i]->finalPosition;
+            }
+            break;
+        case Inner:
+#ifdef MOSEK_AVAILABLE
+            float *ba=new float[vars.size()];
+            float *coords=new float[vars.size()];
+            for(unsigned i=0;i<vars.size();i++) {
+                ba[i]=vars[i]->desiredPosition;
+            }
+            mosek_quad_solve_sep(menv,ba,coords);
+            for(unsigned i=0;i<vars.size();i++) {
+                //printf("x[%d]=%f\n",i,coords[i]);
+                result[i]=coords[i];
+            }
+            delete [] ba;
+            delete [] coords;
+            break;
+#endif
+        default:
+            break;
+    }
+    return activeConstraints;
+}
+
 /*
  * 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;
+unsigned GradientProjection::solve(
+        valarray<double> const &linearCoefficients, 
+        valarray<double> &x) {
+    COLA_ASSERT(linearCoefficients.size()==x.size());
+    COLA_ASSERT(x.size()==denseSize);
+    COLA_ASSERT(numStaticVars>=denseSize);
+    COLA_ASSERT(sparseQ==NULL || 
+                (sparseQ!=NULL && (vars.size()==sparseQ->rowSize())) );
 
-	bool converged=false;
+    if(max_iterations==0) return 0;
 
-    IncSolver* solver=NULL;
+    bool converged=false;
 
     solver = setupVPSC();
-    //cerr << "in gradient projection: n=" << n << endl;
-    for (i=0;i<n;i++) {
-        assert(!IS_NAN(place[i]));
-        assert(IS_FINITE(place[i]));
-        vars[i]->desiredPosition=place[i];
-    }
-    try {
-        solver->satisfy();
-    } catch (char const *str) {
-        dumpVPSCException(str,solver);
-	}
-
-    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;
+#ifdef MOSEK_AVAILABLE
+    if(solveWithMosek==Outer) {
+        float* ba=new float[vars.size()];
+        float* xa=new float[vars.size()];
+        for(unsigned i=0;i<vars.size();i++) {
+            ba[i]=-linearCoefficients[i];
+        }
+        mosek_quad_solve_sep(menv,ba,xa);
+        for(unsigned i=0;i<vars.size();i++) {
+            //printf("mosek result x[%d]=%f\n",i,xa[i]);
+            x[i]=xa[i];
+        }
+        delete [] ba;
+        delete [] xa;
+        return 1;
+    }
+#endif
+    // it may be that we have to consider dummy vars, which the caller didn't know
+    // about.  Thus vars.size() may not equal x.size()
+    unsigned n = vars.size();
+    valarray<double> b(n);
+    result.resize(n);
+    
+    // load desired positions into vars, note that we keep desired positions 
+    // already calculated for dummy vars
+    for (unsigned i=0;i<x.size();i++) {
+        COLA_ASSERT(!isNaN(x[i]));
+        COLA_ASSERT(isFinite(x[i]));
+        b[i]=i<linearCoefficients.size()?linearCoefficients[i]:0;
+        result[i]=x[i];
+        if(scaling) {
+            b[i]*=vars[i]->scale;
+            result[i]/=vars[i]->scale;
+        } 
+        if(!vars[i]->fixedDesiredPosition) vars[i]->desiredPosition=result[i];
+    }
+    runSolver(result);
+        
+    valarray<double> g(n); /* gradient */
+    valarray<double> previous(n); /* stored positions */
+    valarray<double> d(n); /* actual descent vector */
+
+#ifdef CHECK_CONVERGENCE_BY_COST
+    double previousCost = DBL_MAX;
+#endif
+    unsigned counter=0;
+    double stepSize;
+    for (; counter<max_iterations&&!converged; counter++) {
+        previous=result;
+        stepSize=0;
+        double alpha=computeSteepestDescentVector(b,result,g);
 
+        //printf("Iteration[%d]\n",counter);
         // move to new unconstrained position
-		for (i=0; i<n; i++) {
-			place[i]-=alpha*g[i];
-            assert(!IS_NAN(place[i]));
-            assert(IS_FINITE(place[i]));
-            vars[i]->desiredPosition=place[i];
-		}
-        for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
-            (*it)->steepestDescent(alpha);
+        for (unsigned i=0; i<n; i++) {
+            // dividing by variable weight is a cheap trick to make these
+            // weights mean something in terms of the descent vector
+            double step=alpha*g[i]/vars[i]->weight;
+            result[i]+=step;
+            //printf("   after unconstrained step: x[%d]=%f\n",i,result[i]);
+            stepSize+=step*step;
+            COLA_ASSERT(!isNaN(result[i]));
+            COLA_ASSERT(isFinite(result[i]));
+            if(!vars[i]->fixedDesiredPosition) vars[i]->desiredPosition=result[i];
         }
 
         //project to constraint boundary
-        try {
-            solver->satisfy();
-        } catch (char const *str) {
-            dumpVPSCException(str,solver);
-        }
-        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];
+        bool constrainedOptimum = false;
+        constrainedOptimum=runSolver(result);
+        stepSize=0;
+        for (unsigned i=0;i<n;i++) {
+            double step = previous[i]-result[i];
+            stepSize+=step*step;
+        }
+        //constrainedOptimum=false;
+        // beta seems, more often than not, to be >1!
+        if(constrainedOptimum) {
+            // The following step limits the step-size in the feasible
+            // direction
+            d = result - previous;
+            const double beta = 0.5*computeStepSize(g, d);
+            // beta > 1.0 takes us back outside the feasible region
+            // beta < 0 clearly not useful and may happen due to numerical imp.
+            //printf("beta=%f\n",beta);
+            if(beta>0&&beta<0.99999) {
+                stepSize=0;
+                for (unsigned i=0; i<n; i++) {
+                    double step=beta*d[i];
+                    result[i]=previous[i]+step;
+                    stepSize+=step*step;
+                }
             }
-            for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
-                (*it)->feasibleDescent(beta);
+        }
+#ifdef CHECK_CONVERGENCE_BY_COST
+        /* This would be the slow way to detect convergence */
+        //if(counter%2) {
+            double cost = computeCost(b,result);
+            printf("     gp[%d] %.15f %.15f\n",counter,previousCost,cost);
+            //COLA_ASSERT(previousCost>cost);
+            if(fabs(previousCost - cost) < tolerance) {
+                converged = true;
             }
+            previousCost = cost;
+        //}
+#else
+        if(stepSize<tolerance) converged = true; 
+#endif
+    }
+    //printf("GP[%d] converged after %d iterations.\n",k,counter);
+    for(unsigned i=0;i<x.size();i++) {
+        x[i]=result[i];
+        if(scaling) {
+            x[i]*=vars[i]->scale;
         }
-		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(solver);
-	return counter;
+    return counter;
 }
 // Setup an instance of the Variable Placement with Separation Constraints
 // for one iteration.
@@ -168,70 +375,109 @@ unsigned GradientProjection::solve(double * b) {
 // global constraint list (including alignment constraints,
 // dir-edge constraints, containment constraints, etc).
 IncSolver* 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);
+    if(nonOverlapConstraints!=None) {
+        if(clusterHierarchy) {
+            //printf("Setup up cluster constraints, dim=%d--------------\n",k);
+            //clusterHierarchy->generateNonOverlapConstraints(k,nonOverlapConstraints,*rs,vars,lcs);
         } else {
-            m=generateYConstraints(n,rs,vs,tmp_cs); 
-        }
-        for(unsigned i=0;i<m;i++) {
-            lcs.push_back(tmp_cs[i]);
+            for(vector<OrthogonalEdgeConstraint*>::iterator i=orthogonalEdges.begin();i!=orthogonalEdges.end();i++) {
+                OrthogonalEdgeConstraint* e=*i;
+                e->generateTopologyConstraints(k,*rs,vars,lcs);
+            }
+            if(k==HORIZONTAL) {
+                // Make rectangles a little bit wider when processing horizontally so that any overlap
+                // resolved horizontally is strictly non-overlapping when processing vertically
+                Rectangle::setXBorder(0.0001);
+                // use rs->size() rather than n because some of the variables may
+                // be dummy vars with no corresponding rectangle
+                generateXConstraints(*rs,vars,lcs,nonOverlapConstraints==Both?true:false); 
+                Rectangle::setXBorder(0);
+            } else {
+                generateYConstraints(*rs,vars,lcs); 
+            }
         }
     }
-    cs = new Constraint*[lcs.size() + gcs.size()];
-    unsigned m = 0 ;
-    for(vector<Constraint*>::iterator ci = lcs.begin();ci!=lcs.end();++ci) {
-        cs[m++] = *ci;
-    }
-    for(vector<Constraint*>::iterator ci = gcs.begin();ci!=gcs.end();++ci) {
-        cs[m++] = *ci;
-    }
-    return new IncSolver(vars.size(),vs,m,cs);
-}
-void GradientProjection::clearDummyVars() {
-    for(DummyVars::iterator i=dummy_vars.begin();i!=dummy_vars.end();++i) {
-        delete *i;
+    cs=gcs;
+    cs.insert(cs.end(),lcs.begin(),lcs.end());
+    switch(solveWithMosek) {
+        case Off:
+            break;
+#ifdef MOSEK_AVAILABLE
+        case Inner:
+            menv = mosek_init_sep_ls(vars.size(),cs);
+            break;
+        case Outer:
+            unsigned n = vars.size();
+            float* lap = new float[n*(n+1)/2];
+            unsigned k=0;
+            for(unsigned i=0;i<n;i++) {
+                for(unsigned j=i;j<n;j++) {
+                    lap[k]=(*denseQ)[i*n+j];
+                    k++;
+                }
+            }
+            menv = mosek_init_sep(lap,n,cs,1);
+            delete [] lap;
+            break;
+#endif
+        default:
+            break;
     }
-    dummy_vars.clear();
+    return new IncSolver(vars,cs);
 }
 void GradientProjection::destroyVPSC(IncSolver *vpsc) {
-    if(acs) {
-        for(AlignmentConstraints::iterator ac=acs->begin(); ac!=acs->end();++ac) {
-            (*ac)->updatePosition();
+    if(ccs) {
+        for(CompoundConstraints::const_iterator c=ccs->begin(); 
+                c!=ccs->end();++c) {
+            (*c)->updatePosition(vpsc::XDIM);
+            (*c)->updatePosition(vpsc::YDIM);
         }
     }
-    unsigned m,n;
-    Constraint** cs = vpsc->getConstraints(m);
-    const Variable* const* vs = vpsc->getVariables(n);
-    delete vpsc;
-    delete [] cs;
-    delete [] vs;
-    for(vector<Constraint*>::iterator i=lcs.begin();i!=lcs.end();++i) {
-            delete *i;
+    if(unsatisfiableConstraints) {
+        unsatisfiableConstraints->clear();
+        for(Constraints::iterator i=cs.begin();i!=cs.end();i++) {
+            Constraint* c=*i;
+            if(c->unsatisfiable) {
+                UnsatisfiableConstraintInfo *ci = new UnsatisfiableConstraintInfo(c);
+                unsatisfiableConstraints->push_back(ci);
+            }
+        }
+    }
+    if(clusterHierarchy) {
+        clusterHierarchy->computeBoundary(*rs);
+    }
+    if(sparseQ) {
+        for(unsigned i=numStaticVars;i<vars.size();i++) {
+            delete vars[i];
+        }
+        vars.resize(numStaticVars);
+        sparseQ=NULL;
+    }
+    for(vector<Constraint*>::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;
+    delete vpsc;
+#ifdef MOSEK_AVAILABLE
+    if(solveWithMosek!=Off) mosek_delete(menv);
+#endif
+}
+void GradientProjection::straighten(
+    cola::SparseMatrix const * Q, 
+    vector<SeparationConstraint*> const & cs,
+    vector<straightener::Node*> const & snodes) 
+{
+    COLA_ASSERT(Q->rowSize()==snodes.size());
+    COLA_ASSERT(vars.size()==numStaticVars);
+    sparseQ = Q;
+    for(unsigned i=numStaticVars;i<snodes.size();i++) {
+        Variable* v=new vpsc::Variable(i,snodes[i]->pos[k],1);
+        COLA_ASSERT(v->desiredPosition==snodes[i]->pos[k]);
+        vars.push_back(v);
+    }
+    COLA_ASSERT(lcs.size()==0);
+    for(vector<SeparationConstraint*>::const_iterator i=cs.begin();i!=cs.end();i++) {
+        (*i)->generateSeparationConstraints(k, vars, lcs, *rs); 
     }
 }
-
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :
+} // namespace cola