[sivasathyaseeelan]

Adds SimpleImplicitTauLeaping, an implicit tau-leaping method for stiff pure-jump problems.

Explicit tau-leaping is limited by the fastest reaction in the system, so a stiff model forces a step size far below the timescale of interest. Each step here instead solves a nonlinear equation for the state at the end of the step, following Rathinam et al. (2003) and Cao et al. (2004), and the leap is then sampled from the propensities at that predicted state.

using JumpProcesses

# stiff reversible isomerization A <-> B
maj = MassActionJump([100.0, 100.0], [[1 => 1], [2 => 1]],
                     [[1 => -1, 2 => 1], [1 => 1, 2 => -1]])
prob = DiscreteProblem([100.0, 100.0], (0.0, 1.0))
jprob = JumpProblem(prob, PureLeaping(), maj)

sol = solve(jprob, SimpleImplicitTauLeaping())                                  # fully implicit
sol = solve(jprob, SimpleImplicitTauLeaping(solver = TrapezoidalImplicitSolver()))

Two formulations are selected with the solver keyword. NewtonImplicitSolver takes the step in its fully implicit form, while TrapezoidalImplicitSolver averages the propensities at the current and new states, which damps the excessive stiffness of the fully implicit step.

Notes on this revision

This has been rebased onto current master, which in the meantime merged SimpleExplicitTauLeaping. The step-size selection (compute_hor, compute_gi, compute_tau), the mass action rate construction and the saveat handling are now shared with it rather than redefined, and validate_pure_leaping_inputs is widened to cover both algorithms instead of gaining a second identical method.

Two problems from the earlier revision are fixed:

• A non-convergent step could spin forever. The retry path halved tau and then continued, but tau is recomputed from compute_tau at the top of each iteration, so the halved value was discarded. The explicit loop survives the same shape because its retries redraw the Poisson counts, but the implicit solve is deterministic in (u_current, tau) and would fail identically on every pass. The bound is now carried across iterations, and once tau reaches dtmin a non-convergent step is reported as an error rather than silently accepted.
• Setfield was declared as a dependency but never used, which Aqua's stale-dependency check flags.

The algorithm type is also now fully parameterized (SimpleImplicitTauLeaping{T, S}) rather than holding an abstractly-typed solver field, and the residual evaluation reuses preallocated propensity caches instead of allocating on each call.

Testing

The existing SIR and SEIR mean-trajectory comparisons in test/regular_jumps.jl are extended to both implicit formulations, checked against the exact SSA solution at every saved timepoint.

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

[ChrisRackauckas]

This is weird, just make it SimpleImplicitTauLeaping and SimpleTrapezoidalLeaping

[ChrisRackauckas]

Pretty much done, I've got it from here.

[ChrisRackauckas]

I took the commit from here and finished it in #643

[ChrisRackauckas]

Thanks! Your commit is credited with your authorship for doing like 99% of the work.

[sivasathyaseeelan]

Thanks Chris, and thanks for finishing it in #643.