Seeded Poissonians SIR mean assertion is process-nondeterministic

[ChrisRackauckas-Claude]

Summary

The seeded SIR comparison in the Poissonians tests can pass or fail in fresh Julia processes on the same checkout and released dependency graph. The observed failure was a 2.747-standard-error fluctuation at one of 250 tested time points. Repeated runs also produced different symbolic and direct path arrays despite the explicit StableRNG(12345) and per-solve seeds.

I did not change or loosen the assertion. I also do not claim this is a regression from the AutoDespecialize work: clean master passes, the PR head both passes and fails across runs, and even two nominally identical direct JumpProcesses controls diverge.

Official test reproduction

The test is the SIR model: compare symbolic @poissonians to direct JumpProcesses testset. It performs 2,000 symbolic and 2,000 direct simulations and compares their means independently at 250 time points with a 10% pointwise threshold.

Clean master:

git fetch upstream
git worktree add --detach ../mtk-poisson-master 3061ef5b1d422d972cace7af58913c55d4a7f3ed
cd ../mtk-poisson-master
JULIA_PKG_SERVER_REGISTRY_PREFERENCE=eager GROUP=InterfaceII \
  julia +1.13 --project=lib/ModelingToolkitBase \
  -e 'using Pkg; Pkg.test(coverage=false)'

Observed locally:

Poissonians Test | 123  123 | 2m15.4s
Testing ModelingToolkitBase tests passed

PR-head checkout used by the failing CI job:

git fetch upstream pull/4919/head:refs/remotes/upstream/pr-4919
git worktree add --detach ../mtk-poisson-pr 0596aa9e30ed7db1370fa4790439e564dac9dc23
cd ../mtk-poisson-pr
JULIA_PKG_SERVER_REGISTRY_PREFERENCE=eager GROUP=InterfaceII \
  julia +1.13 --project=lib/ModelingToolkitBase \
  -e 'using Pkg; Pkg.test(coverage=false)'

The released-graph rerun of that head passed with 1,040 passes and 7 pre-existing broken tests. The linked Julia-pre CI job failed the SIR assertion on the same commit. Repeating the command in fresh processes is therefore the direct probabilistic reproducer; it is not deterministic on every invocation.

Exact path evidence

For diagnosis I copied the official test's relevant order into a recorder: 2,000 pure-Poisson solves with seeds 1111:3110, 1,000 birth-death solves with seeds 2222:3221, then the 2,000 symbolic and direct SIR paths with seeds 3333:5332. Each SIR path was sampled at 1.0:250.0 and the complete 2000 × 250 arrays were SHA-256 hashed.

Checkout/process	Symbolic SHA-256	Direct SHA-256	Result
master, run 1	0a9cc4e03f052ad7a1494baa10e7f53d674f5f9a6bc3d9adc48b47ba9071b829	b0e7b849ac295ef2f5eb11b0ea6e4cdba41528451da09f6edd5657174a6ff319	pass; max metric 0.04650
master, run 2	c09840ec3ebea26174bb92302bd91d74b0d110549b7178a9f8daebf96d00e631	2e60a96c18846e9362ebe199fe072400cbe85ab7d809a5bd7b515a3667492c4f	pass; max metric 0.04657
PR head, run 1	315277003b4680563dc4a68bb41d5595fdf4add34379870af6d8503a6b680ec4	938f1338c24d1ecf6fabeaa1f4f7729d638dc002171b9a943d0d9aa4a6f6b322	pass; max metric 0.06450
PR head, run 2	ff9e69384b5ff97bb083f58e6245d73114fe8e6657d8f3d4acd36a86af820139	8dbc9c31658a4c1f010d3f5f57dec1e746bfe6cd7b06be39f3657e0d90be23c2	fail at time 29; metric 0.10287

Between the two clean-master processes, the first symbolic divergence was row 15 / seed 3347 / time 1 (0 versus 1), and the first direct divergence was row 130 / seed 3462 / time 1 (0 versus 1). Between the two PR-head processes, the first symbolic divergence was row 74 / seed 3406 / time 1 (1 versus 0), and the first direct divergence was row 68 / seed 3400 / time 1 (0 versus 1). The numerical SIR parameters were identical in all cases: (β, γ) = (0.0001, 0.01).

