Merge pull request #480 from isaacsas/neg_massaction_rate_fix

fix neg mass action rates with FP types

Author: isaacsas

src/massaction_rates.jl

@@ -3,9 +3,8 @@
 # stochiometric coefficient.
 ###############################################################################
 
-@inline function evalrxrate(speciesvec::AbstractVector{T}, rxidx::S,
-        majump::MassActionJump{U, V, W, X})::R where
-        {T, S, R, U <: AbstractVector{R}, V, W, X}
+@inline function evalrxrate(speciesvec::AbstractVector{T}, rxidx,
+        majump::MassActionJump{U})::R where {T <: Integer, R, U <: AbstractVector{R}}
     val = one(T)
     @inbounds for specstoch in majump.reactant_stoch[rxidx]
         specpop = speciesvec[specstoch[1]]
@@ -19,6 +18,25 @@
     @inbounds return val * majump.scaled_rates[rxidx]
 end
 
+@inline function evalrxrate(speciesvec::AbstractVector{T}, rxidx,
+        majump::MassActionJump{U})::R where {T <: Real, R, U <: AbstractVector{R}}
+    val = one(T)
+    @inbounds for specstoch in majump.reactant_stoch[rxidx]
+        specpop = speciesvec[specstoch[1]]
+        val *= specpop
+        @inbounds for k in 2:specstoch[2]
+            specpop -= one(specpop)
+            val *= specpop
+        end
+        # we need to check the smallest rate law term is positive
+        # i.e. for an order k reaction: x - k + 1 > 0 
+        (specpop <= 0) && return zero(R)
+    end
+
+    @inbounds return val * majump.scaled_rates[rxidx]
+end
+
+
 @inline function executerx!(speciesvec::AbstractVector{T}, rxidx::S,
         majump::M) where {T, S, M <: AbstractMassActionJump}
     @inbounds net_stoch = majump.net_stoch[rxidx]

test/fp_unknowns.jl

@@ -0,0 +1,52 @@
+using Test, JumpProcesses, StableRNGs
+
+rng = StableRNG(12345)
+
+# test for https://github.com/SciML/JumpProcesses.jl/issues/479
+#    A, ∅ --> X
+#    1, 2X + Y --> 3X
+#    B, X --> Y
+#    1, X --> ∅
+function test(rng)
+    # dep graphs
+    dg =  [[1, 2, 3, 4],
+            [2, 3, 4],
+            [2, 3, 4],
+            [2, 3, 4]]
+    vtoj = [[2, 3, 4],
+            [2]]
+    jtov = [[1],
+            [1, 2],
+            [1, 2],
+            [1]]
+
+    # reaction as MassActionJump
+    sr = [1.0, .5, 4.0, 1.0]
+    rs = [Pair{Int,Int}[],[1 => 2, 2 => 1], [1 => 1], [1 => 1]]
+    ns = [[1 => 1], [1 => 1, 2 => -1], [1 => -1, 2 => 1], [1 => -1]]
+    maj = MassActionJump(sr, rs, ns; scale_rates = false)
+
+    Nsims = 10000
+    u0 = [1.0, 10.0]
+    tspan = (0.0, 50.0)    
+    dprob = DiscreteProblem(u0, tspan)
+    SSAalgs = JumpProcesses.JUMP_AGGREGATORS
+    Xmeans = zeros(length(SSAalgs))
+    Ymeans = zeros(length(SSAalgs))
+    for (j, agg) in enumerate(SSAalgs)
+        jprob = JumpProblem(dprob, agg, maj; save_positions = (false, false), rng,
+            vartojumps_map = vtoj, jumptovars_map = jtov, dep_graph = dg, 
+            scale_rates = false)
+        for i in 1:Nsims
+            sol = solve(jprob, SSAStepper())
+            Xmeans[j] += sol[1,end]
+            Ymeans[j] += sol[2,end]
+        end
+    end
+    Xmeans ./= Nsims
+    Ymeans ./= Nsims
+    # for i in 2:length(SSAalgs)
+    #     @test abs(Xmeans[i] - Xmeans[1]) < (.1 * Xmeans[1])
+    #     @test abs(Ymeans[i] - Ymeans[1]) < (.1 * Ymeans[1])
+    # end
+end

test/runtests.jl

@@ -16,6 +16,7 @@ using JumpProcesses, DiffEqBase, SafeTestsets
     @time @safetestset "Mass Action Jump Tests; Gene Expr Model" begin include("geneexpr_test.jl") end
     @time @safetestset "Mass Action Jump Tests; Nonlinear Rx Model" begin include("bimolerx_test.jl") end
     @time @safetestset "Mass Action Jump Tests; Special Cases" begin include("degenerate_rx_cases.jl") end
+    @time @safetestset "Mass Action Jump Tests; Floating Point Inputs" begin include("fp_unknowns.jl") end
     @time @safetestset "Direct allocations test" begin include("allocations.jl") end
     @time @safetestset "Bracketing Tests" begin include("bracketing.jl") end
     @time @safetestset "Composition-Rejection Table Tests" begin include("table_test.jl") end