Make Units.waveform_unit configurable

CHANGELOG.md

@@ -2,7 +2,13 @@
 
 ## PyNWB 4.1.1 (Unreleased)
 
+### Changed
+- Added support for NWB Schema 2.10.1
+  - The `unit` attribute of `Units.waveform_mean`, `Units.waveform_sd`, and `Units.waveforms` now has a default value of `"volts"` instead of a fixed value of `"volts"`.
+  - `Units.waveform_mean`, `Units.waveform_sd`, and `Units.waveforms` have a new optional `time_before_peak_in_ms` attribute, exposed as the `waveform_time_before_peak_in_ms` argument and field of `Units`. It holds the time, in milliseconds, from the start of each waveform to the spike peak, i.e., the alignment point used during spike sorting. @rly [#2237](https://github.com/NeurodataWithoutBorders/pynwb/pull/2237)
+
 ### Fixed
+- Fixed `Units.waveform_unit` having no effect on the written file. The `waveform_unit` passed to `Units` is now written to the `waveform_mean`, `waveform_sd`, and `waveforms` columns, and `"volts"` remains the default. PyNWB now also warns when the waveform columns of a file being read carry different `unit` or `sampling_rate` attributes, since only one value per attribute is kept on the `Units` container. @rly [#2162](https://github.com/NeurodataWithoutBorders/pynwb/issues/2162)
 - Fixed `mock_DeviceModel` defaulting `manufacturer` to `None`. The mock now defaults it to `"manufacturer"`. @HugoFara [#2232](https://github.com/NeurodataWithoutBorders/pynwb/pull/2232)
 - Fixed reading a file whose dates carry a sub-minute UTC offset (e.g. `1900-10-01T00:00:00-05:50:36`). @h-mayorquin [#2230](https://github.com/NeurodataWithoutBorders/pynwb/pull/2230)
 - Fixed wide pandas DataFrames in the tutorials spilling out of the content column and into the right margin. @bendichter [#2236](https://github.com/NeurodataWithoutBorders/pynwb/pull/2236)

docs/gallery/domain/ecephys.py

@@ -358,6 +358,24 @@
 # The :py:class:`~pynwb.misc.Units` table can contain simply the spike times of sorted units, or you can also include
 # individual and mean waveform information in some of the optional, predefined :py:class:`~pynwb.misc.Units` table
 # columns: ``waveform_mean``, ``waveform_sd``, or ``waveforms``.
+#
+# The sampling rate and unit of measurement of those three columns are set with the ``waveform_rate`` and
+# ``waveform_unit`` arguments of :py:class:`~pynwb.misc.Units`, which default to ``None`` and ``"volts"``.
+# The ``waveform_time_before_peak_in_ms`` argument records where the spike peak sits within each waveform,
+# that is, the time in milliseconds from the first sample to the alignment point used during spike sorting.
+# Together with ``waveform_rate`` and the number of samples, it locates every waveform sample relative to the
+# spike event. All three are constructor arguments, so set them when you build the
+# :py:class:`~pynwb.misc.Units` table and assign it to :py:attr:`.NWBFile.units`::
+#
+#     from pynwb.misc import Units
+#
+#     nwbfile.units = Units(
+#         name="units",
+#         description="units table",
+#         waveform_rate=30000.0,
+#         waveform_unit="microvolts",
+#         waveform_time_before_peak_in_ms=1.0,
+#     )
 
 nwbfile.units.to_dataframe()
 

src/pynwb/io/core.py

@@ -49,6 +49,8 @@ def get_attr_value(self, **kwargs):
                     return container.parent.waveform_rate
                 if spec.name == 'unit':
                     return container.parent.waveform_unit
+                if spec.name == 'time_before_peak_in_ms':
+                    return container.parent.waveform_time_before_peak_in_ms
             if container.name == 'spike_times':
                 if spec.name == 'resolution':
                     return container.parent.resolution

src/pynwb/io/misc.py

@@ -1,6 +1,9 @@
+import warnings
+
 from hdmf.common.io.table import DynamicTableMap
 
 from .. import register_map
+from .utils import NO_OVERRIDE
 from pynwb.misc import Units
 
 
@@ -10,8 +13,10 @@ class UnitsMap(DynamicTableMap):
     @DynamicTableMap.constructor_arg('resolution')
     def resolution_carg(self, builder, manager):
         if 'spike_times' in builder:
-            return builder['spike_times'].attributes.get('resolution')
-        return None
+            resolution = builder['spike_times'].attributes.get('resolution')
+            if resolution is not None:
+                return resolution
+        return NO_OVERRIDE
 
     @DynamicTableMap.constructor_arg('waveform_rate')
     def waveform_rate_carg(self, builder, manager):
@@ -21,18 +26,33 @@ def waveform_rate_carg(self, builder, manager):
     def waveform_unit_carg(self, builder, manager):
         return self._get_waveform_stat(builder, 'unit')
 
+    @DynamicTableMap.constructor_arg('waveform_time_before_peak_in_ms')
+    def waveform_time_before_peak_in_ms_carg(self, builder, manager):
+        return self._get_waveform_stat(builder, 'time_before_peak_in_ms')
+
     def _get_waveform_stat(self, builder, attribute):
+        """Get the value of an attribute shared by the waveform columns of a Units table.
+
+        The `Units` container holds one `waveform_rate`, one `waveform_unit`, and one
+        `waveform_time_before_peak_in_ms` for the whole table, while the file stores a `sampling_rate`, `unit`,
+        and `time_before_peak_in_ms` attribute on each waveform column. When the columns disagree, the value of
+        the first populated column is used and a warning is raised.
+        """
         waveform_columns = ('waveform_mean', 'waveform_sd', 'waveforms')
-        stats = [builder[column].attributes.get(attribute) for column in waveform_columns if column in builder]
-        if not stats:
-            return None
-        populated_stats = [stat for stat in stats if stat is not None]
-        if len(set(populated_stats)) > 1:
-            # throw warning
-            pass
-        if populated_stats:
-            return populated_stats[0]
-        return None
+        stats = {column: builder[column].attributes.get(attribute)
+                 for column in waveform_columns if column in builder}
+        populated_stats = {column: value for column, value in stats.items() if value is not None}
+        if not populated_stats:
+            return NO_OVERRIDE
+        first_column, first_value = next(iter(populated_stats.items()))
+        if len(set(populated_stats.values())) > 1:
+            warnings.warn(
+                f"The '{attribute}' attribute differs across the waveform columns of Units "
+                f"'{builder.name}': {populated_stats}. Using the value of '{first_column}'.",
+                UserWarning,
+                stacklevel=2
+            )
+        return first_value
 
     @DynamicTableMap.object_attr("electrodes")
     def electrodes_column(self, container, manager):

src/pynwb/misc.py

@@ -154,6 +154,7 @@ class Units(DynamicTable):
     __fields__ = (
         'waveform_rate',
         'waveform_unit',
+        'waveform_time_before_peak_in_ms',
         'resolution'
     )
 
@@ -180,15 +181,24 @@ class Units(DynamicTable):
             {'name': 'electrode_table', 'type': DynamicTable,
              'doc': 'the table that the *electrodes* column indexes', 'default': None},
             {'name': 'waveform_rate', 'type': float,
-             'doc': 'Sampling rate of the waveform means', 'default': None},
+             'doc': 'Sampling rate of the data in the waveform_mean, waveform_sd, and waveforms columns',
+             'default': None},
             {'name': 'waveform_unit', 'type': str,
-             'doc': 'Unit of measurement of the waveform means', 'default': 'volts'},
+             'doc': 'Unit of measurement of the data in the waveform_mean, waveform_sd, and waveforms columns',
+             'default': 'volts'},
+            {'name': 'waveform_time_before_peak_in_ms', 'type': float,
+             'doc': ('Time, in milliseconds, from the start of each waveform in the waveform_mean, waveform_sd, and '
+                     'waveforms columns to the spike peak, i.e., the alignment point used during spike sorting. The '
+                     'same value applies to every unit in the table.'),
+             'default': None},
             {'name': 'resolution', 'type': float,
              'doc': 'The smallest possible difference between two spike times', 'default': None},
             allow_positional=AllowPositional.WARNING,
             )
     def __init__(self, **kwargs):
-        args_to_set = popargs_to_dict(("waveform_rate", "waveform_unit", "resolution"), kwargs)
+        args_to_set = popargs_to_dict(
+            ("waveform_rate", "waveform_unit", "waveform_time_before_peak_in_ms", "resolution"), kwargs
+        )
         electrode_table = popargs("electrode_table", kwargs)
         if kwargs['description'] is None:
             kwargs['description'] = "data on spiking units"

tests/integration/hdf5/test_misc.py

@@ -2,9 +2,10 @@
 import numpy as np
 
 from hdmf.common import VectorData, DynamicTableRegion
-from pynwb import TimeSeries
+from pynwb import NWBHDF5IO, TimeSeries
 from pynwb.misc import Units, DecompositionSeries, FrequencyBandsTable
-from pynwb.testing import NWBH5IOMixin, AcquisitionH5IOMixin, TestCase
+from pynwb.testing import NWBH5IOMixin, AcquisitionH5IOMixin, TestCase, remove_test_file
+from pynwb.testing.mock.file import mock_NWBFile
 from pynwb.ecephys import ElectrodeGroup, ElectrodesTable
 from pynwb.device import Device
 
@@ -123,6 +124,140 @@ def test_waveforms_attributes_written(self):
             self.assertEqual(unit, 'volts')
 
 
+class TestUnitsCustomWaveformUnitIO(AcquisitionH5IOMixin, TestCase):
+    """Test roundtripping a waveform unit other than the default 'volts'."""
+
+    def setUpContainer(self):
+        ut = Units(
+            name='UnitsCustomWaveformUnitTest',
+            description='a simple table for testing a custom Units waveform unit',
+            waveform_rate=40000.,
+            waveform_unit='microvolts',
+        )
+        ut.add_unit(
+            spike_times=[0., 1., 2.],
+            waveform_mean=[1., 2., 3.],
+            waveform_sd=[4., 5., 6.],
+            waveforms=[
+                [  # elec 1
+                    [1, 2, 3],
+                    [1, 2, 3]
+                ], [  # elec 2
+                    [1, 2, 3],
+                    [1, 2, 3]
+                ]
+            ],
+        )
+        return ut
+
+    def test_waveform_unit_roundtrip(self):
+        ut = self.roundtripContainer()
+        self.assertEqual(ut.waveform_unit, 'microvolts')
+
+    def test_waveform_unit_written(self):
+        self.roundtripContainer()
+        with h5py.File(self.filename, 'r') as infile:
+            units = infile['acquisition'][self.container.name]
+            for column in ('waveform_mean', 'waveform_sd', 'waveforms'):
+                unit = units[column].attrs['unit']
+                if isinstance(unit, bytes):
+                    unit = unit.decode('utf-8')
+                self.assertEqual(unit, 'microvolts')
+
+
+class TestUnitsWaveformTimeBeforePeakIO(AcquisitionH5IOMixin, TestCase):
+    """Test roundtripping the time from the start of a waveform to the spike peak."""
+
+    def setUpContainer(self):
+        ut = Units(
+            name='UnitsWaveformTimeBeforePeakTest',
+            description='a simple table for testing the Units waveform peak alignment',
+            waveform_rate=40000.,
+            waveform_time_before_peak_in_ms=1.5,
+        )
+        ut.add_unit(
+            spike_times=[0., 1., 2.],
+            waveform_mean=[1., 2., 3.],
+            waveform_sd=[4., 5., 6.],
+            waveforms=[
+                [  # elec 1
+                    [1, 2, 3],
+                    [1, 2, 3]
+                ], [  # elec 2
+                    [1, 2, 3],
+                    [1, 2, 3]
+                ]
+            ],
+        )
+        return ut
+
+    def test_waveform_time_before_peak_roundtrip(self):
+        ut = self.roundtripContainer()
+        self.assertEqual(ut.waveform_time_before_peak_in_ms, 1.5)
+
+    def test_waveform_time_before_peak_written(self):
+        self.roundtripContainer()
+        with h5py.File(self.filename, 'r') as infile:
+            units = infile['acquisition'][self.container.name]
+            for column in ('waveform_mean', 'waveform_sd', 'waveforms'):
+                self.assertEqual(units[column].attrs['time_before_peak_in_ms'], 1.5)
+
+
+class TestUnitsWaveformTimeBeforePeakOmitted(TestCase):
+    """Test a Units table whose waveform peak alignment is unset."""
+
+    def setUp(self):
+        self.filename = 'test_units_waveform_time_before_peak_omitted.nwb'
+        nwbfile = mock_NWBFile()
+        ut = Units(name='units', description='a table without waveform peak alignment')
+        ut.add_unit(spike_times=[0., 1., 2.], waveform_mean=[1., 2., 3.])
+        nwbfile.units = ut
+        with NWBHDF5IO(self.filename, 'w') as io:
+            io.write(nwbfile)
+
+    def tearDown(self):
+        remove_test_file(self.filename)
+
+    def test_attribute_not_written(self):
+        with h5py.File(self.filename, 'r') as infile:
+            self.assertNotIn('time_before_peak_in_ms', infile['units']['waveform_mean'].attrs)
+
+    def test_read_as_none(self):
+        with NWBHDF5IO(self.filename, 'r') as io:
+            nwbfile = io.read()
+            self.assertIsNone(nwbfile.units.waveform_time_before_peak_in_ms)
+
+
+class TestUnitsMismatchedWaveformUnit(TestCase):
+    """Test reading a file whose waveform columns carry different unit attributes."""
+
+    def setUp(self):
+        self.filename = 'test_units_mismatched_waveform_unit.nwb'
+        nwbfile = mock_NWBFile()
+        ut = Units(name='units', description='a table for testing mismatched waveform units')
+        ut.add_unit(
+            spike_times=[0., 1., 2.],
+            waveform_mean=[1., 2., 3.],
+            waveform_sd=[4., 5., 6.],
+        )
+        nwbfile.units = ut
+        with NWBHDF5IO(self.filename, 'w') as io:
+            io.write(nwbfile)
+        with h5py.File(self.filename, 'r+') as infile:
+            infile['units']['waveform_sd'].attrs['unit'] = 'microvolts'
+
+    def tearDown(self):
+        remove_test_file(self.filename)
+
+    def test_warn_on_mismatched_unit(self):
+        msg = ("The 'unit' attribute differs across the waveform columns of Units 'units': "
+               "{'waveform_mean': 'volts', 'waveform_sd': 'microvolts'}. Using the value of 'waveform_mean'.")
+        with self.assertWarnsWith(UserWarning, msg):
+            with NWBHDF5IO(self.filename, 'r') as io:
+                nwbfile = io.read()
+                self.assertEqual(nwbfile.units.waveform_unit, 'volts')
+
+
 class TestUnitsFileIO(NWBH5IOMixin, TestCase):
 
     def setUpContainer(self):

tests/unit/test_misc.py

@@ -322,6 +322,15 @@ def test_waveform_attrs(self):
         ut = Units(waveform_rate=40000.)
         self.assertEqual(ut.waveform_rate, 40000.)
         self.assertEqual(ut.waveform_unit, 'volts')
+        self.assertIsNone(ut.waveform_time_before_peak_in_ms)
+
+    def test_custom_waveform_unit(self):
+        ut = Units(waveform_unit='microvolts')
+        self.assertEqual(ut.waveform_unit, 'microvolts')
+
+    def test_waveform_time_before_peak_in_ms(self):
+        ut = Units(waveform_time_before_peak_in_ms=1.5)
+        self.assertEqual(ut.waveform_time_before_peak_in_ms, 1.5)
 
     def test_get_starting_time(self):
         """Test get_starting_time returns the earliest spike time across units."""