extend for vector valued obsex

also add docs

Author: hexaeder

docs/src/API.md

@@ -82,6 +82,7 @@ VIndex
 EIndex
 VPIndex
 EPIndex
+@obsex
 ```
 
 ### Index generators

docs/src/symbolic_indexing.md

@@ -114,3 +114,20 @@ Observables can be accessed like any other state, for example, the flows in the
 ```@example si
 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.
+
+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.
+```
+@obsex mag = sqrt(VIndex(1, :u_r)^2 + VIndex(2, :u_i)^2)
+```
+

ext/SymbolicsExt.jl

@@ -15,31 +15,38 @@ function NetworkDynamics.generate_observable_expression(ex::Expr)
     # manually hygen on "mapping", otherwise we cannot escape x
     mapping_sym = gensym("mapping")
     symbolic_expr = postwalk(content) do x
-        :(NetworkDynamics.collect_symbol!($mapping_sym, $x))
+        if x ∈ (:+, :-, :^, :/, :.+, :.-, :.^, :./, Symbol(":")) || x isa QuoteNode
+            x
+        else
+            :(NetworkDynamics.collect_symbol!($mapping_sym, $x))
+        end
     end
     quote
         $(esc(mapping_sym)) = Dict()
         expr = $(esc(symbolic_expr))
-        ObservableExpression($(esc(mapping_sym)), expr, $(Meta.quot(name)))
+        NetworkDynamics.ObservableExpression($(esc(mapping_sym)), expr, $(Meta.quot(name)))
     end
 end
 
+function NetworkDynamics.ObservableExpression(mapping, expr::Vector, names)
+    if names isa Symbol
+        names = [Symbol(names, NetworkDynamics.subscript(i)) for i in 1:length(expr)]
+    elseif isnothing(names)
+        names = Iterators.repeated(nothing, length(expr))
+    end
+    [NetworkDynamics.ObservableExpression(mapping, e, n) for (e, n) in zip(expr, names)]
+end
+
 function NetworkDynamics.ObservableExpression(mapping, expr, name)
     nwidx = collect(keys(mapping))
     syms = collect(values(mapping))
     f = Symbolics.build_function(expr, syms; expression=Val(false))
-    ObservableExpression(nwidx, f, expr, name)
+    NetworkDynamics.ObservableExpression(nwidx, f, expr, name)
 end
 
 NetworkDynamics.collect_symbol!(_, x) = x
-function NetworkDynamics.collect_symbol!(mapping, x::AbstractVector)
-    if length(x) == 1 && only(x) isa Union{SymbolicVertexIndex,SymbolicEdgeIndex}
-        return NetworkDynamics.collect_symbol!(mapping, only(x))
-    elseif any(el -> el isa Union{SymbolicVertexIndex,SymbolicEdgeIndex}, x)
-        throw(ArgumentError("Cannot handle vector expressions in @obsex! Encountered $x."))
-    else
-        x
-    end
+function NetworkDynamics.collect_symbol!(mapping, x::Union{AbstractVector, Tuple})
+    map(el -> NetworkDynamics.collect_symbol!(mapping, el), x)
 end
 function NetworkDynamics.collect_symbol!(mapping, nwindex::Union{SymbolicVertexIndex,SymbolicEdgeIndex})
     if haskey(mapping, nwindex)

src/symbolicindexing.jl

@@ -1148,6 +1148,15 @@ function SII.observed(nw::Network, obsex::ObservableExpression)
         obsex.f(input)
     end
 end
+function SII.observed(nw::Network, obsexs::AbstractVector{<:ObservableExpression})
+    inputfs = map(obsex -> SII.observed(nw, obsex.inputs), obsexs)
+    (u, p, t) -> begin
+        map(obsexs, inputfs) do obsex, inputf
+            input = inputf(u, p, t)
+            obsex.f(input)
+        end
+    end
+end
 
 Base.getindex(s::NWState, idx::ObservableExpression) = SII.getu(s, idx)(s)
 Base.getindex(s::NWParameter, idx::ObservableExpression) = SII.getp(s, idx)(s)

test/symbolicindexing_test.jl

@@ -505,4 +505,8 @@ end
 
     @test SII.getname(@obsex(VIndex(1,:δ) + VIndex(2,:δ))) == Symbol("v1₊δ+v2₊δ")
     @test SII.getname(@obsex(δ²=VIndex(1,:δ)^2)) == :δ²
+
+    @test s[@obsex(vidxs(s, :, :δ) .- VIndex(1, :δ))] == [0, s[VIndex(2,:δ)] - s[VIndex(1,:δ)], s[VIndex(3,:δ)] - s[VIndex(1,:δ)]]
+
+    obsex = @obsex(δ_rel = vidxs(s, :, :δ) .- VIndex(1, :δ))
 end