src/libcola/cola.cpp

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