Symbolic Indexing

Using SciML's SymbolicIndexingInterface.jl, ND.jl provides numerous methods to access and change variables and parameters.

using NetworkDynamics
using Graphs
using OrdinaryDiffEqTsit5
using Plots

Provide Symbol Names

When constructing component models, you can pass symbolic names using the sym and psym keywords.

function _edgef!(e, v_s, v_d, (K,), t)
    e .= K * (v_s[1] .- v_d[1])
end
edgef = EdgeModel(;g=AntiSymmetric(_edgef!), outsym=[:flow], psym=[:K=>1])

EdgeModel :StaticEdgeM PureFeedForward()
 ├─   0 states:  []  
 ├─ 1/1 outputs: src=[₋flow] dst=[flow]
 └─   1 param:   [K=1]

Here we created a static diffusion edge with suitable variable and parameter names. Similarly, we define the diffusion vertex with symbolic names.

function _vertexf!(dv, v, esum, p, t)
    dv[1] = esum[1]
end
vertexf = VertexModel(f=_vertexf!, g=1, sym=[:storage])

VertexModel :VertexM PureStateMap()
 ├─ 1 state:  [storage]
 └─ 1 output: [storage]

Fundamental Symbolic Indices

The default types for this access are the types VIndex, EIndex, VPIndex and EPIndex. Each of those symbolic indices consists of 2 elements: a reference to the network component and a reference to the symbol within that component. As such, VIndex(2, :x) refers to variable with symbolic name :x in vertex number 2. EPIndex(4, 2) would refer to the second parameter of the edge component for the 4th edge.

g = wheel_graph(5)
nw = Network(g, vertexf, edgef)
s = NWState(nw)
s.v[:,:storage] .= randn(5)
prob = ODEProblem(nw, uflat(s), (0,2), pflat(s))
sol = solve(prob, Tsit5())

Those fundamental indices can be used in a lot of scenarios. Most importantly you can use them to

sol(sol.t; idxs=VIndex(1, :storage))   # extract timeseries out ouf solution object
plot(sol; idxs=[VIndex(1, :storage), VIndex(5,:storage)]) # plot storage of two nodes

Generate Symbolic Indices

Often, you need many individual symbolic indices. For that there are the helper methods vidxs, eidxs, vpidxs and epidxs. With the help of those methods you can generate arrays of symbolic indices:

vidxs(nw, :, :storage) # get variable "storage" for all vertices

5-element Vector{VIndex}:
 VIndex(1, :storage)
 VIndex(2, :storage)
 VIndex(3, :storage)
 VIndex(4, :storage)
 VIndex(5, :storage)

plot(sol; idxs=vidxs(nw, :, :storage))

NWState and NWParameter Objects

Internally, both state and parameters of a Network are represented using flat arrays. To access the state or parameters of a network, you can use the NWState and NWParameter objects.

p = NWParameter(nw)

Parameter{Vector{Float64}} of Network (5 vertices, 8 edges)
  ├─ EPIndex(1, :K) => 1.0
  ├─ EPIndex(2, :K) => 1.0
  ├─ EPIndex(3, :K) => 1.0
  ├─ EPIndex(4, :K) => 1.0
  ├─ EPIndex(5, :K) => 1.0
  ├─ EPIndex(6, :K) => 1.0
  ├─ EPIndex(7, :K) => 1.0
  └─ EPIndex(8, :K) => 1.0

creates a NWParameter object for the network nw. It essentially creates a new flat parameter array and fills it with the default parameter values define in the component. The parameters in the NWParameter object can be accessed using symbolic indices.

p[EPIndex(5, :K)] = 2.0 # change the parameter K of the 5th edge

Similarly, you can create a NWState object for the network nw using

s = NWState(nw)

