diff --git a/docs/HISTORY.md b/docs/HISTORY.md index 318834d6..02c53ef5 100644 --- a/docs/HISTORY.md +++ b/docs/HISTORY.md @@ -9,6 +9,15 @@ In 2018, Francesco started neuroseries, a Python package built on Pandas. It was In 2021, Guillaume and other trainees in Adrien's lab decided to fork from neuroseries and started *pynapple*. The core of pynapple is largely built upon neuroseries. Some of the original changes to TSToolbox made by Luke were included in this package, especially the *time_support* property of all ts/tsd objects. +0.6.2 (Soon) +------------------ + +- `smooth` now takes standard deviation in time units +- Fixed `TsGroup` saving method. +- `__getattr__` of `BaseTsd` allow numpy functions to be attached as attributes of Tsd objects +- Added `get` method for `TsGroup` + + 0.6.1 (2024-03-03) ------------------ diff --git a/pynapple/core/base_class.py b/pynapple/core/base_class.py index b0ef6a65..cca9f95c 100644 --- a/pynapple/core/base_class.py +++ b/pynapple/core/base_class.py @@ -282,8 +282,8 @@ def count(self, *args, **kwargs): And bincount automatically inherit ep as time support: >>> bincount.time_support - >>> start end - >>> 0 100.0 800.0 + start end + 0 100.0 800.0 """ bin_size = None if "bin_size" in kwargs: @@ -358,8 +358,8 @@ def restrict(self, iset): The time support of newts automatically inherit the epochs defined by ep. >>> newts.time_support - >>> start end - >>> 0 0.0 500.0 + start end + 0 0.0 500.0 """ assert isinstance(iset, IntervalSet), "Argument should be IntervalSet" diff --git a/pynapple/core/time_series.py b/pynapple/core/time_series.py index 5b4c24c4..ffaa33d5 100644 --- a/pynapple/core/time_series.py +++ b/pynapple/core/time_series.py @@ -453,31 +453,71 @@ def convolve(self, array, ep=None, trim="both"): return self.__class__(t=time_array, d=new_data_array, time_support=ep) - def smooth(self, std, size): - """Smooth a time series with a gaussian kernel. std is the standard deviation and size is the number of point of the window. + def smooth(self, std, time_units="s", size_factor=100, norm=True): + """Smooth a time series with a gaussian kernel. - See the scipy documentation : https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.gaussian.html + `std` is the standard deviation of the gaussian kernel in units of time. + If only `std` is passed, the function will compute the standard deviation and size in number + of time points automatically based on the sampling rate of the time series. + For example, if the time series `tsd` has a sample rate of 100 Hz and `std` is 20 ms, + the standard deviation will be converted to an integer through + `tsd.rate * std = int(100 * 0.05) = 5`. + By default, the function will select a kernel size as 100 times the std in number of time points. + `size_factor` can be used to resize the gaussian kernel. + + `norm` set to True normalizes the gaussian kernel to sum to 1. + + In the following example, a time series `tsd` with a sampling rate of 100 Hz + is convolved with a 50 ms non-normalized gaussian kernel. + + >>> tsd.smooth(50, time_units='ms', norm=False) + + This line is equivalent to : + + >>> from scipy.signal.windows import gaussian + >>> kernel = gaussian(M = 500, std=int(tsd.rate*0.05)) + >>> tsd.convolve(window) + + It is generally a good idea to visualize the kernel before applying any convolution. + + See the scipy documentation for the [gaussian window](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.gaussian.html) Parameters ---------- - std : int - Standard deviation - size : int - Description + std : Number + Standard deviation in units of time + time_units : str, optional + The time units in which std is specified ('us', 'ms', 's' [default]). + size_factor : int, optional + How long should be the kernel size as a function of the standard deviation. Default is 100. + norm : bool, optional + Whether to normalized the gaussian kernel or not. Default is `True`. Returns ------- Tsd, TsdFrame, TsdTensor Time series convolved with a gaussian kernel + """ - assert isinstance(std, int), "std should be type int" - assert isinstance(size, int), "size should be type int" - window = signal.windows.gaussian(size, std=std) - window = window / window.sum() + assert isinstance(std, (int, float)), "std should be type int or float" + assert isinstance(size_factor, int), "size_factor should be of type int" + assert isinstance(norm, bool), "norm should be of type boolean" + assert isinstance(time_units, str), "time_units should be of type str" + + std = TsIndex.format_timestamps(np.array([std]), time_units)[0] + + std_size = int(self.rate * std) + M = std_size * size_factor + window = signal.windows.gaussian(M=M, std=std_size) + + if norm: + window = window / window.sum() + return self.convolve(window) def interpolate(self, ts, ep=None, left=None, right=None): - """Wrapper of the numpy linear interpolation method. See https://numpy.org/doc/stable/reference/generated/numpy.interp.html for an explanation of the parameters. + """Wrapper of the numpy linear interpolation method. See [numpy interpolate](https://numpy.org/doc/stable/reference/generated/numpy.interp.html) + for an explanation of the parameters. The argument ts should be Ts, Tsd, TsdFrame, TsdTensor to ensure interpolating from sorted timestamps in the right unit, Parameters diff --git a/tests/test_time_series.py b/tests/test_time_series.py index b437ecf6..d662b666 100755 --- a/tests/test_time_series.py +++ b/tests/test_time_series.py @@ -2,7 +2,7 @@ # @Author: gviejo # @Date: 2022-04-01 09:57:55 # @Last Modified by: Guillaume Viejo -# @Last Modified time: 2024-04-01 15:45:02 +# @Last Modified time: 2024-04-03 10:23:28 #!/usr/bin/env python """Tests of time series for `pynapple` package.""" @@ -454,25 +454,59 @@ def test_convolve(self, tsd): def test_smooth(self, tsd): if not isinstance(tsd, nap.Ts): from scipy import signal - tsd2 = tsd.smooth(2, 20) + tsd2 = tsd.smooth(1) + tmp = tsd.values.reshape(tsd.shape[0], -1) tmp2 = np.zeros_like(tmp) - window = signal.windows.gaussian(20, std=2) + std = int(tsd.rate * 1) + M = std*100 + window = signal.windows.gaussian(M, std=std) window = window / window.sum() for i in range(tmp.shape[-1]): - tmp2[:,i] = np.convolve(tmp[:,i], window, mode='full')[10:-9] + tmp2[:,i] = np.convolve(tmp[:,i], window, mode='full')[M//2:1-M//2] + np.testing.assert_array_almost_equal( + tmp2, + tsd2.values.reshape(tsd2.shape[0], -1) + ) + + tsd2 = tsd.smooth(1000, time_units='ms') + np.testing.assert_array_almost_equal(tmp2, + tsd2.values.reshape(tsd2.shape[0], -1)) + + tsd2 = tsd.smooth(1000000, time_units='us') + np.testing.assert_array_almost_equal(tmp2, + tsd2.values.reshape(tsd2.shape[0], -1)) + + tsd2 = tsd.smooth(1, size_factor=200, norm=False) + tmp = tsd.values.reshape(tsd.shape[0], -1) + tmp2 = np.zeros_like(tmp) + std = int(tsd.rate * 1) + M = std*200 + window = signal.windows.gaussian(M, std=std) + for i in range(tmp.shape[-1]): + tmp2[:,i] = np.convolve(tmp[:,i], window, mode='full')[M//2:1-M//2] np.testing.assert_array_almost_equal( tmp2, tsd2.values.reshape(tsd2.shape[0], -1) ) + def test_smooth_raise_error(self, tsd): + if not isinstance(tsd, nap.Ts): + with pytest.raises(AssertionError) as e_info: + tsd.smooth('a') + assert str(e_info.value) == "std should be type int or float" + + with pytest.raises(AssertionError) as e_info: + tsd.smooth(1, size_factor='b') + assert str(e_info.value) == "size_factor should be of type int" + with pytest.raises(AssertionError) as e_info: - tsd.smooth('a', 20) - assert str(e_info.value) == "std should be type int" + tsd.smooth(1, norm=1) + assert str(e_info.value) == "norm should be of type boolean" with pytest.raises(AssertionError) as e_info: - tsd.smooth(2, 'b') - assert str(e_info.value) == "size should be type int" + tsd.smooth(1, time_units = 0) + assert str(e_info.value) == "time_units should be of type str"