Enable NaN-safe mode by default

[devmotion]

Fixes #857.

This PR enables NaN-safe mode by default. The previous behaviour can still be selected by setting the nansafe_mode preference to false. The documentation section is rewritten accordingly, and the outdated "~5%-10%" performance claim is replaced with the numbers below. Since the CI matrix already sets the preference explicitly to true and false, both modes remain tested.

With this PR, by default:

julia> ForwardDiff.gradient(v -> sqrt(v[1]) + v[2], [0.0, 1.0])
2-element Vector{Float64}:
 Inf
  1.0

julia> log(ForwardDiff.Dual{:tag}(0.0, 0.0))
Dual{:tag}(-Inf,0.0)

(previously [Inf, NaN] and Dual{:tag}(-Inf,NaN)).

I bumped the version to 1.5.0 since the default behaviour changes. Happy to revert that if you'd prefer a patch release.

Timings

Median times (Chairmarks @be, seconds=3). Each mode ran in a fresh process, and the runs alternated twice (false, true, false, true). Both runs agree within a few percent, so the table shows the second run.

Benchmark	NaN-safe disabled	NaN-safe enabled	Ratio
gradient(rosenbrock, rand(1000))	532 μs	770 μs	1.45
hessian(rosenbrock, rand(10))	3.42 μs	5.35 μs	1.56
gradient(x -> sum(exp, x), rand(100))	8.03 μs	8.13 μs	1.01
jacobian(cumprod, rand(10))	624 ns	614 ns	0.98
derivative(x -> sin(x) * exp(x) / x, 1.3)	16.2 ns	16.6 ns	1.03

Overhead shows up in code that does little more than multiply and add (Rosenbrock), where the compare and select are a sizeable fraction of the work. It is within noise once transcendental functions or memory traffic dominate.

Benchmark script

using ForwardDiff, Chairmarks, Random

function rosenbrock(x::AbstractVector)
    a = one(eltype(x))
    b = 100 * a
    result = zero(eltype(x))
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end
sumexp(x) = sum(exp, x)
scalar(x) = sin(x) * exp(x) / x

Random.seed!(1)
x1000 = rand(1000); x100 = rand(100); x10 = rand(10)

@be ForwardDiff.gradient($rosenbrock, $x1000) seconds=3
@be ForwardDiff.gradient($sumexp, $x100) seconds=3
@be ForwardDiff.jacobian($cumprod, $x10) seconds=3
@be ForwardDiff.hessian($rosenbrock, $x10) seconds=3
@be Ref(1.3) ForwardDiff.derivative($scalar, _[]) seconds=3

versioninfo()

Julia Version 1.13.1
Commit 96ca370cf0e (2026-09-25 19:34 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin27.0.0)
  CPU: 10 × Apple M2 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, apple-m2)
Threads: 1 default, 1 interactive, 1 GC (on 10 virtual cores)

Generated code

In NaN-safe mode, each 2-wide vector of partials gets one extra compare and one extra mask instruction (fcmeq + bic on aarch64, fcmp une + select in LLVM IR). The code remains branch-free and fully vectorized.

@code_native debuginfo=:none ForwardDiff._mul_partials(p, p, 2.0, 1.0) with p::Partials{4,Float64}

NaN-safe mode disabled:

	stp	x29, x30, [sp, #-16]!
	mov	x29, sp
	ldp	q2, q3, [x0]
	fmul.2d	v2, v2, v0[0]
	ldp	q4, q5, [x1]
	fmul.2d	v4, v4, v1[0]
	fadd.2d	v2, v2, v4
	fmul.2d	v0, v3, v0[0]
	fmul.2d	v1, v5, v1[0]
	fadd.2d	v0, v0, v1
	stp	q2, q0, [x8]
	ldp	x29, x30, [sp], #16
	ret

NaN-safe mode enabled:

	stp	x29, x30, [sp, #-16]!
	mov	x29, sp
	ldp	q2, q3, [x0]
	fmul.2d	v4, v2, v0[0]
	fcmeq.2d	v2, v2, #0.0
	bic.16b	v2, v4, v2
	ldp	q4, q5, [x1]
	fmul.2d	v6, v4, v1[0]
	fcmeq.2d	v4, v4, #0.0
	bic.16b	v4, v6, v4
	fadd.2d	v2, v2, v4
	fmul.2d	v0, v3, v0[0]
	fcmeq.2d	v3, v3, #0.0
	bic.16b	v0, v0, v3
	fmul.2d	v1, v5, v1[0]
	fcmeq.2d	v3, v5, #0.0
	bic.16b	v1, v1, v3
	fadd.2d	v0, v0, v1
	stp	q2, q0, [x8]
	ldp	x29, x30, [sp], #16
	ret

@code_llvm debuginfo=:none x * x with x::Dual{Nothing,Float64,4} (function body)

NaN-safe mode disabled:

  %1 = load <2 x double>, ptr %"x::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %"y::Dual.unbox", i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  %5 = load <2 x double>, ptr %"y::Dual.partials_ptr", align 8
  %6 = insertelement <2 x double> poison, double %"x::Dual.unbox", i64 0
  %7 = shufflevector <2 x double> %6, <2 x double> poison, <2 x i32> zeroinitializer
  %8 = fmul <2 x double> %7, %5
  %9 = fadd <2 x double> %4, %8
  store <2 x double> %9, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

NaN-safe mode enabled:

  %1 = load <2 x double>, ptr %"x::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %"y::Dual.unbox", i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  %5 = fcmp une <2 x double> %1, zeroinitializer
  %6 = select <2 x i1> %5, <2 x double> %4, <2 x double> zeroinitializer
  %7 = load <2 x double>, ptr %"y::Dual.partials_ptr", align 8
  %8 = insertelement <2 x double> poison, double %"x::Dual.unbox", i64 0
  %9 = shufflevector <2 x double> %8, <2 x double> poison, <2 x i32> zeroinitializer
  %10 = fmul <2 x double> %9, %7
  %11 = fcmp une <2 x double> %7, zeroinitializer
  %12 = select <2 x i1> %11, <2 x double> %10, <2 x double> zeroinitializer
  %13 = fadd <2 x double> %6, %12
  store <2 x double> %13, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

@code_llvm debuginfo=:none exp(x) with x::Dual{Nothing,Float64,4} (function body)

NaN-safe mode disabled:

  %0 = call double @j_exp_2207(double %"d::Dual.unbox")
  ...
  %1 = load <2 x double>, ptr %"d::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %0, i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  store <2 x double> %4, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

NaN-safe mode enabled:

  %0 = call double @j_exp_2244(double %"d::Dual.unbox")
  ...
  %1 = load <2 x double>, ptr %"d::Dual.partials_ptr", align 8
  %2 = insertelement <2 x double> poison, double %0, i64 0
  %3 = shufflevector <2 x double> %2, <2 x double> poison, <2 x i32> zeroinitializer
  %4 = fmul <2 x double> %3, %1
  %5 = fcmp une <2 x double> %1, zeroinitializer
  %6 = select <2 x i1> %5, <2 x double> %4, <2 x double> zeroinitializer
  store <2 x double> %6, ptr %"new::Dual.sroa.2.0.sret_return.sroa_idx", align 8
  ; (same for partials 3 and 4)

Tests

Pkg.test() passes locally in both modes: 9723 tests with the new default and 9641 with nansafe_mode = false.

🤖 Generated with Claude Code

[KristofferC]

Check your benchmarks with #859 perhaps

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 91.23%. Comparing base (a3c0f4f) to head (da4575c).

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #858   +/-   ##
=======================================
  Coverage   91.23%   91.23%           
=======================================
  Files          11       11           
  Lines        1072     1072           
=======================================
  Hits          978      978           
  Misses         94       94           

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[KristofferC]

I don't think we should merge this until something like #859 is in. I didn't fully get the objection to that PR, it was a very long response but feel free to modify as you see fit.

[devmotion]

I benchmarked #859 on more examples and studied a few other alternatives, but none of them was able to consistently outperform this PR here.

Regarding #859 specifically, I think master/this PR is better since

• Conceptually, it implements exactly the underlying math (a zero partial means unperturbed, so its contribution is zero) for every number type, without relying on properties of the factor's type or implementation details of its arithmetic.
• Its lowering is predictable: The per-partial mask is branch-free and SIMD-friendly (fcmeq + bic on aarch64). There's no branch that LLVM could if-convert, speculate, or keep, and no additional code that could push operations over the inlining budget. Its cost is moderate and consistent: In none of the benchmarks above it's much slower than NaN-safe mode disabled (at most 1.59×), whereas the gains of this PR (ie Reduce the cost of NaN-safe mode by checking the scalar factor once #859) come with regressions of up to 3.8× that depend on N, the call site, and LLVM's heuristics.

[devmotion]

As I mentioned in #857, my strong opinion is that the correct mode should be the default, regardless of its performance impact. Unsafe performance optimization must be an opt-in, similar to @fastmath etc. One can always be fast by returning incorrect results - as ForwardDiff does by default. IMO this PR is a bugfix.

Any user can still opt into the faster but unsafe version by setting the preference.

[KristofferC]

regardless of its performance impact.

This is obviously false. You can test by putting a sleep(100) in the operation here and see what the effect will be on the issue tracker.

[KristofferC]

Conceptually,

Who cares about "conceptually", though? The point is to get the result and to get it quickly.

come with regressions of up to 3.8×

This is from benchmarking a completely useless a * x which says absolutely nothing about what the actual performance change will be in real scenarios.

[devmotion]

I disagree with both points.

There's one thing we want to achieve: zero partials out for zero partials in. The approach on master clearly enforces this, and it's easy to check and validate for any reader. The approach in #859 is indirect and does not work in general, hence ::Real fallbacks are needed, division fast path is limited to single-level Duals of IEEEFloats and multiplication fast path to possibly nested Duals of IEEEFloats, the correctness of which are crucially relying on the implementation detail of * and / for Duals. Hence my assessment that the master approach implements the mathematical property directly, without relying on properties of specific number types and implementation details.

Moreover, on my computer I just couldn't see a consistent advantage of #859 in the benchmarks. Neither in microbenchmarks nor in full derivative/gradient calls. Based on my experiments, I just can't reproduce that it is clearly superior.

On the other hand, what my experiments and study of some alternatives revealed was that the compiler much more consistently and reliably vectorizes the approach on master whereas the approach on #859 is compiled in a much less predictable way, even in the alternatives I tried. This seemed to explain also some of the performance regressions of #859.

All in all, in my experiments #859 doesn't seem to be a clear winner but come with a few caveats and disadvantages.

[devmotion]

I put together a larger set of hopefully somewhat realistic comparisons of NaN-safe mode on master (which #858 enables by default) and in #859. I deliberately included cases where I expected #859 to be slower.

All workloads are hand-written without dependencies besides ForwardDiff and use the public API with preallocated configs and default chunk sizes: Rosenbrock; logistic regression, Weibull survival regression, sigmoid Emax PK-PD model, MLP, polynomial and Michaelis-Menten least squares; ODE right-hand side Jacobians (ROBER, 2D Brusselator, 1D Fisher-KPP); gradients through RK4 solves of PK models; Newton's method for Kepler's equation; and the benchmarks from #858 and #859.

Results

Times are relative to master with NaN-safe mode disabled; in the last column, < 1 means #859 is faster than master with NaN-safe mode enabled (with bootstrap 95% CI).

Workload	NaN-safe off, master	NaN-safe on, master	NaN-safe on, #859	#859 / master (NaN-safe on)
Hessians
Rosenbrock, n = 10	2.33 μs	1.91×	0.92×	0.48× (0.46–0.50)
Rosenbrock, n = 20	17.40 μs	2.03×	0.91×	0.45× (0.42–0.46)
logistic regression NLL, 500 × 8	115.3 μs	1.70×	1.73×	1.02× (0.97–1.04)
Weibull regression NLL, 300 × 7	95.90 μs	3.07×	1.34×	0.44× (0.42–0.44)
sigmoid Emax PK-PD NLL, 200 obs	60.27 μs	2.04×	1.30×	0.64× (0.60–0.65)
Gradients, Jacobians, derivatives
gradient: Rosenbrock, n = 1000	552.5 μs	1.45×	1.15×	0.79× (0.75–0.81)
gradient: logistic regression NLL (X * β), 500 × 8	13.42 μs	1.14×	1.04×	0.91× (0.87–0.94)
gradient: logistic regression NLL (loops), 500 × 8	11.62 μs	1.27×	1.51×	1.19× (1.16–1.21)
gradient: Weibull regression NLL, 300 × 7	22.69 μs	1.03×	1.04×	1.01× (0.98–1.02)
gradient: MLP 4-8-1 loss, 49 params, batch 64	43.46 μs	1.18×	1.00×	0.85× (0.79–0.91)
gradient: sigmoid Emax PK-PD NLL, 200 obs	23.52 μs	1.05×	1.03×	0.98× (0.94–1.01)
gradient: quadratic least squares (broadcast), 1000 points	5.57 μs	1.61×	1.53×	0.95× (0.93–0.97)
gradient: Michaelis-Menten least squares (broadcast), 1000 points	5.46 μs	1.62×	2.14×	1.32× (1.28–1.39)
gradient: RK4 PK 1-compartment, 3 params	31.27 μs	1.14×	1.07×	0.93× (0.88–1.01)
gradient: RK4 PK 2-compartment, 5 params	50.38 μs	1.52×	1.31×	0.87× (0.82–0.91)
gradient: sum(exp, x), n = 100	7.62 μs	1.04×	0.86×	0.83× (0.81–0.87)
gradient: sum(xi -> xi^x[1], x), n = 100	40.56 μs	1.05×	1.35×	1.28× (1.24–1.34)
Jacobian: ROBER RHS, 3 states	24.6 ns	1.14×	1.01×	0.88× (0.84–0.92)
Jacobian: 1D Fisher-KPP RHS, 16 states	539.0 ns	1.08×	1.07×	0.99× (0.96–1.02)
Jacobian: 2D Brusselator RHS, 128 states	45.46 μs	1.09×	1.01×	0.92× (0.88–1.02)
Jacobian: cumprod, n = 10	313.4 ns	1.04×	1.01×	0.97× (0.94–1.01)
derivative: Newton for Kepler's equation, 1000 solves	144.4 μs	1.02×	0.98×	0.97× (0.95–0.98)
derivative: sin(x) * exp(x) / x	16.4 ns	1.06×	1.13×	1.06× (1.04–1.11)

With NaN-safe mode disabled, #859 is within 0.98×–1.02× of master everywhere.

Summary

• Hessians: NaN-safe mode on master costs 1.7×–2.0×, and about 3× for the Weibull Hessian (3.0×–3.2× in repeated runs with other data). Reduce the cost of NaN-safe mode by checking the scalar factor once #859 cuts this to 0.44×–0.64× of master, except for logistic regression (1.02×). The overhead on master comes mostly from masking nested duals, not from branch mispredictions.
• First-order code: Reduce the cost of NaN-safe mode by checking the scalar factor once #859 is mostly faster, but slower for logistic regression written with loops (1.19×), Michaelis-Menten least squares (1.32×), the xi^x[1] gradient (1.28×) and the scalar derivative (1.06×). These are not caused by the inlining failure I described in Reduce the cost of NaN-safe mode by checking the scalar factor once #859. The same logistic regression model is 9% faster with Reduce the cost of NaN-safe mode by checking the scalar factor once #859 when written with X * β and broadcasting, and 19% slower when written with loops.
• Exploit symmetry in the Hessian #844 doesn't change this: in a separate run of the Hessians, both the overhead of NaN-safe mode on master and Reduce the cost of NaN-safe mode by checking the scalar factor once #859's advantage are the same with Exploit symmetry in the Hessian #844.

Overall, #859 is not uniformly faster. It is clearly worthwhile for Hessians, but for first-order code the gains are modest, some realistic workloads regress by up to 1.3×, and whether #859 is faster or slower can flip just based on how the same model is written (X * β vs. loops). Since the gains are concentrated in nested duals, one option might be a hybrid that keeps the per-partial mask for first-order partials and uses the factor check only for nested duals; I haven't tested whether that keeps the Hessian gains without the first-order regressions. Independent of performance, I still think the factor check should be restricted to IEEEFloats (and Duals of them), for the reasons in #859.

Setup and caveats

• Apple M2 Pro, Julia 1.13.1, Chairmarks; aarch64 only, x86 is untested.
• Every configuration ran in a fresh process, 10 interleaved rounds in random order. Per process, each workload was sampled with @be for 1 s (0.5 s for the 7 workloads added in a second run: loop logistic regression, Weibull, Fisher-KPP, quadratic and Michaelis-Menten least squares, Kepler); the table shows the median over rounds of the per-process medians. The Exploit symmetry in the Hessian #844 comparison used 8 rounds at 0.5 s.
• Results agree between master and Reduce the cost of NaN-safe mode by checking the scalar factor once #859 up to signed zeros, except for a last-bit difference in Michaelis-Menten: the per-element terms are identical, but the final sum is vectorized differently (as also happens between NaN-safe mode on and off on master).

Benchmark script (suite.jl)

Run once per configuration in a fresh process, with nansafe_mode set in LocalPreferences.toml: julia --project=<env> suite.jl <outfile> 0.5. It writes, per workload, the median and minimum time and a hash of the result.

# Realistic ForwardDiff workloads for comparing NaN-safe mode implementations.
# Usage: julia --project=<env> suite.jl <outfile> [seconds per benchmark]
# Writes one tab-separated line per benchmark: name, median ns, minimum ns, result hash.

using ForwardDiff, DiffResults, LinearAlgebra, Random, Chairmarks

const SECONDS = length(ARGS) >= 2 ? parse(Float64, ARGS[2]) : 0.5

# Hash of a result, normalizing signed zeros (NaN-safe mode implementations may differ in
# the sign of zero partials)
_normalize(x::Real) = iszero(x) ? zero(x) : x
result_hash(x::Real) = hash(_normalize(x))
result_hash(x::AbstractArray) = hash(map(_normalize, x))
result_hash(x::DiffResults.DiffResult) = hash(map(result_hash, (DiffResults.value(x), DiffResults.derivative(x))))
result_hash(x::Tuple) = hash(map(result_hash, x))

const BENCHMARKS = Pair{String,Any}[]
add!(name, f) = push!(BENCHMARKS, name => f)

# Extended Rosenbrock (Optim/CUTEst test problem), as in #858
function rosenbrock(x::AbstractVector)
    a = one(eltype(x))
    b = 100 * a
    result = zero(eltype(x))
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end
let x = rand(MersenneTwister(1), 1000), out = similar(x), cfg = ForwardDiff.GradientConfig(rosenbrock, x)
    add!("gradient: Rosenbrock, n = 1000", () -> ForwardDiff.gradient!(out, rosenbrock, x, cfg))
end
for n in (10, 20)
    let x = rand(MersenneTwister(1), n), out = zeros(n, n), cfg = ForwardDiff.HessianConfig(rosenbrock, x)
        add!("Hessian: Rosenbrock, n = $n", () -> ForwardDiff.hessian!(out, rosenbrock, x, cfg))
    end
end

# Logistic regression negative log-likelihood (GLM fitting, Newton/IRLS, Laplace)
function logreg_nll(β, X, y)
    η = X * β
    return sum(log1p.(exp.(η)) .- y .* η)
end
let rng = MersenneTwister(2), X = randn(rng, 500, 8), y = Float64.(rand(rng, Bool, 500)),
    β = 0.1 * randn(rng, 8), f = β -> logreg_nll(β, X, y),
    gout = similar(β), gcfg = ForwardDiff.GradientConfig(f, β),
    hout = zeros(8, 8), hcfg = ForwardDiff.HessianConfig(f, β)
    add!("gradient: logistic regression NLL (`X * β`), 500 × 8", () -> ForwardDiff.gradient!(gout, f, β, gcfg))
    add!("Hessian: logistic regression NLL, 500 × 8", () -> ForwardDiff.hessian!(hout, f, β, hcfg))
end

# Small MLP (4 → 8 → 1, tanh), mean squared error over a batch of 64 samples,
# gradient w.r.t. all 49 parameters
function mlp_loss(θ, X, y)
    W1 = reshape(view(θ, 1:32), 8, 4)
    b1 = view(θ, 33:40)
    w2 = view(θ, 41:48)
    b2 = θ[49]
    h = tanh.(W1 * X .+ b1)
    ŷ = vec(w2' * h) .+ b2
    return sum(abs2, ŷ .- y) / length(y)
end
let rng = MersenneTwister(3), X = randn(rng, 4, 64), y = randn(rng, 64), θ = 0.3 * randn(rng, 49),
    f = θ -> mlp_loss(θ, X, y), out = similar(θ), cfg = ForwardDiff.GradientConfig(f, θ)
    add!("gradient: MLP 4-8-1 loss, 49 params, batch 64", () -> ForwardDiff.gradient!(out, f, θ, cfg))
end

# Sigmoid Emax (Hill) PK-PD model with power error model, negative log-likelihood
# θ = (E0, Emax, EC50, γ, σ); `C^γ` and `EC50^γ` (`Real^Dual`, `Dual^Dual`), `f^0.8` (`Dual^Real`)
function emax_nll(θ, C, y)
    E0, Emax, EC50, γ, σ = θ
    nll = zero(eltype(θ))
    for i in eachindex(C, y)
        Cγ = C[i]^γ
        f = E0 + Emax * Cγ / (EC50^γ + Cγ)
        sd = σ * f^0.8
        nll += log(sd) + ((y[i] - f) / sd)^2 / 2
    end
    return nll
end
let rng = MersenneTwister(4), C = 10 .^ range(-1, 2; length = 200), θ = [1.0, 10.0, 5.0, 1.5, 0.1],
    y = [θ[1] + θ[2] * c^θ[4] / (θ[3]^θ[4] + c^θ[4]) for c in C] .* (1 .+ 0.1 .* randn(rng, 200)),
    f = θ -> emax_nll(θ, C, y), gout = similar(θ), gcfg = ForwardDiff.GradientConfig(f, θ),
    hout = zeros(5, 5), hcfg = ForwardDiff.HessianConfig(f, θ)
    add!("gradient: sigmoid Emax PK-PD NLL, 200 obs", () -> ForwardDiff.gradient!(gout, f, θ, gcfg))
    add!("Hessian: sigmoid Emax PK-PD NLL, 200 obs", () -> ForwardDiff.hessian!(hout, f, θ, hcfg))
end

# ODE right-hand side Jacobians, as computed by stiff solvers
# ROBER (3 states)
function rober!(du, u)
    k1, k2, k3 = 0.04, 3e7, 1e4
    y1, y2, y3 = u
    du[1] = -k1 * y1 + k3 * y2 * y3
    du[2] = k1 * y1 - k2 * y2^2 - k3 * y2 * y3
    du[3] = k2 * y2^2
    return du
end
let u = [1.0, 1e-4, 0.5], du = similar(u), J = zeros(3, 3), cfg = ForwardDiff.JacobianConfig(rober!, du, u)
    add!("Jacobian: ROBER RHS, 3 states", () -> ForwardDiff.jacobian!(J, rober!, du, u, cfg))
end
# 2D Brusselator (SciML benchmark), 8 × 8 grid, 128 states
const BRUSS_N = 8
function brusselator!(du, u)
    N = BRUSS_N
    A, B, α = 3.4, 1.0, 10.0 * (N - 1)^2
    xs = range(0, 1; length = N)
    U = reshape(view(u, 1:N^2), N, N)
    V = reshape(view(u, N^2+1:2N^2), N, N)
    dU = reshape(view(du, 1:N^2), N, N)
    dV = reshape(view(du, N^2+1:2N^2), N, N)
    lim(k) = k == 0 ? N : (k == N + 1 ? 1 : k)
    @inbounds for j in 1:N, i in 1:N
        ip, im, jp, jm = lim(i + 1), lim(i - 1), lim(j + 1), lim(j - 1)
        u2v = U[i, j]^2 * V[i, j]
        forcing = ((xs[i] - 0.3)^2 + (xs[j] - 0.6)^2 <= 0.01) ? 5.0 : 0.0
        dU[i, j] = α * (U[im, j] + U[ip, j] + U[i, jp] + U[i, jm] - 4U[i, j]) + B + u2v - (A + 1) * U[i, j] + forcing
        dV[i, j] = α * (V[im, j] + V[ip, j] + V[i, jp] + V[i, jm] - 4V[i, j]) + A * U[i, j] - u2v
    end
    return du
end
let rng = MersenneTwister(5), u = rand(rng, 2BRUSS_N^2), du = similar(u), J = zeros(length(u), length(u)),
    cfg = ForwardDiff.JacobianConfig(brusselator!, du, u)
    add!("Jacobian: 2D Brusselator RHS, 128 states", () -> ForwardDiff.jacobian!(J, brusselator!, du, u, cfg))
end

# Forward sensitivities through a fixed-step RK4 solve of PK models,
# sum of squared residuals at 12 observation times
@inline function rk4_step(f, u, p, dt)
    k1 = f(u, p)
    k2 = f(u .+ (dt / 2) .* k1, p)
    k3 = f(u .+ (dt / 2) .* k2, p)
    k4 = f(u .+ dt .* k3, p)
    return u .+ (dt / 6) .* (k1 .+ 2 .* k2 .+ 2 .* k3 .+ k4)
end
# One-compartment model with first-order absorption, p = (ka, CL, V)
pk1cmt(u, p) = (-p[1] * u[1], p[1] * u[1] - p[2] / p[3] * u[2])
# Two-compartment model with first-order absorption, p = (ka, CL, V1, Q, V2)
pk2cmt(u, p) = (-p[1] * u[1],
                p[1] * u[1] - (p[2] + p[4]) / p[3] * u[2] + p[4] / p[5] * u[3],
                p[4] / p[3] * u[2] - p[4] / p[5] * u[3])
function pk_sse(f, p, u0, obs_steps, obs, dt, cmt)
    u = map(x -> x * one(eltype(p)), u0)
    sse = zero(eltype(p))
    step = 0
    for (k, n) in enumerate(obs_steps)
        while step < n
            u = rk4_step(f, u, p, dt)
            step += 1
        end
        sse += (obs[k] - u[cmt] / p[3])^2
    end
    return sse
end
let obs_steps = [5, 10, 20, 40, 60, 80, 120, 160, 240, 320, 400, 480], dt = 0.05
    for (name, f, p, u0) in (("1-compartment, 3 params", pk1cmt, [1.2, 3.0, 20.0], (100.0, 0.0)),
                             ("2-compartment, 5 params", pk2cmt, [1.2, 3.0, 20.0, 2.0, 40.0], (100.0, 0.0, 0.0)))
        obs = [5 * exp(-0.01 * n) for n in obs_steps]
        g = p -> pk_sse(f, p, u0, obs_steps, obs, dt, 2)
        out = similar(p)
        cfg = ForwardDiff.GradientConfig(g, p)
        add!("gradient: RK4 PK $name, 480 steps", () -> ForwardDiff.gradient!(out, g, p, cfg))
    end
end

# The benchmarks from #858 and from the review in #859
let x = rand(MersenneTwister(1), 100), out = similar(x), f = x -> sum(exp, x), cfg = ForwardDiff.GradientConfig(f, x)
    add!("gradient: sum(exp, x), n = 100", () -> ForwardDiff.gradient!(out, f, x, cfg))
end
let x = rand(MersenneTwister(1), 10), J = zeros(10, 10), cfg = ForwardDiff.JacobianConfig(cumprod, x)
    add!("Jacobian: cumprod, n = 10", () -> ForwardDiff.jacobian!(J, cumprod, x, cfg))
end
let f = x -> sin(x) * exp(x) / x, x = Ref(1.3)
    add!("derivative: sin(x) * exp(x) / x", () -> ForwardDiff.derivative(f, x[]))
end
let x = rand(MersenneTwister(1), 100), out = similar(x), f = x -> sum(xi -> xi^x[1], x), cfg = ForwardDiff.GradientConfig(f, x)
    add!("gradient: sum(xi -> xi^x[1], x), n = 100", () -> ForwardDiff.gradient!(out, f, x, cfg))
end

# Workloads exercising the code paths where #859 regressed in microbenchmarks

# Logistic regression NLL written with explicit loops (`Float64 * Dual{8}` per term)
function logreg_nll_loop(β, X, y)
    nll = zero(eltype(β))
    for i in axes(X, 1)
        η = zero(eltype(β))
        for j in axes(X, 2)
            η += X[i, j] * β[j]
        end
        nll += log1p(exp(η)) - y[i] * η
    end
    return nll
end
let rng = MersenneTwister(2), X = randn(rng, 500, 8), y = Float64.(rand(rng, Bool, 500)),
    β = 0.1 * randn(rng, 8), f = β -> logreg_nll_loop(β, X, y), out = similar(β), cfg = ForwardDiff.GradientConfig(f, β)
    add!("gradient: logistic regression NLL (loops), 500 × 8", () -> ForwardDiff.gradient!(out, f, β, cfg))
end

# Weibull survival regression NLL, 7 covariates + shape (8 params);
# `(t/λ)^k` is `Dual^Dual`, the linear predictor is `Float64 * Dual{8}`
function weibull_nll(θ, X, t, δ)
    k = θ[end]
    nll = zero(eltype(θ))
    for i in eachindex(t, δ)
        η = zero(eltype(θ))
        for j in axes(X, 2)
            η += X[i, j] * θ[j]
        end
        z = t[i] / exp(η)
        nll -= δ[i] * (log(k) - η + (k - 1) * log(z)) - z^k
    end
    return nll
end
let rng = MersenneTwister(6), X = randn(rng, 300, 7), t = 0.1 .+ 2 .* rand(rng, 300), δ = Float64.(rand(rng, 300) .< 0.7),
    θ = [0.1 * randn(rng, 7); 1.5], f = θ -> weibull_nll(θ, X, t, δ),
    gout = similar(θ), gcfg = ForwardDiff.GradientConfig(f, θ), hout = zeros(8, 8), hcfg = ForwardDiff.HessianConfig(f, θ)
    add!("gradient: Weibull regression NLL, 300 × 7", () -> ForwardDiff.gradient!(gout, f, θ, gcfg))
    add!("Hessian: Weibull regression NLL, 300 × 7", () -> ForwardDiff.hessian!(hout, f, θ, hcfg))
end

# Jacobian of a 1D reaction-diffusion (Fisher-KPP) RHS, 16 states (chunk size 8)
function fisher_kpp!(du, u)
    n = length(u)
    D, r = 0.1 * (n + 1)^2, 1.0
    for i in 1:n
        left = i == 1 ? zero(u[i]) : u[i-1]
        right = i == n ? zero(u[i]) : u[i+1]
        du[i] = D * (left - 2 * u[i] + right) + r * u[i] * (1 - u[i])
    end
    return du
end
let u = rand(MersenneTwister(7), 16), du = similar(u), J = zeros(16, 16), cfg = ForwardDiff.JacobianConfig(fisher_kpp!, du, u)
    add!("Jacobian: 1D Fisher-KPP RHS, 16 states", () -> ForwardDiff.jacobian!(J, fisher_kpp!, du, u, cfg))
end

# Polynomial least squares with broadcasting, 3 params, 1000 points
poly_sse(θ, t, y) = sum(abs2, y .- (θ[1] .+ θ[2] .* t .+ θ[3] .* t .^ 2))
let rng = MersenneTwister(8), t = range(0, 1; length = 1000), y = 1 .+ 2 .* t .- t .^ 2 .+ 0.1 .* randn(rng, 1000),
    θ = [0.5, 1.0, 0.0], f = θ -> poly_sse(θ, t, y), out = similar(θ), cfg = ForwardDiff.GradientConfig(f, θ)
    add!("gradient: quadratic least squares (broadcast), 1000 points", () -> ForwardDiff.gradient!(out, f, θ, cfg))
end

# Michaelis-Menten least squares with broadcasting, 2 params, 1000 points (division-heavy)
mm_sse(θ, S, v) = sum(abs2, v .- θ[1] .* S ./ (θ[2] .+ S))
let rng = MersenneTwister(9), S = range(0.1, 10; length = 1000), v = 2 .* S ./ (1.5 .+ S) .+ 0.05 .* randn(rng, 1000),
    θ = [1.0, 1.0], f = θ -> mm_sse(θ, S, v), out = similar(θ), cfg = ForwardDiff.GradientConfig(f, θ)
    add!("gradient: Michaelis-Menten least squares (broadcast), 1000 points", () -> ForwardDiff.gradient!(out, f, θ, cfg))
end

# Newton's method for Kepler's equation `E - e sin(E) = M`, derivatives via ForwardDiff
function kepler(M, e)
    g(E) = E - e * sin(E) - M
    E = M
    for _ in 1:6
        E -= g(E) / ForwardDiff.derivative(g, E)
    end
    return E
end
let Ms = range(0, 2π; length = 1000)
    add!("derivative: Newton for Kepler's equation, 1000 solves", () -> sum(M -> kepler(M, 0.3), Ms))
end

#######
# Run #
#######

function run_suite(outfile)
    open(outfile, "w") do io
        for (name, f) in BENCHMARKS
            h = result_hash(f())
            b = @be f() seconds = SECONDS
            println(io, name, '\t', Chairmarks.median(b).time * 1e9, '\t', Chairmarks.minimum(b).time * 1e9, '\t', h)
            flush(io)
        end
    end
end

run_suite(ARGS[1])

Driver (run.jl): environments and interleaved rounds

# Runs suite.jl for every configuration in fresh, interleaved processes.
# Usage: julia run.jl <master checkout> <PR checkout> <rounds> [seconds per benchmark]

using Pkg, Random

const SOURCES = Dict("master" => ARGS[1], "pr" => ARGS[2])
const CONFIGS = [(src, nansafe) for src in ("master", "pr") for nansafe in (false, true)]
config_name((src, nansafe)) = "$(src)-nansafe=$(nansafe)"

function setup_env(cfg)
    src, nansafe = cfg
    dir = joinpath(@__DIR__, "envs", config_name(cfg))
    mkpath(dir)
    write(joinpath(dir, "LocalPreferences.toml"), "[ForwardDiff]\nnansafe_mode = $nansafe\n")
    Pkg.activate(dir) do
        Pkg.develop(path = SOURCES[src])
        Pkg.add(["Chairmarks", "DiffResults"])
        Pkg.precompile()
    end
    return dir
end

function main(rounds, seconds)
    envs = Dict(cfg => setup_env(cfg) for cfg in CONFIGS)
    resdir = joinpath(@__DIR__, "results")
    mkpath(resdir)
    rng = MersenneTwister(42)
    for round in 1:rounds, cfg in shuffle(rng, CONFIGS)
        out = joinpath(resdir, "$(config_name(cfg))-round=$(round).tsv")
        run(`$(Base.julia_cmd()) --startup-file=no --project=$(envs[cfg]) $(joinpath(@__DIR__, "suite.jl")) $out $seconds`)
    end
end

main(parse(Int, ARGS[3]), length(ARGS) >= 4 ? ARGS[4] : "0.5")

[devmotion]

That being said, I want to repeat once more what I already wrote repeatedly: In my opinion this PR is a long overdue bug fix and hence the performance implications are a side effect but not relevant for whether the bug should be fixed or not and this PR should be merged or not. I can't count how often I ran into this correctness issue, in e.g. Turing, Distributions and Pumas. At least I know that I could try turning on the correct mode, in contrast to many users who don't. We're only having this discussion since for whatever reason what should have been a simple bugfix resolving a correctness issue ended up as a secondary non-default opt-in mode.