src/libcola/cola.cpp

   1 #include "cola.h"
   2 #include "conjugate_gradient.h"
   3 #include "straightener.h"
   4 #include "shortest_paths.h"
   5
   6 namespace cola {
   7
   8 /**
   9  * Find the euclidean distance between a pair of dummy variables
  10  */
  11 inline double dummy_var_euclidean_dist(GradientProjection* gpx, GradientProjection* gpy, unsigned i) {
  12     double dx = gpx->dummy_vars[i]->place_r - gpx->dummy_vars[i]->place_l,
  13         dy = gpy->dummy_vars[i]->place_r - gpy->dummy_vars[i]->place_l;
  14     return sqrt(dx*dx + dy*dy);
  15 }
  16 ConstrainedMajorizationLayout
  17 ::ConstrainedMajorizationLayout(
  18         vector<Rectangle*>& rs,
  19         vector<Edge>& es,
  20         double* eweights,
  21         double idealLength,
  22         TestConvergence& done)
  23     : constrainedLayout(false),
  24       n(rs.size()),
  25       lapSize(n), lap2(new double*[lapSize]),
  26       Q(lap2), Dij(new double*[lapSize]),
  27       tol(0.0001),
  28       done(done),
  29       X(new double[n]),
  30       Y(new double[n]),
  31       clusters(NULL),
  32       linearConstraints(NULL),
  33       gpX(NULL),
  34       gpY(NULL),
  35       straightenEdges(NULL)
  36 {
  37     assert(rs.size()==n);
  38     boundingBoxes = new Rectangle*[rs.size()];
  39     copy(rs.begin(),rs.end(),boundingBoxes);
  40
  41     done.reset();
  42
  43     double** D=new double*[n];
  44     for(unsigned i=0;i<n;i++) {
  45         D[i]=new double[n];
  46     }
  47     shortest_paths::johnsons(n,D,es,eweights);
  48     edge_length = idealLength;
  49     // Lij_{i!=j}=1/(Dij^2)
  50     //
  51     for(unsigned i = 0; i<n; i++) {
  52         X[i]=rs[i]->getCentreX();
  53         Y[i]=rs[i]->getCentreY();
  54         double degree = 0;
  55         lap2[i]=new double[n];
  56         Dij[i]=new double[n];
  57         for(unsigned j=0;j<n;j++) {
  58             double w = edge_length * D[i][j];
  59             Dij[i][j]=w;
  60             if(i==j) continue;
  61             degree+=lap2[i][j]=w>1e-30?1.f/(w*w):0;
  62         }
  63         lap2[i][i]=-degree;
  64         delete [] D[i];
  65     }
  66     delete [] D;
  67 }
  68
  69 void
  70 ConstrainedMajorizationLayout
  71 ::setupDummyVars() {
  72     if(clusters==NULL) return;
  73     double* coords[2]={X,Y};
  74     GradientProjection* gp[2]={gpX,gpY};
  75     for(unsigned k=0;k<2;k++) {
  76         gp[k]->clearDummyVars();
  77         if(clusters) {
  78             for(Clusters::iterator cit=clusters->begin();
  79                     cit!=clusters->end(); ++cit) {
  80                 Cluster *c = *cit;
  81                 DummyVarPair* p = new DummyVarPair(edge_length);
  82                 gp[k]->dummy_vars.push_back(p);
  83                 double minPos=DBL_MAX, maxPos=-DBL_MAX;
  84                 for(Cluster::iterator vit=c->begin();
  85                         vit!=c->end(); ++vit) {
  86                     double pos = coords[k][*vit];
  87                     minPos=min(pos,minPos);
  88                     maxPos=max(pos,maxPos);
  89                     p->leftof.push_back(make_pair(*vit,0));
  90                     p->rightof.push_back(make_pair(*vit,0));
  91                 }
  92                 p->place_l = minPos;
  93                 p->place_r = maxPos;
  94             }
  95         }
  96     }
  97     for(unsigned k=0;k<2;k++) {
  98         unsigned n_d = gp[k]->dummy_vars.size();
  99         if(n_d > 0) {
 100             for(unsigned i=0; i<n_d; i++) {
 101                 gp[k]->dummy_vars[i]->computeLinearTerm(dummy_var_euclidean_dist(gpX, gpY, i));
 102             }
 103         }
 104     }
 105 }
 106 void ConstrainedMajorizationLayout::majlayout(
 107         double** Dij, GradientProjection* gp, double* coords)
 108 {
 109     double b[n];
 110     fill(b,b+n,0);
 111     majlayout(Dij,gp,coords,b);
 112 }
 113 void ConstrainedMajorizationLayout::majlayout(
 114         double** Dij, GradientProjection* gp, double* coords, double* b)
 115 {
 116     double L_ij,dist_ij,degree;
 117     /* compute the vector b */
 118     /* multiply on-the-fly with distance-based laplacian */
 119     for (unsigned i = 0; i < n; i++) {
 120         degree = 0;
 121         if(i<lapSize) {
 122             for (unsigned j = 0; j < lapSize; j++) {
 123                 if (j == i) continue;
 124                 dist_ij = euclidean_distance(i, j);
 125                 if (dist_ij > 1e-30 && Dij[i][j] > 1e-30) {     /* skip zero distances */
 126                     /* calculate L_ij := w_{ij}*d_{ij}/dist_{ij} */
 127                     L_ij = 1.0 / (dist_ij * Dij[i][j]);
 128                     degree -= L_ij;
 129                     b[i] += L_ij * coords[j];
 130                 }
 131             }
 132             b[i] += degree * coords[i];
 133         }
 134         assert(!isnan(b[i]));
 135     }
 136     if(constrainedLayout) {
 137         setupDummyVars();
 138         gp->solve(b);
 139     } else {
 140         conjugate_gradient(lap2, coords, b, n, tol, n);
 141     }
 142     moveBoundingBoxes();
 143 }
 144 inline double ConstrainedMajorizationLayout
 145 ::compute_stress(double **Dij) {
 146     double sum = 0, d, diff;
 147     for (unsigned i = 1; i < lapSize; i++) {
 148         for (unsigned j = 0; j < i; j++) {
 149             d = Dij[i][j];
 150             diff = d - euclidean_distance(i,j);
 151             sum += diff*diff / (d*d);
 152         }
 153     }
 154     if(clusters!=NULL) {
 155         for(unsigned i=0; i<gpX->dummy_vars.size(); i++) {
 156             sum += gpX->dummy_vars[i]->stress(dummy_var_euclidean_dist(gpX, gpY, i));
 157         }
 158     }
 159     return sum;
 160 }
 161 /*
 162 void ConstrainedMajorizationLayout
 163 ::addLinearConstraints(LinearConstraints* linearConstraints) {
 164     n=lapSize+linearConstraints->size();
 165     Q=new double*[n];
 166     X=new double[n];
 167     Y=new double[n];
 168     for(unsigned i = 0; i<n; i++) {
 169         X[i]=rs[i]->getCentreX();
 170         Y[i]=rs[i]->getCentreY();
 171         Q[i]=new double[n];
 172         for(unsigned j=0; j<n; j++) {
 173             if(i<lapSize&&j<lapSize) {
 174                 Q[i][j]=lap2[i][j];
 175             } else {
 176                 Q[i][j]=0;
 177             }
 178         }
 179     }
 180     for(LinearConstraints::iterator i=linearConstraints->begin();
 181            i!= linearConstraints->end();i++) {
 182         LinearConstraint* c=*i;
 183         Q[c->u][c->u]+=c->w*c->duu;
 184         Q[c->u][c->v]+=c->w*c->duv;
 185         Q[c->u][c->b]+=c->w*c->dub;
 186         Q[c->v][c->u]+=c->w*c->duv;
 187         Q[c->v][c->v]+=c->w*c->dvv;
 188         Q[c->v][c->b]+=c->w*c->dvb;
 189         Q[c->b][c->b]+=c->w*c->dbb;
 190         Q[c->b][c->u]+=c->w*c->dub;
 191         Q[c->b][c->v]+=c->w*c->dvb;
 192     }
 193 }
 194 */
 195
 196 bool ConstrainedMajorizationLayout::run() {
 197     /*
 198     for(unsigned i=0;i<n;i++) {
 199         for(unsigned j=0;j<n;j++) {
 200             cout << lap2[i][j] << " ";
 201         }
 202         cout << endl;
 203     }
 204     */
 205     do {
 206         /* Axis-by-axis optimization: */
 207         if(straightenEdges) {
 208             straighten(*straightenEdges,HORIZONTAL);
 209             straighten(*straightenEdges,VERTICAL);
 210         } else {
 211             majlayout(Dij,gpX,X);
 212             majlayout(Dij,gpY,Y);
 213         }
 214     } while(!done(compute_stress(Dij),X,Y));
 215     return true;
 216 }
 217 static bool straightenToProjection=true;
 218 void ConstrainedMajorizationLayout::straighten(vector<straightener::Edge*>& sedges, Dim dim) {
 219         vector<straightener::Node*> snodes;
 220         for (unsigned i=0;i<lapSize;i++) {
 221                 snodes.push_back(new straightener::Node(i,boundingBoxes[i]));
 222         }
 223     SimpleConstraints cs;
 224     straightener::generateConstraints(snodes,sedges,cs,dim);
 225     n=snodes.size();
 226     Q=new double*[n];
 227     delete [] X;
 228     delete [] Y;
 229     X=new double[n];
 230     Y=new double[n];
 231     for(unsigned i = 0; i<n; i++) {
 232         X[i]=snodes[i]->x;
 233         Y[i]=snodes[i]->y;
 234         Q[i]=new double[n];
 235         for(unsigned j=0; j<n; j++) {
 236             if(i<lapSize&&j<lapSize) {
 237                 Q[i][j]=lap2[i][j];
 238             } else {
 239                 Q[i][j]=0;
 240             }
 241         }
 242     }
 243     LinearConstraints linearConstraints;
 244     for(unsigned i=0;i<sedges.size();i++) {
 245         sedges[i]->nodePath(snodes);
 246         vector<unsigned>& path=sedges[i]->path;
 247         // take u and v as the ends of the line
 248         //unsigned u=path[0];
 249         //unsigned v=path[path.size()-1];
 250         double total_length=0;
 251         for(unsigned j=1;j<path.size();j++) {
 252             unsigned u=path[j-1], v=path[j];
 253             total_length+=euclidean_distance(u,v);
 254         }
 255         for(unsigned j=1;j<path.size()-1;j++) {
 256             // find new u and v for each line segment
 257             unsigned u=path[j-1];
 258             unsigned b=path[j];
 259             unsigned v=path[j+1];
 260             double weight=-0.01;
 261             double wub=euclidean_distance(u,b)/total_length;
 262             double wbv=euclidean_distance(b,v)/total_length;
 263             linearConstraints.push_back(new cola::LinearConstraint(u,v,b,weight,wub,wbv,X,Y));
 264         }
 265     }
 266     //cout << "Generated "<<linearConstraints.size()<< " linear constraints"<<endl;
 267     assert(snodes.size()==lapSize+linearConstraints.size());
 268     double b[n],*coords=dim==HORIZONTAL?X:Y,dist_ub,dist_bv;
 269     fill(b,b+n,0);
 270     for(LinearConstraints::iterator i=linearConstraints.begin();
 271            i!= linearConstraints.end();i++) {
 272         LinearConstraint* c=*i;
 273         if(straightenToProjection) {
 274             Q[c->u][c->u]+=c->w*c->duu;
 275             Q[c->u][c->v]+=c->w*c->duv;
 276             Q[c->u][c->b]+=c->w*c->dub;
 277             Q[c->v][c->u]+=c->w*c->duv;
 278             Q[c->v][c->v]+=c->w*c->dvv;
 279             Q[c->v][c->b]+=c->w*c->dvb;
 280             Q[c->b][c->b]+=c->w*c->dbb;
 281             Q[c->b][c->u]+=c->w*c->dub;
 282             Q[c->b][c->v]+=c->w*c->dvb;
 283         } else {
 284             double wub=edge_length*c->frac_ub;
 285             double wbv=edge_length*c->frac_bv;
 286             dist_ub=euclidean_distance(c->u,c->b)*wub;
 287             dist_bv=euclidean_distance(c->b,c->v)*wbv;
 288             wub=max(wub,0.00001);
 289             wbv=max(wbv,0.00001);
 290             dist_ub=max(dist_ub,0.00001);
 291             dist_bv=max(dist_bv,0.00001);
 292             wub=1/(wub*wub);
 293             wbv=1/(wbv*wbv);
 294             Q[c->u][c->u]-=wub;
 295             Q[c->u][c->b]+=wub;
 296             Q[c->v][c->v]-=wbv;
 297             Q[c->v][c->b]+=wbv;
 298             Q[c->b][c->b]-=wbv+wub;
 299             Q[c->b][c->u]+=wub;
 300             Q[c->b][c->v]+=wbv;
 301
 302             b[c->u]+=(coords[c->b]-coords[c->u]) / dist_ub;
 303             b[c->v]+=(coords[c->b]-coords[c->v]) / dist_bv;
 304             b[c->b]+=coords[c->u] / dist_ub + coords[c->v] / dist_bv
 305                    - coords[c->b] / dist_ub - coords[c->b] / dist_bv;
 306         }
 307     }
 308         GradientProjection gp(dim,n,Q,coords,tol,100,
 309             (AlignmentConstraints*)NULL,false,(vpsc::Rectangle**)NULL,(PageBoundaryConstraints*)NULL,&cs);
 310     constrainedLayout = true;
 311     majlayout(Dij,&gp,coords,b);
 312     for(unsigned i=0;i<sedges.size();i++) {
 313         sedges[i]->createRouteFromPath(X,Y);
 314         sedges[i]->dummyNodes.clear();
 315         sedges[i]->path.clear();
 316     }
 317     for(unsigned i=0;i<cs.size();i++) {
 318         delete cs[i];
 319     }
 320     for(unsigned i=0;i<linearConstraints.size();i++) {
 321         delete linearConstraints[i];
 322     }
 323     for(unsigned i=0;i<snodes.size();i++) {
 324         delete snodes[i];
 325     }
 326     for(unsigned i = 0; i<n; i++) {
 327         delete [] Q[i];
 328     }
 329     delete [] Q;
 330     snodes.resize(lapSize);
 331 }
 332
 333 void ConstrainedMajorizationLayout::setupConstraints(
 334         AlignmentConstraints* acsx, AlignmentConstraints* acsy,
 335         bool avoidOverlaps,
 336         PageBoundaryConstraints* pbcx, PageBoundaryConstraints* pbcy,
 337         SimpleConstraints* scx, SimpleConstraints* scy,
 338         Clusters* cs,
 339         vector<straightener::Edge*>* straightenEdges) {
 340     constrainedLayout = true;
 341     this->avoidOverlaps = avoidOverlaps;
 342     if(cs) {
 343         clusters=cs;
 344     }
 345         gpX=new GradientProjection(
 346             HORIZONTAL,n,Q,X,tol,100,acsx,avoidOverlaps,boundingBoxes,pbcx,scx);
 347         gpY=new GradientProjection(
 348             VERTICAL,n,Q,Y,tol,100,acsy,avoidOverlaps,boundingBoxes,pbcy,scy);
 349     this->straightenEdges = straightenEdges;
 350 }
 351 } // namespace cola
 352 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=4:softtabstop=4