1 /** \file
2 * \brief Implementation of the stress minimization via majorization
3 * algorithm.
4 *
5 * \author Mark Ortmann, University of Konstanz
6 *
7 * \par License:
8 * This file is part of the Open Graph Drawing Framework (OGDF).
9 *
10 * \par
11 * Copyright (C)<br>
12 * See README.md in the OGDF root directory for details.
13 *
14 * \par
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * Version 2 or 3 as published by the Free Software Foundation;
18 * see the file LICENSE.txt included in the packaging of this file
19 * for details.
20 *
21 * \par
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * \par
28 * You should have received a copy of the GNU General Public
29 * License along with this program; if not, see
30 * http://www.gnu.org/copyleft/gpl.html
31 */
32
33 #include <ogdf/energybased/StressMinimization.h>
34
35
36 namespace ogdf {
37
38 const double StressMinimization::EPSILON = 10e-4;
39
40 const int StressMinimization::DEFAULT_NUMBER_OF_PIVOTS = 50;
41
42
call(GraphAttributes & GA)43 void StressMinimization::call(GraphAttributes& GA)
44 {
45 const Graph& G = GA.constGraph();
46 // if the graph has at most one node nothing to do
47 if (G.numberOfNodes() <= 1) {
48 // make it exception save
49 for(node v : G.nodes)
50 {
51 GA.x(v) = 0;
52 GA.y(v) = 0;
53 }
54 return;
55 }
56 // Separate component layout cant be applied to a non-connected graph
57 OGDF_ASSERT(!m_componentLayout || isConnected(G));
58 NodeArray<NodeArray<double> > shortestPathMatrix(G);
59 NodeArray<NodeArray<double> > weightMatrix(G);
60 initMatrices(G, shortestPathMatrix, weightMatrix);
61 // if the edge costs are defined by the attribute copy it to an array and
62 // construct the proper shortest path matrix
63 if (m_hasEdgeCostsAttribute) {
64 OGDF_ASSERT(GA.has(GraphAttributes::edgeDoubleWeight));
65 m_avgEdgeCosts = dijkstra_SPAP(GA, shortestPathMatrix);
66 // compute shortest path all pairs
67 } else {
68 m_avgEdgeCosts = m_edgeCosts;
69 bfs_SPAP(G, shortestPathMatrix, m_edgeCosts);
70 }
71 call(GA, shortestPathMatrix, weightMatrix);
72 }
73
74
call(GraphAttributes & GA,NodeArray<NodeArray<double>> & shortestPathMatrix,NodeArray<NodeArray<double>> & weightMatrix)75 void StressMinimization::call(
76 GraphAttributes& GA,
77 NodeArray<NodeArray<double> >& shortestPathMatrix,
78 NodeArray<NodeArray<double> >& weightMatrix)
79 {
80 // compute the initial layout if necessary
81 if (!m_hasInitialLayout) {
82 computeInitialLayout(GA);
83 }
84 const Graph& G = GA.constGraph();
85 // replace infinity distances by sqrt(n) and compute weights.
86 // Note isConnected is only true during calls triggered by the
87 // ComponentSplitterLayout.
88 if (!m_componentLayout && !isConnected(G)) {
89 replaceInfinityDistances(shortestPathMatrix,
90 m_avgEdgeCosts * sqrt((double)(G.numberOfNodes())));
91 }
92 // calculate the weights
93 calcWeights(G, shortestPathMatrix, weightMatrix);
94 // minimize the stress
95 minimizeStress(GA, shortestPathMatrix, weightMatrix);
96 }
97
98
computeInitialLayout(GraphAttributes & GA)99 void StressMinimization::computeInitialLayout(GraphAttributes& GA)
100 {
101 PivotMDS* pivMDS = new PivotMDS();
102 pivMDS->setNumberOfPivots(DEFAULT_NUMBER_OF_PIVOTS);
103 pivMDS->useEdgeCostsAttribute(m_hasEdgeCostsAttribute);
104 pivMDS->setEdgeCosts(m_edgeCosts);
105 if (!m_componentLayout) {
106 // the graph might be disconnected therefore we need
107 // the component layouter
108 //design decision: should that parameter be passed to CSL?
109 ComponentSplitterLayout compLayouter;//(m_hasEdgeCostsAttribute);
110 compLayouter.setLayoutModule(pivMDS);
111 compLayouter.call(GA);
112 } else {
113 pivMDS->call(GA);
114 delete pivMDS;
115 }
116 }
117
118
replaceInfinityDistances(NodeArray<NodeArray<double>> & shortestPathMatrix,double newVal)119 void StressMinimization::replaceInfinityDistances(
120 NodeArray<NodeArray<double> >& shortestPathMatrix,
121 double newVal)
122 {
123 const Graph &G = *shortestPathMatrix.graphOf();
124
125 for(node v : G.nodes) {
126 for(node w : G.nodes) {
127 if (v != w && isinf(shortestPathMatrix[v][w])) {
128 shortestPathMatrix[v][w] = newVal;
129 }
130 }
131 }
132 }
133
134
calcWeights(const Graph & G,NodeArray<NodeArray<double>> & shortestPathMatrix,NodeArray<NodeArray<double>> & weightMatrix)135 void StressMinimization::calcWeights(
136 const Graph& G,
137 NodeArray<NodeArray<double> >& shortestPathMatrix,
138 NodeArray<NodeArray<double> >& weightMatrix)
139 {
140 for (node v : G.nodes) {
141 for (node w : G.nodes) {
142 if (v != w) {
143 // w_ij = d_ij^-2
144 weightMatrix[v][w] = 1
145 / (shortestPathMatrix[v][w] * shortestPathMatrix[v][w]);
146 }
147 }
148 }
149 }
150
151
calcStress(const GraphAttributes & GA,NodeArray<NodeArray<double>> & shortestPathMatrix,NodeArray<NodeArray<double>> & weightMatrix)152 double StressMinimization::calcStress(
153 const GraphAttributes& GA,
154 NodeArray<NodeArray<double> >& shortestPathMatrix,
155 NodeArray<NodeArray<double> >& weightMatrix)
156 {
157 double stress = 0;
158 for (node v = GA.constGraph().firstNode(); v != nullptr; v = v->succ()) {
159 for (node w = v->succ(); w != nullptr; w = w->succ()) {
160 double xDiff = GA.x(v) - GA.x(w);
161 double yDiff = GA.y(v) - GA.y(w);
162 double zDiff = 0.0;
163 if (GA.has(GraphAttributes::threeD))
164 {
165 zDiff = GA.z(v) - GA.z(w);
166 }
167 double dist = sqrt(xDiff * xDiff + yDiff * yDiff + zDiff * zDiff);
168 if (dist != 0) {
169 stress += weightMatrix[v][w] * (shortestPathMatrix[v][w] - dist)
170 * (shortestPathMatrix[v][w] - dist);//
171 }
172 }
173 }
174 return stress;
175 }
176
177
copyLayout(const GraphAttributes & GA,NodeArray<double> & newX,NodeArray<double> & newY)178 void StressMinimization::copyLayout(
179 const GraphAttributes& GA,
180 NodeArray<double>& newX,
181 NodeArray<double>& newY)
182 {
183 // copy the layout
184 for(node v : GA.constGraph().nodes)
185 {
186 newX[v] = GA.x(v);
187 newY[v] = GA.y(v);
188 }
189 }
190
191
copyLayout(const GraphAttributes & GA,NodeArray<double> & newX,NodeArray<double> & newY,NodeArray<double> & newZ)192 void StressMinimization::copyLayout(
193 const GraphAttributes& GA,
194 NodeArray<double>& newX,
195 NodeArray<double>& newY,
196 NodeArray<double>& newZ)
197 {
198 // copy the layout
199 for(node v : GA.constGraph().nodes)
200 {
201 newX[v] = GA.x(v);
202 newY[v] = GA.y(v);
203 newZ[v] = GA.z(v);
204 }
205 }
206
207
minimizeStress(GraphAttributes & GA,NodeArray<NodeArray<double>> & shortestPathMatrix,NodeArray<NodeArray<double>> & weightMatrix)208 void StressMinimization::minimizeStress(
209 GraphAttributes& GA,
210 NodeArray<NodeArray<double> >& shortestPathMatrix,
211 NodeArray<NodeArray<double> >& weightMatrix)
212 {
213 const Graph& G = GA.constGraph();
214 int numberOfPerformedIterations = 0;
215
216 double prevStress = std::numeric_limits<double>::max();
217 double curStress = std::numeric_limits<double>::max();
218
219 if (m_terminationCriterion == TerminationCriterion::Stress) {
220 curStress = calcStress(GA, shortestPathMatrix, weightMatrix);
221 }
222
223 NodeArray<double> newX;
224 NodeArray<double> newY;
225 NodeArray<double> newZ;
226
227 if (m_terminationCriterion == TerminationCriterion::PositionDifference) {
228 newX.init(G);
229 newY.init(G);
230 if (GA.has(GraphAttributes::threeD))
231 newZ.init(G);
232 }
233 do {
234 if (m_terminationCriterion == TerminationCriterion::PositionDifference) {
235 if (GA.has(GraphAttributes::threeD))
236 copyLayout(GA, newX, newY, newZ);
237 else copyLayout(GA, newX, newY);
238 }
239 nextIteration(GA, shortestPathMatrix, weightMatrix);
240 if (m_terminationCriterion == TerminationCriterion::Stress) {
241 prevStress = curStress;
242 curStress = calcStress(GA, shortestPathMatrix, weightMatrix);
243 }
244 } while (!finished(GA, ++numberOfPerformedIterations, newX, newY, prevStress, curStress));
245
246 Logger::slout() << "Iteration count:\t" << numberOfPerformedIterations
247 << "\tStress:\t" << calcStress(GA, shortestPathMatrix, weightMatrix) << std::endl;
248 }
249
250
nextIteration(GraphAttributes & GA,NodeArray<NodeArray<double>> & shortestPathMatrix,NodeArray<NodeArray<double>> & weights)251 void StressMinimization::nextIteration(
252 GraphAttributes& GA,
253 NodeArray<NodeArray<double> >& shortestPathMatrix,
254 NodeArray<NodeArray<double> >& weights)
255 {
256 const Graph& G = GA.constGraph();
257
258 for (node v : G.nodes)
259 {
260 double newXCoord = 0.0;
261 double newYCoord = 0.0;
262 double newZCoord = 0.0;
263 double& currXCoord = GA.x(v);
264 double& currYCoord = GA.y(v);
265 double totalWeight = 0;
266 for (node w : G.nodes)
267 {
268 if (v == w) {
269 continue;
270 }
271 // calculate euclidean distance between both points
272 double xDiff = currXCoord - GA.x(w);
273 double yDiff = currYCoord - GA.y(w);
274 double zDiff = (GA.has(GraphAttributes::threeD)) ? GA.z(v) - GA.z(w) : 0.0;
275 double euclideanDist = sqrt(xDiff * xDiff + yDiff * yDiff + zDiff * zDiff);
276 // get the weight
277 double weight = weights[v][w];
278 // get the desired distance
279 double desDistance = shortestPathMatrix[v][w];
280 // reset the voted x coordinate
281 // if x is not fixed
282 if (!m_fixXCoords) {
283 double voteX = GA.x(w);
284 if (euclideanDist != 0) {
285 // calc the vote
286 voteX += desDistance * (currXCoord - voteX) / euclideanDist;
287 }
288 // add the vote
289 newXCoord += weight * voteX;
290 }
291 // reset the voted y coordinate
292 // y is not fixed
293 if (!m_fixYCoords) {
294 double voteY = GA.y(w);
295 if (euclideanDist != 0) {
296 // calc the vote
297 voteY += desDistance * (currYCoord - voteY) / euclideanDist;
298 }
299 newYCoord += weight * voteY;
300 }
301 if (GA.has(GraphAttributes::threeD) && !m_fixZCoords) {
302 // reset the voted z coordinate
303 double voteZ = GA.z(w);
304 if (euclideanDist != 0) {
305 // calc the vote
306 voteZ += desDistance * (GA.z(v) - voteZ) / euclideanDist;
307 }
308 newZCoord += weight * voteZ;
309 }
310 // sum up the weights
311 totalWeight += weight;
312 }
313 // update the positions
314 if (totalWeight != 0) {
315 if (!m_fixXCoords) {
316 currXCoord = newXCoord / totalWeight;
317 }
318 if (!m_fixYCoords) {
319 currYCoord = newYCoord / totalWeight;
320 }
321 if (GA.has(GraphAttributes::threeD) && !m_fixZCoords) {
322 GA.z(v) = newZCoord / totalWeight;
323 }
324 }
325 }
326 }
327
328
finished(GraphAttributes & GA,int numberOfPerformedIterations,NodeArray<double> & prevXCoords,NodeArray<double> & prevYCoords,const double prevStress,const double curStress)329 bool StressMinimization::finished(
330 GraphAttributes& GA,
331 int numberOfPerformedIterations,
332 NodeArray<double>& prevXCoords,
333 NodeArray<double>& prevYCoords,
334 const double prevStress,
335 const double curStress)
336 {
337 if (numberOfPerformedIterations == m_numberOfIterations) {
338 return true;
339 }
340
341 switch (m_terminationCriterion)
342 {
343 case TerminationCriterion::PositionDifference:
344 {
345 double eucNorm = 0;
346 double dividend = 0;
347 // compute the translation of all nodes between
348 // the consecutive layouts
349 for (node v : GA.constGraph().nodes)
350 {
351 double diffX = prevXCoords[v] - GA.x(v);
352 double diffY = prevYCoords[v] - GA.y(v);
353 dividend += diffX * diffX + diffY * diffY;
354 eucNorm += prevXCoords[v] * prevXCoords[v] + prevYCoords[v] * prevYCoords[v];
355 }
356 return sqrt(dividend) / sqrt(eucNorm) < EPSILON;
357 }
358 case TerminationCriterion::Stress:
359 return curStress == 0 || prevStress - curStress < prevStress * EPSILON;
360
361 default:
362 return false;
363 }
364 }
365
366
initMatrices(const Graph & G,NodeArray<NodeArray<double>> & shortestPathMatrix,NodeArray<NodeArray<double>> & weightMatrix)367 void StressMinimization::initMatrices(const Graph& G,
368 NodeArray<NodeArray<double> >& shortestPathMatrix,
369 NodeArray<NodeArray<double> >& weightMatrix)
370 {
371 // init shortest path matrix by infinity distances
372 for (node v : G.nodes)
373 {
374 shortestPathMatrix[v].init(G, std::numeric_limits<double>::infinity());
375 shortestPathMatrix[v][v] = 0;
376 weightMatrix[v].init(G, 0);
377 }
378 }
379
380 }
381