From 3fb0339db9c638eab4f997b64028f4798d554597 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 8 Jun 2024 12:16:41 -0400
Subject: [PATCH 1/4] save progress

---
 .../bifurcation_kit_extension.jl              |   2 +-
 .../homotopy_continuation_extension.jl        |   2 +-
 src/reactionsystem_conversions.jl             |   2 +-
 .../serialise_reactionsystem.jl               |   6 +-
 src/steady_state_stability.jl                 |   2 +-
 test/dsl/dsl_options.jl                       |   4 +-
 test/extensions/structural_identifiability.jl | 139 +++++++++---------
 .../stability_computation.jl                  |   4 +-
 .../parameter_type_designation.jl             |   2 +-
 9 files changed, 83 insertions(+), 80 deletions(-)

diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
index af2c0bb47e..be0becb3a9 100644
--- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
+++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
@@ -4,7 +4,7 @@
 function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
         plot_var = nothing, record_from_solution = BK.record_sol_default, jac = true, u0 = [], kwargs...)
     if !isautonomous(rs)
-        error("Attempting to create a `BifurcationProblem` for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+        error("Attempting to create a `BifurcationProblem` for a non-autonomous system (e.g. where some rate depend on $(get_iv(rs))). This is not possible.")
     end
 
     # Converts symbols to symbolics.
diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
index 91e868c25e..cc7ab65795 100644
--- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
+++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
@@ -38,7 +38,7 @@ Notes:
 function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative = true,
         neg_thres = -1e-20, u0 = [], kwargs...)
     if !isautonomous(rs)
-        error("Attempting to compute steady state for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+        error("Attempting to compute steady state for a non-autonomous system (e.g. where some rate depend on $(get_iv(rs))). This is not possible.")
     end
     ss_poly = steady_state_polynomial(rs, ps, u0)
     sols = HC.real_solutions(HC.solve(ss_poly; kwargs...))
diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index 99fceb6b58..c98254d442 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -520,7 +520,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
     if !isautonomous(rs)
-        error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(rs.iv)) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
+        error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(get_iv(rs))) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
     end
 
     # Generates system equations.
diff --git a/src/reactionsystem_serialisation/serialise_reactionsystem.jl b/src/reactionsystem_serialisation/serialise_reactionsystem.jl
index 7c647ef2a2..e7ecd69eee 100644
--- a/src/reactionsystem_serialisation/serialise_reactionsystem.jl
+++ b/src/reactionsystem_serialisation/serialise_reactionsystem.jl
@@ -147,13 +147,13 @@ function make_reaction_system_call(
     has_connection_type && (@string_append! reaction_system_string ", connection_type")
 
     # Potentially appends a combinatoric_ratelaws statement.
-    if !Symbolics.unwrap(rs.combinatoric_ratelaws)
+    if !Symbolics.unwrap(combinatoric_ratelaws(rs))
         @string_append! reaction_system_string ", combinatoric_ratelaws = false"
     end
 
     # Potentially appends `ReactionSystem` metadata value(s). Weird composite types are not supported.
-    if !isnothing(rs.metadata)
-        @string_append! reaction_system_string ", metadata = $(x_2_string(rs.metadata))"
+    if !isnothing(MT.get_metadata(rs))
+        @string_append! reaction_system_string ", metadata = $(x_2_string(MT.get_metadata(rs)))"
     end
 
     # Finalises the call. Appends potential annotation. If the system is complete, add a call for this. 
diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl
index 8103d4c61f..42ffec2d0e 100644
--- a/src/steady_state_stability.jl
+++ b/src/steady_state_stability.jl
@@ -48,7 +48,7 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps;
         tol = 10 * sqrt(eps(ss_val_type(u))), ss_jac = steady_state_jac(rs; u0 = u))
     # Warning checks.
     if !isautonomous(rs)
-        error("Attempting to compute stability for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+        error("Attempting to compute stability for a non-autonomous system (e.g. where some rate depend on $(get_iv(rs))). This is not possible.")
     end
 
     # If `u` is a vector of values, we convert it to a map. Also, if there are conservation laws,
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index 93ce90b541..d12cbd4a45 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -574,8 +574,8 @@ let
         end
     end
     V,W = getfield.(observed(rn), :lhs)
-    @test isequal(arguments(ModelingToolkit.unwrap(V)), Any[rn.iv, rn.sivs[1], rn.sivs[2]])
-    @test isequal(arguments(ModelingToolkit.unwrap(W)), Any[rn.iv, rn.sivs[2]])
+    @test isequal(arguments(ModelingToolkit.unwrap(V)), Any[Catalyst.get_iv(rs), Catalyst.get_sivs(rn)[1], Catalyst.get_sivs(rn)[2]])
+    @test isequal(arguments(ModelingToolkit.unwrap(W)), Any[Catalyst.get_iv(rs), Catalyst.get_sivs(rn)[2]])
 end
 
 # Checks that metadata is written properly.
diff --git a/test/extensions/structural_identifiability.jl b/test/extensions/structural_identifiability.jl
index 8ef35c293b..f3617c3409 100644
--- a/test/extensions/structural_identifiability.jl
+++ b/test/extensions/structural_identifiability.jl
@@ -1,7 +1,10 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, StructuralIdentifiability, Test
+using Catalyst, Logging, StructuralIdentifiability, Test
+
+# Sets the `loglevel`.
+loglevel = Logging.Error
 
 # Helper function for checking that results are correct identifiability calls from different packages.
 # Converts the output dicts from StructuralIdentifiability functions from "weird symbol => stuff" to "symbol => stuff" (the output have some strange meta data which prevents equality checks, this enables this).
@@ -28,12 +31,12 @@ let
         (pₑ*M,dₑ), 0 <--> E
         (pₚ*E,dₚ), 0 <--> P
     end
-    gi_1 = assess_identifiability(goodwind_oscillator_catalyst; measured_quantities=[:M])
-    li_1 = assess_local_identifiability(goodwind_oscillator_catalyst; measured_quantities=[:M])
-    ifs_1 = find_identifiable_functions(goodwind_oscillator_catalyst; measured_quantities=[:M])
+    gi_1 = assess_identifiability(goodwind_oscillator_catalyst; measured_quantities = [:M], loglevel)
+    li_1 = assess_local_identifiability(goodwind_oscillator_catalyst; measured_quantities = [:M], loglevel)
+    ifs_1 = find_identifiable_functions(goodwind_oscillator_catalyst; measured_quantities = [:M], loglevel)
 
     # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model.
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M])
     gi_2 = assess_identifiability(si_catalyst_ode)
     li_2 = assess_local_identifiability(si_catalyst_ode)
     ifs_2 = find_identifiable_functions(si_catalyst_ode)
@@ -45,9 +48,9 @@ let
         P'(t) = -dₚ*P(t) + pₚ*E(t),
         y1(t) = M(t)
     )
-    gi_3 = assess_identifiability(goodwind_oscillator_si)
-    li_3 = assess_local_identifiability(goodwind_oscillator_si)
-    ifs_3 = find_identifiable_functions(goodwind_oscillator_si)
+    gi_3 = assess_identifiability(goodwind_oscillator_si; loglevel)
+    li_3 = assess_local_identifiability(goodwind_oscillator_si; loglevel)
+    ifs_3 = find_identifiable_functions(goodwind_oscillator_si; loglevel)
 
     # Check outputs.
     @test sym_dict(gi_1) == sym_dict(gi_2) == sym_dict(gi_3)
@@ -74,15 +77,15 @@ let
         d, X4 --> 0
     end
     @unpack X2, X3 = rs_catalyst    
-    gi_1 = assess_identifiability(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f])
-    li_1 = assess_local_identifiability(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f])
-    ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f])
+    gi_1 = assess_identifiability(rs_catalyst; measured_quantities = [X2, X3], known_p = [:k2f], loglevel)
+    li_1 = assess_local_identifiability(rs_catalyst; measured_quantities = [X2, X3], known_p = [:k2f], loglevel)
+    ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities = [X2, X3], known_p = [:k2f], loglevel)
     
     # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model.
-    rs_ode = make_si_ode(rs_catalyst; measured_quantities=[X2, X3], known_p=[:k2f])
-    gi_2 = assess_identifiability(rs_ode)
-    li_2 = assess_local_identifiability(rs_ode)
-    ifs_2 = find_identifiable_functions(rs_ode)
+    rs_ode = make_si_ode(rs_catalyst; measured_quantities = [X2, X3], known_p = [:k2f])
+    gi_2 = assess_identifiability(rs_ode; loglevel)
+    li_2 = assess_local_identifiability(rs_ode; loglevel)
+    ifs_2 = find_identifiable_functions(rs_ode; loglevel)
 
     # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.).
     rs_si = @ODEmodel(
@@ -94,9 +97,9 @@ let
         y2(t) = X3,
         y3(t) = k2f
     )
-    gi_3 = assess_identifiability(rs_si)
-    li_3 = assess_local_identifiability(rs_si)
-    ifs_3 = find_identifiable_functions(rs_si)
+    gi_3 = assess_identifiability(rs_si; loglevel)
+    li_3 = assess_local_identifiability(rs_si; loglevel)
+    ifs_3 = find_identifiable_functions(rs_si; loglevel)
 
     # Check outputs.
     @test sym_dict(gi_1) == sym_dict(gi_2) == sym_dict(gi_3)
@@ -126,15 +129,15 @@ let
         (kA*X3, kD), Yi <--> Ya
     end
     @unpack X1, X2, X3, X4, k1, k2, Yi, Ya, k1, kD = rs_catalyst
-    gi_1 = assess_identifiability(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD])
-    li_1 = assess_local_identifiability(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD])
-    ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD])
+    gi_1 = assess_identifiability(rs_catalyst; measured_quantities = [X1 + Yi, Ya], known_p = [k1, kD], loglevel)
+    li_1 = assess_local_identifiability(rs_catalyst; measured_quantities = [X1 + Yi, Ya], known_p = [k1, kD], loglevel)
+    ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities = [X1 + Yi, Ya], known_p = [k1, kD], loglevel)
 
     # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model.
     rs_ode = make_si_ode(rs_catalyst; measured_quantities=[X1 + Yi, Ya], known_p=[k1, kD], remove_conserved=false)
-    gi_2 = assess_identifiability(rs_ode)
-    li_2 = assess_local_identifiability(rs_ode)
-    ifs_2 = find_identifiable_functions(rs_ode)
+    gi_2 = assess_identifiability(rs_ode; loglevel)
+    li_2 = assess_local_identifiability(rs_ode; loglevel)
+    ifs_2 = find_identifiable_functions(rs_ode; loglevel)
 
     # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.).
     rs_si = @ODEmodel(
@@ -150,9 +153,9 @@ let
         y3(t) = k1,
         y4(t) = kD
     )
-    gi_3 = assess_identifiability(rs_si)
-    li_3 = assess_local_identifiability(rs_si)
-    ifs_3 = find_identifiable_functions(rs_si)
+    gi_3 = assess_identifiability(rs_si; loglevel)
+    li_3 = assess_local_identifiability(rs_si; loglevel)
+    ifs_3 = find_identifiable_functions(rs_si; loglevel)
     
     # Check outputs.
     @test sym_dict(gi_1) == sym_dict(gi_2) == sym_dict(gi_3)
@@ -174,31 +177,31 @@ let
         (pₚ*E,dₚ), 0 <--> P
     end
     @unpack M, E, P, pₑ, pₚ, pₘ = goodwind_oscillator_catalyst
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p=[:pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M], known_p=[:pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M, :E], known_p=[:pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M], known_p=[:pₑ, :pₚ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[:M, :E], known_p=[:pₑ, :pₚ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p=[pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M], known_p=[pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M, E], known_p=[pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M], known_p=[pₑ, pₚ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M, E], known_p=[pₑ, pₚ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M + pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[M + E, pₑ*M], known_p=[:pₑ])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities=[pₑ, pₚ], known_p=[pₑ])    
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p = [:pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M], known_p = [:pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M, :E], known_p = [:pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M], known_p = [:pₑ, :pₚ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M, :E], known_p = [:pₑ, :pₚ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p = [pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M], known_p = [pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M, E], known_p = [pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M], known_p = [pₑ, pₚ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M, E], known_p = [pₑ, pₚ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M + pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M + E, pₑ*M], known_p = [:pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [pₑ, pₚ], known_p = [pₑ])    
 
     # Tests using model.component style (have to make system complete first).
     gw_osc_complt = complete(goodwind_oscillator_catalyst)
-    @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M]) isa ODE
-    @test make_si_ode(gw_osc_complt; known_p=[gw_osc_complt.pₑ]) isa ODE
-    @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M], known_p=[gw_osc_complt.pₑ]) isa ODE
-    @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M, gw_osc_complt.E], known_p=[gw_osc_complt.pₑ]) isa ODE
-    @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M], known_p=[gw_osc_complt.pₑ, gw_osc_complt.pₚ]) isa ODE
-    @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M], known_p = [:pₚ]) isa ODE
-    @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M*gw_osc_complt.E]) isa ODE
+    @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M]) isa ODE
+    @test make_si_ode(gw_osc_complt; known_p = [gw_osc_complt.pₑ]) isa ODE
+    @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M], known_p = [gw_osc_complt.pₑ]) isa ODE
+    @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M, gw_osc_complt.E], known_p = [gw_osc_complt.pₑ]) isa ODE
+    @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M], known_p = [gw_osc_complt.pₑ, gw_osc_complt.pₚ]) isa ODE
+    @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M], known_p = [:pₚ]) isa ODE
+    @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M*gw_osc_complt.E]) isa ODE
 end
 
 # Tests for hierarchical model with conservation laws at both top and internal levels.
@@ -213,15 +216,15 @@ let
     @named rs_catalyst = compose(rs1, [rs2])
     rs_catalyst = complete(rs_catalyst)
     @unpack X1, X2, k1, k2 = rs1
-    gi_1 = assess_identifiability(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1])
-    li_1 = assess_local_identifiability(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1])
-    ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1])
+    gi_1 = assess_identifiability(rs_catalyst; measured_quantities = [X1, X2, rs2.X3], known_p = [k1], loglevel)
+    li_1 = assess_local_identifiability(rs_catalyst; measured_quantities = [X1, X2, rs2.X3], known_p = [k1], loglevel)
+    ifs_1 = find_identifiable_functions(rs_catalyst; measured_quantities = [X1, X2, rs2.X3], known_p = [k1], loglevel)
 
     # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model.
-    rs_ode = make_si_ode(rs_catalyst; measured_quantities=[X1, X2, rs2.X3], known_p=[k1])
-    gi_2 = assess_identifiability(rs_ode)
-    li_2 = assess_local_identifiability(rs_ode)
-    ifs_2 = find_identifiable_functions(rs_ode)
+    rs_ode = make_si_ode(rs_catalyst; measured_quantities = [X1, X2, rs2.X3], known_p = [k1])
+    gi_2 = assess_identifiability(rs_ode; loglevel)
+    li_2 = assess_local_identifiability(rs_ode; loglevel)
+    ifs_2 = find_identifiable_functions(rs_ode; loglevel)
 
     # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.).
     rs_si = @ODEmodel(
@@ -234,9 +237,9 @@ let
         y3(t) = rs2₊X3,
         y4(t) = k1
     )
-    gi_3 = assess_identifiability(rs_si)
-    li_3 = assess_local_identifiability(rs_si)
-    ifs_3 = find_identifiable_functions(rs_si)
+    gi_3 = assess_identifiability(rs_si; loglevel)
+    li_3 = assess_local_identifiability(rs_si; loglevel)
+    ifs_3 = find_identifiable_functions(rs_si; loglevel)
         
     # Check outputs.
     @test sym_dict(gi_1) == sym_dict(gi_3)
@@ -326,9 +329,9 @@ let
     @unpack p, d, k, c1, c2 = rs
     
     # Tests identifiability assessment when all unknowns are measured.
-    gi_1 = assess_identifiability(rs; measured_quantities=[:X, :V, :C])
-    li_1 = assess_local_identifiability(rs; measured_quantities=[:X, :V, :C])
-    ifs_1 = find_identifiable_functions(rs; measured_quantities=[:X, :V, :C])
+    gi_1 = assess_identifiability(rs; measured_quantities = [:X, :V, :C], loglevel)
+    li_1 = assess_local_identifiability(rs; measured_quantities = [:X, :V, :C], loglevel)
+    ifs_1 = find_identifiable_functions(rs; measured_quantities = [:X, :V, :C], loglevel)
     @test sym_dict(gi_1) == Dict([:X => :globally, :C => :globally, :V => :globally, :k => :globally, 
                       :c1 => :nonidentifiable, :c2 => :nonidentifiable, :p => :globally, :d => :globally])
     @test sym_dict(li_1) == Dict([:X => 1, :C => 1, :V => 1, :k => 1, :c1 => 0, :c2 => 0, :p => 1, :d => 1])
@@ -336,9 +339,9 @@ let
     
     # Tests identifiability assessment when only variables are measured. 
     # Checks that a parameter in an equation can be set as known.
-    gi_2 = assess_identifiability(rs; measured_quantities=[:V, :C], known_p = [:c1])
-    li_2 = assess_local_identifiability(rs; measured_quantities=[:V, :C], known_p = [:c1])
-    ifs_2 = find_identifiable_functions(rs; measured_quantities=[:V, :C], known_p = [:c1])
+    gi_2 = assess_identifiability(rs; measured_quantities = [:V, :C], known_p = [:c1], loglevel)
+    li_2 = assess_local_identifiability(rs; measured_quantities = [:V, :C], known_p = [:c1], loglevel)
+    ifs_2 = find_identifiable_functions(rs; measured_quantities = [:V, :C], known_p = [:c1], loglevel)
     @test sym_dict(gi_2) == Dict([:X => :nonidentifiable, :C => :globally, :V => :globally, :k => :nonidentifiable, 
                       :c1 => :globally, :c2 => :nonidentifiable, :p => :nonidentifiable, :d => :globally])
     @test sym_dict(li_2) == Dict([:X => 0, :C => 1, :V => 1, :k => 0, :c1 => 1, :c2 => 0, :p => 0, :d => 1])
@@ -354,7 +357,7 @@ let
     measured_quantities = [:X]
     
     # Computes bifurcation diagram.
-    @test_throws Exception assess_identifiability(incomplete_network; measured_quantities)
-    @test_throws Exception assess_local_identifiability(incomplete_network; measured_quantities)
-    @test_throws Exception find_identifiable_functions(incomplete_network; measured_quantities)
+    @test_throws Exception assess_identifiability(incomplete_network; measured_quantities, loglevel)
+    @test_throws Exception assess_local_identifiability(incomplete_network; measured_quantities, loglevels)
+    @test_throws Exception find_identifiable_functions(incomplete_network; measured_quantities, loglevel)
 end
\ No newline at end of file
diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl
index 36f6125f02..c786ec68cf 100644
--- a/test/miscellaneous_tests/stability_computation.jl
+++ b/test/miscellaneous_tests/stability_computation.jl
@@ -29,7 +29,7 @@ let
               :d => 0.5 + rand(rng))
 
         # Computes stability using various jacobian options.
-        sss = hc_steady_states(rn, ps)
+        sss = hc_steady_states(rn, ps; show_progress=false)
         stabs_1 = [steady_state_stability(ss, rn, ps) for ss in sss]
         stabs_2 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss]
 
@@ -71,7 +71,7 @@ let
     ps_3 = [rn.k1 => 8.0, rn.k2 => 2.0, rn.k3 => 1.0, rn.k4 => 1.5, rn.kD1 => 0.5, rn.kD2 => 4.0]
     
     # Computes stability using various input forms, and checks that the output is correct. 
-    sss = hc_steady_states(rn, ps_1; u0 = u0_1)
+    sss = hc_steady_states(rn, ps_1; u0 = u0_1, show_progress=false)
     for u0 in [u0_1, u0_2, u0_3, u0_4], ps in [ps_1, ps_2, ps_3]
         stab_1 =  [steady_state_stability(ss, rn, ps) for ss in sss]
         ss_jac = steady_state_jac(rn; u0 = u0)
diff --git a/test/reactionsystem_core/parameter_type_designation.jl b/test/reactionsystem_core/parameter_type_designation.jl
index d6f7a7d13d..ca208a926e 100644
--- a/test/reactionsystem_core/parameter_type_designation.jl
+++ b/test/reactionsystem_core/parameter_type_designation.jl
@@ -64,7 +64,7 @@ let
     for p in p_alts
         oprob = ODEProblem(rs, u0, (0.0, 1000.0), p; abstol = 1e-10, reltol = 1e-10)
         sol = solve(oprob, Tsit5())
-        @test all(sol[end] .≈ 1.0)
+        @test all(sol.u[end] .≈ 1.0)
     end
 end
 

From 300d5e09961b8c7be06730208325f77463a0689b Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 8 Jun 2024 14:02:10 -0400
Subject: [PATCH 2/4] up

---
 .../homotopy_continuation_extension.jl                    | 2 +-
 src/reactionsystem.jl                                     | 2 +-
 src/reactionsystem_serialisation/serialisation_support.jl | 2 +-
 test/dsl/dsl_options.jl                                   | 4 ++--
 test/reactionsystem_core/coupled_equation_crn_systems.jl  | 8 ++++----
 test/reactionsystem_core/reactionsystem.jl                | 2 +-
 6 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
index 6b70922e45..db7a8ff0c4 100644
--- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
+++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
@@ -67,7 +67,7 @@ function make_int_exps(expr)
     wrap(Rewriters.Postwalk(Rewriters.PassThrough(___make_int_exps))(unwrap(expr))).val
 end
 function ___make_int_exps(expr)
-    !istree(expr) && return expr
+    !iscall(expr) && return expr
     if (operation(expr) == ^)
         if isinteger(arguments(expr)[2])
             return arguments(expr)[1]^Int64(arguments(expr)[2])
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 8373bf4990..bf20c94113 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -1479,4 +1479,4 @@ function validate(rs::ReactionSystem, info::String = "")
 end
 
 # Checks if a unit consist of exponents with base 1 (and is this unitless).
-unitless_exp(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)
+unitless_exp(u) = iscall(u) && (operation(u) == ^) && (arguments(u)[1] == 1)
diff --git a/src/reactionsystem_serialisation/serialisation_support.jl b/src/reactionsystem_serialisation/serialisation_support.jl
index 6d8b2152df..18a13c0d80 100644
--- a/src/reactionsystem_serialisation/serialisation_support.jl
+++ b/src/reactionsystem_serialisation/serialisation_support.jl
@@ -240,7 +240,7 @@ const SKIPPED_METADATA = [ModelingToolkit.MTKVariableTypeCtx, Symbolics.Variable
 # Potentially strips the call for a symbolics. E.g. X(t) becomes X (but p remains p). This is used 
 # when variables are written to files, as in code they are used without the call part.
 function strip_call(sym)
-    return istree(sym) ? Sym{Real}(Symbolics.getname(sym)) : sym
+    return iscall(sym) ? Sym{Real}(Symbolics.getname(sym)) : sym
 end
 
 # For an vector of symbolics, creates a dictionary taking each symbolics to each call-stripped form.
diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl
index d12cbd4a45..85dd068aec 100644
--- a/test/dsl/dsl_options.jl
+++ b/test/dsl/dsl_options.jl
@@ -574,8 +574,8 @@ let
         end
     end
     V,W = getfield.(observed(rn), :lhs)
-    @test isequal(arguments(ModelingToolkit.unwrap(V)), Any[Catalyst.get_iv(rs), Catalyst.get_sivs(rn)[1], Catalyst.get_sivs(rn)[2]])
-    @test isequal(arguments(ModelingToolkit.unwrap(W)), Any[Catalyst.get_iv(rs), Catalyst.get_sivs(rn)[2]])
+    @test isequal(arguments(ModelingToolkit.unwrap(V)), Any[Catalyst.get_iv(rn), Catalyst.get_sivs(rn)[1], Catalyst.get_sivs(rn)[2]])
+    @test isequal(arguments(ModelingToolkit.unwrap(W)), Any[Catalyst.get_iv(rn), Catalyst.get_sivs(rn)[2]])
 end
 
 # Checks that metadata is written properly.
diff --git a/test/reactionsystem_core/coupled_equation_crn_systems.jl b/test/reactionsystem_core/coupled_equation_crn_systems.jl
index 33db205dd9..27a77af97d 100644
--- a/test/reactionsystem_core/coupled_equation_crn_systems.jl
+++ b/test/reactionsystem_core/coupled_equation_crn_systems.jl
@@ -47,7 +47,7 @@ let
     # Checks that the correct steady state is found through ODEProblem.
     oprob = ODEProblem(coupled_rs, u0, tspan, ps)
     osol = solve(oprob, Vern7(); abstol = 1e-8, reltol = 1e-8)
-    @test osol[end] ≈ [2.0, 1.0]
+    @test osol.u[end] ≈ [2.0, 1.0]
 
     # Checks that the correct steady state is found through NonlinearProblem.
     nlprob = NonlinearProblem(coupled_rs, u0, ps)
@@ -116,7 +116,7 @@ let
     for coupled_rs in [coupled_rs_prog, coupled_rs_extended, coupled_rs_dsl]
         oprob = ODEProblem(coupled_rs, u0, tspan, ps)
         osol = solve(oprob, Vern7(); abstol = 1e-8, reltol = 1e-8)
-        osol[end] ≈ [10.0, 10.0, 11.0, 11.0]
+        osol.u[end] ≈ [10.0, 10.0, 11.0, 11.0]
     end
 end
 
@@ -160,7 +160,7 @@ let
     # Checks that the correct steady state is found through ODEProblem.
     oprob = ODEProblem(coupled_rs, u0, tspan, ps; structural_simplify = true)
     osol = solve(oprob, Rosenbrock23(); abstol = 1e-8, reltol = 1e-8)
-    @test osol[end] ≈ [2.0, 3.0]
+    @test osol.u[end] ≈ [2.0, 3.0]
 
     # Checks that the correct steady state is found through NonlinearProblem.
     nlprob = NonlinearProblem(coupled_rs, u0, ps)
@@ -207,7 +207,7 @@ let
     osol = solve(oprob, Rosenbrock23(); abstol = 1e-8, reltol = 1e-8)
     sssol = solve(ssprob, DynamicSS(Rosenbrock23()); abstol = 1e-8, reltol = 1e-8)
     nlsol = solve(nlprob; abstol = 1e-8, reltol = 1e-8)
-    @test osol[end] ≈ sssol ≈ nlsol
+    @test osol.u[end] ≈ sssol ≈ nlsol
 end
 
 
diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl
index 6876026e0b..65e42180cc 100644
--- a/test/reactionsystem_core/reactionsystem.jl
+++ b/test/reactionsystem_core/reactionsystem.jl
@@ -356,7 +356,7 @@ let
     oprob1 = ODEProblem(osys, u0map, tspan, pmap)
     sts = [B, D, E, C]
     syms = [:B, :D, :E, :C]
-    ofun = ODEFunction(f!; syms)
+    ofun = ODEFunction(f!; sys = ModelingToolkit.SymbolCache(syms))
     oprob2 = ODEProblem(ofun, u0, tspan, p)
     saveat = tspan[2] / 50
     abstol = 1e-10

From 1f8d805fd01e2cccf1fdc168d72f98cfa431b229 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 8 Jun 2024 16:17:11 -0400
Subject: [PATCH 3/4] fix unnecessary logging

---
 Project.toml                                      |  3 ++-
 test/extensions/structural_identifiability.jl     | 14 +++++++-------
 test/miscellaneous_tests/stability_computation.jl |  4 ++--
 3 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/Project.toml b/Project.toml
index 225d305c09..b9371b3656 100644
--- a/Project.toml
+++ b/Project.toml
@@ -67,6 +67,7 @@ DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -83,4 +84,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["BifurcationKit", "DiffEqCallbacks", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]
+test = ["BifurcationKit", "DiffEqCallbacks", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "Logging", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]
diff --git a/test/extensions/structural_identifiability.jl b/test/extensions/structural_identifiability.jl
index f3617c3409..6b411cc9db 100644
--- a/test/extensions/structural_identifiability.jl
+++ b/test/extensions/structural_identifiability.jl
@@ -37,9 +37,9 @@ let
 
     # Identifiability analysis for Catalyst converted to StructuralIdentifiability.jl model.
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M])
-    gi_2 = assess_identifiability(si_catalyst_ode)
-    li_2 = assess_local_identifiability(si_catalyst_ode)
-    ifs_2 = find_identifiable_functions(si_catalyst_ode)
+    gi_2 = assess_identifiability(si_catalyst_ode; loglevel)
+    li_2 = assess_local_identifiability(si_catalyst_ode; loglevel)
+    ifs_2 = find_identifiable_functions(si_catalyst_ode; loglevel)
 
     # Identifiability analysis for StructuralIdentifiability.jl model (declare this overwrites e.g. X2 variable etc.).
     goodwind_oscillator_si = @ODEmodel(
@@ -264,14 +264,14 @@ let
         k1, x1 --> x2
     end
     # Measure the source
-    id_report = assess_identifiability(rs, measured_quantities = [:x1])
+    id_report = assess_identifiability(rs; measured_quantities = [:x1], loglevel)
     @test sym_dict(id_report) == Dict(
         :x1 => :globally,
         :x2 => :nonidentifiable,
         :k1 => :globally
     )
     # Measure the target instead
-    id_report = assess_identifiability(rs, measured_quantities = [:x2])
+    id_report = assess_identifiability(rs; measured_quantities = [:x2], loglevel)
     @test sym_dict(id_report) == Dict(
         :x1 => :globally,
         :x2 => :globally,
@@ -303,13 +303,13 @@ let
         1, x1 --> x2
         1, x2 --> x3
     end
-    id_report = assess_identifiability(rs, measured_quantities = [:x3])
+    id_report = assess_identifiability(rs; measured_quantities = [:x3], loglevel)
     @test sym_dict(id_report) == Dict(
         :x1 => :globally,
         :x2 => :globally,
         :x3 => :globally,
     )
-    @test length(find_identifiable_functions(rs, measured_quantities = [:x3])) == 1
+    @test length(find_identifiable_functions(rs; measured_quantities = [:x3], loglevel)) == 1
 end
 
 
diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl
index c786ec68cf..542fcd4e51 100644
--- a/test/miscellaneous_tests/stability_computation.jl
+++ b/test/miscellaneous_tests/stability_computation.jl
@@ -29,7 +29,7 @@ let
               :d => 0.5 + rand(rng))
 
         # Computes stability using various jacobian options.
-        sss = hc_steady_states(rn, ps; show_progress=false)
+        sss = hc_steady_states(rn, ps; show_progress = false)
         stabs_1 = [steady_state_stability(ss, rn, ps) for ss in sss]
         stabs_2 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss]
 
