Help with Julia multithreading

I am trying to multithread an estimation (optimization) problem over a set of parameters, while looping over another set of parameters. The optimization function uses a vector of random numbers. Here is the approximate code i am running in Julia 1.6.5:

pars1 = Vector( 1:8 )
lP = length(pars1)
Tup(x,y) = y>x ? (x,y) : return
Indices = [ Tup(x,y) for x in 1:(lP-1) for y in (2:lP)]
deleteat!( Indices, Indices .== nothing )
Indices = Vector{Tuple{Int64, Int64}}( Indices )
lI = length( Indices )
nChunks = Threads.nthreads();
lChunk = round( Int64, lI/nChunks, RoundUp )
rng = MersenneTwister(123)
Siter = 2000

function innerloop( Indices::Vector{Tuple{Int64, Int64}}, ch::Int64, lChunk::Int64, p1::Vector{Int64}, Siter::Int64, draws::Vector{Float64}, data::Vector{Float64}, arguments::Tuple{…} )
lI = length( Indices )
dfm = zeros( Float64, lChunk, 3 );
for j = 1 : lChunk
i = Int64( ch * lChunk - lChunk + j );
if i <= lI
p11 = p1[ Indices[ i ][ 1 ] ]
p12 = p1[ Indices[ i ][ 2 ] ]
res = Estimation(data, p11, p12, Siter, draws, arguments)

        dfm[ j, 1 ] = p11;
        dfm[ j, 2 ] = p12;
        dfm[ j, 3 ] = res[ 5 ];          
    end        
end
return dfm

end

cdfm = Vector{Any}(undef, nChunks)
Threads.@sync for ch = 1 : nChunks
Threads.@spawn begin
draws = rand(Siter)
cdfm[ ch ] = innerloop( Indices, ch, lChunk, pars1, Siter, draws, data, arguments )
end
end

The results I get are not what i expected based on the single thread run. res[5] is the maximum absolute error of the estimation and can be in some cases 100 orders of magnitude higher than in the single thread case.
Is there anything in the above code that could interfere with the multithread execution?
Many thanks.

Just a guess but

cdfm = Vector{Any}(undef, nChunks)
Threads.@sync for ch = 1 : nChunks
    Threads.@spawn begin
        draws = rand(Siter)
        cdfm[ch] = innerloop( Indices, ch, lChunk, pars1, Siter, draws, data, arguments )
    end
end

Looks like a race condition - a thread might be trying to write to cdfm before another thread has finished writing to it. Given your vector is undef you then might end up with random garbage in those entries of the array which haven’t been written to.

Thanks for the reply. I am rather new to Julia and previously work with Matlab, where this would not be an issue.
My understanding of data race was that it occurs when i’m writing and reading the same array from different threads. I thought i was not doing that here

I disagree with @nilshg . This part should not contain a race condition. In fact I can’t find a race condition in the code you posted. Maybe something inside Estimation mutates the input?

You could try and just copy everything you put into innerloop and see whether that fixes things.

cdfm = Vector{Any}(undef, nChunks)
Threads.@sync for ch = 1 : nChunks
  Threads.@spawn begin
    draws = rand(Siter)
    cdfm[ ch ] = innerloop( deepcopy(Indices), ch, lChunk, deepcoopy(pars1), Siter, draws, deepcopy(data), deepcopy(arguments) )
  end
end

Also please try and format your code a bit better. For starters enclose it in triple back ticks ``` like so
```
code here
```

EDIT: Actually maybe there is a race with the RNG. Try the slight variation to generate the random numbers in the main thread:

cdfm = Vector{Any}(undef, nChunks)
Threads.@sync for ch = 1 : nChunks
  draws = rand(Siter)
  Threads.@spawn begin
        cdfm[ ch ] = innerloop( Indices, ch, lChunk, pars1, Siter, draws, data, arguments )
  end
end

Hm you’re right I just checked, this bit seems to work fine with dummy data. I guess no threads will ever access the same index in cdfm so that bit works fine.

Can you try running the code in a new julia version?

The thread-safety of random number generation got overhauled in 1.7.

Also, your rng = MersenneTwister(123) is unused. Just delete that line.

If it is used in Estimation as a global variable:

  1. Don’t do that to us – post code for Estimation. (no hiding of bugs in code we can’t see)
  2. Don’t do that – pass the rng through as a parameter! (use globals with care, i.e. only under very specific circumstances)
  3. Don’t do that – MersenneTwister is not threadsafe and will blow up if you access it from multiple threads (possibly even to the point of corrupting memory, i.e. to the point of being a security vulnerability. I don’t remember the state of locks in dsfmt).

The deepcopy was the problem!
Many, many thanks

Thanks for the MersenneTwister advice. What should i be using instead?

i appreciate the help very much, but i cannot post the code for Estimation and I cannot run it in older versions of Julia either

You can just use the default rng.

In julia 1.6 this is a wrapper construction that is threadsafe and internally uses MersenneTwister (a separate one for every thread). This is not super fast and is not reproducible (via seed!) in context of multi-threading.

From julia 1.7 on, the default RNG is based on xoshiro and is reproducible in most multi-threaded contexts. Some of the details changed after 1.7, but I’m not up to date on them.

If you cannot use julia 1.7+, and you don’t need reproducibility in view of multi-threading, go with default. If you need reproducibility in view of multi-threading then you will need to use some library with small-state (e.g. counter-based) RNG and need to carefully split the RNGs for multi-threading. This is generally a pain, and one reason we changed the thing in 1.7.

Do you want advice on the “generally a pain” path? The requisite packages should still be maintained and in use, as that pain is still required when using Distributed (e.g. for clusters) instead of multithreading (same machine, same process, same address-space).

Glad that this fixed it. Btw, that confirms that you mutate the inputs inside Estimation and thus the race condition is not included in the code you posted. So in essence it was just a lucky guess to find the problem. Please next time post all of your code (if possible simplify it beforehand). Ideally you can even make it runnable (i.e. by filling arrays with some random numbers and setting all variables). That helps a lot in debugging :slight_smile: