[MRG] Simulate spikes / action potentials

neurodsp/sim/cycles.py

@@ -1,13 +1,14 @@
 """Simulating individual cycles, of different types."""
 
 import numpy as np
-from scipy.signal import gaussian, sawtooth
+from scipy.signal import sawtooth
+from scipy.stats import norm
 
 from neurodsp.sim.info import get_sim_func
 from neurodsp.utils.data import compute_nsamples
 from neurodsp.utils.checks import check_param_range, check_param_options
 from neurodsp.utils.decorators import normalize
-from neurodsp.sim.transients import sim_synaptic_kernel
+from neurodsp.sim.transients import sim_synaptic_kernel, sim_action_potential
 
 ###################################################################################################
 ###################################################################################################
@@ -22,15 +23,18 @@ def sim_cycle(n_seconds, fs, cycle_type, phase=0, **cycle_params):
         This is NOT the period of the cycle, but the length of the returned array of the cycle.
     fs : float
         Sampling frequency of the cycle simulation.
-    cycle_type : {'sine', 'asine', 'sawtooth', 'gaussian', 'exp', '2exp'} or callable
+    cycle_type : {'sine', 'asine', 'sawtooth', 'gaussian', 'exp', '2exp', 'ap'} or callable
         What type of cycle to simulate. Options:
 
         * sine: a sine wave cycle
         * asine: an asymmetric sine wave
         * sawtooth: a sawtooth wave
         * gaussian: a gaussian cycle
+        * skewed_gaussian: a skewed gaussian cycle
         * exp: a cycle with exponential decay
         * 2exp: a cycle with exponential rise and decay
+        * ap: an action potential
+
 
     phase : float or {'min', 'max'}, optional, default: 0
         If non-zero, applies a phase shift by rotating the cycle.
@@ -43,8 +47,10 @@ def sim_cycle(n_seconds, fs, cycle_type, phase=0, **cycle_params):
         * asine: `rdsym`, rise-decay symmetry, from 0-1
         * sawtooth: `width`, width of the rising ramp as a proportion of the total cycle
         * gaussian: `std`, standard deviation of the gaussian kernel, in seconds
+        * skewed_gaussian: `center`, `std`, `alpha`, `height`
         * exp: `tau_d`, decay time, in seconds
         * 2exp: `tau_r` & `tau_d` rise time, and decay time, in seconds
+        * ap: `centers`, `stds`, `alphas`, `heights`
 
     Returns
     -------
@@ -219,7 +225,7 @@ def sim_sawtooth_cycle(n_seconds, fs, width):
     return cycle
 
 
-def sim_gaussian_cycle(n_seconds, fs, std):
+def sim_gaussian_cycle(n_seconds, fs, std, center=.5):
     """Simulate a cycle of a gaussian.
 
     Parameters
@@ -230,6 +236,8 @@ def sim_gaussian_cycle(n_seconds, fs, std):
         Sampling frequency of the cycle simulation.
     std : float
         Standard deviation of the gaussian kernel, in seconds.
+    center : float, optional, default: 0.5
+        The center of the gaussian.
 
     Returns
     -------
@@ -243,11 +251,61 @@ def sim_gaussian_cycle(n_seconds, fs, std):
     >>> cycle = sim_gaussian_cycle(n_seconds=0.2, fs=500, std=0.025)
     """
 
-    cycle = gaussian(compute_nsamples(n_seconds, fs), std * fs)
+    xs = np.linspace(0, 1, compute_nsamples(n_seconds, fs))
+    cycle = np.exp(-(xs-center)**2 / (2*std**2))
+
+    return cycle
+
+
+def sim_skewed_gaussian_cycle(n_seconds, fs, center, std, alpha, height=1):
+    """Simulate a cycle of a skewed gaussian.
+
+    Parameters
+    ----------
+    n_seconds : float
+        Length of cycle window in seconds.
+    fs : float
+        Sampling frequency of the cycle simulation.
+    center : float
+        The center of the skewed gaussian.
+    std : float
+        Standard deviation of the gaussian kernel, in seconds.
+    alpha : float
+        Magnitiude and direction of the skew.
+    height : float, optional, default: 1.
+        Maximum value of the cycle.
+
+    Returns
+    -------
+    cycle  : 1d array
+        Output values for skewed gaussian function.
+    """
+
+    n_samples = compute_nsamples(n_seconds, fs)
+
+    # Gaussian distribution
+    cycle = sim_gaussian_cycle(n_seconds, fs, std, center)
+
+    # Skewed cumulative distribution function.
+    #   Assumes time are centered around 0. Adjust to center around non-zero.
+    times = np.linspace(-1, 1, n_samples)
+    cdf = norm.cdf(alpha * ((times - ((center * 2) -1 )) / std))
+
+    # Skew the gaussian
+    cycle = cycle * cdf
+
+    # Rescale height
+    cycle = (cycle / np.max(cycle)) * height
 
     return cycle
 
 
+# Alias action potential from `sim_action_potential`
+def sim_ap_cycle(n_seconds, fs, centers, stds, alphas, heights):
+    return sim_action_potential(n_seconds, fs, centers, stds, alphas, heights)
+sim_ap_cycle.__doc__ = sim_action_potential.__doc__
+
+
 # Alias single exponential cycle from `sim_synaptic_kernel`
 def sim_exp_cycle(n_seconds, fs, tau_d):
     return sim_synaptic_kernel(n_seconds, fs, tau_r=0, tau_d=tau_d)

