From a4b1ee80a13ea9100ce3867b27e0e53620123b29 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Sun, 14 Apr 2024 22:47:03 -0400 Subject: [PATCH 1/7] add ac fit funcs --- doc/api.rst | 2 + neurodsp/aperiodic/autocorr.py | 149 ++++++++++++++++++++++ neurodsp/tests/aperiodic/test_autocorr.py | 20 +++ 3 files changed, 171 insertions(+) diff --git a/doc/api.rst b/doc/api.rst index 503478fe..ce45b9b2 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -206,6 +206,8 @@ Auto-correlation Analyses :toctree: generated/ compute_autocorr + compute_decay_time + fit_autocorr Fluctuation Analyses ~~~~~~~~~~~~~~~~~~~~ diff --git a/neurodsp/aperiodic/autocorr.py b/neurodsp/aperiodic/autocorr.py index 42efd625..8844c14f 100644 --- a/neurodsp/aperiodic/autocorr.py +++ b/neurodsp/aperiodic/autocorr.py @@ -1,6 +1,7 @@ """Autocorrelation analyses of time series.""" import numpy as np +from scipy.optimize import curve_fit ################################################################################################### ################################################################################################### @@ -45,3 +46,151 @@ def compute_autocorr(sig, max_lag=1000, lag_step=1, demean=True): timepoints = np.arange(0, max_lag+1, lag_step) return timepoints, autocorrs + + +def compute_decay_time(timepoints, autocorrs, fs, level=0): + """Compute autocorrelation decay time, from precomputed autocorrelation. + + Parameters + ---------- + timepoints : 1d array + Timepoints for the computed autocorrelations. + autocorrs : 1d array + Autocorrelation values. + fs : int + Sampling rate of the signal. + level : float + Autocorrelation decay threshold. + + Returns + ------- + result : float + Autocorrelation decay time. + If decay time value found, returns nan. + + Notes + ----- + The autocorrelation decay time is computed as the time delay for the autocorrelation + to drop to (or below) the decay time threshold. + """ + + val_checks = autocorrs <= level + + if np.any(val_checks): + # Get the first value to cross the threshold, and convert to time value + result = timepoints[np.argmax(val_checks)] / fs + else: + result = np.nan + + return result + + +def fit_autocorr(timepoints, autocorrs, fs=None, fit_function='single_exp'): + """Fit autocorrelation function, returning timescale estimate. + + Parameters + ---------- + timepoints : 1d array + Timepoints for the computed autocorrelations. + autocorrs : 1d array + Autocorrelation values. + fs : int, optional + Sampling rate of the signal. + If provided, timepoints are converted to time values. + fit_func : {'single_exp', 'double_exp'} + Which fitting function to use to fit the autocorrelation results. + + Returns + ------- + popts + Fit parameters. + """ + + if fs: + timepoints = timepoints / fs + + if fit_function == 'single_exp': + p_bounds =([0, 0, 0], [10, np.inf, np.inf]) + elif fit_function == 'double_exp': + p_bounds =([0, 0, 0, 0, 0], [10, np.inf, np.inf, np.inf, np.inf]) + + popts, _ = curve_fit(AC_FIT_FUNCS[fit_function], timepoints, autocorrs, bounds=p_bounds) + + return popts + + +## AC FITTING + +def exp_decay_func(timepoints, tau, scale, offset): + """Exponential decay fit function. + + Parameters + ---------- + timepoints : 1d array + Time values. + tau : float + Timescale value. + scale : float + Scaling factor, which captures the start value of the function. + offset : float + Offset factor, which captures the end value of the function. + + Returns + ------- + ac_fit : 1d array + Results of fitting the function to the autocorrelation. + """ + + return scale * (np.exp(-timepoints / tau) + offset) + + +def double_exp_decay_func(timepoints, tau1, tau2, scale1, scale2, offset): + """Exponential decay fit funtion with two timescales. + + Parameters + ---------- + timepoints : 1d array + Time values. + tau1, tau2 : float + Timescale values, for the 1st and 2nd timescale. + scale1, scale2 : float + Scaling factors, for the 1st and 2nd timescale. + offset : float + Offset factor. + + Returns + ------- + ac_fit : 1d array + Results of fitting the function to the autocorrelation. + """ + + return scale1 * np.exp(-timepoints / tau1) + scale2 * np.exp(-timepoints / tau2) + offset + + +AC_FIT_FUNCS = { + 'single_exp' : exp_decay_func, + 'double_exp' : double_exp_decay_func, +} + + +def compute_ac_fit(timepoints, *popts, fit_function='single_exp'): + """Regenerate values of the exponential decay fit. + + Parameters + ---------- + timepoints : 1d array + Time values. + *popts + Fit parameters. + fit_func : {'single_exp', 'double_exp'} + Which fitting function to use to fit the autocorrelation results. + + Returns + ------- + fit_values : 1d array + Values of the fit to the autocorrelation values. + """ + + fit_func = AC_FIT_FUNCS[fit_function] + + return fit_func(timepoints, *popts) diff --git a/neurodsp/tests/aperiodic/test_autocorr.py b/neurodsp/tests/aperiodic/test_autocorr.py index 04551581..e500cd6f 100644 --- a/neurodsp/tests/aperiodic/test_autocorr.py +++ b/neurodsp/tests/aperiodic/test_autocorr.py @@ -1,5 +1,7 @@ """Tests for neurodsp.aperiodic.autocorr.""" +from neurodsp.tests.settings import FS + from neurodsp.aperiodic.autocorr import * ################################################################################################### @@ -10,3 +12,21 @@ def test_compute_autocorr(tsig): max_lag = 500 timepoints, autocorrs = compute_autocorr(tsig, max_lag=max_lag) assert len(timepoints) == len(autocorrs) == max_lag + 1 + +def test_compute_decay_time(tsig): + + timepoints, autocorrs = compute_autocorr(tsig, max_lag=500) + decay_time = compute_decay_time(timepoints, autocorrs, FS) + assert isinstance(decay_time, float) + +def test_fit_autocorr(tsig): + + timepoints, autocorrs = compute_autocorr(tsig, max_lag=500) + + popts1 = fit_autocorr(timepoints, autocorrs, fit_function='single_exp') + fit_vals1 = compute_ac_fit(timepoints, *popts1, fit_function='single_exp') + assert np.all(fit_vals1) + + popts2 = fit_autocorr(timepoints, autocorrs, fit_function='double_exp') + fit_vals2 = compute_ac_fit(timepoints, *popts2, fit_function='double_exp') + assert np.all(fit_vals2) From f262da8557ba3b4a09b9e7adc30a9853ee9eb30f Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Thu, 18 Jul 2024 21:24:07 -0400 Subject: [PATCH 2/7] minor docstring fixes --- neurodsp/aperiodic/autocorr.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/neurodsp/aperiodic/autocorr.py b/neurodsp/aperiodic/autocorr.py index 85aaf165..cabc3d0c 100644 --- a/neurodsp/aperiodic/autocorr.py +++ b/neurodsp/aperiodic/autocorr.py @@ -1,4 +1,4 @@ -"""Autocorrelation analyses of time series.""" +"""Autocorrelation related analyses of time series.""" import numpy as np from scipy.optimize import curve_fit @@ -28,7 +28,7 @@ def compute_autocorr(sig, max_lag=1000, lag_step=1, demean=True): timepoints : 1d array Time points, in samples, at which autocorrelations are computed. autocorrs : array - Autocorrelation values, for across time lags. + Autocorrelation values, across time lags. Examples -------- @@ -69,12 +69,12 @@ def compute_decay_time(timepoints, autocorrs, fs, level=0): ------- result : float Autocorrelation decay time. - If decay time value found, returns nan. + If decay time value not found, returns nan. Notes ----- - The autocorrelation decay time is computed as the time delay for the autocorrelation - to drop to (or below) the decay time threshold. + The autocorrelation decay time is computed as the time delay for the + autocorrelation to drop to (or below) the decay time threshold. """ val_checks = autocorrs <= level @@ -141,14 +141,14 @@ def exp_decay_func(timepoints, tau, scale, offset): Returns ------- ac_fit : 1d array - Results of fitting the function to the autocorrelation. + Result of fitting the function to the autocorrelation. """ return scale * (np.exp(-timepoints / tau) + offset) def double_exp_decay_func(timepoints, tau1, tau2, scale1, scale2, offset): - """Exponential decay fit funtion with two timescales. + """Exponential decay fit function with two timescales. Parameters ---------- @@ -164,7 +164,7 @@ def double_exp_decay_func(timepoints, tau1, tau2, scale1, scale2, offset): Returns ------- ac_fit : 1d array - Results of fitting the function to the autocorrelation. + Result of fitting the function to the autocorrelation. """ return scale1 * np.exp(-timepoints / tau1) + scale2 * np.exp(-timepoints / tau2) + offset @@ -186,7 +186,7 @@ def compute_ac_fit(timepoints, *popts, fit_function='single_exp'): *popts Fit parameters. fit_func : {'single_exp', 'double_exp'} - Which fitting function to use to fit the autocorrelation results. + Which fit function to use to fit the autocorrelation results. Returns ------- From 4ec8226051d9cf0d40d3254d7ba36b3e33e5fe34 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Sun, 1 Sep 2024 16:40:44 -0400 Subject: [PATCH 3/7] add compute_ac_fit to api list --- doc/api.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/api.rst b/doc/api.rst index e42e1f28..b6438c76 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -209,6 +209,7 @@ Auto-correlation Analyses compute_autocorr compute_decay_time fit_autocorr + compute_ac_fit Fluctuation Analyses ~~~~~~~~~~~~~~~~~~~~ From 1fed3d5f298877047a9b0b23e7961a21d6eb059d Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Sun, 1 Sep 2024 16:41:05 -0400 Subject: [PATCH 4/7] drop optional fs inpiut to fit_ac --- neurodsp/aperiodic/autocorr.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/neurodsp/aperiodic/autocorr.py b/neurodsp/aperiodic/autocorr.py index cabc3d0c..08a8bd2e 100644 --- a/neurodsp/aperiodic/autocorr.py +++ b/neurodsp/aperiodic/autocorr.py @@ -88,13 +88,13 @@ def compute_decay_time(timepoints, autocorrs, fs, level=0): return result -def fit_autocorr(timepoints, autocorrs, fs=None, fit_function='single_exp'): +def fit_autocorr(timepoints, autocorrs, fit_function='single_exp'): """Fit autocorrelation function, returning timescale estimate. Parameters ---------- timepoints : 1d array - Timepoints for the computed autocorrelations. + Timepoints for the computed autocorrelations, in samples or seconds. autocorrs : 1d array Autocorrelation values. fs : int, optional @@ -107,15 +107,17 @@ def fit_autocorr(timepoints, autocorrs, fs=None, fit_function='single_exp'): ------- popts Fit parameters. - """ - if fs: - timepoints = timepoints / fs + Notes + ----- + The values / units of the returned parameters are dependent on the units of samples. + For example, if the timepoints input is in samples, the fit tau value is too. + """ if fit_function == 'single_exp': - p_bounds =([0, 0, 0], [10, np.inf, np.inf]) + p_bounds = ([0, 0, 0], [np.inf, np.inf, np.inf]) elif fit_function == 'double_exp': - p_bounds =([0, 0, 0, 0, 0], [10, np.inf, np.inf, np.inf, np.inf]) + p_bounds = ([0, 0, 0, 0, 0], [np.inf, np.inf, np.inf, np.inf, np.inf]) popts, _ = curve_fit(AC_FIT_FUNCS[fit_function], timepoints, autocorrs, bounds=p_bounds) @@ -182,7 +184,7 @@ def compute_ac_fit(timepoints, *popts, fit_function='single_exp'): Parameters ---------- timepoints : 1d array - Time values. + Time values, in samples or seconds. *popts Fit parameters. fit_func : {'single_exp', 'double_exp'} From 62fc3b060929e677b013821abf35b179eb5b0633 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Sun, 1 Sep 2024 16:45:07 -0400 Subject: [PATCH 5/7] add param input and udpate docstring --- neurodsp/aperiodic/autocorr.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/neurodsp/aperiodic/autocorr.py b/neurodsp/aperiodic/autocorr.py index 08a8bd2e..09773821 100644 --- a/neurodsp/aperiodic/autocorr.py +++ b/neurodsp/aperiodic/autocorr.py @@ -88,7 +88,7 @@ def compute_decay_time(timepoints, autocorrs, fs, level=0): return result -def fit_autocorr(timepoints, autocorrs, fit_function='single_exp'): +def fit_autocorr(timepoints, autocorrs, fit_function='single_exp', bounds=None): """Fit autocorrelation function, returning timescale estimate. Parameters @@ -102,22 +102,30 @@ def fit_autocorr(timepoints, autocorrs, fit_function='single_exp'): If provided, timepoints are converted to time values. fit_func : {'single_exp', 'double_exp'} Which fitting function to use to fit the autocorrelation results. + bounds : tuple of list + Parameter bounds for fitting. + Organized as ([min_p1, min_p1, ...], [max_p1, max_p2, ...]). Returns ------- popts - Fit parameters. + Fit parameters. Parameters depend on the fitting function. + If `fit_func` is 'single_exp', fit parameters are: tau, scale, offset + If `fit_func` is 'douple_exp', fit parameters are: tau1, tau2, scale1, scale2, offset + See fit function for more details. Notes ----- The values / units of the returned parameters are dependent on the units of samples. For example, if the timepoints input is in samples, the fit tau value is too. + If providing parameter bounds, these also need to match the unit of timepoints. """ - if fit_function == 'single_exp': - p_bounds = ([0, 0, 0], [np.inf, np.inf, np.inf]) - elif fit_function == 'double_exp': - p_bounds = ([0, 0, 0, 0, 0], [np.inf, np.inf, np.inf, np.inf, np.inf]) + if not bounds: + if fit_function == 'single_exp': + p_bounds = ([0, 0, 0], [np.inf, np.inf, np.inf]) + elif fit_function == 'double_exp': + p_bounds = ([0, 0, 0, 0, 0], [np.inf, np.inf, np.inf, np.inf, np.inf]) popts, _ = curve_fit(AC_FIT_FUNCS[fit_function], timepoints, autocorrs, bounds=p_bounds) From 6af8c61918b8f7665985c2ae602310e72e904448 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Sun, 1 Sep 2024 16:52:58 -0400 Subject: [PATCH 6/7] fix bounds in fit_ac and add test --- neurodsp/aperiodic/autocorr.py | 6 +++--- neurodsp/tests/aperiodic/test_autocorr.py | 6 ++++++ 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/neurodsp/aperiodic/autocorr.py b/neurodsp/aperiodic/autocorr.py index 09773821..bbde523e 100644 --- a/neurodsp/aperiodic/autocorr.py +++ b/neurodsp/aperiodic/autocorr.py @@ -123,11 +123,11 @@ def fit_autocorr(timepoints, autocorrs, fit_function='single_exp', bounds=None): if not bounds: if fit_function == 'single_exp': - p_bounds = ([0, 0, 0], [np.inf, np.inf, np.inf]) + bounds = ([0, 0, 0], [np.inf, np.inf, np.inf]) elif fit_function == 'double_exp': - p_bounds = ([0, 0, 0, 0, 0], [np.inf, np.inf, np.inf, np.inf, np.inf]) + bounds = ([0, 0, 0, 0, 0], [np.inf, np.inf, np.inf, np.inf, np.inf]) - popts, _ = curve_fit(AC_FIT_FUNCS[fit_function], timepoints, autocorrs, bounds=p_bounds) + popts, _ = curve_fit(AC_FIT_FUNCS[fit_function], timepoints, autocorrs, bounds=bounds) return popts diff --git a/neurodsp/tests/aperiodic/test_autocorr.py b/neurodsp/tests/aperiodic/test_autocorr.py index e500cd6f..3137e1ba 100644 --- a/neurodsp/tests/aperiodic/test_autocorr.py +++ b/neurodsp/tests/aperiodic/test_autocorr.py @@ -27,6 +27,12 @@ def test_fit_autocorr(tsig): fit_vals1 = compute_ac_fit(timepoints, *popts1, fit_function='single_exp') assert np.all(fit_vals1) + # Test with bounds passed in + bounds = ([0, 0, 0], [10, np.inf, np.inf]) + popts1 = fit_autocorr(timepoints, autocorrs, 'single_exp', bounds) + fit_vals1 = compute_ac_fit(timepoints, *popts1, fit_function='single_exp') + assert np.all(fit_vals1) + popts2 = fit_autocorr(timepoints, autocorrs, fit_function='double_exp') fit_vals2 = compute_ac_fit(timepoints, *popts2, fit_function='double_exp') assert np.all(fit_vals2) From 000322b7443a3e279365cadcb1574735936586b1 Mon Sep 17 00:00:00 2001 From: Tom Donoghue Date: Sun, 1 Sep 2024 16:59:16 -0400 Subject: [PATCH 7/7] add some basic accuracy checking to ac fit --- neurodsp/tests/aperiodic/test_autocorr.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/neurodsp/tests/aperiodic/test_autocorr.py b/neurodsp/tests/aperiodic/test_autocorr.py index 3137e1ba..5db31cca 100644 --- a/neurodsp/tests/aperiodic/test_autocorr.py +++ b/neurodsp/tests/aperiodic/test_autocorr.py @@ -20,6 +20,7 @@ def test_compute_decay_time(tsig): assert isinstance(decay_time, float) def test_fit_autocorr(tsig): + # This is a smoke test - check it runs with no accuracy checking timepoints, autocorrs = compute_autocorr(tsig, max_lag=500) @@ -36,3 +37,21 @@ def test_fit_autocorr(tsig): popts2 = fit_autocorr(timepoints, autocorrs, fit_function='double_exp') fit_vals2 = compute_ac_fit(timepoints, *popts2, fit_function='double_exp') assert np.all(fit_vals2) + +def test_fit_autocorr_acc(): + # This test includes some basic accuracy checking + + fs = 100 + tau = 0.015 + lags = np.linspace(0, 100, 1000) + + # Test can fit and recompute the tau value + corrs1 = exp_decay_func(lags, tau, 1, 0) + params1 = fit_autocorr(lags, corrs1) + assert np.isclose(params1[0], tau, 0.001) + + # Test can fit and recompute the tau value - timepoints as time values + lags = lags / fs + corrs2 = exp_decay_func(lags, tau, 1, 0) + params2 = fit_autocorr(lags, corrs2) + assert np.isclose(params2[0], tau, 0.001)