Evaluation writes iguesser.idx_prev, so threads sharing an interpolant contend on it #618

# Evaluation cost with and without the `idx_prev` write, by grid size, knot spacing and query order.
# Run once per DataInterpolations checkout (master, then the no-write branch), e.g.
#   julia --project=<env with that checkout dev'd> -t 8 bench_sizes.jl
# Clustered knots (cubed) don't look linear, so the hint is the cached `idx_prev`; jittered knots look
# linear, so the hint is computed from `t` and `idx_prev` is never read.
using DataInterpolations, Random, Printf

const NQ = 1_000_000
work(A, x) = (s = 0.0; for xi in x; s += A(xi); end; s)
best(f) = minimum(@elapsed(f()) for _ in 1:5)
function threaded(interps, xs)
    out = zeros(length(xs))
    Threads.@threads :static for k in eachindex(xs)
        out[k] = work(interps[k], xs[k])
    end
    return out
end

println("DataInterpolations ", pkgversion(DataInterpolations), " at ", pkgdir(DataInterpolations), ", ", Threads.nthreads(), " threads")
for n in (12, 100, 1_000, 100_000), spacing in (:clustered, :jittered), method in (:Linear, :CubicSpline)
    t = spacing === :clustered ? 100 .* range(0.0, 1.0; length = n) .^ 3 :
        sort!(100 .* (range(0.0, 1.0; length = n) .+ (0.2 / n) .* rand(Xoshiro(1), n)))
    u = sin.(t ./ 7)
    A = method === :Linear ? LinearInterpolation(u, t) : CubicSpline(u, t)
    lo, hi = first(t), last(t)
    random = lo .+ (hi - lo) .* rand(Xoshiro(2), NQ)
    sorted = sort(random)
    xs = [lo .+ (hi - lo) .* rand(Xoshiro(10 + k), NQ) for k in 1:Threads.nthreads()]
    shared, copies = fill(A, length(xs)), [deepcopy(A) for _ in xs]
    work(A, sorted); work(A, random); threaded(shared, xs); threaded(copies, xs)
    @printf("%-7d %-9s %-11s %-18s lookup=%-5s | 1 thread: sorted %7.2f ms, random %7.2f ms | %d threads random: shared %7.2f ms, copies %7.2f ms\n",
        n, spacing, method, A.kind, A.iguesser.linear_lookup,
        1e3 * best(() -> work(A, sorted)), 1e3 * best(() -> work(A, random)),
        length(xs), 1e3 * best(() -> threaded(shared, xs)), 1e3 * best(() -> threaded(copies, xs)))
end

# What the idx_prev write buys and costs. VARIANT=nowrite replaces get_idx with the same lookup minus the
# write (the hint is still read); VARIANT=changed writes only when the index differs from the cached one.
# Run: julia -t 8 experiment.jl,  VARIANT=nowrite julia -t 8 experiment.jl,  VARIANT=changed julia -t 8 experiment.jl
using DataInterpolations, Random
const Guesser = DataInterpolations.Guesser   # FindFirstFunctions.Guesser, as DataInterpolations imports it
const VARIANT = get(ENV, "VARIANT", "current")
if VARIANT == "nowrite"
    @eval DataInterpolations function get_idx(A::AbstractInterpolation, t, iguess::Guesser; lb = 1, ub_shift = -1, idx_shift = 0, side = :last)
        raw = _dispatch_search(A, A.t, t, iguess(t), side)
        return clamp(raw + idx_shift, lb, length(A.t) + ub_shift)
    end
elseif VARIANT == "changed"
    @eval DataInterpolations function get_idx(A::AbstractInterpolation, t, iguess::Guesser; lb = 1, ub_shift = -1, idx_shift = 0, side = :last)
        raw = _dispatch_search(A, A.t, t, iguess(t), side)
        idx = clamp(raw + idx_shift, lb, length(A.t) + ub_shift)
        iguess.idx_prev[] == idx || (iguess.idx_prev[] = idx)
        return idx
    end
end
t = [0.5, 1.0, 2.0, 3.0, 5.0, 7.0, 10.0, 12.0, 15.0, 20.0, 25.0, 30.0]
u = 0.03 .+ 0.01 .* (1 .- exp.(-t ./ 5))
A = LinearInterpolation(u, t; cache_parameters = true)
const N = 2_000_000
xs = [0.5 .+ 29.5 .* rand(Xoshiro(k), N) for k in 1:Threads.nthreads()]
sorted = sort(xs[1])
work(A, x) = (s = 0.0; for xi in x; s += A(xi); end; s)
function run(interps)
    out = zeros(Threads.nthreads())
    Threads.@threads :static for k in 1:Threads.nthreads()
        out[k] = work(interps[k], xs[k])
    end
    return out
end
shared = fill(A, Threads.nthreads()); copies = [deepcopy(A) for _ in 1:Threads.nthreads()]
run(shared); run(copies); work(A, sorted); work(A, xs[1])
ms(f) = round(1e3 * minimum(@elapsed(f()) for _ in 1:5); digits = 1)
println(rpad(VARIANT, 8),
    " | threaded random: shared ", ms(() -> run(shared)), " ms, copies ", ms(() -> run(copies)),
    " ms | one thread: sorted ", ms(() -> work(A, sorted)), " ms, random ", ms(() -> work(A, xs[1])), " ms")