Event handler

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?

How does newevent relate to the other variables?

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);

For example, something like this seems to work

function handler(x, T, f)
    d = length(x)
    event_type = zeros(Int, d)
    t = zeros(d)
    Q = queue(randexp(d))
    
    while true
        i, t′ = peek(Q)
        t′ > T && break 
        j = event_type[i] + 1
        f[j](i, t′, t, x, event_type, Q)
    end
    t, x
end
@time handler(zeros(d), T, (f1!, f2!));

but the small union

  %19 = Base.getindex(f, j)::Union{typeof(f1!), typeof(f2!)}

relies on the compiler

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 have many unknown event types, “JIT” approach may be useful. See, e.g., Catwalk.jl: [ANN] Catwalk.jl - With dynamic dispatch to the moon! (an adaptive optimizer, aka JIT compiler)

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.

2 Likes

So with GitHub - tisztamo/FunctionWranglers.jl: Fast, inlined execution of arrays of functions and some assistance of @tisztamo who added sindex for this use case I arrived at something neat

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.

PS: The scheduler is starting to take shape https://github.com/mschauer/ZigZagBoomerang.jl/…/src/scheduler.jl. It looks like something like this might find use in billiard problems too. @Datseris Somehow I just understood that https://github.com/mschauer/ZigZagBoomerang.jl and https://github.com/JuliaDynamics/DynamicalBilliards.jl are somewhat related.

1 Like

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.

Quiz: Is this a DynamicalBilliard or a Zig-Zag sampler of a density defined on the unit square:

billiard

@tisztamo

Absolutely nice:

with

  0.814835 seconds (33 allocations: 15.234 KiB

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).

Showing phase diagram of two coordinates out of 80.

1 Like

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).

well given that the “trajectory” changes directions in the middle of the square, it isn’t a billiard.

Well, yes, no,… it’s ray-splitting
at the contour lines of a density

1 Like