Fix data race in BSpline and QuadraticSpline evaluation

src/derivatives.jl

@@ -256,7 +256,8 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
     n = length(A.t)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), n)
     spline_coefficients!(sc, A.d - 1, A.k, t_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
@@ -280,7 +281,8 @@ function _derivative(
     n = length(A.t)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), n)
     spline_coefficients!(sc, A.d - 1, A.k, t_)
     ducum = zeros(size(A.u)[1:(end - 1)]...)
     if t == A.t[1]
@@ -302,7 +304,8 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     idx = get_idx(A, t, iguess)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), A.h)
     spline_coefficients!(sc, A.d - 1, A.k, t_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
@@ -325,7 +328,8 @@ function _derivative(
     idx = get_idx(A, t, iguess)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), A.h)
     spline_coefficients!(sc, A.d - 1, A.k, t_)
     ducum = zeros(size(A.u)[1:(end - 1)]...)
     if t == A.t[1]

src/interpolation_caches.jl

@@ -646,7 +646,7 @@ function QuadraticSpline(
     k, A = quadratic_spline_params(t, sc)
     c = A \ u
 
-    p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
+    p = QuadraticSplineParameterCache(u, t, k, c, cache_parameters)
     A = QuadraticSpline(
         u, t, nothing, p, k, c, sc, extrapolation_left,
         extrapolation_right, cache_parameters, assume_linear_t
@@ -686,7 +686,7 @@ function QuadraticSpline(
         end
     end
 
-    p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
+    p = QuadraticSplineParameterCache(u, t, k, c, cache_parameters)
     A = QuadraticSpline(
         u, t, nothing, p, k, c, sc, extrapolation_left,
         extrapolation_right, cache_parameters, assume_linear_t

src/interpolation_methods.jl

@@ -1151,7 +1151,8 @@ function _interpolate(
     idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
     n = length(A.t)
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), n)
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zero(eltype(A.u))
     for i in nonzero_coefficient_idxs
@@ -1172,7 +1173,8 @@ function _interpolate(
     idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
     n = length(A.t)
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), n)
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...)
     for i in nonzero_coefficient_idxs
@@ -1188,7 +1190,8 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     # change t into param [0 1]
     idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), A.h)
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zero(eltype(A.u))
     for i in nonzero_coefficient_idxs
@@ -1206,7 +1209,8 @@ function _interpolate(
     # change t into param [0 1]
     idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
-    sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
+    # Per-call scratch buffer: evaluation must be reentrant for thread safety (#532)
+    sc = zeros(eltype(t), A.h)
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...)
     for i in nonzero_coefficient_idxs

src/interpolation_utils.jl

@@ -275,7 +275,7 @@ function get_parameters(A::QuadraticSpline, idx)
     return if A.cache_parameters
         A.p.α[idx], A.p.β[idx]
     else
-        quadratic_spline_parameters(A.u, A.t, A.k, A.c, A.sc, idx)
+        quadratic_spline_parameters(A.u, A.t, A.k, A.c, idx)
     end
 end
 

src/parameter_caches.jl

@@ -131,27 +131,40 @@ struct QuadraticSplineParameterCache{pType}
     β::pType
 end
 
-function QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
+function QuadraticSplineParameterCache(u, t, k, c, cache_parameters)
     return if cache_parameters
         parameters = quadratic_spline_parameters.(
-            Ref(u), Ref(t), Ref(k), Ref(c), Ref(sc), 1:(length(t) - 1)
+            Ref(u), Ref(t), Ref(k), Ref(c), 1:(length(t) - 1)
         )
         α, β = collect.(eachrow(stack(collect.(parameters))))
         QuadraticSplineParameterCache(α, β)
     else
         # Compute parameters once to infer types
-        α, β = quadratic_spline_parameters(u, t, k, c, sc, 1)
+        α, β = quadratic_spline_parameters(u, t, k, c, 1)
         QuadraticSplineParameterCache(typeof(α)[], typeof(β)[])
     end
 end
 
-function quadratic_spline_parameters(u, t, k, c, sc, idx)
+function quadratic_spline_parameters(u, t, k, c, idx)
     tᵢ₊ = (t[idx] + t[idx + 1]) / 2
-    nonzero_coefficient_idxs = spline_coefficients!(sc, 2, k, tᵢ₊)
+    # Value of the spline at the segment midpoint via the (degree 2) B-spline basis.
+    # The three nonzero basis values are evaluated into local variables rather than a
+    # shared scratch buffer so that evaluation is reentrant / thread-safe (#532).
+    # `tᵢ₊` is always interior, so only the non-boundary branch of the Cox-de Boor
+    # recursion is needed (cf. `spline_coefficients!`).
+    i = findfirst(x -> x > tᵢ₊, k)::Int - 1
+    w₁ = (k[i + 1] - tᵢ₊) / (k[i + 1] - k[i])
+    w₂ = (tᵢ₊ - k[i]) / (k[i + 1] - k[i])
+    N₁ = (k[i + 1] - tᵢ₊) / (k[i + 1] - k[i - 1]) * w₁
+    N₂ = (tᵢ₊ - k[i - 1]) / (k[i + 1] - k[i - 1]) * w₁ +
+        (k[i + 2] - tᵢ₊) / (k[i + 2] - k[i]) * w₂
+    N₃ = (tᵢ₊ - k[i]) / (k[i + 2] - k[i]) * w₂
+    # Seed the accumulator with `zero(first(u))` (as the buffer-based version did) so
+    # `uᵢ₊` keeps the element type of `u` for vector-valued data.
     uᵢ₊ = zero(first(u))
