Fix benchmark bugs in StiffSDE, NonStiffSDE, HybridJumps, Jumps (StochasticHeat, Oval2Timings, BasicSDEWeak, Synapse, EGFR)

[ChrisRackauckas-Claude]

Please ignore until reviewed by @ChrisRackauckas.

This PR fixes benchmark bugs found by the prose audit of StiffDDE, StiffSDE, NonStiffSDE, HybridJumps, and Jumps. There is one commit per page. Each page was re-run from the worktree root with
~/.juliaup/bin/julia +<ver> --threads=auto --project=. benchmark.jl benchmarks/<Folder>/<file>.jmd, where <ver> is the folder Manifest's Julia version. All five runs exit 0.
The machine was heavily loaded during the runs (load average 300+ on 128 cores), so absolute timings are higher and noisier than in the published output.

StiffSDE/StochasticHeat: implicit errors did not decrease with dt

• Root cause: ImplicitRKMil, ImplicitEM, and ImplicitEulerHeun are adaptive by default, so dt= only set the initial step. Every run took about 75 adaptive steps at the default tolerances whatever dt was. A standalone check (100 samples each) confirms it: adaptive ImplicitRKMil takes 75/77/75 steps for dt = 0.1/0.01/0.001, and with adaptive=false it takes 4/26/251.
• Fix: Add adaptive=false to the implicit runs in the "Highest dt" and "Simple Error Analysis" sections. Update the conclusion to state the new errors.
• Evidence (mean(simple_error(...)), 400 samples; Julia 1.11):

method	dt	before (published)	after
ImplicitRKMil (Strat.)	0.1	0.018185517528369	1.0004903139814554
ImplicitRKMil (Strat.)	0.01	0.018418466765605716	0.014410479368842741
ImplicitRKMil (Strat.)	0.001	0.018451613923516377	0.00012493520216962087
ImplicitEulerHeun	0.1	0.009393747399189691	0.8318866546863489
ImplicitEulerHeun	0.01	0.007250951082943943	0.01423630485660903
ImplicitEulerHeun	0.001	0.007449617244075225	0.0001577490469593784
EulerHeun	5e-5	3.294e-9	3.083e-9

StiffSDE/Oval2Timings: the "fixed dt" loops ran adaptively

• Root cause: The SRIW1, RKMil, ImplicitEM, and ImplicitRKMil loops did not pass adaptive=false. At the default tolerances with an initial dt = 2^-16, adaptive SRIW1 ends with ReturnCode.DtLessThanMin on every trajectory (50/50 in a standalone check), which accounts for the 975–984/1000 "fails". RKMil's time was flat in dt. Two more bugs in the same plot:
  • best_adaptive_time = numfails != 0 ? adaptive_time : min(...) overwrote the best time with a failing run's time.
  • The red "stable" markers used hard-coded indices.
• Fix: Add adaptive=false to those four loops. Keep best_adaptive_time unchanged when a run fails. Compute the markers as the largest dt with zero failures. Add an intro and a conclusion (the page had no prose).
• Evidence:

method	before: fails (dt=2^-16..)	after: fails	after: time (s)
SRIW1	984, 975, 978, 975	8, 4, 2, 0	4.0 → 30.1
RKMil	5, 3, 6, 7, 7, 6 (time flat, 0.22–0.31 s)	10, 1, 1, 0, 0, 0	4.4 → 79.4
EM (unchanged code)	10, 1, 1, 0, 0, 0	10, 1, 1, 0, 0, 0	5.4 → 44.6

Before the fix, ImplicitEM and ImplicitRKMil took about 11 s and 20 s at every dt (flat). After it, ImplicitEM takes 41 → 574 s and ImplicitRKMil 43 → 640 s, with no failures. The best adaptive run without failures (SOSRI, abstol 2^-7, reltol 2^-4) takes 1.08 s.

NonStiffSDE/BasicSDEWeakWorkPrecision: "Sample Error: 1000" plotted the 10^4-sample value

• Root cause: sample_size = Int[10; 1e2; 1e3; 1e4], and the line plotted se[end].
• Fix: Plot se[3] in all six plots. Summary: on the linear problem, "several times more time" becomes "more time", because with the corrected line EM reaches it at roughly twice SRIW1's cost.
• Evidence: In the additive problem, the line moves from 4.7557374180683486e-7 (published se[end]) to se[3], which is 4.774260360097123e-6 in this run. The other claims in the summary still hold against the new figures.
• Separately, DiffEqDevTools.get_sample_errors divides by sqrt(N) twice, so its values scale like 1/N instead of 1/sqrt(N). The published numbers show it: 4.36e-6 at N = 1000 vs 4.76e-7 at N = 10^4. That library fix is in progress; this page does not work around it.

