Merge pull request #205 from JuliaDynamics/hw/unused_p

keep track of unused parameters

Author: hexaeder

.github/workflows/benchmark.yml

@@ -27,8 +27,10 @@ jobs:
                            "https://api.github.com/repos/${{ github.repository }}/pulls/${{ github.event.pull_request.number }}")
           PR_TITLE=$(echo "$PR_DATA" | jq -r '.title')
           PR_BODY=$(echo "$PR_DATA" | jq -r '.body')
+          # Ensure PR_BODY is single-line and sanitized for GitHub env vars
+          PR_BODY_ESCAPED=$(echo "$PR_BODY" | tr '\n' ' ' | tr '\r' ' ')
           echo "PR_TITLE=$PR_TITLE" >> $GITHUB_ENV
-          echo "PR_BODY=$PR_BODY" >> $GITHUB_ENV
+          echo "PR_BODY=$PR_BODY_ESCAPED" >> $GITHUB_ENV
       - name: Check for '[no benchmark]'
         id: check
         run: |

ext/MTKExt.jl

@@ -55,7 +55,14 @@ function VertexModel(sys::ODESystem, inputs, outputs; verbose=false, name=getnam
     sym = [s => _get_metadata(sys, s) for s in _sym]
 
     _psym = getname.(gen.params)
-    psym = [s => _get_metadata(sys, s) for s in _psym]
+    psym = map(gen.params) do p
+        pname = getname(p)
+        md = _get_metadata(sys, pname)
+        if p in gen.unused_params
+            md = (; md..., unused=true)
+        end
+        pname => md
+    end
 
     _obssym = getname.(gen.obsstates)
     obssym = [s => _get_metadata(sys, s) for s in _obssym]
@@ -152,7 +159,14 @@ function EdgeModel(sys::ODESystem, srcin, dstin, srcout, dstout; verbose=false,
     sym = [s => _get_metadata(sys, s) for s in _sym]
 
     _psym = getname.(gen.params)
-    psym = [s => _get_metadata(sys, s) for s in _psym]
+    psym = map(gen.params) do p
+        pname = getname(p)
+        md = _get_metadata(sys, pname)
+        if p in gen.unused_params
+            md = (; md..., unused=true)
+        end
+        pname => md
+    end
 
     _obssym = getname.(gen.obsstates)
     obssym = [s => _get_metadata(sys, s) for s in _obssym]
@@ -373,12 +387,11 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     fftype = _determine_fftype(out_deps, states, allinputs, params, iv)
 
     # filter out unnecessary parameters
-    # var_deps = _all_rhs_symbols(eqs)
-    # used_params = params ∩ (var_deps ∪ obs_deps ∪ out_deps)
-    # if Set(params) != Set(used_params)
-    #     verbose && @info "Remove parameters $(collect(setdiff(params, used_params))) which arn't used in the equations."
-    #     params = used_params
-    # end
+    var_deps = _all_rhs_symbols(eqs)
+    unused_params = Set(setdiff(params, (var_deps ∪ out_deps))) # do not exclud obs_deps
+    if verbose && !isempth(unused_params)
+        @info "Parameters $(unused_params) do not appear in equations of f and g and will be marked as unused."
+    end
 
     # TODO: explore Symbolcs/SymbolicUtils CSE
     # now generate the actual functions
@@ -412,19 +425,22 @@ function generate_io_function(_sys, inputss::Tuple, outputss::Tuple;
     obsformulas = [eq.rhs for eq in obseqs]
     _, obsf_ip = build_function(obsformulas, states, inputss..., params, iv; cse=true, expression)
 
-    return (;f=f_ip, g=g_ip,
-            mass_matrix,
-            states,
-            inputss,
-            outputss,
-            obsstates,
-            fftype,
-            obsf = obsf_ip,
-            equations=formulas,
-            outputeqs=Dict(Iterators.flatten(outputss) .=> gformulas),
-            observed=Dict(getname.(obsstates) .=> obsformulas),
-            odesystem=sys,
-            params)
+    return (;
+        f=f_ip, g=g_ip,
+        mass_matrix,
+        states,
+        inputss,
+        outputss,
+        obsstates,
+        fftype,
+        obsf = obsf_ip,
+        equations=formulas,
+        outputeqs=Dict(Iterators.flatten(outputss) .=> gformulas),
+        observed=Dict(getname.(obsstates) .=> obsformulas),
+        odesystem=sys,
+        params,
+        unused_params
+    )
 end
 
 function _determine_fftype(deps, states, allinputs, params, t)

src/callbacks.jl

@@ -37,7 +37,7 @@ Creates a callback condition for a [`ComponentCallback`].
 Consider a component model with states `[:u1, :u2]`, inputs `[:i]`, outputs
 `[:o]` and parameters `[:p1, :p2]`.
 
-    ComponentCondition([:u1, :o], [:p1]) do (u, p, t)
+    ComponentCondition([:u1, :o], [:p1]) do u, p, t
         # access states symbolicially or via int index
         u[:u1] == u[1]
         u[:o] == u[2]
@@ -81,7 +81,7 @@ Creates a callback condition for a [`ComponentCallback`].
 Consider a component model with states `[:u1, :u2]`, inputs `[:i]`, outputs
 `[:o]` and parameters `[:p1, :p2]`.
 
-    ComponentAffect([:u1, :o], [:p1]) do (u, p, ctx)
+    ComponentAffect([:u1, :o], [:p1]) do u, p, ctx
         u[:u1] = 0 # change the state
         p[:p1] = 1 # change the parameter
         @info "Changed :u1 and :p1 on vertex \$(ctx.vidx)" # access context

src/initialization.jl

@@ -60,7 +60,8 @@ function initialization_problem(cf::T; t=NaN, apply_bound_transformation, verbos
     infree_ms = Tuple((!).(map(s -> has_default(cf, s), sv)) for sv in insym_normalized(cf))
     infixs = Tuple(Float64[has_default(cf, s) ? get_default(cf, s) : NaN for s in sv] for sv in insym_normalized(cf))
 
-    pfree_m = (!).(map(Base.Fix1(has_default, cf), psym(cf)))
+    is_free_p = s -> !is_unused(cf, s) && !has_default(cf, s) # free p are not unused and have no default
+    pfree_m = map(is_free_p, psym(cf))
     pfix = Float64[has_default(cf, s) ? get_default(cf, s) : NaN for s in psym(cf)]
 
     # count free variables and equations
@@ -74,7 +75,7 @@ function initialization_problem(cf::T; t=NaN, apply_bound_transformation, verbos
                    psym(cf)[pfree_m])
     @assert length(freesym) == Nfree
     if Neqs < Nfree
-        throw(ArgumentError("Initialization problem underconstraint. $(Neqs) Equations and  $(Nfree) variables."))
+        throw(ArgumentError("Initialization problem underconstraint. $(Neqs) Equations and  $(Nfree) variables: $freesym"))
     end
 
     # which parts of the nonlinear u (unl) correspond to which free parameters
@@ -282,7 +283,7 @@ function initialize_component!(cf; verbose=true, apply_bound_transformation=true
 end
 
 function isinitialized(cf::ComponentModel)
-    all(has_default_or_init(cf, s) for s in vcat(sym(cf), psym(cf)))
+    all(has_default_or_init(cf, s) || is_unused(cf, s) for s in vcat(sym(cf), psym(cf)))
 end
 
 """
@@ -304,7 +305,7 @@ function init_residual(cf::T; t=NaN, recalc=false) where {T<:ComponentModel}
         outs = Tuple(Float64[get_default_or_init(cf, s) for s in sv] for sv in outsym_normalized(cf))
         u = Float64[get_default_or_init(cf, s) for s in sym(cf)]
         ins = Tuple(Float64[get_default_or_init(cf, s) for s in sv] for sv in insym_normalized(cf))
-        p = Float64[get_default_or_init(cf, s) for s in psym(cf)]
+        p = Float64[is_unused(cf, s) ? NaN : get_default_or_init(cf, s) for s in psym(cf)]
         res = zeros(dim(cf) + sum(outdim_normalized(cf)))
 
         res_fg  = @views res[1:dim(cf)]

src/metadata.jl

@@ -72,6 +72,18 @@ for md in [:default, :guess, :init, :bounds]
 end
 set_graphelement!(c::EdgeModel, p::Pair) = set_graphelement!(c, (;src=p.first, dst=p.second))
 
+"""
+    is_unused(c::ComponentModel, sym::Symbol)
+
+Checks if symbol `sym` is marked as unused (i.e. it does not appear in the equations of f and g explicitly).
+"""
+function is_unused(c::ComponentModel, sym::Symbol)
+    if has_metadata(c, sym, :unused)
+        get_metadata(c, sym, :unused)
+    else
+        false
+    end
+end
 
 #### default or init
 """
@@ -318,7 +330,13 @@ end
 function _append_states!(lns, cf, syms; sigdigits)
     fidx = length(lns)+1
     for sym in syms
-        str = "  &" * string(sym) * " &&= "
+        str = "  &"
+        if is_unused(cf, sym)
+            str *= styled"{gray:$(string(sym)) (unused)}"
+        else
+            str *= string(sym)
+        end
+        str *= " &&= "
         if has_default_or_init(cf, sym)
             val = get_default_or_init(cf, sym)
             val_str = str_significant(val; sigdigits, phantom_minus=true)

test/initialization_test.jl

@@ -48,6 +48,7 @@ end
             V=sqrt(u_r^2 + u_i^2), [description="Voltage magnitude"]
             ω_ref=0, [description="Reference frequency"]
             Pm, [guess=0.1,description="Mechanical Power"]
+            aux, [description="Auxiliary, unused parameter"]
         end
         @equations begin
             Dt(θ) ~ ω - ω_ref
@@ -66,6 +67,8 @@ end
     @test vf.symmetadata[:θ][:bounds] == [-π, π]
     @test vf.symmetadata[:i_r][:default] == 1
     @test vf.symmetadata[:i_i][:default] == 0.1
+    @test filter(p->NetworkDynamics.is_unused(vf, p), psym(vf)) == [:aux]
+    @test NetworkDynamics.has_default_or_init(vf, :aux) == false
 
     NetworkDynamics.initialize_component!(vf; verbose=true)
     @test NetworkDynamics.init_residual(vf) < 1e-8
@@ -77,6 +80,7 @@ end
     set_default!(vf, :ω, get_init(vf, :ω))
 
     NetworkDynamics.initialize_component!(vf; verbose=true)
+    dump_initial_state(vf)
 end
 
 @testset "test component initialization with bounds" begin

test/runtests.jl

@@ -13,6 +13,8 @@ using ExplicitImports
 
 (isinteractive() ? includet : include)(joinpath(pkgdir(NetworkDynamics, "test", "testutils.jl")))
 
+haskey(ENV, "BUILDKITE") && @test CUDA.functional() # fail early in buildkite if cuda is not available
+
 @testset "NetworkDynamics Tests" begin
     @testset "Package Quality Tests" begin
         # print_explicit_imports(NetworkDynamics)