diff options
Diffstat (limited to 'src/libcola/gradient_projection.cpp')
| -rw-r--r-- | src/libcola/gradient_projection.cpp | 483 |
1 files changed, 0 insertions, 483 deletions
diff --git a/src/libcola/gradient_projection.cpp b/src/libcola/gradient_projection.cpp deleted file mode 100644 index 50b8d27d7..000000000 --- a/src/libcola/gradient_projection.cpp +++ /dev/null @@ -1,483 +0,0 @@ -/* - * 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 - * - **********************************************************/ - -#include <iostream> -#include <cmath> -#include <ctime> - -#include <libvpsc/solve_VPSC.h> -#include <libvpsc/variable.h> -#include <libvpsc/isnan.h> -#include <libvpsc/constraint.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; -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); -} - -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. - */ -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())) ); - - if(max_iterations==0) return 0; - - bool converged=false; - - solver = setupVPSC(); -#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 (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 - 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; - } - } - } -#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; - } - } - destroyVPSC(solver); - 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). -IncSolver* GradientProjection::setupVPSC() { - if(nonOverlapConstraints!=None) { - if(clusterHierarchy) { - //printf("Setup up cluster constraints, dim=%d--------------\n",k); - //clusterHierarchy->generateNonOverlapConstraints(k,nonOverlapConstraints,*rs,vars,lcs); - } else { - 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=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; - } - return new IncSolver(vars,cs); -} -void GradientProjection::destroyVPSC(IncSolver *vpsc) { - if(ccs) { - for(CompoundConstraints::const_iterator c=ccs->begin(); - c!=ccs->end();++c) { - (*c)->updatePosition(vpsc::XDIM); - (*c)->updatePosition(vpsc::YDIM); - } - } - 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(); - 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); - } -} -} // namespace cola |