NetworkSystem

[ChrisRackauckas]

It would be nice to have the following work on any AbstractSystem. You supply a network given by a LightGraphs AbstractGraph, you provide a prototype component, and you have a syntax for unidirectional edges. Here's an example. Let g be the graph and suppose we've already defined a component c. Then we might want to construct:

c = ODESystem(...)
connection(cfrom,cto) = [cfrom.x + cto.x ~ 0,cfrom.y - cto.y ~ 0]
ODESystem(connection,g,c)

that clones c into c_i many versions, one for each node, and for every edge makes that connection equation. For more complex examples, we could have

ODESystem([(connection1,g1),(connection2,g2)],c)

or something like that so that way you can have different graphs specifying different edges and different equations associated with those edge graphs. The constraint would be that the nodes would have to be the same between all of the graphs, since you're just changing the associated edges.

I think this could cover neuron simulations, building RDME reaction networks, graph-based epidemiological models, and power systems models. Thoughts @isaacsas @asinghvi17 ?

[mforets]

you may be interested in https://github.com/blegat/HybridSystems.jl which implements hybrid automata.

[ChrisRackauckas]

I'm not sure that's related at all. I plan to add RegimeSwitchingSystems, but that's a fairly different problem from network dynamics.

[mforets]

oops, sorry.. i thought this was about an api for hybrid systems.

[ChrisRackauckas]

Maybe I should clarify a bit. It's for network dynamics where you have for example Lotka-Volterra on a network, so you have the Lotka-Volterra equations for x_i and y_i, and then you have terms in the differential equations that connect ODE sets i and j.

[asinghvi17]

Let me rerun some simulations - as I recall, I kind of dropped this effort because (a) the DAE solve was quite a bit slower than the handwritten ODE solve, and there were some differences between the DAE and ODE solutions. If the DAE solution is accurate, then this approach makes sense to me; otherwise, I was going to try and build up the graph "by hand", by modifying the equations before creating a system.

[ChrisRackauckas]

The DAE should be as accurate: if not there's an error in the model form. But indeed, DAEs are as fast as stiff ODE solves (in mass matrix form since it uses the same methods), but if it's a non-stiff ODE then you need #340 in order to convert back down to an ODE for explicit methods.

[asinghvi17]

Also, I think we should change the required arguments of connection to the following:

connection(
    sys_from, 
    sys_to, 
    edge_weight, # can be any arbitrary metadata, even a function
)::Vector{Equation}

This should take care of weighted connections, which are particularly necessary in neuron simulation. It occurs to me that we could have two constructors - a high level one which constructs the graph with systems for you, and a lower level one where you provide the system for each vertex. (We can't carry this information around in the graph, so we may as well carry it separately in a vector indexed by vertex number). This would allow anisotropic graphs, to selectively introduce defects into systems.

[asinghvi17]

ODESystem([(connection1,g1),(connection2,g2)],c)

How would we unify the graphs? As far as I'm aware, there's no mechanism for graph fusion in LightGraphs. I imagine that a MetaDiGraph with connection functions or some user dispatch mechanisms on edges would be the best way forward for this.

[ChrisRackauckas]

I don't think we'd unify graphs. But we'd somehow need a canonical numbering system for the nodes.

[asinghvi17]

As I understand it, LightGraphs handles the node numbering. For example, the constructor for a directed, weighted graph is:

SimpleWeightedDiGraph(12)

whose vertices are then 1:12. As I understand it, this is the canonical system for all LightGraphs graphs.

[ChrisRackauckas]

And if we have a different graph but with the same number of nodes and different edges, how can we node which nodes are the same?

[asinghvi17]

Wouldn't that just be a degenerate graph, then? I'm not sure if we need to know which nodes are the same, since presumably if you pass two graphs with two coupling mechanisms, you want all of the coupling specified...unless I'm misinterpreting the purpose of this.

[ChrisRackauckas]

Well I want to say "this is node 10", and in graph 1, node 10 has an edge to node 11 so add these equations, in graph 2 it does not so don't add these equations.

[asinghvi17]

That's just the vertex number, no? Here's some example code which might clarify this.

using LightGraphs, SimpleWeightedGraphs
graph = SimpleWeightedDiGraph(2)
e1 = SimpleWeightedEdge(1, 2, 1)
e2 = SimpleWeightedEdge(2, 1, 4)
add_edge!(graph, e1)
add_edge!(graph, e2)

looks like this:

[ChrisRackauckas]

Oh, didn't know there was a global index. I thought all indexing was node local because of the edgelists being vector of vector. Okay, then I think we're good.