From bcba92ff5721822855739fbb5778f8b0d6223968 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 22 Jul 2026 15:25:33 -0400 Subject: [PATCH] Retire unreachable symbolic compatibility ranges Co-Authored-By: Chris Rackauckas --- Project.toml | 16 ++++++++------ src/FiniteStateProjection.jl | 21 +++++++++--------- src/build_rhs.jl | 8 +++---- src/build_rhs_ss.jl | 2 +- src/fspsystem.jl | 32 +++++++++++++++++++-------- src/indexhandlers.jl | 3 ++- test/symbolic_api.jl | 42 ++++++++++++++++++++++++++++++++++++ 7 files changed, 91 insertions(+), 33 deletions(-) create mode 100644 test/symbolic_api.jl diff --git a/Project.toml b/Project.toml index 7bf647c..cf1e695 100644 --- a/Project.toml +++ b/Project.toml @@ -14,22 +14,23 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" [compat] -Catalyst = "15, 16.2" +Catalyst = "16.2" Distributions = "0.25" LinearAlgebra = "1" -ModelingToolkitBase = "1" +ModelingToolkitBase = "1.17" OrdinaryDiffEq = "6, 7" Random = "1" RuntimeGeneratedFunctions = "0.5" SafeTestsets = "0.1, 1" -SciMLBase = "2, 3" +SciMLBase = "2.144, 3" SciMLTesting = "2.4" SparseArrays = "1" SteadyStateDiffEq = "2.5" Sundials = "4, 5, 6" -SymbolicIndexingInterface = "0.3" -Symbolics = "5, 6, 7" -SymbolicUtils = "1, 2, 3, 4" +SymbolicIndexingInterface = "0.3.43" +Symbolics = "7.13" +SymbolicUtils = "4.18" +TOML = "1" Test = "1" julia = "1.10" @@ -42,7 +43,8 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" +TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Distributions", "LinearAlgebra", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLTesting", "SteadyStateDiffEq", "Sundials", "Test"] +test = ["Distributions", "LinearAlgebra", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLTesting", "SteadyStateDiffEq", "Sundials", "TOML", "Test"] diff --git a/src/FiniteStateProjection.jl b/src/FiniteStateProjection.jl index 02697b8..5399e43 100644 --- a/src/FiniteStateProjection.jl +++ b/src/FiniteStateProjection.jl @@ -1,8 +1,8 @@ module FiniteStateProjection -using Catalyst: Reaction, ReactionSystem, jumpratelaw, netstoichmat, numspecies, parameters, - reactions, species -using ModelingToolkitBase: equations, get_systems +using Catalyst: Reaction, ReactionSystem, jumpratelaw, netstoichmat, numspecies, reactions, + species +using ModelingToolkitBase: equations, get_systems, parameters import RuntimeGeneratedFunctions using RuntimeGeneratedFunctions: @RuntimeGeneratedFunction import SciMLBase @@ -10,9 +10,8 @@ using SciMLBase: ODEFunction, ODEProblem, SteadyStateProblem import SparseArrays using SparseArrays: sparse using SymbolicIndexingInterface: getname -using Symbolics: @variables, value -using SymbolicUtils: scalarize, substitute -using SymbolicUtils.Code: toexpr +using Symbolics: @variables, build_function +using SymbolicUtils: substitute import Base: LinearIndices, vec @@ -21,11 +20,7 @@ RuntimeGeneratedFunctions.init(@__MODULE__) export FSPSystem, DefaultIndexHandler, SteadyState _unresolve1(x) = x -@static if VERSION >= v"1.12-" - _unresolve1(f::Function) = nameof(f) -else - _unresolve1(f::Function) = methods(f).mt.name -end +_unresolve1(f::Function) = nameof(f) _prewalk(f, x) = _prewalk_children(f(x), f) _prewalk_children(x, f) = x @@ -33,6 +28,10 @@ _prewalk_children(x::Expr, f) = Expr(x.head, map(arg -> _prewalk(f, arg), x.args function _striplines(ex) ex isa Expr || return ex + if ex.head === :macrocall + # A macro call requires its LineNumberNode in the second slot. + return Expr(:macrocall, ex.args[1:2]..., _striplines.(ex.args[3:end])...) + end args = (arg for arg in ex.args if !(arg isa LineNumberNode)) return Expr(ex.head, _striplines.(args)...) end diff --git a/src/build_rhs.jl b/src/build_rhs.jl index 5aa69f8..c53b202 100644 --- a/src/build_rhs.jl +++ b/src/build_rhs.jl @@ -9,7 +9,7 @@ the variable `psym` represents an `AbstractVector`. See also: [`build_rhs_header`](@ref), [`build_rhs`](@ref) """ function unpackparams(sys::FSPSystem, psym::Symbol) - param_names = Expr(:tuple, map(par -> par.name, parameters(sys.rs))...) + param_names = Expr(:tuple, getname.(parameters(sys.rs))...) return quote $(param_names) = $(psym) @@ -47,8 +47,8 @@ See also: [`build_rhs`](@ref) function build_rhs_firstpass(sys::FSPSystem) isempty(sys.rfs) && return quote end - first_line = :(du[idx_in] = -u[idx_in] * $(sys.rfs[1].body)) - other_lines = (:(du[idx_in] -= u[idx_in] * $(rf.body)) for rf in sys.rfs[2:end]) + first_line = :(du[idx_in] = -u[idx_in] * $(sys.rfs[1].expression)) + other_lines = (:(du[idx_in] -= u[idx_in] * $(rf.expression)) for rf in sys.rfs[2:end]) return quote for idx_in in singleindices($(sys.ih), u) @@ -81,7 +81,7 @@ function build_rhs_secondpass(sys::FSPSystem) ex = quote for (idx_in, idx_out) in pairedindices($(sys.ih), u, $(CartesianIndex(S[:, i]...))) - du[idx_out] += u[idx_in] * $(rf.body) + du[idx_out] += u[idx_in] * $(rf.expression) end end diff --git a/src/build_rhs_ss.jl b/src/build_rhs_ss.jl index a43bb76..eb4a8a5 100644 --- a/src/build_rhs_ss.jl +++ b/src/build_rhs_ss.jl @@ -40,7 +40,7 @@ function build_rhs_singlepass_ss(sys::FSPSystem) ex = quote for (idx_in, idx_out) in pairedindices($(sys.ih), u, $(CartesianIndex(S[:, i]...))) - rate = u[idx_in] * $(rf.body) + rate = u[idx_in] * $(rf.expression) du[idx_in] -= rate du[idx_out] += rate end diff --git a/src/fspsystem.jl b/src/fspsystem.jl index 8944954..c433335 100644 --- a/src/fspsystem.jl +++ b/src/fspsystem.jl @@ -28,6 +28,13 @@ struct FSPSystem{IHT <: AbstractIndexHandler, RT} rfs::RT end +struct RateFunction{F, E} + callable::F + expression::E +end + +@inline (rf::RateFunction)(args...) = rf.callable(args...) + function FSPSystem( rs::ReactionSystem, ih::AbstractIndexHandler = DefaultIndexHandler{length(species(rs))}(); @@ -58,36 +65,43 @@ function build_ratefuncs( rs::ReactionSystem, ih::AbstractIndexHandler; state_sym::Symbol, combinatoric_ratelaw::Bool = true ) + nspecs = numspecies(rs) + state = (@variables ($state_sym)[1:nspecs])[1] + @variables t + params = parameters(rs) substitutions = getsubstitutions(ih, rs, state_sym = state_sym) return map(reactions(rs)) do reac jrl = jumpratelaw(reac; combinatoric_ratelaw) - jrl_s = substitute(jrl, substitutions) - toexpr(jrl_s) + rate = substitute(jrl, substitutions) + ex = build_function(rate, state, t, params...; expression = Val{true}) + ex isa Expr && ex.head === :function && length(ex.args) == 2 || + throw(ArgumentError("Symbolics.build_function returned an unsupported expression")) + ex.args[2] end end function create_ratefuncs(rs::ReactionSystem, ih::AbstractIndexHandler; combinatoric_ratelaw::Bool = true) - paramsyms = Symbol.(parameters(rs)) + params = getname.(parameters(rs)) return tuple( map( - ex -> compile_ratefunc(ex, paramsyms), + rate -> compile_ratefunc(rate, params), build_ratefuncs(rs, ih; state_sym = :idx_in, combinatoric_ratelaw) )... ) end -function compile_ratefunc(ex_rf, params) - ex = _flatten(:((idx_in, t, $(params...)) -> $(ex_rf))) - return @RuntimeGeneratedFunction(ex) +function compile_ratefunc(rate, params) + ex = _flatten(:((idx_in, t, $(params...)) -> $(rate))) + return RateFunction(@RuntimeGeneratedFunction(ex), rate) end _parameter_symbol(key::Symbol) = key -_parameter_symbol(key) = Symbol(value(key)) +_parameter_symbol(key) = getname(key) function pmap_to_p(sys::FSPSystem, pmap) pmap isa SciMLBase.NullParameters && return pmap values_by_parameter = Dict(_parameter_symbol(k) => v for (k, v) in pmap) - return [values_by_parameter[p] for p in Symbol.(parameters(sys.rs))] + return [values_by_parameter[p] for p in getname.(parameters(sys.rs))] end diff --git a/src/indexhandlers.jl b/src/indexhandlers.jl index 602aa3d..8e6a7ce 100644 --- a/src/indexhandlers.jl +++ b/src/indexhandlers.jl @@ -151,7 +151,8 @@ Defines the abundance of species ``S_i`` to be `state_sym[i] - offset`. """ function getsubstitutions(ih::DefaultIndexHandler, rs::ReactionSystem; state_sym::Symbol) nspecs = numspecies(rs) - state_sym_vec = value.(scalarize((@variables ($state_sym)[1:nspecs])[1])) + state_array = (@variables ($state_sym)[1:nspecs])[1] + state_sym_vec = [state_array[i] for i in 1:nspecs] species_orig = species(rs) species_perm = [species_orig[ih.perm[i]] for i in 1:nspecs] diff --git a/test/symbolic_api.jl b/test/symbolic_api.jl new file mode 100644 index 0000000..fccf6a0 --- /dev/null +++ b/test/symbolic_api.jl @@ -0,0 +1,42 @@ +using FiniteStateProjection +using Catalyst +using SciMLBase: ODEProblem +using SparseArrays: SparseMatrixCSC +using Symbolics: @variables +using Test +using TOML + +@testset "Symbolic API and compatibility floors" begin + @variables p + + @test FiniteStateProjection._parameter_symbol(p) == :p + + project = TOML.parsefile(joinpath(pkgdir(FiniteStateProjection), "Project.toml")) + compat = project["compat"] + @test compat["Catalyst"] == "16.2" + @test compat["ModelingToolkitBase"] == "1.17" + @test compat["SciMLBase"] == "2.144, 3" + @test compat["SymbolicIndexingInterface"] == "0.3.43" + @test compat["Symbolics"] == "7.13" + @test compat["SymbolicUtils"] == "4.18" + + rs = @reaction_network begin + birth, 0 --> A + death, A --> 0 + end + sys = FSPSystem(rs) + @test all( + ex -> ex isa Expr, + FiniteStateProjection.build_ratefuncs(rs, sys.ih; state_sym = :idx_in) + ) + @test sys.rfs[1]([2], 0.0, 3.0, 4.0) == 3.0 + @test sys.rfs[2]([2], 0.0, 3.0, 4.0) == 4.0 + + pmap = [:birth => 3.0, :death => 4.0] + u = [0.25, 0.5, 0.25] + du = similar(u) + prob = ODEProblem(sys, u, (0.0, 1.0), pmap) + prob.f(du, u, last.(pmap), 0.0) + A = SparseMatrixCSC(sys, (length(u),), pmap, 0) + @test du == A * u +end