6 * $Date: 2012-07-07 17:14:54 +0200 (Sa, 07. Jul 2012) $
7 ***************************************************************/
10 * \brief Implementation of stress-majorization layout algorithm.
12 * \author Karsten Klein
15 * This file is part of the Open Graph Drawing Framework (OGDF).
19 * See README.txt in the root directory of the OGDF installation for details.
22 * This program is free software; you can redistribute it and/or
23 * modify it under the terms of the GNU General Public License
24 * Version 2 or 3 as published by the Free Software Foundation;
25 * see the file LICENSE.txt included in the packaging of this file
29 * This program is distributed in the hope that it will be useful,
30 * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 * GNU General Public License for more details.
35 * You should have received a copy of the GNU General Public
36 * License along with this program; if not, write to the Free
37 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
38 * Boston, MA 02110-1301, USA.
40 * \see http://www.gnu.org/copyleft/gpl.html
41 ***************************************************************/
43 #include <ogdf/energybased/StressMajorizationSimple.h>
46 #include <ogdf/energybased/FMMMLayout.h>
48 #include <ogdf/basic/SList.h>
51 #include <ogdf/basic/simple_graph_alg.h>
54 const double StressMajorization::startVal
= DBL_MAX
- 1.0;
55 const double StressMajorization::minVal
= DBL_MIN
;
56 const double StressMajorization::desMinLength
= 0.0001;
58 void StressMajorization::initialize(
60 const EdgeArray
<double>& eLength
,
61 NodeArray
< NodeArray
<double> >& oLength
,
62 NodeArray
< NodeArray
<double> >& weights
,
67 const Graph
&G
= GA
.constGraph();
69 m_prevEnergy
= startVal
;
70 m_prevLEnergy
= startVal
;
72 // all edges straight-line
77 //the shortest path lengths
78 forall_nodes (v
, G
) oLength
[v
].init(G
, DBL_MAX
);
79 forall_nodes (v
, G
) weights
[v
].init(G
, 0.0);
81 //-------------------------------------
82 //computes shortest path distances d_ij
83 //-------------------------------------
86 //we use simply BFS n times
87 //TODO experimentally compare speed, also with bintree dijkstra
92 maxDist
= allpairsspBFS(G
, oLength
, weights
);
94 timeUsed
= usedTime(timeUsed
);
95 cout
<< "\n******APSP BFS runtime: \n";
100 EdgeArray
<double> adaptedLength(G
);
101 adaptLengths(G
, GA
, eLength
, adaptedLength
);
102 //we use simply the BFM n times or Floyd instead, leading to cubic runtime
103 //TODO experimentally compare speed, also with bintree dijkstra
104 maxDist
= allpairssp(G
, adaptedLength
, oLength
, weights
, DBL_MAX
);
109 //TODO: Also allow different ways to compute centrality (closeness,...)
110 // and incorporate node sizes
112 computeRadii(G
, oLength
, maxDist
);
114 //Todo: Here we could add distance values depending on given node sizes
118 //uses closeness to compute the radii
119 void StressMajorization::computeRadii(const Graph
& G
, const NodeArray
< NodeArray
<double> >& distances
,
122 m_radii
.init(G
, 1.0);
123 //compute a center (simple version)
125 double minMax
= DBL_MAX
;
127 int numCentralNodes
= 0;
128 double maxCloseness
= 0.0;
129 double minCloseness
= DBL_MAX
;
131 // inverse of sum of shortest path distances
132 NodeArray
<double> closeness(G
, 0.0);
134 //Compute closeness values and min/max
137 double maxDist
= 0.0;
140 if (v
!= w
) closeness
[v
] += distances
[v
][w
];
141 if (distances
[v
][w
] > maxDist
) maxDist
= distances
[v
][w
];
143 if (maxDist
< minMax
)
148 closeness
[v
] = (G
.numberOfNodes()-1)/closeness
[v
]; //was 1/
150 //Check for min/max closeness
151 if (DIsGreater(closeness
[v
], maxCloseness
))
153 maxCloseness
= closeness
[v
];
156 else if (DIsEqual(closeness
[v
], maxCloseness
)) numCentralNodes
++;
157 if (DIsGreater(minCloseness
, closeness
[v
]))
158 minCloseness
= closeness
[v
];
161 //see paper by Brandes, Kenis, Wagner
162 double coincideOfs
= min(0.5, double(numCentralNodes
)/double((G
.numberOfNodes()-1)));
164 OGDF_ASSERT(center
!= 0);
167 //cout << "Diameter: "<<diameter<<"\n";
168 m_radii
[v
] = diameter
/double(2.0) * (1 - (closeness
[v
]-minCloseness
)/(maxCloseness
-minCloseness
+coincideOfs
));
169 //cout << "Radius: "<<m_radii[v]<<"\n";
175 void StressMajorization::mainStep(GraphAttributes
& GA
,
176 NodeArray
< NodeArray
<double> >& oLength
,
177 NodeArray
< NodeArray
<double> >& weights
,
178 const double maxDist
)
180 const Graph
&G
= GA
.constGraph();
183 //preparation for upward constraints
185 NodeArray
< NodeArray
<double> > upWeight(G
);
186 forall_nodes(v
, G
) upWeight
[v
].init(G
, 0.0);
187 //try to add upward constraints by looking at the adjacent edges
190 upWeight
[e
->source()][e
->target()] = 1.0;
191 upWeight
[e
->target()][e
->source()] = -1.0;
195 NodeArray
<int> nodeCounts(G
, 0);
197 NodeArray
<double> invPosNorm(G
);
198 NodeArray
< NodeArray
<double> > invDistb(G
); //bij in paper Pich Brandes
200 //this value is needed quite often
201 NodeArray
<double> wijSum(G
, 0.0); //sum of w_ij for each i over j!= i
209 if (v
!= w
) wijSum
[v
]+= weights
[v
][w
];
213 //--------------------------------------------------
214 // Iterations can either be done until a predefined
215 // number m_numSteps+m_itFac*|G| is reached or a break condition
217 //--------------------------------------------------
218 //for iteration checking
219 double weightsum
= 1.0;
220 if (m_upward
) weightsum
= 0.93;
222 OGDF_ASSERT(DIsGreater(m_numSteps
+m_itFac
*G
.numberOfNodes(), 0.0));
223 double kinv
= weightsum
/double(m_numSteps
+m_itFac
*G
.numberOfNodes());
225 NodeArray
<double> newX(G
);
226 NodeArray
<double> newY(G
);
227 //weighting for different optimization criteria
230 for (double tt
= 0; tt
<= weightsum
; tt
+= kinv
)
232 double iterStressSum
= 0.0;
233 if (m_radial
|| m_upward
) t
= tt
;
235 //compute inverse norm of position for radial distance constraints
240 //set values a_i as in paper pseudo code (simplified)
241 double tmp
= px
*px
+py
*py
;
242 if (DIsGreater(tmp
, 0.0))
243 invPosNorm
[v
] = 1/double(sqrt(tmp
));
244 else invPosNorm
[v
] = 0.0;
245 //cout << "invPosNorm: "<<invPosNorm[v]<<"\n";
247 //set values b_ij as in paper pseudo code (simplified)
251 double tmpij
= (px
-GA
.x(w
))*(px
-GA
.x(w
))+(py
-GA
.y(w
))*(py
-GA
.y(w
));
252 invDistb
[v
][w
] = (DIsGreater(tmpij
, 0.0) ? 1/double(sqrt(tmpij
)) : 0.0);
253 //cout << "invdistb: "<<invDistb[v][w]<<"\n";
260 //value corresponding to radial constraint
261 double radOfsX
= 0.0;
262 double radOfsY
= 0.0;
263 double tmpinvsq
= 0.0;
265 double upWeightSum
= 0.0;
269 upWeightSum
= 0.05+v
->degree()/100.0;
274 double tmp
= radius(v
);
275 OGDF_ASSERT(DIsGreater(tmp
, 0.0))
276 tmpinvsq
= 1/(tmp
*tmp
);
277 double tmppart
= t
*tmpinvsq
*tmp
*invPosNorm
[v
];
279 radOfsX
= tmppart
*GA
.x(v
);
280 radOfsY
= tmppart
*GA
.y(v
);
284 double upOfs
= 0.0; //only for y coordinate
286 //sum over all other nodes
288 double stressSumX
= 0.0;
289 double stressSumY
= 0.0;
294 stressSumX
+= weights
[v
][w
]*(GA
.x(w
)+oLength
[v
][w
]*(GA
.x(v
)-GA
.x(w
))*invDistb
[v
][w
]);
295 stressSumY
+= weights
[v
][w
]*(GA
.y(w
)+oLength
[v
][w
]*(GA
.y(v
)-GA
.y(w
))*invDistb
[v
][w
]);
296 //also add the influence of adjacent nodes that are not placed correctly
299 double val
= upWeight
[v
][w
];
300 //adjacent nodes influence each other
301 if (!DIsEqual(val
, 0.0))
303 if (DIsGreater(val
, 0.0)) //v is source
305 if (GA
.y(v
) > GA
.y(w
)-1.0)
307 upOfs
-= val
*(GA
.y(w
)-1.0*invDistb
[v
][w
]);
310 // else //v is target
312 // if (GA.y(w) > GA.y(v)-1.0)
314 // upOfs -= val*(GA.y(w)+(fabs(GA.y(v)-GA.y(w)))*invDistb[v][w]);
322 upOfs
*= upWeightSum
;
323 if (m_radial
|| m_upward
)
330 newX
[v
] = (stressSumX
+radOfsX
)/double(((1-t
)*wijSum
[v
]+t
*tmpinvsq
));
331 //experimental, only valid if disjoint constraints
332 newY
[v
] = (stressSumY
+radOfsY
+0.2*upOfs
)/double(((1-t
)*wijSum
[v
]+t
*tmpinvsq
+0.2*upWeightSum
));
335 iterStressSum
+= stressSumX
;
336 iterStressSum
+= stressSumY
;
338 //Alternatively, we can do the update after the computation step
339 // forall_nodes(v, G)
344 if (finished(iterStressSum
))
346 cout
<< iterStressSum
<<"\n";
349 }//outer loop: iterations or threshold
354 void StressMajorization::doCall(GraphAttributes
& GA
, const EdgeArray
<double>& eLength
, bool simpleBFS
)
356 const Graph
& G
= GA
.constGraph();
357 double maxDist
; //maximum distance between nodes
358 NodeArray
< NodeArray
<double> > oLength(G
);//first distance, then original length
359 NodeArray
< NodeArray
<double> > weights(G
);//standard weights as in MCGee,Kamada/Kawai
362 OGDF_ASSERT(isConnected(G
));
364 //compute relevant values
365 initialize(GA
, eLength
, oLength
, weights
, maxDist
, simpleBFS
);
367 //main loop with node movement
368 mainStep(GA
, oLength
, weights
, maxDist
);
370 if (simpleBFS
) scale(GA
);
374 void StressMajorization::call(GraphAttributes
& GA
)
376 const Graph
&G
= GA
.constGraph();
377 if(G
.numberOfEdges() < 1)
380 EdgeArray
<double> eLength(G
);//, 1.0);is not used
381 doCall(GA
, eLength
, true);
385 void StressMajorization::call(GraphAttributes
& GA
, const EdgeArray
<double>& eLength
)
387 const Graph
&G
= GA
.constGraph();
388 if(G
.numberOfEdges() < 1)
391 doCall(GA
, eLength
, false);
392 }//call with edge lengths
395 //changes given edge lengths (interpreted as weight factors)
396 //according to additional parameters like node size etc.
397 void StressMajorization::adaptLengths(
399 const GraphAttributes
& GA
,
400 const EdgeArray
<double>& eLengths
,
401 EdgeArray
<double>& adaptedLengths
)
403 //we use the edge lengths as factor and try to respect
404 //the node sizes such that each node has enough distance
406 //adapt to node sizes
409 double smax
= max(GA
.width(e
->source()), GA
.height(e
->source()));
410 double tmax
= max(GA
.width(e
->target()), GA
.height(e
->target()));
412 adaptedLengths
[e
] = (1+eLengths
[e
])*((smax
+tmax
));///2.0);
413 else adaptedLengths
[e
] = 5.0*eLengths
[e
];
418 void StressMajorization::shufflePositions(GraphAttributes
& GA
)
420 //random layout? FMMM? classical MDS? see Paper of Pich and Brandes
421 //just hope that we have low sigma values (distance error)
429 * Initialise the original estimates from nodes and edges.
432 //we could speed this up by not using nested NodeArrays and
433 //by not doing the fully symmetrical computation on undirected graphs
434 //All Pairs Shortest Path Floyd, initializes the whole matrix
435 //returns maximum distance. Does not detect negative cycles (lead to neg. values on diagonal)
436 //threshold is the value for the distance of non-adjacent nodes, distance has to be
438 // The weight parameter here is just for the stress majorization
439 // and directly set here for speedup
440 double StressMajorization::allpairssp(
442 const EdgeArray
<double>& eLengths
,
443 NodeArray
< NodeArray
<double> >& distance
,
444 NodeArray
< NodeArray
<double> >& weights
,
445 const double threshold
)
449 double maxDist
= -threshold
;
453 distance
[v
][v
] = 0.0f
;
454 weights
[v
][v
] = 0.0f
;
457 //TODO: Experimentally compare this against
458 // all nodes and incident edges (memory access) on huge graphs
461 distance
[e
->source()][e
->target()] = eLengths
[e
];
462 distance
[e
->target()][e
->source()] = eLengths
[e
];
466 // * And run the main loop of the algorithm.
475 if ((distance
[u
][v
] < threshold
) && (distance
[v
][w
] < threshold
))
477 distance
[u
][w
] = min( distance
[u
][w
], distance
[u
][v
] + distance
[v
][w
] );
478 weights
[u
][w
] = 1/(distance
[u
][w
]*distance
[u
][w
]);
479 //distance[w][u] = distance[u][w]; //is done anyway afterwards
481 if (distance
[u
][w
] < threshold
)
482 maxDist
= max(maxDist
,distance
[u
][w
]);
490 if (distance
[v
][v
] < 0.0) cerr
<< "\n###Error in shortest path computation###\n\n";
492 cout
<< "Maxdist: "<<maxDist
<<"\n";
497 cout
<< "Distance " << u
->index() << " -> "<<w
->index()<<" "<<distance
[u
][w
]<<"\n";
505 //the same without weights, i.e. all pairs shortest paths with BFS
507 //for compatibility, distances are double
508 // The weight parameter here is just for the stress majorization
509 // and directly set here for speedup
510 double StressMajorization::allpairsspBFS(
512 NodeArray
< NodeArray
<double> >& distance
,
513 NodeArray
< NodeArray
<double> >& weights
)
520 distance
[v
][v
] = 0.0f
;
525 //start in each node once
529 NodeArray
<bool> mark(G
, true);
536 node w
= bfs
.popFrontRet();
538 double d
= distance
[v
][w
]+1.0f
;
539 forall_adj_edges(e
,w
)
541 node u
= e
->opposite(w
);
547 weights
[v
][u
] = 1/(d
*d
);
548 maxDist
= max(maxDist
,d
);
555 //check for negative cycles
558 if (distance
[v
][v
] < 0.0) cerr
<< "\n###Error in shortest path computation###\n\n";
564 cout
<< "Maxdist: "<<maxDist
<<"\n";
569 cout
<< "Distance " << u
->index() << " -> "<<w
->index()<<" "<<distance
[u
][w
]<<"\n";
577 void StressMajorization::scale(GraphAttributes
& GA
)
579 //Simple version: Just scale to max needed
580 //We run over all nodes, find the largest distance needed and scale
581 //the edges uniformly
585 forall_edges(e
, GA
.constGraph())
587 double w1
= sqrt(GA
.width(e
->source())*GA
.width(e
->source())+
588 GA
.height(e
->source())*GA
.height(e
->source()));
589 double w2
= sqrt(GA
.width(e
->target())*GA
.width(e
->target())+
590 GA
.height(e
->target())*GA
.height(e
->target()));
591 w2
= (w1
+w2
)/2.0; //half length of both diagonals
592 double xdist
= GA
.x(e
->source())-GA
.x(e
->target());
593 double ydist
= GA
.y(e
->source())-GA
.y(e
->target());
594 double elength
= sqrt(xdist
*xdist
+ydist
*ydist
);
595 w2
= m_distFactor
* w2
/ elength
;//relative to edge length
602 forall_nodes(v
, GA
.constGraph())
604 GA
.x(v
) = GA
.x(v
)*maxFac
;
605 GA
.y(v
) = GA
.y(v
)*maxFac
;
608 cout
<< "Scaled by factor "<<maxFac
<<"\n";