update project.toml

Author: sivasathyaseeelan

Project.toml

@@ -13,13 +13,13 @@ FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
-NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 PoissonRandom = "e409e4f3-bfea-5376-8464-e040bb5c01ab"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"

src/simple_regular_solve.jl

@@ -263,7 +263,7 @@ function implicit_tau_step(u_prev, t_prev, tau, rate_cache, counts, nu, p, rate,
 
     # Solve the nonlinear system
     prob = NonlinearProblem(f, u_new, nothing)
-    sol = solve(prob, NewtonRaphson())
+    sol = solve(prob, SimpleNewtonRaphson())
 
     # Check for convergence and numerical stability
     if sol.retcode != ReturnCode.Success || any(isnan.(sol.u)) || any(isinf.(sol.u))

test/regular_jumps.jl

@@ -1,41 +1,90 @@
 using JumpProcesses, DiffEqBase
-using Test, LinearAlgebra, Statistics
-using StableRNGs
+using Test, LinearAlgebra
+using StableRNGs, Plots
 rng = StableRNG(12345)
 
-function regular_rate(out, u, p, t)
-    out[1] = (0.1 / 1000.0) * u[1] * u[2]
-    out[2] = 0.01u[2]
-end
+Nsims = 8000
 
-function regular_c(dc, u, p, t, mark)
-    dc[1, 1] = -1
-    dc[2, 1] = 1
-    dc[2, 2] = -1
-    dc[3, 2] = 1
-end
+# SIR model with influx
+let
+    β = 0.1 / 1000.0
+    ν = 0.01
+    influx_rate = 1.0
+    p = (β, ν, influx_rate)
+
+    regular_rate = (out, u, p, t) -> begin
+        out[1] = p[1] * u[1] * u[2]  # β*S*I (infection)
+        out[2] = p[2] * u[2]         # ν*I (recovery)
+        out[3] = p[3]                # influx_rate
+    end
+
+    regular_c = (dc, u, p, t, counts, mark) -> begin
+        dc .= 0.0
+        dc[1] = -counts[1] + counts[3]  # S: -infection + influx
+        dc[2] = counts[1] - counts[2]   # I: +infection - recovery
+        dc[3] = counts[2]               # R: +recovery
+    end
 
-dc = zeros(3, 2)
+    u0 = [999.0, 10.0, 0.0]  # S, I, R
+    tspan = (0.0, 250.0)
 
-rj = RegularJump(regular_rate, regular_c, dc; constant_c = true)
-jumps = JumpSet(rj)
+    prob_disc = DiscreteProblem(u0, tspan, p)
+    rj = RegularJump(regular_rate, regular_c, 3)
+    jump_prob = JumpProblem(prob_disc, Direct(), rj)
 
-prob = DiscreteProblem([999.0, 1.0, 0.0], (0.0, 250.0))
-jump_prob = JumpProblem(prob, Direct(), rj; rng = rng)
-sol = solve(jump_prob, SimpleTauLeaping(); dt = 1.0)
+    sol = solve(EnsembleProblem(jump_prob), SimpleTauLeaping(), EnsembleSerial(); trajectories = Nsims, dt = 1.0)
+    mean_simple = mean(sol.u[i][1,end] for i in 1:Nsims)
 
-const _dc = zeros(3, 2)
-dc[1, 1] = -1
-dc[2, 1] = 1
-dc[2, 2] = -1
-dc[3, 2] = 1
+    sol = solve(EnsembleProblem(jump_prob), SimpleImplicitTauLeaping(), EnsembleSerial(); trajectories = Nsims)
+    mean_implicit = mean(sol.u[i][1,end] for i in 1:Nsims)
 
-function regular_c(du, u, p, t, counts, mark)
-    mul!(du, dc, counts)
+    sol = solve(EnsembleProblem(jump_prob), SimpleAdaptiveTauLeaping(), EnsembleSerial(); trajectories = Nsims)
+    mean_adaptive = mean(sol.u[i][1,end] for i in 1:Nsims)
+
+    @test isapprox(mean_simple, mean_implicit, rtol=0.05)
+    @test isapprox(mean_simple, mean_adaptive, rtol=0.05)
 end
 
-rj = RegularJump(regular_rate, regular_c, 2)
-jumps = JumpSet(rj)
-prob = DiscreteProblem([999, 1, 0], (0.0, 250.0))
-jump_prob = JumpProblem(prob, Direct(), rj; rng = rng)
-sol = solve(jump_prob, SimpleTauLeaping(); dt = 1.0)
+
+# SEIR model with exposed compartment
+let
+    β = 0.3 / 1000.0
+    σ = 0.2
+    ν = 0.01
+    p = (β, σ, ν)
+
+    regular_rate = (out, u, p, t) -> begin
+        out[1] = p[1] * u[1] * u[3]  # β*S*I (infection)
+        out[2] = p[2] * u[2]         # σ*E (progression)
+        out[3] = p[3] * u[3]         # ν*I (recovery)
+    end
+
+    regular_c = (dc, u, p, t, counts, mark) -> begin
+        dc .= 0.0
+        dc[1] = -counts[1]           # S: -infection
+        dc[2] = counts[1] - counts[2] # E: +infection - progression
+        dc[3] = counts[2] - counts[3] # I: +progression - recovery
+        dc[4] = counts[3]            # R: +recovery
+    end
+
+    # Initial state
+    u0 = [999.0, 0.0, 10.0, 0.0]  # S, E, I, R
+    tspan = (0.0, 250.0)
+
+    # Create JumpProblem
+    prob_disc = DiscreteProblem(u0, tspan, p)
+    rj = RegularJump(regular_rate, regular_c, 3)
+    jump_prob = JumpProblem(prob_disc, Direct(), rj; rng=StableRNG(12345))
+
+    sol = solve(EnsembleProblem(jump_prob), SimpleTauLeaping(), EnsembleSerial(); trajectories = Nsims, dt = 1.0)
+    mean_simple = mean(sol.u[i][end,end] for i in 1:Nsims)
+
+    sol = solve(EnsembleProblem(jump_prob), SimpleImplicitTauLeaping(), EnsembleSerial(); trajectories = Nsims)
+    mean_implicit = mean(sol.u[i][end,end] for i in 1:Nsims)
+
+    sol = solve(EnsembleProblem(jump_prob), SimpleAdaptiveTauLeaping(), EnsembleSerial(); trajectories = Nsims)
+    mean_adaptive = mean(sol.u[i][end,end] for i in 1:Nsims)
+
+    @test isapprox(mean_simple, mean_implicit, rtol=0.05)
+    @test isapprox(mean_simple, mean_adaptive, rtol=0.05)
+end