Merge pull request #256 from pynapple-org/dev

Changing smooth function

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)
 ------------------
 

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"

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

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"