Unified graph data type and edge directions

[lindnemi]

This is the state of my thought on edge directions:

Mathematically speaking, an undirected graph is equivalent to a directed graph if for every directed edge i => j there is another edge j => i.

NetworkDynamics.jl graph data structures are derived from LightGraphs.jl, which implements SimpleGraph and SimpleDiGraph differently. For a SimpleGraph edges(g) returns a list of ordered Pairs i => j where i is smaller or equal j. For a directed Graph every Pair of indices in edges(g) represents a single edge with the ordering of i and j corresponding to the direction of that edge.

Thus, if one has an undirected graph g every edge will be saved only once if the undirected graph is represented as a SimpleGraph, but twice (with different direction) if the graph is represented as a SimpleDiGraph.

At the moment network_dynamics works equally well with SimpleGraph and SimpleDiGraph. However, due to the differences in the LightGraphs data structures it will not behave the same way. For an undirected network every coupling function might be evaluated once or twice depending on the data structure chosen as input. This is far from obvious to the user and happens without warning.

From a modelist point of view the SimpleDiGraph data structre is desirable since it is closer to how one would mathematically specify the differential equations defining a complex network. However, SimpleGraph offers a slight performance advantage since for symmetric coupling functions on undirected graphs it is possible to reuse edge values. This is why ND is based on SimpleGraph at the moment.

I would like to make SimpleDiGraph the new default, since the functions for edges and coupling terms are then easier to write and there are less complications due to the fiducial orientation.

My first idea was to convert any graph input by default into a SimpleDiGraph[1]. However, this leads to the issue that a SimpleGraph with ne(g) = N might be converted into a SimpleDiGraph with ne(g) = 2N. In this case the input Arrays of EdgeFunctions, edge parameters, and even edge initial conditions(for ODEEdges) would have to be adjusted as well. While for Edge Functions and parameters automatic conversion is imaginable, for edge initial conditions it would be unclear how to handle the discrepancy of dimension.

At the moment I am quite puzzled on how to untie all those intertwined requirements.

One idea is to have a layer structure where each layer consists of a graph topology and a single edge function. However, this might require a new way to specify edge parameters, since then callable structs carrying their individual parameters are no longer possible.

Another idea is to implement our own graph data type as a subtye of AbstractGraph and introduce a fixed index for every edge. Then we could overload edges(g) as well as ne(g) in order to have a unified interface and to de-couple performance optimization from the interface.

---

[1] except if the user specifies an additional keyword orientation=:fiducial and thereby asserts that she wants to user the harder to write but slightly more performant SimpleGraph version of ND.

[lindnemi]

I guess that in dynamical system on networks the idea of an "undirected edge" is an illusion anyways. Even when you write the equations like

the edge from i to j appears twice, once as g(i,j) at node i and once as g(j,i) and node j. The only thing that you might have are symmetric coupling functions. If the user gives us a hint that her EdgeFunctions are symmetric we might perform some index magic under the hood and save computations even when everything is based on a SimpleDiGraph-like data structure.

[hexaeder]

As I far as understand there are three types of edges: symmetric, asymmetric and anti-symmetric.

# symmetric means
vᵢ = ... + f(vᵢ, vⱼ) + ...
vⱼ = ... + f(vᵢ, vⱼ) + ...
# antisymmetric means
vᵢ = ... + f(vᵢ, vⱼ) + ...
vⱼ = ... - f(vᵢ, vⱼ) + ...
# asymmetric means
vᵢ = ... + h(vᵢ, vⱼ) + ...
vⱼ = ... + g(vⱼ, vᵢ) + ...

Symmetric edges are easy: there is no concept of src and dst. They are just connecting two nodes. In a network with only symmetric edges one could introduce something like a symmetric vertex function which stacks the outgoing and incoming edges on to of each other and call it edges. It is okay to use a SimpleGraph for this because you never exploit the underlying direction of the implementation. It would be also viable to use a SimpleDiGraph and define f1 for i=>j and f2 for j=>. In this case the vertices should only pay attention either to the outgoing or incoming edges (which one is up to convention, from now on in choose the incoming).

For completely asymmetric edges it is easy as well, one has to implement h and g separate on directed edges using a DiGraph. In this case one would again ignore all the outgoing edges in the vertex.

