Add @poissonians macro

[isaacsas]

Note: this is very WIP as I've spent a good amount of time iterating with Claude on this, but haven't yet done extensive code review of what it has generated. I'm sure there are places that need updating / changes for various reasons. I'm putting this up now in case anyone has feedback on the basic design (I will stop working on it if there issues with this approach). I will go through this more carefully and drop the WIP once I have carefully reviewed all the code changes.

Longer-term we can expand from this to do better jump classification (i.e. this does the simple thing, and only has jumps with parameters as their rate become a ConstantRateJump, with everything else VariableRateJump). We can also expand this to support marks and/or Poisson random measures, both of which would build from this anyways.

Add @poissonians Macro for Inline Poisson Process Differentials

Summary

This PR introduces the @poissonians macro, enabling inline Poisson process differentials in ModelingToolkitBase equations. Poissonian variables represent Poisson counting process increments and are automatically converted to ConstantRateJump or VariableRateJump objects during system compilation.

API

Declaring Poissonians

@parameters β γ
@variables S(t) I(t) R(t)

# Inline syntax
@poissonians dN_inf(β * S * I) dN_rec(γ * I)

# Block syntax
@poissonians begin
    dN_inf(β * S * I)
    dN_rec(γ * I)
end

Using in Equations

eqs = [
    D(S) ~ -dN_inf,
    D(I) ~ dN_inf - dN_rec,
    D(R) ~ dN_rec
]

@named sys = System(eqs, t)
compiled_sys = mtkcompile(sys)

Querying Poissonians

ispoissonian(dN)              # true/false
getpoissonianrate(dN)         # returns the rate expression
ModelingToolkitBase.poissonians(sys)  # all poissonians in system

Controlling save_positions for VRJs

# Pass save_positions through mtkcompile to control VRJ saving behavior
compiled_sys = mtkcompile(sys; save_positions = (true, true))

Implementation Details

Phase 1: Macro and Metadata

• Added POISSONIAN to VariableType enum
• Added PoissonianRateCtx metadata for storing rate expressions
• Implemented @poissonians macro with topoissonian helper
• Added ispoissonian and getpoissonianrate accessors

Phase 2: System Storage

• Added poissonians field to System struct
• Two-argument constructor auto-detects poissonians from equations
• Variables/parameters from rate expressions are automatically extracted
• Added poissonians(sys) accessor with hierarchical subsystem support

Phase 3: Compilation to Jumps

• _poissonians_to_jumps extracts coefficients and builds jump affects
• Automatic classification: ConstantRateJump (parameter-only rates) vs VariableRateJump (state/time-dependent rates)
• Same poissonian in multiple equations creates single jump with multiple affects
• Mixed continuous/jump dynamics supported (continuous terms retained)
• Merges with explicitly declared jumps
• Errors on nonlinear poissonian usage (e.g., dN * dN)
• save_positions kwarg flows through mtkcompile to VRJs from poissonians
• User-provided VRJs retain their own save_positions settings
• check_no_poissonians validation ensures mtkcompile is called before JumpProblem

Phase 4: Mathematical Correctness

• Integration tests verify simulation statistics match analytical formulas
• Comparison tests against direct JumpProcesses implementations
• Tests cover pure jump, mixed continuous/jump, and jump-diffusion systems

Example: SIR Model

@parameters β γ
@variables S(t) I(t) R(t)
@poissonians dN_inf(β * S * I) dN_rec(γ * I)

eqs = [
    D(S) ~ -dN_inf,
    D(I) ~ dN_inf - dN_rec,
    D(R) ~ dN_rec
]

@named sir = System(eqs, t)
compiled = mtkcompile(sir)

# Results in 2 VariableRateJumps with proper affects:
# Jump 1: rate = β*S*I, affects = [S ~ Pre(S) - 1, I ~ Pre(I) + 1]
# Jump 2: rate = γ*I, affects = [I ~ Pre(I) - 1, R ~ Pre(R) + 1]

# Create and solve the JumpProblem
jprob = JumpProblem(compiled, [S => 999.0, I => 1.0, R => 0.0, β => 0.1/1000, γ => 0.01], (0.0, 250.0))
sol = solve(jprob, Tsit5())

Example: Jump-Diffusion

@parameters μ σ λ γ
@variables X(t)
@brownians dW
@poissonians dN(λ)

# dX = μX dt + σX dW + γ dN
eqs = [D(X) ~ μ * X + σ * X * dW + γ * dN]

@named sys = System(eqs, t)
compiled = mtkcompile(sys)

# Results in:
# - 1 SDE equation: D(X) ~ μ*X with noise σ*X
# - 1 ConstantRateJump: rate = λ, affect = [X ~ Pre(X) + γ]

[isaacsas]

One question is if we prefer

@poissonians a(rate_expr)

or

@poissonians a [rate = expr]

i.e. explicitly typing the rate metadata. The first seems cleaner to me and more inline with what we write mathematically, but the second has the advantage that we can reuse parse_vars so don't need any new macro parsing code.

[isaacsas]

@AayushSabharwal I think this is good to go assuming the non-broken tests pass. But let me know if you have any feedback. Note that while the amount of code is long, most of it is tests. The MTKBase changes are not too much.

[AayushSabharwal]

Just minor stuff. Thanks for all the code comments!

[isaacsas]

@AayushSabharwal I addressed your comments. I had apparently forgotten to put the new test file in runtests.jl, so I've now added it to InterfaceII right after the jump tests. It passes locally, so hopefully it will now pass on CI. Once that successfully finishes this should be good.

[isaacsas]

Thanks!