Merge pull request #887 from SciML/symbolics_fixes

TorkelE · web-flow · commit 4fcedea2c0bf · 2024-06-08T14:01:53.000-04:00
[WIP] Fix Symbolics related stuff
diff --git a/Project.toml b/Project.toml
@@ -22,7 +22,6 @@ Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
@@ -38,8 +37,8 @@ CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability"
 
 [compat]
 BifurcationKit = "0.3"
-DataStructures = "0.18"
 Combinatorics = "1.0.2"
+DataStructures = "0.18"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
 DynamicPolynomials = "0.5"
@@ -50,15 +49,14 @@ JumpProcesses = "9.3.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "9.11.0"
+ModelingToolkit = "9.16.0"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
 Setfield = "1"
-StructuralIdentifiability = "0.5.1"
-SymbolicUtils = "1.0.3"
-Symbolics = "5.27"
+StructuralIdentifiability = "0.5.8"
+Symbolics = "5.30.1"
 Unitful = "1.12.4"
 julia = "1.10"
 
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -51,7 +51,7 @@ IncompleteLU = "0.2"
 JumpProcesses = "9.11"
 Latexify = "0.16"
 LinearSolve = "2.30"
-ModelingToolkit = "9.15"
+ModelingToolkit = "9.16.0"
 NonlinearSolve = "3.12"
 Optim = "1.9"
 Optimization = "3.25"
@@ -68,5 +68,5 @@ SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 SteadyStateDiffEq = "2.2"
 StochasticDiffEq = "6.65"
-StructuralIdentifiability = "0.5.7"
-Symbolics = "5.28"
+StructuralIdentifiability = "0.5.8"
+Symbolics = "5.30.1"
diff --git a/ext/CatalystHomotopyContinuationExtension.jl b/ext/CatalystHomotopyContinuationExtension.jl
@@ -6,7 +6,8 @@ import DynamicPolynomials
 import ModelingToolkit as MT
 import HomotopyContinuation as HC
 import Setfield: @set
-import Symbolics: unwrap, wrap, Rewriters, symtype, issym, istree
+import Symbolics: unwrap, wrap, Rewriters, symtype, issym
+using Symbolics: iscall
 
 # Creates and exports hc_steady_states function.
 include("CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl")
diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
@@ -80,7 +80,7 @@ end
 # If the input is a fraction, removes the denominator.
 function remove_denominators(expr)
     s_expr = simplify_fractions(expr)
-    !istree(expr) && return expr
+    !iscall(expr) && return expr
     if operation(s_expr) == /
         return remove_denominators(arguments(s_expr)[1])
     end
