1 #include "cola.h"
2 #include "conjugate_gradient.h"
3 #include "straightener.h"
4 #include "shortest_paths.h"
6 namespace cola {
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);
41 done.reset();
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 }
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 */
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;
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 }
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