Symbolic Indexing

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

Provide Symbol Names

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

using NetworkDynamics, Graphs, OrdinaryDiffEqTsit5, Plots
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 create 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]

When constructing component models using ModelingToolkit, the variable names are extracted automatically.

The main types for symbolic indexing are VIndex and EIndex for vertices and edges respectively. Each symbolic index consists of 2 elements: a reference to the network component and a reference to the symbol within that component. Indices may reference states and parameters, but also things like outputs, inputs and observables which do not directly appear in the state/parameter vector.

For accessing by symbol:

VIndex(2, :x) refers to variable with symbolic name :x in vertex number 2.
EIndex(4, :K) refers to parameter :K of the edge component for the 4th edge.

For numeric indexing, use StateIdx and ParamIdx wrappers:

VIndex(2, StateIdx(1)) refers to the first state of vertex 2.
EIndex(4, ParamIdx(2)) refers to the second parameter of edge 4.

Setup code to make following examples work

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 extract points or timeseries from a solution object

ts = 0:0.1:1
sol(ts; idxs=VIndex(1,:storage)) # value of VIndex(1,:storage) at each t in ts

t: 0.0:0.1:1.0
u: 11-element Vector{Float64}:
 -0.21457861470827608
 -0.2285843750805454
 -0.23707930710349903
 -0.24223171646581404
 -0.2453568122720667
 -0.24725230445006693
 -0.24840205216986005
 -0.24909934378380527
 -0.24952214324743344
 -0.24977886132132718
 -0.24993421546857925

Alternatively, you can use them directly in specialized plotting recipes:

plot(sol; idxs=[VIndex(1, :storage), VIndex(5,:storage)]) # plot storage of vertex 1 and vertex 5

It is often advised to choose your timesteps for plotting directly, i.e.

ts = range(0, 1; length=1000)
plot(ts, sol(ts; idxs=VIndex(1,:storage)).u)

gives you more control over how many points are used within the range $(0,1)$.

Generate Symbolic Indices

Often, you need many individual symbolic indices. NetworkDynamics provides several approaches:

Quick Access with Helper Functions

The helper methods vidxs, eidxs, vpidxs and epidxs provide shortcuts for common patterns

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

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

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

Advanced Generation with `generate_indices`

For more complex filtering, use generate_indices, which provides the underlying functionality:

# All edge parameters containing "K" in their name
generate_indices(nw, EIndex(:), "K"; s=false, p=true, out=false)

8-element Vector{NetworkDynamics.SymbolicIndex}:
 EIndex(1, :K)
 EIndex(2, :K)
 EIndex(3, :K)
 EIndex(4, :K)
 EIndex(5, :K)
 EIndex(6, :K)
 EIndex(7, :K)
 EIndex(8, :K)

The helper functions are actually just shortcuts to generate_indices calls:

vidxs(nw, cf, vf) = generate_indices(nw, VIndex(cf), vf; s=true, p=false, out=true, obs=true)

`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)
  ├─ EIndex(1, :K) => 1.0
  ├─ EIndex(2, :K) => 1.0
  ├─ EIndex(3, :K) => 1.0
  ├─ EIndex(4, :K) => 1.0
  ├─ EIndex(5, :K) => 1.0
  ├─ EIndex(6, :K) => 1.0
  ├─ EIndex(7, :K) => 1.0
  └─ EIndex(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 defined in the component. The parameters in the NWParameter object can be accessed using symbolic indices.

p[EIndex(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.

We can set those values like this:

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.0674305185931667
  ├─ VIndex(2, :storage) => -0.37676730636833494
  ├─ VIndex(3, :storage) => -0.6812306797591
  ├─ VIndex(4, :storage) => -0.5720812042830913
  └─ VIndex(5, :storage) => -0.3949489094745284
 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 using the FilteringProxy interface. The filtering proxy can be accessed by calling .v or .e on a state:

s.v

FilteringProxy for NWState()
  Component filter: AllVertices() <- filter further by obj[1], obj["name"], ...
  State filter:     none
  Types:  states ✓  parameters ✓  inputs ✓  outputs ✓  observables ✓
Matching Indices:
  ● VIndex(1, :storage)  -1.0674305   :VertexM
  ● VIndex(2, :storage)  -0.37676731  :VertexM
  ● VIndex(3, :storage)  -0.68123068  :VertexM
  ● VIndex(4, :storage)  -0.5720812   :VertexM
  ● VIndex(5, :storage)  -0.39494891  :VertexM

You can then subsequently filter the list by "indexing" into the FilteringProxy object:

s.v[1:2]

FilteringProxy for NWState()
  Component filter: VIndex([1, 2])
  State filter:     none <- filter states by obj["δ"], obj[:x], ...
  Types:  states ✓  parameters ✓  inputs ✓  outputs ✓  observables ✓
Matching Indices:
  ● VIndex(1, :storage)  -1.0674305   :VertexM
  ● VIndex(2, :storage)  -0.37676731  :VertexM

... until you've filtered on both components and variables and the indexing returns the actual values:

s.v[1:2][:storage]

2-element Vector{Float64}:
 -1.0674305185931667
 -0.37676730636833494

Check out the docstring of FilteringProxy for an in-depth explanation.

The NWState and NWParameter objects are mutable, thus changing them will also change the underlying wrapped flat arrays.

For example, we can use the syntax introduced above to update the values within the state

s.v[1:2][:storage] = [1,2]

We can confirm the update by inspecting the values again:

s.v

FilteringProxy for NWState()
  Component filter: AllVertices() <- filter further by obj[1], obj["name"], ...
  State filter:     none
  Types:  states ✓  parameters ✓  inputs ✓  outputs ✓  observables ✓
Matching Indices:
  ● VIndex(1, :storage)   1           :VertexM
  ● VIndex(2, :storage)   2           :VertexM
  ● VIndex(3, :storage)  -0.68123068  :VertexM
  ● VIndex(4, :storage)  -0.5720812   :VertexM
  ● VIndex(5, :storage)  -0.39494891  :VertexM

You can always access the flat representations by calling uflat and pflat. The ordering of elements in these flat arrays corresponds exactly to the order returned by variable_symbols and parameter_symbols respectively.

Note

The NWState and NWParameter wrappers can be constructed from various objects. For example, within a callback you might construct p = NWParameter(integrator) to then change the parameters of the network within the callback.

Interactive Filtering with FilteringProxy

The FilteringProxy system provides an intuitive way to explore and access network variables through progressive filtering. When you access .v, .e, or .p properties of NWState/NWParameter objects, you get a filtering proxy that can be refined step by step.

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 automatically as observables.

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(nw, 1:2, :storage))

2-element Vector{Int64}:
 1
 2

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

true

Analogous with parameters:

idxs = SII.parameter_index(nw, [EIndex(1, :K), EIndex(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.SymbolicIndex{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.SymbolicIndex{Int64, Symbol}}:
 EIndex(1, :K)
 EIndex(2, :K)
 EIndex(3, :K)
 EIndex(4, :K)
 EIndex(5, :K)
 EIndex(6, :K)
 EIndex(7, :K)
 EIndex(8, :K)

These functions return the symbolic indices in the exact order they appear in the flat arrays returned by uflat and pflat, making them essential when you need to map between flat array indices and symbolic representations.

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