-
Notifications
You must be signed in to change notification settings - Fork 22
/
stress.jl
216 lines (183 loc) · 6.43 KB
/
stress.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
"""
Compute graph layout using stress majorization
Inputs:
δ: Matrix of pairwise distances
p: Dimension of embedding (default: 2)
weights: Matrix of weights. If not specified, defaults to
weights[i,j] = δ[i,j]^-2 if δ[i,j] is nonzero, or 0 otherwise
X0: Initial guess for the layout. Coordinates are given in rows.
If not specified, default to random matrix of Gaussians
Additional optional keyword arguments control the convergence of the algorithm
and the additional output as requested:
iterations: Maximum number of iterations. Default: 400size(X0, 1)^2
abstols: Absolute tolerance for convergence of stress.
The iterations terminate if the difference between two
successive stresses is less than abstol.
Default: √(eps(eltype(X0))
reltols: Relative tolerance for convergence of stress.
The iterations terminate if the difference between two
successive stresses relative to the current stress is less than
reltol. Default: √(eps(eltype(X0))
abstolx: Absolute tolerance for convergence of layout.
The iterations terminate if the Frobenius norm of two successive
layouts is less than abstolx. Default: √(eps(eltype(X0))
Output:
The final layout positions.
Reference:
The main equation to solve is (8) of:
@incollection{
author = {Emden R Gansner and Yehuda Koren and Stephen North},
title = {Graph Drawing by Stress Majorization}
year={2005},
isbn={978-3-540-24528-5},
booktitle={Graph Drawing},
seriesvolume={3383},
series={Lecture Notes in Computer Science},
editor={Pach, J\'anos},
doi={10.1007/978-3-540-31843-9_25},
publisher={Springer Berlin Heidelberg},
pages={239--250},
}
"""
module Stress
using GeometryTypes, Compat, StaticArrays
import Base: start, next, done, *
function (*){T<:LinAlg.BlasFloat,S<:StaticArray}(A::StridedMatrix{T}, x::StridedVector{S})
A_mul_B!(similar(x, S, size(A,1)), A, x)
end
immutable Layout{M1<:AbstractMatrix, M2<:AbstractMatrix, VP<:AbstractVector, FT<:AbstractFloat}
δ::M1
weights::M2
positions::VP
pinvLw::Matrix{Float64}
iterations::Int
abstols::FT
reltols::FT
abstolx::FT
end
function initialweights(D)
map(D) do d
x = Float64(d^(-2.0))
isfinite(x) ? x : zero(Float64)
end
end
function Layout{N, T}(
δ, PT::Type{Point{N, T}}=Point{2, Float64};
startpositions=rand(PT, size(δ,1)), weights=initialweights(δ),
iterations=400*size(δ,1)^2, abstols=√(eps(T)),
reltols=√(eps(T)), abstolx=√(eps(T))
)
@assert size(startpositions, 1)==size(δ, 1)==size(δ, 2)==size(weights, 1)==size(weights, 2)
Lw = weightedlaplacian(weights)
pinvLw = pinv(Lw)
return Layout(δ, weights, startpositions, pinvLw, iterations, abstols, reltols, abstolx)
end
layout(δ, dim::Int; kw_args...) = layout(δ, Point{dim,Float64}; kw_args...)
function layout{N, T}(
δ, PT::Type{Point{N, T}}=Point{2, Float64};
startpositions=rand(PT, size(δ,1)), kw_args...
)
layout!(δ, startpositions; kw_args...)
end
function layout!{N, T}(
δ, startpositions::AbstractVector{Point{N, T}};
iterations=400*size(δ,1)^2, kw_args...
)
iter = Layout(δ, Point{N,T}; startpositions=startpositions, kw_args...)
state = start(iter)
while !done(iter, state)
iter, state = next(iter, state)
end
state[end] > iterations && warn("Maximum number of iterations reached without convergence")
iter.positions
end
function start(network::Layout)
s = stress(network.positions, network.δ, network.weights)
return (s, s, network.positions, 0)
end
function next(iter::Layout, state)
newstress, oldstress, X0, i = state
δ, weights, pinvLw, positions, X0 = iter.δ, iter.weights, iter.pinvLw, iter.positions, copy(iter.positions)
#TODO the faster way is to drop the first row and col from the iteration
t = LZ(X0, δ, weights)
positions = pinvLw * (t*X0)
@assert all(x->all(map(isfinite, x)), positions)
newstress, oldstress = stress(positions, δ, weights), newstress
iter.positions[:] = positions
return iter, (newstress, oldstress, X0, (i+1))
end
function done(iter::Layout, state)
newstress, oldstress, X0, i = state
iterations, reltols = iter.iterations, iter.reltols
positions, abstols, abstolx = iter.positions, iter.abstols, iter.abstolx
(i == 0) && (i<iterations) && return false #special case 0th iteration
return (
i > iterations ||
abs(newstress - oldstress) < reltols * newstress ||
abs(newstress - oldstress) < abstols ||
vecnorm(positions - X0) < abstolx
)
end
"""
Stress function to majorize
Input:
positions: A particular layout (coordinates in rows)
d: Matrix of pairwise distances
weights: Weights for each pairwise distance
See (1) of Reference
"""
function stress{T, N}(
positions::AbstractArray{Point{T, N}},
d=ones(T, length(positions), length(positions)), weights=initialweights(d, T)
)
s = zero(T); n = length(positions)
@assert n==size(d, 1)==size(d, 2)==size(weights, 1)==size(weights, 2)
for j=1:n, i=1:j-1
s += weights[i, j] * (norm(positions[i] - positions[j]) - d[i,j])^2
end
@assert isfinite(s)
s
end
"""
Compute weighted Laplacian given ideal weights weights
Lʷ defined in (4) of the Reference
"""
function weightedlaplacian{T}(weights::AbstractMatrix{T})
n = Compat.LinAlg.checksquare(weights)
Lw = zeros(T, n, n)
for i=1:n
D = zero(T)
for j=1:n
i==j && continue
Lw[i, j] = -weights[i, j]
D += weights[i, j]
end
Lw[i, i] = D
end
Lw
end
"""
Computes L^Z defined in (5) of the Reference
Input: Z: current layout (coordinates)
d: Ideal distances (default: all 1)
weights: weights (default: d.^-2)
"""
function LZ{N,T}(Z::AbstractVector{Point{N,T}}, d, weights)
n = length(Z)
L = zeros(T, n, n)
for i=1:n
D = zero(T)
for j=1:n
i==j && continue
nrmz = norm(Z[i] - Z[j])
nrmz==0 && continue
δ = weights[i, j] * d[i, j]
L[i, j] = -δ/nrmz
D -= -δ/nrmz
end
@assert isfinite(D)
L[i, i] = D
end
L
end
end # end of module