Merge pull request #354 from isaacsas/jumpsys-maj-fix

Fix MassActionJump parameter evaluation

Author: ChrisRackauckas

Project.toml

@@ -15,6 +15,7 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

src/ModelingToolkit.jl

@@ -13,6 +13,8 @@ import GeneralizedGenerated
 using DocStringExtensions
 using Base: RefValue
 
+using RecursiveArrayTools
+
 import SymbolicUtils
 import SymbolicUtils: to_symbolic, FnType
 

src/systems/jumps/jumpsystem.jl

@@ -43,10 +43,10 @@ function assemble_crj(js, crj, statetoid)
 end
 
 function assemble_maj(js, maj::MassActionJump{U,Vector{Pair{V,W}},Vector{Pair{V2,W2}}},
-                      statetoid, ptoid, parammap) where {U,V,W,V2,W2}
+                      statetoid, parammap) where {U,V,W,V2,W2}
     sr = maj.scaled_rates
     if sr isa Operation
-        pval = substitute(sr,parammap)
+        pval = simplify(substitute(sr,parammap)).value
     elseif sr isa Variable
         pval = Dict(parammap)[sr()]
     else
@@ -101,7 +101,7 @@ function DiffEqJump.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
     majs = Vector{MassActionJump}()
     pvars = parameters(js)
     statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
-    parammap = map(Pair,pvars,prob.p)
+    parammap = map((x,y)->Pair(x(),y),pvars,prob.p)
 
     for j in equations(js)
         if j isa ConstantRateJump

test/jumpsystem.jl

@@ -44,7 +44,7 @@ integrator   = TestInt(u,p,tf)
 @test abs(mtjump2.rate(u,p,tf) - jump2.rate(u,p,tf)) < 10*eps()
 mtjump1.affect!(mtintegrator)
 jump1.affect!(integrator)
-@test all(integrator.u .== mtintegrator.u) 
+@test all(integrator.u .== mtintegrator.u)
 mtintegrator.u .= u; integrator.u .= u
 mtjump2.affect!(mtintegrator)
 jump2.affect!(integrator)
@@ -60,7 +60,7 @@ u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => .1/1000, γ => .01]
 dprob = DiscreteProblem(js2, u₀map, tspan, parammap)
 jprob = JumpProblem(js2, dprob, Direct(), save_positions=(false,false))
-Nsims = 10000
+Nsims = 30000
 function getmean(jprob,Nsims)
   m = 0.0
   for i = 1:Nsims
@@ -79,7 +79,7 @@ mtjumps = jprob.discrete_jump_aggregation
 @test abs(mtjumps.rates[2](u,p,tf) - jump2.rate(u,p,tf)) < 10*eps()
 mtjumps.affects![1](mtintegrator)
 jump1.affect!(integrator)
-@test all(integrator.u .== mtintegrator.u) 
+@test all(integrator.u .== mtintegrator.u)
 mtintegrator.u .= u; integrator.u .= u
 mtjumps.affects![2](mtintegrator)
 jump2.affect!(integrator)
@@ -107,7 +107,6 @@ m2 = getmean(jprob,Nsims)
 # test JumpSystem solution agrees with direct version
 @test abs(m-m2)/m < .01
 
-
 # mass action jump tests for SIR model
 maj1 = MassActionJump(2*β/2, [S => 1, I => 1], [S => -1, I => 1])
 maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
@@ -130,7 +129,7 @@ jprob = JumpProblem(js4, dprob, Direct())
 m4 = getmean(jprob,Nsims)
 @test abs(m4 - 2.0/.01)*.01/2.0 < .01
 
-# test second order rx runs 
+# test second order rx runs
 maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
 maj2 = MassActionJump(γ, [S => 2], [S => -1])
 js4   = JumpSystem([maj1,maj2], t, [S], [β,γ])

test/runtests.jl

@@ -12,6 +12,7 @@ using SafeTestsets, Test
 @safetestset "NonlinearSystem Test" begin include("nonlinearsystem.jl") end
 @safetestset "OptimizationSystem Test" begin include("optimizationsystem.jl") end
 @safetestset "ReactionSystem Test" begin include("reactionsystem.jl") end
+@safetestset "JumpSystem Test" begin include("jumpsystem.jl") end
 @safetestset "Build Targets Test" begin include("build_targets.jl") end
 @safetestset "Domain Test" begin include("domains.jl") end
 @safetestset "Constraints Test" begin include("constraints.jl") end