diff --git a/docs/examples/getting_started_with_network_dynamics.jl b/docs/examples/getting_started_with_network_dynamics.jl index f1957fb20..20ad86c86 100644 --- a/docs/examples/getting_started_with_network_dynamics.jl +++ b/docs/examples/getting_started_with_network_dynamics.jl @@ -1,177 +1,330 @@ #= # [Network Diffusion](@id Getting Started) -This introductory example explains the use of the basic types and constructors -in NetworkDynamics.jl by modeling a simple diffusion on an undirected network. +This example explains the use of the basic functions and constructors in NetworkDynamics.jl by modeling a simple +diffusion on an undirected network. -This example can be dowloaded as a normal Julia script [here](@__NAME__.jl). #md +This page will: + * walk you through the theoretical background of a simple diffusion propagating in an undirected network + * explain the network dynamics of the system + * explain how to program them -## Theoretical background +!!! note + An Undirected Network is a network where all the connections between the nodes are bidirectional. + +This example can be downloaded as a normal Julia script [here](@__NAME__.jl). #md +=# -Diffusion processes are relevant for phenomena as diverse as heat conduction, electrical currents, and random walks. Generally speaking they describe the tendency of systems to evolve towards a state of equally distributed heat, charge or concentration. In such system the local temperature (or concentration) changes according to its difference with its neighborhood, i.e. the temperature gradient. +#= +## Theoretical background -Let $g$ be a graph with $N$ nodes and adjacency matrix $A$. Let $v = (v_1, \dots, v_n)$ be a vector of (abstract) temperatures or concentrations at each node $i = 1, \dots, N$. Then the rate of change of state $v_i$ is described by its difference with its neighbors and we obtain the following ordinary differential equation +Diffusion processes appear in phenomena as diverse as heat conduction, electrical currents, and random walks. +Generally speaking they describe the tendency of systems to evolve towards a state of equally distributed entropy +(the entropy being in the form of e.g. heat, charge or concentration). If we assume a thermal system, the temperature +of a specific spot changes depending on the temperature gradient between the temperature of the spot itself and that of +its neighborhood. +We will build a graph $g$ with $N$ nodes and an adjacency matrix $A$, where $v = (v_1, \dots, v_n)$ is the vector of the +(abstract) temperatures at each node $i = 1, \dots, N$. The rate of change of state $v_i$ will be described by the +difference between the temperature of the node and that of its neighbours. From the above we obtain the following ordinary +differential equation: ```math \dot v_i = \sum_{j=1}^N A_{ji} (v_j - v_i). ``` -The sum on the right hand side plays the role of a (discrete) gradient. If the temperature at node $i$ is higher than at its neighboring node $j$ it will decrease along that edge. +The sum on the right hand side plays the role of a (discrete) gradient. If the temperature at node $i$ is higher than +at its neighboring node $j$ it will decrease along that edge. -## Modeling diffusion in NetworkDynamics.jl -We begin by loading the necessary packages. -=# -using Graphs -using NetworkDynamics -using OrdinaryDiffEqTsit5 -using StableRNGs -using Plots -nothing #hide -#= +If a node were to be disconnected from the rest of the network its state would never change, because its adjacency +matrix would be $A_{ji} = 0 \; \forall j$ and hence $\dot v_i = 0$. So because its state would remain unchanged the +model will be comprised of nodes that have no internal dynamics. This means that the evolution of a node will depend +only on its interaction with its neighbors. In NetworkDynamics.jl, interactions between a node and its neighbors are +described using edge equations. + +In order to bring this equation into the form required by NetworkDynamics.jl we need to split the dynamics into edge +dynamics and vertex dynamics and bring them into the correct input-output formulation. -From the above considerations we see that in this model the nodes do not have any internal dynamics - if a node was disconnected from the rest of the network its state would never change, since then $A_{ji} = 0 \; \forall j$ and hence $\dot v_i = 0$. This means that the evolution of a node depends only on the interaction with its neighbors. In NetworkDynamics.jl, interactions with neighbors are described by equations for the edges. -In order to bring this equation into the form required by NetworkDynamics.jl we need split the dynamics into edge and vertex parts and bring them into the correct input-output formulation. -The vertices have one internal state $v$ which is also the output. The input is -the sum over all flows of connected edges. This directly correspons to the component model definition outlined in [Mathematical Model](@ref): +### Vertex Dynamics: +All vertices have one internal state $v$, which is also the vertex output. It is the sum over all incoming edge +flows connected to the vertex. This directly corresponds to the component model definition outlined in [Mathematical Model](@ref): ```math \begin{aligned} \dot x^\mathrm{v} &= f^{\mathrm v}(u^{\mathrm v}, \sum_k^{\text{incident}} y^{\mathrm e}_k, p^{\mathrm v}, t) &&= \sum_k^\mathrm{incident} y^\mathrm{e}_k \\ y^\mathrm{v} &= g^{\mathrm v}(u^{\mathrm v}, \sum_k^{\text{incident}} y^{\mathrm e}_k, p^{\mathrm v}, t) &&= x^\mathrm{v} \end{aligned} ``` -The edge dynamics on the other hand do not have any internal states. Thus we -only define the output as the difference between the source and destination -vertex: +=# + +#= +### Edge Dynamics: +The edge dynamics on the other hand do not have any internal states. Thus we can define the edge output as the +difference between the source and destination vertices: ```math \begin{aligned} y^{\mathrm e}_{\mathrm{dst}} &= g_\mathrm{dst}^{\mathrm e}(u^{\mathrm e}, y^{\mathrm v}_{\mathrm{src}}, y^{\mathrm v}_{\mathrm{dst}}, p^{\mathrm e}, t) &&= y^{\mathrm v}_{\mathrm{src}} - y^{\mathrm v}_{\mathrm{dst}}\\ y^{\mathrm e}_{\mathrm{src}} &= g_\mathrm{src}^{\mathrm e}(u^{\mathrm e}, y^{\mathrm v}_{\mathrm{src}}, y^{\mathrm v}_{\mathrm{dst}}, p^{\mathrm e}, t) &&= y^{\mathrm v}_{\mathrm{dst}} - y^{\mathrm v}_{\mathrm{src}} \end{aligned} ``` +=# + +#= +## Modelling dynamics in NetworkDynamics.jl +To model the vertex and edge dynamics we need to create a `VertexModel` and an `EdgeModel`, respectively. + +As a first step we need to install Julia and the necessary packages (`Graphs`, NetworkDynamics`, `OrdinaryDiffEqTsit5`, +`StableRNGs` and `Plots`). Then we need to load them: +=# + +using Graphs +using NetworkDynamics +using OrdinaryDiffEqTsit5 +using StableRNGs +using Plots +nothing #hide -### Definition of `EdgeModel` + +#= + +### Defining the `EdgeModel` + +Then we must define the `EdgeModel`. To define it we use the function `diffusionedge_g!` which takes as inputs +the current state of the edge `e`, its source vertex `v_src`, its destination vertex `v_dst`, a vector of parameters `p` +and the time `t`. In order to comply with the syntax of NetworkDynamics.jl we must always define functions for edges +with exactly these arguments. (In the case of this example, the values for `p` and `t` are not used since there are no +parameters and there is no explicit time dependency in the system). + +After the function call the edge's output value `e` equals the difference between its source and its destination vertex + (i.e. the discrete gradient along that edge). + +!!! note + `diffusionedge_g!` is called a **mutating** function, because it mutates (modifies) the edge state `e` (which is the + first of its inputs). (In Julia, names of mutating functions end with an `!`.) We use mutating functions because + they reduce the computational resource allocations and therefore, speed up the computations. =# + function diffusionedge_g!(e_dst, v_src, v_dst, p, t) - ## e_dst, v_src, v_dst are arrays, hence we use the broadcasting operator + ## e_dst, v_src, v_dst are arrays, so we use the broadcasting operator . e_dst .= v_src .- v_dst nothing end nothing #hide #md #= -The function `diffusionedge_g!` takes as inputs the current state of the edge `e`, its source vertex `v_src`, its destination vertex `v_dst`, a vector of parameters `p` and the time `t`. In order to comply with the syntax of NetworkDynamics.jl we always have to define functions for edges with exactly these arguments, even though we do not need `p` and `t` for the diffusion example. - -`diffusionedge_g!` is called a **mutating** function, since it modifies (or *mutates*) one of its inputs, namely the edge state `e`. As a convention in Julia names of mutating functions end with an `!`. The use of mutating functions reduces allocations and thereby speeds up computations. After the function call the edge's output value `e` equals the difference between its source and its destination vertex (i.e. the discrete gradient along that edge). - -Notably, this function only models $g_\mathrm{dst}$. However we can wrap this single-sided output function in an [`AntiSymmetric`](@ref) output wrapper to construct the [`EdgeModel`](@ref): +Notably, this function only models $g_\mathrm{dst}$. However we can wrap this single-sided output function in an +[`AntiSymmetric`](@ref) output wrapper to construct the [`EdgeModel`](@ref): =# + nd_diffusion_edge = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow]) + #= -### Definition of `VertexModel` -For undirected graphs, the `edgefunction!` specifies the coupling from a source- to a destination vertex. The contributions of the connected edges to a single vertex are "aggregated". Default aggregation is the summation of all incident edge states. The aggregated edge state is made available via the `esum` argument of the vertex function. +### Defining the `VertexModel` +Next we need to define the `VertexModel`. For undirected graphs, the `edgefunction!` specifies the coupling from a +source to a destination vertex. The contributions of the connected edges to a single vertex are "aggregated". +By default, the aggregation is the summation of all incident edge states. The aggregated edge state is made available +via the `esum` argument of the vertex function. =# + function diffusionvertex_f!(dv, v, esum, p, t) - ## dv, v and esum are arrays, hence we use the broadcasting operator . + ## dv, v and esum are arrays, so we use the broadcasting operator . dv .= esum nothing end nothing #hide #md #= -Just like above the input arguments `v, esum, p, t` are mandatory for the syntax of vertex functions. The additional input `dv` corresponding to the derivative of the vertex' state is mandatory for vertices described by ordinary differential equations. - -The output function `g` is just taking part of the internal states. For that we can use the [`StateMask`](@ref) helper function `g = StateMask(1:1)` +Just like above, the input arguments `v, esum, p, t` are mandatory for the syntax of vertex functions. The additional + input `dv` corresponding to the derivative of the vertex' state is mandatory for vertices described by ordinary + differential equations. The outputs of the vertex function are just a subset of the internal states. Therefore, we + can use the [`StateMask`](@ref) helper function `g = StateMaks(1:1)` to access them. =# -nd_diffusion_vertex = VertexModel(; f=diffusionvertex_f!, g=StateMask(1:1), dim=1) +nd_diffusion_vertex = VertexModel(; f=diffusionvertex_f!, g=StateMask(1:1), dim=1) #= + ## Constructing the network -With the components defined, we can define the topology and assemble the network dynamics. +With the components defined, we can now define the topology and assemble the network dynamics. =# +# The first step is to define the number of nodes the network is comprised of: N = 20 # number of nodes -k = 4 # average degree -g = barabasi_albert(N, k) # a little more exciting than a bare random graph - +k = 4 # average degree distribution nothing #hide #md -# The [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model) generates a scale-free random graph. +# Then we define the network model. We will use the +# [Barabási–Albert model](https://en.wikipedia.org/wiki/Barab%C3%A1si%E2%80%93Albert_model) which generates a +# scale-free random graph. +g = barabasi_albert(N, k) +nothing #hide #md +#= +We then create the network of nodes and edges by using the constructor `Network`. It combines the component +model with the topological information contained in the graph **`g`** and returns an `Network` compatible with the +solvers of `DifferentialEquations.jl`. +=# nd = Network(g, nd_diffusion_vertex, nd_diffusion_edge) +rng = StableRNG(1) + +# We then generate random initial conditions +x0 = randn(rng, N) #= -The constructor `Network` combines the component model with the topological information contained in the graph **`g`** and returns an `Network` compatible with the solvers of `DifferentialEquations.jl`. +We then solve the diffusion problem on the time interval $[0, 2]$ with the `Tsit5()` algorithm, which is recommended +by the authors of `DifferentialEquations.jl` for most non-stiff problems. We first create an `ODEProblem`: =# - -rng = StableRNG(1) -x0 = randn(rng, N) # random initial conditions ode_prob = ODEProblem(nd, x0, (0.0, 2.0)) -Main.test_execution_styles(ode_prob) # testing all ex styles #src + +# We test all ex styles #src +Main.test_execution_styles(ode_prob) #src + +# And then solve the network using ´solve´: sol = solve(ode_prob, Tsit5()); nothing #hide #md #= -We are solving the diffusion problem on the time interval $[0, 2]$ with the `Tsit5()` algorithm, which is recommended by the authors of `DifferentialEquations.jl` for most non-stiff problems. +We plot the results using the `plot` command, which takes three parameters as input `sol`, `idxs` and `fmt`. +The first parameter is the solved network `sol`. +The second parameter is `idxs`, which provides a set of indices. We can use "symbolic" indices, which specify specific +components and their symbols directly (see [Symbolic Indexing](@ref) for more details). In order to collect multiple +indices we can use the helper function [`vidxs`](@ref) and [`eidxs`](@ref), which help us collect all symbolic indices +matching specific criteria. The third and last parameter is `fmt`, which provides the format of the generated plot, in +this case that format being .png. =# - plot(sol; idxs=vidxs(nd, :, :), fmt=:png) #= -The plotting is straightforward. The **`idxs`** keyword allows us to pass a list of indices. Indices can be also "symbolic" indices which specify components and their symbols directly. For example `idxs = VIndex(1, :v)` acesses state `:v` of vertex 1. See [Symbolic Indexing](@ref) for more details. - -In oder to collect multiple indices we can use the helper function [`vidxs`](@ref) and [`eidxs`](@ref), which help to collect all symbolic indices matching a certain criteria. +(@Hans can you provide a description of what we see in this plot? The y axis is missing so it is unclear to me.) ## Two Dimensional Extension -To illustrate a very simple multi-dimensional case, in the following we simulate two independent diffusions on an identical graph. The first uses the symbol `x` and is started with initial conditions drawn from the standard normal distribution $N(0,1)$, the second uses the symbol `ϕ` with squared standard normal inital conditions. +To illustrate a very simple multi-dimensional case, we simulate two independent diffusions on an identical graph. +The first diffusion uses the symbol `x` and is initiated with initial conditions drawn from the standard normal +distribution $N(0,1)$, while the second diffusion uses the symbol `ϕ` with squared standard normal initial conditions. The symbols have to be passed with the keyword **`sym`** to `VertexModel`. + +Once again we define the number of nodes and the type of network we want to generate: =# N = 10 # number of nodes -k = 4 # average degree -g = barabasi_albert(N, k) # a little more exciting than a bare random graph +k = 4 # average degrees distribution +g = barabasi_albert(N, k) -## We will have two independent diffusions on the network, hence dim = 2 +# We have two independent diffusions on the network, hence dim = 2. We now define the `VertexModel` and the `EdgeModel` nd_diffusion_vertex_2 = VertexModel(; f=diffusionvertex_f!, g=1:2, dim=2, sym=[:x, :ϕ]) nd_diffusion_edge_2 = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow_x, :flow_ϕ]) nd_2 = Network(g, nd_diffusion_vertex_2, nd_diffusion_edge_2) -x0_2 = vec(transpose([randn(rng, N) .^ 2 randn(rng, N)])) # x ~ N(0,1)^2; ϕ ~ N(0,1) +#= We define the first diffusion as x ~ $N(0,1)^\mathrm{2}$ and the second diffusion as ϕ ~ $N(0,1)$. So the propagation +of the diffusion from node 0 to node 2 is given by creating a vector: +=# +x0_2 = vec(transpose([randn(rng, N) .^ 2 randn(rng, N)])) +#(@Hans: is my explanation correct?) + +# We then define the [ODEProblem](@ref) ode_prob_2 = ODEProblem(nd_2, x0_2, (0.0, 3.0)) -Main.test_execution_styles(ode_prob_2) # testing all ex styles #src +Main.test_execution_styles(ode_prob_2) #src + +# Then solve the network using: sol_2 = solve(ode_prob_2, Tsit5()); -# Try plotting the variables ϕ_i yourself. [To write ϕ type \phi and press TAB] +# To plot the evolution of variables ϕ_i over time we use the command: plot(sol_2; idxs=vidxs(nd_2, :, :x), fmt=:png) +# [To write ϕ in the terminal type \phi and press TAB] -# Using the `eidxs` helper function we can also plot the flow variables +# Using the `eidxs` helper function we can also plot the evolution of the flow variables over time: plot(sol_2; idxs=eidxs(nd_2, :, :flow_x), fmt=:png) + #= -## Appendix: The network Laplacian $L$ +## Putting it all together -The diffusion equation on a network can be rewritten as +(@Hans: when I am trying to run the script in the "Putting it all together" in a separate .jl file I am getting a +UndefVarError: `AntiSymmetric` not defined in `Main` +Suggestion: check for spelling errors or missing imports. +Stacktrace: + [1] top-level scope + @ REPL[7]:1) +and I have no idea why. Can you please check it? Because if someone copy-pastes it, it will not work for them) +=# + +using Graphs +using NetworkDynamics +using OrdinaryDiffEqTsit5 +using StableRNGs +using Plots +nothing #hide + +function diffusionedge_g!(e_dst, v_src, v_dst, p, t) + ## e_dst, v_src, v_dst are arrays, hence we use the broadcasting operator + e_dst .= v_src .- v_dst + nothing +end +nothing #hide #md + +nd_diffusion_edge = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow]) + +function diffusionvertex_f!(dv, v, esum, p, t) + ## dv, v and esum are arrays, hence we use the broadcasting operator . + dv .= esum + nothing +end +nothing #hide #md + +nd_diffusion_vertex = VertexModel(; f=diffusionvertex_f!, g=StateMask(1:1), dim=1) + +N = 20 +k = 4 +g = barabasi_albert(N, k) +nothing #hide #md + +nd = Network(g, nd_diffusion_vertex, nd_diffusion_edge) +rng = StableRNG(1) +x0 = randn(rng, N) +ode_prob = ODEProblem(nd, x0, (0.0, 2.0)) +Main.test_execution_styles(ode_prob) #src +sol = solve(ode_prob, Tsit5()); +nothing #hide #md +plot(sol; idxs=vidxs(nd, :, :), fmt=:png) +N = 10 # +k = 4 # +g = barabasi_albert(N, k) +nd_diffusion_vertex_2 = VertexModel(; f=diffusionvertex_f!, g=1:2, dim=2, sym=[:x, :ϕ]) +nd_diffusion_edge_2 = EdgeModel(; g=AntiSymmetric(diffusionedge_g!), outsym=[:flow_x, :flow_ϕ]) +nd_2 = Network(g, nd_diffusion_vertex_2, nd_diffusion_edge_2) +x0_2 = vec(transpose([randn(rng, N) .^ 2 randn(rng, N)])) +ode_prob_2 = ODEProblem(nd_2, x0_2, (0.0, 3.0)) +Main.test_execution_styles(ode_prob_2) #src +sol_2 = solve(ode_prob_2, Tsit5()); +plot(sol_2; idxs=vidxs(nd_2, :, :x), fmt=:png, xlabel = "x", ylabel = "nd_2",) +plot(sol_2; idxs=eidxs(nd_2, :, :flow_x), fmt=:png) + + + +#= +## Appendix: The network Laplacian $L$ + +The diffusion equation on a network can be rewritten as: ```math -\dot v_i = \sum_{j=1}^N A_{ji} v_j - d_i v_i = e_i^T A v - d_i v_i +\dot v_i = \sum_{j=1}^N A_{ji} v_j - d_i v_i = e_i^T A v - d_i v_i ``` +where $d_i$ is the degree distribution of node $i$ and $e_i^T$ is the $i$-th standard basis vector. -where $d_i$ is the degree of node $i$ and $e_i^T$ is the $i$-th standard basis vector. Introducing the diagonal matrix $D$ that has the degree of node $i$ in its $i$-th row and the Laplacian matrix $L = D - A$ we arrive at - +By introducing the diagonal matrix $D$, that has the degree of node $i$ in its $i$-th row and the Laplacian matrix +$L = D - A$ we arrive at: ```math \dot v = e_i^T(A - D) v ``` - and finally - ```math \dot v = - L v ``` -This is a linear system of ODEs and its solution is a matrix exponential. To study the asymptotic behaviour of the system it suffices to analyze the eigenspectrum of $L$. For this reason $L$ is an important construction in network science. +This is a linear system of ordinary differential equations (ODEs) and its solution is a matrix exponential. +To study the asymptotic behaviour of the system it suffices to analyze the eigenspectrum of $L$. + For this reason $L$ is an important construction in network science. =# diff --git a/docs/examples/heterogeneous_system.jl b/docs/examples/heterogeneous_system.jl index 0c9cdc7e7..2f7112ecd 100644 --- a/docs/examples/heterogeneous_system.jl +++ b/docs/examples/heterogeneous_system.jl @@ -1,13 +1,13 @@ #= # Modeling a heterogeneous system -This example can be dowloaded as a normal Julia script [here](@__NAME__.jl). #md +This example can be downloaded as a normal Julia script [here](@__NAME__.jl). #md -One of the main purposes of NetworkDynamics.jl is to facilitate -modeling coupled systems with heterogenities. This means that -components can differ in their parameters as well as in their dynamics. +One of the main purposes of NetworkDynamics.jl is to facilitate the +modeling systems whose components can differ in their parameters as well as their dynamics. +These are known as Coupled systems with heterogeneities. -## Heterogenous parameters +## Heterogeneous parameters We start by setting up a simple system of Kuramoto oscillators. =# diff --git a/docs/src/assets/edgemodel.svg b/docs/src/assets/edgemodel.svg new file mode 100644 index 000000000..e876a8ed8 --- /dev/null +++ b/docs/src/assets/edgemodel.svg @@ -0,0 +1,1246 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/edgemodel.tex b/docs/src/assets/edgemodel.tex new file mode 100644 index 000000000..b4c8c19fb --- /dev/null +++ b/docs/src/assets/edgemodel.tex @@ -0,0 +1,24 @@ +% \documentclass[tikz,convert={outfile=\jobname.svg}]{standalone} +\documentclass[tikz,convert=pdf2svg]{standalone} +%\usetikzlibrary{...}% tikz package already loaded by 'tikz' option +\usetikzlibrary{arrows, arrows.meta} +% \tikzset{ +% >=stealth', +% } +\usepackage{amsmath,amssymb} +\usepackage{tikz} +\usepackage{xcolor} + +\begin{document} +\begin{tikzpicture} + \node[draw, minimum height=1.5cm](n){ + Edge Model + }; + \draw[->](n.east)++(0,0.4)--++(1.0,0) node[right]{$y^{\mathrm e}_{\mathrm{src}}$}; \draw[->](n.east)++(0,-0.4)--++(1.0,0) node[right]{$y^{\mathrm e}_{\mathrm{dst}}$}; + + \draw[<-](n.west)++(0,0.4)--++(-1.0,0) node[left]{$y^{\mathrm v}_{i}=i^{\mathrm e}_{\mathrm{src}}$}; + \draw[<-](n.west)++(0,-0.4)--++(-1.0,0) node[left]{$y^{\mathrm v}_{j}=i^{\mathrm e}_{\mathrm{dst}}$}; + + \node[text width=3.5cm, align=center, font=\scriptsize] at ([xshift=-2.7cm, yshift=-1.2cm]n) {mapping of adjecent node outputs as edge inputs}; +\end{tikzpicture} +\end{document} diff --git a/docs/src/assets/makelatexfigs b/docs/src/assets/makelatexfigs index de9feb1d6..4c0d7a00a 100755 --- a/docs/src/assets/makelatexfigs +++ b/docs/src/assets/makelatexfigs @@ -18,4 +18,6 @@ if ! command_exists pdf2svg; then fi latexmk -pdflua -shell-escape mathmodel.tex +latexmk -pdflua -shell-escape edgemodel.tex +latexmk -pdflua -shell-escape nodemodel.tex latexmk -C diff --git a/docs/src/assets/mathmodel.tex b/docs/src/assets/mathmodel.tex index 525421c75..b27c6ff40 100644 --- a/docs/src/assets/mathmodel.tex +++ b/docs/src/assets/mathmodel.tex @@ -8,7 +8,6 @@ \usepackage{amsmath,amssymb} \usepackage{tikz} \usepackage{xcolor} - \makeatletter \def\mathcolor#1#{\@mathcolor{#1}} \def\@mathcolor#1#2#3{% @@ -18,14 +17,13 @@ \endgroup } \makeatother - \definecolor{wblue}{HTML}{0072b2} \definecolor{worange}{HTML}{e69f00} \definecolor{wbblue}{HTML}{56b4e9} \definecolor{wgreen}{HTML}{009e73} \definecolor{wyellow}{HTML}{f0e442} - \begin{document} +\pagecolor{white} \begin{tikzpicture} \node[draw, circle](n1){1}; \node[draw, circle] at (10cm, 1cm)(n2){2}; @@ -39,7 +37,6 @@ \draw[] (n3) --++(1,-1); \draw[] (n1) --++(-1,-2.5); \draw[] (n2) --++(-2,-3.5); - \node[text width=5cm, xshift=-0.5cm, yshift=1.25cm] at (n1){% \small\begin{align*} M\,\dot{x} &= f^{\mathrm{v}}_{1}\left(x, \sum_{k} \mathcolor{worange}{e_{k1}}, p, t\right)\\ @@ -52,7 +49,6 @@ \mathcolor{wblue}{v_{2}} &= g^{\mathrm{v}}_{2}\left(x, p_{\mathrm{v},2}, t\right)\\ \end{align*}% }; - \path (n1)--(n2) node[pos=.5,text width=5cm, xshift=0cm, yshift=1.25cm]{% \small\begin{align*} M\,\dot{x} &= f^{\mathrm{e}}_{12}\left(x, \mathcolor{worange}{v_{\mathrm{src}}}, \mathcolor{wblue}{v_{\mathrm{dst}}}, p, t\right)\\ @@ -60,6 +56,5 @@ \mathcolor{wblue}{e_{12}} &= g^{\mathrm{e}}_{\mathrm{dst}}\left(x, \mathcolor{worange}{v_{\mathrm{src}}}, \mathcolor{wblue}{v_{\mathrm{dst}}}, p, t\right)\\ \end{align*}% }; - \end{tikzpicture} \end{document} diff --git a/docs/src/assets/nodemodel.svg b/docs/src/assets/nodemodel.svg new file mode 100644 index 000000000..f56ac2d0c --- /dev/null +++ b/docs/src/assets/nodemodel.svg @@ -0,0 +1,983 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/nodemodel.tex b/docs/src/assets/nodemodel.tex new file mode 100644 index 000000000..0e2ce7c56 --- /dev/null +++ b/docs/src/assets/nodemodel.tex @@ -0,0 +1,35 @@ +% \documentclass[tikz,convert={outfile=\jobname.svg}]{standalone} +\documentclass[tikz,convert=pdf2svg,border=20pt]{standalone} +%\usetikzlibrary{...}% tikz package already loaded by 'tikz' option +\usetikzlibrary{arrows, arrows.meta} +% \tikzset{ +% >=stealth', +% } + +\usepackage{amsmath,amssymb} +\usepackage{tikz} +\usepackage{xcolor} + +\begin{document} +\begin{tikzpicture} + % Write "Node model" + \node[draw, font=\tiny](n){Node Model}; + + % Write y^v + \draw[->](n.east)--++(1.0,0) node[above, pos=.5, font=\tiny]{$y^{\mathrm v}$}; + + % Draw the accumulator circle + \node[font=\tiny, draw, circle] (acc) at ([xshift=-1.5cm]n.west) {$\mathrm{acc}$}; + + % Connect acc to the node + \draw[->] (acc) -- (n.west) node[above, pos=.5, font=\tiny]{$i^{v}$}; + + % Add three arrows to acc with rectangular corners and write y^e_i, y^e_j and ... + \draw[<-] (acc) |- ++(-1.0,0.7) node[left, font=\tiny](in){$y^e_i$}; + \draw[<-] (acc) -- ++(-1.0,0) node[left, font=\tiny]{$y^e_j$}; + \draw[<-] (acc) |- ++(-1.0,-0.7) node[left, font=\tiny]{$\ldots$}; + + % Add a label for the accumulation function + \node[text width=1.5cm, align=center, font=\itshape\tiny] at ([xshift=-0.8cm, yshift=-1.2cm]acc) {accumulation of adjacent edge outputs}; +\end{tikzpicture} +\end{document} \ No newline at end of file diff --git a/docs/src/data_structure.md b/docs/src/data_structure.md index b26c721ab..fdbd40e19 100644 --- a/docs/src/data_structure.md +++ b/docs/src/data_structure.md @@ -4,19 +4,17 @@ A [`Network`](@ref) contains a list of vertex and edge models along with a graph However, in tight numerical loops, it will never access these lists of models directly. Instead, the network maintains an internal representation that tracks all symbolic indices, defining the precise ordering of states and parameters in a flat array representation. To optimize performance, especially for heterogeneous networks, the network employs specialized data structures that batch identical models together. -This disconnect between the explicit lists and the internal data structures -can be confusing. +This disconnect between the explicit lists and the internal data structures can be confusing. ## Flat Parameter and State Arrays -The vertex and edge models may contain metadata, such as initial values for states and parameters. -Crucially, this metadata is **only** for building and initialization of the -simulation. +The vertex and edge models may contain metadata, such as the initial values for states and parameters. +Crucially, this metadata is **only** for the building and initialization of the simulation. During actual simulation, the state and parameters are handled as **flat arrays**, i.e., plain `Vector{Float64}` objects. [`NWState`](@ref) and [`NWParameter`](@ref) serve as wrappers around flat arrays and the [`Network`](@ref) objects, allowing you to inspect and modify those flat arrays by addressing vertices and edges directly. -A typical workflow is: +A typical workflow is the following: 1. Set default values in the models using the metadata (see [Metadata](@ref)). 2. Create a network (see [Network Construction](@ref)). @@ -44,7 +42,7 @@ You can access the models using `getindex`/`[]` with `VIndex` or `EIndex`: ```@example data_structure v1 === nw[VIndex(1)] ``` -This can be important when changing metadata of components. i.e., both lines below are equivalent: +This can be important when changing the metadata of components. i.e., both lines below are equivalent: ```@example data_structure set_position!(v1, (1,0)) set_position!(nw[VIndex(1)], (1,0)) @@ -52,18 +50,18 @@ nothing #hide ``` !!! note "Aliasing of component models" - Since components are not copied, multiple entries in vertex and edge lists might point to the same instance of a model. + Since components are not copied, multiple entries in the vertex and edge lists might point to the same instance of a model. ```@example data_structure nw = Network(complete_graph(3), [v1,v2,v1], e) v1 === nw[VIndex(1)] === nw[VIndex(3)] ``` Consequently, metadata set for one model might affect another model. This behavior can be beneficial for performance reasons. - To force copying of components, use the `dealias` keyword: + To force the copying of components, use the `dealias` keyword: ```@example data_structure nw = Network(complete_graph(3), [v1,v2,v1], e; dealias=true) nw[VIndex(1)] === nw[VIndex(3)] # neither of them === v1 ``` -## Extracting `Network`-object from Containers -`NetworkDynamics.jl` provides a [`extract_nw`](@ref) function, to get a reference to the wrapped `Network` object from different containers, such as solution objects or integrator objects. +## Extracting a `Network`-object from Containers +`NetworkDynamics.jl` provides a [`extract_nw`](@ref) function, to get a reference to the wrapped `Network` object from different containers, such as solution objects or [integrator objects](@extref DiffEq.integrator). diff --git a/docs/src/index.md b/docs/src/index.md index d45bf9d73..02df1311b 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,40 +1,79 @@ # NetworkDynamics -A package for working with dynamical systems on complex networks. NetworkDynamics.jl provides an interface between [Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl) and [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). -It allows for high performant simulation of dynamic networks by describing the local dynamics on the edges and vertices of the graph. - -The behavior of a node or an edge can be described by algebraic equations, by differential algebraic equation (DAEs) in mass matrix form or by ordinary differential equations (ODE). - -The central construction is the function [`Network`](@ref) that receives functions describing the local dynamics on the edges and nodes of -a graph `g` as inputs, and returns a composite function compatible with the -DifferentialEquations.jl calling syntax. - -```julia -nd = Network(g, vertex_dynamics, edge_dynamics) -nd(dx, x, p, t) +The *NetworkDynamics.jl* package simulates the dynamics of complex networks. It provides an interface +between the [Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl) and the +[DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) packages and facilitates the simulation of +highly efficient dynamic networks by describing the local dynamics of the edges and vertices of the network. + +!!! note + **Complex Networks in a glance** + Complex network systems are composed by the entities that comprise them (the nodes) and the relationships that + connect each entity with one another (the edges). The mathematical structure is called + [Graph](https://en.wikipedia.org/wiki/Graph_theory) and it is used more or less interchangeably with the word + Network). + +In this graph of a simple 5-node network +```@example +using Graphs, GraphMakie, CairoMakie #hide +using GraphMakie.NetworkLayout #hide +CairoMakie.activate!(type="svg") #hide +g = smallgraph(:bull) #hide +fig, ax, p = graphplot(g; ilabels=["v$i" for i in 1:nv(g)], #hide + elabels=["e$i: $(e.src) ↦ $(e.dst)" for (i, e) in enumerate(edges(g))], #hide + layout=Align(Stress()), figure=(;resolution=(830,450))) #hide +ymin, ymax = extrema(last.(p[:node_pos][])) #hide +ylims!(ax, (ymin-0.11*(ymax-ymin), ymax+0.11*(ymax-ymin)))#hide +hidespines!(ax) #hide +hidedecorations!(ax) #hide +fig #hide ``` +we can see that it is composed of nodes (v1 to v5) who are connected to each other. The lines connecting the nodes with +each other ( e1: 1-->2, e2: 1-->3, e3: 2-->3, e4: 2-->4, e5: 3-->5) are called edges. Complex networks are composed of +multiple nodes and edges, with most nodes connected to multiple other nodes with multiple edges. + +(@Hans after rereading the text I realised that the information about the core of the package and the behaviours of the +nodes and edges does not belong in the introduction but rather in the mathematical model, so I moved it. If you are ok +with this just delete this comment) Main features: -- Clear separation of local dynamics and topology: you can easily change topology of your system or switching out dynamical components. -- High performance when working with heterogeneous models (which means heaving different local dynamics in different parts of your network). -- [Symbolic Indexing](@ref) into solutions and states: NetworkDynamics keeps track of the states of the individual subsystems. -- Different execution schemes: NetworkDynamics exploits the known inter-dependencies between components to auto parallelize execution, even on GPUs! -- Equation based models: use [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/dev/) to define local dynamics, use `NetworkDynamics.jl` to combine them into large networks! +- Clear separation of local dynamics and topology: you can easily change the topology of your system or switch out + dynamic components) +- High performance when working with heterogeneous models: you can have different local dynamics in different parts of + your network) +- [Symbolic Indexing](@ref) into solutions and states: NetworkDynamics keeps track of the states of each individual + subsystem. +- Diverse execution schemes: NetworkDynamics exploits the known interdependencies between components to auto + parallelize execution, even on GPUs! +- Equation based models: you can model local dynamics using + [ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/dev/) and them combine them into larger networks using + `NetworkDynamics.jl`! + ## Where to begin? -Check out the [Mathematical Model](@ref) to understand the underlying modelling ideas of NetworkDynamics followed by the page on [Network Construction](@ref) to learn how to implement you own models. +To learn how to implement your own models and understand the underlying modelling ideas of NetworkDynamics you should +first read the [Mathematical Model](@ref) section, followed by section [Network Construction](@ref). If you prefer to look at some concrete code first check out the [Getting Started](@ref) tutorial! ## Installation -Installation is straightforward with Julia's package manager. +1. Install Julia: +- [Julia Installation](https://julialang.org/install/) +- Find your OS and follow the instructions for the installation +2. Install NetworkDynamics.jl with Julia's package manager: ```julia-repl (v1.11) pkg> add NetworkDynamics ``` +3. Install the Julia package LiveServer: +```julia-repl +import Pkg; Pkg.add("LiveServer") +``` + +To learn more about how to use Julia you can visit: [Modern Julia Workflows](https://modernjuliaworkflows.org/) + ## Reproducibility @@ -78,7 +117,8 @@ Pkg.status(; mode = PKGMODE_MANIFEST) # hide ``` ## Funding -Development of this project was in part funded by the *German Federal Ministry for Economic Affairs and Climate Action* as part of the *OpPoDyn*-Project (Project ID 01258425/1, 2024-2027). +Development of this project was in part funded by the *German Federal Ministry for Economic Affairs and Climate Action* +as part of the *OpPoDyn*-Project (Project ID 01258425/1, 2024-2027). ```@raw html diff --git a/docs/src/mathematical_model.md b/docs/src/mathematical_model.md index 26ddff58b..8580f150a 100644 --- a/docs/src/mathematical_model.md +++ b/docs/src/mathematical_model.md @@ -1,49 +1,76 @@ # Mathematical Model -The basic mathematical model of `NetworkDynamics.jl` splits up the system into two parts: vertex and edge components. -The main goal of `NetworkDynamics.jl` is, to express the overall network dynamics as a Differential-Algebraic-Equation (DAE) +The core of the `NetworkDynamics.jl` package is the [`Network`](@ref) function. It accepts the functions describing the +local dynamics on the edges and nodes of the graph `g` as inputs, and returns a composite function compatible with the +DifferentialEquations.jl syntax as output. + +```julia +nd = Network(g, vertex_dynamics, edge_dynamics) +nd(dx, x, p, t) +``` + +In general the local dynamics on the edges and nodes of a graph can be described through the use of (a) algebraic +equations, (b) differential algebraic equation (DAEs) in mass matrix form or (c) ordinary differential equations (ODE). +The `NetworkDynamics.jl` package uses +[Differential-Algebraic-Equation (DAE)](https://mathworld.wolfram.com/Differential-AlgebraicEquation.html) to express +the overall network dynamics: ```math M\,\frac{\mathrm{d}}{\mathrm{d}t}u = f^{\mathrm{nw}}(u, p, t) ``` -where M is a (possibly singular) mass matrix, $u$ is the internal state vector of the system, $p$ are the parameters and $t$ is the time. -To make this compatible with the solvers for `OrdinaryDiffEq.jl`, the created [`Network`](@ref) object is a callable object +where $M$ is a (possibly singular) mass matrix, $u$ is the internal state vector of the system, $p$ are the parameters +and $t$ is the time. To make this compatible with the solvers used in `OrdinaryDiffEq.jl`, the generated +[`Network`](@ref) object is callable ``` -nw(du, u, p, t) # mutates du +nw(du, u, p, t) # mutates du as an "output" ``` -with stored mass matrix information to build an `ODEProblem` based on the `Network`. - -Instead of defining $f^{\mathrm{nw}}$ by hand, `ND.jl` helps you to build it automatically based on a list of decentralized nodal and edge dynamics, so-called `VertexModel` and `EdgeModel` objects. -Each component model $\mathrm c$ is modeled as general input-output-system +and represents the right-hand-side (RHS) of the equation above. The mass-matrix $m$ is stored in the `Network` object +as well. +## Modelling the Dynamics of the System +Each component model $\mathrm c$ is modeled as a general input-output-system: ```math \begin{aligned} M_{\mathrm c}\,\frac{\mathrm{d}}{\mathrm{d}t}x_{\mathrm c} &= f^{\mathrm c}(x^{\mathrm c}, i_{\mathrm c}, p_{\mathrm c}, t)\\ y^{\mathrm c} &= g^{\mathrm c}(x^\mathrm{c}, i_{\mathrm c}, p_{\mathrm c}, t) \end{aligned} ``` +where $M_{\mathrm{c}}$ is the component mass matrix, $x^{\mathrm c}$ are the component states, $i^{\mathrm c}$ are the +inputs of the component and $y^{\mathrm c}$ is the output of the component. If +$\mathrm{dim}(x^{\mathrm{c}}) = 0$, the number of internal states is 0. -where $M_{\mathrm{c}}$ is the component mass matrix, $x^{\mathrm c}$ are the component states, $i^{\mathrm c}$ are the *inputs* of the component and $y^{\mathrm c}$ is the *output* of the component. It is possible to have $\mathrm{dim}(x^{\mathrm{c}}) = 0$ and thus no internal states. +The mathematical model of `NetworkDynamics.jl` splits the network system in two parts: the vertex and +the edge components (the nodes and edges, respectively). Instead of defining the $f^{\mathrm{nw}}$ by hand, `ND.jl` +builds it automatically based on a list of decentralized nodal and edge dynamics that the user provides (the +`VertexModel` and `EdgeModel` objects). -In the network context, the **output of the edges are flow variables**. The **outputs of vertices are potential variables**. In interconnection, the *flow* on the edges depends on the *potentials* at both ends as inputs. The *potentials* of the nodes depend on the incoming *flows* from all connected edges as an input. (Here, flow and potentials are meant in a conceptional and not necessarily physical way.) -```@raw html - -``` +In the context of the network, the **output of the edges are flow variables** and the **outputs of vertices are +potential variables**. When the node and edge models are placed on a graph, the inputs and outputs are connected: +the nodes receive the output of the adjacent edges as inputs and the edges receive the output of the adjacent nodes as +inputs. Thus, the *flow* on the edges depends on the *potentials* at both ends as inputs. The *potentials* of the nodes +depend on the incoming *flows* from all connected edges as an input. (Here, flow and potentials are meant in a +conceptional and not necessarily physical way.) -## Vertex Models -Specifically, a (single-layer) vertex model has one input, and one output. -The input is an aggregation/reduction over all *incident edge outputs*, -```math -i^{\mathrm v} = \mathop{\mathrm{agg}}\limits_k^{\text{incident}} y^{\mathrm e}_k \qquad\text{often}\qquad -i^{\mathrm v} = \sum_k^{\text{incident}} y^{\mathrm e}_k +In this graphical representation of a partial network graph +```@raw html + ``` -The full vertex model +three nodes are visible (node 1, node 2 and node 3) as well as the edges connecting node 1 and node 2 ($e_{\mathrm 12}$, +$e_{\mathrm 21}$). Above the network, the mass matrix equations on node 1 and node 2 ($M_{\mathrm{c}} x_{\mathrm c}$), +the equations on the connecting edges ($e_{\mathrm 12}$, $e_{\mathrm 21}$), as well as the internal state vector +equations of node1 and node2 ($u_1$ and $u_2$) are also shown. + +(@Hans: the .svg graphics in this page looks black with dark letters for the most part in some browsers (e.g. firefox, +duckduckgo). I tried to add a white background to them, but it does not seem to have worked) + +### Vertex Models +The equations of a (single-layer) full vertex model are: ```math \begin{aligned} -M^{\mathrm v}\,\frac{\mathrm{d}}{\mathrm{d}t}x^{\mathrm v} &= f^{\mathrm v}(u^{\mathrm v}, i^{\mathrm v}, p^{\mathrm v}, t)\\ -y^{\mathrm v} &= g^{\mathrm v}(u^{\mathrm v}, i^{\mathrm v}, p^{\mathrm v}, t) +M^{\mathrm v}\,\frac{\mathrm{d}}{\mathrm{d}t}x^{\mathrm v} &= f^{\mathrm v}(x^{\mathrm v}, i^{\mathrm v}, p^{\mathrm v}, t)\\ +y^{\mathrm v} &= g^{\mathrm v}(x^{\mathrm v}, i^{\mathrm v}, p^{\mathrm v}, t) \end{aligned} ``` -corresponds to the Julia functions +and they correspond to the Julia functions: ```julia function fᵥ(dxᵥ, xᵥ, e_aggr, pᵥ, t) # mutate dxᵥ @@ -56,16 +83,36 @@ end vertf = VertexModel(; f=fᵥ, g=gᵥ, mass_matrix=Mᵥ, ...) ``` -## Edge Models -In contrast to vertex models, edge models in general have *two* inputs and *two* outputs, both for source and destination end of the edge. -We commonly use `src` and `dst` to describe the source and destination end of an edge respectively. +A (single-layer) full vertex model has one input, and one output. Its input is an aggregation/reduction over all the +*incident edge outputs* which is calculated using: +```math +i^{\mathrm v} = \mathop{\mathrm{agg}}\limits_k^{\text{incident}} y^{\mathrm e}_k \qquad\text{often}\qquad +i^{\mathrm v} = \sum_k^{\text{incident}} y^{\mathrm e}_k +``` + +The graphical representation of such a model is: +```@raw html + +``` +where $y^e_i$ and $y^e_j$ are two of the $n$ incident edge outputs that are aggregated to produce the model input +$i^u$ and the model output $y^v$ (the vertex model output). -!!! note "On the directionality of edges" - Mathematically, in a system defined on an undirected graph there is no difference between the edge $(1,2)$ and $(2,1)$, the edge has no direction. However, from an implementation point of view we always need to have some kind of ordering for function arguments, state order and so on. For undirected graphs, `Graphs.jl` chooses the direction of an edge `v1->v2` such that `v1 < v2`. + +### Edge Models +In contrast to the vertex models, edge models in general have *two* inputs and *two* outputs, for both the source and +the destination end of the edge. We commonly use `src` and `dst` to describe the source and destination end of an edge, +respectively. -The *inputs* of the edge are just the outputs of the two nodes at both ends. The output is split into two: the `dst` output goes to the input of the vertex at the destination end, the `src` output goes to the input of the vertex at the `src` end. +!!! note "Source and Destination definitions " + A source node (`src`) is the node from which the flow exits. + A destination node (`dst`) is the node into which the flow enters. -The full model of an edge +!!! note "On the directionality of edges" + Mathematically, in a system defined on an undirected graph there is no difference between edge $(1,2)$ and + edge $(2,1)$, because the edge has no direction. However, from an implementation point of view we always need to + have some kind of ordering, which is why we introduce the source (`src`) and destination (`dst`) terminology. + +The full edge model equations are: ```math \begin{aligned} M^{\mathrm e}\,\frac{\mathrm{d}}{\mathrm{d}t}x^{\mathrm e} &= f^{\mathrm e}(u^{\mathrm e}, y^{\mathrm v}_{\mathrm{src}}, y^{\mathrm v}_{\mathrm{dst}}, p^{\mathrm e}, t)\\ @@ -73,7 +120,7 @@ y^{\mathrm e}_{\mathrm{dst}} &= g_\mathrm{dst}^{\mathrm e}(u^{\mathrm e}, y^{\ma y^{\mathrm e}_{\mathrm{src}} &= g_\mathrm{src}^{\mathrm e}(u^{\mathrm e}, y^{\mathrm v}_{\mathrm{src}}, y^{\mathrm v}_{\mathrm{dst}}, p^{\mathrm e}, t) \end{aligned} ``` -corresponds to the Julia functions +and they correspond to the Julia functions: ```julia function fₑ(dxₑ, xₑ, v_src, v_dst, pₑ, t) # mutate dxᵥ @@ -86,26 +133,50 @@ end vertf = EdgeModel(; f=fₑ, g=gₑ, mass_matrix=Mₑ, ...) ``` -The sign convention for both outputs of an edge must be identical, -typically, a positive flow represents a flow *into* the connected vertex. -This is important, because the vertex only receives the flows, it does not know whether the flow was produce by the source or destination end of an edge. -``` - y_src y_dst - V_src o───←─────────→───o V_dst +The *inputs* of an edge connecting two nodes are the *outputs* of the nodes at both their ends. The output of each node +is split into two parts: +1. the *`dst`* output which is used as the input of the vertex at the destination end +2. the `src` output which is used as the input of the vertex at the `src` end. + +Because a Vertex only receives flows from the edges connected to it, it does not know whether the flow it +received was produced by the source or the destination end of an edge. To solve this problem a sign convention has been +introduced. A positive flow represents a flow *into* the connected vertex, while a negative flow represents a flow *out* +of the connected vertex. +When we consider $n_1$ as the source and $n_2$ as the destination, flows out of $n_1$ and into $n_2$ are negative, +so $y_dst < 0$. Whereas flows out of $n_2$ and into $n_1$ are positive, so $y_src > 0$. ``` + y_dst < 0 + n_1 (src) o───────→────────o n_2 (dst) + y_src > 0 + n_1 (src) o───────←────────o n_2 (dst) +``` +But when we consider $n_2$ as the source and $n_1$ as the destination, flows out of $n_2$ and into $n_1$ are negative, +so $y_dst < 0$. Whereas flows out of $n_1$ and into $n_2$ are positive, so $y_src > 0$. +``` + y_src > 0 + $n_1$ (dst) o───────→────────o $n_2$ (src) + y_dst < 0 + $n_1$ (dst) o───────→────────o $n_2$ (src) +``` + +For undirected graphs, `Graphs.jl` chooses the direction of an edge (so which node is the `src` and which the `dst`) +$v_1->v_2$ such that $v_1 < v_2$, so the edge between vertices 16 and 12 will always be an edge with source +`src=12` and destination `dst=16`. + ### Single Sided Edge Outputs -Often, edge outputs will possess some symmetry which makes it more convenient to define "single sided" edge output functions +To simplify the calculations, we split the output flows of edges into "single sided" edge output functions. That way we +can split the output flows of edges into "single sided" edge output functions which simplifies the calculations: ```julia function g_single(y, xᵥ, v_src, v_dst, pₑ, t) # mutate y nothing end ``` -There are multiple wrappers available to automaticially convert them into double-sided edge output functions: +There are multiple wrappers available to automatically convert them into double-sided edge output functions: - `Directed(g_single)` builds a double-sided function *which only couples* to the destination side. - `Symmetric(g_single)` builds a double-sided function in which both ends receive `y`. - `AntiSymmetric(g_single)` builds a double-sided function where the destination receives `y` and the source receives `-y`. @@ -113,9 +184,17 @@ There are multiple wrappers available to automaticially convert them into double ## Feed Forward Behavior -The most general version of the component models can contain direct feed forwards from the input, i.e. the edge output might depend directly on the connected vertices or the vertex output might depend directly on the aggregated edge input. +!!! warning "Feed Forward Vertices" + As of 11/2024, vertices with feed forward behaviour (FF) are not supported at all. Use [`ff_to_constraint`](@ref) to + transform them into vertex models without FF. + +Component models can have a so-called Feed Forward behaviour, which provides a direct link between the input and the +output. -Whenever possible, you should define output functions without feed forwards, i.e. +The most generic version of the component models can contain direct FFs from the input to the output. This means that +the output function $g$ depends directly on the component inputs $i$ rather than just on the component state $x$. + +Whenever possible, you should define output functions without FFs in the following way: ```julia gᵥ_noff(yᵥ, xᵥ, pᵥ, t) gₑ_noff([y_src,] y_dst, xᵥ, pₑ, t) @@ -126,16 +205,16 @@ gᵥ(yᵥ, xᵥ, e_aggr, pᵥ, t) gₑ([y_src], y_dst, xᵥ, v_src, v_dst, pₑ, t) ``` -NetworkDynamics cannot couple two components with feed forward to each other. -It is always possible to transform feed forward behavior to an internal state `x` with mass matrix entry zero to circumvent this problem. This transformation can be performed automatically by using [`ff_to_constraint`](@ref). +NetworkDynamics cannot couple two components with FFs to each other. But, it is always possible to transform +feed forward behaviour to an internal state `x` with mass matrix entry zero to circumvent this problem. This +transformation can be performed automatically using [`ff_to_constraint`](@ref). -!!! warning "Feed Forward Vertices" - As of 11/2024, vertices with feed forward are not supported at all. Use [`ff_to_constraint`](@ref) to transform them into vertex model without FF. +Concretely, NetworkDynamics distinguishes between 4 types of feed forward behaviours of `g` functions based on the +[`FeedForwardType`](@ref) trait. +The feed forward type is inferred automatically based on the provided function `g` (by calculating the signature of the +methods used over the number of arguments it is given.) -Concretely, NetworkDynamics distinguishes between 4 types of feed forward behaviors of `g` functions based on the [`FeedForwardType`](@ref) trait. -The different types the signature of provided function `g`. -Based on the signatures avaialable, ND.jl will try to find the correct type automaticially. Using the `ff` -keyword in the constructors, the user can enforce a specific type. +The code bloke below presents the different `g` signatures for the different feed forward types: **[`PureFeedForward()`](@ref)** ```julia @@ -166,3 +245,5 @@ g!(out_src, out_dst, x) # double-sided edge g!(v_out, x) # single layer vertex ``` +If the automatic inference of feed forward type fails, the user may specify it explicitly using the `ff` keyword +argument of the Edge/VertexModel constructor. diff --git a/docs/src/network_construction.md b/docs/src/network_construction.md index f4d75f4ba..da096f4d6 100644 --- a/docs/src/network_construction.md +++ b/docs/src/network_construction.md @@ -2,9 +2,9 @@ ## Building a Network The main type of `NetworkDynamics.jl` is a [`Network`](@ref). -A network bundles various component models (edge and vertex models) together with a graph to form a callable object which represents the RHS of the overall dynamical system, see [Mathematical Model](@ref). +A network bundles various component models (edge and vertex models) together with a graph to form a callable object which represents the right hand side (RHS) of the overall dynamical system, see [Mathematical Model](@ref). -A `Network` is build by passing a graph `g`, vertex models `vertexm` and edge models `edgem`. +A `Network` is build by passing a graph `g`, vertex models `vertexm` and edge models `edgem` to the [`Network`](@ref) constructor:. ```julia nw = Network(g, vertexm, edgem; kwargs...) ``` @@ -13,20 +13,17 @@ Two important keywords for the [`Network`](@ref) constructor are: - `execution`: Defines the [`ExecutionStyle`](@ref) of the coreloop, e.g. `SequentialExecution{true}()`. - A execution style is a special struct which tells the backend how to parallelize for example. + A execution style is a special Julia object, which tells the backend how to parallelize (e.g. `ThreadedExecution{true}()` will use native Julia threads to parallelize the RHS call). A list of available executions styles can be found under [Execution Types](@ref) in the API. -- `aggregator`: - Tells the backend how to aggregate and which aggregation function to use. - Aggregation is the process of creating a single vertex input by reducing over - the outputs of adjecent edges of said vertex. The `aggregator` contains both the - function and the algorithm. E.g. `SequentialAggregator(+)` is a sequential - aggregation by summation. A list of availabe Aggregators can be found under - [`Aggregators`](@ref) in the API. +- `aggregator`: + Instructs the backend how to perform the aggregation and which aggregation function to use. + Aggregation is the process of creating a single vertex input by reducing over the outputs of adjecent edges of said vertex. The `aggregator` contains both the function and the algorithm. E.g. `SequentialAggregator(+)` is a sequential aggregation by summation. A list of availabe Aggregators can be found under [`Aggregators`](@ref) in the API. ## Building `VertexModel`s -This chapter walks through the most important aspects when defining custom vertex model. For a list of all keyword arguments please check out the docstring of [`VertexModel`](@ref). -As an example, we'll construct an second order kuramoto model, because that's what we do. +This chapter will walk you through the most important aspects of defining a custom vertex model. For a list of all keyword arguments please check out the docstring of [`VertexModel`](@ref). + +As an example, we'll construct an second order kuramoto model, because that is what the package does. ```@example construction using NetworkDynamics #hide function kuramoto_f!(dv, v, esum, p, t) @@ -51,20 +48,21 @@ function kuramoto_g_noff!(y, v, p, t) end VertexModel(; f=kuramoto_f!, g=kuramoto_g_noff!, dim=2, pdim=3, outdim=1) ``` -It is still annoying to explicitly write this trivial output function. You can prevent this by using [`StateMask`](@ref). + +To simplify your programming and avoid explicitly writing the above trivial output function you can use [`StateMask`](@ref). By writing ```@example construction VertexModel(; f=kuramoto_f!, g=StateMask(1:1), dim=2, pdim=3) ``` -we told the vertex model, that the output is part of the states `x[1:1]`. -This enables a few things: -- `outdim` is not needed anymore, can be inferred from `StateMask` +we are instructing the vertex model, that the output is part of the states `x[1:1]`. +This results in the following changes: +- `outdim` is removed because it can be inferred from `StateMask` - `outsym` is not a generic `:o` any more but inferred from the state symbols. We can be even less verbose by writing `g=1:1` or just `g=1`. -In a last we define better names for our states and parameters as well as assigning a position in the graph to enable the graphless network construction. -Whenever you provide `sym` keyword the corresponding `dim` keyword is not neccessary anymore. We end up with a relatively short definition +Lastly, we define improved names for our states and parameters as well as assigning a position in the graph to enable the graphless network construction. +Whenever you provide a `sym` keyword the corresponding `dim` keyword stops being neccessary. So, we end up with a relatively short definition ```@example construction VertexModel(; f=kuramoto_f!, g=1, sym=[:θ, :ω], psym=[:M=>1, :P=>0.1, :D=>0], @@ -72,9 +70,9 @@ VertexModel(; f=kuramoto_f!, g=1, ``` ## Building `EdgeModel`s -This chapter walks through the most important aspects when defining custom edge models. For a list of all keyword arguments please check out the docstring of [`EdgeModel`](@ref). +This chapter walks you through the most important aspects when defining custom edge models. For a list of all keyword arguments please check the docstring of [`EdgeModel`](@ref). -As an example edge model we want to define standard sinusoidal coupling between the vertices in our network. The full definition looks like this: +As an example edge model we define a standard sinusoidal coupling between the vertices in our network. The full definition is: ```@example construction function edge_f!(de, e, vsrc, vdst, p, t) @@ -95,15 +93,15 @@ function edge_g_ff!(ysrc, ydst, vsrc, vdst, p, t) end EdgeModel(;g=edge_g_ff!, pdim=1, outdim=1) ``` -which no classifies as a `PureFeedForward` edge. -In cases like this, where the edge is actually anti symmetric we can alternatively define a single sided output function and wrapping it in an `AntiSymmetric` object +which classifies as a `PureFeedForward` edge. +In cases like this, where the edge is actually anti-symmetrical we can define a single sided output function and wrap it in an `AntiSymmetric` object: ```@example construction function edge_g_s!(ydst, vsrc, vdst, p, t) ydst[1] = p[1] * sin(vsrc[1] - vdst[1]) end EdgeModel(;g=AntiSymmetric(edge_g_ff!), pdim=1, outdim=1) ``` -which can also lead to briefer output naming. Available single sided wrappers are +This can also lead to briefer output naming. Available single sided wrappers are: - [`Directed`](@ref) (no coupling at `src`), - [`AntiSymmetric`](@ref) (same coupling at `src` and `dst`), - [`Symmetric`](@ref) (inverse coupling at `dst`) and diff --git a/docs/src/symbolic_indexing.md b/docs/src/symbolic_indexing.md index c3cc8cbb1..f2d275f35 100644 --- a/docs/src/symbolic_indexing.md +++ b/docs/src/symbolic_indexing.md @@ -1,24 +1,21 @@ # Symbolic Indexing -Using SciML's [`SymbolicIndexingInterface.jl`](https://github.com/SciML/SymbolicIndexingInterface.jl), `ND.jl` provides numerous methods to access and change variables and parameters. - -!!! details "Setup code to make following examples work" - ```@example si - using NetworkDynamics - using Graphs - using OrdinaryDiffEqTsit5 - using Plots - ``` +By using SciML's [`SymbolicIndexingInterface.jl`](https://github.com/SciML/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. ```@example si +using NetworkDynamics +using Graphs +using OrdinaryDiffEqTsit5 +using 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]) ``` -Here we created a static diffusion edge with suitable variable and parameter names. +Here we create a static diffusion edge with suitable variable and parameter names. Similarly, we define the diffusion vertex with symbolic names. ```@example si function _vertexf!(dv, v, esum, p, t) @@ -30,9 +27,11 @@ vertexf = VertexModel(f=_vertexf!, g=1, sym=[:storage]) ## Fundamental Symbolic Indices The default types for this access are the types [`VIndex`](@ref), [`EIndex`](@ref), [`VPIndex`](@ref) and [`EPIndex`](@ref). -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. +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)` refers to the *second* parameter of the edge component for the 4th edge. !!! details "Setup code to make following examples work" ```@example si @@ -47,13 +46,13 @@ As such, `VIndex(2, :x)` refers to variable with symbolic name `:x` in vertex nu Those fundamental indices can be used in a lot of scenarios. Most importantly you can use them to ```@example si -sol(sol.t; idxs=VIndex(1, :storage)) # extract timeseries out ouf solution object +sol(sol.t; idxs=VIndex(1, :storage)) # extract timeseries out of a 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`](@ref), [`eidxs`](@ref), [`vpidxs`](@ref) and [`epidxs`](@ref). -With the help of those methods you can generate arrays of symbolic indices: +Often, you need many individual symbolic indices. To achieve this you can use the helper methods [`vidxs`](@ref), +[`eidxs`](@ref), [`vpidxs`](@ref) and [`epidxs`](@ref). With their help you can generate arrays of symbolic indices: ```@example si vidxs(nw, :, :storage) # get variable "storage" for all vertices @@ -98,7 +97,8 @@ in these flat arrays corresponds exactly to the order returned by [`variable_sym !!! note The `NWState` and `NWParameter` wrappers can be constructed from various objects. - Fore example, within a callback you might construct `p = NWParameter(integrator)` to then change the parameters of the network within the callback. + For example, within a callback you might construct `p = NWParameter(integrator)` to then change the parameters of +the network within the callback. ## Observables @@ -109,10 +109,11 @@ 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 +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: ```@example si plot(sol; idxs=eidxs(nw, :, :flow)) ``` @@ -128,7 +129,8 @@ For example, we can directly plot the storage difference with respect to storage 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. +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) ``` @@ -150,7 +152,7 @@ idxs = SII.parameter_index(nw, eidxs(1:2, :K)) pflat(s)[idxs] == s.p.e[1:2, :K] ``` -If you need the symbols of all the states/parameters in order, you can use +If you need the symbols of all the states/parameters in order, you can use: ```@example si SII.variable_symbols(nw) ```