Retire unreachable symbolic compatibility ranges

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

src/FiniteStateProjection.jl

@@ -1,18 +1,17 @@
 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
 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,18 +20,18 @@ 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
 _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

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
 

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

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

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]

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