src/libcola/gradient_projection.cpp

   1 /**********************************************************
   2  *
   3  * Solve a quadratic function f(X) = X' A X + b X
   4  * subject to a set of separation constraints cs
   5  *
   6  * Tim Dwyer, 2006
   7  **********************************************************/
   8
   9 #include <math.h>
  10 #include <stdlib.h>
  11 #include <time.h>
  12 #include <stdio.h>
  13 #include <float.h>
  14 #include <cassert>
  15 #include <libvpsc/solve_VPSC.h>
  16 #include <libvpsc/variable.h>
  17 #include <libvpsc/constraint.h>
  18 #include "gradient_projection.h"
  19 #include <iostream>
  20
  21 using namespace std;
  22 //#define CONMAJ_LOGGING 1
  23
  24 static void dumpVPSCException(char const *str, IncVPSC* vpsc) {
  25     cerr<<str<<endl;
  26     unsigned m;
  27     Constraint** cs = vpsc->getConstraints(m);
  28     for(unsigned i=0;i<m;i++) {
  29         cerr << *cs[i] << endl;
  30     }
  31 }
  32 /*
  33  * Use gradient-projection to solve an instance of
  34  * the Variable Placement with Separation Constraints problem.
  35  * Uses sparse matrix techniques to handle pairs of dummy
  36  * vars.
  37  */
  38 unsigned GradientProjection::solve(double * b) {
  39         unsigned i,j,counter;
  40         if(max_iterations==0) return 0;
  41
  42         bool converged=false;
  43
  44     IncVPSC* vpsc=NULL;
  45
  46     vpsc = setupVPSC();
  47     //cerr << "in gradient projection: n=" << n << endl;
  48     for (i=0;i<n;i++) {
  49         assert(!isnan(place[i]));
  50         assert(!isinf(place[i]));
  51         vars[i]->desiredPosition=place[i];
  52     }
  53     try {
  54         vpsc->satisfy();
  55     } catch (char const *str) {
  56         dumpVPSCException(str,vpsc);
  57         }
  58
  59     for (i=0;i<n;i++) {
  60         place[i]=vars[i]->position();
  61     }
  62     for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
  63         (*it)->updatePosition();
  64     }
  65
  66         for (counter=0; counter<max_iterations&&!converged; counter++) {
  67                 converged=true;
  68                 // find steepest descent direction
  69         //  g = 2 ( b - Ax )
  70                 for (i=0; i<n; i++) {
  71                         old_place[i]=place[i];
  72                         g[i] = b[i];
  73                         for (j=0; j<n; j++) {
  74                                 g[i] -= A[i][j]*place[j];
  75                         }
  76             g[i] *= 2.0;
  77                 }
  78         for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
  79             (*it)->computeDescentVector();
  80         }
  81         // compute step size: alpha = ( g' g ) / ( 2 g' A g )
  82         //   g terms for dummy vars cancel out so don't consider
  83                 double numerator = 0, denominator = 0, r;
  84                 for (i=0; i<n; i++) {
  85                         numerator += g[i]*g[i];
  86                         r=0;
  87                         for (j=0; j<n; j++) {
  88                                 r += A[i][j]*g[j];
  89                         }
  90                         denominator -= 2.0 * r*g[i];
  91                 }
  92                 double alpha = numerator/denominator;
  93
  94         // move to new unconstrained position
  95                 for (i=0; i<n; i++) {
  96                         place[i]-=alpha*g[i];
  97             assert(!isnan(place[i]));
  98             assert(!isinf(place[i]));
  99             vars[i]->desiredPosition=place[i];
 100                 }
 101         for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
 102             (*it)->steepestDescent(alpha);
 103         }
 104
 105         //project to constraint boundary
 106         try {
 107             vpsc->satisfy();
 108         } catch (char const *str) {
 109             dumpVPSCException(str,vpsc);
 110         }
 111         for (i=0;i<n;i++) {
 112             place[i]=vars[i]->position();
 113         }
 114         for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
 115             (*it)->updatePosition();
 116         }
 117         // compute d, the vector from last pnt to projection pnt
 118                 for (i=0; i<n; i++) {
 119                         d[i]=place[i]-old_place[i];
 120                 }
 121                 // now compute beta, optimal step size from last pnt to projection pnt
 122         //   beta = ( g' d ) / ( 2 d' A d )
 123                 numerator = 0, denominator = 0;
 124                 for (i=0; i<n; i++) {
 125                         numerator += g[i] * d[i];
 126                         r=0;
 127                         for (j=0; j<n; j++) {
 128                                 r += A[i][j] * d[j];
 129                         }
 130                         denominator += 2.0 * r * d[i];
 131                 }
 132         for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
 133             (*it)->betaCalc(numerator,denominator);
 134         }
 135                 double beta = numerator/denominator;
 136
 137         // beta > 1.0 takes us back outside the feasible region
 138         // beta < 0 clearly not useful and may happen due to numerical imp.
 139         if(beta>0&&beta<1.0) {
 140             for (i=0; i<n; i++) {
 141                 place[i]=old_place[i]+beta*d[i];
 142             }
 143             for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
 144                 (*it)->feasibleDescent(beta);
 145             }
 146         }
 147                 double test=0;
 148                 for (i=0; i<n; i++) {
 149                         test += fabs(place[i]-old_place[i]);
 150                 }
 151         for (DummyVars::iterator it=dummy_vars.begin();it!=dummy_vars.end();++it){
 152             test += (*it)->absoluteDisplacement();
 153         }
 154                 if(test>tolerance) {
 155                         converged=false;
 156                 }
 157         }
 158     destroyVPSC(vpsc);
 159         return counter;
 160 }
 161 // Setup an instance of the Variable Placement with Separation Constraints
 162 // for one iteration.
 163 // Generate transient local constraints --- such as non-overlap constraints
 164 // --- that are only relevant to one iteration, and merge these with the
 165 // global constraint list (including alignment constraints,
 166 // dir-edge constraints, containment constraints, etc).
 167 IncVPSC* GradientProjection::setupVPSC() {
 168     Constraint **cs;
 169     //assert(lcs.size()==0);
 170
 171     for(DummyVars::iterator dit=dummy_vars.begin();
 172             dit!=dummy_vars.end(); ++dit) {
 173         (*dit)->setupVPSC(vars,lcs);
 174     }
 175     Variable** vs = new Variable*[vars.size()];
 176     for(unsigned i=0;i<vars.size();i++) {
 177         vs[i]=vars[i];
 178     }
 179     if(nonOverlapConstraints) {
 180         Constraint** tmp_cs=NULL;
 181         unsigned m=0;
 182         if(k==HORIZONTAL) {
 183             Rectangle::setXBorder(0.0001);
 184             m=generateXConstraints(n,rs,vs,tmp_cs,true);
 185             Rectangle::setXBorder(0);
 186         } else {
 187             m=generateYConstraints(n,rs,vs,tmp_cs);
 188         }
 189         for(unsigned i=0;i<m;i++) {
 190             lcs.push_back(tmp_cs[i]);
 191         }
 192     }
 193     cs = new Constraint*[lcs.size() + gcs.size()];
 194     unsigned m = 0 ;
 195     for(Constraints::iterator ci = lcs.begin();ci!=lcs.end();++ci) {
 196         cs[m++] = *ci;
 197     }
 198     for(Constraints::iterator ci = gcs.begin();ci!=gcs.end();++ci) {
 199         cs[m++] = *ci;
 200     }
 201     return new IncVPSC(vars.size(),vs,m,cs);
 202 }
 203 void GradientProjection::clearDummyVars() {
 204     for(DummyVars::iterator i=dummy_vars.begin();i!=dummy_vars.end();++i) {
 205         delete *i;
 206     }
 207     dummy_vars.clear();
 208 }
 209 void GradientProjection::destroyVPSC(IncVPSC *vpsc) {
 210     if(acs) {
 211         for(AlignmentConstraints::iterator ac=acs->begin(); ac!=acs->end();++ac) {
 212             (*ac)->updatePosition();
 213         }
 214     }
 215     unsigned m,n;
 216     Constraint** cs = vpsc->getConstraints(m);
 217     const Variable* const* vs = vpsc->getVariables(n);
 218     delete vpsc;
 219     delete [] cs;
 220     delete [] vs;
 221     for(Constraints::iterator i=lcs.begin();i!=lcs.end();i++) {
 222             delete *i;
 223     }
 224     lcs.clear();
 225     //cout << " Vars count = " << vars.size() << " Dummy vars cnt=" << dummy_vars.size() << endl;
 226     vars.resize(vars.size()-dummy_vars.size()*2);
 227     for(DummyVars::iterator i=dummy_vars.begin();i!=dummy_vars.end();++i) {
 228         DummyVarPair* p = *i;
 229         delete p->left;
 230         delete p->right;
 231     }
 232 }
 233
 234 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :