diff --git a/src/spikeinterface/extractors/cbin_ibl.py b/src/spikeinterface/extractors/cbin_ibl.py index 2a53b999e3..1756049dea 100644 --- a/src/spikeinterface/extractors/cbin_ibl.py +++ b/src/spikeinterface/extractors/cbin_ibl.py @@ -1,4 +1,5 @@ from pathlib import Path +import re import warnings import numpy as np @@ -9,6 +10,89 @@ from spikeinterface.core.core_tools import define_function_from_class +def _parse_spikeglx_meta_table(meta_value): + """Return parenthesized entries from a SpikeGLX meta value. + + Accepts either the raw string form ``"(a)(b)(c)"`` or an already-parsed list. + """ + if isinstance(meta_value, str): + return re.findall(r"\(([^()]*)\)", meta_value) + return list(meta_value) + + +def _parse_saved_channel_subset(subset_text): + if subset_text is None or subset_text == "all": + return None + channels = [] + for entry in subset_text.split(","): + if ":" in entry: + start, stop = entry.split(":") + channels.extend(np.arange(int(start), int(stop) + 1)) + else: + channels.append(int(entry)) + return np.asarray(channels, dtype="int64") + + +def _get_saved_channel_indices(meta): + subset_text = meta.get("snsSaveChanSubset_orig") + if subset_text is None: + subset_text = meta.get("snsSaveChanSubset") + return _parse_saved_channel_subset(subset_text) + + +def _read_cbin_probe(meta_file, meta): + """Build a Probe for a SpikeGLX meta, honoring IBL shank-split subsets. + This is a specific to IBL cbin files, that sometimes do not have `~` in front of the .meta file fields. + """ + if "snsSaveChanSubset_orig" not in meta: + return probeinterface.read_spikeglx(meta_file) + + # probeinterface removed the single-shot `_read_imro_string` helper used + # previously (in 0.3.2). The equivalent flow is: build the full geometric + # probe, parse the IMRO table to a dict, resolve which contact IDs are + # active, and slice the full probe down to those contacts. + from probeinterface.neuropixels_tools import ( + build_neuropixels_probe, + _parse_imro_string, + _get_imro_active_contact_ids, + _annotate_probe_with_adc_sampling_info, + ) + + imDatPrb_pn = meta.get("imDatPrb_pn") + # neo's read_meta_file leaves imroTbl as a raw string when the meta field + # has no leading '~' (the IBL cbin case noted above). _parse_imro_string + # requires the string form, so rebuild it defensively if neo split it. + imro_table_string = meta["imroTbl"] + if isinstance(imro_table_string, list): + imro_table_string = "(" + ")(".join(imro_table_string) + ")" + + full_probe = build_neuropixels_probe(probe_part_number=imDatPrb_pn) + imro_per_channel = _parse_imro_string(imro_table_string) + active_contact_ids = _get_imro_active_contact_ids(imro_per_channel) + contact_id_to_index = {cid: i for i, cid in enumerate(full_probe.contact_ids)} + selected = np.array([contact_id_to_index[cid] for cid in active_contact_ids]) + probe = full_probe.get_slice(selected) + + # Attach ADC sampling annotations (num_channels_per_adc and per-contact + # adc_sample_order) needed by get_neuropixels_sample_shifts_from_probe for + # phase_shift preprocessing. Must run while contacts are still in readout + # channel order, i.e. before the snsSaveChanSubset_orig slice below. + adc_sampling_table = probe.annotations.get("adc_sampling_table") + _annotate_probe_with_adc_sampling_info(probe, adc_sampling_table) + + saved = _get_saved_channel_indices(meta) + if saved is not None: + saved = saved[saved < probe.get_contact_count()] + if saved.size != probe.get_contact_count(): + probe = probe.get_slice(saved) + probe.annotate(serial_number=meta.get("imDatPrb_sn", meta.get("imProbeSN", None))) + probe.annotate(part_number=meta.get("imDatPrb_pn", None)) + probe.annotate(port=meta.get("imDatPrb_port", None)) + probe.annotate(slot=meta.get("imDatPrb_slot", None)) + probe.set_device_channel_indices(np.arange(probe.get_contact_count())) + return probe + + class CompressedBinaryIblExtractor(BaseRecording): """Load IBL data as an extractor object. @@ -99,7 +183,7 @@ def __init__(self, folder_path=None, load_sync_channel=False, stream_name="ap", self.set_channel_offsets(offsets) if not load_sync_channel: - probe = probeinterface.read_spikeglx(meta_file) + probe = _read_cbin_probe(meta_file, meta) if probe.shank_ids is not None: self.set_probe(probe, in_place=True, group_mode="by_shank") @@ -204,9 +288,31 @@ def extract_stream_info(meta_file, meta): info["stream_kind"] = stream_kind info["stream_name"] = stream_name info["units"] = units - info["channel_names"] = [txt.split(";")[0] for txt in meta["snsChanMap"]] + + chan_map_entries = _parse_spikeglx_meta_table(meta["snsChanMap"]) + # First entry is the header tuple like "(384,0,1)"; drop it when present. + if chan_map_entries and "," in chan_map_entries[0] and ";" not in chan_map_entries[0]: + chan_map_entries = chan_map_entries[1:] + full_channel_names = [entry.split(";")[0] for entry in chan_map_entries] + + saved_indices = _get_saved_channel_indices(meta) + if saved_indices is not None and len(full_channel_names) >= saved_indices.max() + 1: + channel_names = [full_channel_names[i] for i in saved_indices] + else: + channel_names = full_channel_names + + channel_offsets = np.zeros(num_chan) + + if len(channel_names) != num_chan or channel_gains.shape[0] != num_chan: + raise ValueError( + "Parsed channel metadata does not match nSavedChans: " + f"channel_names={len(channel_names)}, channel_gains={channel_gains.shape[0]}, " + f"num_chan={num_chan}" + ) + + info["channel_names"] = channel_names info["channel_gains"] = channel_gains - info["channel_offsets"] = np.zeros(info["num_chan"]) + info["channel_offsets"] = channel_offsets info["has_sync_trace"] = has_sync_trace return info diff --git a/src/spikeinterface/extractors/tests/test_cbin_ibl_extractors.py b/src/spikeinterface/extractors/tests/test_cbin_ibl_extractors.py index 905994197e..bc508979b9 100644 --- a/src/spikeinterface/extractors/tests/test_cbin_ibl_extractors.py +++ b/src/spikeinterface/extractors/tests/test_cbin_ibl_extractors.py @@ -1,8 +1,17 @@ import pytest import unittest +import numpy as np + from spikeinterface.extractors.extractor_classes import CompressedBinaryIblExtractor, read_cbin_ibl +from spikeinterface.extractors.cbin_ibl import ( + _get_saved_channel_indices, + _parse_saved_channel_subset, + _parse_spikeglx_meta_table, + _read_cbin_probe, + extract_stream_info, +) from spikeinterface.extractors.tests.common_tests import RecordingCommonTestSuite, SortingCommonTestSuite @@ -13,6 +22,201 @@ class CompressedBinaryIblExtractorTest(RecordingCommonTestSuite, unittest.TestCa entities = [] +# -------------------------------------------------------------------------------- +# NP2.4 shank-split meta handling +# +# When IBL splits an NP2.4 recording into four single-shank recordings, the +# generated .meta file has two peculiarities: +# +# * none of the list-valued fields carry the leading "~" that SpikeGLX writes, +# so neo.rawio.spikeglxrawio.read_meta_file leaves them as raw strings; +# * the true acquisition channels live in "snsSaveChanSubset_orig" while +# "snsSaveChanSubset" only holds the local, post-split renumbering. +# +# probeinterface.read_spikeglx uses "snsSaveChanSubset" and therefore returns a +# silently wrong probe for these files (sync counted as a contact, contacts +# spread over the wrong shanks). Correcting that is what _read_cbin_probe does. +# These tests need only a .meta file, so they run without any downloaded data. +# -------------------------------------------------------------------------------- + +# 16 acquired channels in blocks of four, alternating shank 0 / shank 1, so the +# saved subset is non-contiguous exactly like a real 4-shank split. Channels +# 4-7 and 12-15 are the shank-1 channels, i.e. electrodes 0-7 on shank 1. +N_ACQUIRED_CHANNELS = 16 +SYNC_CHANNEL_INDEX = N_ACQUIRED_CHANNELS +SAVED_CHANNEL_SUBSET_ORIG = "4:7,12:15,16" +EXPECTED_CHANNEL_NAMES = ["AP4", "AP5", "AP6", "AP7", "AP12", "AP13", "AP14", "AP15", "SY0"] + + +def _build_shank_split_meta_text(include_subset_orig=True): + """Build a minimal NP2.4 shank-split .meta text, deliberately without any "~". + + Parameters + ---------- + include_subset_orig : bool, default: True + If False, omit "snsSaveChanSubset_orig" so the meta looks like a normal + SpikeGLX file and _read_cbin_probe should defer to probeinterface. + + Returns + ------- + str + The full contents of the .meta file. + """ + imro_entries, geom_entries, chan_map_entries = [], [], [] + for channel in range(N_ACQUIRED_CHANNELS): + block, position = divmod(channel, 4) + shank = block % 2 + electrode = (block // 2) * 4 + position + # type-24 (NP2.4) IMRO entry: (channel shank bank refid electrode) + imro_entries.append(f"({channel} {shank} 0 1 {electrode})") + # NP2.4 geometry: two columns at 27/59 um from the shank edge, 15 um row pitch + x_um = 27 if electrode % 2 == 0 else 59 + geom_entries.append(f"({shank}:{x_um}:{15 * (electrode // 2)}:1)") + chan_map_entries.append(f"(AP{channel};{channel}:{channel})") + chan_map_entries.append(f"(SY0;{SYNC_CHANNEL_INDEX}:{SYNC_CHANNEL_INDEX})") + + num_saved_channels = len(EXPECTED_CHANNEL_NAMES) + fields = { + "typeThis": "imec", + "acqApLfSy": f"{N_ACQUIRED_CHANNELS},0,1", + "snsApLfSy": f"{num_saved_channels - 1},0,1", + "nSavedChans": str(num_saved_channels), + # the local, post-split numbering that read_spikeglx wrongly relies on + "snsSaveChanSubset": f"0:{num_saved_channels - 1}", + "imSampRate": "30000", + "fileSizeBytes": str(num_saved_channels * 2 * 3000), + "imAiRangeMax": "0.5", + "imMaxInt": "8192", + "imChan0apGain": "80", + "imDatPrb_type": "24", + "imDatPrb_pn": "NP2010", + "imDatPrb_sn": "20472319942", + "imDatPrb_port": "1", + "imDatPrb_slot": "3", + "imroTbl": f"(24,{N_ACQUIRED_CHANNELS})" + "".join(imro_entries), + "snsChanMap": f"({N_ACQUIRED_CHANNELS},0,1)" + "".join(chan_map_entries), + "snsGeomMap": "(NP2010,4,250,70)" + "".join(geom_entries), + } + if include_subset_orig: + fields["snsSaveChanSubset_orig"] = SAVED_CHANNEL_SUBSET_ORIG + + # NOTE: written without the leading "~" that SpikeGLX puts on list-valued + # fields. That omission is the whole point of this fixture. + return "".join(f"{key}={value}\n" for key, value in fields.items()) + + +def _read_meta(meta_file): + """Parse a .meta file the same way CompressedBinaryIblExtractor does.""" + from neo.rawio.spikeglxrawio import read_meta_file + + return read_meta_file(str(meta_file)) + + +@pytest.fixture +def shank_split_meta_file(tmp_path): + """Path to a minimal NP2.4 shank-split meta file.""" + meta_file = tmp_path / "np24_shank_split.ap.meta" + meta_file.write_text(_build_shank_split_meta_text()) + return meta_file + + +def test_parse_spikeglx_meta_table_accepts_string_and_list(): + """The parser must cope with both the raw "~"-less string and a parsed list.""" + assert _parse_spikeglx_meta_table("(a)(b)(c)") == ["a", "b", "c"] + assert _parse_spikeglx_meta_table(["a", "b", "c"]) == ["a", "b", "c"] + + +def test_parse_saved_channel_subset(): + """Ranges, singletons and the "all" sentinel.""" + assert _parse_saved_channel_subset("all") is None + assert _parse_saved_channel_subset(None) is None + np.testing.assert_array_equal(_parse_saved_channel_subset("0:3"), [0, 1, 2, 3]) + np.testing.assert_array_equal(_parse_saved_channel_subset("4:7,12:15,16"), [4, 5, 6, 7, 12, 13, 14, 15, 16]) + + +def test_saved_channel_indices_prefer_orig(shank_split_meta_file): + """snsSaveChanSubset_orig must take precedence over the local snsSaveChanSubset.""" + meta = _read_meta(shank_split_meta_file) + np.testing.assert_array_equal(_get_saved_channel_indices(meta), [4, 5, 6, 7, 12, 13, 14, 15, 16]) + + +def test_meta_without_tilde_stays_a_string(shank_split_meta_file): + """Guards the premise of the fix: neo cannot list-parse these fields.""" + meta = _read_meta(shank_split_meta_file) + assert isinstance(meta["imroTbl"], str) + assert isinstance(meta["snsChanMap"], str) + + +def test_read_cbin_probe_selects_the_correct_shank(shank_split_meta_file): + """The probe must hold only the 8 shank-1 contacts, at the right positions.""" + meta = _read_meta(shank_split_meta_file) + probe = _read_cbin_probe(str(shank_split_meta_file), meta) + + assert probe.get_contact_count() == 8 + assert set(probe.shank_ids) == {"1"} + # shank pitch 250 um, two columns 32 um apart, four rows at 15 um pitch + np.testing.assert_array_equal(np.unique(probe.contact_positions[:, 0]), [250.0, 282.0]) + np.testing.assert_array_equal(np.unique(probe.contact_positions[:, 1]), [0.0, 15.0, 30.0, 45.0]) + np.testing.assert_array_equal(probe.device_channel_indices, np.arange(8)) + assert probe.annotations["serial_number"] == "20472319942" + assert probe.annotations["part_number"] == "NP2010" + + +def test_read_cbin_probe_differs_from_plain_read_spikeglx(shank_split_meta_file): + """Regression guard: plain read_spikeglx returns a wrong probe for this meta. + + It applies the local "snsSaveChanSubset" to the geometry table, so the sync + channel becomes a contact and the contacts land on two different shanks. + If probeinterface ever learns about "snsSaveChanSubset_orig" this test will + fail, which is the signal that _read_cbin_probe can be removed. + """ + import probeinterface + + meta = _read_meta(shank_split_meta_file) + wrong_probe = probeinterface.read_spikeglx(str(shank_split_meta_file)) + correct_probe = _read_cbin_probe(str(shank_split_meta_file), meta) + + assert wrong_probe.get_contact_count() == 9 # 8 contacts + sync wrongly counted as a contact + assert set(wrong_probe.shank_ids) == {"0", "1"} # contacts wrongly spread over two shanks + assert correct_probe.get_contact_count() != wrong_probe.get_contact_count() + + +def test_read_cbin_probe_falls_back_when_no_subset_orig(tmp_path): + """A normal SpikeGLX meta must be passed through to probeinterface untouched.""" + import probeinterface + + meta_file = tmp_path / "plain.ap.meta" + meta_file.write_text(_build_shank_split_meta_text(include_subset_orig=False)) + meta = _read_meta(meta_file) + + probe = _read_cbin_probe(str(meta_file), meta) + reference_probe = probeinterface.read_spikeglx(str(meta_file)) + np.testing.assert_array_equal(probe.contact_positions, reference_probe.contact_positions) + + +def test_extract_stream_info_uses_the_original_subset(shank_split_meta_file): + """Channel names must follow snsSaveChanSubset_orig, not the local numbering.""" + meta = _read_meta(shank_split_meta_file) + info = extract_stream_info(str(shank_split_meta_file), meta) + + assert info["num_chan"] == len(EXPECTED_CHANNEL_NAMES) + assert info["channel_names"] == EXPECTED_CHANNEL_NAMES + assert info["channel_gains"].shape[0] == info["num_chan"] + assert info["channel_offsets"].shape[0] == info["num_chan"] + assert info["has_sync_trace"] + + +def test_extract_stream_info_raises_on_inconsistent_meta(tmp_path): + """The nSavedChans consistency guard must fire on a corrupted meta.""" + meta_text = _build_shank_split_meta_text().replace("nSavedChans=9", "nSavedChans=42") + meta_file = tmp_path / "bad.ap.meta" + meta_file.write_text(meta_text) + meta = _read_meta(meta_file) + + with pytest.raises(ValueError, match="does not match nSavedChans"): + extract_stream_info(str(meta_file), meta) + + # ~ def test_read_cbin_ibl(): # ~ base_folder = '/media/samuel/dataspikesorting/DataSpikeSorting/olivier_destripe/' # ~ data_folder = base_folder + '4c04120d-523a-4795-ba8f-49dbb8d9f63a'