At the failing PR-head point, time 29:

symbolic mean = 1.501
direct mean   = 1.361
difference    = 0.140
symbolic variance = 2.835416708354179
direct variance   = 2.358858429214608
SE(difference) = sqrt((2.835416708354179 + 2.358858429214608) / 2000)
               = 0.05096211895893256
z = 0.140 / 0.05096211895893256
  = 2.747138518961856

The assertion tests 250 correlated time points and fails if any point crosses the threshold. This evidence makes the CI failure consistent with a statistically brittle assertion, while the differing full-array hashes show that the explicit seeds do not currently make these paths process-reproducible.

Separate JumpProblem specialization observation

This does not explain the stochastic assertion failure, but it is a separate coverage gap exposed while checking the feature. On one compiled MTK system, I constructed a Float64-parameter JumpProblem, replaced only the tunable parameter through SciMLStructures.replace(Tunable(), ..., Float32[0.1]), and used supported remake to create the second problem. The results were:

parameter_types_equal=false
outer_jump_problem_types_equal=false
inner_problem_types_equal=false
inner_function_specialization=(SciMLBase.AutoSpecialize, SciMLBase.AutoSpecialize)
inner_function_types_equal=true
solved_problem_types_equal=false
solved_function_types_equal=false

Thus the same symbolic function has a stable function type before solving, but the outer JumpProblem, its concrete inner problem, and the solved problem/function types still vary with the concrete MTK parameter layout. MTK builds the underlying problem and then delegates to the JumpProcesses constructor; JumpProcesses stores the concrete inner problem type in JumpProblem and creates a new concrete problem in its constructor. This should be addressed separately from the stochastic test.

No code PR is attached because I do not have a deterministic fail-before test, and I did not loosen the existing threshold.

Environment

Julia Version 1.13.0-rc2
Commit 4f599a0fa94 (2026-08-07 08:32 UTC)
Linux x86_64; 128 × AMD EPYC 7502; LLVM 20.1.8; 1 Julia thread

ModelingToolkitBase 1.64.0 (local checkout at the hashes above)
SciMLBase 3.46.0
DiffEqBase 7.14.0
OrdinaryDiffEq 7.6.0
OrdinaryDiffEqCore 4.14.2
JumpProcesses 9.29.2
StableRNGs 1.0.4

Links

• Failing Julia-pre CI job: https://github.com/SciML/ModelingToolkit.jl/actions/runs/31519448170/job/93874566648
• Test source on clean master:

  ModelingToolkit.jl/lib/ModelingToolkitBase/test/poissonians.jl

  Lines 760 to 846 in 3061ef5

  	# Should match analytical steady state
  	@test abs(mean_sym - E_X_ss) < 0.1 * E_X_ss
  	end
  	@testset "SIR model: compare symbolic @poissonians to direct JumpProcesses" begin
  	# Classic SIR: S + I --(β*S*I)--> 2I, I --(γ*I)--> R
  	# State-dependent rates → VariableRateJumps → needs ODE solver
  	@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 sys = System(eqs, t)
  	compiled_sys = mtkcompile(sys)
  	# Both rates depend on unknowns → both are VariableRateJumps
  	sys_jumps = ModelingToolkitBase.jumps(compiled_sys)
  	@test length(sys_jumps) == 2
  	@test all(j -> j isa VariableRateJump, sys_jumps)
  	β_val = 0.1 / 1000
  	γ_val = 0.01
  	S0, I0, R0 = 999.0, 1.0, 0.0
  	T = 250.0
  	Nsims = 2000
  	save_times = 1.0:1.0:T
  	Ntimes = length(save_times)
  	jprob_sym = JumpProblem(
  	compiled_sys,
  	[S => S0, I => I0, R => R0, β => β_val, γ => γ_val], (0.0, T);
  	save_positions = (false, false), rng
  	)
  	seed = 3333
  	R_sym = zeros(Nsims, Ntimes)
  	for i in 1:Nsims
  	sol = solve(jprob_sym, Tsit5(); seed)
  	for (k, t_) in enumerate(save_times)
  	R_sym[i, k] = sol(t_; idxs = R)
  	end
  	seed += 1
  	end
  	# Build directly with JumpProcesses using VariableRateJumps for fair comparison
  	f_direct(du, u, p, t) = (du .= 0) # No continuous dynamics
  	oprob_direct = ODEProblem(f_direct, [S0, I0, R0], (0.0, T), (β_val, γ_val))
  	r1(u, p, t) = p[1] * u[1] * u[2]
  	function a1!(integ)
  	integ.u[1] -= 1
  	integ.u[2] += 1
  	end
  	j1 = VariableRateJump(r1, a1!)
  	r2(u, p, t) = p[2] * u[2]
  	function a2!(integ)
  	integ.u[2] -= 1
  	integ.u[3] += 1
  	end
  	j2 = VariableRateJump(r2, a2!)
  	jprob_direct = JumpProblem(
  	oprob_direct, Direct(), j1, j2;
  	rng, save_positions = (false, false)
  	)
  	seed = 3333
  	R_direct = zeros(Nsims, Ntimes)
  	for i in 1:Nsims
  	sol = solve(jprob_direct, Tsit5(); seed)
  	for (k, t_) in enumerate(save_times)
  	R_direct[i, k] = sol(t_)[3]
  	end
  	seed += 1
  	end
  	# Compare means at every saved time point
  	mean_sym = vec(mean(R_sym; dims = 1))
  	mean_direct = vec(mean(R_direct; dims = 1))
  	# All time points should match within tolerance
  	@test all(
  	abs(mean_sym[k] - mean_direct[k]) < 0.1 * max(mean_direct[k], 1.0)
  	for k in 1:Ntimes
  	)

