# ============================================================================
# Data race in `BSplineInterpolation` under concurrent evaluation
# ============================================================================
#
# WHAT: Evaluating a single, shared `BSplineInterpolation` from several
# threads silently returns WRONG values (no error is thrown).
#
# WHY: `BSplineInterpolation` keeps a preallocated coefficient buffer in its
# `sc` field and OVERWRITES it in place on every evaluation
# (src/interpolation_methods.jl):
#
# # reuse the shared buffer A.sc (only the Dual branch allocates fresh)
# sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
# spline_coefficients!(sc, A.d, A.k, t) # writes basis(t) into A.sc
# for i in nonzero_coefficient_idxs
# ucum += sc[i] * A.c[i] # then reads A.sc back
# end
#
# So a "read" (evaluation) mutates state inside the object, and it is
# not reentrant. Two threads at different points race on `A.sc`: one
# thread's `spline_coefficients!` overwrites the buffer between the
# other thread's write and read, producing a garbage result.
#
# FIX: Take the same fresh-buffer path the `ForwardDiff.Dual` branch already
# uses for the plain-number case (`sc = zeros(eltype(t), n)`), or
# document `BSplineInterpolation` as non-reentrant for concurrent eval.
#
# NOTE: Local interpolants (`CubicSpline`, `LinearInterpolation`, PCHIP, ...)
# recompute each segment's parameters from read-only fields per call and
# keep no such scratch buffer, so they are unaffected — see Part 3.
#
# RUN: julia -t auto di_threadsafety_mwe.jl
# ----------------------------------------------------------------------------
using DataInterpolations
using Base.Threads
println("DataInterpolations v", pkgversion(DataInterpolations),
" | Julia ", VERSION, " | threads = ", nthreads())
# Minimal toy data — the bug is general, nothing here is domain-specific.
# The two probe points (2.5, 5.5) have very different values, so a corrupted
# result is unmistakable rather than off in the last bit.
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]
itp = BSplineInterpolation(u, t, 3, :Uniform, :Average)
# ── Part 1 — the root cause, with NO concurrency ────────────────────────────
# Evaluate at one point, then another: the SAME internal buffer `itp.sc`
# changes. In other words, evaluating the curve mutates hidden state.
itp(2.5); buffer_after_2_5 = copy(itp.sc)
itp(5.5); buffer_after_5_5 = copy(itp.sc)
println("\n[1] itp.sc is one buffer, overwritten on every evaluation: ",
buffer_after_2_5 != buffer_after_5_5) # => true
# ── Part 2 — the consequence under threads ──────────────────────────────────
# Each point has one correct value (computed single-threaded just below). A
# thread-safe curve would always return it. Evaluating the two points
# concurrently, the shared buffer gets clobbered and some results are wrong.
@assert nthreads() > 1 "Start Julia with >1 thread, e.g. `julia -t auto`."
correct_2_5, correct_5_5 = itp(2.5), itp(5.5)
wrong = Atomic{Int}(0)
@threads for _ in 1:1_000_000
itp(2.5) == correct_2_5 || atomic_add!(wrong, 1)
itp(5.5) == correct_5_5 || atomic_add!(wrong, 1)
end
println("[2] BSplineInterpolation: ", wrong[], " wrong / 2,000,000 evaluations") # thousands
# ── Part 3 — control: a local interpolant keeps no scratch buffer ───────────
safe = CubicSpline(u, t)
ok_2_5, ok_5_5 = safe(2.5), safe(5.5)
wrong_safe = Atomic{Int}(0)
@threads for _ in 1:1_000_000
safe(2.5) == ok_2_5 || atomic_add!(wrong_safe, 1)
safe(5.5) == ok_5_5 || atomic_add!(wrong_safe, 1)
end
println("[3] CubicSpline (control): ", wrong_safe[], " wrong / 2,000,000 evaluations") # 0# running the above produces
DataInterpolations v8.10.0 | Julia 1.12.5 | threads = 10
[1] itp.sc is one buffer, overwritten on every evaluation: true
[2] BSplineInterpolation: 1591 wrong / 2,000,000 evaluations
[3] CubicSpline (control): 0 wrong / 2,000,000 evaluations
Error & Stacktrace⚠️