diff --git a/synapse/cli/offline_hdf5_plotter.py b/synapse/cli/offline_hdf5_plotter.py index a4709f86..e5494be9 100644 --- a/synapse/cli/offline_hdf5_plotter.py +++ b/synapse/cli/offline_hdf5_plotter.py @@ -105,14 +105,20 @@ def load_h5_data(data_file, console, time_range=None): # List immediate groups (top-level only) print_tree(f, console) - # Get channel information + # Get channel information. `id` is one id per entry per frame, in frame + # order — physical electrode ids for recordings that could resolve them + # (see id_source), logical channel ids or row indices otherwise. channels = f["/general/extracellular_ephys/electrodes/"] channel_ids = channels["id"][:].tolist() number_of_channels = len(channel_ids) sample_rate = float(attributes["sample_rate_hz"]) console.print(f"Sample rate: {sample_rate} Hz") - console.print(f"Found {number_of_channels} channels") + id_source = channels.attrs.get("id_source") + if id_source is not None: + console.print(f"Found {number_of_channels} channels (id source: {id_source})") + else: + console.print(f"Found {number_of_channels} channels") # Get frame data info frame_data = f["/acquisition/ElectricalSeries"] @@ -204,8 +210,11 @@ def load_h5_data(data_file, console, time_range=None): actual_samples_per_channel, number_of_channels ) - # Create DataFrame - df = pd.DataFrame(reshaped_data, columns=range(number_of_channels)) + # Create DataFrame, labeling each column with the channel id of the + # row it came from. filter_channels selects by label (--channels + # takes channel ids), so positional column labels quietly selected + # the wrong traces for any recording whose ids aren't 0..N-1. + df = pd.DataFrame(reshaped_data, columns=channel_ids) return PlotData(data=df, sample_rate=sample_rate, channel_ids=channel_ids) diff --git a/synapse/cli/streaming.py b/synapse/cli/streaming.py index dc758d57..0283822f 100644 --- a/synapse/cli/streaming.py +++ b/synapse/cli/streaming.py @@ -16,6 +16,12 @@ from synapse.api.node_pb2 import NodeType from synapse.client.taps import Tap from synapse.utils.proto import load_device_config +from synapse.utils.electrode_ids import ( + GPIO_CHANNEL_ID_OFFSET, + ElectrodeRowIds, + build_channel_to_electrode_map, + derive_electrode_row_ids, +) from synapse.api.datatype_pb2 import BroadbandFrame @@ -252,7 +258,7 @@ def get_stats(self): def set_attributes( self, sample_rate_hz: float, - channels: list, + electrode_rows: ElectrodeRowIds, broadband_lsb_uv: float, session_description: str = "", ): @@ -269,7 +275,22 @@ def set_attributes( electrodes_group = self.file.create_group( "general/extracellular_ephys/electrodes" ) - electrodes_group.create_dataset("id", data=channels, dtype="uint32") + # `id` holds physical electrode ids where we could resolve them (see + # derive_electrode_row_ids), and degrades to logical channel ids or row + # indices otherwise — `id_source` says which. The logical channel id and + # channel type are written alongside it, as equal-length datasets, so + # downstream analysis can reconstruct the full mapping. + electrodes_group.create_dataset( + "id", data=electrode_rows.ids, dtype="uint32" + ) + electrodes_group.create_dataset( + "channel_id", data=electrode_rows.channel_ids, dtype="uint32" + ) + electrodes_group.create_dataset( + "channel_type", data=electrode_rows.types, dtype="uint32" + ) + electrodes_group.attrs["id_source"] = electrode_rows.source + electrodes_group.attrs["gpio_channel_id_offset"] = GPIO_CHANNEL_ID_OFFSET def start(self): """Start the writer thread""" @@ -685,8 +706,15 @@ def list_available_taps(args, device, console): ) -def detect_stream_parameters(broadband_tap, console): - """Detect sample rate and available channels from the first message""" +def detect_stream_parameters(broadband_tap, console, channel_to_electrode=None): + """Detect sample rate and channel identity from the first message + + Returns (sample_rate, electrode_rows, first_frame). `electrode_rows` carries + one id per frame_data entry: the physical electrode where the device + configuration lets us resolve it, and a documented fallback otherwise. It + replaces the old range(num_channels), which labeled every row by its + position and so mislabeled any sparse or non-ascending configuration. + """ console.log("[cyan]Detecting stream parameters from first message...[/cyan]") try: @@ -704,15 +732,23 @@ def detect_stream_parameters(broadband_tap, console): # Extract parameters sample_rate = first_frame.sample_rate_hz - num_channels = len(first_frame.frame_data) - available_channels = list(range(num_channels)) + electrode_rows = derive_electrode_row_ids(first_frame, channel_to_electrode) + if electrode_rows is None: + console.print( + "[bold red]First message carried no channel data[/bold red]" + ) + return None, None, None + num_channels = len(electrode_rows.ids) console.log(f"[green]Detected sample rate: {sample_rate} Hz[/green]") console.log( - f"[green]Detected {num_channels} channels (0-{num_channels - 1})[/green]" + f"[green]Detected {num_channels} channels " + f"(id source: {electrode_rows.source})[/green]" ) + for warning in electrode_rows.warnings: + console.print(f"[yellow]{warning}[/yellow]") - return sample_rate, available_channels, first_frame + return sample_rate, electrode_rows, first_frame except Exception as e: console.print(f"[bold red]Error detecting stream parameters: {e}[/bold red]") @@ -944,19 +980,24 @@ def read(args): console.print("[bold red]Failed to get broadband tap[/bold red]") return + # Get the latest info on the device. The stream identifies its rows only by + # logical channel id, so the applied configuration is the only place the + # physical electrode behind each row is recorded. + device_info = device.info() + channel_to_electrode = build_channel_to_electrode_map(device_info) + # Detect stream parameters from the first message - sample_rate, available_channels, first_frame = detect_stream_parameters( - broadband_tap, console + sample_rate, electrode_rows, first_frame = detect_stream_parameters( + broadband_tap, console, channel_to_electrode ) if sample_rate is None: console.print("[bold red]Failed to detect stream parameters[/bold red]") return + available_channels = electrode_rows.ids # Setup our HDF5 writer if output is requested writer = None if args.output: - # Get the latest info on the device - device_info = device.info() broadband_node = get_broadband_node_status(device_info) if broadband_node is None: console.print( @@ -967,7 +1008,7 @@ def read(args): writer = BroadbandFrameWriter(args.output) writer.set_attributes( sample_rate_hz=sample_rate, - channels=available_channels, + electrode_rows=electrode_rows, broadband_lsb_uv=broadband_lsb_uv, ) writer.start() diff --git a/synapse/tests/cli/test_broadband_frame_writer.py b/synapse/tests/cli/test_broadband_frame_writer.py new file mode 100644 index 00000000..bd56c673 --- /dev/null +++ b/synapse/tests/cli/test_broadband_frame_writer.py @@ -0,0 +1,76 @@ +import h5py + +from synapse.api.channel_pb2 import ChannelRange, ChannelType +from synapse.api.datatype_pb2 import BroadbandFrame +from synapse.cli.streaming import BroadbandFrameWriter +from synapse.utils.electrode_ids import ( + GPIO_CHANNEL_ID_OFFSET, + derive_electrode_row_ids, +) + + +def write_attributes(tmp_path, frame, channel_to_electrode=None): + """Run set_attributes for a frame and hand back the resulting file.""" + rows = derive_electrode_row_ids(frame, channel_to_electrode) + writer = BroadbandFrameWriter(str(tmp_path)) + try: + writer.set_attributes( + sample_rate_hz=30000.0, + electrode_rows=rows, + broadband_lsb_uv=0.195, + ) + finally: + writer.file.close() + return writer.filename + + +def test_writes_physical_electrode_ids(tmp_path): + frame = BroadbandFrame( + frame_data=[0, 0, 0, 0, 0], + channel_ranges=[ + ChannelRange( + type=ChannelType.ELECTRODE, count=5, channel_ids=[0, 1, 2, 3, 4] + ) + ], + ) + filename = write_attributes( + tmp_path, frame, {0: 6, 1: 16, 2: 10, 3: 14, 4: 20} + ) + + with h5py.File(filename, "r") as f: + electrodes = f["general/extracellular_ephys/electrodes"] + assert electrodes["id"][:].tolist() == [6, 16, 10, 14, 20] + assert electrodes["channel_id"][:].tolist() == [0, 1, 2, 3, 4] + assert electrodes["channel_type"][:].tolist() == [0, 0, 0, 0, 0] + assert electrodes.attrs["id_source"] == "electrode_map" + assert electrodes.attrs["gpio_channel_id_offset"] == GPIO_CHANNEL_ID_OFFSET + assert f.attrs["sample_rate_hz"] == 30000.0 + assert f.attrs["lsb_uv"] == 0.195 + + +def test_electrode_table_length_matches_the_per_frame_entry_count(tmp_path): + """GPIO rows are part of ElectricalSeries, so they need table rows too.""" + frame = BroadbandFrame( + frame_data=[0, 0, 0], + channel_ranges=[ + ChannelRange(type=ChannelType.ELECTRODE, count=2, channel_ids=[0, 1]), + ChannelRange(type=ChannelType.GPIO, count=1, channel_ids=[1]), + ], + ) + filename = write_attributes(tmp_path, frame, {0: 6, 1: 16}) + + with h5py.File(filename, "r") as f: + electrodes = f["general/extracellular_ephys/electrodes"] + assert len(electrodes["id"]) == len(frame.frame_data) + assert electrodes["id"][:].tolist() == [6, 16, GPIO_CHANNEL_ID_OFFSET + 1] + assert electrodes["channel_type"][:].tolist() == [0, 0, 1] + + +def test_records_the_legacy_fallback_in_id_source(tmp_path): + frame = BroadbandFrame(frame_data=[0, 0, 0]) + filename = write_attributes(tmp_path, frame) + + with h5py.File(filename, "r") as f: + electrodes = f["general/extracellular_ephys/electrodes"] + assert electrodes["id"][:].tolist() == [0, 1, 2] + assert electrodes.attrs["id_source"] == "positional" diff --git a/synapse/tests/utils/test_electrode_ids.py b/synapse/tests/utils/test_electrode_ids.py new file mode 100644 index 00000000..6fe17d7e --- /dev/null +++ b/synapse/tests/utils/test_electrode_ids.py @@ -0,0 +1,168 @@ +from synapse.api.channel_pb2 import Channel, ChannelRange, ChannelType +from synapse.api.datatype_pb2 import BroadbandFrame +from synapse.api.device_pb2 import DeviceInfo, Peripheral +from synapse.api.node_pb2 import NodeType +from synapse.utils.electrode_ids import ( + GPIO_CHANNEL_ID_OFFSET, + build_channel_to_electrode_map, + derive_electrode_row_ids, +) + + +def make_frame(num_rows, ranges=None): + frame = BroadbandFrame(frame_data=list(range(num_rows))) + for channel_range in ranges or []: + frame.channel_ranges.append(channel_range) + return frame + + +def electrode_range(count, channel_ids=None): + return ChannelRange( + type=ChannelType.ELECTRODE, count=count, channel_ids=channel_ids or [] + ) + + +def make_device_info(channels, peripheral=None): + info = DeviceInfo() + if peripheral is not None: + info.peripherals.append(peripheral) + node = info.configuration.nodes.add() + node.type = NodeType.kBroadbandSource + node.broadband_source.peripheral_id = ( + peripheral.peripheral_id if peripheral is not None else 0 + ) + node.broadband_source.signal.electrode.channels.extend(channels) + return info + + +# Horacio's report: 5 channels referenced to ground, configured out of ascending +# order, two of them on non-functional electrodes. +HORACIO_ELECTRODES = [6, 16, 10, 14, 20] +HORACIO_CHANNELS = [ + Channel(id=i, electrode_id=electrode_id, reference_id=520) + for i, electrode_id in enumerate(HORACIO_ELECTRODES) +] + + +def test_returns_none_for_an_empty_frame(): + assert derive_electrode_row_ids(make_frame(0)) is None + + +def test_maps_logical_channel_ids_to_physical_electrode_ids(): + frame = make_frame(5, [electrode_range(5, [0, 1, 2, 3, 4])]) + rows = derive_electrode_row_ids(frame, {0: 6, 1: 16, 2: 10, 3: 14, 4: 20}) + + assert rows.ids == HORACIO_ELECTRODES + assert rows.channel_ids == [0, 1, 2, 3, 4] + assert rows.source == "electrode_map" + assert rows.warnings == [] + + +def test_preserves_the_configured_electrode_order(): + """The configured order is the recorded order — no ascending re-sort.""" + frame = make_frame(5, [electrode_range(5, [0, 1, 2, 3, 4])]) + rows = derive_electrode_row_ids( + frame, build_channel_to_electrode_map(make_device_info(HORACIO_CHANNELS)) + ) + + assert rows.ids == [6, 16, 10, 14, 20] + assert rows.ids != sorted(rows.ids) + + +def test_falls_back_to_logical_ids_without_a_map(): + frame = make_frame(3, [electrode_range(3, [7, 8, 9])]) + rows = derive_electrode_row_ids(frame) + + assert rows.ids == [7, 8, 9] + assert rows.source == "channel_ranges" + assert any("No channel-to-electrode map" in w for w in rows.warnings) + + +def test_falls_back_to_logical_ids_when_the_map_is_incomplete(): + frame = make_frame(3, [electrode_range(3, [0, 1, 2])]) + rows = derive_electrode_row_ids(frame, {0: 6, 1: 16}) + + assert rows.ids == [0, 1, 2] + assert rows.source == "channel_ranges" + assert any("does not cover every streamed channel" in w for w in rows.warnings) + + +def test_rekeys_gpio_rows_and_preserves_the_row_count(): + frame = make_frame( + 4, + [ + electrode_range(2, [0, 1]), + ChannelRange(type=ChannelType.GPIO, count=2, channel_ids=[1, 3]), + ], + ) + rows = derive_electrode_row_ids(frame, {0: 6, 1: 16}) + + assert rows.ids == [6, 16, GPIO_CHANNEL_ID_OFFSET + 1, GPIO_CHANNEL_ID_OFFSET + 3] + assert rows.types == [ + ChannelType.ELECTRODE, + ChannelType.ELECTRODE, + ChannelType.GPIO, + ChannelType.GPIO, + ] + assert len(rows.ids) == len(frame.frame_data) + assert rows.source == "electrode_map" + + +def test_uses_contiguous_ids_when_a_range_omits_channel_ids(): + """Pre-fix Nixel512 firmware sends a bare count.""" + frame = make_frame(3, [electrode_range(3)]) + rows = derive_electrode_row_ids(frame) + + assert rows.ids == [0, 1, 2] + assert rows.source == "positional" + assert any("omit channel_ids" in w for w in rows.warnings) + + +def test_resolves_electrodes_for_a_contiguous_range_when_mapped(): + frame = make_frame(3, [electrode_range(3)]) + rows = derive_electrode_row_ids(frame, {0: 6, 1: 16, 2: 10}) + + assert rows.ids == [6, 16, 10] + assert rows.source == "electrode_map" + + +def test_falls_back_to_positional_ids_without_channel_ranges(): + rows = derive_electrode_row_ids(make_frame(3), {0: 6, 1: 16, 2: 10}) + + assert rows.channel_ids == [0, 1, 2] + assert rows.ids == [6, 16, 10] + assert any("no channel_ranges" in w for w in rows.warnings) + + +def test_falls_back_to_positional_ids_when_ranges_disagree_with_frame_size(): + frame = make_frame(3, [electrode_range(5, [0, 1, 2, 3, 4])]) + rows = derive_electrode_row_ids(frame) + + assert rows.channel_ids == [0, 1, 2] + assert rows.source == "positional" + assert any("describe 5 channels" in w for w in rows.warnings) + + +def test_builds_the_map_from_the_device_configuration(): + assert build_channel_to_electrode_map(make_device_info(HORACIO_CHANNELS)) == { + 0: 6, + 1: 16, + 2: 10, + 3: 14, + 4: 20, + } + + +def test_skips_peripherals_with_placeholder_electrode_ids(): + virtual = Peripheral( + name="SciFi Virtual Recording Peripheral", + vendor="Science Corporation", + peripheral_id=3, + ) + channels = [Channel(id=i, electrode_id=2 * i) for i in range(3)] + + assert build_channel_to_electrode_map(make_device_info(channels, virtual)) == {} + + +def test_builds_an_empty_map_for_a_missing_device_info(): + assert build_channel_to_electrode_map(None) == {} diff --git a/synapse/utils/electrode_ids.py b/synapse/utils/electrode_ids.py new file mode 100644 index 00000000..0de69f94 --- /dev/null +++ b/synapse/utils/electrode_ids.py @@ -0,0 +1,211 @@ +"""Per-row channel identity for broadband recordings. + +The broadband stream carries one entry per channel per frame and identifies +those entries only by *logical* channel id — the 0-based index of the channel +within the BroadbandSource configuration. The physical electrode behind each +logical channel lives in the device configuration, not in the stream, so +recording the stream alone loses the mapping from row to electrode. + +This module reconstructs that identity: `build_channel_to_electrode_map` reads +the logical -> physical mapping out of a DeviceInfo, and +`derive_electrode_row_ids` combines it with a frame's `channel_ranges` to +produce one id per frame_data entry. + +Mirrors `deriveElectrodeRowIds` in nexus-desktop +(src/ipc/main/hdf5-writer.ts); the two write the same HDF5 layout and must stay +in sync. +""" + +from dataclasses import dataclass +from typing import Dict, List, Optional + +from synapse.api.channel_pb2 import ChannelType +from synapse.api.node_pb2 import NodeType + +# The device emits GPIO channels in the broadband stream keyed by their GPIO +# line index, which collides with the contiguous electrode channel ids — +# electrode channel 1 and GPIO line 1 are both "1". GPIO rows are re-keyed into +# this high namespace so every id in the electrode table is unique. Mirrors +# GPIO_CHANNEL_ID_OFFSET in nexus-desktop. +GPIO_CHANNEL_ID_OFFSET = 1_000_000 + +# Peripherals whose configuration carries placeholder electrode ids rather than +# real ones. The virtual recording peripheral generates synthetic data and +# derives electrode_id = 2 * id, so there is no physical electrode behind those +# ids and remapping the stream through them only corrupts the labels (0, 1, 2 -> +# 0, 2, 4). Matched on the (name, vendor) the firmware reports. +_PLACEHOLDER_ELECTRODE_PERIPHERALS = frozenset( + {("SciFi Virtual Recording Peripheral", "Science Corporation")} +) + +# Where the per-row ids in general/extracellular_ephys/electrodes/id came from. +ID_SOURCE_ELECTRODE_MAP = "electrode_map" # real physical electrode ids +ID_SOURCE_CHANNEL_RANGES = "channel_ranges" # device-reported logical channel ids +ID_SOURCE_POSITIONAL = "positional" # row index; no identity at all (legacy) + + +@dataclass +class ElectrodeRowIds: + """Per-frame_data-entry channel identity. + + Every list has exactly len(frame.frame_data) entries, so the length of + electrodes/id keeps matching the number of entries per frame in + ElectricalSeries (which may include non-electrode GPIO channels). + """ + + ids: List[int] + """One id per frame_data entry — what goes into electrodes/id.""" + + channel_ids: List[int] + """Logical stream channel id per frame_data entry.""" + + types: List[int] + """synapse.ChannelType per frame_data entry.""" + + source: str + """One of the ID_SOURCE_* constants above.""" + + warnings: List[str] + """Set when identity had to degrade; the caller reports these.""" + + +def build_channel_to_electrode_map(device_info) -> Dict[int, int]: + """Build the logical channel id -> physical electrode id lookup. + + Reads every BroadbandSource node's configured channels, each of which + carries both an `id` (the logical/stream id) and an `electrode_id` (the + real electrode). Nodes fed by a placeholder-electrode peripheral are + skipped so their channels stay unmapped and fall back to logical ids. + """ + mapping: Dict[int, int] = {} + if device_info is None: + return mapping + + peripherals_by_id = { + peripheral.peripheral_id: peripheral + for peripheral in device_info.peripherals + } + + for node in device_info.configuration.nodes: + if node.type != NodeType.kBroadbandSource: + continue + peripheral = peripherals_by_id.get(node.broadband_source.peripheral_id) + if peripheral is not None and ( + peripheral.name, + peripheral.vendor, + ) in _PLACEHOLDER_ELECTRODE_PERIPHERALS: + continue + for channel in node.broadband_source.signal.electrode.channels: + mapping[channel.id] = channel.electrode_id + + return mapping + + +def derive_electrode_row_ids( + frame, channel_to_electrode: Optional[Dict[int, int]] = None +) -> Optional[ElectrodeRowIds]: + """Derive per-row channel identity for a BroadbandFrame. + + The frame describes its own layout in `channel_ranges`: each range gives a + ChannelType, a count, and — when the peripheral populates it — the logical + channel ids in frame order. Ranges that omit channel_ids are contiguous, + with ids running in frame order. + + Returns None for an empty frame. Otherwise the result always has one entry + per frame_data entry, degrading through the ID_SOURCE_* levels (and + recording a warning) whenever real identity isn't available. + """ + num_rows = len(frame.frame_data) if frame is not None else 0 + if num_rows <= 0: + return None + + channel_ids = [0] * num_rows + types = [ChannelType.ELECTRODE] * num_rows + warnings: List[str] = [] + has_explicit_ids = False + used_ranges = False + + ranges = list(frame.channel_ranges) + if ranges: + declared = sum(r.count for r in ranges) + if declared != num_rows: + warnings.append( + f"BroadbandFrame.channel_ranges describe {declared} channels but the " + f"frame has {num_rows} entries; falling back to positional channel ids." + ) + else: + row = 0 + next_contiguous_id = 0 + for channel_range in ranges: + explicit = list(channel_range.channel_ids) + for i in range(channel_range.count): + if i < len(explicit): + channel_ids[row] = explicit[i] + has_explicit_ids = True + else: + # Contiguous range: ids run in frame order from where + # the previous range left off. + channel_ids[row] = next_contiguous_id + next_contiguous_id += 1 + types[row] = channel_range.type + row += 1 + used_ranges = True + else: + warnings.append( + "BroadbandFrame carries no channel_ranges; falling back to positional " + "channel ids (all entries assumed to be electrodes)." + ) + + if not used_ranges: + channel_ids = list(range(num_rows)) + types = [ChannelType.ELECTRODE] * num_rows + elif not has_explicit_ids: + warnings.append( + "BroadbandFrame.channel_ranges omit channel_ids; channel identity falls " + "back to positional ordering within each range." + ) + + # Resolve logical channel ids to physical electrode ids, but only when + # *every* electrode row resolves: a partial map would mix the two id spaces + # in one dataset, which is worse than an honest logical-id dataset. + electrode_rows = [i for i in range(num_rows) if types[i] != ChannelType.GPIO] + all_mapped = bool(channel_to_electrode) and all( + channel_ids[i] in channel_to_electrode for i in electrode_rows + ) + use_map = all_mapped and len(electrode_rows) > 0 + + ids: List[int] = [] + for i in range(num_rows): + if types[i] == ChannelType.GPIO: + ids.append(GPIO_CHANNEL_ID_OFFSET + channel_ids[i]) + elif use_map: + ids.append(channel_to_electrode[channel_ids[i]]) + else: + ids.append(channel_ids[i]) + + if not use_map and electrode_rows: + if channel_to_electrode: + warnings.append( + "Channel-to-electrode map does not cover every streamed channel; " + "writing logical channel ids instead of physical electrode ids." + ) + else: + warnings.append( + "No channel-to-electrode map supplied; writing logical channel ids " + "instead of physical electrode ids." + ) + + if use_map: + source = ID_SOURCE_ELECTRODE_MAP + elif has_explicit_ids: + source = ID_SOURCE_CHANNEL_RANGES + else: + source = ID_SOURCE_POSITIONAL + + return ElectrodeRowIds( + ids=ids, + channel_ids=channel_ids, + types=[int(t) for t in types], + source=source, + warnings=warnings, + )