# Appropiate use of Threads.@threads?

**URL:** <https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480>\
**Category:** New to Julia\
**Tags:** question, multithreading\
**Created:** [June 27, 2017, 1:35am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480 "2017-06-27T01:35:45Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![zorp](https://avatars.discourse-cdn.com/v4/letter/z/c2a13f/32.png) [@zorp](https://discourse.julialang.org/u/zorp)\
**Post date:** [June 27, 2017, 1:35am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480/1 "2017-06-27T01:35:45Z")

</div>

I’ve written the following:

```julia
using DataFrames
using Query

plots = readtable("plots.csv")
species = readtable("plot_species.csv")

# remove taxa without a species or those that have [ or {'s in their name
species = @from i in species begin
	@where (search.(i.currentTaxonName, '[') == 0) &&
		(search.(i.currentTaxonName, '{') == 0) &&
		(search.(i.currentTaxonName, ' ') != 0)
	@select i
	@collect DataFrame
end

if !isfile("occurrence_counts.csv")
	occurrence_count = by(species, :currentTaxonName, nrow)
	rename!(occurrence_count, :currentTaxonName, :taxon)
	rename!(occurrence_count, :x1. :n)
	writetable("occurrence_counts.csv", occurrence_count)
else
	occurrence_counts = readtable("occurrence_counts.csv")
	unique_names = unique(convert(Array, occurrence_counts[:taxon]))
end

if !isfile("plot_alphas.csv")
	plot_alphas = by(species, :OBSERVATION_ID, df -> DataFrame(α=nrow))
	writetable("plot_alphas.csv", plot_alphas)
else
	plot_alphas = readtable("plot_alphas.csv")
end

###
### Calculating theta values
###
const N = 25
const ITER_COUNT = 100

function select_N_rand_plots(nums)
	s = Set{Int}()
	len = length(nums)

	while length(s) < N
		push!(s, nums[rand(1:len)])
	end

	return s
end

function calculate_mean_alpha(nums)
	cum_sum = 0
	
	for id in nums
		cum_sum += plot_alphas[plot_alphas[:id] .== id, :][:α][1]
	end
	
	cum_sum / length(nums)
end

function calculate_gamma(nums)
	spps = Set{String}()

	for id in nums
		occs = species[species[:OBSERVATION_ID] .== id, :][:currentTaxonName]
		for o in occs
			push!(spps, o)
		end
	end

	return length(spps)
end

unique_species = convert(Array, unique(species[:currentTaxonName]))
thetas = Float64[]
taxa = String[]
i = 1

for current_spp in unique_species
	occs = species[species[:currentTaxonName] .== current_spp, :]
	nrow(occs) < N && continue
	push!(taxa, current_spp)
	Θ = 0

	 Threads.@threads for a=1:ITER_COUNT
		sample_plots = select_N_rand_plots(convert(Array, occs[:OBSERVATION_ID]))

		γ = calculate_gamma(sample_plots)
		ᾱ = calculate_mean_alpha(sample_plots)
		Θ += (γ - ᾱ)
	end

	Θ /= ITER_COUNT
	push!(thetas, Θ)
	i += 1
	println("* finished $(current_spp) ($(i))")
end

results = DataFrame()
results[:taxa] = taxa
results[:Θ] = thetas
writetable("theta_species.csv", results)

```

…and I’m worried that my use of @threads is not correct. While working in Juno, if I run the entire script, with the @time macro placed on the outermost for loop, the amount of memory allocated seems **quite large** :

```julia
80.600184 seconds (210.53 M allocations: 45.076 GiB, 9.65% gc time)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 27, 2017, 1:46am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480/2 "2017-06-27T01:46:10Z")

</div>

> [@zorp](#):
>
> Θ = 0
> 
> ```
> Threads.@threads for a=1:ITER_COUNT
> sample_plots = select_N_rand_plots(convert(Array, occs[:OBSERVATION_ID]))
> 
> γ = calculate_gamma(sample_plots)
> ᾱ = calculate_mean_alpha(sample_plots)
> Θ += (γ - ᾱ)
> end
> 
> ```

that’s not thread-safe. Think about what each thread is doing. They will all grab the `Θ` and update it. What if two threads grab the `Θ` at the same time? They will both add and then write, but the second write won’t include the add of the first, so the value will be off. You should instead make an array with a separate `Θ` for each `threadid()`, reduce on each of those, and do a final reduce afterwards. The more expensive the computation the more rare this will be, but this code is open for this to occur and it will skew the results.

> [@zorp](#):
>
> function calculate\_gamma(nums)  
> spps = Set{String}()
> 
> ```
> for id in nums
> occs = species[species[:OBSERVATION_ID] .== id, :][:currentTaxonName]
> for o in occs
> push!(spps, o)
> end
> end
> 
> return length(spps)
> 
> ```
> 
> end

That’s a very expensive way to do that calculation. You don’t need `spps`, just the length, right? So instead of `push!`ing values, just keep a running counter for the length.

```julia
function calculate_gamma(nums)
	spps = 0

	for id in nums
		occs = species[species[:OBSERVATION_ID] .== id, :][:currentTaxonName]
                spps += length(occs)
	end

	return spps
end

```

But even then, `occs = species[species[:OBSERVATION_ID] .== id, :][:currentTaxonName]` that computation is going to allocate a lot. `species[:OBSERVATION_ID] .== id` will create a temporary array, and `species[species[:OBSERVATION_ID] .== id, :]` will then generate a temporary array for the slice. That could be a ton if the data set is decently sized! The first can be fixed by having a cache array to write into, then use a view:

```julia
cache = BitArray(length(species[:OBSERVATION_ID] ))
for id in nums
   cache .= species[:OBSERVATION_ID] .== id
   @view(species[cache,:])[:currentTaxonName]
end

```

I believe that should work, though likely you might want to handle the data slightly differently if you’re always cutting it by OBSERVATION\_ID.

That gives you an idea for the kinds of optimizations you can do. I would highly recommend optimizing this before multithreading, since there’s a ton that can be improved.

---

<div class="post-metadata">

**Author:** ![zorp](https://avatars.discourse-cdn.com/v4/letter/z/c2a13f/32.png) [@zorp](https://discourse.julialang.org/u/zorp)\
**Post date:** [June 27, 2017, 2:07am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480/3 "2017-06-27T02:07:22Z")

</div>

Thanks for the feedback! I have some questions about what you wrote:

> [@ChrisRackauckas](#):
>
> You should instead make an array with a separate Θ for each threadid(), reduce on each of those, and do a final reduce afterwards.

Is something like this:

```julia
Θ = Array{Float64}(Threads.nthreads())

	Threads.@threads for a=1:ITER_COUNT
		sample_plots = select_N_rand_plots(convert(Array, occs[:OBSERVATION_ID]))

		γ = calculate_gamma(sample_plots)
		ᾱ = calculate_mean_alpha(sample_plots)
		Θ[Threads.threadid()] = (γ - ᾱ)
	end

	Θ = reduce(+, Θ)
	Θ /= ITER_COUNT
	push!(thetas, Θ)

```

along the lines of what you were talking about?

> [@ChrisRackauckas](#):
>
> I believe that should work, though likely you might want to handle the data slightly differently if you’re always cutting it by OBSERVATION\_ID.
> 
> That’s a very expensive way to do that calculation. You don’t need spps, just the length, right? So instead of push!ing values, just keep a running counter for the length.
> 
> ```julia
> function calculate_gamma(nums)
> spps = 0
> 
> for id in nums
> occs = species[species[:OBSERVATION_ID] .== id, :][:currentTaxonName]
> spps += length(occs)
> end
> 
> return spps
> end
> 
> ```

I’m not sure that does what I was hoping to anymore. The function is supposed to return the number of unique species contained within the plots that have nums as their ID numbers. By adding the length of occs each iteration, I may very well be counting the same species (if it occurs at different plots in nums) several times, right?

> [@ChrisRackauckas](#):
>
> That gives you an idea for the kinds of optimizations you can do. I would highly recommend optimizing this before multithreading, since there’s a ton that can be improved.

Are there other areas that you see as low-hanging optimization fruits? I’m rather new to this and would really appreciate any additional input!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 27, 2017, 2:18am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480/4 "2017-06-27T02:18:12Z")

</div>

> [@zorp](#):
>
> The function is supposed to return the number of unique species contained within the plots that have nums as their ID numbers. By adding the length of occs each iteration, I may very well be counting the same species (if it occurs at different plots in nums) several times, right?

Oh, it’s supposed to be the unique species? Then yes, you’ll need to do something different, but there should be a way to re-use cache vectors to cut allocations down. I am not entirely sure why you need to plot it either. Can’t you get similar information directly from the data and not have to go through the plot?

> [@zorp](#):
>
> Are there other areas that you see as low-hanging optimization fruits?

With that amount of allocations, I would profile and find out ways to cut down the allocations. The best way to optimize is to target what is actually taking a long time.

> **[GitHub - timholy/ProfileView.jl: Visualization of Julia profiling data](https://github.com/timholy/ProfileView.jl)**
>
> Visualization of Julia profiling data. Contribute to timholy/ProfileView.jl development by creating an account on GitHub.

ProfileView.jl will give you a nice chart which shows where time is spent: use that to find out which lines you need to optimize. From experience I know that it will be something to with temporary arrays given how many allocations you have. There must be some point in there that is allocating a lot and clogging the whole system.

---

<div class="post-metadata">

**Author:** ![zorp](https://avatars.discourse-cdn.com/v4/letter/z/c2a13f/32.png) [@zorp](https://discourse.julialang.org/u/zorp)\
**Post date:** [June 27, 2017, 2:46am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480/5 "2017-06-27T02:46:43Z")

</div>

I installed ProfileView and think I have identified the culprit:

```julia
function calculate_gamma(nums)
    for id in nums
		---> occs = species[species[:OBSERVATION_ID] .== id, :][:currentTaxonName]
		...
	return length(spps)
end

```

Can you suggest any way in which I might reduce the time spent executing that line?

> [@ChrisRackauckas](#):
>
> I am not entirely sure why you need to plot it either. Can’t you get similar information directly from the data and not have to go through the plot?

Sorry for the confusion. When I said “plot,” I was talking about locations at which the data that I’m working with were sampled.

---

<div class="post-metadata">

**Author:** ![yuyichao](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yuyichao/32/20_2.png) [@yuyichao](https://discourse.julialang.org/u/yuyichao)\
**Post date:** [June 27, 2017, 3:02am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480/6 "2017-06-27T03:02:22Z")

</div>

> [@zorp](#):
>
> Θ = Array{Float64}(Threads.nthreads())

From what you are doing, you want `zeros` and `+=` in the loop.

---

<div class="post-metadata">

**Author:** ![zorp](https://avatars.discourse-cdn.com/v4/letter/z/c2a13f/32.png) [@zorp](https://discourse.julialang.org/u/zorp)\
**Post date:** [June 27, 2017, 3:33am UTC](https://discourse.julialang.org/t/appropiate-use-of-threads-threads/4480/7 "2017-06-27T03:33:40Z")

</div>

> [@yuyichao](#):
>
> From what you are doing, you want zeros and += in the loop.

You mean like this?

```julia
@time for current_spp in unique_species
	occs = species[species[:currentTaxonName] .== current_spp, :]
	nrow(occs) < N && continue
	push!(taxa, current_spp)
	Θ = zeros(Threads.nthreads())

	Threads.@threads for a=1:ITER_COUNT
		sample_plots = select_N_rand_plots(occs[:OBSERVATION_ID])

		γ = calculate_gamma(sample_plots)
		ᾱ = calculate_mean_alpha(sample_plots)
		Θ[Threads.threadid()] += (γ - ᾱ)
	end

	Θ = reduce(+, Θ)
	Θ /= ITER_COUNT
	push!(thetas, Θ)
	i += 1
	println("* finished $(current_spp) ($(i))")
end

```

I’ve reduced my time quite a bit thanks to the suggestion that I profile. I rewrote my gamma calculating function by removing the most expensive computation. Here’s my time now:

```julia
10.806910 seconds (45.52 M allocations: 2.283 GiB, 7.51% gc time)

```