From f16c9505942c5edf7c5fd7efdac75b771b7d5e68 Mon Sep 17 00:00:00 2001 From: Alec Loudenback Date: Fri, 29 May 2026 17:09:18 -0600 Subject: [PATCH 1/3] Fix data race in BSpline and QuadraticSpline evaluation (#532) BSplineInterpolation, BSplineApprox, and QuadraticSpline (with the default cache_parameters=false) kept a scratch buffer in their `sc` field and overwrote it in place on every evaluation. Because a "read" (evaluation) mutated shared state inside the object, evaluation was not reentrant: evaluating a single shared interpolant from multiple threads raced on `A.sc` and silently returned wrong values. - BSpline{Interpolation,Approx} (interpolation_methods.jl, derivatives.jl): allocate the spline-coefficient scratch per call instead of reusing `A.sc`. This matches the existing ForwardDiff.Dual branch and makes eval reentrant. - QuadraticSpline (parameter_caches.jl): compute the three nonzero degree-2 B-spline basis values in local variables instead of writing them into the shared `A.sc` buffer. Keeps evaluation allocation-free (verified for scalar and SVector data) while making it reentrant; also covers derivatives/integrals. - Add thread-safety regression tests: a deterministic no-mutation-on-eval check plus a concurrent stress test. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/derivatives.jl | 12 +++-- src/interpolation_caches.jl | 4 +- src/interpolation_methods.jl | 12 +++-- src/interpolation_utils.jl | 2 +- src/parameter_caches.jl | 29 ++++++++---- test/Methods/thread_safety_tests.jl | 68 +++++++++++++++++++++++++++++ 6 files changed, 108 insertions(+), 19 deletions(-) create mode 100644 test/Methods/thread_safety_tests.jl diff --git a/src/derivatives.jl b/src/derivatives.jl index 268e3ec0..40c74c1f 100644 --- a/src/derivatives.jl +++ b/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] diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index cdcd6995..35da14d4 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -668,7 +668,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, t_props @@ -710,7 +710,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, t_props diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index aecdeec9..00c9da1f 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1264,7 +1264,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 @@ -1285,7 +1286,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 @@ -1301,7 +1303,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 @@ -1319,7 +1322,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 diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index ceea6aa6..0807a346 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -301,7 +301,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 diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl index d30930a4..cccc4e7d 100644 --- a/src/parameter_caches.jl +++ b/src/parameter_caches.jl @@ -131,21 +131,21 @@ 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) uᵢ₊ = if length(t) == 2 # For 2 data points the knot vector has boundary multiplicity 2 < degree + 1, # so the B-spline basis cannot be evaluated at interior points; the spline @@ -153,11 +153,24 @@ function quadratic_spline_parameters(u, t, k, c, sc, idx) (u[1] + u[2]) / 2 else 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] uᵢ₊ end α = 2 * (u[idx + 1] + u[idx]) - 4uᵢ₊ diff --git a/test/Methods/thread_safety_tests.jl b/test/Methods/thread_safety_tests.jl new file mode 100644 index 00000000..3920b354 --- /dev/null +++ b/test/Methods/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)), + ("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 From a60ca024c1d3980a909326e23025ead3151ede70 Mon Sep 17 00:00:00 2001 From: Alec Loudenback Date: Fri, 29 May 2026 19:08:50 -0600 Subject: [PATCH 2/3] Cover degree-1 (linear) BSplineInterpolation in thread-safety tests A degree-1 B-spline goes through the same _interpolate method as higher degrees (dispatch is on type, not degree), so it was equally subject to the shared-buffer race; add it explicitly to document the coverage. Co-Authored-By: Claude Opus 4.8 (1M context) --- test/Methods/thread_safety_tests.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/Methods/thread_safety_tests.jl b/test/Methods/thread_safety_tests.jl index 3920b354..f7a327ec 100644 --- a/test/Methods/thread_safety_tests.jl +++ b/test/Methods/thread_safety_tests.jl @@ -19,6 +19,8 @@ function thread_safety_interpolants() ] # 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)), From d60cf288410669330ec71486c0d2fc198bbf6928 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 29 Jun 2026 07:09:13 -0400 Subject: [PATCH 3/3] Make BSpline evaluation allocation-free as well as thread-safe MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The thread-safety fix for #532 replaced the shared `sc` scratch field with a per-call `zeros(eltype(t), n)` buffer in BSpline{Interpolation,Approx} evaluation and differentiation. That made eval reentrant but reintroduced a heap allocation on every call (848 B for a length-99 scalar spline; QuadraticSpline was already allocation-free via its local-variable rewrite). A degree-`d` B-spline has only `d + 1` nonzero basis functions at any point, so a buffer sized to the full knot vector is unnecessary. `bspline_nonzero_coefficients` computes those `d + 1` values into a stack-allocated `MVector` (StaticArrays) and returns them as an `SVector` window plus the control-point offset. This is both allocation-free and reentrant — no shared or heap state per call. All four `_interpolate` and four `_derivative` BSpline methods use it. Because the guarantee requires a single statically-allocation-free path (AllocCheck sees every branch reachable from the call operator, including the extrapolation paths that dispatch to `_derivative`), the heap fallback is removed and the degree is bounded at construction (`d < 16`). The largest degree used anywhere in the repo is 3; degree ≥16 global B-splines are numerically degenerate. StaticArrays moves from a test-only dependency to a hard dependency. Allocations per call (length-99 scalar spline): value 848 B → 0 B, derivative 848 B → 0 B; BSplineApprox 128 B → 0 B. Array-valued eval still allocates its output array, which is inherent. Adds AllocCheck regression coverage for BSpline value and derivative evaluation to test/qa/alloc_tests.jl. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01TpSbnpEmw2Z2z5rxvHWLGn --- Project.toml | 4 ++-- src/DataInterpolations.jl | 1 + src/derivatives.jl | 44 ++++++++++++++++++++---------------- src/interpolation_caches.jl | 8 +++++++ src/interpolation_methods.jl | 36 +++++++++++++---------------- src/interpolation_utils.jl | 40 ++++++++++++++++++++++++++++++++ test/qa/alloc_tests.jl | 28 +++++++++++++++++++++++ 7 files changed, 119 insertions(+), 42 deletions(-) diff --git a/Project.toml b/Project.toml index 894332f7..cb763806 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [weakdeps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -73,11 +74,10 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "AllocCheck", "BenchmarkTools", "SafeTestsets", "SciMLTesting", "ChainRulesCore", "Mooncake", "Optim", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "StaticArrays", "Symbolics", "Unitful", "Zygote", "SparseConnectivityTracer"] +test = ["Aqua", "AllocCheck", "BenchmarkTools", "SafeTestsets", "SciMLTesting", "ChainRulesCore", "Mooncake", "Optim", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Unitful", "Zygote", "SparseConnectivityTracer"] diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl index ffaf8157..73f0c34b 100644 --- a/src/DataInterpolations.jl +++ b/src/DataInterpolations.jl @@ -9,6 +9,7 @@ using RecipesBase: RecipesBase, @recipe, @series using PrettyTables: PrettyTables, pretty_table using ForwardDiff: ForwardDiff using EnumX: EnumX, @enumx +using StaticArrays: MVector, SVector import FindFirstFunctions import FindFirstFunctions: Guesser diff --git a/src/derivatives.jl b/src/derivatives.jl index 40c74c1f..68bb015e 100644 --- a/src/derivatives.jl +++ b/src/derivatives.jl @@ -256,15 +256,16 @@ 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 - # 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] ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2]) else - for i in 1:(n - 1) - ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) + # Stack-allocated basis window: differentiation must be reentrant (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, n) + @inbounds for l in 1:m + i = offset + l - 1 + (1 <= i <= n - 1) || continue + ducum += vals[l] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end return ducum * A.d * scale @@ -281,16 +282,17 @@ 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 - # 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] ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2]) else - for i in 1:(n - 1) + # Stack-allocated basis window: differentiation must be reentrant (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, n) + @inbounds for l in 1:m + i = offset + l - 1 + (1 <= i <= n - 1) || continue ducum = ducum + - sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / + vals[l] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end @@ -304,15 +306,16 @@ 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 - # 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] ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2]) else - for i in 1:(A.h - 1) - ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) + # Stack-allocated basis window: differentiation must be reentrant (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, A.h) + @inbounds for l in 1:m + i = offset + l - 1 + (1 <= i <= A.h - 1) || continue + ducum += vals[l] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end return ducum * A.d * scale @@ -328,15 +331,16 @@ 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 - # 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] ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2]) else - for i in 1:(A.h - 1) - ducum += sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / + # Stack-allocated basis window: differentiation must be reentrant (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, A.h) + @inbounds for l in 1:m + i = offset + l - 1 + (1 <= i <= A.h - 1) || continue + ducum = ducum + vals[l] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / (A.k[i + A.d + 1] - A.k[i + 1]) end end diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl index 35da14d4..ea59cb85 100644 --- a/src/interpolation_caches.jl +++ b/src/interpolation_caches.jl @@ -1020,6 +1020,8 @@ function BSplineInterpolation( t_props = something(search_properties, FindFirstFunctions.SearchProperties(t)) n = length(t) n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d + 1) points.") + d ≥ BSPLINE_STACK_MAXLEN && + error("BSplineInterpolation supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.") s = zero(eltype(u)) p = zero(t) k = zeros(eltype(t), n + d + 1) @@ -1101,6 +1103,8 @@ function BSplineInterpolation( t_props = something(search_properties, FindFirstFunctions.SearchProperties(t)) n = length(t) n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d + 1) points.") + d ≥ BSPLINE_STACK_MAXLEN && + error("BSplineInterpolation supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.") s = zero(eltype(u)) p = zero(t) k = zeros(eltype(t), n + d + 1) @@ -1274,6 +1278,8 @@ function BSplineApprox( t_props = something(search_properties, FindFirstFunctions.SearchProperties(t)) n = length(t) h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d + 1) control points.") + d ≥ BSPLINE_STACK_MAXLEN && + error("BSplineApprox supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.") s = zero(eltype(u)) p = zero(t) k = zeros(eltype(t), h + d + 1) @@ -1376,6 +1382,8 @@ function BSplineApprox( t_props = something(search_properties, FindFirstFunctions.SearchProperties(t)) n = length(t) h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d + 1) control points.") + d ≥ BSPLINE_STACK_MAXLEN && + error("BSplineApprox supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.") s = zero(eltype(u)) p = zero(t) k = zeros(eltype(t), h + d + 1) diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl index 00c9da1f..1929a0cf 100644 --- a/src/interpolation_methods.jl +++ b/src/interpolation_methods.jl @@ -1264,12 +1264,11 @@ 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) - # 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) + # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, n) ucum = zero(eltype(A.u)) - for i in nonzero_coefficient_idxs - ucum += sc[i] * A.c[i] + @inbounds for l in 1:m + ucum += vals[l] * A.c[offset + l] end return ucum end @@ -1286,12 +1285,11 @@ 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) - # 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) + # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, n) ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...) - for i in nonzero_coefficient_idxs - ucum = ucum + (sc[i] * A.c[ax_u..., i]) + @inbounds for l in 1:m + ucum = ucum + (vals[l] * A.c[ax_u..., offset + l]) end return ucum end @@ -1303,12 +1301,11 @@ 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]) - # 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) + # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, A.h) ucum = zero(eltype(A.u)) - for i in nonzero_coefficient_idxs - ucum += sc[i] * A.c[i] + @inbounds for l in 1:m + ucum += vals[l] * A.c[offset + l] end return ucum end @@ -1322,12 +1319,11 @@ 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]) - # 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) + # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532) + vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, A.h) ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...) - for i in nonzero_coefficient_idxs - ucum = ucum + (sc[i] * A.c[ax_u..., i]) + @inbounds for l in 1:m + ucum = ucum + (vals[l] * A.c[ax_u..., offset + l]) end return ucum end diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl index 0807a346..72bf7c6f 100644 --- a/src/interpolation_utils.jl +++ b/src/interpolation_utils.jl @@ -67,6 +67,46 @@ function spline_coefficients!(N, d, k, u::AbstractVector) return nothing end +# A degree-`d` B-spline has only `d + 1` nonzero basis functions at any point, so +# evaluation needs a scratch buffer of just `d + 1` entries — not one sized to the +# full knot/control-point vector. `bspline_nonzero_coefficients` computes those +# values into a stack-allocated `MVector` (no heap allocation, hence reentrant / +# thread-safe, #532) for degrees up to `BSPLINE_STACK_MAXLEN - 1`; callers fall back +# to the heap `spline_coefficients!` for larger degrees, which essentially never occur. +const BSPLINE_STACK_MAXLEN = 16 + +# Returns `(vals, offset, m)`: `vals[l]` is the basis weight of control point +# `offset + l` for `l in 1:m`. `ncp` is the number of control points (needed only to +# place the single nonzero weight at the right boundary). Mirrors the locator and +# recurrence of `spline_coefficients!`, but writes into a local window indexed `1:m` +# instead of absolute knot indices. +@inline function bspline_nonzero_coefficients(d::Integer, k, u::T, ncp::Integer) where {T} + N = zero(MVector{BSPLINE_STACK_MAXLEN, T}) + if u == k[1] + @inbounds N[1] = one(u) + return SVector(N), 0, 1 + elseif u == k[end] + @inbounds N[1] = one(u) + return SVector(N), ncp - 1, 1 + end + i = searchsortedlast(k, u) + # Local index of global knot index `g` is `g + off` (so `i - d → 1`, `i → d + 1`). + off = d + 1 - i + @inbounds begin + N[i + off] = one(u) + for deg in 1:d + N[i - deg + off] = (k[i + 1] - u) / (k[i + 1] - k[i - deg + 1]) * + N[i - deg + 1 + off] + for j in (i - deg + 1):(i - 1) + N[j + off] = (u - k[j]) / (k[j + deg] - k[j]) * N[j + off] + + (k[j + deg + 1] - u) / (k[j + deg + 1] - k[j + 1]) * N[j + 1 + off] + end + N[i + off] = (u - k[i]) / (k[i + deg] - k[i]) * N[i + off] + end + end + return SVector(N), i - d - 1, d + 1 +end + function quadratic_spline_params(t::AbstractVector, sc::AbstractVector) # Duplicate time points make the collocation system singular if any(i -> t[i] == t[i + 1], 1:(length(t) - 1)) diff --git a/test/qa/alloc_tests.jl b/test/qa/alloc_tests.jl index 4235f843..a5e187b4 100644 --- a/test/qa/alloc_tests.jl +++ b/test/qa/alloc_tests.jl @@ -1,8 +1,10 @@ using DataInterpolations +using DataInterpolations: derivative using AllocCheck: @check_allocs using StaticArrays: SVector @check_allocs(test_allocs(itp, x) = itp(x)) # Reuse function definition to save on compilation time +@check_allocs(test_derivative_allocs(itp, x) = derivative(itp, x)) @testset "Allocation-free interpolation tests" begin @testset "LinearInterpolation" begin @@ -53,6 +55,32 @@ using StaticArrays: SVector @test_nowarn test_allocs(A_s, 0.7) end + # B-spline evaluation must be allocation-free *and* reentrant: it computes the + # d+1 nonzero basis functions into a stack buffer instead of a shared scratch + # field, which is what makes it thread-safe (#532). Scalar data keeps the return + # value itself off the heap. + @testset "BSplineInterpolation" begin + t = collect(range(0.0, 2π; length = 12)) + u = sin.(t) + for d in (1, 2, 3) + A = BSplineInterpolation( + u, t, d, :Uniform, :Average; extrapolation = ExtrapolationType.Extension + ) + @test_nowarn test_allocs(A, 3.1) + @test_nowarn test_derivative_allocs(A, 3.1) + end + end + + @testset "BSplineApprox" begin + t = collect(range(0.0, 2π; length = 12)) + u = sin.(t) + A = BSplineApprox( + u, t, 3, 8, :Uniform, :Average; extrapolation = ExtrapolationType.Extension + ) + @test_nowarn test_allocs(A, 3.1) + @test_nowarn test_derivative_allocs(A, 3.1) + end + @testset "CubicSpline" begin t = [-1.0, 0.0, 1.0] u_ = [0.0, 1.0, 3.0]' .* ones(4)