Skip to content

Instantly share code, notes, and snippets.

@cscheid
Created Sep 10, 2014
Embed
What would you like to do?
bare-bones stress majorization in julia
# You should totally ignore this because it's the first piece of Julia I've ever written.
# 2-clause BSD, blah.
function make_cycle(n)
result = Dict{Int32, Array{Int32}}()
for i = 1:n
ii = i - 1
push!(result, i, [1 + ((i+n-2) % n), 1 + i % n])
end
result
end
function apsp_floyd_warshall(g)
n = length(g)
result = Array(Int32, n, n)
fill!(result, 1000000)
for i in 1:n
result[i,i] = 0
end
for (from, to_list) in g
for to in to_list
result[from, to] = 1
end
end
for i in 1:n
for j in 1:n
for k in 1:n
result[i,j] = min(result[i,k] + result[k,j], result[i,j])
end
end
end
result
end
function make_laplacian!(m)
n = size(m, 1)
for i in 1:n
m[i, i] = 0.0
end
for i in 1:n
for j in 1:n
if i != j
m[i, i] -= m[i, j]
end
end
end
m
end
function stress_majorization(g, dim = 2, max_iter = 100)
n = length(g)
d = apsp_floyd_warshall(g)
w = make_laplacian!(-map(x -> if x == 0 0.0 else 1.0/(x ^ 2) end, d))
function pseudo_inv_diag(v)
[if x > 1e-10 1/x else x end for x in v]
end
svd = svdfact(w)
# setting w_ij = 1/d_ij^2, delta_ij := w_ij d_ij = 1/d_ij
delta = map(x -> if x == 0 0.0 else 1.0/x end, d)
# Horribly inefficient
function invert_vector(v)
svd[:Vt]' * diagm(pseudo_inv_diag(svd[:S])) * svd[:U]' * v
end
result = rand(n, dim)
function dist(v1, v2)
d = sum((v1 - v2) .^ 2)
if (d > 0)
d ^ -0.5
else
0
end
end
relative_movement = 1
iter = 0
while relative_movement > 1e-4 && iter < max_iter
Lz = make_laplacian!(Float32[ - delta[i,j] * dist(result[i,:], result[j,:]) for i in 1:n, j in 1:n ])
new_result = invert_vector(Lz * result)
tl = sum(result .^ 2) # total_length
tdl = sum((result - new_result) .^ 2) # total_delta_length
relative_movement = (tdl / tl) ^ 0.5 / (n * dim)
iter += 1
println("relative change: $relative_movement, iteration $iter")
result = new_result
end
result
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment