Merge 0f6cab6 into ea2175c

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"]

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.

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...))
@@ -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])

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.

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.

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)

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.

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.

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. 

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,

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(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.