# Why is this code slow?

**URL:** <https://discourse.julialang.org/t/why-is-this-code-slow/24403>\
**Category:** General Usage\
**Created:** [May 20, 2019, 2:21pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403 "2019-05-20T14:21:17Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [May 20, 2019, 2:21pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/1 "2019-05-20T14:21:17Z")

</div>

I am trying to draw a random index from a vector of probabilities. Here is my code:

```julia
import Random
using Random: AbstractRNG

struct FastCategorical
    _cdf::Vector{Float64}
    function FastCategorical(pdf)
        cdf = cumsum(pdf)
        @assert last(cdf) ≈ 1
        new(cdf)
    end
end

@inline function Random.rand(rng::AbstractRNG, o::FastCategorical)
    p = rand(rng)
    searchsortedfirst(o._cdf, p)
end

using BenchmarkTools
using LinearAlgebra
pdf = normalize!(rand(10^5), 1)
d = FastCategorical(pdf)
rng = Random.GLOBAL_RNG

@btime rand($rng, $d) # 128.708 ns (0 allocations: 0 bytes)
@btime searchsortedfirst($(d._cdf), x) setup=(x = rand()) # 39.083 ns (0 allocations: 0 bytes)
@btime rand($(rng)) # 2.865 ns (0 allocations: 0 bytes)

```

Observe that `rand` is much slower then running both its lines individually.

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 20, 2019, 3:43pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/2 "2019-05-20T15:43:42Z")

</div>

Apparently this problem has little to do with your structure of code

```julia
using BenchmarkTools, LinearAlgebra
function test()
    cdf = cumsum(normalize!(rand(10^5), 1))
    println("Summed times")
    @btime searchsortedfirst($cdf, rand())
    println("Single times")
    @btime searchsortedfirst($cdf, x) setup=(x=rand())
    @btime rand()
end
test()

--------

Summed times
  187.690 ns (0 allocations: 0 bytes)
Single times
  56.369 ns (0 allocations: 0 bytes)
  4.125 ns (0 allocations: 0 bytes)

```

:: Thanks to @giordano for making the testing rigorous

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [May 20, 2019, 3:56pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/3 "2019-05-20T15:56:28Z")

</div>

```julia
julia> function test3()
           cdf = cumsum(normalize!(rand(10^5), 1))
           f=rand(); 
           @btime searchsortedfirst($cdf, $(rand()))
           @btime searchsortedfirst($cdf, $f) 
           @btime searchsortedfirst($cdf, $(pi/10))
           @btime rand()
       end
test3 (generic function with 1 method)

julia> test3()
  44.017 ns (0 allocations: 0 bytes)
  44.198 ns (0 allocations: 0 bytes)
  44.163 ns (0 allocations: 0 bytes)
  3.012 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [May 20, 2019, 3:58pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/4 "2019-05-20T15:58:12Z")

</div>

CF [PSA: Microbenchmarks remember branch history](https://discourse.julialang.org/t/psa-microbenchmarks-remember-branch-history/17436)

It is even worse here, since your cdf does not fit into L1. You can probably get a significant speedup by using a cache-oblivious layout of the implicit search tree, instead of a linearly ordered vector.

Unfortunately, I am not aware of any julia packages implementing search-optimized layout of sorted lists.

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 20, 2019, 4:01pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/5 "2019-05-20T16:01:15Z")

</div>

@giordano thanks! my mistake, pi/3 also must be interpolated,  
but rand() was intentionally uninterpolated, in order for it to be calculated every time

Edited:

@foobar_lv2 explanatory post is right on the money.

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [May 20, 2019, 4:06pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/6 "2019-05-20T16:06:02Z")

</div>

> [@francesco.alemanno](#):
>
> but rand() was intentionally uninterpolated, in order for it to be calculated every time

You can use the `setup` keyword to `@btime` in that case, see [GitHub - JuliaCI/BenchmarkTools.jl: A benchmarking framework for the Julia language](https://github.com/JuliaCI/BenchmarkTools.jl#quick-start)

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 20, 2019, 4:12pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/7 "2019-05-20T16:12:56Z")

</div>

i was oversimplifying too much, and ended up making the test pointless, i edited according to your suggestions, thanks

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 20, 2019, 4:54pm UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/8 "2019-05-20T16:54:12Z")

</div>

Sidestepping the whole optimization question: if you want a lot of draws from a discrete/categorical distribution, it is worth building an [alias table](https://en.wikipedia.org/wiki/Alias_method).

It is implemented in Distributions.jl.

---

<div class="post-metadata">

**Author:** ![jw3126](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jw3126/32/3086_2.png) [@jw3126](https://discourse.julialang.org/u/jw3126)\
**Post date:** [May 21, 2019, 4:40am UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/9 "2019-05-21T04:40:47Z")

</div>

But why is the `@btime searchsortedfirst($cdf, x) setup=(x=rand())` fast? You tell me that the CPU learns to predict `searchsortedfirst` branches from the `rand` call??

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [May 21, 2019, 8:42am UTC](https://discourse.julialang.org/t/why-is-this-code-slow/24403/10 "2019-05-21T08:42:47Z")

</div>

> [@jw3126](#):
>
> But why is the `@btime searchsortedfirst($cdf, x) setup=(x=rand())` fast?

Because you also need to set `evals=1`. This is imo a design flaw in the BenchmarkTools API, but currently that’s how it works.

```julia
julia> @benchmark searchsortedfirst($cdf, x) setup=(x=rand())
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 66.402 ns (0.00% GC)
  median time: 71.518 ns (0.00% GC)
  mean time: 71.189 ns (0.00% GC)
  maximum time: 152.288 ns (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 978

julia> @benchmark searchsortedfirst($cdf, x) setup=(x=rand()) evals=1
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 149.000 ns (0.00% GC)
  median time: 283.000 ns (0.00% GC)
  mean time: 297.651 ns (0.00% GC)
  maximum time: 11.499 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 1

julia> g(a)=searchsorted(a, rand())
g (generic function with 2 methods)

julia> @benchmark g($cdf)
BenchmarkTools.Trial: 
  memory estimate: 0 bytes
  allocs estimate: 0
  --------------
  minimum time: 238.247 ns (0.00% GC)
  median time: 245.880 ns (0.00% GC)
  mean time: 249.613 ns (0.00% GC)
  maximum time: 5.811 μs (0.00% GC)
  --------------
  samples: 10000
  evals/sample: 425

```