ImplicitTauLeaping solver and Kernels

[sivasathyaseeelan]

using JumpProcesses, DiffEqBase
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

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

dc = zeros(3, 2)

rj = RegularJump(regular_rate, regular_c, dc; constant_c = true)
jumps = JumpSet(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, ImplicitTauLeaping())

const _dc = zeros(3, 2)
dc[1, 1] = -1
dc[2, 1] = 1
dc[2, 2] = -1
dc[3, 2] = 1

function regular_c(du, u, p, t, counts, mark)
    mul!(du, dc, counts)
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, ImplicitTauLeaping())
plot(sol)

using JumpProcesses, DiffEqBase
using Test, LinearAlgebra
using StableRNGs, Plots
rng = StableRNG(12345)

# 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

u0 = [10000.0, 0.0, 0.0]  # S1, S2, S3
tspan = (0.0, 4.0)

# 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, ImplicitTauLeaping())
plot(sol)

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[ChrisRackauckas]

SimpleNonlinearSolve.jl