Feat request: integrate `by_axis_local_stress_graph` network layout from GraphRecipes about graphmakie.jl HOT 4 CLOSED

using Pkg
Pkg.activate(temp=true)
Pkg.add("GLMakie")
Pkg.add("LightGraphs")
Pkg.add("GraphMakie")

using GLMakie
using LightGraphs
using GraphMakie
using GraphMakie.NetworkLayout
using SparseArrays

adj_matrix = [0 1 1 0 0 0 0 0 0 0;
              0 0 0 0 1 1 0 0 0 0;
              0 0 0 1 0 0 1 0 1 0;
              0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 1 0 1;
              0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0;
              0 0 0 0 0 0 0 0 0 0]

g = SimpleGraph(SimpleDiGraph(adj_matrix))
graphplot(g; layout=Stress())

# code copied and adapted from GraphRecipes
# https://github.com/JuliaPlots/GraphRecipes.jl/blob/5559436a981089eb842f4c7b7a96a4dbf9abe466/src/graph_layouts.jl#L157-L204
function by_axis_local_stress_graph(g)
    adjmat = adjacency_matrix(g)

    # parameters
    node_weights = ones(size(adjmat,1))
    dim = 2
    free_dims = 1:dim
    x = rand(length(node_weights))
    y = rand(length(node_weights))
    z = rand(length(node_weights))
    maxiter = 1000

    n = length(node_weights)

    # graph-theoretical distance between node i and j (i.e. shortest path distance)
    # TODO: calculate a real distance
    dist = estimate_distance(adjmat)
    # @show dist

    # also known as kᵢⱼ in "axis-by-axis stress minimization".  the -2 could also be 0 or -1?
    w = dist .^ -2

    # in each iteration, we update one dimension/node at a time, reducing the total stress with each update
    X = dim == 2 ? (x, y) : (x, y, z)
    laststress = stressf(X, dist, w)
    for k in 1:maxiter
        for p in free_dims
            for i=1:n
                numer, denom = 0.0, 0.0
                for j=1:n
                    i==j && continue
                    numer += w[i,j] * (X[p][j] + dist[i,j] * (X[p][i] - X[p][j]) / norm_ij(X, i, j))
                    denom += w[i,j]
                end
                if denom != 0
                    X[p][i] = numer / denom
                end
            end
        end

        # check for convergence of the total stress
        thisstress = stressf(X, dist, w)
        if abs(thisstress - laststress) / abs(laststress) < 1e-6
            # info("converged. numiter=$k last=$laststress this=$thisstress")
            break
        end
        laststress = thisstress
    end

    GraphMakie.Pointf.(zip(X...))
end

function estimate_distance(adjmat::AbstractMatrix)
    g = LightGraphs.Graph(adjmat)
    dists = convert(Matrix{Float64}, hcat(map(i->LightGraphs.dijkstra_shortest_paths(g, i).dists, LightGraphs.vertices(g))...))
    tot = 0.0; cnt = 0
    for (i,d) in enumerate(dists)
        if d < 1e10
            tot += d
            cnt += 1
        end
    end
    avg = cnt > 0 ? tot / cnt : 1.0
    for (i,d) in enumerate(dists)
        if d > 1e10
            dists[i] = 3avg
        end
    end
    dists
end

norm_ij(X, i, j) = sqrt(sum(Float64[(v[i]-v[j])^2 for v in X]))
stressf(X, dist, w, i, j) = w[i,j] * (norm_ij(X, i, j) - dist[i,j])^2
function stressf(X, dist, w)
    tot = 0.0
    for i=1:size(X,1), j=1:i-1
        tot += stressf(X, dist, w, i, j)
    end
    tot
end

graphplot(g; layout=by_axis_local_stress_graph)