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/JumpProcesses.jl

@@ -4,7 +4,7 @@ using Reexport
 @reexport using DiffEqBase
 
 using LinearAlgebra, Markdown, DocStringExtensions
-using NonlinearSolve
+using SimpleNonlinearSolve
 using DataStructures, PoissonRandom, Random, ArrayInterface
 using FunctionWrappers, UnPack
 using Graphs

src/simple_regular_solve.jl

@@ -251,7 +251,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,90 +1,35 @@
 using JumpProcesses, DiffEqBase
-using Test, LinearAlgebra, Statistics
-using StableRNGs
+using Test, LinearAlgebra
+using StableRNGs, Plots
 rng = StableRNG(12345)
 
-Nsims = 8000
-
-# 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
-
-    u0 = [999.0, 10.0, 0.0]  # S, I, R
-    tspan = (0.0, 250.0)
-
-    prob_disc = DiscreteProblem(u0, tspan, p)
-    rj = RegularJump(regular_rate, regular_c, 3)
-    jump_prob = JumpProblem(prob_disc, Direct(), rj)
-
-    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)
-
-    sol = solve(EnsembleProblem(jump_prob), SimpleImplicitTauLeaping(), EnsembleSerial(); trajectories = Nsims)
-    mean_implicit = mean(sol.u[i][1,end] for i in 1:Nsims)
-
-    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)
+# Parameters
+c1 = 1.0      # S1 -> 0
+c2 = 10.0     # S1 + S1 <- S2
+c3 = 1000.0   # S1 + S1 -> S2
+c4 = 0.1      # S2 -> S3
+p = (c1, c2, c3, c4)
+
+regular_rate_implicit = (out, u, p, t) -> begin
+    out[1] = p[1] * u[1]          # S1 -> 0
+    out[2] = p[2] * u[2]          # S1 + S1 <- S2
+    out[3] = p[3] * u[1] * (u[1] - 1) / 2  # S1 + S1 -> S2
+    out[4] = p[4] * u[2]          # S2 -> S3
 end
 
+regular_c_implicit = (dc, u, p, t, counts, mark) -> begin
+    dc .= 0.0
+    dc[1] = -counts[1] - 2 * counts[3] + 2 * counts[2]  # S1: -decay - 2*forward + 2*backward
+    dc[2] = counts[3] - counts[2] - counts[4]           # S2: +forward - backward - decay
+    dc[3] = counts[4]                                   # S3: +decay
+end
 
-# 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)
+u0 = [10000.0, 0.0, 0.0]  # S1, S2, S3
+tspan = (0.0, 4.0)
 
-    @test isapprox(mean_simple, mean_implicit, rtol=0.05)
-    @test isapprox(mean_simple, mean_adaptive, rtol=0.05)
-end
+# Create JumpProblem with proper parameter passing
+prob_disc = DiscreteProblem(u0, tspan, p)
+rj = RegularJump(regular_rate_implicit, regular_c_implicit, 4)
+jump_prob = JumpProblem(prob_disc, Direct(), rj; rng=StableRNG(12345))
+sol = solve(jump_prob, SimpleAdaptiveTauLeaping())
+plot(sol)