Merge pull request #422 from Neuroblox/ho/firing_rate

Add firing rate calculator and recipe

ext/MakieExtension.jl

@@ -4,14 +4,16 @@ isdefined(Base, :get_extension) ? using Makie : using ..Makie
 
 using Neuroblox
 using Neuroblox: AbstractBlox, AbstractNeuronBlox, CompositeBlox, VLState, VLSetup
-using Neuroblox: meanfield_timeseries, voltage_timeseries, detect_spikes, get_neurons
+using Neuroblox: meanfield_timeseries, voltage_timeseries, detect_spikes, firing_rate, get_neurons
 using Neuroblox: powerspectrum
 using SciMLBase: AbstractSolution
 using LinearAlgebra: diag
-using DSP, SparseArrays
+using SparseArrays
+using DSP
 using Statistics: mean
 
-import Neuroblox: meanfield, meanfield!, rasterplot, rasterplot!, stackplot, stackplot!, voltage_stack, effectiveconnectivity, effectiveconnectivity!, ecbarplot, freeenergy, freeenergy!
+
+import Neuroblox: meanfield, meanfield!, rasterplot, rasterplot!, stackplot, stackplot!, frplot, frplot!, voltage_stack, effectiveconnectivity, effectiveconnectivity!, ecbarplot, freeenergy, freeenergy!
 import Neuroblox: powerspectrumplot, powerspectrumplot!
 
 @recipe(FreeEnergy, spDCMresults) do scene
@@ -95,20 +97,6 @@ end
 
 argument_names(::Type{<: RasterPlot}) = (:blox, :sol)
 
-# function Makie.plot!(p::RasterPlot)
-#     sol = p.sol[]
-#     t = sol.t
-#     blox = p.blox[]
-#     neurons = get_neurons(blox)
-
-#     for (i, n) in enumerate(neurons)
-#         spike_idxs = detect_spikes(n, sol)
-#         scatter!(p, t[spike_idxs], fill(i, length(spike_idxs)); color=p.color[])
-#     end
-
-#     return p
-# end
-
 function Makie.plot!(p::RasterPlot)
     sol = p.sol[]
     t = sol.t
@@ -120,8 +108,8 @@ function Makie.plot!(p::RasterPlot)
     ax.ylabel = p.Axis.ylabel[]
 
     spikes = detect_spikes(blox, sol; threshold=threshold)
-    neuron_indices, spike_times = findnz(spikes)
-    scatter!(p, spike_times, neuron_indices; color=p.color[])
+    spike_times, neuron_indices = findnz(spikes)
+    scatter!(p, sol.t[spike_times], neuron_indices; color=p.color[])
 
     return p
 end
@@ -159,6 +147,39 @@ function Makie.plot!(p::StackPlot)
     return p
 end
 
+@recipe(FRPlot, blox, sol) do scene
+    Theme(
+        color = :black,
+        Axis = (
+            ylabel = "Frequency (Hz)",
+            xlabel = "Time (s)"
+        ),
+        win_size = 10, # ms
+        overlap = 0,
+        transient = 0
+    )
+end
+
+argument_names(::Type{<: FRPlot}) = (:blox, :sol)
+
+function Makie.plot!(p::FRPlot)
+    sol = p.sol[]
+    blox = p.blox[]
+
+    ax = current_axis()
+    ax.xlabel = p.Axis.xlabel[]
+    ax.ylabel = p.Axis.ylabel[]
+
+    hideydecorations!(ax)
+
+    fr = firing_rate(blox, sol; win_size = p.win_size[], overlap = p.overlap[], transient = p.transient[])
+
+    t = range(p.transient[], stop = last(sol.t), length = length(fr))
+    lines!(p, t .* 1e-3, fr; color = p.color[])
+
+    return p
+end
+
 function Makie.convert_arguments(::Makie.PointBased, blox::AbstractNeuronBlox, sol::AbstractSolution)
     V = voltage_timeseries(blox, sol)
 

src/Neuroblox.jl

@@ -196,6 +196,9 @@ function rasterplot! end
 function stackplot end
 function stackplot! end
 
+function frplot end
+function frplot! end
+
 function voltage_stack end
 
 function ecbarplot end
@@ -249,8 +252,8 @@ export get_weights, get_dynamic_states, get_idx_tagged_vars, get_eqidx_tagged_va
 export BalloonModel,LeadField, boldsignal_endo_balloon
 export PINGNeuronExci, PINGNeuronInhib
 export PYR_Izh, QIF_PING_NGNMM
-export meanfield, meanfield!, rasterplot, rasterplot!, stackplot, stackplot!, voltage_stack, effectiveconnectivity, effectiveconnectivity!, ecbarplot, freeenergy, freeenergy!
+export meanfield, meanfield!, rasterplot, rasterplot!, stackplot, stackplot!, frplot, frplot!, voltage_stack, effectiveconnectivity, effectiveconnectivity!, ecbarplot, freeenergy, freeenergy!
 export powerspectrumplot, powerspectrumplot!, welch_pgram, periodogram, hanning, hamming
-export detect_spikes, mean_firing_rate
+export detect_spikes, mean_firing_rate, firing_rate
 export voltage_timeseries, meanfield_timeseries, state_timeseries
 end

src/blox/blox_utilities.jl

@@ -442,12 +442,29 @@ function detect_spikes(
     neurons = get_neurons(blox)
 
     S = mapreduce(sparse_hcat, neurons) do neuron
-        detect_spikes(neuron, sol; threshold)
+        detect_spikes(neuron, sol; threshold, ts)
     end
 
     return S
 end
 
+function firing_rate(
+    blox, sol::SciMLBase.AbstractSolution; 
+    win_size = last(sol.t), win_resolution = 1e-3, 
+    transient = 0, overlap = 0, threshold = nothing)
+
+    ts = sol.t
+    t_win_start = transient:(win_size - win_size*overlap):(last(ts) - win_size)
+
+    fr = map(t_win_start) do tws
+        spikes = detect_spikes(blox, sol; threshold, ts = tws:win_resolution:(tws + win_size))
+        N_neurons = size(spikes, 2)
+        1000.0 * (nnz(spikes) / N_neurons) / win_size
+    end
+
+    return fr
+end
+
 function mean_firing_rate(spikes::SparseMatrixCSC, sol; trim_transient = 0,
                  firing_rate_Δt = last(sol.t) - trim_transient,)