still bug hunting

Author: isaacsas

test/jumpsystem.jl

@@ -1,7 +1,7 @@
 using ModelingToolkit, DiffEqBase, DiffEqJump, Test, LinearAlgebra
 MT = ModelingToolkit
 
-# basic SIR model with tweaks
+# basic MT SIR model with tweaks
 @parameters β γ t
 @variables S I R
 rate₁   = β*S*I
@@ -32,25 +32,24 @@ jump2 = VariableRateJump(rate2,affect2!)
 u = [100, 9, 5]
 p = (0.1/1000,0.01)
 tf = 1.0
-mutable struct TestInt
-    u
-    p
-    t
+mutable struct TestInt{U,V,T}
+    u::U
+    p::V
+    t::T
 end
 mtintegrator = TestInt(u,p,tf)
 integrator   = TestInt(u,p,tf)
 @test abs(mtjump1.rate(u,p,tf) - jump1.rate(u,p,tf)) < 10*eps()
 @test abs(mtjump2.rate(u,p,tf) - jump2.rate(u,p,tf)) < 10*eps()
 mtjump1.affect!(mtintegrator)
 jump1.affect!(integrator)
-@test norm(integrator.u - mtintegrator.u) < 10*eps()
+@test all(integrator.u .== mtintegrator.u) 
 mtintegrator.u .= u; integrator.u .= u
 mtjump2.affect!(mtintegrator)
 jump2.affect!(integrator)
-@test norm(integrator.u - mtintegrator.u) < 10*eps()
+@test all(integrator.u .== mtintegrator.u)
 
-
-# test can make and solve a jump problem
+# test MT can make and solve a jump problem
 rate₃   = γ*I
 affect₃ = [I ~ I - 1, R ~ R + 1]
 j₃ = ConstantRateJump(rate₃,affect₃)
@@ -60,5 +59,38 @@ dprob = DiscreteProblem(u₀,tspan,p)
 jprob = JumpProblem(js2, dprob, Direct())
 sol = solve(jprob, SSAStepper())
 
-using Plots
-plot(sol)
+# test the MT JumpProblem rates/affects are correct
+rate2(u,p,t) = 0.01u[2]
+jump2 = ConstantRateJump(rate2,affect2!)
+mtjumps = jprob.discrete_jump_aggregation
+@test abs(mtjumps.rates[1](u,p,tf) - jump1.rate(u,p,tf)) < 10*eps()
+@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) 
+mtintegrator.u .= u; integrator.u .= u
+mtjumps.affects![2](mtintegrator)
+jump2.affect!(integrator)
+@test all(integrator.u .== mtintegrator.u)
+
+# # direct vers
+# p = (0.1/1000,0.01)
+# prob = DiscreteProblem([999,1,0],(0.0,250.0),p)
+# r1(u,p,t) = (0.1/1000.0)*u[1]*u[2]
+# function a1!(integrator)
+#   integrator.u[1] -= 1
+#   integrator.u[2] += 1
+# end
+# j1 = ConstantRateJump(r1,a1!)
+# r2(u,p,t) = 0.01u[2]
+# function a2!(integrator)
+#   integrator.u[2] -= 1
+#   integrator.u[3] += 1
+# end
+# j2 = ConstantRateJump(r2,a2!)
+# jset = JumpSet((),(j1,j2),nothing,nothing)
+# jprob = JumpProblem(prob,Direct(),jset)
+# sol = solve(jprob, SSAStepper())
+
+# using Plots
+# plot(sol)