-    for j in nonzero_coefficient_idxs
-        uᵢ₊ += sc[j] * c[j]
-    end
+    uᵢ₊ += N₁ * c[i - 2]
+    uᵢ₊ += N₂ * c[i - 1]
+    uᵢ₊ += N₃ * c[i]
     α = 2 * (u[idx + 1] + u[idx]) - 4uᵢ₊
     β = 4 * (uᵢ₊ - u[idx]) - (u[idx + 1] - u[idx])
     return α, β

test/runtests.jl

@@ -15,6 +15,7 @@ const GROUP = get(ENV, "GROUP", "All")
     if GROUP == "All" || GROUP == "Methods"
         @testset "Methods" begin
             @safetestset "Derivative Tests" include("derivative_tests.jl")
+            @safetestset "Thread Safety Tests" include("thread_safety_tests.jl")
             @safetestset "Integral Tests" include("integral_tests.jl")
             @safetestset "Integral Inverse Tests" include("integral_inverse_tests.jl")
             @safetestset "Online Tests" include("online_tests.jl")

test/thread_safety_tests.jl

@@ -0,0 +1,68 @@
+using DataInterpolations
+using DataInterpolations: derivative
+using Test
+using Base.Threads
+
+# Regression tests for issue #532: a `BSplineInterpolation`, `BSplineApprox` or
+# (default, uncached) `QuadraticSpline` kept a scratch buffer in its `sc` field
+# and overwrote it in place on every evaluation. Because a "read" (evaluation)
+# mutated shared state, evaluating one shared interpolant from several threads
+# silently returned wrong values. Evaluation (and differentiation) must be
+# reentrant: it must not write to any state stored in the interpolant.
+
+function thread_safety_interpolants()
+    u = [0.0, 1.0, 0.0, 1.0, 0.0, 1.0]
+    t = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
+    umat = [0.0 1.0 0.0 1.0 0.0 1.0
+            1.0 0.0 1.0 0.0 1.0 0.0]            # 2 × length(t): last axis indexes t
+    return [
+        ("BSplineInterpolation", BSplineInterpolation(u, t, 3, :Uniform, :Average)),
+        # degree 1 == "linear" B-spline: same `_interpolate` method, so equally affected
+        ("BSplineInterpolation (degree 1)", BSplineInterpolation(u, t, 1, :Uniform, :Average)),
+        ("BSplineApprox", BSplineApprox(u, t, 3, 4, :Uniform, :Average)),
+        ("BSplineInterpolation (array)", BSplineInterpolation(umat, t, 3, :Uniform, :Average)),
+        ("BSplineApprox (array)", BSplineApprox(umat, t, 3, 4, :Uniform, :Average)),
+        ("QuadraticSpline", QuadraticSpline(u, t)),                       # cache_parameters = false
+        ("QuadraticSpline (cache_parameters)", QuadraticSpline(u, t; cache_parameters = true)),
+    ]
+end
+
+@testset "Thread safety / reentrancy (issue #532)" begin
+    pts = collect(range(1.0, 6.0; length = 41))   # includes the knots themselves
+
+    # Deterministic guard (catches the bug even on a single thread): repeated and
+    # interleaved evaluation must leave the interpolant's scratch buffer — and its
+    # results — unchanged.
+    @testset "evaluation does not mutate the interpolant: $name" for (name, A) in
+                                                                     thread_safety_interpolants()
+        sc_before = copy(A.sc)
+        ref = [A(x) for x in pts]
+        dref = [derivative(A, x) for x in pts]
+        for _ in 1:200, x in pts
+            A(x)
+            derivative(A, x)
+        end
+        @test A.sc == sc_before
+        @test [A(x) for x in pts] == ref
+        @test [derivative(A, x) for x in pts] == dref
+    end
+
+    # End-to-end check: many threads hammering one shared interpolant must all
+    # agree with the single-threaded reference. Only meaningful with >1 thread.
+    if Threads.nthreads() > 1
+        @testset "concurrent evaluation is race-free: $name" for (name, A) in
+                                                                 thread_safety_interpolants()
+            ref = [A(x) for x in pts]
+            wrong = Threads.Atomic{Int}(0)
+            Threads.@threads for _ in 1:50_000
+                for (j, x) in enumerate(pts)
+                    A(x) == ref[j] || Threads.atomic_add!(wrong, 1)
+                end
+            end
+            @test wrong[] == 0
+        end
+    else
+        @info "issue #532: concurrent stress test skipped; start Julia with more " *
+              "than one thread (e.g. `julia -t auto`) to exercise it."
+    end
+end