diff options
Diffstat (limited to 'src/libcola/gradient_projection.cpp')
| -rw-r--r-- | src/libcola/gradient_projection.cpp | 234 |
1 files changed, 234 insertions, 0 deletions
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 : |