neurodsp/sim/transients.py

@@ -2,9 +2,11 @@
 
 from warnings import warn
 
+from itertools import repeat
+
 import numpy as np
 
-from neurodsp.utils.data import create_times
+from neurodsp.utils.data import create_times, compute_nsamples
 
 ###################################################################################################
 ###################################################################################################
@@ -71,3 +73,64 @@ def sim_synaptic_kernel(n_seconds, fs, tau_r, tau_d):
     kernel = kernel / np.sum(kernel)
 
     return kernel
+
+
+def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights):
+    """Simulate an action potential as the sum of skewed gaussians.
+
+    Parameters
+    ----------
+    n_seconds : float
+        Length of cycle window in seconds.
+    fs : float
+        Sampling frequency of the cycle simulation.
+    centers : array-like or float
+        Times where the peak occurs in the pre-skewed gaussian.
+    stds : array-like or float
+        Standard deviations of the gaussian kernels, in seconds.
+    alphas : array-like or float
+        Magnitiude and direction of the skew.
+    heights : array-like or float
+        Maximum value of the cycles.
+
+    Returns
+    -------
+    cycle : 1d array
+        Simulated spike cycle.
+    """
+
+    # Prevent circular import
+    from neurodsp.sim.cycles import sim_skewed_gaussian_cycle
+
+    # Determine number of parameters and repeat if necessary
+    params = []
+    n_params = []
+
+    for param in [centers, stds, alphas, heights]:
+
+        if isinstance(param, (tuple, list, np.ndarray)):
+            n_params.append(len(param))
+        else:
+            param = repeat(param)
+
+        params.append(param)
+
+    # Parameter checking
+    if len(n_params) != 0 and len(set(n_params)) != 1:
+        raise ValueError('Unequal lengths between two or more of {centers, stds, alphas, heights}.')
+
+    # Simulate
+    elif len(n_params) == 0:
+        # Single gaussian
+        cycle = sim_skewed_gaussian_cycle(n_seconds, fs, centers, stds, alphas, heights)
+
+    else:
+        # Multiple gaussians
+        cycle = np.zeros((n_params[0], compute_nsamples(n_seconds, fs)))
+
+        for idx, (center, std, alpha, height) in enumerate(zip(*params)):
+            cycle[idx] = sim_skewed_gaussian_cycle(n_seconds, fs,  center, std, alpha,height)
+
+        cycle = np.sum(cycle, axis=0)
+
+    return cycle

neurodsp/tests/sim/test_cycles.py

@@ -95,6 +95,17 @@ def test_sim_gaussian_cycle():
     cycle = sim_gaussian_cycle(N_SECONDS_CYCLE, FS_ODD, 2)
     check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD)
 
+def test_sim_skewed_gaussian_cycle():
+
+    cycle = sim_skewed_gaussian_cycle(N_SECONDS_CYCLE, FS, 0.5, 0.25, 2)
+    check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE)
+
+    cycle = sim_skewed_gaussian_cycle(N_SECONDS_ODD, FS, 0.5, 0.25, 2)
+    check_sim_output(cycle, n_seconds=N_SECONDS_ODD)
+
+    cycle = sim_skewed_gaussian_cycle(N_SECONDS_CYCLE, FS_ODD, 0.5, 0.25, 2)
+    check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD)
+
 def test_create_cycle_time():
 
     times = create_cycle_time(N_SECONDS_CYCLE, FS)

neurodsp/tests/sim/test_transients.py

@@ -3,8 +3,7 @@
 import numpy as np
 
 from neurodsp.tests.tutils import check_sim_output
-from neurodsp.tests.settings import N_SECONDS, FS
-
+from neurodsp.tests.settings import N_SECONDS, N_SECONDS_ODD, FS, FS_ODD
 from neurodsp.sim.transients import *
 
 ###################################################################################################
@@ -18,3 +17,22 @@ def test_sim_synaptic_kernel():
         kernel = sim_synaptic_kernel(N_SECONDS, FS, tau_r, tau_d)
         check_sim_output(kernel)
         assert np.all(kernel >= 0.)
+
+def test_sim_action_potential():
+
+    stds = (.25, .2)
+    alphas = (8, .2)
+    centers = (.25, .5)
+    heights = (15, 2.5)
+
+    cycle = sim_action_potential(N_SECONDS, FS, centers, stds, alphas, heights)
+    check_sim_output(cycle, n_seconds=N_SECONDS)
+
+    cycle = sim_action_potential(N_SECONDS_ODD, FS, centers, stds, alphas, heights)
+    check_sim_output(cycle, n_seconds=N_SECONDS_ODD)
+
+    cycle = sim_action_potential(N_SECONDS, FS_ODD, centers, stds, alphas, heights)
+    check_sim_output(cycle, n_seconds=N_SECONDS, fs=FS_ODD)
+
+    cycle = sim_action_potential(N_SECONDS, FS, centers, stds, alphas, heights)
+    check_sim_output(cycle, n_seconds=N_SECONDS)