fix(pchip): close the flat-data holes in the slope guard (Dual) and its adjoint (Float64)

[omeneghini-pf]

Two related bugs, one trigger: two adjacent control points holding the same value, which
makes a secant exactly zero. Both derivative paths then fall into the degenerate branch they
were meant to avoid. Values are fine — only the derivatives break.

Exact ties are not exotic: clamped or saturated data, quantized/rounded inputs, a plateau in an
otherwise monotone profile, or optimizer variables resting on a shared bound all produce them.

1. Forward: _pchip_harmonic_mean returns Dual(NaN, NaN)

The guard for flat data (both secants zero → 0*0/0, added with the single-division form in
#168) is iszero(den). On a ForwardDiff.Dual, iszero is partial-aware, so a seeded
derivative skips it:

julia> w1 = w2 = 3.0;
julia> δp = Dual{:t}(0.0, 8e-9); δc = Dual{:t}(0.0, 0.0);
julia> den = w1*δc + w2*δp        # Dual(0.0, 2.4e-8)  →  iszero(den) == false
julia> (w1+w2)*δp*δc/den
Dual{:t}(NaN,NaN)                 # NaN in the VALUE, not just the partial

End to end:

julia> x  = collect(0.0:0.1:0.6);
julia> xq = collect(range(0.0, 0.6; length=13));
julia> g(p) = pchip_interp(x, [0.0, 1.0, 2.0 + p[1], 2.0 + p[2], 2.0, 2.0, 5.0];
                           extrap=ExtendExtrap()).(xq);
julia> any(isnan, ForwardDiff.jacobian(g, [0.0, 0.0]))
true

_constant_kernel already handles this correctly (dL_primal = _extract_primal(dL)); this is
the one place the convention was missed.

2. Adjoint: pchip_adjoint returns NaN on plain Float64

Same trigger, no AD required. sign(δ_prev) != sign(δ_curr) is false for 0 vs 0, so the
interior loop takes the active branch and divides by δ²:

julia> y = [0.0, 1.0, 2.0, 2.0, 2.0, 2.0, 5.0];              # plateau between two ramps
julia> pchip_adjoint(x, y, xq; extrap=ExtendExtrap())(ones(13))
7-element Vector{Float64}:
 1.3125, 2.25, NaN, NaN, NaN, NaN, 1.3125

julia> y2 = [0.0, 1.0, 2.0, 2.001, 2.002, 2.003, 5.0];       # same shape, no exact ties
julia> pchip_adjoint(x, y2, xq; extrap=ExtendExtrap())(ones(13))
7-element Vector{Float64}:
 1.3125, 2.25, 1.9375, 2.0, 1.9375, 2.25, 1.3125

(D = w1/0 + w2/0 = Inf, then S*w1/(D²*δ²) = Inf/NaN.) After the fix the tied case returns
that second vector — the no-ties limit — which is one of the added tests.

Why this is easy to miss

• The forward case needs two live seeds. Seed only one of two tied control points and the
  two secant partials cancel on a uniform grid, restoring iszero(den) — so single-variable
  AD/FD spot-checks pass while a full ForwardDiff.jacobian is NaN-poisoned.
• It can surface as a zero derivative rather than a NaN. Any caller that treats a NaN
  evaluation as an invalid point and substitutes a constant fallback hands ForwardDiff a
  partial-free value, so the whole seeded chunk comes back as exactly-zero Jacobian columns —
  a structurally singular system with nothing to warn on. That failure mode is silent, which
  is why it seems worth fixing at the source rather than leaving to callers.

Fix

• Forward: iszero(_extract_primal(den)), matching _constant_kernel.
• Adjoint: skip the node when the forward pass takes its zero-denominator branch, via a shared
  _flat_secants(δp, δc) predicate, so forward and reverse select the same branch. (The
  adjoint builds D = w1/δp + w2/δc rather than den, hence stating the test in terms of the
  secants; on the monotone branch the two are equivalent.) Applied to both the interior loop
  and the periodic kernel; the endpoint blocks never divide by a secant and are unaffected.

At exactly-flat data the Fritsch–Carlson slope is genuinely non-differentiable (positively
homogeneous there), so zero is a subgradient choice — and it is the value the Float64 forward
path already returns, which is what keeps the three paths consistent.

No behavior change for any input whose denominator is nonzero.

Testing

New test/test_pchip_flat_data.jl (13 tests): the kernel guard under a seeded Dual; forward
ForwardDiff.jacobian across a plateau (no NaN, values unchanged, and agreement with finite
differences where the data is smooth); and the adjoint on flat data (finite, equal to the
no-ties limit, and satisfying the dot-product identity ⟨W·y, ȳ⟩ == ⟨y, Wᵀ·ȳ⟩ on uniform and
non-uniform grids).

Full suite on this branch: 70486 pass, 24 broken (pre-existing), 0 failures.

[github-actions]

FastInterpolations.jl Benchmarks

🚨 Regression baseline is a different CPU

	CPU
This PR run	`sapphirerapids
Master baseline (latest commit)	`znver3

The Previous / Imm. Ratio / Grad. Ratio columns compare against master history, but master's most recent commit was benchmarked on a different CPU — so a naive comparison is cross-CPU and unreliable. The ratios below use only same-sapphirerapids|4c history and are blank if master has none on this CPU yet.

All benchmarks (56 total, click to expand)

Benchmark	Current: 60946c0	Previous	Imm. Ratio	Grad. Ratio
10_nd_construct/bicubic_2d	34345 ns	35156.0 ns	0.977	0.902
10_nd_construct/bilinear_2d	767.26 ns	794.3 ns	0.966	0.944
10_nd_construct/tricubic_3d	269582 ns	267086.0 ns	1.009	0.937
10_nd_construct/trilinear_3d	2354.66 ns	2401.0 ns	0.981	0.938
11_nd_eval/bicubic_2d_batch	1303.1 ns	1305.7 ns	0.998	0.93
11_nd_eval/bicubic_2d_scalar	16.89 ns	17.0 ns	0.995	0.904
11_nd_eval/bilinear_2d_scalar	7.9 ns	7.9 ns	0.999	0.878
11_nd_eval/tricubic_3d_batch	2888.5 ns	2896.4 ns	0.997	0.939
11_nd_eval/tricubic_3d_scalar	32.74 ns	32.8 ns	0.998	0.904
11_nd_eval/trilinear_3d_scalar	13.77 ns	13.8 ns	0.997	0.904
12_cubic_eval_gridquery/range_random	3904.88 ns	3908.8 ns	0.999	0.944
12_cubic_eval_gridquery/range_sorted	3897.46 ns	3901.8 ns	0.999	0.938
12_cubic_eval_gridquery/vec_random	7035.5 ns	6246.3 ns	1.126	0.867
12_cubic_eval_gridquery/vec_sorted	2906.26 ns	2909.3 ns	0.999	0.931
13_nd_oneshot_gridquery/bicubic_2d_rand_rand	72214.5 ns	73595.2 ns	0.981	0.923
13_nd_oneshot_gridquery/bicubic_2d_sort_rand	51526.8 ns	50511.4 ns	1.02	0.953
13_nd_oneshot_gridquery/bicubic_2d_sort_sort	49158.6 ns	48382.9 ns	1.016	0.949
13_nd_oneshot_gridquery/bilinear_2d_rand_rand	28417.8 ns	31740.0 ns	0.895	1.074
13_nd_oneshot_gridquery/bilinear_2d_sort_rand	7249.46 ns	8053.3 ns	0.9	0.848
13_nd_oneshot_gridquery/bilinear_2d_sort_sort	4754.34 ns	4755.0 ns	1.0	0.935
14_series_oneshot_batch/constant_inplace_vec_k8_q1000_rand	17117 ns	22362.5 ns	0.765	0.706
14_series_oneshot_batch/linear_inplace_vec_k8_q1000_rand	12959.3 ns	10770.2 ns	1.203	0.955
15_gridded_query/cubic_oneshot_2d_40x32	38792 ns	38839.4 ns	0.999	0.894
15_gridded_query/cubic_persistent_2d_40x32	7073.1 ns	7065.9 ns	1.001	0.93
15_gridded_query/linear_oneshot_2d_40x32	1290 ns	1291.0 ns	0.999	0.926
15_gridded_query/linear_oneshot_3d_12x10x8	1695.34 ns	1694.9 ns	1.0	0.929
15_gridded_query/linear_persistent_2d_40x32	1270.1 ns	1273.7 ns	0.997	0.943
15_gridded_query/linear_persistent_3d_12x10x8	1680.38 ns	1671.9 ns	1.005	0.938
1_cubic_oneshot/q00001	400.54 ns	400.1 ns	1.001	0.902
1_cubic_oneshot/q10000	42049.7 ns	42085.4 ns	0.999	0.932
2_cubic_construct/g0100	1289.3 ns	1292.9 ns	0.997	0.918
2_cubic_construct/g1000	13054.7 ns	13013.8 ns	1.003	0.934
3_cubic_eval/q00001	15.27 ns	18.4 ns	0.83	0.759
3_cubic_eval/q00100	412.58 ns	413.1 ns	0.999	0.925
3_cubic_eval/q10000	39938.1 ns	40006.9 ns	0.998	0.936
4_linear_oneshot/q00001	16.77 ns	16.9 ns	0.993	0.877
4_linear_oneshot/q10000	20457.8 ns	20440.0 ns	1.001	0.936
5_linear_construct/g0100	42.07 ns	44.3 ns	0.95	0.907
5_linear_construct/g1000	347.05 ns	352.3 ns	0.985	0.978
6_linear_eval/q00001	6.67 ns	6.9 ns	0.972	0.878
6_linear_eval/q00100	215.18 ns	215.2 ns	1.0	0.932
6_linear_eval/q10000	20341.3 ns	20389.0 ns	0.998	0.931
7_cubic_range/scalar_query	5.97 ns	5.8 ns	1.029	0.904
7_cubic_vec/scalar_query	8.14 ns	8.3 ns	0.981	0.856
8_cubic_multi/construct_s001_q100	463.56 ns	457.2 ns	1.014	0.915
8_cubic_multi/construct_s010_q100	4078.72 ns	4075.7 ns	1.001	0.937
8_cubic_multi/construct_s100_q100	38662.4 ns	38718.2 ns	0.999	0.944
8_cubic_multi/eval_s001_q100	450.34 ns	452.4 ns	0.995	0.898
8_cubic_multi/eval_s010_q100	1223.22 ns	1226.0 ns	0.998	0.888
8_cubic_multi/eval_s010_q100_scalar_loop	1863.62 ns	1862.2 ns	1.001	0.942
8_cubic_multi/eval_s100_q100	8928.9 ns	9253.3 ns	0.965	0.891
8_cubic_multi/eval_s100_q100_scalar_loop	2981 ns	3029.6 ns	0.984	0.911
9_nd_oneshot/bicubic_2d	33194.9 ns	33695.5 ns	0.985	0.938
9_nd_oneshot/bilinear_2d	505.44 ns	504.7 ns	1.001	0.945
9_nd_oneshot/tricubic_3d	268174.8 ns	265182.1 ns	1.011	0.94
9_nd_oneshot/trilinear_3d	964.3 ns	963.0 ns	1.001	0.911

⚠️ Performance Regression Confirmed ⚠️

After re-running 3 flagged benchmark(s) 10 time(s), 2 regression(s) confirmed.

Benchmark	Current	Previous	Imm. Ratio	Grad. Ratio	Tier
12_cubic_eval_gridquery/vec_random	7035.5 ns	6246.3 ns	1.126	0.867	immediate
14_series_oneshot_batch/linear_inplace_vec_k8_q1000_rand	12959.3 ns	10770.2 ns	1.203	0.955	immediate

Thresholds: immediate > 1.1x (vs latest master), gradual > 1.1x (vs sliding window)

Runner: sapphirerapids|4c — INTEL(R) XEON(R) PLATINUM 8573C, julia 1.12.6. Times are min-merged and compared only against this same machine's history.

This comment was automatically generated by Benchmark workflow.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.31%. Comparing base (071fbdd) to head (60946c0).

Additional details and impacted files

@@           Coverage Diff           @@
##           master     #203   +/-   ##
=======================================
  Coverage   96.31%   96.31%           
=======================================
  Files         164      164           
  Lines       13965    13967    +2     
=======================================
+ Hits        13450    13452    +2     
  Misses        515      515           

Files with missing lines	Coverage Δ
src/pchip/pchip_adjoint.jl	100.00% <100.00%> (ø)
src/pchip/pchip_slopes.jl	93.61% <100.00%> (+0.13%)	⬆️

🚀 New features to boost your workflow:

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