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