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