diff --git a/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl b/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl
@@ -73,7 +73,7 @@ function SI.assess_local_identifiability(rs::ReactionSystem, args...;
     funcs_to_check = make_ftc(funcs_to_check, conseqs, vars)
 
     # Computes identifiability and converts it to a easy to read form.
-    out = SI.assess_local_identifiability(sys, args...; measured_quantities,
+    out = SI.assess_local_identifiability(osys, args...; measured_quantities,
         funcs_to_check, kwargs...)
     return make_output(out, funcs_to_check, reverse.(conseqs))
 end
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
@@ -23,8 +23,8 @@ using RuntimeGeneratedFunctions
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
 import Symbolics: BasicSymbolic
-import SymbolicUtils
-using ModelingToolkit: Symbolic, value, istree, get_unknowns, get_ps, get_iv, get_systems,
+using Symbolics: iscall
+using ModelingToolkit: Symbolic, value, get_unknowns, get_ps, get_iv, get_systems,
                        get_eqs, get_defaults, toparam, get_var_to_name, get_observed,
                        getvar
 
@@ -44,6 +44,7 @@ import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertice
 import DataStructures: OrderedDict, OrderedSet
 import Parameters: @with_kw_noshow
 import Symbolics: occursin, wrap
+import Symbolics.RewriteHelpers: hasnode, replacenode
 
 # globals for the modulate
 function default_time_deriv()
diff --git a/src/latexify_recipes.jl b/src/latexify_recipes.jl
@@ -193,7 +193,7 @@ function make_stoich_str(spec, stoich, subber; mathrm = true, kwargs...)
     if isequal(stoich, one(stoich))
         prestr * latexraw(subber(spec); kwargs...) * poststr
     else
-        if (stoich isa Symbolic) && istree(stoich)
+        if (stoich isa Symbolic) && iscall(stoich)
             LaTeXString("(") *
             latexraw(subber(stoich); kwargs...) *
             LaTeXString(")") *
diff --git a/src/reaction.jl b/src/reaction.jl
@@ -260,7 +260,7 @@ function print_rxside(io::IO, specs, stoich)
                      spec : MT.operation(spec)
             if isequal(stoich[i], one(stoich[i]))
                 print(io, prspec)
-            elseif istree(stoich[i])
+            elseif iscall(stoich[i])
                 print(io, "(", stoich[i], ")*", prspec)
             else
                 print(io, stoich[i], "*", prspec)
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
@@ -1464,7 +1464,7 @@ function validate(rs::ReactionSystem, info::String = "")
         # Needs additional checks because for cases: (1.0^n) and (1.0^n1)*(1.0^n2).
         # These are not considered (be default) considered equal to `1.0` for unitless reactions.
         isequal(rxunits, rateunits) && continue
-        if istree(rxunits)
+        if iscall(rxunits)
             unitless_exp(rxunits) && continue
             (operation(rxunits) == *) &&
                 all(unitless_exp(arg) for arg in arguments(rxunits)) && continue
diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
@@ -440,8 +440,8 @@ end
 
 # Finds and differentials in an expression, and sets these to 0.
 function remove_diffs(expr)
-    if Symbolics._occursin(Symbolics.is_derivative, expr)
-        return Symbolics.replace(expr, diff_2_zero)
+    if hasnode(Symbolics.is_derivative, expr)
+        return replacenode(expr, diff_2_zero)
     else
         return expr
     end
@@ -556,8 +556,8 @@ function nonlinear_convert_differentials_check(rs::ReactionSystem)
         # If the lhs upper level function is not a differential w.r.t. time.
         # If the contenct of the differential is not a variable (and nothing more).
         # If either of this is a case, throws the warning.
-        if Symbolics._occursin(Symbolics.is_derivative, eq.rhs) ||
-           !Symbolics.istree(eq.lhs) ||
+        if hasnode(Symbolics.is_derivative, eq.rhs) ||
+           !Symbolics.is_derivative(eq.lhs) ||
            !isequal(Symbolics.operation(eq.lhs), Differential(get_iv(rs))) ||
            (length(arguments(eq.lhs)) != 1) ||
            !any(isequal(arguments(eq.lhs)[1]), nonspecies(rs))
@@ -892,7 +892,7 @@ function to_multivariate_poly(polyeqs::AbstractVector{Symbolics.BasicSymbolic{Re
 
     pvar2sym, sym2term = SymbolicUtils.get_pvar2sym(), SymbolicUtils.get_sym2term()
     ps = map(polyeqs) do x
-        if istree(x) && operation(x) == (/)
+        if iscall(x) && operation(x) == (/)
             error("We should not be able to get here, please contact the package authors.")
         else
             PolyForm(x, pvar2sym, sym2term).p
diff --git a/src/reactionsystem_serialisation/serialisation_support.jl b/src/reactionsystem_serialisation/serialisation_support.jl
@@ -105,9 +105,9 @@ function sym_2_declaration_string(sym; multiline_format = false)
     dec_string = "$(sym)"
 
     # If the symbol has a non-default type, appends the declaration of this.
-    # Assumes that the type is on the form `SymbolicUtils.BasicSymbolic{X}`. Contain error checks
+    # Assumes that the type is on the form `BasicSymbolic{X}`. Contain error checks
     # to ensure that this is the case.
-    if !(sym isa SymbolicUtils.BasicSymbolic{Real})
+    if !(sym isa BasicSymbolic{Real})
         sym_type = String(Symbol(typeof(Symbolics.unwrap(sym))))
         if (get_substring(sym_type, 1, 28) != "SymbolicUtils.BasicSymbolic{") ||
            (get_char_end(sym_type, 0) != '}')
@@ -142,7 +142,7 @@ end
 # and metadata values (which hopefully almost exclusively) have simple, supported, types. Ideally,
 # more supported types can be added here.
 x_2_string(x::Num) = expression_2_string(x)
-x_2_string(x::SymbolicUtils.BasicSymbolic{<:Real}) = expression_2_string(x)
+x_2_string(x::BasicSymbolic{<:Real}) = expression_2_string(x)
 x_2_string(x::Bool) = string(x)
 x_2_string(x::String) = "\"$x\""
 x_2_string(x::Char) = "\'$x\'"
@@ -289,6 +289,6 @@ end
 function complicated_declaration(sym)
     isempty(get_metadata_to_declare(sym)) || (return true)
     ModelingToolkit.hasdefault(sym) && (return true)
-    (sym isa SymbolicUtils.BasicSymbolic{Real}) || (return true)
+    (sym isa BasicSymbolic{Real}) || (return true)
     return false
 end
diff --git a/src/reactionsystem_serialisation/serialise_fields.jl b/src/reactionsystem_serialisation/serialise_fields.jl
@@ -445,7 +445,7 @@ function discrete_event_string(discrete_event, strip_call_dict)
     # Creates the string corresponding to the conditions. The special check is if the condition is
     # an expression like `X > 5.0`. Here, "(...)" is added for purely aesthetic reasons.
     condition_string = x_2_string(discrete_event.condition)
-    if discrete_event.condition isa SymbolicUtils.BasicSymbolic
+    if discrete_event.condition isa BasicSymbolic
         @string_prepend! "(" condition_string
         @string_append! condition_string ")"
     end
diff --git a/src/registered_functions.jl b/src/registered_functions.jl
@@ -117,7 +117,7 @@ expand_registered_functions(expr)
 Takes an expression, and expands registered function expressions. E.g. `mm(X,v,K)` is replaced with v*X/(X+K). Currently supported functions: `mm`, `mmr`, `hill`, `hillr`, and `hill`.
 """
 function expand_registered_functions(expr)
-    istree(expr) || return expr
+    iscall(expr) || return expr
     args = arguments(expr)
     if operation(expr) == Catalyst.mm
         return args[2] * args[1] / (args[1] + args[3])
diff --git a/src/spatial_reaction_systems/utility.jl b/src/spatial_reaction_systems/utility.jl
@@ -212,7 +212,7 @@ end
 # If all parameters the rate depend on are uniform all edges, this becomes a length 1 vector.
 # Else a vector with each value corresponding to the rate at one specific edge.
 function compute_transport_rates(rate_law::Num,
-        p_val_dict::Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Float64}}, num_edges::Int64)
+        p_val_dict::Dict{BasicSymbolic{Real}, Vector{Float64}}, num_edges::Int64)
     # Finds parameters involved in rate and create a function evaluating the rate law.
     relevant_ps = Symbolics.get_variables(rate_law)
     rate_law_func = drop_expr(@RuntimeGeneratedFunction(build_function(
diff --git a/test/dsl/dsl_basic_model_construction.jl b/test/dsl/dsl_basic_model_construction.jl
@@ -2,8 +2,9 @@
 
 # Fetch packages.
 using DiffEqBase, Catalyst, Random, Test
-using ModelingToolkit: operation, istree, get_unknowns, get_ps, get_eqs, get_systems,
+using ModelingToolkit: operation, get_unknowns, get_ps, get_eqs, get_systems,
                        get_iv, nameof
+using Symbolics: iscall
 
 # Sets stable rng number.
 using StableRNGs
@@ -22,7 +23,7 @@ function unpacksys(sys)
     get_eqs(sys), get_iv(sys), get_unknowns(sys), get_ps(sys), nameof(sys), get_systems(sys)
 end
 
-opname(x) = istree(x) ? nameof(operation(x)) : nameof(x)
+opname(x) = iscall(x) ? nameof(operation(x)) : nameof(x)
 alleq(xs, ys) = all(isequal(x, y) for (x, y) in zip(xs, ys))
 
 # Gets all the reactants in a set of equations.
diff --git a/test/miscellaneous_tests/reactionsystem_serialisation.jl b/test/miscellaneous_tests/reactionsystem_serialisation.jl
@@ -4,10 +4,11 @@
 using Catalyst, Test
 using Catalyst: get_rxs
 using ModelingToolkit: getdefault, getdescription, get_metadata
+using Symbolics: getmetadata
 
 # Creates missing getters for MTK metadata (can be removed once added to MTK).
-getmisc(x) = SymbolicUtils.getmetadata(Symbolics.unwrap(x), ModelingToolkit.VariableMisc, nothing)
-getinput(x) = SymbolicUtils.getmetadata(Symbolics.unwrap(x), ModelingToolkit.VariableInput, nothing)
+getmisc(x) = getmetadata(Symbolics.unwrap(x), ModelingToolkit.VariableMisc, nothing)
+getinput(x) = getmetadata(Symbolics.unwrap(x), ModelingToolkit.VariableInput, nothing)
 
 # Sets the default `t` and `D` to use.
 t = default_t()
@@ -191,16 +192,16 @@ let
     @test isequal(getdefault(rs_loaded.rs2.W2), float_md)
 
     # Checks that `Reaction` metadata fields are correct.
-    @test isequal(getmetadata(get_rxs(rs_loaded)[1], :misc), bool_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded)[2], :misc), int_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded)[3], :misc), sym_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded)[4], :misc), str_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded)[5], :misc), nothing_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded.rs2)[1], :misc), expr_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded.rs2)[2], :misc), tup_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded.rs2)[3], :misc), vec_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded.rs2)[4], :misc), dict_md)
-    @test isequal(getmetadata(get_rxs(rs_loaded.rs2)[5], :misc), mat_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded)[1]), bool_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded)[2]), int_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded)[3]), sym_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded)[4]), str_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded)[5]), nothing_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded.rs2)[1]), expr_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded.rs2)[2]), tup_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded.rs2)[3]), vec_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded.rs2)[4]), dict_md)
+    @test isequal(Catalyst.getmisc(get_rxs(rs_loaded.rs2)[5]), mat_md)
 
     # Checks that `ReactionSystem` metadata fields are correct.
     @test isequal(get_metadata(rs_loaded), mat_md)
diff --git a/test/reactionsystem_core/parameter_type_designation.jl b/test/reactionsystem_core/parameter_type_designation.jl
@@ -2,7 +2,7 @@
 
 # Fetch packages.
 using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, StochasticDiffEq, Test
-using Symbolics: unwrap
+using Symbolics: BasicSymbolic, unwrap
 
 # Sets stable rng number.
 using StableRNGs
@@ -47,16 +47,16 @@ end
 
 # Tests that parameters stored in the system have the correct type.
 let
-    @test Symbolics.unwrap(rs.p1) isa SymbolicUtils.BasicSymbolic{Real}
-    @test Symbolics.unwrap(rs.d1) isa SymbolicUtils.BasicSymbolic{Real}
-    @test Symbolics.unwrap(rs.p2) isa SymbolicUtils.BasicSymbolic{Float64}
-    @test Symbolics.unwrap(rs.d2) isa SymbolicUtils.BasicSymbolic{Float64}
-    @test Symbolics.unwrap(rs.p3) isa SymbolicUtils.BasicSymbolic{Int64}
-    @test Symbolics.unwrap(rs.d3) isa SymbolicUtils.BasicSymbolic{Int64}
-    @test Symbolics.unwrap(rs.p4) isa SymbolicUtils.BasicSymbolic{Float32}
-    @test Symbolics.unwrap(rs.d4) isa SymbolicUtils.BasicSymbolic{Rational{Int64}}
-    @test Symbolics.unwrap(rs.p5) isa SymbolicUtils.BasicSymbolic{Rational{Int64}}
-    @test Symbolics.unwrap(rs.d5) isa SymbolicUtils.BasicSymbolic{Float32}
+    @test Symbolics.unwrap(rs.p1) isa BasicSymbolic{Real}
+    @test Symbolics.unwrap(rs.d1) isa BasicSymbolic{Real}
+    @test Symbolics.unwrap(rs.p2) isa BasicSymbolic{Float64}
+    @test Symbolics.unwrap(rs.d2) isa BasicSymbolic{Float64}
+    @test Symbolics.unwrap(rs.p3) isa BasicSymbolic{Int64}
+    @test Symbolics.unwrap(rs.d3) isa BasicSymbolic{Int64}
+    @test Symbolics.unwrap(rs.p4) isa BasicSymbolic{Float32}
+    @test Symbolics.unwrap(rs.d4) isa BasicSymbolic{Rational{Int64}}
+    @test Symbolics.unwrap(rs.p5) isa BasicSymbolic{Rational{Int64}}
+    @test Symbolics.unwrap(rs.d5) isa BasicSymbolic{Float32}
 end
 
 # Tests that simulations with differentially typed variables yields correct results.
diff --git a/test/reactionsystem_core/symbolic_stoichiometry.jl b/test/reactionsystem_core/symbolic_stoichiometry.jl
@@ -2,7 +2,7 @@
 
 # Fetch packages.
 using Catalyst, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Statistics, Test
-using Symbolics: unwrap
+using Symbolics: BasicSymbolic, unwrap
 
 # Sets stable rng number.
 using StableRNGs
@@ -46,8 +46,8 @@ let
     @test rs1 == rs2 == rs3
     @test issetequal(unknowns(rs1), [X, Y])
     @test issetequal(parameters(rs1), [p, k, d, n1, n2, n3])
-    @test unwrap(d) isa SymbolicUtils.BasicSymbolic{Float64}
-    @test unwrap(n1) isa SymbolicUtils.BasicSymbolic{Int64}
+    @test unwrap(d) isa BasicSymbolic{Float64}
+    @test unwrap(n1) isa BasicSymbolic{Int64}
 end
 
 # Declares a network, parameter values, and initial conditions, to be used for the next couple of tests.
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -18,7 +18,7 @@ using SafeTestsets, Test
         @time @safetestset "Symbolic Stoichiometry" begin include("reactionsystem_core/symbolic_stoichiometry.jl") end
         @time @safetestset "Parameter Type Designation" begin include("reactionsystem_core/parameter_type_designation.jl") end
         @time @safetestset "Custom CRN Functions" begin include("reactionsystem_core/custom_crn_functions.jl") end
-        # @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
+        @time @safetestset "Coupled CRN/Equation Systems" begin include("reactionsystem_core/coupled_equation_crn_systems.jl") end
         @time @safetestset "Events" begin include("reactionsystem_core/events.jl") end
     
         # Tests model creation via the @reaction_network DSL.
@@ -74,9 +74,9 @@ using SafeTestsets, Test
         end
 
         # Tests extensions.
-        # @time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end
-        # @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end
-        # @time @safetestset "Structural Identifiability Extension" begin include("extensions/structural_identifiability.jl") end
+        @time @safetestset "BifurcationKit Extension" begin include("extensions/bifurcation_kit.jl") end
+        @time @safetestset "HomotopyContinuation Extension" begin include("extensions/homotopy_continuation.jl") end
+        @time @safetestset "Structural Identifiability Extension" begin include("extensions/structural_identifiability.jl") end
     #end
 
-end # @time
+end # @time
diff --git a/test/spatial_modelling/lattice_reaction_systems.jl b/test/spatial_modelling/lattice_reaction_systems.jl