• PR-head MTK JumpProblem construction:

  ModelingToolkit.jl/lib/ModelingToolkitBase/src/problems/jumpproblem.jl

  Lines 2 to 138 in 0596aa9

  	@fallback_iip_specialize function JumpProcesses.JumpProblem{iip, spec}(
  	sys::System, op, tspan::Union{Tuple, Nothing};
  	check_compatibility = true, eval_expression = false, eval_module = @__MODULE__,
  	checkbounds = false, aggregator = JumpProcesses.NullAggregator(),
  	callback = nothing, rng = nothing, save_positions = (true, true), kwargs...
  	) where {iip, spec}
  	check_complete(sys, JumpProblem)
  	check_compatibility && check_compatible_system(JumpProblem, sys)
  	if haskey(kwargs, :tstops)
  	throw(
  	ArgumentError(
  	"Passing `tstops` directly to `JumpProblem(::System, ...)` is not supported. " *
  	"Define tstops on the `System` via the `tstops` keyword instead."
  	)
  	)
  	end
  	has_vrjs = any(x -> x isa VariableRateJump, jumps(sys))
  	has_eqs = !isempty(equations(sys))
  	has_noise = get_noise_eqs(sys) !== nothing || !isempty(brownians(sys))
  	if (has_vrjs || has_eqs)
  	if has_eqs && has_noise
  	prob = SDEProblem{iip, spec}(
  	sys, op, tspan; check_compatibility = false,
  	build_initializeprob = false, checkbounds, check_length = false,
  	_skip_events = true, _skip_tstops = true, kwargs...
  	)
  	elseif has_eqs
  	prob = ODEProblem{iip, spec}(
  	sys, op, tspan; check_compatibility = false,
  	build_initializeprob = false, checkbounds, check_length = false,
  	_skip_events = true, _skip_tstops = true, kwargs...
  	)
  	else
  	_, u0,
  	p = process_SciMLProblem(
  	EmptySciMLFunction{iip}, sys, op;
  	t = tspan === nothing ? nothing : tspan[1],
  	check_length = false, build_initializeprob = false, kwargs...
  	)
  	observedfun = ObservedFunctionCache(
  	sys; eval_expression, eval_module,
  	checkbounds
  	)
  	f = (du, u, p, t) -> (du .= 0; nothing)
  	df = ODEFunction{true, spec}(f; sys, observed = observedfun)
  	prob = ODEProblem{true}(df, u0, tspan, p; kwargs...)
  	end
  	else
  	_f, u0,
  	p = process_SciMLProblem(
  	EmptySciMLFunction{iip}, sys, op;
  	t = tspan === nothing ? nothing : tspan[1], check_length = false, build_initializeprob = false, kwargs...
  	)
  	f = DiffEqBase.DISCRETE_INPLACE_DEFAULT
  	observedfun = ObservedFunctionCache(
  	sys; eval_expression, eval_module, checkbounds
  	)
  	df = DiscreteFunction{true, true}(
  	f; sys = sys, observed = observedfun,
  	initialization_data = get(_f.kwargs, :initialization_data, nothing)
  	)
  	prob = DiscreteProblem(df, u0, tspan, p; kwargs...)
  	end
  	# Create SymbolicTstops for all paths and forward via JumpProblem kwargs.
  	# Inner problems (SDEProblem/ODEProblem) are created with _skip_tstops = true
  	# to avoid duplication.
  	tstops = SymbolicTstops(
  	sys, GeneratedFunctionOptions(; expression = Val{false}, eval_expression, eval_module)
  	)
  	dvs = unknowns(sys)
  	unknowntoid = Dict(value(unknown) => i for (i, unknown) in enumerate(dvs))
  	js = jumps(sys)
  	invttype = prob.tspan[1] === nothing ? Float64 : typeof(1 / prob.tspan[2])
  	# handling parameter substitution and empty param vecs
  	p = (prob.p isa DiffEqBase.NullParameters || prob.p === nothing) ? Num[] : prob.p
  	majpmapper = JumpSysMajParamMapper(sys, p; jseqs = js, rateconsttype = invttype)
  	_majs = Vector{MassActionJump}(filter(x -> x isa MassActionJump, js))
  	_crjs = Vector{ConstantRateJump}(filter(x -> x isa ConstantRateJump, js))
  	vrjs = Vector{VariableRateJump}(filter(x -> x isa VariableRateJump, js))
  	majs = isempty(_majs) ? nothing : assemble_maj(_majs, unknowntoid, majpmapper)
  	crjs = ConstantRateJump[
  	assemble_crj(sys, j, unknowntoid; eval_expression, eval_module)
  	for j in _crjs
  	]
  	vrjs = VariableRateJump[
  	assemble_vrj(sys, j, unknowntoid; eval_expression, eval_module)
  	for j in vrjs
  	]
  	jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs)
  	# dep graphs are only for constant rate jumps
  	nonvrjs = ArrayPartition(_majs, _crjs)
  	if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator) ||
  	(aggregator isa JumpProcesses.NullAggregator)
  	jdeps = asgraph(sys; eqs = nonvrjs)
  	vdeps = variable_dependencies(sys; eqs = nonvrjs)
  	vtoj = jdeps.badjlist
  	jtov = vdeps.badjlist
  	jtoj = needs_depgraph(aggregator) ? eqeq_dependencies(jdeps, vdeps).fadjlist :
  	nothing
  	else
  	vtoj = nothing
  	jtov = nothing
  	jtoj = nothing
  	end
  	# handle events, making sure to reset aggregators in the generated affect functions
  	# preprocess op to convert Symbol keys to Symbolic using main system before passing
  	# to process_events (which may create ImplicitDiscreteProblems for affect subsystems)
  	op_processed = operating_point_preprocess(sys, op)
  	cbs = process_events(
  	sys; callback, eval_expression, eval_module, op = op_processed, reset_jumps = true,
  	tspan
  	)
  	if rng !== nothing
  	kwargs = (; kwargs..., rng)
  	end
  	if tstops !== nothing
  	kwargs = (; kwargs..., tstops)
  	end
  	# MTK requires pre-scaled rate expressions; never ask JumpProcesses to rescale.
  	return JumpProblem(
  	prob, aggregator, jset; dep_graph = jtoj, vartojumps_map = vtoj,
  	jumptovars_map = jtov, scale_rates = false, nocopy = true,
  	callback = cbs, save_positions, kwargs...
  	)
  	end

• JumpProcesses 9.29.2 problem construction: https://github.com/SciML/JumpProcesses.jl/blob/5dc1d35a6277c4e2d7eaf4727faa2378ea87e7d2/src/problem.jl#L73-L317
• JumpProcesses solve path: https://github.com/SciML/JumpProcesses.jl/blob/5dc1d35a6277c4e2d7eaf4727faa2378ea87e7d2/src/solve.jl#L1-L65