From 2957868768281a0406fa12aa9ec2caad78a88a07 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Mon, 19 Apr 2021 16:39:22 -0700 Subject: [PATCH 1/6] add skewed_gaussian and ap cycles --- neurodsp/sim/cycles.py | 98 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 94 insertions(+), 4 deletions(-) diff --git a/neurodsp/sim/cycles.py b/neurodsp/sim/cycles.py index 67c2d2b4..c43ef941 100644 --- a/neurodsp/sim/cycles.py +++ b/neurodsp/sim/cycles.py @@ -1,7 +1,8 @@ """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 @@ -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, height=1.): """Simulate a cycle of a gaussian. Parameters @@ -230,6 +236,10 @@ 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. + height : float, optional, default: 1. + Maximum value of the cycle. Returns ------- @@ -243,7 +253,87 @@ 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 = height * 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. + center : float, optional, default: 0.5 + The time where the peak occurs. + height : float, optional, default: 1. + Maximum value of the cycle. + + Returns + ------- + ys : 1d array + Output values for skewed gaussian function. + """ + + # Gaussian distribution + cycle = sim_gaussian_cycle(n_seconds, fs, std/2, center, height) + + # Skewed cumulative distribution function. + # Assumes time are centered around 0. Adjust to center around 0.5. + times = np.linspace(-1, 1, fs) + 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 + + +def sim_ap_cycle(n_seconds, fs, centers, stds, alphas, heights, max_extrema='peak'): + """Simulate an action potential as the sum of two inverse, skewed gaussians. + + Parameters + ---------- + n_seconds : float + Length of cycle window in seconds. + fs : float + Sampling frequency of the cycle simulation. + centers : tuple of (float, float) + The centers of the skewed gaussians. + stds : tuple of (float, float) + Standard deviations of the gaussian kernels, in seconds. + heights : tuple of (float, float) + Maximum value of the cycles. + max_extrema : {'peak', 'trough'} + Defines the larger spike deflection as either peak or trough. + + Returns + ------- + cycle : 1d array + Simulated spike cycle. + """ + + polar = sim_skewed_gaussian_cycle(n_seconds, fs, centers[0], stds[0], + alphas[0], height=heights[0]) + + repolar = sim_skewed_gaussian_cycle(n_seconds, fs, centers[1], stds[1], + alphas[1], height=heights[1]) + + cycle = polar - repolar + + if max_extrema == 'trough': + cycle = -cycle return cycle From 9cb015892b0e200a9ac906d43e8102547455fb1f Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Mon, 19 Apr 2021 16:52:58 -0700 Subject: [PATCH 2/6] tests updated --- neurodsp/sim/cycles.py | 4 +++- neurodsp/tests/sim/test_cycles.py | 31 +++++++++++++++++++++++++++++++ 2 files changed, 34 insertions(+), 1 deletion(-) diff --git a/neurodsp/sim/cycles.py b/neurodsp/sim/cycles.py index c43ef941..a97b463c 100644 --- a/neurodsp/sim/cycles.py +++ b/neurodsp/sim/cycles.py @@ -283,12 +283,14 @@ def sim_skewed_gaussian_cycle(n_seconds, fs, center, std, alpha, height=1): Output values for skewed gaussian function. """ + n_samples = compute_nsamples(n_seconds, fs) + # Gaussian distribution cycle = sim_gaussian_cycle(n_seconds, fs, std/2, center, height) # Skewed cumulative distribution function. # Assumes time are centered around 0. Adjust to center around 0.5. - times = np.linspace(-1, 1, fs) + times = np.linspace(-1, 1, n_samples) cdf = norm.cdf(alpha * ((times - ((center * 2) -1 )) / std)) # Skew the gaussian diff --git a/neurodsp/tests/sim/test_cycles.py b/neurodsp/tests/sim/test_cycles.py index af3b89f6..303f328d 100644 --- a/neurodsp/tests/sim/test_cycles.py +++ b/neurodsp/tests/sim/test_cycles.py @@ -95,6 +95,37 @@ 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_sim_ap_cycle(): + + centers=(.25, .5) + stds=(.25, .2) + alphas=(8, .2) + heights=(15, 2.5) + + cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS, centers, stds, alphas, heights) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) + + cycle = sim_ap_cycle(N_SECONDS_ODD, FS, centers, stds, alphas, heights) + check_sim_output(cycle, n_seconds=N_SECONDS_ODD) + + cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS_ODD, centers, stds, alphas, heights) + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) + + cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS, centers, stds, alphas, heights, 'trough') + check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) + + def test_create_cycle_time(): times = create_cycle_time(N_SECONDS_CYCLE, FS) From 08337cbde8f63d98d2324594cbd08740b07cd451 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Wed, 21 Apr 2021 19:09:28 -0700 Subject: [PATCH 3/6] skewed gaussian std fix --- neurodsp/sim/cycles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neurodsp/sim/cycles.py b/neurodsp/sim/cycles.py index a97b463c..6a8b6913 100644 --- a/neurodsp/sim/cycles.py +++ b/neurodsp/sim/cycles.py @@ -286,7 +286,7 @@ def sim_skewed_gaussian_cycle(n_seconds, fs, center, std, alpha, height=1): n_samples = compute_nsamples(n_seconds, fs) # Gaussian distribution - cycle = sim_gaussian_cycle(n_seconds, fs, std/2, center, height) + cycle = sim_gaussian_cycle(n_seconds, fs, std, center, height) # Skewed cumulative distribution function. # Assumes time are centered around 0. Adjust to center around 0.5. From 489896d8a4e893d51f69b1fe99741dae58ba2bd8 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Mon, 26 Apr 2021 20:21:49 -0700 Subject: [PATCH 4/6] review updates --- neurodsp/sim/cycles.py | 89 ++++---------------------------------- neurodsp/sim/transients.py | 88 +++++++++++++++++++++++++++++++++++++ 2 files changed, 97 insertions(+), 80 deletions(-) diff --git a/neurodsp/sim/cycles.py b/neurodsp/sim/cycles.py index 6a8b6913..453a686d 100644 --- a/neurodsp/sim/cycles.py +++ b/neurodsp/sim/cycles.py @@ -8,7 +8,7 @@ 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_skewed_gaussian, sim_action_potential ################################################################################################### ################################################################################################### @@ -225,7 +225,7 @@ def sim_sawtooth_cycle(n_seconds, fs, width): return cycle -def sim_gaussian_cycle(n_seconds, fs, std, center=.5, height=1.): +def sim_gaussian_cycle(n_seconds, fs, std, center=.5): """Simulate a cycle of a gaussian. Parameters @@ -238,8 +238,6 @@ def sim_gaussian_cycle(n_seconds, fs, std, center=.5, height=1.): Standard deviation of the gaussian kernel, in seconds. center : float, optional, default: 0.5 The center of the gaussian. - height : float, optional, default: 1. - Maximum value of the cycle. Returns ------- @@ -254,90 +252,21 @@ def sim_gaussian_cycle(n_seconds, fs, std, center=.5, height=1.): """ xs = np.linspace(0, 1, compute_nsamples(n_seconds, fs)) - cycle = height * np.exp(-(xs-center)**2 / (2*std**2)) + cycle = np.exp(-(xs-center)**2 / (2*std**2)) return cycle +# Alias skewed gaussian cycle from `sim_skewed_gaussian` 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. - center : float, optional, default: 0.5 - The time where the peak occurs. - height : float, optional, default: 1. - Maximum value of the cycle. - - Returns - ------- - ys : 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, height) - - # Skewed cumulative distribution function. - # Assumes time are centered around 0. Adjust to center around 0.5. - 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 + return sim_skewed_gaussian(n_seconds, fs, center, std, alpha, height) +sim_skewed_gaussian_cycle.__doc__ = sim_skewed_gaussian.__doc__ +# Alias action potential from `sim_action_potential` def sim_ap_cycle(n_seconds, fs, centers, stds, alphas, heights, max_extrema='peak'): - """Simulate an action potential as the sum of two inverse, skewed gaussians. - - Parameters - ---------- - n_seconds : float - Length of cycle window in seconds. - fs : float - Sampling frequency of the cycle simulation. - centers : tuple of (float, float) - The centers of the skewed gaussians. - stds : tuple of (float, float) - Standard deviations of the gaussian kernels, in seconds. - heights : tuple of (float, float) - Maximum value of the cycles. - max_extrema : {'peak', 'trough'} - Defines the larger spike deflection as either peak or trough. - - Returns - ------- - cycle : 1d array - Simulated spike cycle. - """ - - polar = sim_skewed_gaussian_cycle(n_seconds, fs, centers[0], stds[0], - alphas[0], height=heights[0]) - - repolar = sim_skewed_gaussian_cycle(n_seconds, fs, centers[1], stds[1], - alphas[1], height=heights[1]) - - cycle = polar - repolar - - if max_extrema == 'trough': - cycle = -cycle - - return cycle + return sim_action_potential(n_seconds, fs, centers, stds, alphas, heights, max_extrema='peak') +sim_ap_cycle.__doc__ = sim_action_potential.__doc__ # Alias single exponential cycle from `sim_synaptic_kernel` diff --git a/neurodsp/sim/transients.py b/neurodsp/sim/transients.py index d6828e62..c98ed0f2 100644 --- a/neurodsp/sim/transients.py +++ b/neurodsp/sim/transients.py @@ -71,3 +71,91 @@ def sim_synaptic_kernel(n_seconds, fs, tau_r, tau_d): kernel = kernel / np.sum(kernel) return kernel + + +def sim_skewed_gaussian(n_seconds, fs, center, std, alpha, height=1): + """Simulate 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. + center : float, optional, default: 0.5 + Time where the peak occurs in the pre-skewed gaussian. + height : float, optional, default: 1. + Maximum value of the cycle. + + Returns + ------- + ys : 1d array + Output values for skewed gaussian function. + """ + + # Prevent circular import + from neurodsp.sim.cycles import sim_gaussian_cycle + + 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 0.5. + 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 + + +def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights, max_extrema='peak'): + """Simulate an action potential as the sum of two inverse, skewed gaussians. + + Parameters + ---------- + n_seconds : float + Length of cycle window in seconds. + fs : float + Sampling frequency of the cycle simulation. + centers : tuple of (float, float) + Times where the peak occurs in the pre-skewed gaussian. + stds : tuple of (float, float) + Standard deviations of the gaussian kernels, in seconds. + alpha : float + Magnitiude and direction of the skew. + heights : tuple of (float, float) + Maximum value of the cycles. + max_extrema : {'peak', 'trough'} + Defines the larger spike deflection as either peak or trough. + + Returns + ------- + cycle : 1d array + Simulated spike cycle. + """ + + polar = sim_skewed_gaussian(n_seconds, fs, centers[0], stds[0], + alphas[0], height=heights[0]) + + repolar = sim_skewed_gaussian(n_seconds, fs, centers[1], stds[1], + alphas[1], height=heights[1]) + + cycle = polar - repolar + + if max_extrema == 'trough': + cycle = -cycle + + return cycle From f157421e6397f786f48de984b048854524096b35 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Mon, 26 Apr 2021 20:26:45 -0700 Subject: [PATCH 5/6] variable number of ap cycle params --- neurodsp/sim/cycles.py | 5 ++- neurodsp/sim/transients.py | 57 +++++++++++++++++++++++-------- neurodsp/tests/sim/test_cycles.py | 2 +- 3 files changed, 45 insertions(+), 19 deletions(-) diff --git a/neurodsp/sim/cycles.py b/neurodsp/sim/cycles.py index 453a686d..10feda2e 100644 --- a/neurodsp/sim/cycles.py +++ b/neurodsp/sim/cycles.py @@ -2,7 +2,6 @@ import numpy as np 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 @@ -264,8 +263,8 @@ def sim_skewed_gaussian_cycle(n_seconds, fs, center, std, alpha, height=1): # Alias action potential from `sim_action_potential` -def sim_ap_cycle(n_seconds, fs, centers, stds, alphas, heights, max_extrema='peak'): - return sim_action_potential(n_seconds, fs, centers, stds, alphas, heights, max_extrema='peak') +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__ diff --git a/neurodsp/sim/transients.py b/neurodsp/sim/transients.py index c98ed0f2..1bd367f5 100644 --- a/neurodsp/sim/transients.py +++ b/neurodsp/sim/transients.py @@ -3,8 +3,9 @@ from warnings import warn import numpy as np +from scipy.stats import norm -from neurodsp.utils.data import create_times +from neurodsp.utils.data import create_times, compute_nsamples ################################################################################################### ################################################################################################### @@ -121,7 +122,7 @@ def sim_skewed_gaussian(n_seconds, fs, center, std, alpha, height=1): return cycle -def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights, max_extrema='peak'): +def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights): """Simulate an action potential as the sum of two inverse, skewed gaussians. Parameters @@ -130,16 +131,14 @@ def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights, max_extr Length of cycle window in seconds. fs : float Sampling frequency of the cycle simulation. - centers : tuple of (float, float) + centers : array-like or float Times where the peak occurs in the pre-skewed gaussian. - stds : tuple of (float, float) + stds : array-like or float Standard deviations of the gaussian kernels, in seconds. - alpha : float + alphas : array-like or float Magnitiude and direction of the skew. - heights : tuple of (float, float) + heights : array-like or float Maximum value of the cycles. - max_extrema : {'peak', 'trough'} - Defines the larger spike deflection as either peak or trough. Returns ------- @@ -147,15 +146,43 @@ def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights, max_extr Simulated spike cycle. """ - polar = sim_skewed_gaussian(n_seconds, fs, centers[0], stds[0], - alphas[0], height=heights[0]) + # Determine number of parameters + n_params = [] + + if isinstance(centers, (tuple, list, np.ndarray)): + n_params.append(len(centers)) + else: + centers = repeat(centers) + + if isinstance(stds, (tuple, list, np.ndarray)): + n_params.append(len(stds)) + else: + stds = repeat(stds) + + if isinstance(heights, (tuple, list, np.ndarray)): + n_params.append(len(heights)) + else: + heights = repeat(heights) + + # Parameter checking + if len(n_params) == 0: + raise ValueError('Array-like expected for one of {centers, stds, heights}.') + + for param_len in n_params[1:]: + if param_len != n_params[0]: + raise ValueError('Unequal lengths between two or more of {centers, stds, heights}') + + # Initialize cycle array + n_samples = compute_nsamples(n_seconds, fs) + + n_params = n_params[0] - repolar = sim_skewed_gaussian(n_seconds, fs, centers[1], stds[1], - alphas[1], height=heights[1]) + cycle = np.zeros((n_params, n_samples)) - cycle = polar - repolar + # Simulate + for idx, (center, std, alpha, height) in enumerate(zip(centers, stds, alphas, heights)): + cycle[idx] = sim_skewed_gaussian(n_seconds, fs, center, std, alpha, height) - if max_extrema == 'trough': - cycle = -cycle + cycle = np.sum(cycle, axis=0) return cycle diff --git a/neurodsp/tests/sim/test_cycles.py b/neurodsp/tests/sim/test_cycles.py index 303f328d..21366381 100644 --- a/neurodsp/tests/sim/test_cycles.py +++ b/neurodsp/tests/sim/test_cycles.py @@ -122,7 +122,7 @@ def test_sim_ap_cycle(): cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS_ODD, centers, stds, alphas, heights) check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) - cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS, centers, stds, alphas, heights, 'trough') + cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS, centers, stds, alphas, heights) check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) From 86e26578d3c9dd1c1046f6254fc3dd5f2b19ff18 Mon Sep 17 00:00:00 2001 From: ryanhammonds Date: Thu, 27 May 2021 11:36:35 -0700 Subject: [PATCH 6/6] reorganization and cleans of skewed gaussian cycle --- neurodsp/sim/cycles.py | 46 +++++++++++- neurodsp/sim/transients.py | 104 +++++++------------------- neurodsp/tests/sim/test_cycles.py | 20 ----- neurodsp/tests/sim/test_transients.py | 22 +++++- 4 files changed, 88 insertions(+), 104 deletions(-) diff --git a/neurodsp/sim/cycles.py b/neurodsp/sim/cycles.py index 10feda2e..21ec960c 100644 --- a/neurodsp/sim/cycles.py +++ b/neurodsp/sim/cycles.py @@ -2,12 +2,13 @@ import numpy as np 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, sim_skewed_gaussian, sim_action_potential +from neurodsp.sim.transients import sim_synaptic_kernel, sim_action_potential ################################################################################################### ################################################################################################### @@ -256,10 +257,47 @@ def sim_gaussian_cycle(n_seconds, fs, std, center=.5): return cycle -# Alias skewed gaussian cycle from `sim_skewed_gaussian` def sim_skewed_gaussian_cycle(n_seconds, fs, center, std, alpha, height=1): - return sim_skewed_gaussian(n_seconds, fs, center, std, alpha, height) -sim_skewed_gaussian_cycle.__doc__ = sim_skewed_gaussian.__doc__ + """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` diff --git a/neurodsp/sim/transients.py b/neurodsp/sim/transients.py index 1bd367f5..9e725200 100644 --- a/neurodsp/sim/transients.py +++ b/neurodsp/sim/transients.py @@ -2,8 +2,9 @@ from warnings import warn +from itertools import repeat + import numpy as np -from scipy.stats import norm from neurodsp.utils.data import create_times, compute_nsamples @@ -74,56 +75,8 @@ def sim_synaptic_kernel(n_seconds, fs, tau_r, tau_d): return kernel -def sim_skewed_gaussian(n_seconds, fs, center, std, alpha, height=1): - """Simulate 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. - center : float, optional, default: 0.5 - Time where the peak occurs in the pre-skewed gaussian. - height : float, optional, default: 1. - Maximum value of the cycle. - - Returns - ------- - ys : 1d array - Output values for skewed gaussian function. - """ - - # Prevent circular import - from neurodsp.sim.cycles import sim_gaussian_cycle - - 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 0.5. - 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 - - def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights): - """Simulate an action potential as the sum of two inverse, skewed gaussians. + """Simulate an action potential as the sum of skewed gaussians. Parameters ---------- @@ -146,43 +99,38 @@ def sim_action_potential(n_seconds, fs, centers, stds, alphas, heights): Simulated spike cycle. """ - # Determine number of parameters + # Prevent circular import + from neurodsp.sim.cycles import sim_skewed_gaussian_cycle + + # Determine number of parameters and repeat if necessary + params = [] n_params = [] - if isinstance(centers, (tuple, list, np.ndarray)): - n_params.append(len(centers)) - else: - centers = repeat(centers) + for param in [centers, stds, alphas, heights]: - if isinstance(stds, (tuple, list, np.ndarray)): - n_params.append(len(stds)) - else: - stds = repeat(stds) + if isinstance(param, (tuple, list, np.ndarray)): + n_params.append(len(param)) + else: + param = repeat(param) - if isinstance(heights, (tuple, list, np.ndarray)): - n_params.append(len(heights)) - else: - heights = repeat(heights) + params.append(param) # Parameter checking - if len(n_params) == 0: - raise ValueError('Array-like expected for one of {centers, stds, heights}.') - - for param_len in n_params[1:]: - if param_len != n_params[0]: - raise ValueError('Unequal lengths between two or more of {centers, stds, heights}') + if len(n_params) != 0 and len(set(n_params)) != 1: + raise ValueError('Unequal lengths between two or more of {centers, stds, alphas, heights}.') - # Initialize cycle array - n_samples = compute_nsamples(n_seconds, fs) - - n_params = n_params[0] + # Simulate + elif len(n_params) == 0: + # Single gaussian + cycle = sim_skewed_gaussian_cycle(n_seconds, fs, centers, stds, alphas, heights) - cycle = np.zeros((n_params, n_samples)) + else: + # Multiple gaussians + cycle = np.zeros((n_params[0], compute_nsamples(n_seconds, fs))) - # Simulate - for idx, (center, std, alpha, height) in enumerate(zip(centers, stds, alphas, heights)): - cycle[idx] = sim_skewed_gaussian(n_seconds, fs, center, std, alpha, height) + 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) + cycle = np.sum(cycle, axis=0) return cycle diff --git a/neurodsp/tests/sim/test_cycles.py b/neurodsp/tests/sim/test_cycles.py index 21366381..68205613 100644 --- a/neurodsp/tests/sim/test_cycles.py +++ b/neurodsp/tests/sim/test_cycles.py @@ -106,26 +106,6 @@ def test_sim_skewed_gaussian_cycle(): 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_sim_ap_cycle(): - - centers=(.25, .5) - stds=(.25, .2) - alphas=(8, .2) - heights=(15, 2.5) - - cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS, centers, stds, alphas, heights) - check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) - - cycle = sim_ap_cycle(N_SECONDS_ODD, FS, centers, stds, alphas, heights) - check_sim_output(cycle, n_seconds=N_SECONDS_ODD) - - cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS_ODD, centers, stds, alphas, heights) - check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE, fs=FS_ODD) - - cycle = sim_ap_cycle(N_SECONDS_CYCLE, FS, centers, stds, alphas, heights) - check_sim_output(cycle, n_seconds=N_SECONDS_CYCLE) - - def test_create_cycle_time(): times = create_cycle_time(N_SECONDS_CYCLE, FS) diff --git a/neurodsp/tests/sim/test_transients.py b/neurodsp/tests/sim/test_transients.py index eb7cc77f..6a4ea2bd 100644 --- a/neurodsp/tests/sim/test_transients.py +++ b/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)