src/libcola/gradient_projection.h

   1 #ifndef _GRADIENT_PROJECTION_H
   2 #define _GRADIENT_PROJECTION_H
   3
   4 #include <libvpsc/solve_VPSC.h>
   5 #include <libvpsc/variable.h>
   6 #include <libvpsc/constraint.h>
   7 #include <libvpsc/generate-constraints.h>
   8 #include <vector>
   9 #include <iostream>
  10 #include <math.h>
  11
  12 using namespace std;
  13
  14 typedef vector<Constraint*> Constraints;
  15 typedef vector<Variable*> Variables;
  16 typedef vector<pair<unsigned,double> > OffsetList;
  17
  18 class SimpleConstraint {
  19 public:
  20     SimpleConstraint(unsigned l, unsigned r, double g)
  21         : left(l), right(r), gap(g)  {}
  22     unsigned left;
  23     unsigned right;
  24     double gap;
  25 };
  26 typedef vector<SimpleConstraint*> SimpleConstraints;
  27 class AlignmentConstraint {
  28 friend class GradientProjection;
  29 public:
  30     AlignmentConstraint(double pos) : position(pos), variable(NULL) {}
  31     void updatePosition() {
  32         position = variable->position();
  33     }
  34     OffsetList offsets;
  35     void* guide;
  36     double position;
  37 private:
  38     Variable* variable;
  39 };
  40 typedef vector<AlignmentConstraint*> AlignmentConstraints;
  41
  42 class PageBoundaryConstraints {
  43 public:
  44     PageBoundaryConstraints(double lm, double rm, double w)
  45         : leftMargin(lm), rightMargin(rm), weight(w) { }
  46     void createVarsAndConstraints(Variables &vs, Constraints &cs) {
  47         Variable* vl, * vr;
  48         // create 2 dummy vars, based on the dimension we are in
  49         vs.push_back(vl=new Variable(vs.size(), leftMargin, weight));
  50         vs.push_back(vr=new Variable(vs.size(), rightMargin, weight));
  51
  52         // for each of the "real" variables, create a constraint that puts that var
  53         // between our two new dummy vars, depending on the dimension.
  54         for(OffsetList::iterator o=offsets.begin(); o!=offsets.end(); ++o)  {
  55             cs.push_back(new Constraint(vl, vs[o->first], o->second));
  56             cs.push_back(new Constraint(vs[o->first], vr, o->second));
  57         }
  58     }
  59     OffsetList offsets;
  60 private:
  61     double leftMargin;
  62     double rightMargin;
  63     double weight;
  64 };
  65
  66 typedef vector<pair<unsigned,double> > CList;
  67 /**
  68  * A DummyVarPair is a pair of variables with an ideal distance between them and which have no
  69  * other interaction with other variables apart from through constraints.  This means that
  70  * the entries in the laplacian matrix for dummy vars and other vars would be 0 - thus, sparse
  71  * matrix techniques can be used in laplacian operations.
  72  * The constraints are specified by a two lists of pairs of variable indexes and required separation.
  73  * The two lists are:
  74  *   leftof: variables to which left must be to the left of,
  75  *   rightof: variables to which right must be to the right of.
  76  */
  77 class DummyVarPair {
  78 public:
  79     DummyVarPair(double desiredDist) : dist(desiredDist), lap2(1.0/(desiredDist*desiredDist)) { }
  80     CList leftof; // variables to which left dummy var must be to the left of
  81     CList rightof; // variables to which right dummy var must be to the right of
  82     double place_l;
  83     double place_r;
  84     void computeLinearTerm(double euclideanDistance) {
  85         if(euclideanDistance > 1e-30) {
  86             b = place_r - place_l;
  87             b /= euclideanDistance * dist;
  88         } else { b=0; }
  89     }
  90     double stress(double euclideanDistance) {
  91         double diff = dist - euclideanDistance;
  92         return diff*diff / (dist*dist);
  93     }
  94 private:
  95 friend class GradientProjection;
  96     /**
  97      * Setup vars and constraints for an instance of the VPSC problem.
  98      * Adds generated vars and constraints to the argument vectors.
  99      */
 100     void setupVPSC(Variables &vars, Constraints &cs) {
 101         double weight=1;
 102         left = new Variable(vars.size(),place_l,weight);
 103         vars.push_back(left);
 104         right = new Variable(vars.size(),place_r,weight);
 105         vars.push_back(right);
 106         for(CList::iterator cit=leftof.begin();
 107                 cit!=leftof.end(); ++cit) {
 108             Variable* v = vars[(*cit).first];
 109             cs.push_back(new Constraint(left,v,(*cit).second));
 110         }
 111         for(CList::iterator cit=rightof.begin();
 112                 cit!=rightof.end(); ++cit) {
 113             Variable* v = vars[(*cit).first];
 114             cs.push_back(new Constraint(v,right,(*cit).second));
 115         }
 116     }
 117     /**
 118      * Extract the result of a VPSC solution to the variable positions
 119      */
 120     void updatePosition() {
 121         place_l=left->position();
 122         place_r=right->position();
 123     }
 124     /**
 125      * Compute the descent vector, also copying the current position to old_place for
 126      * future reference
 127      */
 128     void computeDescentVector() {
 129         old_place_l=place_l;
 130         old_place_r=place_r;
 131         g = 2.0 * ( b + lap2 * ( place_l - place_r ) );
 132     }
 133     /**
 134      * move in the direction of steepest descent (based on g computed by computeGradient)
 135      * alpha is the step size.
 136      */
 137     void steepestDescent(double alpha) {
 138         place_l -= alpha*g;
 139         place_r += alpha*g;
 140         left->desiredPosition=place_l;
 141         right->desiredPosition=place_r;
 142     }
 143     /**
 144      * add dummy vars' contribution to numerator and denominator for
 145      * beta step size calculation
 146      */
 147     void betaCalc(double &numerator, double &denominator) {
 148         double dl = place_l-old_place_l,
 149                dr = place_r-old_place_r,
 150                r = 2.0 * lap2 * ( dr - dl );
 151         numerator += g * ( dl - dr );
 152         denominator += r*dl - r * dr;
 153     }
 154     /**
 155      * move by stepsize beta from old_place to place
 156      */
 157     void feasibleDescent(double beta) {
 158         left->desiredPosition = place_l = old_place_l + beta*(place_l - old_place_l);
 159         right->desiredPosition = place_r = old_place_r + beta*(place_r - old_place_r);
 160     }
 161     double absoluteDisplacement() {
 162         return fabs(place_l - old_place_l) + fabs(place_r - old_place_r);
 163     }
 164     double dist; // ideal distance between vars
 165     double b; // linear coefficient in quad form for left (b_right = -b)
 166     Variable* left; // Variables used in constraints
 167     Variable* right;
 168     double lap2; // laplacian entry
 169     double g; // descent vec for quad form for left (g_right = -g)
 170     double old_place_l; // old_place is where the descent vec g was computed
 171     double old_place_r;
 172 };
 173 typedef vector<DummyVarPair*> DummyVars;
 174
 175 enum Dim { HORIZONTAL, VERTICAL };
 176
 177 class GradientProjection {
 178 public:
 179         GradientProjection(
 180         const Dim k,
 181                 unsigned n,
 182                 double** A,
 183                 double* x,
 184                 double tol,
 185                 unsigned max_iterations,
 186         AlignmentConstraints* acs=NULL,
 187         bool nonOverlapConstraints=false,
 188         Rectangle** rs=NULL,
 189         PageBoundaryConstraints *pbc = NULL,
 190         SimpleConstraints *sc = NULL)
 191             : k(k), n(n), A(A), place(x), rs(rs),
 192               nonOverlapConstraints(nonOverlapConstraints),
 193               tolerance(tol), acs(acs), max_iterations(max_iterations),
 194               g(new double[n]), d(new double[n]), old_place(new double[n]),
 195               constrained(false)
 196     {
 197         for(unsigned i=0;i<n;i++) {
 198             vars.push_back(new Variable(i,1,1));
 199         }
 200         if(acs) {
 201             for(AlignmentConstraints::iterator iac=acs->begin();
 202                     iac!=acs->end();++iac) {
 203                 AlignmentConstraint* ac=*iac;
 204                 Variable *v=ac->variable=new Variable(vars.size(),ac->position,0.0001);
 205                 vars.push_back(v);
 206                 for(OffsetList::iterator o=ac->offsets.begin();
 207                         o!=ac->offsets.end();
 208                         o++) {
 209                     gcs.push_back(new Constraint(v,vars[o->first],o->second,true));
 210                 }
 211             }
 212         }
 213         if (pbc)  {
 214             pbc->createVarsAndConstraints(vars,gcs);
 215         }
 216         if (sc) {
 217             for(SimpleConstraints::iterator c=sc->begin(); c!=sc->end();++c) {
 218                 gcs.push_back(new Constraint(
 219                         vars[(*c)->left],vars[(*c)->right],(*c)->gap));
 220             }
 221         }
 222         if(!gcs.empty() || nonOverlapConstraints) {
 223             constrained=true;
 224         }
 225         }
 226     ~GradientProjection() {
 227         delete [] g;
 228         delete [] d;
 229         delete [] old_place;
 230         for(Constraints::iterator i(gcs.begin()); i!=gcs.end(); i++) {
 231             delete *i;
 232         }
 233         gcs.clear();
 234         for(unsigned i=0;i<vars.size();i++) {
 235             delete vars[i];
 236         }
 237     }
 238     void clearDummyVars();
 239         unsigned solve(double* b);
 240     DummyVars dummy_vars; // special vars that must be considered in Lapl.
 241 private:
 242     IncVPSC* setupVPSC();
 243     void destroyVPSC(IncVPSC *vpsc);
 244     Dim k;
 245         unsigned n; // number of actual vars
 246         double** A; // Graph laplacian matrix
 247     double* place;
 248         Variables vars; // all variables
 249                           // computations
 250     Constraints gcs; /* global constraints - persist throughout all
 251                                 iterations */
 252     Constraints lcs; /* local constraints - only for current iteration */
 253     Rectangle** rs;
 254     bool nonOverlapConstraints;
 255     double tolerance;
 256     AlignmentConstraints* acs;
 257     unsigned max_iterations;
 258         double* g; /* gradient */
 259         double* d;
 260         double* old_place;
 261     bool constrained;
 262 };
 263
 264 #endif /* _GRADIENT_PROJECTION_H */
 265
 266 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4 :