Thread-safety

[cortner]

Is it possible that Dierckx.jl is not thread safe? evaluation is fine but derivatives seem to fail; see below for an explanation.

# file mtsplines.jl
using Dierckx, Random, .Threads
let 
   xdat = range(0, 1, length=100)
   spl = Spline1D(xdat, exp.(xdat))
   Random.seed!(0)
   r = rand(10_000)
   s = zeros(nthreads())
   @threads for n = 1:length(r)
      s[threadid()] += Dierckx.derivative(spl, r[n])
   end
   @show sum(s) / length(r)
end

This script yields the following output:

(base) Fuji-2:temp ortner$ j16 --threads=1 mtsplines.jl 
sum(s) / length(r) = 1.7170295968665945
(base) Fuji-2:temp ortner$ j16 --threads=1 mtsplines.jl 
sum(s) / length(r) = 1.7170295968665945
(base) Fuji-2:temp ortner$ j16 --threads=4 mtsplines.jl 
sum(s) / length(r) = -792.7283310922177
(base) Fuji-2:temp ortner$ j16 --threads=4 mtsplines.jl 
sum(s) / length(r) = -104735.52438361687

[cortner]

Looking at the code, it seems the work array might be the issue. For the record, here is a temporary work-around.

let 
   xdat = range(0, 1, length=100)
   spl = Spline1D(xdat, exp.(xdat))
   Random.seed!(0)
   r = rand(10_000)
   s = zeros(nthreads())
   wrk = [ zeros(length(spl.t)) for _=1:nthreads() ]
   @threads for n = 1:length(r)
      s[threadid()] += Dierckx._derivative(spl.t, spl.c, spl.k, r[n], 
                                           1, spl.bc, wrk[threadid()])
   end
   @show sum(s) / length(r)
end

[cortner]

I guess that as a workaround, I can pass in a separate work array for each thread, but how would people feelabout just allocating a separate array for each thread within Dierckx.jl?

The line
https://github.com/kbarbary/Dierckx.jl/blob/4eae3fc4af1acbfcef2f92eca7acacd65fab2ceb/src/Dierckx.jl#L72

would just become an array or work arrays then.

[cortner]

for the record, I'm happy to make a PR, especially since I was the one introducing the pre-allocation of the wrk array :); but would be good to get feedback first.

[kbarbary]

It seems like the only other option is to go back to allocating the work array in each function call, rather than in the constructor of Spline1D. If you think the "array of work arrays" option would be better, that sounds good to me!

ps: sorry for the ridiculous delay here. Industry job and life and all....

[cortner]

no worries at all, I understand :) - not in industry but still busy.... TBH, I'm not sure. the allocating approach will be much easier to maintain of course.

Yet another idea would be to implement a simple Object Pool. I've been using this extensively in another code of mine and love how well this works and how transparent things become.

[cortner]

I should also admit that I just switched my main code to Interpolations, but I might still use Dierckx in the future when I need the fitting. So from my end going back to allocating is fine and if I or anybody else has a performance problem then they can raise an issue and then one could do something more performant.

[kbarbary]

Yeah I think long term, Interpolations is probably the way to go. This package is probably nice for checking consistency and/or for people migrating code from scipy.interpolate.

[kbarbary]

Ease of maintenance is probably my top priority at this point given limited time for it!