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