2023-08-09__Poisson_cquantile_upperbound¶

using WithFeedback

@withfb using Revise

using Revise … ✔ (1.9 s)

@withfb using GlobalMacros
@withfb using Units

using GlobalMacros … ✔
using Units … ✔

@withfb using PyPlot
@withfb using Sciplotlib
@withfb using PhDPlots

using PyPlot … ✔ (16.4 s)
using Sciplotlib … ✔ (11.1 s)
using PhDPlots … ✔

@withfb using Distributions

using Distributions … ✔ (1.5 s)

f(x) = cquantile(Poisson(x), 1E-14)

f (generic function with 1 method)

x = 1:10:(4*600)
g(x) = 300+1.05*x
Sciplotlib.plot(x, f.(x))
Sciplotlib.plot(x, g.(x));

../_images/2023-08-09__Poisson_cquantile_upperbound_8_0.png

x = 1:10:(4*60000)
all(f.(x) .< g.(x))

true

using Nto1AdEx

out = Nto1AdEx.sim(6500, 10*seconds);

PhDPlots.plotsig(out.V / mV, [0, 1000], ms);

../_images/2023-08-09__Poisson_cquantile_upperbound_12_0.png

@time Nto1AdEx.sim(6500, 10*seconds);

  0.067683 seconds (13.03 k allocations: 31.572 MiB, 21.20% gc time)

🤯

Measuring compilation time¶

re-pre-compiling the module
first time running the func

@time using Nto1AdEx

