I am aware that the two versions have to be different somewhere. But where? I Can’t see it. Also, please excuse the “Not Minimal Enough” Minimal Working Example
Summary
So I am writing a very nice example that highlights the strength of Julia and multiple dispatch. This will be used in conjuction to this discourse post: Idea: repository for examples highlighting Julia's strengths
In essense I am making a smaller version of DynamicalBilliards. In this tutorial, a billiard bd is a tuple of obstacles, that are of different type. When I am accessing this tuple with some index , e.g. o = bd[i], Julia allocates as this is a type unstable operation. It is not known in advance what type is obstacle bd[i]. However, with some magic that I don’t understand, it is actually possible to make this operation non allocating. Because, in DynamicalBilliards this is actually non-allocating, see below.
MWE
# barebone definitions
using StaticArrays
const SV = SVector{2,Float64}
mutable struct Particle; pos::SV; vel::SV; end
Particle(x0, y0, φ0) = Particle(SV(x0, y0), SV(cos(φ0), sin(φ0)))
struct Wall; sp::SV; ep::SV; normal::SV; end
struct Disk; c::SV; r::Float64; end
collisiontime(p::Particle, w::Wall) = cos(p.pos[1])
collisiontime(p::Particle, w::Disk) = sin(p.pos[1])
struct Billiard{L, O<:Tuple}
obstacles::O
end
Billiard(stuff::O) where {O<:Tuple} = Billiard{length(stuff), O}(stuff)
@inline Base.getindex(bd::Billiard, i) = bd.obstacles[i]
@inline Base.eachindex(bd::Billiard) = eachindex(bd.obstacles)
@inline Base.iterate(bd::Billiard, state = 1) = iterate(bd.obstacles, state)
@inline Base.length(bd::Billiard{L}) where {L} = L
# only relevant function
@generated function next_collision(p::Particle, bd::Billiard{L}) where {L}
out = :(j = 0; ct = Inf)
for i=1:L
push!(out.args, quote
let o = bd[$i]
t = collisiontime(p, o)
if t < ct
ct = t
j = $i
end
end
end
)
end
push!(out.args, :(return j, ct))
return out
end
# function where the allocation actually happens
@inline function bounce!(p::Particle, bd::Billiard)
i::Int, ct::Float64 = next_collision(p, bd)
o = bd[i] # <------ ALLOCATION HERE!
x = collisiontime(p, o) # <- irrelevant line, just so that `o` is used.
return i, x
end
##### Instantiate and benchmark:
x, y, r = 1.0, 1.0, 0.3
sp = [0.0,y]; ep = [0.0, 0.0]; n = [x,0.0]
leftw = Wall(sp, ep, n)
sp = [x,0.0]; ep = [x, y]; n = [-x,0.0]
rightw = Wall(sp, ep, n)
sp = [x,y]; ep = [0.0, y]; n = [0.0,-y]
topw = Wall(sp, ep, n)
sp = [0.0,0.0]; ep = [x, 0.0]; n = [0.0,y]
botw = Wall(sp, ep, n)
disk = Disk([x/2, y/2], r)
bd = Billiard((botw, rightw, topw, leftw, disk))
p = Particle(0.2, 0.1, 2π*rand())
using BenchmarkTools
@btime bounce!($p, $bd)
@profiler for i in 1:1000000; bounce!(p, bd); end
83.746 ns (2 allocations: 256 bytes)
The profiler results show that half of the function time is spent in the line o = bd[i].
Okay, I can understand that, it is a type unstable operation after all.
Same code, no allocation.
Exactly the process described above, is done in the DynamicalBilliards package. The function next_collision is identical (link) and the individual struct and collisiontime differences are irrelevant. To the best of my knowledge even the use of @inline is identical. The only differences are:
- The
Billiardstruct has an extra but irrelevant field. link - After the line
o = bd[i],bounce!calls some other functions that useoin a meaningful way. link
Yet, how can it be that (in Julia 1.0):
julia> using BenchmarkTools, DynamicalBilliards
julia> bd = billiard_sinai(); p = randominside(bd);
julia> @btime bounce!($p, $bd);
34.408 ns (0 allocations: 0 bytes)
??? (I promise you the line bd = billiard_sinai(); p = randominside(bd); replicates identical setting as the MWE). Doing profile on the code version of DynamicalBilliards shows that there is nearly 0 time spend on the line o = bd[i].
Why is this happening? :S :S I am super confused.
I’d really appreciate some help, as I have finished my nice example since two weeks ago, but I can’t post it since I don’t know why is this happening… 
p.s.: This has nothing to do with the MWE not being inside a module. I have already checked that. It also has nothing to do with collisiontime not using the identical functions as in DynamicalBilliards, I have also checked that. (in fact I just modified the MWE to not use them to save space)