I have a number of functions, all looking a bit like that:
while true
i = peel(Q)
if i == 1 && mode = 1
a, b = dosomething1(a, b)
elseif i == 2
b, c = dosomething2(b, c)
elseif i == 3
c, d = dosomething2(c, d)
end
enqueue!(Q, newevent)
end
so these are event handlers based on a priority queue.
To me it looks like I want to be a bit more meta here and have something like an “engine”, but I have to be a bit careful as these are hot loops. Context ZigZagBoomerang.jl. Any pointers?
Let me give a more serious minimal working example with two event types:
So the Queue keeps only track of event times, and I have a vector for each event type for each element of a state
using ZigZagBoomerang: SPriorityQueue, enqueue!
using Random
function queue(τ)
Q = SPriorityQueue{Int,eltype(τ)}()
for i in eachindex(τ)
enqueue!(Q, i => τ[i])
end
Q
end
T = 100.0
d = 10000
state = rand(d)
function handler(x, T)
d = length(x)
event_type = zeros(Int, d)
t = zeros(d)
Q = queue(randexp(d))
while true
i, t′ = peek(Q)
t′ > T && break
if event_type[i] == 0 # fast
x[i] += (t′ - t[i])
t[i] = t′
event_type[i] = 1 - event_type[i] # next slow
Q[i] = t[i] + randexp()
else # slow
x[i] += 0.5*(t′ - t[i])
t[i] = t′
event_type[i] = 1 - event_type[i] # next fast
Q[i] = t[i] + 0.5*randexp()
end
end
t, x
end
@time handler(zeros(d), T);
If you only have a couple of known events, I think a sequence of if is a decent approach. If the computation inside state transition is relatively small and the transitions between events have some structure (not all states jump to arbitrary state), maybe you can squeeze out some performance with @goto.
If you know the possible concrete types of f[j] statically, you can use a manual union splitting hack like this. (This is a manual version of Catwalk.jl etc.). This is equivalent to plain ifs, but it can be useful for code organization.
function handle!(t, x, θ, m, f!, next, event_type, Q)
i, t′ = peek(Q)
j = event_type[i]
if j > 0
sindex(f!, j, i, t′, t, x, θ, m, event_type, Q)
end
τ, jnext = sfindmin(next, i, t′, t, x, θ, m)
Q[i] = τ
event_type[i] = jnext
(t′, i, x[i], j)
end
Getting the hang of it I added sfindmin to schedule the next relevant event for each particle.
Hi, Sorry I don’t see the connection between the two; I’m not familiar with the contents of ZigZagBoomerang.jl by having a look at its page its the first time I see/read any of it.
for 200000 calls of the event loop computing the random event times and the wall collision times for all 3 “walls” for each of 80 correlated coordinates (hard wall at 1, soft or hard wall at 0 and “sticky” barrier at 0.25).
Great to hear that FunctionWranglers works for you!
I do not really understand what ZigZagBoomerang does, but I see the three large parts on the flame graph and that there is not much extra time between them. If the number of event dispatches is 200_000 * 3 * 80 during this run, then your gain must be very significant to other solutions I think (Not counting compilation).
