Support self-loops

[lindnemi]

Maybe by taking adjacency inputs?[However explicit self loops might pose some issues if they accidentally get looped over twice) Think about it first.

At the moment the best way is to pull the self-interaction into the Vertex Function, e.g.

[...]
for e in edges
  dv[2] += e[1]
end
dv[2] += sin(v[1] - v[1] + alpha) # placeholder for a self coupling term
[...]

[hexaeder]

Is there a problem with self edges currently? If so we should check the constructors and add errors. On first glance it seems to work just fine and the graph structure looks about right. Even though testing self loops with diffusion edges is kinda pointless.

using NetworkDynamics
using Graphs
using OrdinaryDiffEq
using Plots

function diffusionedge!(e, v_s, v_d, p, t)
    e .= v_s - v_d
    nothing
end

function diffusionvertex!(dv, v, edges, p, t)
    dv .= 0.0
    for e in edges
        dv .+= e
    end
    nothing
end

g = SimpleGraph(2)
add_edge!(g, 1, 2)
add_edge!(g, 1, 1)
add_edge!(g, 2, 2)

nd_diffusion_vertex = ODEVertex(; f=diffusionvertex!, dim=1)
nd_diffusion_edge = StaticEdge(; f=diffusionedge!, dim=1)

nd = network_dynamics(nd_diffusion_vertex, nd_diffusion_edge, g)

x0 = randn(2) # random initial conditions
ode_prob = ODEProblem(nd, x0, (0.0, 4.0))
sol = solve(ode_prob, Tsit5());

plot(sol)

structure = nd.f.graph_structure
@test structure.num_e == 3
@test structure.e_dims == [2, 2, 2]
@test structure.e_offs == [0, 2, 4]

[lindnemi]

The problem is: For an undirected coupling every edge is evaluated twics: s => d and d => s. This means a self loop would have twice the coupling strength of any other edge. An easy way out might be to allow self loops only in directed graphs?

[lindnemi]

And, tbh, i thought SimpleGraphs does not support self-loops

[hexaeder]

Interesting enough, thats how Graphs inteprets the adj matrix of an undirected graph with a self edge anyway (the double counting).

julia> g = SimpleGraph(2)
julia> add_edge!(g,1,2)
julia> add_edge!(g,1,1)
julia> adjacency_matrix(g)
2×2 SparseArrays.SparseMatrixCSC{Int64, Int64} with 3 stored entries:
 2  1
 1  ⋅

Well, i'd still say in our concept of graphs there really are no undirected edges, only pairs of directed edges going in both directions (which is reflected by the symmetry of the adj matrix in traditional matrix formulation of the network dynamics), a fact which is some what hidden for the user right now because we do this transformation in reconstruct_edge.... i think at some point this should be reflected in the public interface... but I digress :D