From 960242f3b756e79cef153ecdf71c3db23831b198 Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Fri, 13 Nov 2026 12:26:26 +0100 Subject: [PATCH 1/3] Fix reading and writing of NWB files not created by Neo NWBIO fails to read most NWB files created by other software: fields that are links to other NWB objects or string datasets are copied into annotations, containers and image series in acquisition are treated as TimeSeries, units without obs_intervals or spike_times and intervals without tags raise errors, and lazy loading fails because the NWB proxies have no parent links. Writing renames every signal and event in place, does not filter nwb: annotations against the target class, ignores offset and channel_conversion, and does not write signal descriptions. The generic BaseTestIO tests never ran for NWBIO because its mode was "one-file". It is now "file", the NeuroConv and tutorial files already in ephy_testing_data are added to the tested entities, and new tests that build NWB files with pynwb run without network access. --- neo/io/nwbio.py | 683 ++++++++++++++++++++++++++-------- neo/test/iotest/test_nwbio.py | 364 +++++++++++++++++- 2 files changed, 881 insertions(+), 166 deletions(-) diff --git a/neo/io/nwbio.py b/neo/io/nwbio.py index 55e5048a3..67362736f 100644 --- a/neo/io/nwbio.py +++ b/neo/io/nwbio.py @@ -18,7 +18,6 @@ import json import logging import os -from collections import defaultdict from itertools import chain from json.decoder import JSONDecodeError @@ -26,6 +25,8 @@ import quantities as pq from neo.core import Segment, SpikeTrain, Epoch, Event, AnalogSignal, IrregularlySampledSignal, Block, ImageSequence +from neo.core.baseneo import _check_annotations +from neo.core.dataobject import ArrayDict from neo.io.baseio import BaseIO from neo.io.proxyobjects import ( AnalogSignalProxy as BaseAnalogSignalProxy, @@ -63,6 +64,10 @@ POSSIBLE_JSON_FIELDS = ("source_script", "description") +# TimeSeries fields that are handled explicitly when writing, so are not taken from `nwb:` annotations. +# The scaling fields are applied to the data when reading, and so must not be applied again. +WRITE_HANDLED_FIELDS = ("conversion", "offset", "channel_conversion", "comments", "description") + prefix_map = {1e9: "giga", 1e6: "mega", 1e3: "kilo", 1: "", 1e-3: "milli", 1e-6: "micro", 1e-9: "nano", 1e-12: "pico"} @@ -75,10 +80,172 @@ def try_json_field(content): """ try: return json.loads(content) - except JSONDecodeError: + except (JSONDecodeError, TypeError): return content +def _parse_json_dict(content): + """ + Return `content` parsed as a JSON object (dict), or None if `content` + is not a string containing a JSON object. + """ + parsed = try_json_field(content) + if isinstance(parsed, dict): + return parsed + return None + + +# Hierarchy information stored by Neo in the `comments` field of a TimeSeries +# always contains these keys +_HIERARCHY_KEYS = ("block", "segment") + +# Maximum number of elements read from a dataset into an annotation +_MAX_ANNOTATION_SIZE = 10000 + + +class _NotAnAnnotation: + """Sentinel type marking a value that cannot be stored as a Neo annotation.""" + + +NOT_AN_ANNOTATION = _NotAnAnnotation() + + +def _to_annotation_value(value, max_size=_MAX_ANNOTATION_SIZE): + """ + Convert a value read from an NWB file to a value that can be stored as a Neo annotation. + + Datasets (h5py, including hdmf's StrDataset) are read into memory, provided they + contain no more than `max_size` elements, and arrays of Python objects + (e.g. variable-length strings) are converted to arrays of str. + + Returns `NOT_AN_ANNOTATION` if the value cannot be stored as an annotation, + e.g. if it is a link to another NWB object (such as a DynamicTableRegion). + """ + from hdmf.container import AbstractContainer + + if isinstance(value, AbstractContainer): + return NOT_AN_ANNOTATION + if hasattr(value, "shape") and hasattr(value, "dtype") and not isinstance(value, (np.ndarray, np.generic)): + # a dataset that has not yet been read from the file + size = int(np.prod(value.shape)) if value.shape else 1 + if size > max_size: + return NOT_AN_ANNOTATION + value = value[()] + if isinstance(value, np.ndarray) and value.dtype.kind == "O": + try: + value = value.astype(str) + except (TypeError, ValueError): + return NOT_AN_ANNOTATION + if isinstance(value, (list, tuple)): + converted = [_to_annotation_value(item, max_size) for item in value] + if any(item is NOT_AN_ANNOTATION for item in converted): + return NOT_AN_ANNOTATION + value = type(value)(converted) + try: + _check_annotations(value) + except (ValueError, TypeError): + return NOT_AN_ANNOTATION + return value + + +def _get_hierarchy(timeseries): + """ + Return the Neo hierarchy information (block and segment names, and optionally the + original Neo object name and the NWB comments) that Neo stores as JSON in the `comments` + field of a TimeSeries, or None if the TimeSeries was not written by Neo. + """ + hierarchy = _parse_json_dict(timeseries.comments) + if hierarchy is not None and all(key in hierarchy for key in _HIERARCHY_KEYS): + return hierarchy + return None + + +def _foreign_comments(timeseries): + """ + Return the `comments` of a TimeSeries, or None if the comments field contains + the default value or Neo hierarchy information. + """ + hierarchy = _get_hierarchy(timeseries) + if hierarchy is not None: + return hierarchy.get("comments", None) + if timeseries.comments in (None, "", "no comments"): + return None + return timeseries.comments + + +def _annotations_equal(a, b): + """Compare two annotation values, which may be NumPy arrays.""" + if isinstance(a, np.ndarray) or isinstance(b, np.ndarray): + return np.array_equal(np.asarray(a), np.asarray(b)) + try: + return bool(a == b) + except ValueError: # e.g. lists containing arrays + return False + + +def _container_fields(container, exclude=()): + """ + Return a dict containing those fields of an NWB container (e.g. an IntracellularElectrode + or a Device) that can be stored as Neo annotations. + """ + fields = {} + for field_name in container.__class__.__nwbfields__: + if field_name in exclude: + continue + value = getattr(container, field_name, None) + if value is None: + continue + value = _to_annotation_value(value) + if value is NOT_AN_ANNOTATION: + logger.debug(f"Unable to store field '{field_name}' of '{container.name}' as an annotation") + else: + fields[field_name] = value + return fields + + +def _ragged_rows(vector_index, rows): + """ + Return the values of the given rows of a ragged column of a DynamicTable, + as a list of arrays, given the VectorIndex for the column. + + We read the index and the data directly, since indexing a VectorIndex with an array of + row numbers does not give the expected result in all versions of hdmf. + """ + ends = np.asarray(vector_index.data[:], dtype=np.int64) + starts = np.concatenate([[0], ends[:-1]]) + data = vector_index.target.data + if len(data) == 0: + return [np.array([], dtype=str) for row in rows] + data = np.asarray(data[:]) + if data.dtype.kind in "OS": + data = data.astype(str) + return [data[starts[row] : ends[row]] for row in rows] + + +def _init_proxy(proxy, name=None, description=None, annotations=None, array_annotations=None, n_channels=None): + """ + Initialise the attributes of an NWB proxy object that would normally be set by + `BaseNeo.__init__()` and `BaseProxy.__init__()`, i.e. name, description, annotations, + array annotations and the links to parent containers. + + The NWB proxies do not call `BaseProxy.__init__()`, since it expects a RawIO. + """ + annotations = annotations or {} + _check_annotations(annotations) + proxy.annotations = annotations + proxy.name = name + proxy.description = description + proxy.file_origin = None + for parent in proxy._parent_containers: + setattr(proxy, parent, None) + if n_channels is not None: + proxy.array_annotations = ArrayDict(n_channels) + if array_annotations: + proxy.array_annotations.update(array_annotations) + else: + proxy.array_annotations = {} + + def get_class(module, name): """ Given a module path and a class name, return the class object @@ -218,7 +385,7 @@ class NWBIO(BaseIO): name = "NeoNWB IO" description = "This IO reads/writes experimental data from/to an .nwb dataset" extensions = ["nwb"] - mode = "one-file" + mode = "file" is_readable = True is_writable = True @@ -253,8 +420,7 @@ def read_all_blocks(self, lazy=False, **kwargs): try: self._file = self._io_nwb.read() except ValueError: - print("Error: Unable to read this version of NWB file.") - print("Please convert to a later NWB format.") + self.logger.error("Unable to read this version of NWB file. Please convert to a later NWB format.") raise self.global_block_metadata = {} @@ -263,6 +429,10 @@ def read_all_blocks(self, lazy=False, **kwargs): if value is not None: if annotation_name in POSSIBLE_JSON_FIELDS: value = try_json_field(value) + value = _to_annotation_value(value) + if value is NOT_AN_ANNOTATION: + self.logger.debug(f"Unable to store NWB file attribute '{annotation_name}' as an annotation") + continue self.global_block_metadata[annotation_name] = value if "session_description" in self.global_block_metadata: self.global_block_metadata["description"] = self.global_block_metadata["session_description"] @@ -278,6 +448,9 @@ def read_all_blocks(self, lazy=False, **kwargs): self._read_units(lazy=lazy) self._read_epochs_group(lazy) + if not self._blocks: + # No data objects could be read, but we still return the file metadata + self._blocks["default"] = Block(name="default", **self.global_block_metadata) return list(self._blocks.values()) def read_block(self, lazy=False, block_index=0, **kargs): @@ -318,6 +491,9 @@ def _read_epochs_group(self, lazy): epoch_names = None if epoch_names is not None: + epoch_names = np.asarray(epoch_names) + segment_names = np.asarray(segment_names) + block_names = np.asarray(block_names) unique_epoch_names = np.unique(epoch_names) for epoch_name in unique_epoch_names: (index,) = np.where((epoch_names == epoch_name)) @@ -342,14 +518,28 @@ def _read_timeseries_group(self, group_name, lazy): group = getattr(self._file, group_name) for timeseries in group.values(): - try: - # NWB files created by Neo store the segment and block names in the comments field - hierarchy = json.loads(timeseries.comments) - except JSONDecodeError: - # For NWB files created with other applications, we put everything in a single - # segment in a single block - # todo: investigate whether there is a reliable way to create multiple segments, - # e.g. using Trial information + if not isinstance(timeseries, pynwb.TimeSeries): + # todo: handle containers such as LFP, Position, BehavioralEvents, + # and neurodata types from extensions + self.logger.info( + f"Skipping '{timeseries.name}' in '{group_name}': " + f"reading {timeseries.__class__.__name__} objects is not yet supported" + ) + continue + if isinstance(timeseries, pynwb.image.ImageSeries) or len(timeseries.data.shape) > 2: + # todo: read ImageSeries as ImageSequence + self.logger.info( + f"Skipping '{timeseries.name}' in '{group_name}': " + f"reading {timeseries.__class__.__name__} objects is not yet supported" + ) + continue + # NWB files created by Neo store the segment and block names in the comments field. + # For NWB files created with other applications, we put everything in a single + # segment in a single block + # todo: investigate whether there is a reliable way to create multiple segments, + # e.g. using Trial information + hierarchy = _get_hierarchy(timeseries) + if hierarchy is None: block_name = "default" segment_name = "default" else: @@ -373,7 +563,18 @@ def _read_timeseries_group(self, group_name, lazy): segment.irregularlysampledsignals.append(signal) def _read_units(self, lazy): - if self._file.units: + if self._file.units is not None and "spike_times" not in self._file.units.colnames: + # e.g. units tables in some ophys files, which hold only cell metadata + self.logger.info("Skipping the units table, since it does not contain spike times") + elif self._file.units: + # Used for units without observation intervals: read the flattened spike times + # of all units once, rather than once per unit. + all_spike_times = self._file.units["spike_times"].target.data[:] + if len(all_spike_times) > 0: + default_t_start = min(0.0, float(np.min(all_spike_times))) + default_t_stop = float(np.max(all_spike_times)) + else: + default_t_start = default_t_stop = 0.0 for id in range(len(self._file.units)): try: # NWB files created by Neo store the segment and block names as extra columns @@ -385,7 +586,9 @@ def _read_units(self, lazy): segment_name = "default" block_name = "default" segment = self._get_segment(block_name, segment_name) - spiketrain = SpikeTrainProxy(self._file.units, id) + spiketrain = SpikeTrainProxy( + self._file.units, id, default_t_start=default_t_start, default_t_stop=default_t_stop + ) if not lazy: spiketrain = spiketrain.load() segment.spiketrains.append(spiketrain) @@ -397,26 +600,24 @@ def _read_stimulus_group(self, lazy): self._read_timeseries_group("stimulus", lazy) def _build_global_annotations(self, blocks): - annotations = defaultdict(set) + annotations = {} for annotation_name in GLOBAL_ANNOTATIONS: if annotation_name in self.annotations: annotations[annotation_name] = self.annotations[annotation_name] else: + # values may be unhashable (e.g. arrays of keywords), so we can't use a set + values = [] for block in blocks: if annotation_name in block.annotations: - try: - annotations[annotation_name].add(block.annotations[annotation_name]) - except TypeError: - if annotation_name in POSSIBLE_JSON_FIELDS: - encoded = json.dumps(block.annotations[annotation_name]) - annotations[annotation_name].add(encoded) - else: - raise - if annotation_name in annotations: - if len(annotations[annotation_name]) > 1: - raise NotImplementedError(f"We don't yet support multiple values for {annotation_name}") - # take single value from set - (annotations[annotation_name],) = annotations[annotation_name] + value = block.annotations[annotation_name] + if annotation_name in POSSIBLE_JSON_FIELDS and isinstance(value, (dict, list)): + value = json.dumps(value) + if not any(_annotations_equal(value, other) for other in values): + values.append(value) + if len(values) > 1: + raise NotImplementedError(f"We don't yet support multiple values for {annotation_name}") + if values: + annotations[annotation_name] = values[0] if "identifier" not in annotations: annotations["identifier"] = str(self.filename) @@ -437,6 +638,11 @@ def write_all_blocks(self, blocks, validate=True, **kwargs): """ import pynwb + if len(blocks) == 0: + raise ValueError("No blocks to write") + # names given to TimeSeries in the NWB file, which must be unique + self._nwb_names = set() + self._dropped_fields = set() global_annotations = self._build_global_annotations(blocks) self._nwbfile = pynwb.NWBFile(**global_annotations) @@ -521,55 +727,147 @@ def _write_electrodes(self, nwbfile, block): def _write_segment(self, nwbfile, segment, electrodes): # maybe use NWB trials to store Segment metadata? + # Note that we do not modify the names of the Neo objects. Where necessary, a different, + # unique name is used in the NWB file, and the original Neo name is stored in the + # hierarchy information in the TimeSeries comments field. for i, signal in enumerate(chain(segment.analogsignals, segment.irregularlysampledsignals)): if signal.segment is not segment: raise TypeError(f"signal.segment must be segment and is {signal.segment}") - if hasattr(signal, "name"): - signal.name = f"{segment.name} {signal.name} {i}" - logging.warning(f"Warning signal name exists. New name: {signal.name}") - else: - signal.name = f"{segment.name} analogsignal{signal.name} {i}" - self._write_signal(self._nwbfile, signal, electrodes) + nwb_name = self._unique_nwb_name(signal.name, f"{segment.name} analogsignal{i}") + self._write_signal(self._nwbfile, signal, electrodes, nwb_name) for i, train in enumerate(segment.spiketrains): if train.segment is not segment: raise TypeError(f"train.segment must be segment and is {train.segment}") - if not train.name: - train.name = f"{segment.name} spiketrain{i}" - self._write_spiketrain(self._nwbfile, train) + # an empty string in the _name column represents a SpikeTrain without a name + self._write_spiketrain(self._nwbfile, train, train.name or "") for i, event in enumerate(segment.events): if event.segment is not segment: - raise TypeError(f"event.segment mst be segment and is {event.segment}") - if hasattr(event, "name"): - event.name = f"{segment.name} {event.name} {i}" - logging.warning(f"Warning event name exists. New name: {event.name}") - else: - event.name = f"{segment.name} event{event.name} {i}" - self._write_event(self._nwbfile, event) + raise TypeError(f"event.segment must be segment and is {event.segment}") + nwb_name = self._unique_nwb_name(event.name, f"{segment.name} event{i}") + self._write_event(self._nwbfile, event, nwb_name) for i, epoch in enumerate(segment.epochs): - if not epoch.name: - epoch.name = f"{segment.name} epoch{i}" - self._write_epoch(self._nwbfile, epoch) + # the _name column is used to group rows of the epochs table into Epochs when reading, + # so Epochs without a name need a unique generated name + self._write_epoch(self._nwbfile, epoch, epoch.name or f"{segment.name} epoch{i}") + + def _unique_nwb_name(self, name, default): + """ + Return a name for a TimeSeries that is unique within the NWB file being written. + + Uses `name` if it is not None and has not already been used, + otherwise `default`, followed by a numerical suffix if necessary. + """ + candidate = name or default + if candidate in self._nwb_names: + base = f"{default} {name}" if name else default + candidate = base + n = 1 + while candidate in self._nwb_names: + candidate = f"{base} ({n})" + n += 1 + self._nwb_names.add(candidate) + return candidate + + @staticmethod + def _accepted_args(timeseries_class): + """Return the names of the arguments accepted by the constructor of an NWB class.""" + from hdmf.utils import get_docval + + return {arg["name"] for arg in get_docval(timeseries_class.__init__)} + + @staticmethod + def _missing_required_args(timeseries_class, provided): + """ + Return the names of the required arguments of the constructor of an NWB class + that are not in `provided`. + """ + from hdmf.utils import get_docval + + return [ + arg["name"] + for arg in get_docval(timeseries_class.__init__) + if "default" not in arg and arg["name"] not in provided + ] + + def _filter_nwb_kwargs(self, timeseries_class, kwargs): + """ + Remove arguments not accepted by the constructor of `timeseries_class`. + + Values taken from `nwb:` annotations may include fields that cannot be set + when constructing the target class (e.g. fields that have fixed values for IZeroClampSeries), + or that belong to a different class. + """ + accepted = self._accepted_args(timeseries_class) + dropped = sorted(set(kwargs) - accepted) + for field_name in dropped: + kwargs.pop(field_name) + key = (timeseries_class.__name__, field_name) + if key not in self._dropped_fields: + self._dropped_fields.add(key) + self.logger.warning( + f"Annotation 'nwb:{field_name}' cannot be written to a {timeseries_class.__name__}, ignoring it" + ) + return kwargs + + def _hierarchy_json(self, obj, segment, nwb_name): + """ + Return the JSON string that Neo stores in the `comments` field of a TimeSeries, + recording the names of the Block and Segment the object belongs to, + the Neo name of the object if this differs from the name used in the NWB file, + and any NWB comments. + """ + hierarchy = {"block": segment.block.name, "segment": segment.name} + if obj.name != nwb_name: + hierarchy["name"] = obj.name + if "nwb:comments" in obj.annotations: + hierarchy["comments"] = obj.annotations["nwb:comments"] + return json.dumps(hierarchy) - def _write_signal(self, nwbfile, signal, electrodes): + def _write_signal(self, nwbfile, signal, electrodes, nwb_name): import pynwb - hierarchy = {"block": signal.segment.block.name, "segment": signal.segment.name} + comments = self._hierarchy_json(signal, signal.segment, nwb_name) + timeseries_class = pynwb.TimeSeries # default if "nwb_neurodata_type" in signal.annotations: - timeseries_class = get_class(*signal.annotations["nwb_neurodata_type"]) - else: - timeseries_class = pynwb.TimeSeries # default - additional_metadata = {name[4:]: value for name, value in signal.annotations.items() if name.startswith("nwb:")} + try: + timeseries_class = get_class(*signal.annotations["nwb_neurodata_type"]) + except (ValueError, AttributeError, ImportError) as err: + self.logger.warning( + f"Unable to find the NWB class {signal.annotations['nwb_neurodata_type']} ({err}), " + f"writing signal '{nwb_name}' as a TimeSeries" + ) + additional_metadata = { + name[4:]: value + for name, value in signal.annotations.items() + if name.startswith("nwb:") and name[4:] not in WRITE_HANDLED_FIELDS + } if "nwb_electrode" in signal.annotations: electrode_name = signal.annotations["nwb_electrode"]["name"] additional_metadata["electrode"] = electrodes[electrode_name] + missing = self._missing_required_args( + timeseries_class, set(additional_metadata) | {"name", "data", "unit", "rate", "starting_time", "timestamps"} + ) + if missing: + # todo: write links to other NWB objects, e.g. the electrodes of an ElectricalSeries + self.logger.warning( + f"Unable to write signal '{nwb_name}' as a {timeseries_class.__name__}, " + f"since the following information is missing: {', '.join(missing)}. Writing it as a TimeSeries." + ) + timeseries_class = pynwb.TimeSeries if timeseries_class != pynwb.TimeSeries: conversion, units = get_units_conversion(signal, timeseries_class) additional_metadata["conversion"] = conversion else: units = signal.units + description = signal.description + if description is None and isinstance(signal.annotations.get("notes", None), dict): + description = json.dumps(signal.annotations["notes"]) + if description is not None: + additional_metadata["description"] = description + additional_metadata = self._filter_nwb_kwargs(timeseries_class, additional_metadata) if hasattr(signal, "proxy_for") and signal.proxy_for in [AnalogSignal, IrregularlySampledSignal]: signal = signal.load() if issubclass(timeseries_class, pynwb.icephys.PatchClampSeries): @@ -582,31 +880,21 @@ def _write_signal(self, nwbfile, signal, electrodes): data = signal.ravel() # convert to 1D else: data = signal + kwargs = dict(name=nwb_name, data=data, unit=units.dimensionality.string, comments=comments) if isinstance(signal, AnalogSignal): - sampling_rate = signal.sampling_rate.rescale("Hz") - tS = timeseries_class( - name=signal.name, - starting_time=time_in_seconds(signal.t_start), - data=data, - unit=units.dimensionality.string, - rate=float(sampling_rate), - comments=json.dumps(hierarchy), - **additional_metadata, - ) + kwargs["starting_time"] = time_in_seconds(signal.t_start) + kwargs["rate"] = float(signal.sampling_rate.rescale("Hz")) # todo: try to add array_annotations via "control" attribute elif isinstance(signal, IrregularlySampledSignal): - tS = timeseries_class( - name=signal.name, - data=data, - unit=units.dimensionality.string, - timestamps=signal.times.rescale("second").magnitude, - comments=json.dumps(hierarchy), - **additional_metadata, - ) + kwargs["timestamps"] = signal.times.rescale("second").magnitude else: raise TypeError( f"signal has type { signal.__class__.__name__}, should be AnalogSignal or IrregularlySampledSignal" ) + if "unit" not in self._accepted_args(timeseries_class): + # some classes have fixed units, e.g. ElectricalSeries (volts) + kwargs.pop("unit") + tS = timeseries_class(**kwargs, **additional_metadata) nwb_group = signal.annotations.get("nwb_group", "acquisition") add_method_map = {"acquisition": self._nwbfile.add_acquisition, "stimulus": self._nwbfile.add_stimulus} if nwb_group in add_method_map: @@ -616,14 +904,18 @@ def _write_signal(self, nwbfile, signal, electrodes): add_time_series(tS) return tS - def _write_spiketrain(self, nwbfile, spiketrain): + def _write_spiketrain(self, nwbfile, spiketrain, name): segment = spiketrain.segment if hasattr(spiketrain, "proxy_for") and spiketrain.proxy_for is SpikeTrain: spiketrain = spiketrain.load() + if "nwb:obs_intervals" in spiketrain.annotations: + obs_intervals = np.asarray(spiketrain.annotations["nwb:obs_intervals"], dtype=float) + else: + obs_intervals = [[float(spiketrain.t_start.rescale("s")), float(spiketrain.t_stop.rescale("s"))]] self._nwbfile.add_unit( spike_times=spiketrain.rescale("s").magnitude, - obs_intervals=[[float(spiketrain.t_start.rescale("s")), float(spiketrain.t_stop.rescale("s"))]], - _name=spiketrain.name, + obs_intervals=obs_intervals, + _name=name, # _description=spiketrain.description, segment=segment.name, block=segment.block.name, @@ -633,36 +925,42 @@ def _write_spiketrain(self, nwbfile, spiketrain): # todo: handle spike waveforms, if any (see SpikeEventSeries) return self._nwbfile.units - def _write_event(self, nwbfile, event): + def _write_event(self, nwbfile, event, nwb_name): import pynwb segment = event.segment + comments = self._hierarchy_json(event, segment, nwb_name) if hasattr(event, "proxy_for") and event.proxy_for == Event: event = event.load() - hierarchy = {"block": segment.block.name, "segment": segment.name} + description = event.description + if description is None and isinstance(event.annotations.get("notes", None), dict): + description = json.dumps(event.annotations["notes"]) tS_evt = pynwb.misc.AnnotationSeries( - name=event.name, + name=nwb_name, data=event.labels, timestamps=event.times.rescale("second").magnitude, - description=event.description or "", - comments=json.dumps(hierarchy), + description=description or "", + comments=comments, ) self._nwbfile.add_acquisition(tS_evt) return tS_evt - def _write_epoch(self, nwbfile, epoch): + def _write_epoch(self, nwbfile, epoch, name): segment = epoch.segment if hasattr(epoch, "proxy_for") and epoch.proxy_for == Epoch: epoch = epoch.load() + labels = epoch.labels + if len(labels) != len(epoch): + # Neo allows Epochs without labels + labels = [""] * len(epoch) for t_start, duration, label in zip( - epoch.rescale("s").magnitude, epoch.durations.rescale("s").magnitude, epoch.labels + epoch.rescale("s").magnitude, epoch.durations.rescale("s").magnitude, labels ): self._nwbfile.add_epoch( - t_start, - t_start + duration, - [label], - [], - _name=epoch.name, + start_time=t_start, + stop_time=t_start + duration, + tags=[label] if label else [], + _name=name, segment=segment.name, block=segment.block.name, ) @@ -687,13 +985,18 @@ class AnalogSignalProxy(BaseAnalogSignalProxy): "timestamps_unit", "electrode", "stream_id", + # scaling fields, which are applied to the data when loading + "conversion", + "offset", + "channel_conversion", ) def __init__(self, timeseries, nwb_group): self._timeseries = timeseries - self.units = timeseries.unit - if timeseries.conversion: - self.units = _recompose_unit(timeseries.unit, timeseries.conversion) + self.shape = self._timeseries.data.shape + if len(self.shape) == 1: + self.shape = (self.shape[0], 1) + self._set_scaling(timeseries) if timeseries.starting_time is not None: self.t_start = timeseries.starting_time * pq.s else: @@ -702,44 +1005,87 @@ def __init__(self, timeseries, nwb_group): self.sampling_rate = timeseries.rate * pq.Hz else: self.sampling_rate = None - self.name = timeseries.name - self.annotations = {"nwb_group": nwb_group} - self.description = try_json_field(timeseries.description) - if isinstance(self.description, dict): - self.annotations["notes"] = self.description - if "name" in self.annotations: - self.annotations.pop("name") - self.description = None - self.shape = self._timeseries.data.shape - if len(self.shape) == 1: - self.shape = (self.shape[0], 1) - metadata_fields = list(timeseries.__nwbfields__) - for field_name in self.__class__.common_metadata_fields: # already handled - try: - metadata_fields.remove(field_name) - except ValueError: - pass + + annotations = {"nwb_group": nwb_group} + description = timeseries.description + notes = _parse_json_dict(description) + if notes is not None: + annotations["notes"] = notes + description = None + comments = _foreign_comments(timeseries) + if comments is not None: + annotations["nwb:comments"] = comments + metadata_fields = [ + field_name + for field_name in timeseries.__nwbfields__ + if field_name not in self.__class__.common_metadata_fields + ] for field_name in metadata_fields: value = getattr(timeseries, field_name) if value is not None: - self.annotations[f"nwb:{field_name}"] = value - self.annotations["nwb_neurodata_type"] = (timeseries.__class__.__module__, timeseries.__class__.__name__) - if hasattr(timeseries, "electrode"): + value = _to_annotation_value(value) + if value is NOT_AN_ANNOTATION: + logger.debug(f"Unable to store field '{field_name}' of '{timeseries.name}' as an annotation") + else: + annotations[f"nwb:{field_name}"] = value + annotations["nwb_neurodata_type"] = (timeseries.__class__.__module__, timeseries.__class__.__name__) + if getattr(timeseries, "electrode", None) is not None: # todo: once the Group class is available, we could add electrode metadata # to a Group containing all signals that share that electrode # This would reduce the amount of redundancy (repeated metadata in every signal) - electrode_metadata = {"device": {}} - metadata_fields = list(timeseries.electrode.__class__.__nwbfields__) + ["name"] - metadata_fields.remove("device") # needs special handling - for field_name in metadata_fields: - value = getattr(timeseries.electrode, field_name) - if value is not None: - electrode_metadata[field_name] = value - for field_name in timeseries.electrode.device.__class__.__nwbfields__: - value = getattr(timeseries.electrode.device, field_name) - if value is not None: - electrode_metadata["device"][field_name] = value - self.annotations["nwb_electrode"] = electrode_metadata + electrode = timeseries.electrode + electrode_metadata = _container_fields(electrode, exclude=("device",)) + electrode_metadata["name"] = electrode.name + electrode_metadata["device"] = _container_fields(electrode.device) + annotations["nwb_electrode"] = electrode_metadata + + name = timeseries.name + hierarchy = _get_hierarchy(timeseries) + if hierarchy is not None and "name" in hierarchy: + name = hierarchy["name"] + _init_proxy(self, name=name, description=description, annotations=annotations, n_channels=self.shape[-1]) + + def _set_scaling(self, timeseries): + """ + Determine the units of the signal, and the gain and offset to be applied to the data + when loading, from the `unit`, `conversion`, `offset` and (for ElectricalSeries) + `channel_conversion` fields of the TimeSeries. + + In NWB, the value in `unit` is data * channel_conversion * conversion + offset. + Where possible, the conversion factor is folded into the units (e.g. a conversion + of 1e-3 with unit "volts" gives units of mV), so that no scaling of the data is needed. + """ + conversion = timeseries.conversion if timeseries.conversion else 1.0 + offset = getattr(timeseries, "offset", None) or 0.0 + try: + self.units = _recompose_unit(timeseries.unit, conversion) + gain = 1.0 + unit_conversion = conversion + except ValueError: + # conversion factor is not a power of ten, so we apply it to the data + self.units = _recompose_unit(timeseries.unit, 1.0) + gain = conversion + unit_conversion = 1.0 + channel_conversion = getattr(timeseries, "channel_conversion", None) + if channel_conversion is not None: + channel_conversion = np.asarray(channel_conversion[:], dtype=float) + if channel_conversion.shape == (self.shape[-1],): + gain = gain * channel_conversion + else: + logger.warning( + f"Ignoring channel_conversion of '{timeseries.name}': it has {channel_conversion.size} values " + f"but the data have {self.shape[-1]} channels" + ) + self._gain = gain + self._offset = offset / unit_conversion + + def _scale(self, data): + """Apply the gain and offset determined by `_set_scaling()` to the data.""" + if np.all(np.asarray(self._gain) == 1.0) and self._offset == 0.0: + return data + if data.ndim == 1: + data = data.reshape(-1, 1) + return data * self._gain + self._offset def load(self, time_slice=None, strict_slicing=True): """ @@ -756,7 +1102,7 @@ def load(self, time_slice=None, strict_slicing=True): i_start, i_stop = np.searchsorted(self._timeseries.timestamps, time_slice) else: i_start, i_stop, sig_t_start = self._time_slice_indices(time_slice, strict_slicing=strict_slicing) - signal = self._timeseries.data[i_start:i_stop] + signal = self._scale(self._timeseries.data[i_start:i_stop]) if self.sampling_rate is None: return IrregularlySampledSignal( self._timeseries.timestamps[i_start:i_stop] * pq.s, @@ -787,13 +1133,21 @@ class EventProxy(BaseEventProxy): def __init__(self, timeseries, nwb_group): self._timeseries = timeseries - self.name = timeseries.name - self.annotations = {"nwb_group": nwb_group} - self.description = try_json_field(timeseries.description) - if isinstance(self.description, dict): - self.annotations.update(self.description) - self.description = None self.shape = self._timeseries.data.shape + annotations = {"nwb_group": nwb_group} + description = timeseries.description + notes = _parse_json_dict(description) + if notes is not None: + annotations["notes"] = notes + description = None + comments = _foreign_comments(timeseries) + if comments is not None: + annotations["nwb:comments"] = comments + name = timeseries.name + hierarchy = _get_hierarchy(timeseries) + if hierarchy is not None and "name" in hierarchy: + name = hierarchy["name"] + _init_proxy(self, name=name, description=description, annotations=annotations, n_channels=self.shape[-1]) def load(self, time_slice=None, strict_slicing=True): """ @@ -821,17 +1175,23 @@ def __init__(self, time_intervals, epoch_name=None, index=None): :param epoch_name: (str) Name of the epoch object :param index: (np.array, slice) - Slice object or array of bool values masking time_intervals to be used. In case of - an array it has to have the same shape as `time_intervals`. + Array of integer indices or of bool values, or slice object, selecting the rows of + time_intervals to be used. In the case of an array of bool values, + it has to have the same shape as `time_intervals`. """ self._time_intervals = time_intervals if index is not None: self._index = index - self.shape = (index.sum(),) + if isinstance(index, slice): + self.shape = (len(range(*index.indices(len(time_intervals)))),) + elif np.asarray(index).dtype == bool: + self.shape = (int(np.sum(index)),) + else: + self.shape = (len(index),) else: self._index = slice(None) self.shape = (len(time_intervals),) - self.name = epoch_name + _init_proxy(self, name=epoch_name, n_channels=self.shape[-1]) def load(self, time_slice=None, strict_slicing=True): """ @@ -845,41 +1205,53 @@ def load(self, time_slice=None, strict_slicing=True): if time_slice: raise NotImplementedError("todo") else: - start_times = self._time_intervals.start_time[self._index] - stop_times = self._time_intervals.stop_time[self._index] + start_times = np.asarray(self._time_intervals.start_time[self._index]) + stop_times = np.asarray(self._time_intervals.stop_time[self._index]) durations = stop_times - start_times - labels = self._time_intervals.tags[self._index] + if "tags" in self._time_intervals.colnames: + # each row may have zero, one or several tags; Neo has a single label per epoch + rows = np.arange(len(self._time_intervals))[self._index] + tags = _ragged_rows(self._time_intervals["tags"], rows) + labels = np.array([", ".join(str(tag) for tag in row_tags) for row_tags in tags], dtype=str) + else: + labels = None return Epoch(times=start_times * pq.s, durations=durations * pq.s, labels=labels, name=self.name) class SpikeTrainProxy(BaseSpikeTrainProxy): - def __init__(self, units_table, id): + def __init__(self, units_table, id, default_t_start=0.0, default_t_stop=0.0): """ :param units_table: A Units table (see https://pynwb.readthedocs.io/en/stable/pynwb.misc.html#pynwb.misc.Units) :param id: the cell/unit ID (integer) + :param default_t_start: start time (in seconds) to use if the unit has no observation intervals + :param default_t_stop: stop time (in seconds) to use if the unit has no observation intervals """ self._units_table = units_table self.id = id self.units = pq.s - obs_intervals = units_table.get_unit_obs_intervals(id) + annotations = {"nwb_group": "acquisition"} + if "obs_intervals" in units_table.colnames: + obs_intervals = np.asarray(units_table.get_unit_obs_intervals(id), dtype=float) + else: + obs_intervals = np.empty((0, 2)) if len(obs_intervals) == 0: - t_start, t_stop = None, None - elif len(obs_intervals) == 1: - t_start, t_stop = obs_intervals[0] + t_start, t_stop = default_t_start, default_t_stop else: - raise NotImplementedError("Can't yet handle multiple observation intervals") + t_start, t_stop = obs_intervals[:, 0].min(), obs_intervals[:, 1].max() + if len(obs_intervals) > 1: + annotations["nwb:obs_intervals"] = obs_intervals self.t_start = t_start * pq.s self.t_stop = t_stop * pq.s - self.annotations = {"nwb_group": "acquisition"} - try: + if "_name" in units_table.colnames: # NWB files created by Neo store the name as an extra column - self.name = units_table._name[id] - except AttributeError: - self.name = None + name = units_table["_name"][id] or None + else: + name = None self.shape = None # no way to get this without reading the data + _init_proxy(self, name=name, annotations=annotations) def load(self, time_slice=None, strict_slicing=True): """ @@ -891,15 +1263,24 @@ def load(self, time_slice=None, strict_slicing=True): (t_start or t_stop) is outside the real time range of the segment. """ interval = None + t_start, t_stop = self.t_start, self.t_stop if time_slice: - interval = (float(t) for t in time_slice) # convert from quantities - spike_times = self._units_table.get_unit_spike_times(self.id, in_interval=interval) + interval = tuple(float(t.rescale(pq.s)) for t in time_slice) + t_start = max(t_start, time_slice[0]) + t_stop = min(t_stop, time_slice[1]) + spike_times = np.asarray(self._units_table.get_unit_spike_times(self.id, in_interval=interval)) + if spike_times.size > 0: + # some files contain spikes outside the observation intervals + if spike_times.min() < float(t_start.rescale(pq.s)) or spike_times.max() > float(t_stop.rescale(pq.s)): + logger.debug(f"Spike times of unit {self.id} extend beyond its observation intervals") + t_start = min(t_start, spike_times.min() * pq.s) + t_stop = max(t_stop, spike_times.max() * pq.s) return SpikeTrain( spike_times * self.units, - self.t_stop, + t_stop, units=self.units, # sampling_rate=array(1.) * Hz, - t_start=self.t_start, + t_start=t_start, # waveforms=None, # left_sweep=None, name=self.name, diff --git a/neo/test/iotest/test_nwbio.py b/neo/test/iotest/test_nwbio.py index 6a96fbf3d..c950b1bc0 100644 --- a/neo/test/iotest/test_nwbio.py +++ b/neo/test/iotest/test_nwbio.py @@ -3,14 +3,12 @@ Tests of neo.io.nwbio """ -import os +import json +import tempfile import unittest -from datetime import datetime +from datetime import datetime, timezone +from pathlib import Path -try: - from urllib.request import urlretrieve -except ImportError: - from urllib import urlretrieve from neo.test.iotest.common_io_test import BaseTestIO from neo.core import AnalogSignal, SpikeTrain, Event, Epoch, IrregularlySampledSignal, Segment, Block @@ -32,12 +30,41 @@ @unittest.skipUnless(HAVE_PYNWB, "requires pynwb") class TestNWBIO(BaseTestIO, unittest.TestCase): + """Generic compliance tests, using NWB files from the ephy_testing_data repository.""" + ioclass = NWBIO entities_to_download = ["nwb"] entities_to_test = [ # Files from Allen Institute: "nwb/H19.29.141.11.21.01.nwb", # 7 MB + # File generated by the PyNWB ecephys tutorial + "nwb/ecephys_tutorial_v2.5.0.nwb", + # Files generated by NeuroConv + "nwb/AbfInterface_3.nwb", + "nwb/AxonaLFPDataInterface.nwb", + "nwb/AxonaRecordingInterface.nwb", + "nwb/BlackrockSortingInterface.nwb", ] + # The generic write tests use the same io_kwargs for reading and writing, whereas NWBIO needs + # mode="w" for writing, and the generated Blocks lack the 'session_start_time' that NWB requires. + # Writing is tested in TestNWBIOInMemory and TestNWBIOForeignFiles instead. + read_and_write_is_bijective = False + + @classmethod + def generate_files_for_io_able_to_write(cls): + cls.files_generated = [] + + +@unittest.skipUnless(HAVE_PYNWB, "requires pynwb") +class TestNWBIOInMemory(unittest.TestCase): + """Tests that create their own NWB files, so do not need downloaded test data.""" + + def setUp(self): + self._tmpdir = tempfile.TemporaryDirectory() + self.tmpdir = Path(self._tmpdir.name) + + def tearDown(self): + self._tmpdir.cleanup() def test_roundtrip(self): @@ -115,7 +142,7 @@ def test_roundtrip(self): seg.events.append(evt) # write to file - test_file_name = "test_round_trip.nwb" + test_file_name = self.tmpdir / "test_round_trip.nwb" iow = NWBIO(filename=test_file_name, mode="w") iow.write_all_blocks(original_blocks) @@ -166,16 +193,25 @@ def test_roundtrip(self): original_epoch_11 = original_blocks[1].segments[1].epochs[0] retrieved_epoch_11 = retrieved_blocks[1].segments[1].epochs[0] - for attr_name in ("name",): - retrieved_attribute = getattr(retrieved_epoch_11, attr_name) - original_attribute = getattr(original_epoch_11, attr_name) - self.assertEqual(retrieved_attribute, original_attribute) + # Epochs without names are given a generated name in the NWB file, + # but the original Epoch should not be modified + self.assertIsNone(original_epoch_11.name) + self.assertEqual(retrieved_epoch_11.name, "Second block segment1 epoch0") assert_array_equal(retrieved_epoch_11.rescale("ms").magnitude, original_epoch_11.rescale("ms").magnitude) assert_allclose( retrieved_epoch_11.durations.rescale("ms").magnitude, original_epoch_11.durations.rescale("ms").magnitude ) assert_array_equal(retrieved_epoch_11.labels, original_epoch_11.labels) - os.remove(test_file_name) + + # writing should not modify the names of the original objects, + # and the original names should be restored when reading + self.assertEqual(original_blocks[0].segments[0].analogsignals[0].name, "Signal_a None") + for original_block, retrieved_block in zip(original_blocks, retrieved_blocks): + for original_seg, retrieved_seg in zip(original_block.segments, retrieved_block.segments): + self.assertEqual( + [sig.name for sig in original_seg.analogsignals], [sig.name for sig in retrieved_seg.analogsignals] + ) + self.assertEqual([ev.name for ev in original_seg.events], [ev.name for ev in retrieved_seg.events]) def test_roundtrip_with_annotations(self): # test with NWB-specific annotations @@ -223,7 +259,7 @@ def test_roundtrip_with_annotations(self): segment.analogsignals = [stimulus, response] assert stimulus.segment is response.segment is segment - test_file_name = "test_round_trip_with_annotations.nwb" + test_file_name = self.tmpdir / "test_round_trip_with_annotations.nwb" iow = NWBIO(filename=test_file_name, mode="w") iow.write_all_blocks([original_block]) @@ -244,10 +280,21 @@ def test_roundtrip_with_annotations(self): self.assertEqual(retrieved_attribute, original_attribute) assert_array_equal(retrieved_response.magnitude, original_response.magnitude) - os.remove(test_file_name) + # NWB-specific annotations should be preserved + for original_signal in (stimulus, response): + retrieved_signal = retrieved_block.segments[0].filter(name=original_signal.name)[0] + for key, value in original_signal.annotations.items(): + if key == "nwb_electrode": + for elec_key, elec_value in value.items(): + if elec_key == "device": + self.assertEqual(retrieved_signal.annotations[key]["device"]["name"], elec_value["name"]) + else: + self.assertEqual(retrieved_signal.annotations[key][elec_key], elec_value) + else: + self.assertEqual(retrieved_signal.annotations[key], value) def test_write_proxy_objects(self): - test_file_name = self.local_test_dir / "test_round_trip_with_annotations.nwb" + test_file_name = self.tmpdir / "test_write_proxy_objects.nwb" # generate dummy IO as basis for ProxyObjects self.proxy_reader = ExampleRawIO(filename="my_filename.fake") @@ -304,6 +351,293 @@ def test_write_proxy_objects(self): assert_array_equal(original_segment.epochs[0].load().magnitude, retrieved_segment.epochs[0].magnitude) +def _make_nwbfile(): + """Create a minimal NWBFile, as it might be created by software other than Neo.""" + return pynwb.NWBFile( + session_description="test file", + identifier="test-file", + session_start_time=datetime(2020, 1, 1, tzinfo=timezone.utc), + ) + + +def _add_electrodes(nwbfile, n_electrodes): + """Add electrodes to an NWBFile, and return a region covering all of them.""" + device = nwbfile.create_device(name="probe") + group = nwbfile.create_electrode_group(name="shank0", description="shank", location="CA1", device=device) + for i in range(n_electrodes): + nwbfile.add_electrode(group=group, location="CA1") + return nwbfile.create_electrode_table_region(list(range(n_electrodes)), "all electrodes") + + +@unittest.skipUnless(HAVE_PYNWB, "requires pynwb") +class TestNWBIOForeignFiles(unittest.TestCase): + """Tests reading (and re-writing) NWB files that were not created by Neo.""" + + def setUp(self): + self._tmpdir = tempfile.TemporaryDirectory() + self.tmpdir = Path(self._tmpdir.name) + + def tearDown(self): + self._tmpdir.cleanup() + + def _write(self, nwbfile, name="foreign.nwb"): + path = self.tmpdir / name + with pynwb.NWBHDF5IO(path, mode="w") as io: + io.write(nwbfile) + return path + + def _read(self, path, lazy=False): + blocks = NWBIO(path, mode="r").read_all_blocks(lazy=lazy) + if lazy: + for block in blocks: + for segment in block.segments: + for proxy in segment.data_children: + self.assertIs(proxy.segment, segment) + return blocks + + def test_read_lazy(self): + nwbfile = _make_nwbfile() + nwbfile.add_acquisition( + pynwb.TimeSeries(name="ts", data=np.arange(10.0), unit="volts", rate=1000.0, starting_time=0.0) + ) + path = self._write(nwbfile) + segment = self._read(path, lazy=True)[0].segments[0] + self.assertEqual(len(segment.analogsignals), 1) + signal = segment.analogsignals[0].load() + assert_array_equal(signal.magnitude.ravel(), np.arange(10.0)) + + def test_fields_that_cannot_be_annotations(self): + # The `electrodes` field of an ElectricalSeries is a link to rows of the electrodes table, + # and string arrays are read as StrDataset; neither can be stored directly as an annotation + nwbfile = _make_nwbfile() + region = _add_electrodes(nwbfile, 4) + nwbfile.keywords = ["hippocampus", "LFP"] + nwbfile.add_acquisition( + pynwb.ecephys.ElectricalSeries(name="es", data=np.ones((20, 4)), electrodes=region, rate=1000.0) + ) + path = self._write(nwbfile) + for lazy in (False, True): + block = self._read(path, lazy=lazy)[0] + assert_array_equal(block.annotations["keywords"], ["hippocampus", "LFP"]) + signal = block.segments[0].analogsignals[0] + if lazy: + signal = signal.load() + self.assertEqual(signal.shape, (20, 4)) + self.assertNotIn("nwb:electrodes", signal.annotations) + + # writing back: we can't yet write the electrode information, + # so the signal should be written as a generic TimeSeries + out_path = self.tmpdir / "rewritten.nwb" + with self.assertLogs("neo.io.nwbio.NWBIO", level="WARNING") as cm: + NWBIO(out_path, mode="w").write_all_blocks([block]) + self.assertIn("ElectricalSeries", cm.output[0]) + with pynwb.NWBHDF5IO(out_path, mode="r") as io: + rewritten = io.read() + self.assertEqual(list(rewritten.keywords), ["hippocampus", "LFP"]) + self.assertEqual(type(rewritten.acquisition["es"]), pynwb.TimeSeries) + + def test_scaling(self): + nwbfile = _make_nwbfile() + raw = np.arange(40, dtype=np.int16).reshape(10, 4) + # conversion factor that is not a power of ten, as used by e.g. Intan + nwbfile.add_acquisition( + pynwb.TimeSeries( + name="not_power_of_ten", data=raw, unit="volts", rate=1000.0, conversion=0.195e-6, offset=-1e-3 + ) + ) + # conversion factor that can be expressed through the units + nwbfile.add_acquisition( + pynwb.TimeSeries(name="power_of_ten", data=raw, unit="volts", rate=1000.0, conversion=1e-3, offset=0.5) + ) + region = _add_electrodes(nwbfile, 4) + nwbfile.add_acquisition( + pynwb.ecephys.ElectricalSeries( + name="channel_conversion", + data=raw, + electrodes=region, + rate=1000.0, + conversion=1e-6, + channel_conversion=[1.0, 2.0, 3.0, 4.0], + ) + ) + path = self._write(nwbfile) + with pynwb.NWBHDF5IO(path, mode="r") as io: + nwbfile_in = io.read() + expected = {name: ts.get_data_in_units() for name, ts in nwbfile_in.acquisition.items()} + for lazy in (False, True): + segment = self._read(path, lazy=lazy)[0].segments[0] + for signal in segment.analogsignals: + if lazy: + signal = signal.load() + assert_allclose(signal.rescale("V").magnitude, expected[signal.name]) + for key in ("nwb:conversion", "nwb:offset", "nwb:channel_conversion"): + self.assertNotIn(key, signal.annotations) + + # the scaling should not be applied a second time when writing + block = self._read(path)[0] + out_path = self.tmpdir / "rewritten.nwb" + NWBIO(out_path, mode="w").write_all_blocks([block]) + segment = self._read(out_path)[0].segments[0] + for signal in segment.analogsignals: + assert_allclose(signal.rescale("V").magnitude, expected[signal.name]) + + def test_skip_unsupported_objects(self): + # containers such as Position, and image data, are not yet supported, + # but should not prevent the rest of the file being read + nwbfile = _make_nwbfile() + position = pynwb.behavior.Position(name="Position") + position.create_spatial_series( + name="pos", data=np.ones((10, 2)), reference_frame="origin", rate=10.0, starting_time=0.0 + ) + nwbfile.add_acquisition(position) + nwbfile.add_acquisition( + pynwb.image.ImageSeries(name="images", data=np.zeros((5, 4, 4)), unit="n.a.", rate=10.0) + ) + nwbfile.add_acquisition(pynwb.TimeSeries(name="ts", data=np.arange(10.0), unit="volts", rate=10.0)) + path = self._write(nwbfile) + for lazy in (False, True): + segment = self._read(path, lazy=lazy)[0].segments[0] + self.assertEqual([sig.name for sig in segment.analogsignals], ["ts"]) + + def test_units_without_obs_intervals(self): + nwbfile = _make_nwbfile() + nwbfile.add_unit(spike_times=[0.5, 1.5, 2.5]) + nwbfile.add_unit(spike_times=[-0.5, 3.5]) + path = self._write(nwbfile) + for lazy in (False, True): + spiketrains = self._read(path, lazy=lazy)[0].segments[0].spiketrains + if lazy: + spiketrains = [st.load() for st in spiketrains] + assert_array_equal(spiketrains[0].magnitude, [0.5, 1.5, 2.5]) + for st in spiketrains: + self.assertEqual(st.t_start, -0.5 * pq.s) + self.assertEqual(st.t_stop, 3.5 * pq.s) + + def test_units_with_multiple_obs_intervals(self): + nwbfile = _make_nwbfile() + nwbfile.add_unit(spike_times=[0.5, 1.5, 4.5], obs_intervals=[[0.0, 2.0], [4.0, 5.0]]) + # spikes outside the observation interval + nwbfile.add_unit(spike_times=[0.5, 6.0], obs_intervals=[[0.0, 5.0]]) + path = self._write(nwbfile) + spiketrains = self._read(path)[0].segments[0].spiketrains + self.assertEqual(spiketrains[0].t_start, 0.0 * pq.s) + self.assertEqual(spiketrains[0].t_stop, 5.0 * pq.s) + assert_array_equal(spiketrains[0].annotations["nwb:obs_intervals"], [[0.0, 2.0], [4.0, 5.0]]) + assert_array_equal(spiketrains[1].magnitude, [0.5, 6.0]) + self.assertEqual(spiketrains[1].t_stop, 6.0 * pq.s) + + # the observation intervals should be preserved when writing + out_path = self.tmpdir / "rewritten.nwb" + block = self._read(path)[0] + NWBIO(out_path, mode="w").write_all_blocks([block]) + with pynwb.NWBHDF5IO(out_path, mode="r") as io: + assert_array_equal(io.read().units.get_unit_obs_intervals(0), [[0.0, 2.0], [4.0, 5.0]]) + + def test_units_without_spike_times(self): + # some ophys files use the units table only for cell metadata + nwbfile = _make_nwbfile() + nwbfile.add_unit_column("cell_id", "cell identifier") + nwbfile.add_unit(cell_id=42) + path = self._write(nwbfile) + block = self._read(path)[0] + self.assertEqual(len(block.segments), 0) + self.assertEqual(block.annotations["identifier"], "test-file") + + def test_epochs_without_tags(self): + nwbfile = _make_nwbfile() + nwbfile.add_epoch_column("condition", "experimental condition") + nwbfile.add_epoch(start_time=0.0, stop_time=1.0, condition="a") + nwbfile.add_epoch(start_time=2.0, stop_time=4.0, condition="b") + nwbfile.add_acquisition(pynwb.TimeSeries(name="ts", data=np.arange(10.0), unit="volts", rate=10.0)) + path = self._write(nwbfile) + for lazy in (False, True): + epoch = self._read(path, lazy=lazy)[0].segments[0].epochs[0] + self.assertEqual(epoch.shape, (2,)) + if lazy: + epoch = epoch.load() + assert_array_equal(epoch.times.magnitude, [0.0, 2.0]) + assert_array_equal(epoch.durations.magnitude, [1.0, 2.0]) + + # writing back + out_path = self.tmpdir / "rewritten.nwb" + NWBIO(out_path, mode="w").write_all_blocks(self._read(path)) + epoch = self._read(out_path)[0].segments[0].epochs[0] + assert_array_equal(epoch.times.magnitude, [0.0, 2.0]) + + def test_description_and_comments(self): + nwbfile = _make_nwbfile() + notes = {"protocol": "step", "amplitude": 5} + nwbfile.add_acquisition( + pynwb.TimeSeries( + name="ts", + data=np.arange(10.0), + unit="volts", + rate=10.0, + description="membrane potential", + comments="recorded on rig 2", + ) + ) + nwbfile.add_acquisition( + pynwb.TimeSeries( + name="ts_notes", data=np.arange(10.0), unit="volts", rate=10.0, description=json.dumps(notes) + ) + ) + nwbfile.add_acquisition( + pynwb.TimeSeries(name="ts_number", data=np.arange(10.0), unit="volts", rate=10.0, description="42") + ) + nwbfile.add_acquisition( + pynwb.misc.AnnotationSeries( + name="events", data=["a", "b"], timestamps=[0.1, 0.2], description=json.dumps(notes) + ) + ) + path = self._write(nwbfile) + + def check(segment): + signals = {sig.name: sig for sig in segment.analogsignals} + self.assertEqual(signals["ts"].description, "membrane potential") + self.assertEqual(signals["ts"].annotations["nwb:comments"], "recorded on rig 2") + self.assertIsNone(signals["ts_notes"].description) + self.assertEqual(signals["ts_notes"].annotations["notes"], notes) + self.assertNotIn("nwb:comments", signals["ts_notes"].annotations) + # a description that happens to be valid JSON, but not a JSON object, is kept as a string + self.assertEqual(signals["ts_number"].description, "42") + event = segment.events[0] + self.assertIsNone(event.description) + self.assertEqual(event.annotations["notes"], notes) + self.assertNotIn("protocol", event.annotations) + + block = self._read(path)[0] + check(block.segments[0]) + + out_path = self.tmpdir / "rewritten.nwb" + NWBIO(out_path, mode="w").write_all_blocks([block]) + check(self._read(out_path)[0].segments[0]) + + def test_fields_not_accepted_by_class(self): + # IZeroClampSeries has fixed values for bias_current etc., which can't be passed to the constructor + nwbfile = _make_nwbfile() + device = nwbfile.create_device(name="amplifier") + electrode = nwbfile.create_icephys_electrode(name="elec0", description="patch pipette", device=device) + nwbfile.add_acquisition( + pynwb.icephys.IZeroClampSeries( + name="izero", data=np.arange(10.0), electrode=electrode, rate=10.0, gain=1.0, sweep_number=0 + ) + ) + path = self._write(nwbfile) + block = self._read(path)[0] + self.assertIn("nwb:bias_current", block.segments[0].analogsignals[0].annotations) + + out_path = self.tmpdir / "rewritten.nwb" + with self.assertLogs("neo.io.nwbio.NWBIO", level="WARNING"): + NWBIO(out_path, mode="w").write_all_blocks([block]) + with pynwb.NWBHDF5IO(out_path, mode="r") as io: + self.assertIsInstance(io.read().acquisition["izero"], pynwb.icephys.IZeroClampSeries) + + def test_write_no_blocks(self): + with self.assertRaises(ValueError): + NWBIO(self.tmpdir / "empty.nwb", mode="w").write_all_blocks([]) + + if __name__ == "__main__": if HAVE_PYNWB: print("pynwb.__version__ = ", pynwb.__version__) From af1d790a543dc130d951534eeb22231cd726688c Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Fri, 13 Nov 2026 13:42:24 +0100 Subject: [PATCH 2/3] Write and read SpikeTrain descriptions in NWB files NWBIO does not store the description of a SpikeTrain. It is now written to a _description column of the units table, with an empty string representing None, and read back from that column when present. This partly addresses #1516. --- neo/io/nwbio.py | 15 ++++++++++----- neo/test/iotest/test_nwbio.py | 21 +++++++++++++++++++++ 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/neo/io/nwbio.py b/neo/io/nwbio.py index 67362736f..74b27cc6b 100644 --- a/neo/io/nwbio.py +++ b/neo/io/nwbio.py @@ -648,8 +648,7 @@ def write_all_blocks(self, blocks, validate=True, **kwargs): if sum(statistics(block)["SpikeTrain"]["count"] for block in blocks) > 0: self._nwbfile.add_unit_column("_name", "the name attribute of the SpikeTrain") - # nwbfile.add_unit_column('_description', - # 'the description attribute of the SpikeTrain') + self._nwbfile.add_unit_column("_description", "the description attribute of the SpikeTrain") self._nwbfile.add_unit_column("segment", "the name of the Neo Segment to which the SpikeTrain belongs") self._nwbfile.add_unit_column("block", "the name of the Neo Block to which the SpikeTrain belongs") @@ -916,7 +915,8 @@ def _write_spiketrain(self, nwbfile, spiketrain, name): spike_times=spiketrain.rescale("s").magnitude, obs_intervals=obs_intervals, _name=name, - # _description=spiketrain.description, + # an empty string represents a SpikeTrain without a description + _description=spiketrain.description or "", segment=segment.name, block=segment.block.name, ) @@ -1250,8 +1250,13 @@ def __init__(self, units_table, id, default_t_start=0.0, default_t_stop=0.0): name = units_table["_name"][id] or None else: name = None + if "_description" in units_table.colnames: + # NWB files created by recent versions of Neo also store the description + description = units_table["_description"][id] or None + else: + description = None self.shape = None # no way to get this without reading the data - _init_proxy(self, name=name, annotations=annotations) + _init_proxy(self, name=name, description=description, annotations=annotations) def load(self, time_slice=None, strict_slicing=True): """ @@ -1285,7 +1290,7 @@ def load(self, time_slice=None, strict_slicing=True): # left_sweep=None, name=self.name, # file_origin=None, - # description=None, + description=self.description, # array_annotations=None, **self.annotations, ) diff --git a/neo/test/iotest/test_nwbio.py b/neo/test/iotest/test_nwbio.py index c950b1bc0..05f475f28 100644 --- a/neo/test/iotest/test_nwbio.py +++ b/neo/test/iotest/test_nwbio.py @@ -350,6 +350,27 @@ def test_write_proxy_objects(self): assert_array_equal(original_segment.events[0].load().magnitude, retrieved_segment.events[0].magnitude) assert_array_equal(original_segment.epochs[0].load().magnitude, retrieved_segment.epochs[0].magnitude) + def test_roundtrip_descriptions(self): + # see https://github.com/NeuralEnsemble/python-neo/issues/1516 + block = Block(name="block", session_start_time=datetime.now().astimezone()) + segment = Segment(name="segment") + block.segments.append(segment) + segment.analogsignals.append( + AnalogSignal( + np.random.randn(10, 2) * pq.mV, sampling_rate=1 * pq.kHz, name="signal", description="a signal" + ) + ) + segment.spiketrains.append(SpikeTrain([1, 2, 3] * pq.s, t_stop=10 * pq.s, description="a spike train")) + segment.spiketrains.append(SpikeTrain([4, 5] * pq.s, t_stop=10 * pq.s)) + test_file_name = self.tmpdir / "test_roundtrip_descriptions.nwb" + NWBIO(filename=test_file_name, mode="w").write_all_blocks([block]) + + retrieved_segment = NWBIO(filename=test_file_name, mode="r").read_all_blocks()[0].segments[0] + self.assertEqual(retrieved_segment.analogsignals[0].description, "a signal") + self.assertEqual(retrieved_segment.spiketrains[0].description, "a spike train") + self.assertIsNone(retrieved_segment.spiketrains[1].description) + self.assertIsNone(retrieved_segment.spiketrains[1].name) + def _make_nwbfile(): """Create a minimal NWBFile, as it might be created by software other than Neo.""" From 0cf37d15eea1650ed63bc635177a9eba43606f7f Mon Sep 17 00:00:00 2001 From: Andrew Davison Date: Fri, 13 Nov 2026 16:09:23 +0100 Subject: [PATCH 3/3] Fix lazy loading of irregularly sampled signals from NWB files NWBIO uses its AnalogSignalProxy for both regularly and irregularly sampled signals, but the class declares itself a proxy for AnalogSignal, so with lazy=True, adding the proxy for an irregularly sampled signal to a Segment raises a TypeError. The proxy now declares IrregularlySampledSignal when the TimeSeries has no sampling rate. --- neo/io/nwbio.py | 2 ++ neo/test/iotest/test_nwbio.py | 7 +++++++ 2 files changed, 9 insertions(+) diff --git a/neo/io/nwbio.py b/neo/io/nwbio.py index 74b27cc6b..5bbbe0ef3 100644 --- a/neo/io/nwbio.py +++ b/neo/io/nwbio.py @@ -1005,6 +1005,8 @@ def __init__(self, timeseries, nwb_group): self.sampling_rate = timeseries.rate * pq.Hz else: self.sampling_rate = None + # this proxy class is used for both regularly and irregularly sampled signals + self.proxy_for = IrregularlySampledSignal annotations = {"nwb_group": nwb_group} description = timeseries.description diff --git a/neo/test/iotest/test_nwbio.py b/neo/test/iotest/test_nwbio.py index 05f475f28..fd9571eff 100644 --- a/neo/test/iotest/test_nwbio.py +++ b/neo/test/iotest/test_nwbio.py @@ -421,11 +421,18 @@ def test_read_lazy(self): nwbfile.add_acquisition( pynwb.TimeSeries(name="ts", data=np.arange(10.0), unit="volts", rate=1000.0, starting_time=0.0) ) + nwbfile.add_acquisition( + pynwb.TimeSeries(name="irregular", data=np.arange(5.0), unit="volts", timestamps=[0.1, 0.2, 0.4, 0.8, 1.6]) + ) path = self._write(nwbfile) segment = self._read(path, lazy=True)[0].segments[0] self.assertEqual(len(segment.analogsignals), 1) signal = segment.analogsignals[0].load() assert_array_equal(signal.magnitude.ravel(), np.arange(10.0)) + self.assertEqual(len(segment.irregularlysampledsignals), 1) + irregular = segment.irregularlysampledsignals[0].load() + self.assertIsInstance(irregular, IrregularlySampledSignal) + assert_array_equal(irregular.times.magnitude, [0.1, 0.2, 0.4, 0.8, 1.6]) def test_fields_that_cannot_be_annotations(self): # The `electrodes` field of an ElectricalSeries is a link to rows of the electrodes table,