HybridJumps/Synapse: swapped axis labels

• Fix: xlabel = "Time", with ylabel = "Voltage" or "N". The re-run figures Synapse_22_1.png and Synapse_23_1.png show time (0–1000) on x.

Jumps/EGFR_Benchmark: y label said "Average" for a median

• Fix: ylabel = "Median Time (s)". The re-run bar chart (NRM 18.4 s, CCNRM 14.1 s, DirectCR 18.0 s, RSSACR 0.85 s) shows the new label.

Not fixed here / not a bug / upstream

• StiffDDE/Oregonator_wpd (flat reference solution): The bug is in DDEProblemLibrary: the sign of k₁ A u₂ in the u1 equation is wrong compared with RADAR5's dr-oregon.f. Upstream fix with test: DDEProblemLibrary: fix sign error in RADAR5 Oregonator (u1 equation) DiffEqProblemLibrary.jl#233. The page change (tolerances scaled to the solution) is a separate draft PR, blocked on that release: StiffDDE/Oregonator_wpd: tolerances for the corrected (oscillating) RADAR5 Oregonator #1905.
• StiffDDE/QuorumSensing (one Tsit5 point at ~1e10 in the L∞ RK diagram): not a benchmark bug, no change. The point is MethodOfSteps(Tsit5()) at abstol=1e-4, reltol=1e-1. It returns ReturnCode.Success with u3 oscillating between ±7.1e10 from t ≈ 36 (72 accepted / 28 rejected steps). The true u3 is ~1e-7, far below abstol, and with reltol=0.1 the error estimator accepts a growing unstable mode. reltol=5e-2 gives ReturnCode.Unstable (NaN'd by DiffEqDevTools), and reltol=1e-2 gives max|u3| = 8.9e-5. With MethodOfSteps(Tsit5(); constrained=true), which keeps steps no longer than the lag τ = 2, max|u3| stays at 1.3e-4 even at reltol=1e-1. So this is a real result: an explicit method on a stiff DDE at a very loose tolerance.
• HybridJumps/MultivariateHawkes (VR_Direct (recursive) 9.987 ms at V=30 vs 3.864 s at V=40): under investigation in JumpProcesses; this PR does not change the page.

🤖 Generated with Claude Code (model: claude-opus-5-5[1m])

https://claude.ai/code/session_01Vrsu4PESdABpNaxYiTBVfd

---

Independent review (Devin CLI / fusion-claude-opus-5-5-high-sidekick-swe-2-medium)

VERDICT: CHANGES

Risk assessment

• Risk: low. The PR changes five benchmark .jmd pages and nothing else: no tolerances loosened, no solvers dropped, no dependency, Project or Manifest changes, and no added code comments.
• Blast radius: the published pages StiffSDE/StochasticHeat, StiffSDE/Oval2Timings, NonStiffSDE/BasicSDEWeakWorkPrecision, HybridJumps/Synapse and Jumps/EGFR_Benchmark.
• Evidence: one of the five pages has a CI weave that passed (Oval2Timings). I re-ran the StochasticHeat numerics locally. The other three changes only touch labels or an index.
• Independent review: this review.
• Merge: needs review. The PR is a draft marked "ignore until reviewed by @ChrisRackauckas".

Full review

VERDICT: CHANGES

Review of #1906 (head 466b4ea, content-identical to fb9edca plus a master merge).

Risk assessment

• Risk: low. The PR changes five benchmark .jmd pages and nothing else: no tolerances loosened, no solvers dropped, no dependency, Project or Manifest changes, and no added code comments.
• Blast radius: the published pages StiffSDE/StochasticHeat, StiffSDE/Oval2Timings, NonStiffSDE/BasicSDEWeakWorkPrecision, HybridJumps/Synapse and Jumps/EGFR_Benchmark.
• Evidence: one of the five pages has a CI weave that passed (Oval2Timings). I re-ran the StochasticHeat numerics locally. The other three changes only touch labels or an index.
• Independent review: this review.
• Merge: needs review. The PR is a draft marked "ignore until reviewed by @ChrisRackauckas".

Blocking findings

1. The new Oval2Timings conclusion states timings that contradict the page CI produces. Confirmed from the CI artifact. The conclusion copies absolute timings from the author's machine, which the PR body itself describes as running at "load average 300+ on 128 cores". The CI weave on the benchmark runner (job 107296517719, artifact benchmark-benchmarks-StiffSDE-Oval2Timings-jmd) prints these numbers, and they are what the published page will show:

   Claim in the new ## Conclusion	CI weave output
   EM, RKMil and SRIW1 at Δt = 2^-19: "about 18 s, 19 s, and 30 s"	3.86 s, 4.26 s, 8.21 s
   Implicit methods "between about 40 s and 640 s"	ImplicitEM 11.0–169.7 s, ImplicitRKMil 11.1–165.9 s
   SOSRI (2^-7, 2^-4) "about 1.1 s"	0.124 s
   "roughly 16 times faster than the cheapest fixed-step run without failures"	3.86 / 0.124 ≈ 31×

   The qualitative claims all hold in CI: all three explicit methods first have zero failures at Δt = 2^-19 (EM and RKMil fail 10/1/1/0, SRIW1 fails 8/4/2/0), the implicit methods have no failures, and SOSRI(2^-7, 2^-4) is the fastest adaptive run. Only the magnitudes are wrong. The ratio is off by 2×, so the shift is not a uniform slowdown. Fix: either use the CI numbers, or state results in terms that survive machine load (Δt thresholds, failure counts, "more than an order of magnitude faster"). A page whose purpose is to correct prose should not publish seconds that disagree with the output printed directly above them.

Non-blocking findings

1. CI has not woven 4 of the 5 touched pages. Confirmed via gh. In run 36085775894 (current head) all five Run: jobs are queued. In the previous run, 35894793886, the StochasticHeat, Synapse, EGFR and BasicSDEWeakWorkPrecision jobs were cancelled after 24 h without ever getting a runner (runner_name empty, no steps). So the claim that "all five runs exit 0" rests only on the author's local runs, and the PR body pastes no logs for them. The risk is small: Synapse and EGFR change labels only, and BasicSDEWeak changes se[end] to se[3] on a length-4 vector. Still, let the queued run finish before merging, at least for StochasticHeat and BasicSDEWeak.
2. BasicSDEWeakWorkPrecision: the index fix is correct, but the threshold it plots is still wrong upstream. Confirmed by reading the code. In DiffEqDevTools src/benchmark.jl (get_sample_errors, around lines 1041–1049), mean_solution_ends is already a mean over numruns[i] samples, and its std is then divided by sqrt(numruns[i]) again. The PR body's claim that the values scale like 1/N is therefore correct, and the "Sample Error: 1000" line sits about √1000 ≈ 32× below the page's own definition E_S = sqrt(V/N). The PR still rewrites the Summary ("several times more time" → "more time"), and "EM and RKMil do not reach the sample error" is judged against that too-low line. Consider leaving the Summary alone until the upstream fix lands, or noting the caveat. The body says the library fix is "in progress", but I found no open DiffEqDevTools PR for it.
3. StochasticHeat: the new numbers make the surviving prose weaker. At dt = 0.1 the fixed-step squared error (~0.96–0.99) is larger than the squared norm of the true solution (0.649 at t = 0.25 with W = 0). "Do not have a stepsize limit" is true only in the stability sense. The sentence "unless we want an error of around 10^-10 we are better off using an implicit method" predates this PR, and no timing on the page supports it. This is out of scope for a bug-fix PR, but worth a follow-up by the prose audit.
4. The Oval2 "after" table in the PR body (for example SRIW1 "4.0 → 30.1" s, EM "5.4 → 44.6" s) also differs from CI (1.04 → 8.21 s, 0.51 → 14.7 s). This is body-only, so it matters less. The failure counts in the body match CI exactly.

Checked and fine:

• adaptive=false on the Oval2 fixed-step loops.
• best_adaptive_time now ignores failing runs. The first chunk seeds it with Inf.
• findfirst(==(0), fails[:,i]) with fails initialised to -1 picks the largest tested Δt with zero failures. The plotted markers in the CI figure land at 2^-19 for all three methods, consistent with the failure counts.
• Synapse axes: res.t is on x.
• EGFR: the bars are median(bm).time, matching the new "Median Time" label.
• One commit per page, no scope creep.

What I ran

• gh pr view/diff/checks 1906, gh run view for runs 35894793886 and 36085775894, gh api .../actions/jobs/<id> for the runner and steps of the cancelled jobs.
• Downloaded the Oval2Timings CI artifact and read the woven Oval2Timings.md output and Oval2Timings_19_1.png. Source of the blocking-finding numbers: confirmed from CI output.
• Cloned the PR head read-only to /tmp/rev1906/repo. Read DiffEqDevTools get_sample_errors at master 70a3759.
• Confirmed by running: a reduced StochasticHeat reproduction (/tmp/rev1906/stochheat_check.jl, 100 samples, julia 1.11.9, --project=benchmarks/StiffSDE instantiated from the PR's Manifest, exit 0). With adaptive=true the error is flat in dt: ImplicitRKMil 0.0184 / 0.0209 / 0.0209 with 75–77 points, ImplicitEulerHeun 0.0073 / 0.0053 / 0.0055 with 105–112 points. With adaptive=false, for dt = 0.1 / 0.01 / 0.001: ImplicitRKMil 0.991 / 0.0161 / 1.33e-4 (4 / 26 / 251 points) and ImplicitEulerHeun 0.963 / 0.0144 / 1.51e-4. EulerHeun at dt = 5e-5 gives 2.99e-9. These match the PR table and the new conclusion ("about 1 → 10^-2 → 10^-4", "about 3×10^-9") within sampling noise.

What I did not verify

• I did not run full weaves of any page. The StochasticHeat "Highest dt" @time chunks, BasicSDEWeakWorkPrecision (including the claim that EM reaches the line at "roughly twice SRIW1's cost"), Synapse and EGFR are unverified beyond reading the code.
• The rendered Synapse and EGFR figures and the BasicSDEWeak plots: unverified.
• The PR body's claims about StiffDDE/Oregonator, QuorumSensing and MultivariateHawkes: this PR does not touch those pages, and I did not check them.

Links:

• PR: Fix benchmark bugs in StiffSDE, NonStiffSDE, HybridJumps, Jumps (StochasticHeat, Oval2Timings, BasicSDEWeak, Synapse, EGFR) #1906
• Oval2 CI weave (passed): https://github.com/SciML/SciMLBenchmarks.jl/actions/runs/35894793886/job/107296517719
• Previous benchmark run (4 jobs cancelled, never started): https://github.com/SciML/SciMLBenchmarks.jl/actions/runs/35894793886
• Current benchmark run (queued): https://github.com/SciML/SciMLBenchmarks.jl/actions/runs/36085775894
• DiffEqDevTools get_sample_errors: https://github.com/SciML/DiffEqDevTools.jl/blob/70a3759be01c0c1eaea0384a8a3acb071db8b30a/src/benchmark.jl#L977-L1051

---

🤖 Posted by an AI agent — orchestrator: Devin Sidekick on behalf of @ChrisRackauckas · reviewer: Devin CLI / fusion-claude-opus-5-5-high-sidekick-swe-2-medium
Conversation: local session; review transcript ~/sandbox/fleet-master/jobs/909/log-SciMLBenchmarks.jl-1906.txt on Chris's Mac

[pankgeorg]

VERDICT: MERGE

No blocking findings in the changed code.

• benchmarks/StiffSDE/Oval2Timings.jmd:43 and benchmarks/StiffSDE/Oval2Timings.jmd:142: failed adaptive runs now leave the best successful time unchanged; the intended fixed-step runs explicitly disable adaptivity. At benchmarks/StiffSDE/Oval2Timings.jmd:190, dts is descending and findfirst(==(0), fails[:,i]) selects the largest tested step with zero failures, while untested -1 entries cannot qualify.
• benchmarks/NonStiffSDE/BasicSDEWeakWorkPrecision.jmd:76: all six markers now use the third entry of [10,100,1000,10000], matching the 1,000-trajectory error experiments. This is an indexing correction; it does not fix or validate the separately reported get_sample_errors normalization issue.
• benchmarks/StiffSDE/StochasticHeat.jmd:62: the implicit comparisons now actually vary the fixed step. benchmarks/HybridJumps/Synapse.jmd:2451 and benchmarks/Jumps/EGFR_Benchmark.jmd:63 correctly label the plotted axes and median statistic.

Validation: inspected the surrounding loops and sampling setup. The Oval2 CI weave passed. I did not reproduce the stochastic timing numbers or run benchmarks locally; the other cancelled weaves do not establish successful execution.

Reviewed head: fb9edcaba7bd3fea56d35c9273f5a46b754198c5.

CI comparison (checked 2026-09-24; each row is a failed PR job):

Job	PR	Default branch
Runic / Runic Format Check	FAILURE	SUCCESS — latest run

All these PR Runic logs report formatting in unchanged cutest_benchmark_utils.jl. That failure also occurred on the September 23 base run, but the latest default-branch Runic run is green; it should not be described as currently red on default.

Cancelled checks are incomplete validation, not passing or red test results:

Job	PR	Latest matching default job
Run: benchmarks/StiffSDE/StochasticHeat.jmd	CANCELLED	CANCELLED; absent from newest workflow run
Run: benchmarks/HybridJumps/Synapse.jmd	CANCELLED	CANCELLED; absent from newest workflow run
Run: benchmarks/Jumps/EGFR_Benchmark.jmd	CANCELLED	CANCELLED; absent from newest workflow run
Run: benchmarks/NonStiffSDE/BasicSDEWeakWorkPrecision.jmd	CANCELLED	CANCELLED; absent from newest workflow run

🤖 Posted for @pankgeorg by an AI agent — harness: Codex CLI 0.155.0 · model: gpt-6-astra
Session: local Codex session 01a0d371-cd25-7b23-9706-a1a1fa644f9f