Threading improvements

[palday]

better safety in nested threading

If you used use_threads=false but nested that within a Threads.@threads for, then you could be on a thread with id > 1 but with the shared resource arrays only 1 unit long. This problem was first detected by @behinger in a different context.

potential performance increase in multithreading.

For this, I replaced the ReentrantLock with a SpinLock. Spin locks are faster, but aren't re-entrant safe. We're not sharing this state anywhere and the RNG call is ideally quick, so we should be able to avoid the deadlock issues that the locks documentation warns about. Nonetheless, I've added a warning to the docstring discussing this.

The "potential" performance improvements are a complicated story, but it looks like we do get a boost for many models (the sleepstudy example is the big outlier). Here are some timing metrics (with progress bars and NLopt warnings cut out) using 4 threads on a machine with 8 cores on Julia 1.6-rc1 and default OpenBLAS and number of OpenBLAS threads:

user@box:~/Work/MixedModels.jl$ JULIA_NUM_THREADS=4 julia --project=.
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.6.0-rc1 (2021-02-06)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

[ Info: using Revise.jl
julia> Threads.nthreads()
4

julia> using BenchmarkTools

julia> using MixedModels

julia> using Random

julia> using MixedModels: dataset

julia> fm = fit(MixedModel, 
                @formula(reaction ~ 1 + days + (1 + days | subj) ),
                dataset(:sleepstudy));

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=false)
BenchmarkTools.Trial: 
  memory estimate:  48.18 MiB
  allocs estimate:  2591111
  --------------
  minimum time:     654.390 ms (0.57% GC)
  median time:      661.538 ms (0.52% GC)
  mean time:        662.146 ms (0.53% GC)
  maximum time:     670.782 ms (0.51% GC)
  --------------
  samples:          8
  evals/sample:     1

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=true)
BenchmarkTools.Trial: 
  memory estimate:  48.43 MiB
  allocs estimate:  2598271
  --------------
  minimum time:     1.579 s (0.00% GC)
  median time:      1.639 s (0.00% GC)
  mean time:        1.631 s (0.30% GC)
  maximum time:     1.668 s (1.19% GC)
  --------------
  samples:          4
  evals/sample:     1

julia> fm = fit(MixedModel, 
                @formula(yield ~ 1 + (1 | batch) ),
                dataset(:dyestuff));

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=false)
BenchmarkTools.Trial: 
  memory estimate:  10.78 MiB
  allocs estimate:  413441
  --------------
  minimum time:     41.938 ms (0.00% GC)
  median time:      43.071 ms (0.00% GC)
  mean time:        45.094 ms (3.35% GC)
  maximum time:     56.119 ms (12.58% GC)
  --------------
  samples:          111
  evals/sample:     1

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=true)
BenchmarkTools.Trial: 
  memory estimate:  10.84 MiB
  allocs estimate:  414659
  --------------
  minimum time:     13.829 ms (0.00% GC)
  median time:      15.324 ms (0.00% GC)
  mean time:        17.112 ms (7.55% GC)
  maximum time:     45.174 ms (56.36% GC)
  --------------
  samples:          293
  evals/sample:     1

julia> fm = fit(MixedModel,
                @formula(rt_trunc ~ 1 + spkr * prec * load + (1 | subj) + (1 | item)),
                dataset(:kb07));

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=false)
BenchmarkTools.Trial: 
  memory estimate:  30.83 MiB
  allocs estimate:  1327351
  --------------
  minimum time:     986.587 ms (0.00% GC)
  median time:      1.056 s (0.54% GC)
  mean time:        1.066 s (0.35% GC)
  maximum time:     1.134 s (0.61% GC)
  --------------
  samples:          5
  evals/sample:     1

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=true)
BenchmarkTools.Trial: 
  memory estimate:  32.29 MiB
  allocs estimate:  1331753
  --------------
  minimum time:     421.722 ms (0.00% GC)
  median time:      496.380 ms (0.00% GC)
  mean time:        488.109 ms (0.66% GC)
  maximum time:     517.069 ms (3.26% GC)
  --------------
  samples:          11
  evals/sample:     1

julia> fm = fit(MixedModel,
                @formula(rt_trunc ~ 1 + spkr * prec * load + (1 + spkr + prec + load | subj) + (1 + spkr + prec + load| item)),
                dataset(:kb07));

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=false)
BenchmarkTools.Trial: 
  memory estimate:  1.49 GiB
  allocs estimate:  91419845
  --------------
  minimum time:     293.846 s (0.04% GC)
  median time:      293.846 s (0.04% GC)
  mean time:        293.846 s (0.04% GC)
  maximum time:     293.846 s (0.04% GC)
  --------------
  samples:          1
  evals/sample:     1

julia> @benchmark parametricbootstrap($(MersenneTwister(1234321)), 1000, fm, use_threads=true)
BenchmarkTools.Trial: 
  memory estimate:  1.49 GiB
  allocs estimate:  91406667
  --------------
  minimum time:     170.185 s (0.15% GC)
  median time:      170.185 s (0.15% GC)
  mean time:        170.185 s (0.15% GC)
  maximum time:     170.185 s (0.15% GC)
  --------------
  samples:          1
  evals/sample:     1

[codecov]

Codecov Report

Merging #493 (1f297c3) into main (b9c73da) will increase coverage by 0.00%.
The diff coverage is 100.00%.

❗ Current head 1f297c3 differs from pull request most recent head 3c914ae. Consider uploading reports for the commit 3c914ae to get more accurate results

@@           Coverage Diff           @@
##             main     #493   +/-   ##
=======================================
  Coverage   93.95%   93.96%           
=======================================
  Files          26       26           
  Lines        2152     2153    +1     
=======================================
+ Hits         2022     2023    +1     
  Misses        130      130           

Impacted Files	Coverage Δ
src/bootstrap.jl	97.84% <100.00%> (+0.01%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b9c73da...3c914ae. Read the comment docs.

[dmbates]

Good catch that the threadid() can be other than 1 when nested.

[palday]

The original error came from @behinger actually using this within a threaded for loop.