Merge pull request #228 from JuliaDynamics/hw/init_tutorial

improve docs

Author: hexaeder

docs/examples/gas_network.jl

@@ -1,5 +1,5 @@
 #=
-# Dynamic Flow in Simple Gas Network
+# [Dynamic Flow in Simple Gas Network](@id gas-example)
 
 This example is based on the paper
 

docs/examples/init_tutorial.jl

@@ -2,7 +2,10 @@
 # [Tutorial on Stepwise Initialization of a Complex Model](@id init-tutorial)
 
 This example demonstrates how to initialize a complex network model with both static
-and dynamic components. We'll create a simple gas network model with three nodes
+and dynamic components.
+The models are closely related to the ones used in the [gas network example](@ref gas-example),
+but greatly simplified for the sake of this tutorial.
+We'll create a gas network model with three nodes
 and pipes connecting them, and show how to:
 
 1. Create static models for initialization
@@ -361,14 +364,3 @@ let
 
     fig
 end
-
-#=
-## Interactive Visualization
-
-You can also visualize the results interactively using NetworkDynamicsInspector:
-
-```julia
-using NetworkDynamicsInspector
-inspect(sol)
-```
-=#

docs/make.jl

@@ -45,6 +45,7 @@ kwargs = (;
         "General" => "index.md",
         "mathematical_model.md",
         "network_construction.md",
+        "data_structure.md",
         "Features" => [
             "symbolic_indexing.md",
             "metadata.md",
@@ -57,14 +58,12 @@ kwargs = (;
         "API.md",
         "Tutorials" => [
             "Getting Started" => "generated/getting_started_with_network_dynamics.md",
-            "Directed and Weighted Graphs" => "generated/directed_and_weighted_graphs.md",
             "Heterogeneous Systems" => "generated/heterogeneous_system.md",
-            # "Stochastic differential equations" => "generated/StochasticSystem.md",
-            # "Delay differential equations" => "generated/kuramoto_delay.md",
+            "Initialization" => "generated/init_tutorial.md",
             "Cascading Failure" => "generated/cascading_failure.md",
-            "Stress on Truss" => "generated/stress_on_truss.md",
             "Gas Network" => "generated/gas_network.md",
-            "Initialization" => "generated/init_tutorial.md",
+            "Stress on Truss" => "generated/stress_on_truss.md",
+            "Directed and Weighted Graphs" => "generated/directed_and_weighted_graphs.md",
         ]
     ],
     draft=false,

docs/src/API.md

@@ -191,6 +191,7 @@ KAAggregator
 save_parameters!
 ff_to_constraint
 Base.copy(::NetworkDynamics.ComponentModel)
+extract_nw
 ```
 
 ## NetworkDynamicsInspector API

docs/src/Multithreading.md

@@ -0,0 +1,69 @@
+# Data Structure
+
+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.
+
+## 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.
+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:
+
+1. Set default values in the models using the metadata (see [Metadata](@ref)).
+2. Create a network (see [Network Construction](@ref)).
+3. Generate a state `s = NWState(nw)` which will be prefilled with the default values from the component metadata (see [SymbolicIndexing](@ref)).
+4. Change the values of `s`, i.e., `s.v[1,:x] = 1.0`: This changes the **underlying flat array** but not the metadata of the models.
+5. Build a problem with the updated flat arrays using `uflat(s)` and `pflat(s)`.
+
+## Accessing Components
+Per default, the models are not copied on Network construction:
+
+```@example data_structure
+using NetworkDynamics # hide
+using Graphs #hide
+include(joinpath(pkgdir(NetworkDynamics), "test", "ComponentLibrary.jl")) # hide
+kuramoto_first = Lib.kuramoto_vertex! # hide
+kuramoto_secnd = Lib.kuramoto_inertia! # hide
+kuramoto_edge = Lib.kuramoto_edge! # hide
+
+v1 = VertexModel(f=kuramoto_first, sym=[:θ], psym=[:ω], g=1)
+v2 = VertexModel(f=kuramoto_secnd, sym=[:δ, :ω], psym=[:M, :D, :Pm], g=1)
+e = EdgeModel(;g=AntiSymmetric(kuramoto_edge), outsym=[:P], psym=[:K])
+nw = Network(complete_graph(2), [v1, v2], e)
+```
+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:
+```@example data_structure
+set_position!(v1, (1,0))
+set_position!(nw[VIndex(1)], (1,0))
+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. 
+    ```@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:
+    ```@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.

docs/src/data_structure.md

@@ -28,7 +28,7 @@ Special cases for symbol metadata are:
 - `default`: Stores default values for states/parameters. In initialization, those are considered fixed.
 - `guess`: Stores a guess for a state/parameter which needs to solved during initialization ("free" variables).
 - `bounds`: Store bounds for variables/parameters
-- `init`: Stores the solution of the "free" variables during initialization.
+- `init`: Stores the solution of the "free" variables, this is rarely set manual but instead when calling [`initialize_component!`](@ref).
 
 Fore those, there are special functions `has_*`, `get_*` and `set_*!`. See [Per Symbol Metadata API](@ref).
 

docs/src/metadata.md

@@ -36,7 +36,7 @@ as tearing of thousands of equations.
 
 ## RC-Circuit Example
 In good [MTK tradition](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/acausal_components/), this feature will be explained along a simple RC circuit example.
-The [Dynamic Flow in Simple Gas Network](@ref) example is another showcase of the MTK constructors.
+The [Gas Network Example](@ref gas-example) or [Initialization Tutorial](@ref init-tutorial) also showcase the MTK constructors.
 
 The system to model is 2 node, 1 edge network. The node output states are the voltage (to ground), the edge output sates are the currents at both ends.
 ```

docs/src/mtk_integration.md

@@ -1,6 +1,6 @@
 # Symbolic Indexing
 
-Using SciML's [`SymblicIndexingInterface.jl`](https://github.com/SciML/SymbolicIndexingInterface.jl), `ND.jl` provides lots of methods to access and change variables and Parameters.
+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
@@ -11,7 +11,7 @@ Using SciML's [`SymblicIndexingInterface.jl`](https://github.com/SciML/SymbolicI
     ```
 
 ## Provide Symbol Names
-When construction component models, you can pass symbolic names using the `sym` and `psym` keywords.
+When constructing component models, you can pass symbolic names using the `sym` and `psym` keywords.
 ```@example si
 function _edgef!(e, v_s, v_d, (K,), t)
     e .= K * (v_s[1] .- v_d[1])
@@ -28,9 +28,9 @@ vertexf = VertexModel(f=_vertexf!, g=1, sym=[:storage])
 ```
 
 
-## Fundamental Symblic Indices
+## 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 componen and a reference to the symbol within that component.
+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.
 
@@ -56,7 +56,7 @@ Often, you need many individual symbolic indices. For that there are the helper
 With the help of those methods you can generate arrays of symbolic indices:
 
 ```@example si
-vidxs(nw, :, :storage) # get variable "storage" for all nodes
+vidxs(nw, :, :storage) # get variable "storage" for all vertices
 ```
 ```@example si
 plot(sol; idxs=vidxs(nw, :, :storage))
@@ -70,7 +70,7 @@ p = NWParameter(nw)
 ```
 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 the symbolic indices.
+The parameters in the `NWParameter` object can be accessed using symbolic indices.
 ```@example si
 p[EPIndex(5, :K)] = 2.0 # change the parameter K of the 5th edge
 nothing #hide
@@ -79,20 +79,20 @@ Similarly, you can create a `NWState` object for the network `nw` using
 ```@example si
 s = NWState(nw)
 ```
-No default values were provided in the network components, so the state array is filled with `NaN`s.
+No default values were provided in the network components, so the state array is filled with `NaN` values.
 ```@example si
-s[VIndex(:, :storage)] .= randn(5) # set the (initial) storage for alle nodes 
+s[VIndex(:, :storage)] .= randn(5) # set the (initial) storage for all vertices 
 s #hide
 ```
-For both `NWState` and `NWParameter` objects, the there is a more convenient way to access the variables and parameters.
+For both `NWState` and `NWParameter` objects, there is a more convenient way to access the variables and parameters.
 ```@example si
 @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 allways access the flat representations by calling [`uflat`](@ref) and [`pflat`](@ref).
+You can always access the flat representations by calling [`uflat`](@ref) and [`pflat`](@ref).
 
 !!! note
     The `NWState` and `NWParameter` wrappers can be constructed from various objects.
@@ -101,8 +101,8 @@ You can allways access the flat representations by calling [`uflat`](@ref) and [
 
 ## 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 those states are called Observables.
+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
@@ -117,17 +117,43 @@ plot(sol; idxs=eidxs(nw, :, :flow))
 
 ## Derived `ObservableExpressions` using `@obsex`
 
-Sometimes it is usefull to plot or observe some simple derived quantity.For that,
-one can used the [`@obsex`](@ref) macro, to define simple derived quantities.
+Sometimes it is useful to plot or observe simple derived quantities. For that,
+one can use the [`@obsex`](@ref) macro to define simple derived quantities.
 
 For example, we can directly plot the storage difference with respect to storage of node 1.
 
 ```@example si
 plot(sol; idxs=@obsex(vidxs(nw,:,:storage) .- VIndex(1,:storage)))
 ```
 
-Other examples are the calculation of magnitude and argument of complex values which are modeld in real and imaginary part.
+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`:
+
+```@example si
+using NetworkDynamics: SII # SII = SymbolicIndexingInterface
+idxs = SII.variable_index(nw, vidxs(1:2, :storage))
+```
+```@example si
+uflat(s)[idxs] == s.v[1:2, :storage]
+```
+Analogous with parmeters:
+```@example si
+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
+```@example si
+SII.variable_symbols(nw)
+```
+and
+```@example si
+SII.parameter_symbols(nw)
+```
+All above examples also work on other "symbolic containers", e.g. `SII.variable_symbols(::NWState)`.