Merge branch 'master' into plot_event_timeline

code/chirpdetection.py

@@ -589,16 +589,16 @@ def main(datapath: str, plot: str) -> None:
     raw_time = np.arange(data.raw.shape[0]) / data.raw_rate
 
     # good chirp times for data: 2022-06-02-10_00
-    window_start_index = (3 * 60 * 60 + 6 * 60 + 43.5 + 5) * data.raw_rate
-    window_duration_index = 60 * data.raw_rate
+    # window_start_index = (3 * 60 * 60 + 6 * 60 + 43.5 + 5) * data.raw_rate
+    # window_duration_index = 60 * data.raw_rate
 
     #     t0 = 0
     #     dt = data.raw.shape[0]
     # window_start_seconds = (23495 + ((28336-23495)/3)) * data.raw_rate
     # window_duration_seconds = (28336 - 23495) * data.raw_rate
 
-    # window_start_index = 0
-    # window_duration_index = data.raw.shape[0]
+    window_start_index = 0
+    window_duration_index = data.raw.shape[0]
 
     # generate starting points of rolling window
     window_start_indices = np.arange(

code/modules/datahandling.py

@@ -1,9 +1,10 @@
 import numpy as np
 from typing import List, Any
 from scipy.ndimage import gaussian_filter1d
+from scipy.stats import gamma, norm
 
 
-def norm(data):
+def scale01(data):
     """
     Normalize data to [0, 1]
 
@@ -209,6 +210,117 @@ def flatten(list: List[List[Any]]) -> List:
     return [item for sublist in list for item in sublist]
 
 
+def causal_kde1d(spikes, time, width, shape=2):
+    """
+    causalkde computes a kernel density estimate using a causal kernel (i.e. exponential or gamma distribution).
+    A shape of 1 turns the gamma distribution into an exponential.
+
+    Parameters
+    ----------
+    spikes : array-like
+        spike times
+    time : array-like
+        sampling time
+    width : float
+        kernel width
+    shape : int, optional
+        shape of gamma distribution, by default 1
+
+    Returns
+    -------
+    rate : array-like
+        instantaneous firing rate
+    """
+
+    # compute dt
+    dt = time[1] - time[0]
+
+    # time on which to compute kernel:
+    tmax = 10 * width
+
+    # kernel not wider than time
+    if 2 * tmax > time[-1] - time[0]:
+        tmax = 0.5 * (time[-1] - time[0])
+
+    # kernel time
+    ktime = np.arange(-tmax, tmax, dt)
+
+    # gamma kernel centered in ktime:
+    kernel = gamma.pdf(
+        x=ktime,
+        a=shape,
+        loc=0,
+        scale=width,
+    )
+
+    # indices of spikes in time array:
+    indices = np.asarray((spikes - time[0]) / dt, dtype=int)
+
+    # binary spike train:
+    brate = np.zeros(len(time))
+    brate[indices[(indices >= 0) & (indices < len(time))]] = 1.0
+
+    # convolution with kernel:
+    rate = np.convolve(brate, kernel, mode="same")
+
+    return rate
+
+
+def acausal_kde1d(spikes, time, width):
+    """
+    causalkde computes a kernel density estimate using a causal kernel (i.e. exponential or gamma distribution).
+    A shape of 1 turns the gamma distribution into an exponential.
+
+    Parameters
+    ----------
+    spikes : array-like
+        spike times
+    time : array-like
+        sampling time
+    width : float
+        kernel width
+    shape : int, optional
+        shape of gamma distribution, by default 1
+
+    Returns
+    -------
+    rate : array-like
+        instantaneous firing rate
+    """
+
+    # compute dt
+    dt = time[1] - time[0]
+
+    # time on which to compute kernel:
+    tmax = 10 * width
+
+    # kernel not wider than time
+    if 2 * tmax > time[-1] - time[0]:
+        tmax = 0.5 * (time[-1] - time[0])
+
+    # kernel time
+    ktime = np.arange(-tmax, tmax, dt)
+
+    # gamma kernel centered in ktime:
+    kernel = norm.pdf(
+        x=ktime,
+        loc=0,
+        scale=width,
+    )
+
+    # indices of spikes in time array:
+    indices = np.asarray((spikes - time[0]) / dt, dtype=int)
+
+    # binary spike train:
+    brate = np.zeros(len(time))
+    brate[indices[(indices >= 0) & (indices < len(time))]] = 1.0
+
+    # convolution with kernel:
+    rate = np.convolve(brate, kernel, mode="same")
+
+    return rate
+
+
 if __name__ == "__main__":
 
     timestamps = [