bare-bones stress majorization in julia

main.jl

# 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