@@ -71,7 +71,7 @@ let
     ps_3 = [rn.k1 => 8.0, rn.k2 => 2.0, rn.k3 => 1.0, rn.k4 => 1.5, rn.kD1 => 0.5, rn.kD2 => 4.0]
     
     # Computes stability using various input forms, and checks that the output is correct. 
-    sss = hc_steady_states(rn, ps_1; u0 = u0_1, show_progress=false)
+    sss = hc_steady_states(rn, ps_1; u0 = u0_1, show_progress = false)
     for u0 in [u0_1, u0_2, u0_3, u0_4], ps in [ps_1, ps_2, ps_3]
         stab_1 =  [steady_state_stability(ss, rn, ps) for ss in sss]
         ss_jac = steady_state_jac(rn; u0 = u0)

From 0f6cab61efdb029e70c128fcb0e446ff2e48f123 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 8 Jun 2024 17:41:35 -0400
Subject: [PATCH 4/4] SI updates

---
 ext/CatalystStructuralIdentifiabilityExtension.jl      |  1 +
 .../structural_identifiability_extension.jl            |  8 ++++----
 test/extensions/structural_identifiability.jl          | 10 +++++-----
 3 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/ext/CatalystStructuralIdentifiabilityExtension.jl b/ext/CatalystStructuralIdentifiabilityExtension.jl
index 0a442affa8..ed80952bd1 100644
--- a/ext/CatalystStructuralIdentifiabilityExtension.jl
+++ b/ext/CatalystStructuralIdentifiabilityExtension.jl
@@ -2,6 +2,7 @@ module CatalystStructuralIdentifiabilityExtension
 
 # Fetch packages.
 using Catalyst