[ Info: Precompiling Nto1AdEx [368485ca-63bc-4029-854f-349d2205662c]

  0.727026 seconds (195.55 k allocations: 13.672 MiB, 2.12% gc time, 26.00% compilation time: 19% of which was recompilation)

@time using Units

  0.000880 seconds (377 allocations: 25.625 KiB)

@time Nto1AdEx.sim(6500, 10*seconds);

  0.041351 seconds (13.03 k allocations: 31.572 MiB, 27.28% gc time)

@time Nto1AdEx.sim(6500, 10*seconds);

  0.039387 seconds (13.03 k allocations: 31.572 MiB, 31.90% gc time)

:OOOOOO holy damn!
Changing laptop battery mode from ‘better battery’ to ‘best performance,
changes sim time from ~0.15 sec to 0.026 sec.

Gotta re-do Brian timings on ‘best performance’.

(Can change setting not only by clicking battery (which doesn’t always work 🙄), but also start > settings > power & bat > power mode.\ (https://support.microsoft.com/en-us/windows/change-the-power-mode-for-your-windows-pc-c2aff038-22c9-f46d-5ca0-78696fdf2de8#Category=Windows_11)

Hm, how can I trigger pkg recomp again? Maybe not using Revise (it’s too fast now, lol (when I add method or change a param).

Ok, that did it. (I restarted nb, didn’t load Revise). Then:

[ Info: Precompiling Nto1AdEx [368485ca-63bc-4029-854f-349d2205662c]

  0.727026 seconds (195.55 k allocations: 13.672 MiB, 2.12% gc time, 26.00% compilation time: 19% of which was recompilation)

First run time of sim:

 0.041351 seconds (13.03 k allocations: 31.572 MiB, 27.28% gc time)

hm, is this cached?

After clearing C:\Users\tfiers\.julia\compiled\v1.9\Nto1AdEx (had to stop this jl kernel)

@time using Nto1AdEx

[ Info: Precompiling Nto1AdEx [368485ca-63bc-4029-854f-349d2205662c]

  0.679249 seconds (193.93 k allocations: 13.567 MiB, 1.45% gc time, 24.78% compilation time: 10% of which was recompilation)

@time using Units

  0.000860 seconds (377 allocations: 25.625 KiB)

@time Nto1AdEx.sim(6500, 10*seconds);

  0.037517 seconds (13.03 k allocations: 31.572 MiB, 30.89% gc time)

@time Nto1AdEx.sim(6500, 10*seconds);

  0.035507 seconds (13.03 k allocations: 31.572 MiB, 28.18% gc time)

Ok. incredible.

Also, there’s just no compilation cost on first function call?

One last try: emptying module.

@time using Nto1AdEx

[ Info: Precompiling Nto1AdEx [368485ca-63bc-4029-854f-349d2205662c]

  0.657741 seconds (193.18 k allocations: 13.418 MiB, 2.26% gc time, 27.71% compilation time: 9% of which was recompilation)

(This was precompiling an empty module (just commented out lines))

Restarting and nb and re-filling in:
(but adding a comment line to not hit cache, maybe)

@time using Nto1AdEx

[ Info: Precompiling Nto1AdEx [368485ca-63bc-4029-854f-349d2205662c]

  0.693661 seconds (194.14 k allocations: 13.590 MiB, 1.51% gc time, 26.59% compilation time: 11% of which was recompilation)

What the hell? Mad? Is this real?

@time Nto1AdEx.sim(6500, 10*seconds);

  0.038055 seconds (13.03 k allocations: 31.572 MiB, 30.97% gc time)

Ok, one last test. A new module. Maybe here?

New module¶

module Nto1AdEx2

using GlobalMacros
using Units
using Random

@typed begin
    # AdEx LIF neuron params (cortical RS)
    C  = 104  * pF
    gₗ = 4.3  * nS
    Eₗ = -65  * mV
    Vₜ = -52  * mV
    Δₜ = 0.8  * mV
    Vₛ =  40  * mV
    Vᵣ = -53  * mV
    a  = -0.8 * nS
    b  =  65  * pA
    τw =  88  * ms
    # Conductance-based synapses
    Eₑ =   0 * mV
    Eᵢ = -80 * mV
    τg =   7 * ms
    # Simulation timestep
    Δt = 0.1 * ms
    # Input firing rate distribution
    μₓ = 4 * Hz
    σ = sqrt(0.6)
    μ = log(μₓ / Hz) - σ^2 / 2
end

const T = Float64

@kwdef mutable struct Neuron
    V    ::T = Eₗ
    w    ::T = 0 * pA
    gₑ   ::T = 0 * nS
    gᵢ   ::T = 0 * nS
    DₜV  ::T = 0 * mV/second
    Dₜw  ::T = 0 * pA/second
    Dₜgₑ ::T = 0 * nS/second
    Dₜgᵢ ::T = 0 * nS/second
end

# Calculate & store derivatives
f!(n::Neuron) = let (; V, w, gₑ, gᵢ) = n
    # Synaptic current
    Iₛ = gₑ*(V - Eₑ) + gᵢ*(V - Eᵢ)
    # Diffeqs
    n.DₜV  = (-gₗ*(V - Eₗ) + gₗ*Δₜ*exp((V-Vₜ)/Δₜ) - Iₛ - w) / C
    n.Dₜw  = (a*(V - Eₗ) - w) / τw
    n.Dₜgₑ = -gₑ / τg
    n.Dₜgᵢ = -gᵢ / τg
end

eulerstep!(n::Neuron) = begin
    n.V  += n.DₜV  * Δt
    n.w  += n.Dₜw  * Δt
    n.gₑ += n.Dₜgₑ * Δt
    n.gᵢ += n.Dₜgᵢ * Δt
end

has_spiked(n::Neuron) = n.V > Vₛ
on_self_spike!(n::Neuron) = begin
    n.V = Vᵣ
    n.w += b
end

include("../pkg/Nto1AdEx/src/poisson.jl")

struct Spike
    time   ::Float64
    source ::Int
end
Base.isless(x::Spike, y::Spike) = x.time < y.time

spikevec(src_id, times) = [Spike(t, src_id) for t in times]

multiplex(spiketrains) = begin
    # Join spiketimes of different trains into one sorted stream of Spikes
    spikevecs = [spikevec(i, times) for (i, times) in enumerate(spiketrains)]
    spikes = reduce(vcat, spikevecs)
    sort!(spikes)
end

sim(N, duration, seed=1, wₑ=14*pS, wᵢ=4*wₑ) = begin
    Random.seed!(seed)
    num_steps = round(Int, duration / Δt)
    t = 0 * second
    n = Neuron()
    V = Vector{T}(undef, num_steps)
    spiketimes = T[]
    Nₑ = round(Int, N * 4/5)
    rates = exp.(randn(N) .* σ .+ μ) .* Hz
    trains = [poisson_spikes(r, duration) for r in rates]
    spikes = multiplex(trains)
    # Index to keep track of input spikes processed
    j = 1
    for i in 1:num_steps
        # Process incoming spikes
        while j < length(spikes) && spikes[j].time < t
            # New spike arrival
            spike = spikes[j]
            if spike.source ≤ Nₑ
                n.gₑ += wₑ
            else
                n.gᵢ += wᵢ
            end
            j += 1
        end
        # Update neuron state variables (calc & apply diffeqs)
        f!(n)
        eulerstep!(n)
        if has_spiked(n)
            on_self_spike!(n)
            push!(spiketimes, t)
        end
        V[i] = n.V
        t += Δt
    end
    spikerate = length(spiketimes) / duration
    return (; V, spiketimes, rates, trains, Nₑ, spikerate, wₑ)
end

end

WARNING: replacing module Nto1AdEx2.

Main.Nto1AdEx2

@time using .Nto1AdEx2

  0.000009 seconds (9 allocations: 416 bytes)

using Units

@time Nto1AdEx2.sim(6500, 10*seconds);

  0.042777 seconds (13.03 k allocations: 31.572 MiB, 23.99% gc time)

@time Nto1AdEx2.sim(6500, 10*seconds);

  0.026069 seconds (13.03 k allocations: 31.572 MiB)

Absolute madness. Didn’t know Julia was this good here.

Voltage-to-Wiring Sim

2023-08-09__Poisson_cquantile_upperbound

Contents

2023-08-09__Poisson_cquantile_upperbound¶

Measuring compilation time¶

New module¶