NWState{Vector{Float64}} of Network (5 vertices, 8 edges)
  ├─ VIndex(1, :storage) => NaN
  ├─ VIndex(2, :storage) => NaN
  ├─ VIndex(3, :storage) => NaN
  ├─ VIndex(4, :storage) => NaN
  └─ VIndex(5, :storage) => NaN
 p = NWParameter([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
 t = nothing

No default values were provided in the network components, so the state array is filled with NaN values.

s[VIndex(:, :storage)] .= randn(5) # set the (initial) storage for all vertices

NWState{Vector{Float64}} of Network (5 vertices, 8 edges)
  ├─ VIndex(1, :storage) => 1.172913653567515
  ├─ VIndex(2, :storage) => -1.6288456218595309
  ├─ VIndex(3, :storage) => 0.3379821353790073
  ├─ VIndex(4, :storage) => 1.2415725230351353
  └─ VIndex(5, :storage) => 0.7647289771919843
 p = NWParameter([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
 t = nothing

For both NWState and NWParameter objects, there is a more convenient way to access the variables and parameters.

@assert s.v[1, :storage] == s[VIndex(1, :storage)] # s.v -> access vertex states
@assert s.e[1, :flow]    == s[EIndex(1, :flow)]    # s.e -> access edge states
@assert s.p.e[1,:K]      == p[EPIndex(1, :K)]      # s.p -> access parameters

The NWState and NWParameter objects are mutable, thus changing them will also change the underlying wrapped flat arrays. You can always access the flat representations by calling uflat and pflat.

Observables

Sometimes, the "states" you're interested in aren't really states in the DAE sense but rather algebraic derivations from DAE states, parameters, and time – in accordance with the naming in the SciML ecosystem, these values are called Observables.

A prime example of Observables are edge/vertex-outputs, such as the flow in the edge model defined above. It is also possible to define additional Observables manually by using the obssym and obsf keyword on the EdgeModel/VertexModel constructors. When building models using ModelingToolkit, the reduced algebraic states will be preserved as observables automatically.

Observables can be accessed like any other state, for example, the flows in the network don't show up in the state array but can be accessed in all the ways discussed above, for example

plot(sol; idxs=eidxs(nw, :, :flow))

Derived ObservableExpressions using @obsex

Sometimes it is useful to plot or observe simple derived quantities. For that, one can use the @obsex macro to define simple derived quantities.

For example, we can directly plot the storage difference with respect to storage of node 1.

plot(sol; idxs=@obsex(vidxs(nw,:,:storage) .- VIndex(1,:storage)))

Other examples include calculating the magnitude and argument of complex values that are modeled using real and imaginary parts.

@obsex mag = sqrt(VIndex(1, :u_r)^2 + VIndex(2, :u_i)^2)

Low-level accessors for flat array indices

Sometimes, you want to know the indices of your states in the flat arrays. For that, you can use the low-level methods defined in SymbolicIndexingInterface.jl:

using NetworkDynamics: SII # SII = SymbolicIndexingInterface
idxs = SII.variable_index(nw, vidxs(1:2, :storage))

2-element Vector{Int64}:
 1
 2

uflat(s)[idxs] == s.v[1:2, :storage]

true

Analogous with parmeters:

idxs = SII.parameter_index(nw, eidxs(1:2, :K))
pflat(s)[idxs] == s.p.e[1:2, :K]

true

If you need the symbols of all the states/parameters in order, you can use

SII.variable_symbols(nw)

5-element Vector{NetworkDynamics.SymbolicStateIndex{Int64, Symbol}}:
 VIndex(1, :storage)
 VIndex(2, :storage)
 VIndex(3, :storage)
 VIndex(4, :storage)
 VIndex(5, :storage)

and

SII.parameter_symbols(nw)

8-element Vector{NetworkDynamics.SymbolicParameterIndex{Int64, Symbol}}:
 EPIndex(1, :K)
 EPIndex(2, :K)
 EPIndex(3, :K)
 EPIndex(4, :K)
 EPIndex(5, :K)
 EPIndex(6, :K)
 EPIndex(7, :K)
 EPIndex(8, :K)

All above examples also work on other "symbolic containers", e.g. SII.variable_symbols(::NWState).