+import DataStructures.OrderedDict
 import StructuralIdentifiability as SI
 
 # Creates and exports hc_steady_states function.
diff --git a/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl b/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl
index d4cc1ead23..c5906eaf6e 100644
--- a/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl
+++ b/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl
@@ -189,8 +189,8 @@ function make_measured_quantities(
         rs::ReactionSystem, measured_quantities::Vector{T}, known_p::Vector{S},
         conseqs; ignore_no_measured_warn = false) where {T, S}
     # Warning if the user didn't give any measured quantities.
-    if ignore_no_measured_warn || isempty(measured_quantities)
-        @warn "No measured quantity provided to the `measured_quantities` argument, any further identifiability analysis will likely fail. You can disable this warning by setting `ignore_no_measured_warn=true`."
+    if !ignore_no_measured_warn && isempty(measured_quantities)
+        @warn "No measured quantity provided to the `measured_quantities` argument, any further identifiability analysis will likely fail. You can disable this warning by setting `ignore_no_measured_warn = true`."
     end
 
     # Appends the known parameters to the measured_quantities vector. Converts any Symbols to symbolics.
@@ -220,9 +220,9 @@ end
 # Sorts the output according to their input order (defaults to the `[unknowns; parameters]` order).
 function make_output(out, funcs_to_check, conseqs)
     funcs_to_check = vector_subs(funcs_to_check, conseqs)
-    out = Dict(zip(vector_subs(keys(out), conseqs), values(out)))
+    out = OrderedDict(zip(vector_subs(keys(out), conseqs), values(out)))
     sortdict = Dict(ftc => i for (i, ftc) in enumerate(funcs_to_check))
-    return sort(out; by = x -> sortdict[x])
+    return sort!(out; by = x -> sortdict[x])
 end
 
 # For a vector of expressions and a conservation law, substitutes the law into every equation.
diff --git a/test/extensions/structural_identifiability.jl b/test/extensions/structural_identifiability.jl
index 6b411cc9db..105d990d48 100644
--- a/test/extensions/structural_identifiability.jl
+++ b/test/extensions/structural_identifiability.jl
@@ -178,13 +178,13 @@ let
     end
     @unpack M, E, P, pₑ, pₚ, pₘ = goodwind_oscillator_catalyst
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p = [:pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p = [:pₑ], ignore_no_measured_warn = true)
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M], known_p = [:pₑ])
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M, :E], known_p = [:pₑ])
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M], known_p = [:pₑ, :pₚ])
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [:M, :E], known_p = [:pₑ, :pₚ])
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M])
-    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p = [pₑ])
+    si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; known_p = [pₑ], ignore_no_measured_warn = true)
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M], known_p = [pₑ])
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M, E], known_p = [pₑ])
     si_catalyst_ode = make_si_ode(goodwind_oscillator_catalyst; measured_quantities = [M], known_p = [pₑ, pₚ])
@@ -196,7 +196,7 @@ let
     # Tests using model.component style (have to make system complete first).
     gw_osc_complt = complete(goodwind_oscillator_catalyst)
     @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M]) isa ODE
-    @test make_si_ode(gw_osc_complt; known_p = [gw_osc_complt.pₑ]) isa ODE
+    @test make_si_ode(gw_osc_complt; known_p = [gw_osc_complt.pₑ], ignore_no_measured_warn = true) isa ODE
     @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M], known_p = [gw_osc_complt.pₑ]) isa ODE
     @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M, gw_osc_complt.E], known_p = [gw_osc_complt.pₑ]) isa ODE
     @test make_si_ode(gw_osc_complt; measured_quantities = [gw_osc_complt.M], known_p = [gw_osc_complt.pₑ, gw_osc_complt.pₚ]) isa ODE
@@ -288,7 +288,7 @@ let
         b, A0 --> 2A2
         c, A0 --> A1 + A2
     end
-    id_report = assess_identifiability(rs, measured_quantities = [:A0, :A1, :A2])
+    id_report = assess_identifiability(rs; measured_quantities = [:A0, :A1, :A2], loglevel)
     @test sym_dict(id_report) == Dict(
         :A0 => :globally,
         :A1 => :globally,
@@ -358,6 +358,6 @@ let
     
     # Computes bifurcation diagram.
     @test_throws Exception assess_identifiability(incomplete_network; measured_quantities, loglevel)
-    @test_throws Exception assess_local_identifiability(incomplete_network; measured_quantities, loglevels)
+    @test_throws Exception assess_local_identifiability(incomplete_network; measured_quantities, loglevel)
     @test_throws Exception find_identifiable_functions(incomplete_network; measured_quantities, loglevel)
 end
\ No newline at end of file