The antisymmetric version makes everything complicated. One could just treat them like asymmetric terms, splitting them up into two parts and only care about either the outgoing or incoming ones in the vertex function. Another approach (which is used in the docs) is to define a antisymmetric function and resolve this antisymmetry through the concept of outgoing and incoming edges. I think on should never use this approach together with SimpleGraph because it abuses the type. The underlying direction of a SimpleGraph is due to implimentation not due to concept and shouldn't be exploited in algorithms. As a first countermeasure one could check all the edges for symmetry if defined on a SimpleGraph and print out a warning if there is any.

The problem seems to concentrate around the vertex functions. There are three corresponding implementations

function f1!(dx, x, out_edges, in_edges) # weird but useful for performant antisymmetric edges
function f2!(dx, x, in_edges) # cleanest for asymmetric edges
function f3!(dx, x, all_edges) # cleanest for symmetric edges

whereat the first one is the one actually used. Personally I think it is still the nicest. It offers the most control over the different scenarios and allow for combinations of edge types. I guess the idea for a multilayer network would be to implement multiple graphs on the same nodes. For each layer one could specialize on the implementation for the edge type and the absolute vertex value will be the sum over all layers? I think this is a very interesting concept but probably even more confusing for the user than the status quo?

Another idea is to implement our own graph data type as a subtye of AbstractGraph and introduce a fixed index for every edge. Then we could overload edges(g) as well as ne(g) in order to have a unified interface and to de-couple performance optimization from the interface.

This is a good idea anyways. The package internals need strict indexing throughout the edges. Right know this indexing is given by the sorting of the edges(g::AbstractGraph iterator. This confused me on the first tries with the package since the ordering is not correlated to the order of add_edge! calls.

[lindnemi]

We should leave the vertex function as they are. In some scenarios one might want to couple via the out-edges and in other scenarios via the incoming edges.

What is strange is that e_s, e_d is not consistent with what it claims to be. If a vertex i is a source of an edge, then i would expect the edge to be i => j and vice versa for destination. On the other hand, i would expect all out-edges to be contained in e_s. However that's not the case. Some out-edges are in e_s and others are in e_d just with the opposite sign.

I just had the idea to add a new field sign to the EdgeData views. If one uses an anti-symmetric coupling function on an undirected SimpleGraph, every EdgeData view would be constructed twice with opposite signs. One could write the same functions as for SimpleDiGraph and couple only via the in-coming edges, but under the hood the actual edge values would be computed only once.

[lindnemi]

I noticed something strange:

As a physicist I am interested in the coupling sum

x_i = f_i + \sum_j A_{ij} g_{ij}

The term A_ij is 1 if there is an edge i => j. That is g_ij computes the flow on the edge from i to j.

However, if I take the perspective of the edge, the flow going from i to j should arrive at the destination vertex j and not at the source vertex i.

EDIT: The following suggestion is not at all more consistent.

---

To be more consistent one could define h_ij = g_ji and write equivalently

x_i = f_i + \sum_j A_{ij} h_{ji}

[lindnemi]

Here is a suggestion for static half-edges. Tomorrow I want to ry to implement it:

There are four types of edges: directed, symmetric, antisymmetric, asymmetric

The user specifies every edge as if it was directed, including its dimension, e.g.

function edge!(e, v_s, v_d, p ,t)
    sin(v_d[1] - v_s[1])
end

ef! = StaticEdge(f!=edge!, dim = 1)

The user may specify a hint if the edge function is symmetric or antisymmetric (we might even @assert that it really is)

Then a dispatch happens when the graph is constructed.

For a SimpleDiGraph everything stays as it is.

For a SimpleGraph the EdgeFunction will be reconstructed. First the dimension is doubled. Then a wrapper for f! is defined. If the Edge Function is...

• ...symmetric the wrapper is: f!(e[1]); e[2] = e[1]
• ...antisymmetric the wrapper is: f!(e[1]); e[2] = -e[1]
• ...asymmetric the wrapper is: f!(e[1], v_s, v_d); f!(e[2], v_d,v_s)

Everything else should stay more or less as it is. Probably the indices for the views need to be adjusted, since the vertex function will only see one array of (incoming) edges, instead of two arrays. I didn't spell that part out in detail.

The only case this doesn't cover in the most efficient manner yet is antisymmetric ODE edges. But i have never heard of them anyways. An orthogonal issue are edges that have internal algebraic equations which could be parsed as expressions instead of as ODE variables. I haven't thought about that any more, the relevant prototype is here: https://github.com/Lima-Lima/NetworkDynamics.jl/commit/5f4ea9115dd33292e47a1af277a5a4a5f10037ce