From 23cef31ac371d56149e82d14620da1198c58992a Mon Sep 17 00:00:00 2001
From: Alessio Buccino <alejoe9187@gmail.com>
Date: Tue, 16 Dec 2025 13:19:20 +0100
Subject: [PATCH 01/46] wip: load_analyzer_from_nwb function

---
 .../extractors/extractor_classes.py           |   2 +
 .../extractors/nwbextractors.py               | 287 +++++++++++++++++-
 2 files changed, 288 insertions(+), 1 deletion(-)

diff --git a/src/spikeinterface/extractors/extractor_classes.py b/src/spikeinterface/extractors/extractor_classes.py
index 4f0f586f18..7a91927ec9 100644
--- a/src/spikeinterface/extractors/extractor_classes.py
+++ b/src/spikeinterface/extractors/extractor_classes.py
@@ -35,6 +35,7 @@
     read_nwb_recording,
     read_nwb_sorting,
     read_nwb_timeseries,
+    load_analyzer_from_nwb,
 )
 
 from .cbin_ibl import CompressedBinaryIblExtractor, read_cbin_ibl
@@ -194,6 +195,7 @@
 __all__.extend(
     [
         "read_nwb",  # convenience function for multiple nwb formats
+        "load_analyzer_from_nwb",
         "recording_extractor_full_dict",
         "sorting_extractor_full_dict",
         "event_extractor_full_dict",
diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 976e752a62..069f014970 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -7,7 +7,14 @@
 import numpy as np
 
 from spikeinterface import get_global_tmp_folder
-from spikeinterface.core import BaseRecording, BaseRecordingSegment, BaseSorting, BaseSortingSegment
+from spikeinterface.core import (
+    BaseRecording,
+    BaseRecordingSegment,
+    BaseSorting,
+    BaseSortingSegment,
+    SortingAnalyzer,
+    get_default_analyzer_extension_params,
+)
 from spikeinterface.core.core_tools import define_function_from_class
 
 
@@ -1259,6 +1266,7 @@ def _fetch_sorting_segment_info_backend(
 
         # need this for later
         self.units_table = units_table
+        self._file = open_file
 
         return unit_ids, spike_times_data, spike_times_index_data
 
@@ -1789,3 +1797,280 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
         outputs = outputs[0]
 
     return outputs
+
+
+def load_analyzer_from_nwb(
+    file_path: str | Path,
+    t_start: float | None = None,
+    sampling_frequency: float | None = None,
+    electrical_series_path: str | None = None,
+    unit_table_path: str | None = None,
+    stream_mode: Literal["fsspec", "remfile", "zarr"] | None = None,
+    stream_cache_path: str | Path | None = None,
+    cache: bool = False,
+    storage_options: dict | None = None,
+    use_pynwb: bool = False,
+    group_name: str | None = None,
+    compute_extra: List[str] | None = ["unit_locations", "correlograms"],
+    compute_extra_params: dict | None = None,
+    verbose: bool = False,
+) -> SortingAnalyzer:
+    import pandas as pd
+    from spikeinterface.metrics.template import ComputeTemplateMetrics
+    from spikeinterface.metrics.quality import ComputeQualityMetrics
+
+    # try to read recording object to get the analyzer
+    try:
+        recording = NwbRecordingExtractor(
+            file_path=file_path,
+            electrical_series_path=electrical_series_path,
+            stream_mode=stream_mode,
+            stream_cache_path=stream_cache_path,
+            cache=cache,
+            storage_options=storage_options,
+            use_pynwb=use_pynwb,
+        )
+    except Exception:
+        if verbose:
+            print("Could not load recording, proceeding without it")
+        recording = None
+
+    t_start_tmp = 0 if t_start is None else t_start
+
+    sorting_tmp = NwbSortingExtractor(
+        file_path=file_path,
+        electrical_series_path=electrical_series_path,
+        unit_table_path=unit_table_path,
+        stream_mode=stream_mode,
+        stream_cache_path=stream_cache_path,
+        cache=cache,
+        storage_options=storage_options,
+        use_pynwb=use_pynwb,
+        t_start=t_start_tmp,
+        sampling_frequency=sampling_frequency,
+    )
+
+    if recording is None and t_start is None:
+        # re-estimate t_start from spike times
+        if verbose:
+            print("Re-estimating t_start from spike_times")
+        t_start_new = np.min(sorting_tmp._sorting_segments[0].spike_times_data) - 0.001
+        if verbose:
+            print(f"Found new t_start: {t_start_new} s")
+        sorting = NwbSortingExtractor(
+            file_path=file_path,
+            electrical_series_path=electrical_series_path,
+            unit_table_path=unit_table_path,
+            stream_mode=stream_mode,
+            stream_cache_path=stream_cache_path,
+            cache=cache,
+            storage_options=storage_options,
+            use_pynwb=use_pynwb,
+            t_start=t_start_new,
+            sampling_frequency=sampling_frequency,
+        )
+    else:
+        sorting = sorting_tmp
+
+    if use_pynwb:
+        units = sorting.units_table
+        colnames = units.colnames
+        units = units.to_dataframe(index=True)
+    else:
+        units_dset = sorting._file["units"]
+        units = make_df(units_dset)
+        colnames = units.columns
+
+    electrodes_indices = None
+    if use_pynwb:
+        electrodes_table = sorting._nwbfile.electrodes.to_dataframe(index=True)
+        if "electrodes" in colnames:
+            electrodes_indices = units["electrodes"]
+    else:
+        electrodes_table = make_df(sorting._file["/general/extracellular_ephys/electrodes"])
+        if "electrodes" in colnames:
+            electrodes_indices = electrodes_indices = units["electrodes"][:]
+
+    if electrodes_indices is not None:
+        # here we assume all groups are the same for each unit, so we just check one.
+        if "group_name" in electrodes_table.columns:
+            group_names = np.array([electrodes_table.iloc[int(ei[0])]["group_name"] for ei in electrodes_indices])
+            if len(np.unique(group_names)) > 0:
+                if group_name is None:
+                    raise Exception(
+                        f"More than one group, use group_name option to select units. Available groups: {np.unique(group_names)}"
+                    )
+                else:
+                    unit_mask = group_names == group_name
+                    if verbose:
+                        print(f"Selecting {sum(unit_mask)} / {len(units)} units from {group_name}")
+                    sorting = sorting.select_units(unit_ids=sorting.unit_ids[unit_mask])
+                    units = units.loc[units.index[unit_mask]]
+                    electrodes_indices = units["electrodes"]
+
+    # TODO: figure out sparsity
+
+    # handle recording if available
+    if recording is not None:
+        # check groups
+        group_names = np.unique(recording.get_channel_groups())
+        if group_name is not None and len(group_names) > 1:
+            recording = recording.split_by("group")[group_name]
+        rec_attributes = None
+    else:
+        recording = None
+        rec_attributes = {}
+
+        # get sliced electrodes table from electrode_indices union
+        electrode_indices_all = []
+        for ei in electrodes_indices:
+            electrode_indices_all.extend(ei)
+        electrode_indices_all = np.sort(np.unique(electrode_indices_all))
+        if verbose:
+            print(f"Found {len(electrode_indices_all)} electrodes")
+        electrodes_table_sliced = electrodes_table.iloc[electrode_indices_all]
+        if "channel_name" in electrodes_table_sliced:
+            channel_ids = electrodes_table_sliced["channel_name"][:]
+        else:
+            channel_ids = electrodes_table_sliced["id"][:]
+        num_samples = [sorting.to_spike_vector()[-1]["sample_index"]]
+        rec_attributes = dict(
+            channel_ids=channel_ids,
+            sampling_frequency=sorting.sampling_frequency,
+            num_channels=len(channel_ids),
+            num_samples=num_samples,
+            is_filtered=True,
+            dtype="float32",
+        )
+        # make a probegroup
+        electrode_colnames = electrodes_table_sliced.columns
+        assert (
+            "rel_x" in electrode_colnames and "rel_y" in electrode_colnames
+        ), "'rel_x' and 'rel_y' should be columns in the electrode name"
+        locations = np.array([electrodes_table_sliced["rel_x"][:], electrodes_table_sliced["rel_y"][:]]).T
+        probegroup = create_dummy_probegroup_from_locations(locations)
+        rec_attributes["probegroup"] = probegroup
+
+    # instantiate analyzer
+    analyzer = SortingAnalyzer.create_memory(
+        sorting=sorting, recording=recording, sparsity=None, rec_attributes=rec_attributes, return_in_uV=True
+    )
+
+    # templates
+    if "waveform_mean" in units:
+        from spikeinterface.core.analyzer_extension_core import ComputeTemplates, ComputeRandomSpikes
+
+        # instantiate templates
+        analyzer.compute("random_spikes", method="all")
+
+        templates_ext = ComputeTemplates(sorting_analyzer=analyzer)
+        templates_avg_data = np.array([t for t in units["waveform_mean"].values]).astype("float")
+        total_ms = templates_avg_data.shape[1] / analyzer.sampling_frequency * 1000
+        template_params = get_default_analyzer_extension_params("templates")
+        if total_ms != template_params["ms_before"] + template_params["ms_after"]:
+            if verbose:
+                print("Guessing correct template cutouts")
+            template_params["ms_before"] = int(1 / 3 * total_ms)
+            template_params["ms_after"] = total_ms - template_params["ms_before"]
+        template_params["operators"] = ["average", "std"]
+        templates_ext.set_params(**template_params)
+        templates_avg_data = np.array([t for t in units["waveform_mean"].values]).astype("float")
+        templates_ext.data["average"] = templates_avg_data
+        if "waveforms_sd" in units:
+            templates_std_data = np.array([t for t in units["waveform_sd"].values]).astype("float")
+        else:
+            templates_std_data = np.zeros_like(templates_avg_data)
+        templates_ext.data["std"] = templates_std_data
+        templates_ext.run_info["run_completed"] = True
+
+        analyzer.extensions["templates"] = templates_ext
+
+    template_metric_columns = ComputeTemplateMetrics.get_metric_columns()
+    quality_metric_columns = ComputeQualityMetrics.get_metric_columns()
+
+    tm = pd.DataFrame(index=sorting.unit_ids)
+    qm = pd.DataFrame(index=sorting.unit_ids)
+
+    for col in units.columns:
+        if col in template_metric_columns:
+            tm.loc[:, col] = units[col].values
+        if col in quality_metric_columns:
+            qm.loc[:, col] = units[col].values
+
+    if len(tm.columns) > 0:
+        if verbose:
+            print("Adding template metrics")
+        tm_ext = ComputeTemplateMetrics(analyzer)
+        tm_ext.data["metrics"] = tm
+        tm_ext.run_info["run_completed"] = True
+        analyzer.extensions["template_metrics"] = tm_ext
+    if len(qm.columns) > 0:
+        if verbose:
+            print("Adding quality metrics")
+        qm_ext = ComputeQualityMetrics(analyzer)
+        qm_ext.data["metrics"] = qm
+        qm_ext.run_info["run_completed"] = True
+        analyzer.extensions["quality_metrics"] = qm_ext
+
+    # compute extra required
+    if compute_extra is not None:
+        if verbose:
+            print(f"Computing extra extensions: {compute_extra}")
+        compute_extra_params = {} if compute_extra_params is None else compute_extra_params
+        analyzer.compute(compute_extra, **compute_extra_params)
+
+    return analyzer
+
+
+def create_dummy_probegroup_from_locations(locations, shape="circle", shape_params={"radius": 1}):
+    """
+    Creates a "dummy" probe based on locations.
+
+    Parameters
+    ----------
+    locations : np.array
+        Array with channel locations (num_channels, ndim) [ndim can be 2 or 3]
+    shape : str, default: "circle"
+        Electrode shapes
+    shape_params : dict, default: {"radius": 1}
+        Shape parameters
+
+    Returns
+    -------
+    probe : Probe
+        The created probe
+    """
+    from probeinterface import Probe, ProbeGroup
+
+    ndim = locations.shape[1]
+    assert ndim == 2
+    probe = Probe(ndim=2)
+    probe.set_contacts(locations, shapes=shape, shape_params=shape_params)
+    probe.set_device_channel_indices(np.arange(len(probe.contact_positions)))
+    probe.create_auto_shape()
+    probegroup = ProbeGroup()
+    probegroup.add_probe(probe)
+
+    return probegroup
+
+
+def make_df(group):
+    """Makes pandas DataFrame from hdf5/zarr NWB group"""
+    import pandas as pd
+
+    colnames = list(group.keys())
+    data = {}
+    for col in colnames:
+        if "_index" in col:
+            continue
+        item = group[col][:]
+        if f"{col}_index" in colnames:
+            item = np.split(item, group[f"{col}_index"][:])[:-1]
+            data[col] = item
+        elif item.ndim > 1:
+            data[col] = [item_flat for item_flat in item]
+        else:
+            data[col] = item
+    df = pd.DataFrame(data=data)
+    df.set_index("id", inplace=True)
+    return df

From fcba7d3394dba85af5ffc8a5d7076125611ca071 Mon Sep 17 00:00:00 2001
From: Alessio Buccino <alejoe9187@gmail.com>
Date: Tue, 23 Dec 2025 15:58:15 +0100
Subject: [PATCH 02/46] Make dtype/is_filtered optional attrs and suggestions
 from code review

---
 src/spikeinterface/core/sortinganalyzer.py    |  4 +-
 src/spikeinterface/extractors/__init__.py     |  1 +
 .../extractors/extractor_classes.py           |  2 -
 .../extractors/nwbextractors.py               | 63 ++++++++++---------
 4 files changed, 37 insertions(+), 33 deletions(-)

diff --git a/src/spikeinterface/core/sortinganalyzer.py b/src/spikeinterface/core/sortinganalyzer.py
index 1870c24e7a..d2a9cbe2df 100644
--- a/src/spikeinterface/core/sortinganalyzer.py
+++ b/src/spikeinterface/core/sortinganalyzer.py
@@ -1501,7 +1501,7 @@ def is_sparse(self) -> bool:
         return self.sparsity is not None
 
     def is_filtered(self) -> bool:
-        return self.rec_attributes["is_filtered"]
+        return self.rec_attributes.get("is_filtered", True)
 
     def get_sorting_provenance(self):
         """
@@ -1602,7 +1602,7 @@ def get_sorting_property(self, key) -> np.ndarray:
         return self.sorting.get_property(key)
 
     def get_dtype(self):
-        return self.rec_attributes["dtype"]
+        return self.rec_attributes.get("dtype")
 
     def get_num_units(self) -> int:
         return self.sorting.get_num_units()
diff --git a/src/spikeinterface/extractors/__init__.py b/src/spikeinterface/extractors/__init__.py
index 216c668e0e..000c3832a5 100644
--- a/src/spikeinterface/extractors/__init__.py
+++ b/src/spikeinterface/extractors/__init__.py
@@ -6,6 +6,7 @@
 from .neuropixels_utils import get_neuropixels_channel_groups, get_neuropixels_sample_shifts
 from .neoextractors import get_neo_num_blocks, get_neo_streams
 from .phykilosortextractors import read_kilosort_as_analyzer
+from .nwbextractors import read_nwb_as_analyzer
 
 from warnings import warn
 
diff --git a/src/spikeinterface/extractors/extractor_classes.py b/src/spikeinterface/extractors/extractor_classes.py
index 7a91927ec9..4f0f586f18 100644
--- a/src/spikeinterface/extractors/extractor_classes.py
+++ b/src/spikeinterface/extractors/extractor_classes.py
@@ -35,7 +35,6 @@
     read_nwb_recording,
     read_nwb_sorting,
     read_nwb_timeseries,
-    load_analyzer_from_nwb,
 )
 
 from .cbin_ibl import CompressedBinaryIblExtractor, read_cbin_ibl
@@ -195,7 +194,6 @@
 __all__.extend(
     [
         "read_nwb",  # convenience function for multiple nwb formats
-        "load_analyzer_from_nwb",
         "recording_extractor_full_dict",
         "sorting_extractor_full_dict",
         "event_extractor_full_dict",
diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 069f014970..2c683088d1 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1799,7 +1799,7 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
     return outputs
 
 
-def load_analyzer_from_nwb(
+def read_nwb_as_analyzer(
     file_path: str | Path,
     t_start: float | None = None,
     sampling_frequency: float | None = None,
@@ -1878,7 +1878,7 @@ def load_analyzer_from_nwb(
         units = units.to_dataframe(index=True)
     else:
         units_dset = sorting._file["units"]
-        units = make_df(units_dset)
+        units = _create_df_from_nwb_table(units_dset)
         colnames = units.columns
 
     electrodes_indices = None
@@ -1887,7 +1887,7 @@ def load_analyzer_from_nwb(
         if "electrodes" in colnames:
             electrodes_indices = units["electrodes"]
     else:
-        electrodes_table = make_df(sorting._file["/general/extracellular_ephys/electrodes"])
+        electrodes_table = _create_df_from_nwb_table(sorting._file["/general/extracellular_ephys/electrodes"])
         if "electrodes" in colnames:
             electrodes_indices = electrodes_indices = units["electrodes"][:]
 
@@ -1939,8 +1939,6 @@ def load_analyzer_from_nwb(
             sampling_frequency=sorting.sampling_frequency,
             num_channels=len(channel_ids),
             num_samples=num_samples,
-            is_filtered=True,
-            dtype="float32",
         )
         # make a probegroup
         electrode_colnames = electrodes_table_sliced.columns
@@ -1948,7 +1946,7 @@ def load_analyzer_from_nwb(
             "rel_x" in electrode_colnames and "rel_y" in electrode_colnames
         ), "'rel_x' and 'rel_y' should be columns in the electrode name"
         locations = np.array([electrodes_table_sliced["rel_x"][:], electrodes_table_sliced["rel_y"][:]]).T
-        probegroup = create_dummy_probegroup_from_locations(locations)
+        probegroup = _create_dummy_probegroup_from_locations(locations)
         rec_attributes["probegroup"] = probegroup
 
     # instantiate analyzer
@@ -1960,18 +1958,22 @@ def load_analyzer_from_nwb(
     if "waveform_mean" in units:
         from spikeinterface.core.analyzer_extension_core import ComputeTemplates, ComputeRandomSpikes
 
-        # instantiate templates
+        # compute random spikes, which is a dependency for templates
+        # since we don't know the spike samples, we compute with method 'all'
         analyzer.compute("random_spikes", method="all")
 
+        # instantiate templates
         templates_ext = ComputeTemplates(sorting_analyzer=analyzer)
         templates_avg_data = np.array([t for t in units["waveform_mean"].values]).astype("float")
         total_ms = templates_avg_data.shape[1] / analyzer.sampling_frequency * 1000
-        template_params = get_default_analyzer_extension_params("templates")
-        if total_ms != template_params["ms_before"] + template_params["ms_after"]:
-            if verbose:
-                print("Guessing correct template cutouts")
-            template_params["ms_before"] = int(1 / 3 * total_ms)
-            template_params["ms_after"] = total_ms - template_params["ms_before"]
+        # estimate ms_before and ms_after from minimum point in the average template
+        nbefore = np.unravel_index(np.argmin(templates_avg_data, axis=1), templates_avg_data.shape)[1]
+        print(nbefore)
+        ms_before = int(nbefore / analyzer.sampling_frequency * 1000)
+        ms_after = int(total_ms - ms_before)
+        template_params = {}
+        template_params["ms_before"] = ms_before
+        template_params["ms_after"] = ms_after
         template_params["operators"] = ["average", "std"]
         templates_ext.set_params(**template_params)
         templates_avg_data = np.array([t for t in units["waveform_mean"].values]).astype("float")
@@ -1988,29 +1990,32 @@ def load_analyzer_from_nwb(
     template_metric_columns = ComputeTemplateMetrics.get_metric_columns()
     quality_metric_columns = ComputeQualityMetrics.get_metric_columns()
 
-    tm = pd.DataFrame(index=sorting.unit_ids)
-    qm = pd.DataFrame(index=sorting.unit_ids)
+    template_metrics_df = pd.DataFrame(index=sorting.unit_ids)
+    quality_metric_df = pd.DataFrame(index=sorting.unit_ids)
 
     for col in units.columns:
         if col in template_metric_columns:
-            tm.loc[:, col] = units[col].values
+            template_metrics_df.loc[:, col] = units[col].values
         if col in quality_metric_columns:
-            qm.loc[:, col] = units[col].values
+            quality_metric_df.loc[:, col] = units[col].values
 
-    if len(tm.columns) > 0:
+    if len(template_metrics_df.columns) > 0:
         if verbose:
             print("Adding template metrics")
-        tm_ext = ComputeTemplateMetrics(analyzer)
-        tm_ext.data["metrics"] = tm
-        tm_ext.run_info["run_completed"] = True
-        analyzer.extensions["template_metrics"] = tm_ext
-    if len(qm.columns) > 0:
+        template_metrics_ext = ComputeTemplateMetrics(analyzer)
+        template_metrics_ext.data["metrics"] = template_metrics_df
+        template_metrics_ext.run_info["run_completed"] = True
+        # cast to correct dtypes
+        template_metrics_ext._cast_metrics()
+        analyzer.extensions["template_metrics"] = template_metrics_ext
+    if len(quality_metric_df.columns) > 0:
         if verbose:
             print("Adding quality metrics")
-        qm_ext = ComputeQualityMetrics(analyzer)
-        qm_ext.data["metrics"] = qm
-        qm_ext.run_info["run_completed"] = True
-        analyzer.extensions["quality_metrics"] = qm_ext
+        quality_metrics_ext = ComputeQualityMetrics(analyzer)
+        quality_metrics_ext.data["metrics"] = quality_metric_df
+        quality_metrics_ext.run_info["run_completed"] = True
+        quality_metrics_ext._cast_metrics()
+        analyzer.extensions["quality_metrics"] = quality_metrics_ext
 
     # compute extra required
     if compute_extra is not None:
@@ -2022,7 +2027,7 @@ def load_analyzer_from_nwb(
     return analyzer
 
 
-def create_dummy_probegroup_from_locations(locations, shape="circle", shape_params={"radius": 1}):
+def _create_dummy_probegroup_from_locations(locations, shape="circle", shape_params={"radius": 1}):
     """
     Creates a "dummy" probe based on locations.
 
@@ -2054,7 +2059,7 @@ def create_dummy_probegroup_from_locations(locations, shape="circle", shape_para
     return probegroup
 
 
-def make_df(group):
+def _create_df_from_nwb_table(group):
     """Makes pandas DataFrame from hdf5/zarr NWB group"""
     import pandas as pd
 

From b72a39d8f38d9ed458693d9b9dc45390533f5d9b Mon Sep 17 00:00:00 2001
From: Alessio Buccino <alejoe9187@gmail.com>
Date: Tue, 23 Dec 2025 16:56:48 +0100
Subject: [PATCH 03/46] fix ibl tests

---
 src/spikeinterface/extractors/tests/test_iblextractors.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/spikeinterface/extractors/tests/test_iblextractors.py b/src/spikeinterface/extractors/tests/test_iblextractors.py
index 972a8e7bb0..56d01e38cf 100644
--- a/src/spikeinterface/extractors/tests/test_iblextractors.py
+++ b/src/spikeinterface/extractors/tests/test_iblextractors.py
@@ -76,8 +76,8 @@ def test_offsets(self):
 
     def test_probe_representation(self):
         probe = self.recording.get_probe()
-        expected_probe_representation = "Probe - 384ch - 1shanks"
-        assert repr(probe) == expected_probe_representation
+        expected_probe_representation = "Probe - 384ch"
+        assert expected_probe_representation in repr(probe)
 
     def test_property_keys(self):
         expected_property_keys = [

From 8e9e12d871406640b07cb082db6f5843b623c33b Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Tue, 7 Jul 2026 16:29:14 +0100
Subject: [PATCH 04/46] fix file references

---
 src/spikeinterface/extractors/nwbextractors.py | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index dd41621e41..23b3472b08 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1780,17 +1780,17 @@ def read_nwb_as_analyzer(
         colnames = units.colnames
         units = units.to_dataframe(index=True)
     else:
-        units_dset = sorting._file["units"]
+        units_dset = sorting.units_table
         units = _create_df_from_nwb_table(units_dset)
         colnames = units.columns
 
     electrodes_indices = None
     if use_pynwb:
-        electrodes_table = sorting._nwbfile.electrodes.to_dataframe(index=True)
+        electrodes_table = sorting._reader.nwbfile.electrodes.to_dataframe(index=True)
         if "electrodes" in colnames:
             electrodes_indices = units["electrodes"]
     else:
-        electrodes_table = _create_df_from_nwb_table(sorting._file["/general/extracellular_ephys/electrodes"])
+        electrodes_table = _create_df_from_nwb_table(sorting._reader.file["/general/extracellular_ephys/electrodes"])
         if "electrodes" in colnames:
             electrodes_indices = electrodes_indices = units["electrodes"][:]
 
@@ -1798,7 +1798,7 @@ def read_nwb_as_analyzer(
         # here we assume all groups are the same for each unit, so we just check one.
         if "group_name" in electrodes_table.columns:
             group_names = np.array([electrodes_table.iloc[int(ei[0])]["group_name"] for ei in electrodes_indices])
-            if len(np.unique(group_names)) > 0:
+            if len(np.unique(group_names)) > 1:
                 if group_name is None:
                     raise Exception(
                         f"More than one group, use group_name option to select units. Available groups: {np.unique(group_names)}"
@@ -1854,7 +1854,13 @@ def read_nwb_as_analyzer(
 
     # instantiate analyzer
     analyzer = SortingAnalyzer.create_memory(
-        sorting=sorting, recording=recording, sparsity=None, rec_attributes=rec_attributes, return_in_uV=True
+        sorting=sorting,
+        recording=recording,
+        sparsity=None,
+        return_in_uV=True,
+        peak_sign="neg",
+        peak_mode="extremum",
+        rec_attributes=rec_attributes,
     )
 
     # templates

From 97a6d23729a960edc1ad2595b9bb0db6bfc61aee Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Tue, 7 Jul 2026 16:51:51 +0100
Subject: [PATCH 05/46] fix main indices

---
 .../extractors/nwbextractors.py               | 96 +++++++++++++------
 1 file changed, 66 insertions(+), 30 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 23b3472b08..23a1407acd 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1852,11 +1852,48 @@ def read_nwb_as_analyzer(
         probegroup = _create_dummy_probegroup_from_locations(locations)
         rec_attributes["probegroup"] = probegroup
 
+    # Per-unit sparsity and channel map from the Units `electrodes` region. Each unit's
+    # `waveform_mean` is stored only on a subset of channels (near its peak); the region gives
+    # which channels those are. We use it both for the analyzer sparsity and to scatter the
+    # waveforms onto their true channel positions instead of stacking the sparse block densely.
+    sparsity = None
+    unit_local_channels = None
+    if electrodes_indices is not None:
+        from spikeinterface.core.sparsity import ChannelSparsity
+
+        if recording is not None:
+            analyzer_channel_ids = list(recording.get_channel_ids())
+        else:
+            analyzer_channel_ids = list(channel_ids)
+        num_channels = len(analyzer_channel_ids)
+
+        if "channel_name" in electrodes_table.columns:
+            electrode_row_to_channel_id = electrodes_table["channel_name"].to_numpy()
+        else:
+            electrode_row_to_channel_id = electrodes_table.index.to_numpy()
+        position_of_channel_id = {channel_id: pos for pos, channel_id in enumerate(analyzer_channel_ids)}
+
+        unit_local_channels = []
+        sparsity_mask = np.zeros((sorting.get_num_units(), num_channels), dtype=bool)
+        for unit_index, region in enumerate(electrodes_indices):
+            positions = [
+                position_of_channel_id[electrode_row_to_channel_id[int(electrode_row)]]
+                for electrode_row in region
+                if electrode_row_to_channel_id[int(electrode_row)] in position_of_channel_id
+            ]
+            unit_local_channels.append(positions)
+            sparsity_mask[unit_index, positions] = True
+        sparsity = ChannelSparsity(
+            mask=sparsity_mask,
+            unit_ids=np.asarray(sorting.unit_ids),
+            channel_ids=np.asarray(analyzer_channel_ids),
+        )
+
     # instantiate analyzer
     analyzer = SortingAnalyzer.create_memory(
         sorting=sorting,
         recording=recording,
-        sparsity=None,
+        sparsity=sparsity,
         return_in_uV=True,
         peak_sign="neg",
         peak_mode="extremum",
@@ -1864,34 +1901,35 @@ def read_nwb_as_analyzer(
     )
 
     # templates
-    if "waveform_mean" in units:
-        from spikeinterface.core.analyzer_extension_core import ComputeTemplates, ComputeRandomSpikes
+    if "waveform_mean" in units and unit_local_channels is not None:
+        from spikeinterface.core.analyzer_extension_core import ComputeTemplates
 
-        # compute random spikes, which is a dependency for templates
-        # since we don't know the spike samples, we compute with method 'all'
+        # random spikes is a dependency for templates; the raw spike samples are not stored so we
+        # select all of them
         analyzer.compute("random_spikes", method="all")
 
-        # instantiate templates
+        waveform_mean = np.array([np.asarray(t, dtype="float") for t in units["waveform_mean"].values])
+        num_samples_template = waveform_mean.shape[1]
+
+        # scatter each unit's real waveform columns onto their channels; the NaN-padded columns of
+        # units with fewer channels are simply never touched
+        dense_templates = np.zeros((sorting.get_num_units(), num_samples_template, num_channels), dtype="float32")
+        for unit_index, positions in enumerate(unit_local_channels):
+            k = len(positions)
+            dense_templates[unit_index][:, positions] = waveform_mean[unit_index][:, :k]
+
+        # scalar alignment (nbefore): the sample where the average peak-channel amplitude is largest
+        peak_amplitude_per_sample = np.abs(np.nan_to_num(waveform_mean)).max(axis=2).mean(axis=0)
+        nbefore = int(np.argmax(peak_amplitude_per_sample))
+
         templates_ext = ComputeTemplates(sorting_analyzer=analyzer)
-        templates_avg_data = np.array([t for t in units["waveform_mean"].values]).astype("float")
-        total_ms = templates_avg_data.shape[1] / analyzer.sampling_frequency * 1000
-        # estimate ms_before and ms_after from minimum point in the average template
-        nbefore = np.unravel_index(np.argmin(templates_avg_data, axis=1), templates_avg_data.shape)[1]
-        print(nbefore)
-        ms_before = int(nbefore / analyzer.sampling_frequency * 1000)
-        ms_after = int(total_ms - ms_before)
-        template_params = {}
-        template_params["ms_before"] = ms_before
-        template_params["ms_after"] = ms_after
-        template_params["operators"] = ["average", "std"]
-        templates_ext.set_params(**template_params)
-        templates_avg_data = np.array([t for t in units["waveform_mean"].values]).astype("float")
-        templates_ext.data["average"] = templates_avg_data
-        if "waveforms_sd" in units:
-            templates_std_data = np.array([t for t in units["waveform_sd"].values]).astype("float")
-        else:
-            templates_std_data = np.zeros_like(templates_avg_data)
-        templates_ext.data["std"] = templates_std_data
+        templates_ext.set_params(
+            ms_before=nbefore / analyzer.sampling_frequency * 1000,
+            ms_after=(num_samples_template - nbefore) / analyzer.sampling_frequency * 1000,
+            operators=["average", "std"],
+        )
+        templates_ext.data["average"] = dense_templates
+        templates_ext.data["std"] = np.zeros_like(dense_templates)
         templates_ext.run_info["run_completed"] = True
 
         analyzer.extensions["templates"] = templates_ext
@@ -1912,18 +1950,16 @@ def read_nwb_as_analyzer(
         if verbose:
             print("Adding template metrics")
         template_metrics_ext = ComputeTemplateMetrics(analyzer)
-        template_metrics_ext.data["metrics"] = template_metrics_df
-        template_metrics_ext.run_info["run_completed"] = True
         # cast to correct dtypes
-        template_metrics_ext._cast_metrics()
+        template_metrics_ext.data["metrics"] = template_metrics_ext._cast_metrics(template_metrics_df)
+        template_metrics_ext.run_info["run_completed"] = True
         analyzer.extensions["template_metrics"] = template_metrics_ext
     if len(quality_metric_df.columns) > 0:
         if verbose:
             print("Adding quality metrics")
         quality_metrics_ext = ComputeQualityMetrics(analyzer)
-        quality_metrics_ext.data["metrics"] = quality_metric_df
+        quality_metrics_ext.data["metrics"] = quality_metrics_ext._cast_metrics(quality_metric_df)
         quality_metrics_ext.run_info["run_completed"] = True
-        quality_metrics_ext._cast_metrics()
         analyzer.extensions["quality_metrics"] = quality_metrics_ext
 
     # compute extra required

From 20d38450cef2fa4b252cd0bc14f6b4da9f944beb Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 8 Jul 2026 11:03:35 +0100
Subject: [PATCH 06/46] Fix performance

---
 .../extractors/nwbextractors.py               | 81 +++++++++++--------
 1 file changed, 48 insertions(+), 33 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index b80daf8862..3f5b343375 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -441,6 +441,38 @@ def _read_column(self, table, name, indices=None):
             values = np.array([v.decode("utf-8") if isinstance(v, bytes) else v for v in values])
         return values
 
+    def _read_table_property(self, table, name, columns, indices=None):
+        # Load one DynamicTable column as a per-row property, or return None if it cannot be one.
+        # Scalar columns are read directly (and optionally reordered to `indices`). Ragged columns
+        # (those with a `<name>_index`) are decided from their index alone.
+        #
+        # Performance, read the (tiny) index before the (potentially huge) data. Whether a ragged
+        # column can be a per-row property is fully determined by its index: unequal per-row lengths
+        # cannot stack into a per-row array, and a ragged column cannot be reordered to a channel
+        # subset. So we inspect the index first and, when the column is unusable, return None WITHOUT
+        # ever reading its data. This matters because ragged columns are frequently per-spike, e.g. a
+        # Units table's spike amplitudes or spike depths, which are as long as the whole spike vector
+        # (tens of millions of values, hundreds of MB), while the index is only `num_rows` values. The
+        # old code read the full data and then discarded it after finding the shapes unequal. On local
+        # disk that waste is nearly invisible; over a network file (remfile / streaming from DANDI) it
+        # is severe, because each of those columns becomes many range requests for bytes we throw away.
+        # For a typical IBL session this avoided reading ~330 MB (two ~20M-long per-spike columns) that
+        # were only ever skipped. Do not reorder this to read `data` before checking the index.
+        index_name = f"{name}_index"
+        if index_name not in columns:
+            return self._read_column(table, name, indices)
+
+        index = table[index_name]
+        index = index.data[:] if hasattr(index, "data") else index[:]
+        if indices is not None or np.unique(np.diff(index, prepend=0)).size != 1:
+            return None
+        data = table[name]
+        if hasattr(data, "data"):  # pynwb resolves the ragged structure when sliced
+            return data[:]
+        data = np.asarray(data[:])
+        starts = np.concatenate([[0], index[:-1]])
+        return [data[start:end] for start, end in zip(starts, index)]
+
     def close(self):
         # Release the open handle (raw hdf5 / zarr store, or the pynwb read IO).
         if self.file is not None:
@@ -607,7 +639,8 @@ def _has_electrical_series(self):
     def _fetch_unit_properties(self):
         # Return the unit-table columns as a name -> per-unit values mapping in a backend-independent
         # format. Skips id / spike-time columns, index columns, nested ragged arrays, and ragged
-        # columns whose per-unit shapes are unequal.
+        # columns whose per-unit shapes are unequal. The last are detected from the (small) index
+        # alone, so their (potentially large, e.g. per-spike) data is never read.
         columns = self.units_column_names
 
         properties_to_skip = ["spike_times", "spike_times_index", "unit_name", "id"]
@@ -618,36 +651,13 @@ def _fetch_unit_properties(self):
 
         properties = dict()
         for property_name in properties_to_add:
-            corresponding_index_name = f"{property_name}_index"
-            not_ragged_array = corresponding_index_name not in columns
-            if not_ragged_array:
-                # _read_column materializes a plain, decoded, DynamicTableRegion-safe per-unit array
-                values = self._read_column(self.units_table, property_name)
-            else:  # TODO if we want we could make this recursive to handle nested ragged arrays
-                data = self.units_table[property_name][:]
-                data_index = self.units_table[corresponding_index_name]
-                if hasattr(data_index, "data"):
-                    # for pynwb we need to get the data from the data attribute
-                    data_index = data_index.data[:]
-                else:
-                    data_index = data_index[:]
-                index_spacing = np.diff(data_index, prepend=0)
-                all_index_spacing_are_the_same = np.unique(index_spacing).size == 1
-                if all_index_spacing_are_the_same:
-                    if hasattr(self.units_table[corresponding_index_name], "data"):
-                        # ragged array indexing is handled by pynwb
-                        values = data
-                    else:
-                        # ravel array based on data_index
-                        start_indices = [0] + list(data_index[:-1])
-                        end_indices = list(data_index)
-                        values = [
-                            data[start_index:end_index] for start_index, end_index in zip(start_indices, end_indices)
-                        ]
-                else:
-                    warnings.warn(f"Skipping {property_name} because of unequal shapes across units")
-                    continue
-            properties[property_name] = values
+            values = self._read_table_property(self.units_table, property_name, columns)
+            # A None result is a ragged column with a variable length per unit (typically per-spike
+            # data such as spike amplitudes): not a per-unit property, so it is skipped and its data is
+            # never read. No warning is emitted, this is expected for normal files and there is nothing
+            # for the user to act on.
+            if values is not None:
+                properties[property_name] = values
 
         return properties
 
@@ -751,7 +761,9 @@ def column_names(self):
 
     def read_electrode_property(self, name):
         # An electrode property is one electrodes-table column, reordered to this series' channels.
-        return self._read_column(self.electrodes_table, name, self.electrodes_indices)
+        # Ragged electrode columns cannot be per-channel scalars, so they are skipped (None returned)
+        # without their data being read.
+        return self._read_table_property(self.electrodes_table, name, self.column_names, self.electrodes_indices)
 
     def _read_ids(self):
         if self.reading_method == "use_pynwb":
@@ -1048,8 +1060,11 @@ def __init__(
             for column in self._reader.column_names:
                 if column in columns_mapped_to_core_fields:
                     continue
+                values = self._reader.read_electrode_property(column)
+                if values is None:  # ragged electrode column, not a per-channel property
+                    continue
                 property_name = "brain_area" if column == "location" else column
-                self.set_property(property_name, self._reader.read_electrode_property(column))
+                self.set_property(property_name, values)
 
         if stream_mode is None and file_path is not None:
             file_path = str(Path(file_path).resolve())

From 8888b60717cf1dd91078b1278b27063c87d1c616 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 8 Jul 2026 13:36:19 +0100
Subject: [PATCH 07/46] more performance fixes

---
 src/spikeinterface/core/sortinganalyzer.py    |  16 ++-
 .../extractors/nwbextractors.py               | 136 +++++++++++++-----
 2 files changed, 111 insertions(+), 41 deletions(-)

diff --git a/src/spikeinterface/core/sortinganalyzer.py b/src/spikeinterface/core/sortinganalyzer.py
index a1521390c7..5be9b67496 100644
--- a/src/spikeinterface/core/sortinganalyzer.py
+++ b/src/spikeinterface/core/sortinganalyzer.py
@@ -529,7 +529,9 @@ def load(cls, folder, recording=None, load_extensions=True, format="auto", backe
         return sorting_analyzer
 
     @classmethod
-    def create_memory(cls, sorting, recording, sparsity, return_in_uV, peak_sign, peak_mode, rec_attributes):
+    def create_memory(
+        cls, sorting, recording, sparsity, return_in_uV, peak_sign, peak_mode, rec_attributes, copy_sorting=True
+    ):
         # used by create and save_as
 
         if rec_attributes is None:
@@ -540,11 +542,17 @@ def create_memory(cls, sorting, recording, sparsity, return_in_uV, peak_sign, pe
             # a copy is done to avoid shared dict between instances (which can block garbage collector)
             rec_attributes = rec_attributes.copy()
 
-        # a copy of sorting is copied in memory for fast access
-        sorting_copy = NumpySorting.from_sorting(sorting, with_metadata=True, copy_spike_vector=True)
+        if copy_sorting:
+            # a copy of sorting is materialized in memory for fast access
+            analyzer_sorting = NumpySorting.from_sorting(sorting, with_metadata=True, copy_spike_vector=True)
+        else:
+            # keep the given (possibly lazy) sorting as-is: its spike times are read on demand rather
+            # than materialized up front. Useful for large streamed sortings where a view may need only
+            # a few units' trains (e.g. curation from stored templates + metrics).
+            analyzer_sorting = sorting
 
         sorting_analyzer = SortingAnalyzer(
-            sorting=sorting_copy,
+            sorting=analyzer_sorting,
             recording=recording,
             rec_attributes=rec_attributes,
             format="memory",
diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 42edc54b78..0d2b0e4f77 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1859,38 +1859,57 @@ def read_nwb_as_analyzer(
         use_pynwb=use_pynwb,
         t_start=t_start_tmp,
         sampling_frequency=sampling_frequency,
+        load_unit_properties=False,  # columns are read deliberately below, into their extensions
     )
 
+    sorting = sorting_tmp
     if recording is None and t_start is None:
-        # re-estimate t_start from spike times
-        if verbose:
-            print("Re-estimating t_start from spike_times")
-        t_start_new = np.min(sorting_tmp._sorting_segments[0].spike_times_data) - 0.001
+        # Estimate t_start from the first stored spike time and set it on the existing sorting in place.
+        # spike_times is stored per unit in ascending order, so the very first element is a good proxy
+        # for the earliest spike. Reading only that one element (instead of np.min over the whole lazy
+        # dataset, or constructing the sorting a second time) avoids streaming the entire spike_times
+        # array just to fix one scalar.
+        t_start_new = float(np.asarray(sorting._sorting_segments[0].spike_times_data[0])) - 0.001
         if verbose:
             print(f"Found new t_start: {t_start_new} s")
-        sorting = NwbSortingExtractor(
-            file_path=file_path,
-            electrical_series_path=electrical_series_path,
-            unit_table_path=unit_table_path,
-            stream_mode=stream_mode,
-            stream_cache_path=stream_cache_path,
-            cache=cache,
-            storage_options=storage_options,
-            use_pynwb=use_pynwb,
-            t_start=t_start_new,
-            sampling_frequency=sampling_frequency,
-        )
-    else:
-        sorting = sorting_tmp
-
+        for segment in sorting._sorting_segments:
+            segment._t_start = t_start_new
+
+    # Read the Units table deliberately. Classify each column and only materialize the ones that
+    # become part of the analyzer: `waveform_mean` (templates), the `electrodes` region (sparsity), the
+    # per-unit numeric metric columns (quality/template metrics), and the per-unit label columns
+    # (identity properties). The large per-spike columns (spike amplitudes, spike depths, each as long
+    # as the whole spike vector) are never read here, which is what makes streaming this table viable.
     if use_pynwb:
-        units = sorting.units_table
-        colnames = units.colnames
-        units = units.to_dataframe(index=True)
+        units_table = sorting.units_table
+        colnames = list(units_table.colnames)
+        units = units_table.to_dataframe(index=True)
+        structural = {"waveform_mean", "electrodes"}
+        metric_colnames = [c for c in units.columns if c not in structural and units[c].dtype.kind in "fiu"]
+        label_colnames = [
+            c
+            for c in units.columns
+            if c not in structural
+            and units[c].dtype.kind == "O"
+            and not isinstance(units[c].iloc[0], (list, np.ndarray))
+        ]
     else:
-        units_dset = sorting.units_table
-        units = _create_df_from_nwb_table(units_dset)
-        colnames = units.columns
+        units_group = sorting.units_table
+        colnames = list(units_group.keys())
+        structural = ("id", "spike_times", "waveform_mean", "electrodes")
+        metric_colnames, label_colnames = [], []
+        for column_name in colnames:
+            if column_name.endswith("_index") or column_name in structural or f"{column_name}_index" in colnames:
+                continue  # index columns, structural columns, and ragged/per-spike columns
+            dataset = units_group[column_name]
+            if dataset.attrs.get("neurodata_type") == "DynamicTableRegion":
+                continue  # e.g. max_electrode: a per-unit reference into electrodes, not a metric
+            if dataset.ndim == 1 and dataset.dtype.kind in "fiu":
+                metric_colnames.append(column_name)
+            elif dataset.ndim == 1 and dataset.dtype.kind in "OSU":
+                label_colnames.append(column_name)
+        needed = [c for c in ("waveform_mean", "electrodes") if c in colnames] + metric_colnames + label_colnames
+        units = _create_df_from_nwb_table(units_group, columns=needed)
 
     electrodes_indices = None
     if use_pynwb:
@@ -1919,7 +1938,10 @@ def read_nwb_as_analyzer(
                     units = units.loc[units.index[unit_mask]]
                     electrodes_indices = units["electrodes"]
 
-    # TODO: figure out sparsity
+    # per-unit label / identity columns (e.g. kilosort2_label, cluster_uuid) -> sorting properties
+    for label_column in label_colnames:
+        if label_column in units.columns and len(units[label_column]) == sorting.get_num_units():
+            sorting.set_property(label_column, units[label_column].values)
 
     # handle recording if available
     if recording is not None:
@@ -1944,7 +1966,11 @@ def read_nwb_as_analyzer(
             channel_ids = electrodes_table_sliced["channel_name"][:]
         else:
             channel_ids = electrodes_table_sliced["id"][:]
-        num_samples = [sorting.to_spike_vector()[-1]["sample_index"]]
+        # num_samples only needs to bound the recording length. Estimate it from the last stored spike
+        # time (one element read) rather than sorting.to_spike_vector(), which would materialize every
+        # spike just to read the maximum sample index.
+        last_spike_time = float(np.asarray(sorting._sorting_segments[0].spike_times_data[-1]))
+        num_samples = [int(last_spike_time * sorting.sampling_frequency) + 1]
         rec_attributes = dict(
             channel_ids=channel_ids,
             sampling_frequency=sorting.sampling_frequency,
@@ -2006,15 +2032,37 @@ def read_nwb_as_analyzer(
         peak_sign="neg",
         peak_mode="extremum",
         rec_attributes=rec_attributes,
+        # recordingless: keep the lazy sorting so spike_times is read on demand, not materialized here
+        copy_sorting=recording is not None,
     )
 
     # templates
     if "waveform_mean" in units and unit_local_channels is not None:
         from spikeinterface.core.analyzer_extension_core import ComputeTemplates
 
-        # random spikes is a dependency for templates; the raw spike samples are not stored so we
-        # select all of them
-        analyzer.compute("random_spikes", method="all")
+        # random_spikes is the analyzer's required extension. With a recording it drives waveform
+        # extraction, so compute it for real. Without a recording it is never used to extract anything
+        # (there are no traces); computing it would call sorting.to_spike_vector(), which materializes
+        # the whole spike_times array (e.g. ~165 MB for a full IBL session) and dominates the cost.
+        # Instead build it from the per-unit spike counts (the spike_times index), so spike_times is
+        # not read here; the lazy sorting still reads a unit's spikes on demand when a view needs them.
+        if recording is not None:
+            analyzer.compute("random_spikes", method="all")
+        else:
+            from spikeinterface.core.analyzer_extension_core import ComputeRandomSpikes
+
+            spike_times_index = sorting._sorting_segments[0].spike_times_index_data
+            total_num_spikes = int(np.asarray(spike_times_index[-1]))
+            max_spikes_per_unit = 500
+            num_random = min(max_spikes_per_unit * sorting.get_num_units(), total_num_spikes)
+            random_spikes_indices = np.sort(
+                np.random.default_rng(0).choice(total_num_spikes, size=num_random, replace=False)
+            )
+            random_spikes_ext = ComputeRandomSpikes(sorting_analyzer=analyzer)
+            random_spikes_ext.set_params(method="uniform", max_spikes_per_unit=max_spikes_per_unit, seed=0)
+            random_spikes_ext.data["random_spikes_indices"] = random_spikes_indices
+            random_spikes_ext.run_info["run_completed"] = True
+            analyzer.extensions["random_spikes"] = random_spikes_ext
 
         waveform_mean = np.array([np.asarray(t, dtype="float") for t in units["waveform_mean"].values])
         num_samples_template = waveform_mean.shape[1]
@@ -2043,15 +2091,20 @@ def read_nwb_as_analyzer(
         analyzer.extensions["templates"] = templates_ext
 
     template_metric_columns = ComputeTemplateMetrics.get_metric_columns()
-    quality_metric_columns = ComputeQualityMetrics.get_metric_columns()
 
     template_metrics_df = pd.DataFrame(index=sorting.unit_ids)
     quality_metric_df = pd.DataFrame(index=sorting.unit_ids)
 
-    for col in units.columns:
+    # Route every per-unit numeric column to its metrics extension. Canonical template-metric names go
+    # to template_metrics; everything else goes to quality_metrics, including pipeline-specific metrics
+    # that are not canonical SpikeInterface names (e.g. IBL's ibl_quality_score, noise_cutoff), so they
+    # are preserved and shown in the metrics view rather than dropped.
+    for col in metric_colnames:
+        if col not in units.columns:
+            continue
         if col in template_metric_columns:
             template_metrics_df.loc[:, col] = units[col].values
-        if col in quality_metric_columns:
+        else:
             quality_metric_df.loc[:, col] = units[col].values
 
     if len(template_metrics_df.columns) > 0:
@@ -2112,17 +2165,26 @@ def _create_dummy_probegroup_from_locations(locations, shape="circle", shape_par
     return probegroup
 
 
-def _create_df_from_nwb_table(group):
-    """Makes pandas DataFrame from hdf5/zarr NWB group"""
+def _create_df_from_nwb_table(group, columns=None):
+    """Makes pandas DataFrame from hdf5/zarr NWB group.
+
+    If `columns` is given, only those columns (plus the `id` index) are read; the data of every other
+    column is never touched. This lets callers avoid materializing large columns, e.g. the per-spike
+    arrays of a Units table (spike amplitudes, spike depths), when only a few columns are needed.
+    """
     import pandas as pd
 
-    colnames = list(group.keys())
+    all_colnames = list(group.keys())
+    if columns is None:
+        colnames = all_colnames
+    else:
+        colnames = [c for c in ["id", *columns] if c in all_colnames]
     data = {}
     for col in colnames:
         if "_index" in col:
             continue
         item = group[col][:]
-        if f"{col}_index" in colnames:
+        if f"{col}_index" in all_colnames:
             item = np.split(item, group[f"{col}_index"][:])[:-1]
             data[col] = item
         elif item.ndim > 1:

From b15964b468b306c412293bd7e5eeff74ddb2dc23 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 8 Jul 2026 15:16:54 +0100
Subject: [PATCH 08/46] spike vector

---
 .../extractors/nwbextractors.py               | 58 +++++++++++++++----
 1 file changed, 48 insertions(+), 10 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index acdfd793fb..6067b88101 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -8,6 +8,7 @@
 
 from spikeinterface import get_global_tmp_folder
 from spikeinterface.core import BaseRecording, BaseRecordingSegment, BaseSorting, BaseSortingSegment
+from spikeinterface.core.base import minimum_spike_dtype
 from spikeinterface.core.core_tools import define_function_from_class
 
 if importlib.util.find_spec("pynwb") is not None:
@@ -1426,6 +1427,25 @@ def __init__(
             "t_start": self.t_start,
         }
 
+    def _compute_and_cache_spike_vector(self) -> None:
+        # An NWB Units table is always single-segment. Spikes are stored as one flat times dataset
+        # grouped by unit, so the whole spike vector is built from a single bulk read instead of one
+        # streamed slice per unit (the generic BaseSorting implementation). This matters when streaming.
+        segment = self.segments[0]
+        sample_indices, unit_indices = segment._get_all_spike_frames_and_unit_indices()
+
+        # Stable sort by sample_index keeps unit_index ordering on ties, matching the way the generic
+        # BaseSorting builder constructs the vector (unit order within a sample).
+        order = np.argsort(sample_indices, stable=True)
+        n = sample_indices.size
+
+        spikes = np.zeros(n, dtype=minimum_spike_dtype)
+        spikes["sample_index"] = sample_indices[order]
+        spikes["unit_index"] = unit_indices[order]
+        # segment_index stays 0 for the single segment.
+        self._cached_spike_vector = spikes
+        self._cached_spike_vector_segment_slices = np.array([[0, n]], dtype="int64")
+
     @staticmethod
     def fetch_available_units_tables(
         file_path: str | Path,
@@ -1468,6 +1488,18 @@ def _frame_to_time(self, frame):
             return self._time_vector[frame]
         return frame / self._sampling_frequency + self._t_start
 
+    def _times_to_frames(self, spike_times) -> np.ndarray:
+        # Map seconds to samples. An irregular clock (per-sample timestamps) maps exactly with a
+        # searchsorted into the time vector; a uniform clock uses frame = (t - t_start) * sampling_frequency.
+        if self._time_vector is not None:
+            # searchsorted needs the timestamps in memory: materialize them once, on first use, and cache.
+            if not isinstance(self._time_vector, np.ndarray):
+                self._time_vector = np.asarray(self._time_vector)
+            frames = np.searchsorted(self._time_vector, spike_times, side="right") - 1
+        else:
+            frames = np.round((spike_times - self._t_start) * self._sampling_frequency)
+        return frames.astype("int64", copy=False)
+
     def get_unit_spike_train(
         self,
         unit_id,
@@ -1481,16 +1513,22 @@ def get_unit_spike_train(
         # Get spike times in seconds
         spike_times = self.get_unit_spike_train_in_seconds(unit_id=unit_id, start_time=start_time, end_time=end_time)
 
-        # Map seconds to samples. An irregular clock (per-sample timestamps) maps exactly with a
-        # searchsorted into the time vector; a uniform clock uses frame = (t - t_start) * sampling_frequency.
-        if self._time_vector is not None:
-            # searchsorted needs the timestamps in memory: materialize them once, on first use, and cache.
-            if not isinstance(self._time_vector, np.ndarray):
-                self._time_vector = np.asarray(self._time_vector)
-            frames = np.searchsorted(self._time_vector, spike_times, side="right") - 1
-        else:
-            frames = np.round((spike_times - self._t_start) * self._sampling_frequency)
-        return frames.astype("int64", copy=False)
+        return self._times_to_frames(spike_times)
+
+    def _get_all_spike_frames_and_unit_indices(self):
+        """Read all units' spike trains (as frames) in one bulk read, plus each spike's unit index.
+
+        NWB stores all spike times in a single flat dataset grouped by unit (with per-unit
+        boundaries in ``spike_times_index_data``), so the whole segment is read at once instead
+        of one streamed slice per unit.
+        """
+        spike_times = np.asarray(self.spike_times_data[:], dtype="float64")
+        # boundaries[u]:boundaries[u + 1] is the slice of spike_times belonging to unit u.
+        boundaries = np.concatenate(([0], np.asarray(self.spike_times_index_data[:], dtype="int64")))
+        counts = np.diff(boundaries)
+        unit_indices = np.repeat(np.arange(counts.size, dtype="int64"), counts)
+        sample_indices = self._times_to_frames(spike_times)
+        return sample_indices, unit_indices
 
     def get_unit_spike_train_in_seconds(
         self,

From a90b018f3ec04420c4afe61a6e7a36d0c4b626df Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 8 Jul 2026 15:25:45 +0100
Subject: [PATCH 09/46] Spike vector alterantive construction on nwb sorting

---
 src/spikeinterface/extractors/nwbextractors.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 6067b88101..0eb5550d8f 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1428,9 +1428,14 @@ def __init__(
         }
 
     def _compute_and_cache_spike_vector(self) -> None:
-        # An NWB Units table is always single-segment. Spikes are stored as one flat times dataset
-        # grouped by unit, so the whole spike vector is built from a single bulk read instead of one
-        # streamed slice per unit (the generic BaseSorting implementation). This matters when streaming.
+        # Performance: this deliberately overrides the generic BaseSorting builder, and the override is
+        # the point. The generic builder reads one spike train per unit (one streamed read per unit),
+        # whereas an NWB Units table stores every unit's spikes in a single flat times dataset grouped by
+        # unit. Reading that dataset once (see NwbSortingSegment._get_all_spike_frames_and_unit_indices)
+        # turns num_units streamed reads into one, which is a large difference for remote/streamed files.
+        # Do not drop this override in favor of the generic path without a reason: doing so silently
+        # reintroduces the per-unit read pattern. An NWB Units table is always single-segment, so the
+        # single segment is built directly.
         segment = self.segments[0]
         sample_indices, unit_indices = segment._get_all_spike_frames_and_unit_indices()
 

From 66af315f0a1075a31c2f548fdc888f6c0c91ff07 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 8 Jul 2026 15:47:42 +0100
Subject: [PATCH 10/46] nwb extractors

---
 .../extractors/nwbextractors.py               | 345 ++++++++++--------
 1 file changed, 197 insertions(+), 148 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 979942d94b..b0b80e725f 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1884,10 +1884,6 @@ def read_nwb_as_analyzer(
     compute_extra_params: dict | None = None,
     verbose: bool = False,
 ) -> SortingAnalyzer:
-    import pandas as pd
-    from spikeinterface.metrics.template import ComputeTemplateMetrics
-    from spikeinterface.metrics.quality import ComputeQualityMetrics
-
     # try to read recording object to get the analyzer
     try:
         recording = NwbRecordingExtractor(
@@ -2001,162 +1997,225 @@ def read_nwb_as_analyzer(
         if label_column in units.columns and len(units[label_column]) == sorting.get_num_units():
             sorting.set_property(label_column, units[label_column].values)
 
-    # handle recording if available
+    # Obtain a recording for the analyzer. With a real recording we use it directly. Without one (the
+    # streaming curation case) we build a lightweight placeholder recording that carries only the probe
+    # geometry, hand it to the standard constructor, and drop it afterwards. This mirrors
+    # `read_kilosort_as_analyzer` and avoids hand-assembling `rec_attributes` and channel ids.
     if recording is not None:
-        # check groups
         group_names = np.unique(recording.get_channel_groups())
         if group_name is not None and len(group_names) > 1:
             recording = recording.split_by("group")[group_name]
-        rec_attributes = None
+        analyzer_recording = recording
+        analyzer_channel_ids = list(recording.get_channel_ids())
     else:
-        recording = None
-        rec_attributes = {}
-
-        # get sliced electrodes table from electrode_indices union
-        electrode_indices_all = []
-        for ei in electrodes_indices:
-            electrode_indices_all.extend(ei)
-        electrode_indices_all = np.sort(np.unique(electrode_indices_all))
-        if verbose:
-            print(f"Found {len(electrode_indices_all)} electrodes")
-        electrodes_table_sliced = electrodes_table.iloc[electrode_indices_all]
-        if "channel_name" in electrodes_table_sliced:
-            channel_ids = electrodes_table_sliced["channel_name"][:]
-        else:
-            channel_ids = electrodes_table_sliced["id"][:]
-        # num_samples only needs to bound the recording length. Estimate it from the last stored spike
-        # time (one element read) rather than sorting.to_spike_vector(), which would materialize every
-        # spike just to read the maximum sample index.
-        last_spike_time = float(np.asarray(sorting._sorting_segments[0].spike_times_data[-1]))
-        num_samples = [int(last_spike_time * sorting.sampling_frequency) + 1]
-        rec_attributes = dict(
-            channel_ids=channel_ids,
-            sampling_frequency=sorting.sampling_frequency,
-            num_channels=len(channel_ids),
-            num_samples=num_samples,
+        analyzer_recording, analyzer_channel_ids = _make_placeholder_recording_from_electrodes(
+            sorting, electrodes_table, electrodes_indices, verbose=verbose
         )
-        # make a probegroup
-        electrode_colnames = electrodes_table_sliced.columns
-        assert (
-            "rel_x" in electrode_colnames and "rel_y" in electrode_colnames
-        ), "'rel_x' and 'rel_y' should be columns in the electrode name"
-        locations = np.array([electrodes_table_sliced["rel_x"][:], electrodes_table_sliced["rel_y"][:]]).T
-        probegroup = _create_dummy_probegroup_from_locations(locations)
-        rec_attributes["probegroup"] = probegroup
-
-    # Per-unit sparsity and channel map from the Units `electrodes` region. Each unit's
-    # `waveform_mean` is stored only on a subset of channels (near its peak); the region gives
-    # which channels those are. We use it both for the analyzer sparsity and to scatter the
-    # waveforms onto their true channel positions instead of stacking the sparse block densely.
-    sparsity = None
-    unit_local_channels = None
-    if electrodes_indices is not None:
-        from spikeinterface.core.sparsity import ChannelSparsity
-
-        if recording is not None:
-            analyzer_channel_ids = list(recording.get_channel_ids())
-        else:
-            analyzer_channel_ids = list(channel_ids)
-        num_channels = len(analyzer_channel_ids)
 
-        if "channel_name" in electrodes_table.columns:
-            electrode_row_to_channel_id = electrodes_table["channel_name"].to_numpy()
-        else:
-            electrode_row_to_channel_id = electrodes_table.index.to_numpy()
-        position_of_channel_id = {channel_id: pos for pos, channel_id in enumerate(analyzer_channel_ids)}
-
-        unit_local_channels = []
-        sparsity_mask = np.zeros((sorting.get_num_units(), num_channels), dtype=bool)
-        for unit_index, region in enumerate(electrodes_indices):
-            positions = [
-                position_of_channel_id[electrode_row_to_channel_id[int(electrode_row)]]
-                for electrode_row in region
-                if electrode_row_to_channel_id[int(electrode_row)] in position_of_channel_id
-            ]
-            unit_local_channels.append(positions)
-            sparsity_mask[unit_index, positions] = True
-        sparsity = ChannelSparsity(
-            mask=sparsity_mask,
-            unit_ids=np.asarray(sorting.unit_ids),
-            channel_ids=np.asarray(analyzer_channel_ids),
-        )
+    # Per-unit sparsity and channel map from the Units `electrodes` region. Each unit's `waveform_mean`
+    # is stored only on a subset of channels (near its peak); the region gives which channels those are.
+    # We use it both for the analyzer sparsity and to scatter the waveforms onto their true channel
+    # positions instead of stacking the sparse block densely.
+    sparsity, unit_local_channels = _make_sparsity_from_electrodes(
+        sorting, electrodes_table, electrodes_indices, analyzer_channel_ids
+    )
 
-    # instantiate analyzer
+    # Instantiate the analyzer through the standard constructor. For the recordingless case keep the
+    # lazy sorting (copy_sorting=False) so spike_times is read on demand rather than materialized here,
+    # and drop the placeholder recording once its geometry has been captured into rec_attributes.
     analyzer = SortingAnalyzer.create_memory(
         sorting=sorting,
-        recording=recording,
+        recording=analyzer_recording,
         sparsity=sparsity,
         return_in_uV=True,
         peak_sign="neg",
         peak_mode="extremum",
-        rec_attributes=rec_attributes,
-        # recordingless: keep the lazy sorting so spike_times is read on demand, not materialized here
+        rec_attributes=None,
         copy_sorting=recording is not None,
     )
+    if recording is None:
+        analyzer._recording = None
 
-    # templates
+    num_channels = len(analyzer_channel_ids)
+
+    # templates, and the required random_spikes extension that must exist alongside them
     if "waveform_mean" in units and unit_local_channels is not None:
-        from spikeinterface.core.analyzer_extension_core import ComputeTemplates
-
-        # random_spikes is the analyzer's required extension. With a recording it drives waveform
-        # extraction, so compute it for real. Without a recording it is never used to extract anything
-        # (there are no traces); computing it would call sorting.to_spike_vector(), which materializes
-        # the whole spike_times array (e.g. ~165 MB for a full IBL session) and dominates the cost.
-        # Instead build it from the per-unit spike counts (the spike_times index), so spike_times is
-        # not read here; the lazy sorting still reads a unit's spikes on demand when a view needs them.
-        if recording is not None:
-            analyzer.compute("random_spikes", method="all")
-        else:
-            from spikeinterface.core.analyzer_extension_core import ComputeRandomSpikes
-
-            spike_times_index = sorting._sorting_segments[0].spike_times_index_data
-            total_num_spikes = int(np.asarray(spike_times_index[-1]))
-            max_spikes_per_unit = 500
-            num_random = min(max_spikes_per_unit * sorting.get_num_units(), total_num_spikes)
-            random_spikes_indices = np.sort(
-                np.random.default_rng(0).choice(total_num_spikes, size=num_random, replace=False)
-            )
-            random_spikes_ext = ComputeRandomSpikes(sorting_analyzer=analyzer)
-            random_spikes_ext.set_params(method="uniform", max_spikes_per_unit=max_spikes_per_unit, seed=0)
-            random_spikes_ext.data["random_spikes_indices"] = random_spikes_indices
-            random_spikes_ext.run_info["run_completed"] = True
-            analyzer.extensions["random_spikes"] = random_spikes_ext
-
-        waveform_mean = np.array([np.asarray(t, dtype="float") for t in units["waveform_mean"].values])
-        num_samples_template = waveform_mean.shape[1]
-
-        # scatter each unit's real waveform columns onto their channels; the NaN-padded columns of
-        # units with fewer channels are simply never touched
-        dense_templates = np.zeros((sorting.get_num_units(), num_samples_template, num_channels), dtype="float32")
-        for unit_index, positions in enumerate(unit_local_channels):
-            k = len(positions)
-            dense_templates[unit_index][:, positions] = waveform_mean[unit_index][:, :k]
-
-        # scalar alignment (nbefore): the sample where the average peak-channel amplitude is largest
-        peak_amplitude_per_sample = np.abs(np.nan_to_num(waveform_mean)).max(axis=2).mean(axis=0)
-        nbefore = int(np.argmax(peak_amplitude_per_sample))
-
-        templates_ext = ComputeTemplates(sorting_analyzer=analyzer)
-        templates_ext.set_params(
-            ms_before=nbefore / analyzer.sampling_frequency * 1000,
-            ms_after=(num_samples_template - nbefore) / analyzer.sampling_frequency * 1000,
-            operators=["average", "std"],
-        )
-        templates_ext.data["average"] = dense_templates
-        templates_ext.data["std"] = np.zeros_like(dense_templates)
-        templates_ext.run_info["run_completed"] = True
+        _make_random_spikes(analyzer, sorting, has_recording=recording is not None)
+        _make_templates(analyzer, units, unit_local_channels, num_channels)
 
-        analyzer.extensions["templates"] = templates_ext
+    # per-unit numeric columns -> template_metrics / quality_metrics
+    _make_metrics(analyzer, units, metric_colnames, verbose=verbose)
 
-    template_metric_columns = ComputeTemplateMetrics.get_metric_columns()
+    # compute extra required
+    if compute_extra is not None:
+        if verbose:
+            print(f"Computing extra extensions: {compute_extra}")
+        compute_extra_params = {} if compute_extra_params is None else compute_extra_params
+        analyzer.compute(compute_extra, **compute_extra_params)
+
+    return analyzer
+
+
+def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, electrodes_indices, verbose=False):
+    """Build a lightweight placeholder recording that carries only the probe geometry of the electrodes
+    a unit's waveforms live on, for the recordingless (streaming) case.
 
-    template_metrics_df = pd.DataFrame(index=sorting.unit_ids)
-    quality_metric_df = pd.DataFrame(index=sorting.unit_ids)
+    Mirrors `read_kilosort_as_analyzer`: `generate_ground_truth_recording` produces a fully lazy
+    recording (noise generated on the fly, no traces materialized) whose only purpose is to feed probe
+    geometry, channel ids, and a bounding length into the standard analyzer constructor. The caller
+    drops the recording (`analyzer._recording = None`) right after construction. Returns the recording
+    and the ordered channel ids.
+    """
+    from probeinterface import Probe
+    from spikeinterface.core import generate_ground_truth_recording
+
+    # union of electrode rows referenced by any unit, in sorted order
+    electrode_indices_all = []
+    for region in electrodes_indices:
+        electrode_indices_all.extend(region)
+    electrode_indices_all = np.sort(np.unique(electrode_indices_all))
+    if verbose:
+        print(f"Found {len(electrode_indices_all)} electrodes")
+    electrodes_table_sliced = electrodes_table.iloc[electrode_indices_all]
+
+    if "channel_name" in electrodes_table_sliced:
+        channel_ids = np.asarray(electrodes_table_sliced["channel_name"][:])
+    else:
+        channel_ids = electrodes_table_sliced.index.to_numpy()
+
+    electrode_colnames = electrodes_table_sliced.columns
+    assert (
+        "rel_x" in electrode_colnames and "rel_y" in electrode_colnames
+    ), "'rel_x' and 'rel_y' should be columns in the electrodes table"
+    locations = np.array([electrodes_table_sliced["rel_x"][:], electrodes_table_sliced["rel_y"][:]]).T
+
+    # The recording length only needs to bound the spikes. Estimate it from the last stored spike time
+    # (one element read) rather than sorting.to_spike_vector(), which would materialize every spike.
+    # spike_times is stored in seconds, so this is already a duration in seconds (no sampling-rate
+    # division, unlike the Kilosort path where spike times are in samples).
+    last_spike_time = float(np.asarray(sorting._sorting_segments[0].spike_times_data[-1]))
+    duration = last_spike_time + 1.0
+
+    probe = Probe(si_units="um")
+    probe.set_contacts(locations, shapes="circle", shape_params={"radius": 1})
+    probe.set_device_channel_indices(np.arange(len(channel_ids)))
+
+    recording, _ = generate_ground_truth_recording(
+        durations=[duration],
+        sampling_frequency=sorting.sampling_frequency,
+        probe=probe,
+        num_units=1,
+        seed=0,
+    )
+    recording = recording.rename_channels(channel_ids)
+    return recording, list(channel_ids)
+
+
+def _make_sparsity_from_electrodes(sorting, electrodes_table, electrodes_indices, analyzer_channel_ids):
+    """Build the analyzer `ChannelSparsity` and the per-unit local channel lists from the Units
+    `electrodes` region, mapping each unit's electrode rows onto positions in `analyzer_channel_ids`.
+    Returns (sparsity, unit_local_channels), or (None, None) if there is no electrodes region."""
+    if electrodes_indices is None:
+        return None, None
+
+    from spikeinterface.core.sparsity import ChannelSparsity
+
+    num_channels = len(analyzer_channel_ids)
+    if "channel_name" in electrodes_table.columns:
+        electrode_row_to_channel_id = electrodes_table["channel_name"].to_numpy()
+    else:
+        electrode_row_to_channel_id = electrodes_table.index.to_numpy()
+    position_of_channel_id = {channel_id: pos for pos, channel_id in enumerate(analyzer_channel_ids)}
+
+    unit_local_channels = []
+    sparsity_mask = np.zeros((sorting.get_num_units(), num_channels), dtype=bool)
+    for unit_index, region in enumerate(electrodes_indices):
+        positions = [
+            position_of_channel_id[electrode_row_to_channel_id[int(electrode_row)]]
+            for electrode_row in region
+            if electrode_row_to_channel_id[int(electrode_row)] in position_of_channel_id
+        ]
+        unit_local_channels.append(positions)
+        sparsity_mask[unit_index, positions] = True
+    sparsity = ChannelSparsity(
+        mask=sparsity_mask,
+        unit_ids=np.asarray(sorting.unit_ids),
+        channel_ids=np.asarray(analyzer_channel_ids),
+    )
+    return sparsity, unit_local_channels
+
+
+def _make_random_spikes(analyzer, sorting, has_recording):
+    """Attach the required `random_spikes` extension.
+
+    With a recording it drives waveform extraction, so compute it for real. Without a recording it is
+    never used to extract anything (there are no traces); computing it would call
+    `sorting.to_spike_vector()`, which materializes the whole spike_times array (e.g. ~165 MB for a
+    full IBL session) and dominates the cost. Instead build it from the per-unit spike counts (the
+    spike_times index), so spike_times is not read here; the lazy sorting still reads a unit's spikes on
+    demand when a view needs them.
+    """
+    if has_recording:
+        analyzer.compute("random_spikes", method="all")
+        return
+
+    from spikeinterface.core.analyzer_extension_core import ComputeRandomSpikes
+
+    spike_times_index = sorting._sorting_segments[0].spike_times_index_data
+    total_num_spikes = int(np.asarray(spike_times_index[-1]))
+    max_spikes_per_unit = 500
+    num_random = min(max_spikes_per_unit * sorting.get_num_units(), total_num_spikes)
+    random_spikes_indices = np.sort(np.random.default_rng(0).choice(total_num_spikes, size=num_random, replace=False))
+    random_spikes_ext = ComputeRandomSpikes(sorting_analyzer=analyzer)
+    random_spikes_ext.set_params(method="uniform", max_spikes_per_unit=max_spikes_per_unit, seed=0)
+    random_spikes_ext.data["random_spikes_indices"] = random_spikes_indices
+    random_spikes_ext.run_info["run_completed"] = True
+    analyzer.extensions["random_spikes"] = random_spikes_ext
+
+
+def _make_templates(analyzer, units, unit_local_channels, num_channels):
+    """Attach the `templates` extension from the Units `waveform_mean` column, scattering each unit's
+    sparse waveform block onto its true channel positions."""
+    from spikeinterface.core.analyzer_extension_core import ComputeTemplates
+
+    waveform_mean = np.array([np.asarray(t, dtype="float") for t in units["waveform_mean"].values])
+    num_samples_template = waveform_mean.shape[1]
+
+    # scatter each unit's real waveform columns onto their channels; the NaN-padded columns of units
+    # with fewer channels are simply never touched
+    dense_templates = np.zeros((analyzer.get_num_units(), num_samples_template, num_channels), dtype="float32")
+    for unit_index, positions in enumerate(unit_local_channels):
+        k = len(positions)
+        dense_templates[unit_index][:, positions] = waveform_mean[unit_index][:, :k]
+
+    # scalar alignment (nbefore): the sample where the average peak-channel amplitude is largest
+    peak_amplitude_per_sample = np.abs(np.nan_to_num(waveform_mean)).max(axis=2).mean(axis=0)
+    nbefore = int(np.argmax(peak_amplitude_per_sample))
+
+    templates_ext = ComputeTemplates(sorting_analyzer=analyzer)
+    templates_ext.set_params(
+        ms_before=nbefore / analyzer.sampling_frequency * 1000,
+        ms_after=(num_samples_template - nbefore) / analyzer.sampling_frequency * 1000,
+        operators=["average", "std"],
+    )
+    templates_ext.data["average"] = dense_templates
+    templates_ext.data["std"] = np.zeros_like(dense_templates)
+    templates_ext.run_info["run_completed"] = True
+    analyzer.extensions["templates"] = templates_ext
 
-    # Route every per-unit numeric column to its metrics extension. Canonical template-metric names go
-    # to template_metrics; everything else goes to quality_metrics, including pipeline-specific metrics
-    # that are not canonical SpikeInterface names (e.g. IBL's ibl_quality_score, noise_cutoff), so they
-    # are preserved and shown in the metrics view rather than dropped.
+
+def _make_metrics(analyzer, units, metric_colnames, verbose=False):
+    """Route every per-unit numeric column to its metrics extension. Canonical template-metric names go
+    to `template_metrics`; everything else goes to `quality_metrics`, including pipeline-specific metrics
+    that are not canonical SpikeInterface names (e.g. IBL's noise_cutoff), so they are preserved and
+    shown in the metrics view rather than dropped."""
+    import pandas as pd
+    from spikeinterface.metrics.template import ComputeTemplateMetrics
+    from spikeinterface.metrics.quality import ComputeQualityMetrics
+
+    template_metric_columns = ComputeTemplateMetrics.get_metric_columns()
+    template_metrics_df = pd.DataFrame(index=analyzer.unit_ids)
+    quality_metric_df = pd.DataFrame(index=analyzer.unit_ids)
     for col in metric_colnames:
         if col not in units.columns:
             continue
@@ -2169,7 +2228,6 @@ def read_nwb_as_analyzer(
         if verbose:
             print("Adding template metrics")
         template_metrics_ext = ComputeTemplateMetrics(analyzer)
-        # cast to correct dtypes
         template_metrics_ext.data["metrics"] = template_metrics_ext._cast_metrics(template_metrics_df)
         template_metrics_ext.run_info["run_completed"] = True
         analyzer.extensions["template_metrics"] = template_metrics_ext
@@ -2181,15 +2239,6 @@ def read_nwb_as_analyzer(
         quality_metrics_ext.run_info["run_completed"] = True
         analyzer.extensions["quality_metrics"] = quality_metrics_ext
 
-    # compute extra required
-    if compute_extra is not None:
-        if verbose:
-            print(f"Computing extra extensions: {compute_extra}")
-        compute_extra_params = {} if compute_extra_params is None else compute_extra_params
-        analyzer.compute(compute_extra, **compute_extra_params)
-
-    return analyzer
-
 
 def _create_dummy_probegroup_from_locations(locations, shape="circle", shape_params={"radius": 1}):
     """

From c53e3366e6f06d2b9d2a0fea394f7b5aa06cbc04 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 8 Jul 2026 16:24:04 +0100
Subject: [PATCH 11/46] stream improvements

---
 .../extractors/nwbextractors.py               | 65 ++++---------------
 1 file changed, 14 insertions(+), 51 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index b0b80e725f..a04436cf06 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1880,7 +1880,7 @@ def read_nwb_as_analyzer(
     storage_options: dict | None = None,
     use_pynwb: bool = False,
     group_name: str | None = None,
-    compute_extra: List[str] | None = ["unit_locations", "correlograms"],
+    compute_extra: List[str] | None = ["unit_locations"],
     compute_extra_params: dict | None = None,
     verbose: bool = False,
 ) -> SortingAnalyzer:
@@ -1917,17 +1917,10 @@ def read_nwb_as_analyzer(
     )
 
     sorting = sorting_tmp
-    if recording is None and t_start is None:
-        # Estimate t_start from the first stored spike time and set it on the existing sorting in place.
-        # spike_times is stored per unit in ascending order, so the very first element is a good proxy
-        # for the earliest spike. Reading only that one element (instead of np.min over the whole lazy
-        # dataset, or constructing the sorting a second time) avoids streaming the entire spike_times
-        # array just to fix one scalar.
-        t_start_new = float(np.asarray(sorting._sorting_segments[0].spike_times_data[0])) - 0.001
-        if verbose:
-            print(f"Found new t_start: {t_start_new} s")
-        for segment in sorting._sorting_segments:
-            segment._t_start = t_start_new
+    # Recordingless case: leave t_start at 0 (set when the sorting was constructed). NWB spike times are
+    # in seconds from session start, so they are always >= 0 and map to non-negative frames with
+    # t_start = 0. Anchoring t_start to the first spike would only remove empty space before it on the
+    # timeline (cosmetic) and would cost a streamed read, so we skip it to keep the build cheap.
 
     # Read the Units table deliberately. Classify each column and only materialize the ones that
     # become part of the analyzer: `waveform_mean` (templates), the `electrodes` region (sparsity), the
@@ -2062,9 +2055,9 @@ def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, elect
 
     Mirrors `read_kilosort_as_analyzer`: `generate_ground_truth_recording` produces a fully lazy
     recording (noise generated on the fly, no traces materialized) whose only purpose is to feed probe
-    geometry, channel ids, and a bounding length into the standard analyzer constructor. The caller
-    drops the recording (`analyzer._recording = None`) right after construction. Returns the recording
-    and the ordered channel ids.
+    geometry, channel ids, and a bounding length into the standard analyzer constructor. The caller drops
+    the recording (`analyzer._recording = None`) right after construction. Returns the recording and the
+    ordered channel ids.
     """
     from probeinterface import Probe
     from spikeinterface.core import generate_ground_truth_recording
@@ -2089,10 +2082,12 @@ def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, elect
     ), "'rel_x' and 'rel_y' should be columns in the electrodes table"
     locations = np.array([electrodes_table_sliced["rel_x"][:], electrodes_table_sliced["rel_y"][:]]).T
 
-    # The recording length only needs to bound the spikes. Estimate it from the last stored spike time
-    # (one element read) rather than sorting.to_spike_vector(), which would materialize every spike.
-    # spike_times is stored in seconds, so this is already a duration in seconds (no sampling-rate
-    # division, unlike the Kilosort path where spike times are in samples).
+    # The recording length only needs to loosely bound the timeline for the recordingless GUI; nothing
+    # about curation depends on it being exact. Estimate it cheaply from the last stored spike time (one
+    # element, i.e. only the last spike_times chunk) rather than scanning the whole array for the true
+    # global maximum. spike_times is in seconds, so this is already a duration in seconds. It is only an
+    # approximation: because spike_times is concatenated per unit and not globally sorted, this is the
+    # last unit's last spike, not necessarily the latest spike overall.
     last_spike_time = float(np.asarray(sorting._sorting_segments[0].spike_times_data[-1]))
     duration = last_spike_time + 1.0
 
@@ -2240,38 +2235,6 @@ def _make_metrics(analyzer, units, metric_colnames, verbose=False):
         analyzer.extensions["quality_metrics"] = quality_metrics_ext
 
 
-def _create_dummy_probegroup_from_locations(locations, shape="circle", shape_params={"radius": 1}):
-    """
-    Creates a "dummy" probe based on locations.
-
-    Parameters
-    ----------
-    locations : np.array
-        Array with channel locations (num_channels, ndim) [ndim can be 2 or 3]
-    shape : str, default: "circle"
-        Electrode shapes
-    shape_params : dict, default: {"radius": 1}
-        Shape parameters
-
-    Returns
-    -------
-    probe : Probe
-        The created probe
-    """
-    from probeinterface import Probe, ProbeGroup
-
-    ndim = locations.shape[1]
-    assert ndim == 2
-    probe = Probe(ndim=2)
-    probe.set_contacts(locations, shapes=shape, shape_params=shape_params)
-    probe.set_device_channel_indices(np.arange(len(probe.contact_positions)))
-    probe.create_auto_shape()
-    probegroup = ProbeGroup()
-    probegroup.add_probe(probe)
-
-    return probegroup
-
-
 def _create_df_from_nwb_table(group, columns=None):
     """Makes pandas DataFrame from hdf5/zarr NWB group.
 

From 469dc9b6e0ac4336ff423ca25a708d5fd5448476 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 8 Jul 2026 17:05:39 +0100
Subject: [PATCH 12/46] fix channel indices

---
 .../extractors/nwbextractors.py               | 19 +++++++++++++++----
 1 file changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index a04436cf06..1e1bf3b913 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2072,7 +2072,12 @@ def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, elect
     electrodes_table_sliced = electrodes_table.iloc[electrode_indices_all]
 
     if "channel_name" in electrodes_table_sliced:
-        channel_ids = np.asarray(electrodes_table_sliced["channel_name"][:])
+        # IBL stores channel_name as variable-length HDF5 byte strings (e.g. b"AP0"); decode to str so
+        # SpikeInterface accepts them as channel ids (object-dtype arrays are rejected). This must match
+        # the decode in _make_sparsity_from_electrodes so channel ids and sparsity keys stay consistent.
+        channel_ids = np.array(
+            [c.decode("utf-8") if isinstance(c, bytes) else c for c in electrodes_table_sliced["channel_name"][:]]
+        )
     else:
         channel_ids = electrodes_table_sliced.index.to_numpy()
 
@@ -2117,7 +2122,10 @@ def _make_sparsity_from_electrodes(sorting, electrodes_table, electrodes_indices
 
     num_channels = len(analyzer_channel_ids)
     if "channel_name" in electrodes_table.columns:
-        electrode_row_to_channel_id = electrodes_table["channel_name"].to_numpy()
+        # decode HDF5 byte-string channel names to str, matching _make_placeholder_recording_from_electrodes
+        electrode_row_to_channel_id = np.array(
+            [c.decode("utf-8") if isinstance(c, bytes) else c for c in electrodes_table["channel_name"].to_numpy()]
+        )
     else:
         electrode_row_to_channel_id = electrodes_table.index.to_numpy()
     position_of_channel_id = {channel_id: pos for pos, channel_id in enumerate(analyzer_channel_ids)}
@@ -2176,8 +2184,11 @@ def _make_templates(analyzer, units, unit_local_channels, num_channels):
     waveform_mean = np.array([np.asarray(t, dtype="float") for t in units["waveform_mean"].values])
     num_samples_template = waveform_mean.shape[1]
 
-    # scatter each unit's real waveform columns onto their channels; the NaN-padded columns of units
-    # with fewer channels are simply never touched
+    # Scatter each unit's waveform onto its channels. This relies on the alignment contract that
+    # waveform_mean[:, i] corresponds to the unit's electrodes-region entry i. The IBL writer guarantees
+    # it by ordering each unit's real channels first (padding last) and building the electrodes region
+    # from exactly those real channels in the same order, so the first k = len(region) columns are the
+    # real channels in region order and any trailing padding columns are correctly dropped by [:, :k].
     dense_templates = np.zeros((analyzer.get_num_units(), num_samples_template, num_channels), dtype="float32")
     for unit_index, positions in enumerate(unit_local_channels):
         k = len(positions)

From 47f9fc782ea9fa415a9dc3d8a24ac96258bc4eea Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 9 Jul 2026 11:27:09 +0100
Subject: [PATCH 13/46] better std behavior

---
 .../extractors/nwbextractors.py               | 27 ++++++++++++++-----
 1 file changed, 21 insertions(+), 6 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 1e1bf3b913..005c655137 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1931,7 +1931,7 @@ def read_nwb_as_analyzer(
         units_table = sorting.units_table
         colnames = list(units_table.colnames)
         units = units_table.to_dataframe(index=True)
-        structural = {"waveform_mean", "electrodes"}
+        structural = {"waveform_mean", "waveform_sd", "electrodes"}
         metric_colnames = [c for c in units.columns if c not in structural and units[c].dtype.kind in "fiu"]
         label_colnames = [
             c
@@ -1943,7 +1943,7 @@ def read_nwb_as_analyzer(
     else:
         units_group = sorting.units_table
         colnames = list(units_group.keys())
-        structural = ("id", "spike_times", "waveform_mean", "electrodes")
+        structural = ("id", "spike_times", "waveform_mean", "waveform_sd", "electrodes")
         metric_colnames, label_colnames = [], []
         for column_name in colnames:
             if column_name.endswith("_index") or column_name in structural or f"{column_name}_index" in colnames:
@@ -1955,7 +1955,9 @@ def read_nwb_as_analyzer(
                 metric_colnames.append(column_name)
             elif dataset.ndim == 1 and dataset.dtype.kind in "OSU":
                 label_colnames.append(column_name)
-        needed = [c for c in ("waveform_mean", "electrodes") if c in colnames] + metric_colnames + label_colnames
+        needed = [
+            c for c in ("waveform_mean", "waveform_sd", "electrodes") if c in colnames
+        ] + metric_colnames + label_colnames
         units = _create_df_from_nwb_table(units_group, columns=needed)
 
     electrodes_indices = None
@@ -2199,13 +2201,26 @@ def _make_templates(analyzer, units, unit_local_channels, num_channels):
     nbefore = int(np.argmax(peak_amplitude_per_sample))
 
     templates_ext = ComputeTemplates(sorting_analyzer=analyzer)
+    templates_ext.data["average"] = dense_templates
+    operators = ["average"]
+
+    # Only expose the "std" operator when the file actually stores waveform_sd. Fabricating a zero std
+    # would falsely imply the unit's spikes have no waveform variability, so when it is absent we declare
+    # only "average" rather than inventing data.
+    if "waveform_sd" in units.columns:
+        waveform_sd = np.array([np.asarray(t, dtype="float") for t in units["waveform_sd"].values])
+        dense_std = np.zeros((analyzer.get_num_units(), num_samples_template, num_channels), dtype="float32")
+        for unit_index, positions in enumerate(unit_local_channels):
+            k = len(positions)
+            dense_std[unit_index][:, positions] = waveform_sd[unit_index][:, :k]
+        templates_ext.data["std"] = dense_std
+        operators.append("std")
+
     templates_ext.set_params(
         ms_before=nbefore / analyzer.sampling_frequency * 1000,
         ms_after=(num_samples_template - nbefore) / analyzer.sampling_frequency * 1000,
-        operators=["average", "std"],
+        operators=operators,
     )
-    templates_ext.data["average"] = dense_templates
-    templates_ext.data["std"] = np.zeros_like(dense_templates)
     templates_ext.run_info["run_completed"] = True
     analyzer.extensions["templates"] = templates_ext
 

From d010d13fd31076d7578654c569627b7ec4ff5fad Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 9 Jul 2026 10:52:38 +0000
Subject: [PATCH 14/46] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 src/spikeinterface/extractors/nwbextractors.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 005c655137..cb10c05a07 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1955,9 +1955,11 @@ def read_nwb_as_analyzer(
                 metric_colnames.append(column_name)
             elif dataset.ndim == 1 and dataset.dtype.kind in "OSU":
                 label_colnames.append(column_name)
-        needed = [
-            c for c in ("waveform_mean", "waveform_sd", "electrodes") if c in colnames
-        ] + metric_colnames + label_colnames
+        needed = (
+            [c for c in ("waveform_mean", "waveform_sd", "electrodes") if c in colnames]
+            + metric_colnames
+            + label_colnames
+        )
         units = _create_df_from_nwb_table(units_group, columns=needed)
 
     electrodes_indices = None

From 2c7ffa54c8ebafaa4d7448bd35169a653fcbfabd Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 9 Jul 2026 12:50:26 +0100
Subject: [PATCH 15/46] better naming

---
 src/spikeinterface/extractors/__init__.py      | 2 +-
 src/spikeinterface/extractors/nwbextractors.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/spikeinterface/extractors/__init__.py b/src/spikeinterface/extractors/__init__.py
index 42b382393a..58f7860d3d 100644
--- a/src/spikeinterface/extractors/__init__.py
+++ b/src/spikeinterface/extractors/__init__.py
@@ -5,7 +5,7 @@
 
 from .neoextractors import get_neo_num_blocks, get_neo_streams
 from .phykilosortextractors import read_kilosort_as_analyzer
-from .nwbextractors import read_nwb_as_analyzer
+from .nwbextractors import read_nwb_sorting_analyzer
 
 from warnings import warn
 
diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 005c655137..43a84ccfe2 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1868,7 +1868,7 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
     return outputs
 
 
-def read_nwb_as_analyzer(
+def read_nwb_sorting_analyzer(
     file_path: str | Path,
     t_start: float | None = None,
     sampling_frequency: float | None = None,

From bdbbe86fc4960034aa8a6b97f38a46e5c37ac3cf Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 9 Jul 2026 13:43:12 +0100
Subject: [PATCH 16/46] zach comments

---
 src/spikeinterface/core/sortinganalyzer.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/spikeinterface/core/sortinganalyzer.py b/src/spikeinterface/core/sortinganalyzer.py
index 8c92040dfd..1f7fa87348 100644
--- a/src/spikeinterface/core/sortinganalyzer.py
+++ b/src/spikeinterface/core/sortinganalyzer.py
@@ -1962,7 +1962,7 @@ def is_sparse(self) -> bool:
         return self.sparsity is not None
 
     def is_filtered(self) -> bool:
-        return self.rec_attributes.get("is_filtered", True)
+        return self.rec_attributes["is_filtered"]
 
     def get_sorting_provenance(self):
         """
@@ -2063,7 +2063,7 @@ def get_sorting_property(self, key) -> np.ndarray:
         return self.sorting.get_property(key)
 
     def get_dtype(self):
-        return self.rec_attributes.get("dtype")
+        return self.rec_attributes["dtype"]
 
     def get_num_units(self) -> int:
         return self.sorting.get_num_units()

From ee97e9c9bd1e9376167e1fcfef089f9cc0b1a889 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 9 Jul 2026 16:13:53 +0100
Subject: [PATCH 17/46] mapping

---
 .../extractors/nwbextractors.py               | 160 +++++++++++++-----
 1 file changed, 117 insertions(+), 43 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index a3361af5dd..a08c91f5d5 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1868,6 +1868,78 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
     return outputs
 
 
+# Where each SortingAnalyzer extension's data lives in an NWB file. Each extension maps to an ordered
+# list of read-sources tried until one resolves: the typed ndx-spikesorting container first, then the
+# generic-NWB convention. A user/dataset override merges over this (see `extension_map` in
+# `read_nwb_sorting_analyzer`); `None` disables an extension. This is only for real analyzer extensions;
+# sparsity, sorting properties, and recording metadata are handled separately by the reader.
+#
+# The "typed_container" source is a hook: the typed reader currently lives in ndx-spikesorting, so for
+# now it never resolves here and the convention source is used. When that logic moves into
+# spikeinterface, this source will resolve against the ndx-spikesorting types.
+DEFAULT_EXTENSION_MAP = {
+    "templates": [
+        {"source": "typed_container", "type": "Templates"},
+        {"source": "column", "column": "waveform_mean", "std_column": "waveform_sd"},
+    ],
+    "quality_metrics": [
+        {"source": "typed_container", "type": "UnitsMetrics"},
+        {"source": "columns", "columns": "canonical_quality"},
+    ],
+    "template_metrics": [
+        {"source": "typed_container", "type": "UnitsMetrics"},
+    ],
+    "unit_locations": [
+        {"source": "typed_container", "type": "UnitLocations"},
+    ],
+}
+
+
+def _resolve_extension_columns(extension_map, scalar_colnames, all_colnames):
+    """Resolve the convention (column-based) sources of an extension map against a Units table's columns.
+
+    Returns a dict with the resolved template column (and optional std column) and the per-extension
+    column lists for the column-backed extensions. The "typed_container" source is skipped here (the
+    typed reader lives in ndx-spikesorting for now), so resolution falls through to the "column"/
+    "columns" sources. `scalar_colnames` are the 1-D per-unit columns eligible to be metrics/properties;
+    `all_colnames` is every column name (used to check a source's target exists).
+    """
+    from spikeinterface.metrics.quality import ComputeQualityMetrics
+
+    resolved = {"templates_column": None, "templates_std_column": None, "quality_metric_colnames": []}
+
+    for extension, sources in extension_map.items():
+        if sources is None:  # explicitly disabled
+            continue
+        for source in sources:
+            kind = source.get("source")
+            if kind == "typed_container":
+                continue  # hook: not resolvable here yet
+            if kind == "column":
+                column = source.get("column")
+                if column in all_colnames:
+                    if extension == "templates":
+                        resolved["templates_column"] = column
+                        std_column = source.get("std_column")
+                        resolved["templates_std_column"] = std_column if std_column in all_colnames else None
+                    break
+            if kind == "columns":
+                spec = source.get("columns")
+                if spec == "canonical_quality":
+                    canonical = set(ComputeQualityMetrics.get_metric_columns())
+                    cols = [c for c in scalar_colnames if c in canonical]
+                else:  # an explicit list of column names
+                    cols = [c for c in spec if c in scalar_colnames]
+                if extension == "quality_metrics":
+                    resolved["quality_metric_colnames"] = cols
+                break
+
+    # every scalar column not claimed as a quality metric becomes a sorting property
+    claimed = set(resolved["quality_metric_colnames"])
+    resolved["property_colnames"] = [c for c in scalar_colnames if c not in claimed]
+    return resolved
+
+
 def read_nwb_sorting_analyzer(
     file_path: str | Path,
     t_start: float | None = None,
@@ -1882,8 +1954,13 @@ def read_nwb_sorting_analyzer(
     group_name: str | None = None,
     compute_extra: List[str] | None = ["unit_locations"],
     compute_extra_params: dict | None = None,
+    extension_map: dict | None = None,
     verbose: bool = False,
 ) -> SortingAnalyzer:
+    # extension_map overrides (per extension) merge over DEFAULT_EXTENSION_MAP; see its docstring.
+    resolved_extension_map = dict(DEFAULT_EXTENSION_MAP)
+    if extension_map is not None:
+        resolved_extension_map.update(extension_map)
     # try to read recording object to get the analyzer
     try:
         recording = NwbRecordingExtractor(
@@ -1955,11 +2032,21 @@ def read_nwb_sorting_analyzer(
                 metric_colnames.append(column_name)
             elif dataset.ndim == 1 and dataset.dtype.kind in "OSU":
                 label_colnames.append(column_name)
-        needed = (
-            [c for c in ("waveform_mean", "waveform_sd", "electrodes") if c in colnames]
-            + metric_colnames
-            + label_colnames
-        )
+
+    # Resolve the extension map against the Units columns: which column holds the templates, and which
+    # scalar columns are quality metrics (canonical names by default) vs plain sorting properties.
+    scalar_colnames = metric_colnames + label_colnames
+    resolved = _resolve_extension_columns(resolved_extension_map, scalar_colnames, colnames)
+    templates_column = resolved["templates_column"]
+    templates_std_column = resolved["templates_std_column"]
+    quality_metric_colnames = resolved["quality_metric_colnames"]
+    property_colnames = resolved["property_colnames"]
+
+    if not use_pynwb:
+        # deliberate read: only the resolved template column(s), the electrodes region, and the scalar
+        # columns that become metrics or properties are materialized; ragged/per-spike columns are not.
+        needed = [c for c in (templates_column, templates_std_column, "electrodes") if c and c in colnames]
+        needed += scalar_colnames
         units = _create_df_from_nwb_table(units_group, columns=needed)
 
     electrodes_indices = None
@@ -1989,10 +2076,12 @@ def read_nwb_sorting_analyzer(
                     units = units.loc[units.index[unit_mask]]
                     electrodes_indices = units["electrodes"]
 
-    # per-unit label / identity columns (e.g. kilosort2_label, cluster_uuid) -> sorting properties
-    for label_column in label_colnames:
-        if label_column in units.columns and len(units[label_column]) == sorting.get_num_units():
-            sorting.set_property(label_column, units[label_column].values)
+    # Every scalar column not claimed as a quality metric becomes a sorting property (labels like
+    # cluster_uuid, and any non-canonical numeric column). Properties show up in the GUI's unit table
+    # alongside metrics, so nothing is lost, without asserting an unknown column is a "quality metric".
+    for property_column in property_colnames:
+        if property_column in units.columns and len(units[property_column]) == sorting.get_num_units():
+            sorting.set_property(property_column, units[property_column].values)
 
     # Obtain a recording for the analyzer. With a real recording we use it directly. Without one (the
     # streaming curation case) we build a lightweight placeholder recording that carries only the probe
@@ -2035,13 +2124,13 @@ def read_nwb_sorting_analyzer(
 
     num_channels = len(analyzer_channel_ids)
 
-    # templates, and the required random_spikes extension that must exist alongside them
-    if "waveform_mean" in units and unit_local_channels is not None:
+    # templates (from the resolved templates column), and the required random_spikes extension
+    if templates_column is not None and templates_column in units and unit_local_channels is not None:
         _make_random_spikes(analyzer, sorting, has_recording=recording is not None)
-        _make_templates(analyzer, units, unit_local_channels, num_channels)
+        _make_templates(analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels)
 
-    # per-unit numeric columns -> template_metrics / quality_metrics
-    _make_metrics(analyzer, units, metric_colnames, verbose=verbose)
+    # the resolved quality-metric columns -> quality_metrics extension
+    _make_metrics(analyzer, units, quality_metric_colnames, verbose=verbose)
 
     # compute extra required
     if compute_extra is not None:
@@ -2180,12 +2269,12 @@ def _make_random_spikes(analyzer, sorting, has_recording):
     analyzer.extensions["random_spikes"] = random_spikes_ext
 
 
-def _make_templates(analyzer, units, unit_local_channels, num_channels):
-    """Attach the `templates` extension from the Units `waveform_mean` column, scattering each unit's
+def _make_templates(analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels):
+    """Attach the `templates` extension from the resolved templates column, scattering each unit's
     sparse waveform block onto its true channel positions."""
     from spikeinterface.core.analyzer_extension_core import ComputeTemplates
 
-    waveform_mean = np.array([np.asarray(t, dtype="float") for t in units["waveform_mean"].values])
+    waveform_mean = np.array([np.asarray(t, dtype="float") for t in units[templates_column].values])
     num_samples_template = waveform_mean.shape[1]
 
     # Scatter each unit's waveform onto its channels. This relies on the alignment contract that
@@ -2206,11 +2295,11 @@ def _make_templates(analyzer, units, unit_local_channels, num_channels):
     templates_ext.data["average"] = dense_templates
     operators = ["average"]
 
-    # Only expose the "std" operator when the file actually stores waveform_sd. Fabricating a zero std
-    # would falsely imply the unit's spikes have no waveform variability, so when it is absent we declare
-    # only "average" rather than inventing data.
-    if "waveform_sd" in units.columns:
-        waveform_sd = np.array([np.asarray(t, dtype="float") for t in units["waveform_sd"].values])
+    # Only expose the "std" operator when the file actually stores a std column (resolved). Fabricating a
+    # zero std would falsely imply the unit's spikes have no waveform variability, so when it is absent we
+    # declare only "average" rather than inventing data.
+    if templates_std_column is not None and templates_std_column in units.columns:
+        waveform_sd = np.array([np.asarray(t, dtype="float") for t in units[templates_std_column].values])
         dense_std = np.zeros((analyzer.get_num_units(), num_samples_template, num_channels), dtype="float32")
         for unit_index, positions in enumerate(unit_local_channels):
             k = len(positions)
@@ -2227,33 +2316,18 @@ def _make_templates(analyzer, units, unit_local_channels, num_channels):
     analyzer.extensions["templates"] = templates_ext
 
 
-def _make_metrics(analyzer, units, metric_colnames, verbose=False):
-    """Route every per-unit numeric column to its metrics extension. Canonical template-metric names go
-    to `template_metrics`; everything else goes to `quality_metrics`, including pipeline-specific metrics
-    that are not canonical SpikeInterface names (e.g. IBL's noise_cutoff), so they are preserved and
-    shown in the metrics view rather than dropped."""
+def _make_metrics(analyzer, units, quality_metric_colnames, verbose=False):
+    """Attach the `quality_metrics` extension from the resolved quality-metric columns. Other per-unit
+    scalar columns are loaded as sorting properties instead (see the caller), and `template_metrics` is
+    never synthesized from value columns because it carries structured data a values-only table lacks."""
     import pandas as pd
-    from spikeinterface.metrics.template import ComputeTemplateMetrics
     from spikeinterface.metrics.quality import ComputeQualityMetrics
 
-    template_metric_columns = ComputeTemplateMetrics.get_metric_columns()
-    template_metrics_df = pd.DataFrame(index=analyzer.unit_ids)
     quality_metric_df = pd.DataFrame(index=analyzer.unit_ids)
-    for col in metric_colnames:
-        if col not in units.columns:
-            continue
-        if col in template_metric_columns:
-            template_metrics_df.loc[:, col] = units[col].values
-        else:
+    for col in quality_metric_colnames:
+        if col in units.columns:
             quality_metric_df.loc[:, col] = units[col].values
 
-    if len(template_metrics_df.columns) > 0:
-        if verbose:
-            print("Adding template metrics")
-        template_metrics_ext = ComputeTemplateMetrics(analyzer)
-        template_metrics_ext.data["metrics"] = template_metrics_ext._cast_metrics(template_metrics_df)
-        template_metrics_ext.run_info["run_completed"] = True
-        analyzer.extensions["template_metrics"] = template_metrics_ext
     if len(quality_metric_df.columns) > 0:
         if verbose:
             print("Adding quality metrics")

From 3a42237c01fecd96076c1a4d423b57a63e1654ab Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 9 Jul 2026 16:19:07 +0100
Subject: [PATCH 18/46] tests

---
 .../extractors/tests/test_nwbextractors.py    | 106 ++++++++++++++++++
 1 file changed, 106 insertions(+)

diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index d3f26f84db..0f0754b8ce 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1266,6 +1266,112 @@ def test_time_series_recording_equality_with_pynwb_and_backend(generate_nwbfile_
     check_recordings_equal(recording_backend, recording_pynwb)
 
 
+def _make_units_nwb(path, n_units=6, n_ch=8, n_samp=30, with_std=False):
+    """Build a minimal NWB file with a Units table for read_nwb_sorting_analyzer tests.
+
+    Pure pynwb, no external writers. The Units table has `waveform_mean` (templates), a canonical
+    quality metric (`snr`), a non-canonical numeric column (`custom_score`), a string label
+    (`ks_label`), and a per-unit `electrodes` region; the electrodes table carries `rel_x`/`rel_y`.
+    """
+    from datetime import datetime
+    from pynwb import NWBFile, NWBHDF5IO
+
+    nwbfile = NWBFile(session_description="t", identifier="id", session_start_time=datetime(2020, 1, 1).astimezone())
+    device = nwbfile.create_device(name="probe")
+    group = nwbfile.create_electrode_group(name="s0", description="d", location="loc", device=device)
+    nwbfile.add_electrode_column(name="rel_x", description="x")
+    nwbfile.add_electrode_column(name="rel_y", description="y")
+    for i in range(n_ch):
+        nwbfile.add_electrode(group=group, location="brain", rel_x=float(i % 4) * 20.0, rel_y=float(i // 4) * 20.0)
+    nwbfile.add_unit_column(name="snr", description="canonical quality metric")
+    nwbfile.add_unit_column(name="custom_score", description="non-canonical per-unit value")
+    nwbfile.add_unit_column(name="ks_label", description="curation label")
+    rng = np.random.default_rng(0)
+    for u in range(n_units):
+        kwargs = dict(
+            spike_times=np.sort(rng.uniform(0, 100, size=int(rng.integers(50, 500)))),
+            electrodes=list(range(n_ch)),
+            waveform_mean=rng.standard_normal((n_samp, n_ch)),
+            snr=float(rng.uniform(2, 10)),
+            custom_score=float(rng.uniform(0, 1)),
+            ks_label="good" if u % 2 == 0 else "mua",
+        )
+        if with_std:
+            kwargs["waveform_sd"] = np.abs(rng.standard_normal((n_samp, n_ch)))
+        nwbfile.add_unit(**kwargs)
+    with NWBHDF5IO(str(path), "w") as io:
+        io.write(nwbfile)
+    return str(path)
+
+
+def test_read_nwb_sorting_analyzer_default(tmp_path):
+    pytest.importorskip("pynwb")
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    path = _make_units_nwb(tmp_path / "units.nwb")
+    analyzer = read_nwb_sorting_analyzer(path, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None)
+
+    # recordingless, sparse, with templates for every unit
+    assert not analyzer.has_recording()
+    assert analyzer.sparsity is not None
+    assert analyzer.has_extension("templates")
+    assert analyzer.get_extension("templates").get_data().shape[0] == analyzer.get_num_units()
+
+    # by default only canonical quality columns become quality_metrics; the rest become properties, and
+    # template_metrics is never synthesized from value columns
+    assert list(analyzer.get_extension("quality_metrics").get_data().columns) == ["snr"]
+    assert not analyzer.has_extension("template_metrics")
+    properties = set(analyzer.sorting.get_property_keys())
+    assert "custom_score" in properties and "ks_label" in properties
+    assert "snr" not in properties
+
+
+def test_read_nwb_sorting_analyzer_extension_map_override(tmp_path):
+    pytest.importorskip("pynwb")
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    path = _make_units_nwb(tmp_path / "units.nwb")
+    analyzer = read_nwb_sorting_analyzer(
+        path,
+        use_pynwb=False,
+        sampling_frequency=30000.0,
+        compute_extra=None,
+        extension_map={"quality_metrics": [{"source": "columns", "columns": ["snr", "custom_score"]}]},
+    )
+    assert set(analyzer.get_extension("quality_metrics").get_data().columns) == {"snr", "custom_score"}
+    # custom_score is now a metric, so it is no longer a property
+    assert "custom_score" not in analyzer.sorting.get_property_keys()
+
+
+def test_read_nwb_sorting_analyzer_extension_map_disable(tmp_path):
+    pytest.importorskip("pynwb")
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    path = _make_units_nwb(tmp_path / "units.nwb")
+    analyzer = read_nwb_sorting_analyzer(
+        path,
+        use_pynwb=False,
+        sampling_frequency=30000.0,
+        compute_extra=None,
+        extension_map={"quality_metrics": None},
+    )
+    assert not analyzer.has_extension("quality_metrics")
+    # with quality_metrics disabled, every scalar column (including snr) becomes a property
+    assert {"snr", "custom_score", "ks_label"} <= set(analyzer.sorting.get_property_keys())
+
+
+def test_read_nwb_sorting_analyzer_waveform_sd(tmp_path):
+    pytest.importorskip("pynwb")
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    path = _make_units_nwb(tmp_path / "units.nwb", with_std=True)
+    analyzer = read_nwb_sorting_analyzer(path, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None)
+    templates = analyzer.get_extension("templates")
+    # the std operator is populated only because the file stores waveform_sd
+    assert "std" in templates.params["operators"]
+    assert "std" in templates.data
+
+
 if __name__ == "__main__":
     tmp_path = Path("tmp")
     if tmp_path.is_dir():

From 8d020740465e5f396da4ab96d1cde797b80bd833 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 9 Jul 2026 17:22:37 +0100
Subject: [PATCH 19/46] review

---
 .../extractors/nwbextractors.py               | 30 +++++++++----------
 1 file changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 0eb5550d8f..45ba91afc5 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1431,25 +1431,25 @@ def _compute_and_cache_spike_vector(self) -> None:
         # Performance: this deliberately overrides the generic BaseSorting builder, and the override is
         # the point. The generic builder reads one spike train per unit (one streamed read per unit),
         # whereas an NWB Units table stores every unit's spikes in a single flat times dataset grouped by
-        # unit. Reading that dataset once (see NwbSortingSegment._get_all_spike_frames_and_unit_indices)
+        # unit. Reading that dataset once (see NwbSortingSegment._get_all_spike_samples_and_unit_indices)
         # turns num_units streamed reads into one, which is a large difference for remote/streamed files.
         # Do not drop this override in favor of the generic path without a reason: doing so silently
         # reintroduces the per-unit read pattern. An NWB Units table is always single-segment, so the
         # single segment is built directly.
         segment = self.segments[0]
-        sample_indices, unit_indices = segment._get_all_spike_frames_and_unit_indices()
+        sample_indices, unit_indices = segment._get_all_spike_samples_and_unit_indices()
 
         # Stable sort by sample_index keeps unit_index ordering on ties, matching the way the generic
         # BaseSorting builder constructs the vector (unit order within a sample).
         order = np.argsort(sample_indices, stable=True)
-        n = sample_indices.size
+        num_spikes = sample_indices.size
 
-        spikes = np.zeros(n, dtype=minimum_spike_dtype)
+        spikes = np.zeros(num_spikes, dtype=minimum_spike_dtype)
         spikes["sample_index"] = sample_indices[order]
         spikes["unit_index"] = unit_indices[order]
         # segment_index stays 0 for the single segment.
         self._cached_spike_vector = spikes
-        self._cached_spike_vector_segment_slices = np.array([[0, n]], dtype="int64")
+        self._cached_spike_vector_segment_slices = np.array([[0, num_spikes]], dtype="int64")
 
     @staticmethod
     def fetch_available_units_tables(
@@ -1493,17 +1493,17 @@ def _frame_to_time(self, frame):
             return self._time_vector[frame]
         return frame / self._sampling_frequency + self._t_start
 
-    def _times_to_frames(self, spike_times) -> np.ndarray:
+    def _times_to_samples(self, spike_times) -> np.ndarray:
         # Map seconds to samples. An irregular clock (per-sample timestamps) maps exactly with a
-        # searchsorted into the time vector; a uniform clock uses frame = (t - t_start) * sampling_frequency.
+        # searchsorted into the time vector; a uniform clock uses sample = (t - t_start) * sampling_frequency.
         if self._time_vector is not None:
             # searchsorted needs the timestamps in memory: materialize them once, on first use, and cache.
             if not isinstance(self._time_vector, np.ndarray):
                 self._time_vector = np.asarray(self._time_vector)
-            frames = np.searchsorted(self._time_vector, spike_times, side="right") - 1
+            samples = np.searchsorted(self._time_vector, spike_times, side="right") - 1
         else:
-            frames = np.round((spike_times - self._t_start) * self._sampling_frequency)
-        return frames.astype("int64", copy=False)
+            samples = np.round((spike_times - self._t_start) * self._sampling_frequency)
+        return samples.astype("int64", copy=False)
 
     def get_unit_spike_train(
         self,
@@ -1518,10 +1518,10 @@ def get_unit_spike_train(
         # Get spike times in seconds
         spike_times = self.get_unit_spike_train_in_seconds(unit_id=unit_id, start_time=start_time, end_time=end_time)
 
-        return self._times_to_frames(spike_times)
+        return self._times_to_samples(spike_times)
 
-    def _get_all_spike_frames_and_unit_indices(self):
-        """Read all units' spike trains (as frames) in one bulk read, plus each spike's unit index.
+    def _get_all_spike_samples_and_unit_indices(self):
+        """Read all units' spike trains (as samples) in one bulk read, plus each spike's unit index.
 
         NWB stores all spike times in a single flat dataset grouped by unit (with per-unit
         boundaries in ``spike_times_index_data``), so the whole segment is read at once instead
@@ -1532,7 +1532,7 @@ def _get_all_spike_frames_and_unit_indices(self):
         boundaries = np.concatenate(([0], np.asarray(self.spike_times_index_data[:], dtype="int64")))
         counts = np.diff(boundaries)
         unit_indices = np.repeat(np.arange(counts.size, dtype="int64"), counts)
-        sample_indices = self._times_to_frames(spike_times)
+        sample_indices = self._times_to_samples(spike_times)
         return sample_indices, unit_indices
 
     def get_unit_spike_train_in_seconds(
@@ -1544,7 +1544,7 @@ def get_unit_spike_train_in_seconds(
         """Get the spike train times for a unit in seconds.
 
         This method returns spike times directly in seconds without conversion
-        to frames, avoiding double conversion for NWB files that already store
+        to samples, avoiding double conversion for NWB files that already store
         spike times as timestamps.
 
         Parameters

From 0ed91530490fc3c83046817aaad4e00cb0d32da7 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 9 Jul 2026 17:26:55 +0100
Subject: [PATCH 20/46] add comments

---
 src/spikeinterface/extractors/nwbextractors.py | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 45ba91afc5..fd45aa07c2 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1439,8 +1439,12 @@ def _compute_and_cache_spike_vector(self) -> None:
         segment = self.segments[0]
         sample_indices, unit_indices = segment._get_all_spike_samples_and_unit_indices()
 
-        # Stable sort by sample_index keeps unit_index ordering on ties, matching the way the generic
-        # BaseSorting builder constructs the vector (unit order within a sample).
+        # NWB stores spike_times grouped by unit, not globally ordered by time: pynwb itself notes that
+        # "the earliest spike may be in any unit" (Units.get_earliest_spike_time), and within-unit ordering
+        # is only a best-practice (checked by nwbinspector), not a schema guarantee. So a global sort is
+        # required to build the time-ordered spike vector. It is stable so unit_index order is preserved on
+        # equal sample_index ties, matching the generic BaseSorting builder; on already-ordered input the
+        # sort is a no-op (argsort returns the identity permutation).
         order = np.argsort(sample_indices, stable=True)
         num_spikes = sample_indices.size
 
@@ -1528,7 +1532,14 @@ def _get_all_spike_samples_and_unit_indices(self):
         of one streamed slice per unit.
         """
         spike_times = np.asarray(self.spike_times_data[:], dtype="float64")
-        # boundaries[u]:boundaries[u + 1] is the slice of spike_times belonging to unit u.
+        # NWB stores a ragged column as one flat data array (spike_times, all units concatenated) plus an
+        # index array of cumulative end offsets (spike_times_index): spike_times_index[u] is where unit u's
+        # spikes end and unit u+1's begin, so unit u is spike_times[spike_times_index[u-1] : spike_times_index[u]].
+        # Prepending a 0 gives every unit an explicit start, so boundaries[u]:boundaries[u+1] is unit u's
+        # slice; the gap between consecutive boundaries (np.diff) is that unit's spike count, and repeating
+        # each unit index by its count labels every spike with its unit.
+        # Example: 3 units with 2, 0, 3 spikes -> spike_times_index = [2, 2, 5],
+        #   boundaries = [0, 2, 2, 5], counts = [2, 0, 3], unit_indices = [0, 0, 2, 2, 2].
         boundaries = np.concatenate(([0], np.asarray(self.spike_times_index_data[:], dtype="int64")))
         counts = np.diff(boundaries)
         unit_indices = np.repeat(np.arange(counts.size, dtype="int64"), counts)

From 8b64779b503c61401c002fb40fa5a2e15231b2b6 Mon Sep 17 00:00:00 2001
From: Alessio Buccino <alejoe9187@gmail.com>
Date: Thu, 16 Jul 2026 13:43:11 +0200
Subject: [PATCH 21/46] wip

---
 src/spikeinterface/core/basesorting.py        | 17 +++++++++++++++
 .../extractors/nwbextractors.py               | 21 ++++++++++++++++++-
 2 files changed, 37 insertions(+), 1 deletion(-)

diff --git a/src/spikeinterface/core/basesorting.py b/src/spikeinterface/core/basesorting.py
index 51ec90bf37..3c023bb6e7 100644
--- a/src/spikeinterface/core/basesorting.py
+++ b/src/spikeinterface/core/basesorting.py
@@ -471,6 +471,23 @@ def get_last_spike_frame(self, segment_index: int | None = None) -> int:
             return 0
         return int(np.max(spikes_in_segment["sample_index"]))
 
+    def get_last_spike_time(self, segment_index: int | None = None) -> float:
+        """Get the time of the last spike in a segment across all units.
+
+        Parameters
+        ----------
+        segment_index : int or None, default: None
+            The segment index (required for multi-segment)
+
+        Returns
+        -------
+        float
+            The time of the last spike in seconds, or 0.0 if no spikes exist.
+        """
+        segment_index = self._check_segment_index(segment_index)
+        last_spike_frame = self.get_last_spike_frame(segment_index=segment_index)
+        return self.sample_index_to_time(last_spike_frame, segment_index=segment_index)
+
     def get_times(
         self,
         segment_index: int | None = None,
diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 7920c64410..7d0d89e364 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1553,6 +1553,25 @@ def _get_all_spike_samples_and_unit_indices(self):
         sample_indices = self._times_to_samples(spike_times)
         return sample_indices, unit_indices
 
+    def get_last_spike_time(self, segment_index: int | None = None) -> float:
+        """Get the time of the last spike in a segment across all units.
+        Overridden from BaseSorting to use spike_times_data.
+
+        Parameters
+        ----------
+        segment_index : int or None, default: None
+            The segment index (required for multi-segment)
+
+        Returns
+        -------
+        float
+            The time of the last spike in seconds, or 0.0 if no spikes exist.
+        """
+        segment = self.segments[segment_index] if segment_index is not None else self.segments[0]
+        if segment.spike_times_data.size == 0:
+            return 0.0
+        return segment.spike_times_data[-1]
+
     def get_unit_spike_train_in_seconds(
         self,
         unit_id,
@@ -2197,7 +2216,7 @@ def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, elect
     # global maximum. spike_times is in seconds, so this is already a duration in seconds. It is only an
     # approximation: because spike_times is concatenated per unit and not globally sorted, this is the
     # last unit's last spike, not necessarily the latest spike overall.
-    last_spike_time = float(np.asarray(sorting._sorting_segments[0].spike_times_data[-1]))
+    last_spike_time = float(np.asarray(sorting.get_last_spike_time(segment_index=0)))
     duration = last_spike_time + 1.0
 
     probe = Probe(si_units="um")

From 31f3e6819f62b69bcf713ca89c524e761d659249 Mon Sep 17 00:00:00 2001
From: Alessio Buccino <alejoe9187@gmail.com>
Date: Thu, 16 Jul 2026 14:47:01 +0200
Subject: [PATCH 22/46] fix: always compute all random spikes

---
 .../extractors/nwbextractors.py               | 35 ++++++++-----------
 1 file changed, 15 insertions(+), 20 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 7d0d89e364..17af691c11 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2156,7 +2156,7 @@ def read_nwb_sorting_analyzer(
 
     # templates (from the resolved templates column), and the required random_spikes extension
     if templates_column is not None and templates_column in units and unit_local_channels is not None:
-        _make_random_spikes(analyzer, sorting, has_recording=recording is not None)
+        _make_random_spikes(analyzer, sorting)
         _make_templates(analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels)
 
     # the resolved quality-metric columns -> quality_metrics extension
@@ -2271,30 +2271,25 @@ def _make_sparsity_from_electrodes(sorting, electrodes_table, electrodes_indices
     return sparsity, unit_local_channels
 
 
-def _make_random_spikes(analyzer, sorting, has_recording):
-    """Attach the required `random_spikes` extension.
+def _make_random_spikes(analyzer, sorting):
+    """Attach the required `random_spikes` extension, always using the "all" method.
 
-    With a recording it drives waveform extraction, so compute it for real. Without a recording it is
-    never used to extract anything (there are no traces); computing it would call
-    `sorting.to_spike_vector()`, which materializes the whole spike_times array (e.g. ~165 MB for a
-    full IBL session) and dominates the cost. Instead build it from the per-unit spike counts (the
-    spike_times index), so spike_times is not read here; the lazy sorting still reads a unit's spikes on
-    demand when a view needs them.
+    If the segment is an NwbSortingSegment sorting segment, it uses the spike_times_index_data to determine the
+    total number of spikes. Otherwise, it counts the total number of spikes from the sorting object.
     """
-    if has_recording:
-        analyzer.compute("random_spikes", method="all")
-        return
-
     from spikeinterface.core.analyzer_extension_core import ComputeRandomSpikes
 
-    spike_times_index = sorting._sorting_segments[0].spike_times_index_data
-    total_num_spikes = int(np.asarray(spike_times_index[-1]))
-    max_spikes_per_unit = 500
-    num_random = min(max_spikes_per_unit * sorting.get_num_units(), total_num_spikes)
-    random_spikes_indices = np.sort(np.random.default_rng(0).choice(total_num_spikes, size=num_random, replace=False))
+    segment = sorting._sorting_segments[0]
+    if isinstance(segment, NwbSortingSegment):
+        # the lazy segment has a spike_times_index_data array that is the same as the NWB spike_times_index
+        # dataset, so we can use it to get the total number of spikes without materializing spike_times.
+        spike_times_index = segment.spike_times_index_data
+        total_num_spikes = int(np.asarray(spike_times_index[-1]))
+    else:
+        total_num_spikes = sorting.count_total_num_spikes()
     random_spikes_ext = ComputeRandomSpikes(sorting_analyzer=analyzer)
-    random_spikes_ext.set_params(method="uniform", max_spikes_per_unit=max_spikes_per_unit, seed=0)
-    random_spikes_ext.data["random_spikes_indices"] = random_spikes_indices
+    random_spikes_ext.set_params(method="all")
+    random_spikes_ext.data["random_spikes_indices"] = np.arange(total_num_spikes, dtype="int64")
     random_spikes_ext.run_info["run_completed"] = True
     analyzer.extensions["random_spikes"] = random_spikes_ext
 

From 3f1c13109d971af33d8bf3bc144b3595e59a64f2 Mon Sep 17 00:00:00 2001
From: Alessio Buccino <alejoe9187@gmail.com>
Date: Thu, 16 Jul 2026 16:12:59 +0200
Subject: [PATCH 23/46] fix for multi-group

---
 .../extractors/nwbextractors.py               | 25 +++++++++++++++++--
 1 file changed, 23 insertions(+), 2 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 17af691c11..3b5001031d 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2047,6 +2047,9 @@ def read_nwb_sorting_analyzer(
             and units[c].dtype.kind == "O"
             and not isinstance(units[c].iloc[0], (list, np.ndarray))
         ]
+        waveform_unit = units_table.waveform_unit
+        if waveform_unit is None:
+            waveform_unit = "volts"  # default to volts if not present
     else:
         units_group = sorting.units_table
         colnames = list(units_group.keys())
@@ -2062,6 +2065,12 @@ def read_nwb_sorting_analyzer(
                 metric_colnames.append(column_name)
             elif dataset.ndim == 1 and dataset.dtype.kind in "OSU":
                 label_colnames.append(column_name)
+        # In NWB the waveform unit is stored as the "unit" attribute of the waveform_mean dataset,
+        # not as a group-level attribute (this mirrors pynwb's Units.waveform_unit).
+        if "waveform_mean" in units_group:
+            waveform_unit = units_group["waveform_mean"].attrs.get("unit", "volts")
+        else:
+            waveform_unit = "volts"  # default to volts if not present
 
     # Resolve the extension map against the Units columns: which column holds the templates, and which
     # scalar columns are quality metrics (canonical names by default) vs plain sorting properties.
@@ -2157,7 +2166,9 @@ def read_nwb_sorting_analyzer(
     # templates (from the resolved templates column), and the required random_spikes extension
     if templates_column is not None and templates_column in units and unit_local_channels is not None:
         _make_random_spikes(analyzer, sorting)
-        _make_templates(analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels)
+        _make_templates(
+            analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels, waveform_unit
+        )
 
     # the resolved quality-metric columns -> quality_metrics extension
     _make_metrics(analyzer, units, quality_metric_colnames, verbose=verbose)
@@ -2294,7 +2305,9 @@ def _make_random_spikes(analyzer, sorting):
     analyzer.extensions["random_spikes"] = random_spikes_ext
 
 
-def _make_templates(analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels):
+def _make_templates(
+    analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels, waveform_unit
+):
     """Attach the `templates` extension from the resolved templates column, scattering each unit's
     sparse waveform block onto its true channel positions."""
     from spikeinterface.core.analyzer_extension_core import ComputeTemplates
@@ -2312,6 +2325,14 @@ def _make_templates(analyzer, units, templates_column, templates_std_column, uni
         k = len(positions)
         dense_templates[unit_index][:, positions] = waveform_mean[unit_index][:, :k]
 
+    # TODO: this would be nice, but there is a bug in NWB. Even when explicitly settting "microvolts" as the waveform
+    # unit, pynwb writes "volts" into the waveform_mean dataset's unit attribute. So we can't rely on it to be correct.
+    # For now, we just assume the waveform is in volts and convert to microvolts.
+    # if waveform_unit == "volts":
+    #     dense_templates *= 1e6  # convert to microvolts
+    # elif waveform_unit == "millivolts":
+    #     dense_templates *= 1e3  # convert to microvolts
+
     # scalar alignment (nbefore): the sample where the average peak-channel amplitude is largest
     peak_amplitude_per_sample = np.abs(np.nan_to_num(waveform_mean)).max(axis=2).mean(axis=0)
     nbefore = int(np.argmax(peak_amplitude_per_sample))

From ca815baa73c35273eb6cb7f234fe72ff745d6f6d Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Fri, 17 Jul 2026 10:48:28 -0600
Subject: [PATCH 24/46] fix groups

---
 .../extractors/nwbextractors.py               | 21 ++++++++++---------
 1 file changed, 11 insertions(+), 10 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 3b5001031d..27ccfa3181 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2206,12 +2206,10 @@ def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, elect
     electrodes_table_sliced = electrodes_table.iloc[electrode_indices_all]
 
     if "channel_name" in electrodes_table_sliced:
-        # IBL stores channel_name as variable-length HDF5 byte strings (e.g. b"AP0"); decode to str so
-        # SpikeInterface accepts them as channel ids (object-dtype arrays are rejected). This must match
-        # the decode in _make_sparsity_from_electrodes so channel ids and sparsity keys stay consistent.
-        channel_ids = np.array(
-            [c.decode("utf-8") if isinstance(c, bytes) else c for c in electrodes_table_sliced["channel_name"][:]]
-        )
+        # channel_name is already decoded to str at the source (_create_df_from_nwb_table), but pandas
+        # holds string columns as object dtype; cast to a numpy unicode array so SpikeInterface accepts
+        # the channel ids (object-dtype ids are rejected).
+        channel_ids = np.asarray(electrodes_table_sliced["channel_name"][:], dtype=str)
     else:
         channel_ids = electrodes_table_sliced.index.to_numpy()
 
@@ -2256,10 +2254,7 @@ def _make_sparsity_from_electrodes(sorting, electrodes_table, electrodes_indices
 
     num_channels = len(analyzer_channel_ids)
     if "channel_name" in electrodes_table.columns:
-        # decode HDF5 byte-string channel names to str, matching _make_placeholder_recording_from_electrodes
-        electrode_row_to_channel_id = np.array(
-            [c.decode("utf-8") if isinstance(c, bytes) else c for c in electrodes_table["channel_name"].to_numpy()]
-        )
+        electrode_row_to_channel_id = electrodes_table["channel_name"].to_numpy()
     else:
         electrode_row_to_channel_id = electrodes_table.index.to_numpy()
     position_of_channel_id = {channel_id: pos for pos, channel_id in enumerate(analyzer_channel_ids)}
@@ -2408,6 +2403,12 @@ def _create_df_from_nwb_table(group, columns=None):
         elif item.ndim > 1:
             data[col] = [item_flat for item_flat in item]
         else:
+            if item.dtype.kind in ("O", "S"):
+                # HDF5 stores string columns as variable-length or fixed-length bytes (e.g. b"AP0",
+                # b"Probe00"). Decode to str at the source so every caller gets plain strings: channel and
+                # group ids must be str (SpikeInterface rejects object-dtype byte-string ids), and group
+                # selection compares the column against a str argument.
+                item = np.array([v.decode("utf-8") if isinstance(v, bytes) else v for v in item])
             data[col] = item
     df = pd.DataFrame(data=data)
     df.set_index("id", inplace=True)

From 61ffd02cf6bdbc012d78561d60e63c1a3ead257c Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Fri, 17 Jul 2026 13:14:02 -0600
Subject: [PATCH 25/46] logic for unit fixing

---
 .../extractors/nwbextractors.py               | 64 +++++++++++++------
 1 file changed, 46 insertions(+), 18 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 27ccfa3181..3f0a5b15f8 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1985,6 +1985,7 @@ def read_nwb_sorting_analyzer(
     compute_extra: List[str] | None = ["unit_locations"],
     compute_extra_params: dict | None = None,
     extension_map: dict | None = None,
+    rescale_templates_to_uV: bool = True,
     verbose: bool = False,
 ) -> SortingAnalyzer:
     # extension_map overrides (per extension) merge over DEFAULT_EXTENSION_MAP; see its docstring.
@@ -2047,9 +2048,6 @@ def read_nwb_sorting_analyzer(
             and units[c].dtype.kind == "O"
             and not isinstance(units[c].iloc[0], (list, np.ndarray))
         ]
-        waveform_unit = units_table.waveform_unit
-        if waveform_unit is None:
-            waveform_unit = "volts"  # default to volts if not present
     else:
         units_group = sorting.units_table
         colnames = list(units_group.keys())
@@ -2065,12 +2063,6 @@ def read_nwb_sorting_analyzer(
                 metric_colnames.append(column_name)
             elif dataset.ndim == 1 and dataset.dtype.kind in "OSU":
                 label_colnames.append(column_name)
-        # In NWB the waveform unit is stored as the "unit" attribute of the waveform_mean dataset,
-        # not as a group-level attribute (this mirrors pynwb's Units.waveform_unit).
-        if "waveform_mean" in units_group:
-            waveform_unit = units_group["waveform_mean"].attrs.get("unit", "volts")
-        else:
-            waveform_unit = "volts"  # default to volts if not present
 
     # Resolve the extension map against the Units columns: which column holds the templates, and which
     # scalar columns are quality metrics (canonical names by default) vs plain sorting properties.
@@ -2167,7 +2159,13 @@ def read_nwb_sorting_analyzer(
     if templates_column is not None and templates_column in units and unit_local_channels is not None:
         _make_random_spikes(analyzer, sorting)
         _make_templates(
-            analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels, waveform_unit
+            analyzer,
+            units,
+            templates_column,
+            templates_std_column,
+            unit_local_channels,
+            num_channels,
+            rescale_to_uV=rescale_templates_to_uV,
         )
 
     # the resolved quality-metric columns -> quality_metrics extension
@@ -2300,8 +2298,34 @@ def _make_random_spikes(analyzer, sorting):
     analyzer.extensions["random_spikes"] = random_spikes_ext
 
 
+def _infer_template_scale_to_uV(dense_templates):
+    """Infer the factor that brings templates into microvolts, SpikeInterface's convention.
+
+    The NWB waveform unit attribute is schema-fixed to "volts" and never reflects the true scale
+    (pynwb issue #2162), so it cannot drive the conversion. The template magnitude can: an extracellular
+    spike peaks at roughly tens to hundreds of microvolts, so pick the factor in {1, 1e3, 1e6} that lands
+    the median per-unit peak in that physiological window. This reads no extra data (it operates on the
+    already-built templates). Returns 1.0 if the magnitude is unusable (all zero).
+    """
+    per_unit_peak = np.abs(dense_templates).max(axis=(1, 2))
+    per_unit_peak = per_unit_peak[per_unit_peak > 0]
+    if per_unit_peak.size == 0:
+        return 1.0
+    median_peak = float(np.median(per_unit_peak))
+    for factor in (1.0, 1e3, 1e6):
+        if 5.0 <= median_peak * factor <= 2000.0:
+            return factor
+    return 1.0
+
+
 def _make_templates(
-    analyzer, units, templates_column, templates_std_column, unit_local_channels, num_channels, waveform_unit
+    analyzer,
+    units,
+    templates_column,
+    templates_std_column,
+    unit_local_channels,
+    num_channels,
+    rescale_to_uV=False,
 ):
     """Attach the `templates` extension from the resolved templates column, scattering each unit's
     sparse waveform block onto its true channel positions."""
@@ -2320,13 +2344,15 @@ def _make_templates(
         k = len(positions)
         dense_templates[unit_index][:, positions] = waveform_mean[unit_index][:, :k]
 
-    # TODO: this would be nice, but there is a bug in NWB. Even when explicitly settting "microvolts" as the waveform
-    # unit, pynwb writes "volts" into the waveform_mean dataset's unit attribute. So we can't rely on it to be correct.
-    # For now, we just assume the waveform is in volts and convert to microvolts.
-    # if waveform_unit == "volts":
-    #     dense_templates *= 1e6  # convert to microvolts
-    # elif waveform_unit == "millivolts":
-    #     dense_templates *= 1e3  # convert to microvolts
+    # The stored waveform unit attribute is schema-fixed to "volts" and never reflects the real scale
+    # (pynwb #2162), so it cannot drive a conversion. When rescaling is requested, infer the factor from
+    # the template magnitude instead and apply the same factor to the std templates below so both stay in
+    # the same unit. On by default (read_nwb_sorting_analyzer) so templates read in microvolts; it can be
+    # turned off to keep the file's native unit, since the scale does not affect peak channel, shape,
+    # location, or sparsity.
+    template_scale = _infer_template_scale_to_uV(dense_templates) if rescale_to_uV else 1.0
+    if template_scale != 1.0:
+        dense_templates *= template_scale
 
     # scalar alignment (nbefore): the sample where the average peak-channel amplitude is largest
     peak_amplitude_per_sample = np.abs(np.nan_to_num(waveform_mean)).max(axis=2).mean(axis=0)
@@ -2345,6 +2371,8 @@ def _make_templates(
         for unit_index, positions in enumerate(unit_local_channels):
             k = len(positions)
             dense_std[unit_index][:, positions] = waveform_sd[unit_index][:, :k]
+        if template_scale != 1.0:
+            dense_std *= template_scale
         templates_ext.data["std"] = dense_std
         operators.append("std")
 

From 41ce6237dd2b7a206dff6c8670059067e7cd82c6 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Sat, 18 Jul 2026 00:38:05 -0600
Subject: [PATCH 26/46] improvements

---
 .../extractors/nwbextractors.py               | 96 ++++++++++---------
 1 file changed, 51 insertions(+), 45 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 3f0a5b15f8..1289cd0cb1 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -563,14 +563,23 @@ def _locate_units_table(self, backend):
         # Resolve unit_table_path (auto-discovering it when the file has exactly one Units table)
         # and return the table handle, raising a helpful error that lists the options on failure.
         if self.unit_table_path is None:
-            available = _find_neurodata_type_from_backend(self.file, neurodata_type="Units", backend=backend)
-            if len(available) != 1:
-                raise ValueError(
-                    "Multiple Units tables found in the file. "
-                    "Please specify the 'unit_table_path' argument:"
-                    f"Available options are: {available}."
-                )
-            self.unit_table_path = available[0]
+            # Fast path: the primary NWB Units table lives at the canonical top-level "units". Check it
+            # directly instead of walking the whole HDF5/zarr tree, which over a stream costs one network
+            # round-trip per object (hundreds on a rich file, dominating the extractor init). Only fall
+            # back to a full search when "units" is absent (a file that stores units in a non-standard
+            # location). When "units" is present it is used as-is, so a file with several Units tables is no
+            # longer flagged here; pass `unit_table_path` explicitly to select a non-canonical one.
+            if "units" in self.file and self.file["units"].attrs.get("neurodata_type") == "Units":
+                self.unit_table_path = "units"
+            else:
+                available = _find_neurodata_type_from_backend(self.file, neurodata_type="Units", backend=backend)
+                if len(available) != 1:
+                    raise ValueError(
+                        "Multiple Units tables found in the file. "
+                        "Please specify the 'unit_table_path' argument:"
+                        f"Available options are: {available}."
+                    )
+                self.unit_table_path = available[0]
         try:
             return self.file[self.unit_table_path]
         except KeyError:
@@ -1333,8 +1342,6 @@ def __init__(
         # unless it is named, because the ElectricalSeries in a file is not necessarily the one the
         # sorting was computed from.
 
-        # First, look for any ElectricalSeries in the file.
-        file_has_electrical_series = self._reader._has_electrical_series()
         electrical_series_path_is_provided = electrical_series_path is not None
         sampling_frequency_is_provided = self.provided_sampling_frequency is not None
         t_start_is_provided = self.t_start is not None
@@ -1357,18 +1364,13 @@ def __init__(
                 samples_for_rate_estimation
             )
         else:
-            # No ElectricalSeries named. An unnamed recording is ambiguous (it may not be the one the
-            # sorting was computed from), so t_start must be given explicitly. With no recording in the
-            # file, anchor the frame grid at 0 (the spike times in seconds are preserved regardless).
-            unnamed_recording_needs_t_start = file_has_electrical_series and not t_start_is_provided
-            if unnamed_recording_needs_t_start:
-                raise ValueError(
-                    "The file contains an ElectricalSeries but 't_start' was not provided. Pass "
-                    "'electrical_series_path' to read the time base from it, or pass 't_start' explicitly."
-                )
-
-            no_recording_to_align_to = not file_has_electrical_series and not t_start_is_provided
-            if no_recording_to_align_to:
+            # No ElectricalSeries named. Spike times are stored as seconds from session start, so anchor
+            # the frame grid at t_start = 0 by default (frames = times * sampling_frequency). Alignment to a
+            # particular recording is an explicit opt-in via `electrical_series_path`; we do not search the
+            # file for an unnamed ElectricalSeries to align to, since that walks the whole HDF5 tree (one
+            # network round-trip per object, dominating the streamed init) and the detected series would not
+            # be used for alignment anyway.
+            if not t_start_is_provided:
                 self.t_start = 0.0
 
             # sampling_frequency: the argument, or the Units.resolution attribute.
@@ -2125,8 +2127,19 @@ def read_nwb_sorting_analyzer(
         analyzer_recording = recording
         analyzer_channel_ids = list(recording.get_channel_ids())
     else:
+        # The placeholder recording only needs to loosely bound the timeline: it is dropped right after the
+        # build and its length feeds the GUI time axis, not any curation quantity. Reading the full
+        # spike_times array for the exact last spike would cost ~130 MB over a stream, so instead read only
+        # the last stored spike (one chunk, ~4 MB) and double it. spike_times is stored per-unit-
+        # concatenated (not globally time-sorted), so the last stored value is a lower bound on the true
+        # last spike; doubling generously covers the gap while over-estimating only adds a harmless empty
+        # tail. Read from sorting_tmp because a group selection wraps `sorting` in a UnitsSelectionSorting
+        # whose segment has no spike_times_data.
+        spike_times_data = sorting_tmp._sorting_segments[0].spike_times_data
+        last_stored_spike_time = 0.0 if spike_times_data.shape[0] == 0 else float(np.asarray(spike_times_data[-1]))
+        placeholder_duration = max(last_stored_spike_time * 2.0, 1.0)
         analyzer_recording, analyzer_channel_ids = _make_placeholder_recording_from_electrodes(
-            sorting, electrodes_table, electrodes_indices, verbose=verbose
+            sorting, electrodes_table, electrodes_indices, duration=placeholder_duration, verbose=verbose
         )
 
     # Per-unit sparsity and channel map from the Units `electrodes` region. Each unit's `waveform_mean`
@@ -2157,7 +2170,15 @@ def read_nwb_sorting_analyzer(
 
     # templates (from the resolved templates column), and the required random_spikes extension
     if templates_column is not None and templates_column in units and unit_local_channels is not None:
-        _make_random_spikes(analyzer, sorting)
+        # random_spikes only needs the total spike count of the (possibly group-selected) sorting. Read it
+        # from the NWB spike_times_index (per-unit cumulative counts, ~KB) on the unwrapped sorting_tmp,
+        # restricted to the selected units, rather than counting through the sorting object: a group
+        # selection wraps `sorting` in a UnitsSelectionSorting whose count materializes the parent's whole
+        # spike vector (~130-480 MB). np.isin keeps the units still present after any group selection.
+        per_unit_spike_counts = np.diff(sorting_tmp._sorting_segments[0].spike_times_index_data[:], prepend=0)
+        selected_units_mask = np.isin(sorting_tmp.unit_ids, sorting.unit_ids)
+        total_num_spikes = int(per_unit_spike_counts[selected_units_mask].sum())
+        _make_random_spikes(analyzer, total_num_spikes)
         _make_templates(
             analyzer,
             units,
@@ -2181,7 +2202,7 @@ def read_nwb_sorting_analyzer(
     return analyzer
 
 
-def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, electrodes_indices, verbose=False):
+def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, electrodes_indices, duration, verbose=False):
     """Build a lightweight placeholder recording that carries only the probe geometry of the electrodes
     a unit's waveforms live on, for the recordingless (streaming) case.
 
@@ -2217,15 +2238,6 @@ def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, elect
     ), "'rel_x' and 'rel_y' should be columns in the electrodes table"
     locations = np.array([electrodes_table_sliced["rel_x"][:], electrodes_table_sliced["rel_y"][:]]).T
 
-    # The recording length only needs to loosely bound the timeline for the recordingless GUI; nothing
-    # about curation depends on it being exact. Estimate it cheaply from the last stored spike time (one
-    # element, i.e. only the last spike_times chunk) rather than scanning the whole array for the true
-    # global maximum. spike_times is in seconds, so this is already a duration in seconds. It is only an
-    # approximation: because spike_times is concatenated per unit and not globally sorted, this is the
-    # last unit's last spike, not necessarily the latest spike overall.
-    last_spike_time = float(np.asarray(sorting.get_last_spike_time(segment_index=0)))
-    duration = last_spike_time + 1.0
-
     probe = Probe(si_units="um")
     probe.set_contacts(locations, shapes="circle", shape_params={"radius": 1})
     probe.set_device_channel_indices(np.arange(len(channel_ids)))
@@ -2275,22 +2287,16 @@ def _make_sparsity_from_electrodes(sorting, electrodes_table, electrodes_indices
     return sparsity, unit_local_channels
 
 
-def _make_random_spikes(analyzer, sorting):
+def _make_random_spikes(analyzer, total_num_spikes):
     """Attach the required `random_spikes` extension, always using the "all" method.
 
-    If the segment is an NwbSortingSegment sorting segment, it uses the spike_times_index_data to determine the
-    total number of spikes. Otherwise, it counts the total number of spikes from the sorting object.
+    With no recording this extension is never used to extract waveforms (there are no traces), but it is a
+    required dependency of templates, so it must exist and be consistent with the sorting's spike count.
+    `method="all"` means the indices are just `arange(total_num_spikes)`; the caller supplies the count
+    (read cheaply from `spike_times_index`) so this never materializes the spike vector.
     """
     from spikeinterface.core.analyzer_extension_core import ComputeRandomSpikes
 
-    segment = sorting._sorting_segments[0]
-    if isinstance(segment, NwbSortingSegment):
-        # the lazy segment has a spike_times_index_data array that is the same as the NWB spike_times_index
-        # dataset, so we can use it to get the total number of spikes without materializing spike_times.
-        spike_times_index = segment.spike_times_index_data
-        total_num_spikes = int(np.asarray(spike_times_index[-1]))
-    else:
-        total_num_spikes = sorting.count_total_num_spikes()
     random_spikes_ext = ComputeRandomSpikes(sorting_analyzer=analyzer)
     random_spikes_ext.set_params(method="all")
     random_spikes_ext.data["random_spikes_indices"] = np.arange(total_num_spikes, dtype="int64")

From 3a1d69e86c9d0b42db441612e037a868883b41be Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Sat, 18 Jul 2026 01:10:46 -0600
Subject: [PATCH 27/46] more fixes

---
 .../extractors/nwbextractors.py               | 21 ++++++++++++++++++-
 1 file changed, 20 insertions(+), 1 deletion(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 1289cd0cb1..bb1fbb255f 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -350,6 +350,25 @@ def _find_neurodata_type_from_backend(group, path="", result=None, neurodata_typ
     return result
 
 
+def _find_electrical_series_paths(file_handle, backend="hdf5"):
+    """Paths of every ElectricalSeries, searching only where they can live instead of the whole file.
+
+    In NWB, ElectricalSeries live under ``/acquisition`` (raw) or ``/processing/ecephys`` (LFP/filtered,
+    possibly nested in an LFP or FilteredEphys container). Restricting the recursive search to those two
+    subtrees avoids walking hundreds of unrelated objects (trials, behavior, stimulus, ...) over a stream,
+    where each object visited is a network round-trip. A file that stores an ElectricalSeries in a
+    non-standard location must be read by passing ``electrical_series_path`` explicitly.
+    """
+    result = []
+    for root in ("acquisition", "processing/ecephys"):
+        try:
+            group = file_handle[root]
+        except (KeyError, TypeError):
+            continue
+        _find_neurodata_type_from_backend(group, root, result, "ElectricalSeries", backend)
+    return result
+
+
 def _retrieve_electrodes_indices_from_electrical_series_backend(open_file, electrical_series, backend="hdf5"):
     """
     Retrieves the indices of the electrodes from the electrical series.
@@ -751,7 +770,7 @@ def _locate_electrical_series(self, backend):
         # Resolve electrical_series_path (auto-discovering it when the file has exactly one series)
         # and return the series handle, raising a helpful error that lists the options on failure.
         if self.electrical_series_path is None:
-            available = _find_neurodata_type_from_backend(self.file, neurodata_type="ElectricalSeries", backend=backend)
+            available = _find_electrical_series_paths(self.file, backend=backend)
             if len(available) != 1:
                 raise ValueError(
                     "Multiple ElectricalSeries found in the file. "

From 37306db6e5a416e605d8eed45f8fd918b26bbf54 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Sat, 18 Jul 2026 11:26:32 -0600
Subject: [PATCH 28/46] remove dead code

---
 src/spikeinterface/extractors/nwbextractors.py | 15 ---------------
 1 file changed, 15 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index bb1fbb255f..88e5e3fcc3 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -657,21 +657,6 @@ def _sampling_frequency_from_units_metadata(self):
             return None
         return 1.0 / resolution
 
-    def _has_electrical_series(self):
-        # Whether the file contains any ElectricalSeries at all. Used to decide, for a sorting, whether
-        # there is a recording to align to: if there is none, t_start has no origin and defaults to 0.
-        if self.reading_method == "use_pynwb":
-            from pynwb.ecephys import ElectricalSeries
-
-            electrical_series = [item for item in self.nwbfile.all_children() if isinstance(item, ElectricalSeries)]
-        else:
-            backend = "zarr" if self.reading_method == "use_zarr" else "hdf5"
-            electrical_series = _find_neurodata_type_from_backend(
-                self.file, neurodata_type="ElectricalSeries", backend=backend
-            )
-        file_has_electrical_series = len(electrical_series) > 0
-        return file_has_electrical_series
-
     def _fetch_unit_properties(self):
         # Return the unit-table columns as a name -> per-unit values mapping in a backend-independent
         # format. Skips id / spike-time columns, index columns, nested ragged arrays, and ragged

From 180a9f022d7ebff07b264c789b91b807771a3146 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Sat, 18 Jul 2026 18:12:39 -0600
Subject: [PATCH 29/46] fix tests

---
 .../extractors/nwbextractors.py               | 29 +++++++++++++++----
 1 file changed, 23 insertions(+), 6 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 88e5e3fcc3..ee89834239 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -657,6 +657,18 @@ def _sampling_frequency_from_units_metadata(self):
             return None
         return 1.0 / resolution
 
+    def _has_electrical_series(self):
+        # Whether the file contains any ElectricalSeries at all. Used to decide, for a sorting, whether
+        # there is a recording to align to: if there is none, t_start has no origin and defaults to 0.
+        if self.reading_method == "use_pynwb":
+            from pynwb.ecephys import ElectricalSeries
+
+            electrical_series = [item for item in self.nwbfile.all_children() if isinstance(item, ElectricalSeries)]
+        else:
+            backend = "zarr" if self.reading_method == "use_zarr" else "hdf5"
+            electrical_series = _find_electrical_series_paths(self.file, backend=backend)
+        return len(electrical_series) > 0
+
     def _fetch_unit_properties(self):
         # Return the unit-table columns as a name -> per-unit values mapping in a backend-independent
         # format. Skips id / spike-time columns, index columns, nested ragged arrays, and ragged
@@ -1368,13 +1380,18 @@ def __init__(
                 samples_for_rate_estimation
             )
         else:
-            # No ElectricalSeries named. Spike times are stored as seconds from session start, so anchor
-            # the frame grid at t_start = 0 by default (frames = times * sampling_frequency). Alignment to a
-            # particular recording is an explicit opt-in via `electrical_series_path`; we do not search the
-            # file for an unnamed ElectricalSeries to align to, since that walks the whole HDF5 tree (one
-            # network round-trip per object, dominating the streamed init) and the detected series would not
-            # be used for alignment anyway.
+            # No ElectricalSeries named. An unnamed recording is ambiguous (it may not be the one the
+            # sorting was computed from), so t_start must be given explicitly. With no recording in the
+            # file, anchor the frame grid at 0 (the spike times in seconds are preserved regardless).
+            # Whether the file has an ElectricalSeries only affects t_start's origin, so when t_start is
+            # already provided we skip that search entirely: it walks the whole HDF5 tree (one network
+            # round-trip per object) and otherwise dominates the extractor init on a rich streamed file.
             if not t_start_is_provided:
+                if self._reader._has_electrical_series():
+                    raise ValueError(
+                        "The file contains an ElectricalSeries but 't_start' was not provided. Pass "
+                        "'electrical_series_path' to read the time base from it, or pass 't_start' explicitly."
+                    )
                 self.t_start = 0.0
 
             # sampling_frequency: the argument, or the Units.resolution attribute.

From 7776d3a9886f88d8ea7421e0182bec324f9328b9 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Tue, 6 Oct 2026 13:15:11 -0600
Subject: [PATCH 30/46] Do not read per-spike Units columns on the pynwb path
 of read_nwb_sorting_analyzer

---
 .../extractors/nwbextractors.py               |  7 ++++-
 .../extractors/tests/test_nwbextractors.py    | 28 ++++++++++++++++++-
 2 files changed, 33 insertions(+), 2 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 1b32c45046..eae4b387c4 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2077,9 +2077,14 @@ def read_nwb_sorting_analyzer(
     # (identity properties). The large per-spike columns (spike amplitudes, spike depths, each as long
     # as the whole spike vector) are never read here, which is what makes streaming this table viable.
     if use_pynwb:
+        from hdmf.common import VectorIndex
+
         units_table = sorting.units_table
         colnames = list(units_table.colnames)
-        units = units_table.to_dataframe(index=True)
+        # ragged columns (per-spike data) are excluded so to_dataframe never reads them; electrodes is
+        # ragged too but is the per-unit region needed for sparsity
+        ragged_colnames = {c for c in colnames if c != "electrodes" and isinstance(units_table[c], VectorIndex)}
+        units = units_table.to_dataframe(index=True, exclude=ragged_colnames)
         structural = {"waveform_mean", "waveform_sd", "electrodes"}
         metric_colnames = [c for c in units.columns if c not in structural and units[c].dtype.kind in "fiu"]
         label_colnames = [
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index 0f0754b8ce..dda3cfbbef 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1266,7 +1266,7 @@ def test_time_series_recording_equality_with_pynwb_and_backend(generate_nwbfile_
     check_recordings_equal(recording_backend, recording_pynwb)
 
 
-def _make_units_nwb(path, n_units=6, n_ch=8, n_samp=30, with_std=False):
+def _make_units_nwb(path, n_units=6, n_ch=8, n_samp=30, with_std=False, with_per_spike=False):
     """Build a minimal NWB file with a Units table for read_nwb_sorting_analyzer tests.
 
     Pure pynwb, no external writers. The Units table has `waveform_mean` (templates), a canonical
@@ -1286,6 +1286,8 @@ def _make_units_nwb(path, n_units=6, n_ch=8, n_samp=30, with_std=False):
     nwbfile.add_unit_column(name="snr", description="canonical quality metric")
     nwbfile.add_unit_column(name="custom_score", description="non-canonical per-unit value")
     nwbfile.add_unit_column(name="ks_label", description="curation label")
+    if with_per_spike:
+        nwbfile.add_unit_column(name="spike_amplitudes_uV", description="per-spike amplitudes", index=True)
     rng = np.random.default_rng(0)
     for u in range(n_units):
         kwargs = dict(
@@ -1298,6 +1300,8 @@ def _make_units_nwb(path, n_units=6, n_ch=8, n_samp=30, with_std=False):
         )
         if with_std:
             kwargs["waveform_sd"] = np.abs(rng.standard_normal((n_samp, n_ch)))
+        if with_per_spike:
+            kwargs["spike_amplitudes_uV"] = rng.standard_normal(kwargs["spike_times"].size)
         nwbfile.add_unit(**kwargs)
     with NWBHDF5IO(str(path), "w") as io:
         io.write(nwbfile)
@@ -1372,6 +1376,28 @@ def test_read_nwb_sorting_analyzer_waveform_sd(tmp_path):
     assert "std" in templates.data
 
 
+@pytest.mark.parametrize("use_pynwb", [True, False])
+def test_read_nwb_sorting_analyzer_skips_per_spike_columns(tmp_path, use_pynwb, monkeypatch):
+    pytest.importorskip("pynwb")
+    import h5py
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    path = _make_units_nwb(tmp_path / "units.nwb", with_per_spike=True)
+
+    read_datasets = set()
+    original_getitem = h5py.Dataset.__getitem__
+
+    def logging_getitem(self, key):
+        read_datasets.add(self.name)
+        return original_getitem(self, key)
+
+    monkeypatch.setattr(h5py.Dataset, "__getitem__", logging_getitem)
+    analyzer = read_nwb_sorting_analyzer(path, use_pynwb=use_pynwb, sampling_frequency=30000.0, compute_extra=None)
+
+    assert "/units/spike_amplitudes_uV" not in read_datasets
+    assert "spike_amplitudes_uV" not in analyzer.sorting.get_property_keys()
+
+
 if __name__ == "__main__":
     tmp_path = Path("tmp")
     if tmp_path.is_dir():

From f04e2b2d908b9853cec814e2ea8b3da28600a1a4 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Tue, 6 Oct 2026 14:03:44 -0600
Subject: [PATCH 31/46] Require electrical_series_path in
 read_nwb_sorting_analyzer and warn when the duration is estimated

---
 .../extractors/nwbextractors.py               | 34 +++++++++++++------
 .../extractors/tests/test_nwbextractors.py    | 31 +++++++++++++++--
 2 files changed, 51 insertions(+), 14 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index eae4b387c4..3627fcf7e6 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2013,9 +2013,10 @@ def _resolve_extension_columns(extension_map, scalar_colnames, all_colnames):
 
 def read_nwb_sorting_analyzer(
     file_path: str | Path,
+    electrical_series_path: str | None,
     t_start: float | None = None,
     sampling_frequency: float | None = None,
-    electrical_series_path: str | None = None,
+    recording_duration: float | None = None,
     unit_table_path: str | None = None,
     stream_mode: Literal["fsspec", "remfile", "zarr"] | None = None,
     stream_cache_path: str | Path | None = None,
@@ -2033,8 +2034,10 @@ def read_nwb_sorting_analyzer(
     resolved_extension_map = dict(DEFAULT_EXTENSION_MAP)
     if extension_map is not None:
         resolved_extension_map.update(extension_map)
-    # try to read recording object to get the analyzer
-    try:
+    # The recording is attached only when the caller names its ElectricalSeries. Guessing it would attach
+    # whatever single series the file holds (e.g. LFP), and a failed load would silently fall back to the
+    # recordingless case. electrical_series_path=None is the explicit opt-out.
+    if electrical_series_path is not None:
         recording = NwbRecordingExtractor(
             file_path=file_path,
             electrical_series_path=electrical_series_path,
@@ -2044,12 +2047,12 @@ def read_nwb_sorting_analyzer(
             storage_options=storage_options,
             use_pynwb=use_pynwb,
         )
-    except Exception:
-        if verbose:
-            print("Could not load recording, proceeding without it")
+        # the ElectricalSeries provides the time base, so t_start is passed through as given and the
+        # sorting raises if the caller also set it
+        t_start_tmp = t_start
+    else:
         recording = None
-
-    t_start_tmp = 0 if t_start is None else t_start
+        t_start_tmp = 0 if t_start is None else t_start
 
     sorting_tmp = NwbSortingExtractor(
         file_path=file_path,
@@ -2179,9 +2182,18 @@ def read_nwb_sorting_analyzer(
         # last spike; doubling generously covers the gap while over-estimating only adds a harmless empty
         # tail. Read from sorting_tmp because a group selection wraps `sorting` in a UnitsSelectionSorting
         # whose segment has no spike_times_data.
-        spike_times_data = sorting_tmp._sorting_segments[0].spike_times_data
-        last_stored_spike_time = 0.0 if spike_times_data.shape[0] == 0 else float(np.asarray(spike_times_data[-1]))
-        placeholder_duration = max(last_stored_spike_time * 2.0, 1.0)
+        if recording_duration is not None:
+            placeholder_duration = recording_duration
+        else:
+            spike_times_data = sorting_tmp._sorting_segments[0].spike_times_data
+            last_stored_spike_time = 0.0 if spike_times_data.shape[0] == 0 else float(np.asarray(spike_times_data[-1]))
+            placeholder_duration = max(last_stored_spike_time * 2.0, 1.0)
+            warnings.warn(
+                "No recording attached and no `recording_duration` given, so the recording duration is estimated. "
+                "Recomputing metrics that depend on it (firing_rate, presence_ratio, isi_violation, rp_violation, "
+                "sliding_rp_violation, firing_range, amplitude_cv, drift) will give wrong values. The metrics "
+                "stored in the file are unaffected. Pass `recording_duration` (in seconds) if it is known."
+            )
         analyzer_recording, analyzer_channel_ids = _make_placeholder_recording_from_electrodes(
             sorting, electrodes_table, electrodes_indices, duration=placeholder_duration, verbose=verbose
         )
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index dda3cfbbef..65408b6c63 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1313,7 +1313,9 @@ def test_read_nwb_sorting_analyzer_default(tmp_path):
     from spikeinterface.extractors import read_nwb_sorting_analyzer
 
     path = _make_units_nwb(tmp_path / "units.nwb")
-    analyzer = read_nwb_sorting_analyzer(path, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None)
+    analyzer = read_nwb_sorting_analyzer(
+        path, electrical_series_path=None, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None
+    )
 
     # recordingless, sparse, with templates for every unit
     assert not analyzer.has_recording()
@@ -1337,6 +1339,7 @@ def test_read_nwb_sorting_analyzer_extension_map_override(tmp_path):
     path = _make_units_nwb(tmp_path / "units.nwb")
     analyzer = read_nwb_sorting_analyzer(
         path,
+        electrical_series_path=None,
         use_pynwb=False,
         sampling_frequency=30000.0,
         compute_extra=None,
@@ -1354,6 +1357,7 @@ def test_read_nwb_sorting_analyzer_extension_map_disable(tmp_path):
     path = _make_units_nwb(tmp_path / "units.nwb")
     analyzer = read_nwb_sorting_analyzer(
         path,
+        electrical_series_path=None,
         use_pynwb=False,
         sampling_frequency=30000.0,
         compute_extra=None,
@@ -1369,7 +1373,9 @@ def test_read_nwb_sorting_analyzer_waveform_sd(tmp_path):
     from spikeinterface.extractors import read_nwb_sorting_analyzer
 
     path = _make_units_nwb(tmp_path / "units.nwb", with_std=True)
-    analyzer = read_nwb_sorting_analyzer(path, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None)
+    analyzer = read_nwb_sorting_analyzer(
+        path, electrical_series_path=None, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None
+    )
     templates = analyzer.get_extension("templates")
     # the std operator is populated only because the file stores waveform_sd
     assert "std" in templates.params["operators"]
@@ -1392,12 +1398,31 @@ def logging_getitem(self, key):
         return original_getitem(self, key)
 
     monkeypatch.setattr(h5py.Dataset, "__getitem__", logging_getitem)
-    analyzer = read_nwb_sorting_analyzer(path, use_pynwb=use_pynwb, sampling_frequency=30000.0, compute_extra=None)
+    analyzer = read_nwb_sorting_analyzer(
+        path, electrical_series_path=None, use_pynwb=use_pynwb, sampling_frequency=30000.0, compute_extra=None
+    )
 
     assert "/units/spike_amplitudes_uV" not in read_datasets
     assert "spike_amplitudes_uV" not in analyzer.sorting.get_property_keys()
 
 
+def test_read_nwb_sorting_analyzer_recording_duration(tmp_path):
+    pytest.importorskip("pynwb")
+    import warnings
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    path = _make_units_nwb(tmp_path / "units.nwb")
+    with pytest.warns(UserWarning, match="duration is estimated"):
+        read_nwb_sorting_analyzer(path, electrical_series_path=None, sampling_frequency=30000.0, compute_extra=None)
+
+    with warnings.catch_warnings():
+        warnings.filterwarnings("error", message=".*duration is estimated.*")
+        analyzer = read_nwb_sorting_analyzer(
+            path, electrical_series_path=None, sampling_frequency=30000.0, recording_duration=120.0, compute_extra=None
+        )
+    assert analyzer.get_total_duration() == pytest.approx(120.0)
+
+
 if __name__ == "__main__":
     tmp_path = Path("tmp")
     if tmp_path.is_dir():

From 56fd24d6daed2a5e1f9fe3d33622f910d2a30b1b Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Tue, 6 Oct 2026 14:10:08 -0600
Subject: [PATCH 32/46] Add a docstring to read_nwb_sorting_analyzer

---
 .../extractors/nwbextractors.py               | 70 +++++++++++++++++++
 1 file changed, 70 insertions(+)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 3627fcf7e6..6c6b014e20 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2030,6 +2030,76 @@ def read_nwb_sorting_analyzer(
     rescale_templates_to_uV: bool = True,
     verbose: bool = False,
 ) -> SortingAnalyzer:
+    """
+    Build a SortingAnalyzer from the Units table of an NWB file.
+
+    The analyzer is assembled from what the file stores, not recomputed: templates come from
+    `waveform_mean` (and `waveform_sd`), sparsity from the per-unit `electrodes` region, quality metrics
+    from the canonical quality-metric columns, and every other per-unit scalar column becomes a sorting
+    property. Per-spike columns (e.g. spike amplitudes) are not read.
+
+    Parameters
+    ----------
+    file_path : str or Path
+        Path to the NWB file.
+    electrical_series_path : str or None
+        Path of the ElectricalSeries holding the recording the sorting was run on. It has no default, so
+        attaching a recording is always an explicit choice. When given, that recording is attached to the
+        analyzer and provides the channel ids, the channel groups (for `group_name`) and the time base
+        (`sampling_frequency` and `t_start`, so neither can be passed as well). When None, the analyzer
+        has no recording: the probe is built from the electrodes table and the recording duration is
+        taken from `recording_duration`, or estimated (see `recording_duration`).
+    t_start : float or None, default: None
+        Time (in seconds, on the NWB session clock) of the recording's first sample. Only used when
+        `electrical_series_path` is None, where it defaults to 0.
+    sampling_frequency : float or None, default: None
+        The sampling frequency in Hz. Only used when `electrical_series_path` is None. If None, it is read
+        from the Units table ``resolution`` attribute.
+    recording_duration : float or None, default: None
+        Duration of the recording in seconds, used when `electrical_series_path` is None. If None, the
+        duration is estimated from the spike times and a warning is raised: the metrics stored in the file
+        are unaffected, but recomputing metrics that depend on the duration (e.g. firing_rate,
+        presence_ratio, isi_violation) will give wrong values.
+    unit_table_path : str or None, default: None
+        The path of the Units table in the NWB file. If None, the canonical ``/units`` table is used.
+    stream_mode : "fsspec" | "remfile" | "zarr" | None, default: None
+        The streaming mode to use. If None it assumes the file is on the local disk.
+    stream_cache_path : str or Path or None, default: None
+        Local path for caching. If None it uses the system temporary directory.
+    cache : bool, default: False
+        If True, the file is cached in the file passed to stream_cache_path.
+    storage_options : dict or None, default: None
+        Additional kwargs (e.g. AWS credentials) passed to the zarr.open convenience function.
+        Only used with the "zarr" stream_mode.
+    use_pynwb : bool, default: False
+        Uses the pynwb library to read the NWB file. If False, h5py (or zarr) is used directly.
+    group_name : str or None, default: None
+        The electrode group whose units are loaded. Required when the units span more than one group.
+    compute_extra : list of str or None, default: ["unit_locations"]
+        Extensions to compute on the analyzer after it is built.
+    compute_extra_params : dict or None, default: None
+        Parameters passed to `analyzer.compute()` for `compute_extra`.
+    extension_map : dict or None, default: None
+        Per-extension overrides of `DEFAULT_EXTENSION_MAP`, which maps each extension to where its data
+        lives in the file. Setting an extension to None disables it.
+    rescale_templates_to_uV : bool, default: True
+        Infer the scale of `waveform_mean` from its magnitude and convert the templates to microvolts.
+    verbose : bool, default: False
+        If True, print progress messages.
+
+    Returns
+    -------
+    sorting_analyzer : SortingAnalyzer
+        An in-memory SortingAnalyzer with the `random_spikes` and `templates` extensions (when the file
+        stores templates), `quality_metrics` (when it stores canonical metric columns) and `compute_extra`.
+
+    Notes
+    -----
+    `random_spikes` is set to cover every spike, because it is a required dependency of `templates`.
+    Without a recording, extensions that need traces (e.g. waveforms) cannot be computed. With a
+    recording they can, but recompute `random_spikes` first to extract waveforms from a subset of spikes,
+    and note that recomputing `templates` replaces the templates read from the file.
+    """
     # extension_map overrides (per extension) merge over DEFAULT_EXTENSION_MAP; see its docstring.
     resolved_extension_map = dict(DEFAULT_EXTENSION_MAP)
     if extension_map is not None:

From 532ae4ee148a11789a43b08fae2e75cf0c805cea Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Tue, 6 Oct 2026 14:27:23 -0600
Subject: [PATCH 33/46] Support a non-ragged electrodes column in
 read_nwb_sorting_analyzer and test the electrical series path

---
 .../extractors/nwbextractors.py               |  5 +-
 .../extractors/tests/test_nwbextractors.py    | 68 ++++++++++++++++++-
 2 files changed, 71 insertions(+), 2 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 6c6b014e20..2b39629614 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2210,6 +2210,9 @@ def read_nwb_sorting_analyzer(
             electrodes_indices = electrodes_indices = units["electrodes"][:]
 
     if electrodes_indices is not None:
+        # the schema also allows a non-ragged electrodes column (one electrode row per unit, no
+        # electrodes_index); make every unit's region an array so both layouts are handled the same way
+        electrodes_indices = [np.atleast_1d(region) for region in electrodes_indices]
         # here we assume all groups are the same for each unit, so we just check one.
         if "group_name" in electrodes_table.columns:
             group_names = np.array([electrodes_table.iloc[int(ei[0])]["group_name"] for ei in electrodes_indices])
@@ -2224,7 +2227,7 @@ def read_nwb_sorting_analyzer(
                         print(f"Selecting {sum(unit_mask)} / {len(units)} units from {group_name}")
                     sorting = sorting.select_units(unit_ids=sorting.unit_ids[unit_mask])
                     units = units.loc[units.index[unit_mask]]
-                    electrodes_indices = units["electrodes"]
+                    electrodes_indices = [np.atleast_1d(region) for region in units["electrodes"]]
 
     # Every scalar column not claimed as a quality metric becomes a sorting property (labels like
     # cluster_uuid, and any non-canonical numeric column). Properties show up in the GUI's unit table
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index 65408b6c63..b5ff41fecc 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1266,7 +1266,9 @@ def test_time_series_recording_equality_with_pynwb_and_backend(generate_nwbfile_
     check_recordings_equal(recording_backend, recording_pynwb)
 
 
-def _make_units_nwb(path, n_units=6, n_ch=8, n_samp=30, with_std=False, with_per_spike=False):
+def _make_units_nwb(
+    path, n_units=6, n_ch=8, n_samp=30, with_std=False, with_per_spike=False, with_electrical_series=False
+):
     """Build a minimal NWB file with a Units table for read_nwb_sorting_analyzer tests.
 
     Pure pynwb, no external writers. The Units table has `waveform_mean` (templates), a canonical
@@ -1283,6 +1285,18 @@ def _make_units_nwb(path, n_units=6, n_ch=8, n_samp=30, with_std=False, with_per
     nwbfile.add_electrode_column(name="rel_y", description="y")
     for i in range(n_ch):
         nwbfile.add_electrode(group=group, location="brain", rel_x=float(i % 4) * 20.0, rel_y=float(i // 4) * 20.0)
+    if with_electrical_series:
+        from pynwb.ecephys import ElectricalSeries
+
+        region = nwbfile.create_electrode_table_region(region=list(range(n_ch)), description="all electrodes")
+        series = ElectricalSeries(
+            name="ElectricalSeries",
+            data=np.zeros((1_000, n_ch), dtype="int16"),
+            electrodes=region,
+            rate=30_000.0,
+            starting_time=5.0,
+        )
+        nwbfile.add_acquisition(series)
     nwbfile.add_unit_column(name="snr", description="canonical quality metric")
     nwbfile.add_unit_column(name="custom_score", description="non-canonical per-unit value")
     nwbfile.add_unit_column(name="ks_label", description="curation label")
@@ -1423,6 +1437,58 @@ def test_read_nwb_sorting_analyzer_recording_duration(tmp_path):
     assert analyzer.get_total_duration() == pytest.approx(120.0)
 
 
+@pytest.mark.parametrize("use_pynwb", [True, False])
+def test_read_nwb_sorting_analyzer_non_ragged_electrodes(tmp_path, use_pynwb):
+    """The schema allows a Units `electrodes` column with one electrode row per unit and no index."""
+    pytest.importorskip("pynwb")
+    import h5py
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    n_units, n_ch = 6, 8
+    path = _make_units_nwb(tmp_path / "units.nwb", n_units=n_units, n_ch=n_ch)
+    peak_electrodes = np.arange(n_units) % n_ch
+    # pynwb always writes Units.electrodes with an index, so rewrite it as one row per unit
+    with h5py.File(path, "r+") as file:
+        units = file["units"]
+        attrs = dict(units["electrodes"].attrs)
+        del units["electrodes"], units["electrodes_index"]
+        units.create_dataset("electrodes", data=peak_electrodes)
+        units["electrodes"].attrs.update(attrs)
+
+    analyzer = read_nwb_sorting_analyzer(
+        path, electrical_series_path=None, use_pynwb=use_pynwb, sampling_frequency=30000.0, compute_extra=None
+    )
+    channel_ids = analyzer.channel_ids
+    for unit_id, peak_electrode in zip(analyzer.unit_ids, peak_electrodes):
+        assert list(analyzer.sparsity.unit_id_to_channel_ids[unit_id]) == [channel_ids[peak_electrode]]
+
+
+def test_read_nwb_sorting_analyzer_electrical_series_path(tmp_path):
+    pytest.importorskip("pynwb")
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    path = _make_units_nwb(tmp_path / "units.nwb", with_electrical_series=True)
+
+    # naming the series attaches it and takes the time base from it
+    analyzer = read_nwb_sorting_analyzer(
+        path, electrical_series_path="acquisition/ElectricalSeries", compute_extra=None
+    )
+    assert analyzer.has_recording()
+    assert analyzer.sampling_frequency == 30_000.0
+
+    # the series is the only time base, so an explicit t_start conflicts with it
+    with pytest.raises(ValueError):
+        read_nwb_sorting_analyzer(
+            path, electrical_series_path="acquisition/ElectricalSeries", t_start=0.0, compute_extra=None
+        )
+
+    # None opts out even when the file holds a series
+    analyzer = read_nwb_sorting_analyzer(
+        path, electrical_series_path=None, sampling_frequency=30_000.0, recording_duration=10.0, compute_extra=None
+    )
+    assert not analyzer.has_recording()
+
+
 if __name__ == "__main__":
     tmp_path = Path("tmp")
     if tmp_path.is_dir():

From 838ab240b3e97bafc7bbc9b70f4ec1d0bc395637 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 12:59:32 -0600
Subject: [PATCH 34/46] another joe fix

---
 .../extractors/nwbextractors.py               | 13 +++++
 .../extractors/tests/test_nwbextractors.py    | 48 +++++++++++++++++--
 2 files changed, 57 insertions(+), 4 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 2b39629614..971027b2eb 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2467,6 +2467,17 @@ def _make_templates(
     from spikeinterface.core.analyzer_extension_core import ComputeTemplates
 
     waveform_mean = np.array([np.asarray(t, dtype="float") for t in units[templates_column].values])
+    # The schema also allows one waveform per unit with no electrode axis (e.g. DANDI 000363). It belongs to
+    # the unit's electrode only when the region has exactly one, otherwise we cannot tell which one it is.
+    single_channel = waveform_mean.ndim == 2
+    if single_channel:
+        if any(len(positions) != 1 for positions in unit_local_channels):
+            warnings.warn(
+                f"`{templates_column}` has one waveform per unit but some units have several electrodes, "
+                "so the electrode of each waveform is unknown. The templates extension is not loaded."
+            )
+            return
+        waveform_mean = waveform_mean[:, :, np.newaxis]
     num_samples_template = waveform_mean.shape[1]
 
     # Scatter each unit's waveform onto its channels. This relies on the alignment contract that
@@ -2502,6 +2513,8 @@ def _make_templates(
     # declare only "average" rather than inventing data.
     if templates_std_column is not None and templates_std_column in units.columns:
         waveform_sd = np.array([np.asarray(t, dtype="float") for t in units[templates_std_column].values])
+        if single_channel:
+            waveform_sd = waveform_sd[:, :, np.newaxis]
         dense_std = np.zeros((analyzer.get_num_units(), num_samples_template, num_channels), dtype="float32")
         for unit_index, positions in enumerate(unit_local_channels):
             k = len(positions)
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index b5ff41fecc..7c4d2927b9 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1267,13 +1267,21 @@ def test_time_series_recording_equality_with_pynwb_and_backend(generate_nwbfile_
 
 
 def _make_units_nwb(
-    path, n_units=6, n_ch=8, n_samp=30, with_std=False, with_per_spike=False, with_electrical_series=False
+    path,
+    n_units=6,
+    n_ch=8,
+    n_samp=30,
+    with_std=False,
+    with_per_spike=False,
+    with_electrical_series=False,
+    single_channel=False,
 ):
     """Build a minimal NWB file with a Units table for read_nwb_sorting_analyzer tests.
 
     Pure pynwb, no external writers. The Units table has `waveform_mean` (templates), a canonical
     quality metric (`snr`), a non-canonical numeric column (`custom_score`), a string label
     (`ks_label`), and a per-unit `electrodes` region; the electrodes table carries `rel_x`/`rel_y`.
+    With `single_channel`, each unit has one electrode (`unit % n_ch`) and a 1D `waveform_mean`.
     """
     from datetime import datetime
     from pynwb import NWBFile, NWBHDF5IO
@@ -1303,17 +1311,18 @@ def _make_units_nwb(
     if with_per_spike:
         nwbfile.add_unit_column(name="spike_amplitudes_uV", description="per-spike amplitudes", index=True)
     rng = np.random.default_rng(0)
+    waveform_shape = (n_samp,) if single_channel else (n_samp, n_ch)
     for u in range(n_units):
         kwargs = dict(
             spike_times=np.sort(rng.uniform(0, 100, size=int(rng.integers(50, 500)))),
-            electrodes=list(range(n_ch)),
-            waveform_mean=rng.standard_normal((n_samp, n_ch)),
+            electrodes=[u % n_ch] if single_channel else list(range(n_ch)),
+            waveform_mean=rng.standard_normal(waveform_shape),
             snr=float(rng.uniform(2, 10)),
             custom_score=float(rng.uniform(0, 1)),
             ks_label="good" if u % 2 == 0 else "mua",
         )
         if with_std:
-            kwargs["waveform_sd"] = np.abs(rng.standard_normal((n_samp, n_ch)))
+            kwargs["waveform_sd"] = np.abs(rng.standard_normal(waveform_shape))
         if with_per_spike:
             kwargs["spike_amplitudes_uV"] = rng.standard_normal(kwargs["spike_times"].size)
         nwbfile.add_unit(**kwargs)
@@ -1463,6 +1472,37 @@ def test_read_nwb_sorting_analyzer_non_ragged_electrodes(tmp_path, use_pynwb):
         assert list(analyzer.sparsity.unit_id_to_channel_ids[unit_id]) == [channel_ids[peak_electrode]]
 
 
+@pytest.mark.parametrize("use_pynwb", [True, False])
+def test_read_nwb_sorting_analyzer_single_channel_waveforms(tmp_path, use_pynwb):
+    """The schema allows a 1D `waveform_mean` per unit, recorded on the unit's single electrode."""
+    pytest.importorskip("pynwb")
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    n_units, n_ch = 6, 8
+    path = _make_units_nwb(tmp_path / "units.nwb", n_units=n_units, n_ch=n_ch, with_std=True, single_channel=True)
+    analyzer = read_nwb_sorting_analyzer(
+        path,
+        electrical_series_path=None,
+        use_pynwb=use_pynwb,
+        sampling_frequency=30000.0,
+        rescale_templates_to_uV=False,
+        compute_extra=None,
+    )
+
+    from pynwb import NWBHDF5IO
+
+    with NWBHDF5IO(path, "r") as io:
+        units = io.read().units
+        stored_mean = [units["waveform_mean"][unit_index] for unit_index in range(n_units)]
+    templates = analyzer.get_extension("templates")
+    assert "std" in templates.params["operators"]
+    average = templates.data["average"]
+    for unit_index in range(n_units):
+        electrode = unit_index % n_ch
+        np.testing.assert_allclose(average[unit_index][:, electrode], stored_mean[unit_index], rtol=1e-6)
+        assert not np.any(np.delete(average[unit_index], electrode, axis=1))
+
+
 def test_read_nwb_sorting_analyzer_electrical_series_path(tmp_path):
     pytest.importorskip("pynwb")
     from spikeinterface.extractors import read_nwb_sorting_analyzer

From 0e5a391c608fc3fef9b84ffaaa46a7790203494d Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 14:14:52 -0600
Subject: [PATCH 35/46] Add nwb_to_analyzer_mapping to
 read_nwb_sorting_analyzer to read files that store sparsity and channel
 locations elsewhere

---
 .../extractors/nwbextractors.py               | 299 ++++++++++++------
 .../extractors/tests/test_nwbextractors.py    |  90 +++++-
 2 files changed, 297 insertions(+), 92 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 971027b2eb..a7120f8499 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1939,23 +1939,34 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
     return outputs
 
 
-# Where each SortingAnalyzer extension's data lives in an NWB file. Each extension maps to an ordered
-# list of read-sources tried until one resolves: the typed ndx-spikesorting container first, then the
-# generic-NWB convention. A user/dataset override merges over this (see `extension_map` in
-# `read_nwb_sorting_analyzer`); `None` disables an extension. This is only for real analyzer extensions;
-# sparsity, sorting properties, and recording metadata are handled separately by the reader.
-#
-# The "typed_container" source is a hook: the typed reader currently lives in ndx-spikesorting, so for
-# now it never resolves here and the convention source is used. When that logic moves into
-# spikeinterface, this source will resolve against the ndx-spikesorting types.
-DEFAULT_EXTENSION_MAP = {
+# Where each piece of a SortingAnalyzer lives in an NWB file. Each key is the analyzer object being loaded
+# and maps to an ordered list of places to look, tried until one resolves. A user mapping merges over this
+# one key by key (see `nwb_to_analyzer_mapping` in `read_nwb_sorting_analyzer`), and `None` disables a key.
+# Each place is a dict whose "source" names the kind of place and whose other keys say which one:
+#   * "units_columns": columns of the Units table
+#   * "electrodes_columns": columns of the electrodes table
+#   * "peak_channel_group": a Units column holding each unit's peak electrode id; the unit covers every
+#     electrode in that electrode's group, ordered by the electrodes column "order_by"
+#   * "typed_container": an ndx-spikesorting object. This is a hook: the typed reader lives in
+#     ndx-spikesorting for now, so it never resolves here and the next place is used.
+DEFAULT_NWB_TO_ANALYZER_MAPPING = {
+    # which electrodes each unit's waveform covers (sparsity, and where the templates are placed)
+    "sparsity": [
+        {"source": "units_columns", "electrodes": "electrodes"},
+    ],
+    # the probe contact positions, used when no recording is attached
+    "channel_locations": [
+        {"source": "electrodes_columns", "x": "rel_x", "y": "rel_y"},
+    ],
     "templates": [
         {"source": "typed_container", "type": "Templates"},
-        {"source": "column", "column": "waveform_mean", "std_column": "waveform_sd"},
+        {"source": "units_columns", "mean": "waveform_mean", "std": "waveform_sd"},
     ],
+    # "columns" is "canonical_quality" (the columns named like a SpikeInterface quality metric), a list of
+    # column names, or a dict from SpikeInterface metric name to column name
     "quality_metrics": [
         {"source": "typed_container", "type": "UnitsMetrics"},
-        {"source": "columns", "columns": "canonical_quality"},
+        {"source": "units_columns", "columns": "canonical_quality"},
     ],
     "template_metrics": [
         {"source": "typed_container", "type": "UnitsMetrics"},
@@ -1966,51 +1977,100 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
 }
 
 
-def _resolve_extension_columns(extension_map, scalar_colnames, all_colnames):
-    """Resolve the convention (column-based) sources of an extension map against a Units table's columns.
+def _mapping_structural_colnames(mapping):
+    """The Units columns that the mapping reads as structure (the electrodes region and the waveforms), so
+    they are not classified as metrics or properties and are read even when ragged."""
+    structural = set()
+    for key in ("sparsity", "templates"):
+        for place in mapping.get(key) or []:
+            if place.get("source") == "units_columns":
+                structural.update(place.get(role) for role in ("electrodes", "mean", "std") if place.get(role))
+    return structural
+
 
-    Returns a dict with the resolved template column (and optional std column) and the per-extension
-    column lists for the column-backed extensions. The "typed_container" source is skipped here (the
-    typed reader lives in ndx-spikesorting for now), so resolution falls through to the "column"/
-    "columns" sources. `scalar_colnames` are the 1-D per-unit columns eligible to be metrics/properties;
-    `all_colnames` is every column name (used to check a source's target exists).
+def _resolve_nwb_to_analyzer_mapping(mapping, scalar_colnames, all_colnames):
+    """Resolve the Units-table places of the mapping against the columns of the Units table.
+
+    Returns which place gives each unit's electrodes, which columns hold the templates, which columns are
+    quality metrics (as a dict from metric name to column) and which scalar columns become sorting
+    properties. "typed_container" places are skipped (see `DEFAULT_NWB_TO_ANALYZER_MAPPING`).
+    `scalar_colnames` are the 1D per-unit columns eligible to be metrics or properties, `all_colnames` is
+    every column name.
     """
     from spikeinterface.metrics.quality import ComputeQualityMetrics
 
-    resolved = {"templates_column": None, "templates_std_column": None, "quality_metric_colnames": []}
-
-    for extension, sources in extension_map.items():
-        if sources is None:  # explicitly disabled
+    resolved = {
+        "unit_electrodes": None,
+        "templates_column": None,
+        "templates_std_column": None,
+        "quality_metric_columns": {},
+    }
+
+    for place in mapping.get("sparsity") or []:
+        kind = place.get("source")
+        if kind == "units_columns" and place.get("electrodes") in all_colnames:
+            resolved["unit_electrodes"] = place
+            break
+        if kind == "peak_channel_group" and place.get("peak_channel") in all_colnames:
+            resolved["unit_electrodes"] = place
+            break
+
+    for place in mapping.get("templates") or []:
+        if place.get("source") == "units_columns" and place.get("mean") in all_colnames:
+            resolved["templates_column"] = place["mean"]
+            std_column = place.get("std")
+            resolved["templates_std_column"] = std_column if std_column in all_colnames else None
+            break
+
+    for place in mapping.get("quality_metrics") or []:
+        if place.get("source") != "units_columns":
             continue
-        for source in sources:
-            kind = source.get("source")
-            if kind == "typed_container":
-                continue  # hook: not resolvable here yet
-            if kind == "column":
-                column = source.get("column")
-                if column in all_colnames:
-                    if extension == "templates":
-                        resolved["templates_column"] = column
-                        std_column = source.get("std_column")
-                        resolved["templates_std_column"] = std_column if std_column in all_colnames else None
-                    break
-            if kind == "columns":
-                spec = source.get("columns")
-                if spec == "canonical_quality":
-                    canonical = set(ComputeQualityMetrics.get_metric_columns())
-                    cols = [c for c in scalar_colnames if c in canonical]
-                else:  # an explicit list of column names
-                    cols = [c for c in spec if c in scalar_colnames]
-                if extension == "quality_metrics":
-                    resolved["quality_metric_colnames"] = cols
-                break
-
-    # every scalar column not claimed as a quality metric becomes a sorting property
-    claimed = set(resolved["quality_metric_colnames"])
+        spec = place.get("columns")
+        if spec == "canonical_quality":
+            canonical = set(ComputeQualityMetrics.get_metric_columns())
+            columns = {c: c for c in scalar_colnames if c in canonical}
+        elif isinstance(spec, dict):
+            columns = {metric: column for metric, column in spec.items() if column in scalar_colnames}
+        else:  # a list of column names, kept under their own names
+            columns = {c: c for c in spec if c in scalar_colnames}
+        resolved["quality_metric_columns"] = columns
+        break
+
+    # every scalar column not read as a quality metric becomes a sorting property
+    claimed = set(resolved["quality_metric_columns"].values())
     resolved["property_colnames"] = [c for c in scalar_colnames if c not in claimed]
     return resolved
 
 
+def _resolve_channel_location_columns(mapping, electrodes_colnames):
+    """The (x, y) electrodes-table columns of the first "channel_locations" place present in the table."""
+    for place in mapping.get("channel_locations") or []:
+        if place.get("source") == "electrodes_columns":
+            if place.get("x") in electrodes_colnames and place.get("y") in electrodes_colnames:
+                return place["x"], place["y"]
+    return None
+
+
+def _electrode_rows_from_peak_channel_group(peak_electrode_ids, electrodes_table, order_by):
+    """For each unit, the electrodes-table rows of every electrode in the group of its peak electrode,
+    ordered by the `order_by` column (table order when `order_by` is None)."""
+    if "group_name" not in electrodes_table.columns:
+        raise ValueError("A 'peak_channel_group' place needs a 'group_name' column in the electrodes table.")
+    if order_by is not None and order_by not in electrodes_table.columns:
+        raise ValueError(f"The 'order_by' column {order_by!r} is not in the electrodes table.")
+    peak_rows = electrodes_table.index.get_indexer(np.asarray(peak_electrode_ids))
+    if np.any(peak_rows < 0):
+        raise ValueError("Some peak electrode ids are not ids of the electrodes table.")
+    group_names = electrodes_table["group_name"].to_numpy()
+    rows_of_group = {}
+    for group in np.unique(group_names):
+        rows = np.flatnonzero(group_names == group)
+        if order_by is not None:
+            rows = rows[np.argsort(electrodes_table[order_by].to_numpy()[rows], kind="stable")]
+        rows_of_group[group] = rows
+    return [rows_of_group[group_names[row]] for row in peak_rows]
+
+
 def read_nwb_sorting_analyzer(
     file_path: str | Path,
     electrical_series_path: str | None,
@@ -2026,7 +2086,7 @@ def read_nwb_sorting_analyzer(
     group_name: str | None = None,
     compute_extra: List[str] | None = ["unit_locations"],
     compute_extra_params: dict | None = None,
-    extension_map: dict | None = None,
+    nwb_to_analyzer_mapping: dict | None = None,
     rescale_templates_to_uV: bool = True,
     verbose: bool = False,
 ) -> SortingAnalyzer:
@@ -2036,7 +2096,8 @@ def read_nwb_sorting_analyzer(
     The analyzer is assembled from what the file stores, not recomputed: templates come from
     `waveform_mean` (and `waveform_sd`), sparsity from the per-unit `electrodes` region, quality metrics
     from the canonical quality-metric columns, and every other per-unit scalar column becomes a sorting
-    property. Per-spike columns (e.g. spike amplitudes) are not read.
+    property. Per-spike columns (e.g. spike amplitudes) are not read. Files that store these elsewhere can
+    be read by passing `nwb_to_analyzer_mapping` (see the example below).
 
     Parameters
     ----------
@@ -2079,9 +2140,11 @@ def read_nwb_sorting_analyzer(
         Extensions to compute on the analyzer after it is built.
     compute_extra_params : dict or None, default: None
         Parameters passed to `analyzer.compute()` for `compute_extra`.
-    extension_map : dict or None, default: None
-        Per-extension overrides of `DEFAULT_EXTENSION_MAP`, which maps each extension to where its data
-        lives in the file. Setting an extension to None disables it.
+    nwb_to_analyzer_mapping : dict or None, default: None
+        Where each analyzer object lives in the file, merged key by key over
+        `DEFAULT_NWB_TO_ANALYZER_MAPPING`. Each key is an analyzer object ("sparsity", "channel_locations",
+        "templates", "quality_metrics", ...) and maps to a list of places tried in order. Setting a key to
+        None disables it.
     rescale_templates_to_uV : bool, default: True
         Infer the scale of `waveform_mean` from its magnitude and convert the templates to microvolts.
     verbose : bool, default: False
@@ -2099,11 +2162,34 @@ def read_nwb_sorting_analyzer(
     Without a recording, extensions that need traces (e.g. waveforms) cannot be computed. With a
     recording they can, but recompute `random_spikes` first to extract waveforms from a subset of spikes,
     and note that recomputing `templates` replaces the templates read from the file.
+
+    Examples
+    --------
+    The Allen Institute Neuropixels files on DANDI (e.g. 000021, 000253) have no Units `electrodes` column.
+    Each unit stores its peak electrode id in `peak_channel_id` and its waveform on every electrode of that
+    probe, the electrode positions are in `probe_horizontal_position` and `probe_vertical_position`, and
+    some metric columns use other names:
+
+    >>> allen_mapping = {
+    ...     "sparsity": [{"source": "peak_channel_group", "peak_channel": "peak_channel_id", "order_by": "local_index"}],
+    ...     "channel_locations": [
+    ...         {"source": "electrodes_columns", "x": "probe_horizontal_position", "y": "probe_vertical_position"}
+    ...     ],
+    ...     "quality_metrics": [{"source": "units_columns", "columns": {
+    ...         "amplitude_cutoff": "amplitude_cutoff", "d_prime": "d_prime", "firing_rate": "firing_rate",
+    ...         "isi_violations_ratio": "isi_violations", "isolation_distance": "isolation_distance",
+    ...         "l_ratio": "l_ratio", "nn_hit_rate": "nn_hit_rate", "nn_miss_rate": "nn_miss_rate",
+    ...         "presence_ratio": "presence_ratio", "silhouette": "silhouette_score", "snr": "snr",
+    ...     }}],
+    ... }
+    >>> analyzer = read_nwb_sorting_analyzer(
+    ...     file_path, electrical_series_path=None, group_name="probeA", nwb_to_analyzer_mapping=allen_mapping
+    ... )  # doctest: +SKIP
     """
-    # extension_map overrides (per extension) merge over DEFAULT_EXTENSION_MAP; see its docstring.
-    resolved_extension_map = dict(DEFAULT_EXTENSION_MAP)
-    if extension_map is not None:
-        resolved_extension_map.update(extension_map)
+    # the user mapping merges key by key over DEFAULT_NWB_TO_ANALYZER_MAPPING; see its comment
+    mapping = dict(DEFAULT_NWB_TO_ANALYZER_MAPPING)
+    if nwb_to_analyzer_mapping is not None:
+        mapping.update(nwb_to_analyzer_mapping)
     # The recording is attached only when the caller names its ElectricalSeries. Guessing it would attach
     # whatever single series the file holds (e.g. LFP), and a failed load would silently fall back to the
     # recordingless case. electrical_series_path=None is the explicit opt-out.
@@ -2156,9 +2242,10 @@ def read_nwb_sorting_analyzer(
         colnames = list(units_table.colnames)
         # ragged columns (per-spike data) are excluded so to_dataframe never reads them; electrodes is
         # ragged too but is the per-unit region needed for sparsity
-        ragged_colnames = {c for c in colnames if c != "electrodes" and isinstance(units_table[c], VectorIndex)}
+        # the structural columns (electrodes region, waveforms) are kept even when ragged
+        structural = _mapping_structural_colnames(mapping)
+        ragged_colnames = {c for c in colnames if c not in structural and isinstance(units_table[c], VectorIndex)}
         units = units_table.to_dataframe(index=True, exclude=ragged_colnames)
-        structural = {"waveform_mean", "waveform_sd", "electrodes"}
         metric_colnames = [c for c in units.columns if c not in structural and units[c].dtype.kind in "fiu"]
         label_colnames = [
             c
@@ -2170,7 +2257,7 @@ def read_nwb_sorting_analyzer(
     else:
         units_group = sorting.units_table
         colnames = list(units_group.keys())
-        structural = ("id", "spike_times", "waveform_mean", "waveform_sd", "electrodes")
+        structural = {"id", "spike_times"} | _mapping_structural_colnames(mapping)
         metric_colnames, label_colnames = [], []
         for column_name in colnames:
             if column_name.endswith("_index") or column_name in structural or f"{column_name}_index" in colnames:
@@ -2183,36 +2270,43 @@ def read_nwb_sorting_analyzer(
             elif dataset.ndim == 1 and dataset.dtype.kind in "OSU":
                 label_colnames.append(column_name)
 
-    # Resolve the extension map against the Units columns: which column holds the templates, and which
-    # scalar columns are quality metrics (canonical names by default) vs plain sorting properties.
+    # Resolve the mapping against the Units columns: where each unit's electrodes come from, which column
+    # holds the templates, and which scalar columns are quality metrics vs plain sorting properties.
     scalar_colnames = metric_colnames + label_colnames
-    resolved = _resolve_extension_columns(resolved_extension_map, scalar_colnames, colnames)
+    resolved = _resolve_nwb_to_analyzer_mapping(mapping, scalar_colnames, colnames)
+    unit_electrodes = resolved["unit_electrodes"]
     templates_column = resolved["templates_column"]
     templates_std_column = resolved["templates_std_column"]
-    quality_metric_colnames = resolved["quality_metric_colnames"]
+    quality_metric_columns = resolved["quality_metric_columns"]
     property_colnames = resolved["property_colnames"]
+    electrodes_column = None
+    if unit_electrodes is not None and unit_electrodes["source"] == "units_columns":
+        electrodes_column = unit_electrodes["electrodes"]
 
     if not use_pynwb:
         # deliberate read: only the resolved template column(s), the electrodes region, and the scalar
         # columns that become metrics or properties are materialized; ragged/per-spike columns are not.
-        needed = [c for c in (templates_column, templates_std_column, "electrodes") if c and c in colnames]
+        needed = [c for c in (templates_column, templates_std_column, electrodes_column) if c and c in colnames]
         needed += scalar_colnames
         units = _create_df_from_nwb_table(units_group, columns=needed)
 
-    electrodes_indices = None
     if use_pynwb:
         electrodes_table = sorting._reader.nwbfile.electrodes.to_dataframe(index=True)
-        if "electrodes" in colnames:
-            electrodes_indices = units["electrodes"]
     else:
         electrodes_table = _create_df_from_nwb_table(sorting._reader.file["/general/extracellular_ephys/electrodes"])
-        if "electrodes" in colnames:
-            electrodes_indices = electrodes_indices = units["electrodes"][:]
 
-    if electrodes_indices is not None:
+    # each unit's electrodes, as rows of the electrodes table
+    electrodes_indices = None
+    if electrodes_column is not None:
         # the schema also allows a non-ragged electrodes column (one electrode row per unit, no
         # electrodes_index); make every unit's region an array so both layouts are handled the same way
-        electrodes_indices = [np.atleast_1d(region) for region in electrodes_indices]
+        electrodes_indices = [np.atleast_1d(region) for region in units[electrodes_column]]
+    elif unit_electrodes is not None:  # "peak_channel_group"
+        electrodes_indices = _electrode_rows_from_peak_channel_group(
+            units[unit_electrodes["peak_channel"]].to_numpy(), electrodes_table, unit_electrodes.get("order_by")
+        )
+
+    if electrodes_indices is not None:
         # here we assume all groups are the same for each unit, so we just check one.
         if "group_name" in electrodes_table.columns:
             group_names = np.array([electrodes_table.iloc[int(ei[0])]["group_name"] for ei in electrodes_indices])
@@ -2227,7 +2321,7 @@ def read_nwb_sorting_analyzer(
                         print(f"Selecting {sum(unit_mask)} / {len(units)} units from {group_name}")
                     sorting = sorting.select_units(unit_ids=sorting.unit_ids[unit_mask])
                     units = units.loc[units.index[unit_mask]]
-                    electrodes_indices = [np.atleast_1d(region) for region in units["electrodes"]]
+                    electrodes_indices = [region for region, keep in zip(electrodes_indices, unit_mask) if keep]
 
     # Every scalar column not claimed as a quality metric becomes a sorting property (labels like
     # cluster_uuid, and any non-canonical numeric column). Properties show up in the GUI's unit table
@@ -2267,8 +2361,26 @@ def read_nwb_sorting_analyzer(
                 "sliding_rp_violation, firing_range, amplitude_cv, drift) will give wrong values. The metrics "
                 "stored in the file are unaffected. Pass `recording_duration` (in seconds) if it is known."
             )
+        if electrodes_indices is None:
+            raise ValueError(
+                "Without a recording, the channels are taken from the electrodes of each unit, but no 'sparsity' "
+                "place of the mapping was found in the Units table (by default the 'electrodes' column). Set "
+                "'sparsity' in `nwb_to_analyzer_mapping` to where the file stores each unit's electrodes."
+            )
+        location_columns = _resolve_channel_location_columns(mapping, electrodes_table.columns)
+        if location_columns is None:
+            raise ValueError(
+                "No 'channel_locations' place of the mapping was found in the electrodes table (by default the "
+                "'rel_x' and 'rel_y' columns). Set 'channel_locations' in `nwb_to_analyzer_mapping` to the "
+                "columns holding the electrode positions."
+            )
         analyzer_recording, analyzer_channel_ids = _make_placeholder_recording_from_electrodes(
-            sorting, electrodes_table, electrodes_indices, duration=placeholder_duration, verbose=verbose
+            sorting,
+            electrodes_table,
+            electrodes_indices,
+            location_columns,
+            duration=placeholder_duration,
+            verbose=verbose,
         )
 
     # Per-unit sparsity and channel map from the Units `electrodes` region. Each unit's `waveform_mean`
@@ -2319,7 +2431,7 @@ def read_nwb_sorting_analyzer(
         )
 
     # the resolved quality-metric columns -> quality_metrics extension
-    _make_metrics(analyzer, units, quality_metric_colnames, verbose=verbose)
+    _make_metrics(analyzer, units, quality_metric_columns, verbose=verbose)
 
     # compute extra required
     if compute_extra is not None:
@@ -2331,7 +2443,9 @@ def read_nwb_sorting_analyzer(
     return analyzer
 
 
-def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, electrodes_indices, duration, verbose=False):
+def _make_placeholder_recording_from_electrodes(
+    sorting, electrodes_table, electrodes_indices, location_columns, duration, verbose=False
+):
     """Build a lightweight placeholder recording that carries only the probe geometry of the electrodes
     a unit's waveforms live on, for the recordingless (streaming) case.
 
@@ -2361,11 +2475,8 @@ def _make_placeholder_recording_from_electrodes(sorting, electrodes_table, elect
     else:
         channel_ids = electrodes_table_sliced.index.to_numpy()
 
-    electrode_colnames = electrodes_table_sliced.columns
-    assert (
-        "rel_x" in electrode_colnames and "rel_y" in electrode_colnames
-    ), "'rel_x' and 'rel_y' should be columns in the electrodes table"
-    locations = np.array([electrodes_table_sliced["rel_x"][:], electrodes_table_sliced["rel_y"][:]]).T
+    x_column, y_column = location_columns
+    locations = np.array([electrodes_table_sliced[x_column][:], electrodes_table_sliced[y_column][:]]).T
 
     probe = Probe(si_units="um")
     probe.set_contacts(locations, shapes="circle", shape_params={"radius": 1})
@@ -2478,6 +2589,15 @@ def _make_templates(
             )
             return
         waveform_mean = waveform_mean[:, :, np.newaxis]
+    # Some writers store each unit's waveform as electrodes x samples (e.g. the Allen Institute files)
+    # instead of the schema's samples x electrodes. We read it in that orientation when every unit's waveform
+    # has exactly one row per electrode of the unit and a different number of columns.
+    electrodes_by_samples = len(unit_local_channels) > 0 and all(
+        waveform_mean.shape[1] == len(positions) and waveform_mean.shape[2] != len(positions)
+        for positions in unit_local_channels
+    )
+    if electrodes_by_samples:
+        waveform_mean = waveform_mean.transpose(0, 2, 1)
     num_samples_template = waveform_mean.shape[1]
 
     # Scatter each unit's waveform onto its channels. This relies on the alignment contract that
@@ -2515,6 +2635,8 @@ def _make_templates(
         waveform_sd = np.array([np.asarray(t, dtype="float") for t in units[templates_std_column].values])
         if single_channel:
             waveform_sd = waveform_sd[:, :, np.newaxis]
+        if electrodes_by_samples:
+            waveform_sd = waveform_sd.transpose(0, 2, 1)
         dense_std = np.zeros((analyzer.get_num_units(), num_samples_template, num_channels), dtype="float32")
         for unit_index, positions in enumerate(unit_local_channels):
             k = len(positions)
@@ -2533,17 +2655,18 @@ def _make_templates(
     analyzer.extensions["templates"] = templates_ext
 
 
-def _make_metrics(analyzer, units, quality_metric_colnames, verbose=False):
-    """Attach the `quality_metrics` extension from the resolved quality-metric columns. Other per-unit
+def _make_metrics(analyzer, units, quality_metric_columns, verbose=False):
+    """Attach the `quality_metrics` extension from the resolved quality-metric columns, a dict from metric
+    name to Units column. Other per-unit
     scalar columns are loaded as sorting properties instead (see the caller), and `template_metrics` is
     never synthesized from value columns because it carries structured data a values-only table lacks."""
     import pandas as pd
     from spikeinterface.metrics.quality import ComputeQualityMetrics
 
     quality_metric_df = pd.DataFrame(index=analyzer.unit_ids)
-    for col in quality_metric_colnames:
-        if col in units.columns:
-            quality_metric_df.loc[:, col] = units[col].values
+    for metric_name, column in quality_metric_columns.items():
+        if column in units.columns:
+            quality_metric_df.loc[:, metric_name] = units[column].values
 
     if len(quality_metric_df.columns) > 0:
         if verbose:
@@ -2570,8 +2693,8 @@ def _create_df_from_nwb_table(group, columns=None):
         colnames = [c for c in ["id", *columns] if c in all_colnames]
     data = {}
     for col in colnames:
-        if "_index" in col:
-            continue
+        if col.endswith("_index") and col[: -len("_index")] in all_colnames:
+            continue  # the index of a ragged column (e.g. spike_times_index), read with that column
         item = group[col][:]
         if f"{col}_index" in all_colnames:
             item = np.split(item, group[f"{col}_index"][:])[:-1]
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index 7c4d2927b9..87ad7ca044 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1355,7 +1355,7 @@ def test_read_nwb_sorting_analyzer_default(tmp_path):
     assert "snr" not in properties
 
 
-def test_read_nwb_sorting_analyzer_extension_map_override(tmp_path):
+def test_read_nwb_sorting_analyzer_mapping_override(tmp_path):
     pytest.importorskip("pynwb")
     from spikeinterface.extractors import read_nwb_sorting_analyzer
 
@@ -1366,14 +1366,14 @@ def test_read_nwb_sorting_analyzer_extension_map_override(tmp_path):
         use_pynwb=False,
         sampling_frequency=30000.0,
         compute_extra=None,
-        extension_map={"quality_metrics": [{"source": "columns", "columns": ["snr", "custom_score"]}]},
+        nwb_to_analyzer_mapping={"quality_metrics": [{"source": "units_columns", "columns": ["snr", "custom_score"]}]},
     )
     assert set(analyzer.get_extension("quality_metrics").get_data().columns) == {"snr", "custom_score"}
     # custom_score is now a metric, so it is no longer a property
     assert "custom_score" not in analyzer.sorting.get_property_keys()
 
 
-def test_read_nwb_sorting_analyzer_extension_map_disable(tmp_path):
+def test_read_nwb_sorting_analyzer_mapping_disable(tmp_path):
     pytest.importorskip("pynwb")
     from spikeinterface.extractors import read_nwb_sorting_analyzer
 
@@ -1384,7 +1384,7 @@ def test_read_nwb_sorting_analyzer_extension_map_disable(tmp_path):
         use_pynwb=False,
         sampling_frequency=30000.0,
         compute_extra=None,
-        extension_map={"quality_metrics": None},
+        nwb_to_analyzer_mapping={"quality_metrics": None},
     )
     assert not analyzer.has_extension("quality_metrics")
     # with quality_metrics disabled, every scalar column (including snr) becomes a property
@@ -1503,6 +1503,88 @@ def test_read_nwb_sorting_analyzer_single_channel_waveforms(tmp_path, use_pynwb)
         assert not np.any(np.delete(average[unit_index], electrode, axis=1))
 
 
+def test_read_nwb_sorting_analyzer_peak_channel_group_mapping(tmp_path):
+    """A file laid out like the Allen Institute ones: no Units `electrodes` column, a peak electrode id per
+    unit, each waveform stored as electrodes x samples over every electrode of the probe, the positions in
+    other electrodes columns and a metric under another name."""
+    pytest.importorskip("pynwb")
+    import h5py
+    from datetime import datetime
+    from pynwb import NWBFile, NWBHDF5IO
+    from spikeinterface.extractors import read_nwb_sorting_analyzer
+
+    n_probes, n_ch, n_samp, n_units = 2, 6, 20, 5
+    nwbfile = NWBFile(session_description="t", identifier="id", session_start_time=datetime(2020, 1, 1).astimezone())
+    device = nwbfile.create_device(name="probe")
+    for column in ("probe_horizontal_position", "probe_vertical_position", "local_index"):
+        nwbfile.add_electrode_column(name=column, description=column)
+    for probe_index in range(n_probes):
+        group = nwbfile.create_electrode_group(name=f"probe{probe_index}", description="d", location="l", device=device)
+        # rows are stored in reverse local_index order, so the reader has to sort by local_index
+        for local_index in reversed(range(n_ch)):
+            nwbfile.add_electrode(
+                group=group,
+                location="brain",
+                probe_horizontal_position=float(local_index % 2) * 16.0,
+                probe_vertical_position=float(local_index) * 20.0,
+                local_index=local_index,
+            )
+    nwbfile.add_unit_column(name="peak_channel_id", description="electrode id of the peak channel")
+    nwbfile.add_unit_column(name="isi_violations", description="isi violations ratio")
+    rng = np.random.default_rng(0)
+    peak_channel_ids = [1, 3, 7, 8, 11]  # units 0-1 on probe0 (ids 0-5), units 2-4 on probe1 (ids 6-11)
+    for unit_index in range(n_units):
+        nwbfile.add_unit(
+            spike_times=np.sort(rng.uniform(0, 100, size=50)),
+            peak_channel_id=peak_channel_ids[unit_index],
+            isi_violations=float(rng.uniform(0, 1)),
+        )
+    path = tmp_path / "allen_like.nwb"
+    with NWBHDF5IO(str(path), "w") as io:
+        io.write(nwbfile)
+    # pynwb writes waveform_mean as samples x electrodes per unit, so write the ragged electrodes x samples
+    # layout with h5py: row j of a unit's block is the electrode with local_index j of its probe
+    blocks = rng.standard_normal((n_units, n_ch, n_samp))
+    with h5py.File(path, "r+") as file:
+        file["units"].create_dataset("waveform_mean", data=blocks.reshape(n_units * n_ch, n_samp))
+        file["units"].create_dataset("waveform_mean_index", data=np.arange(1, n_units + 1) * n_ch)
+
+    allen_mapping = {
+        "sparsity": [{"source": "peak_channel_group", "peak_channel": "peak_channel_id", "order_by": "local_index"}],
+        "channel_locations": [
+            {"source": "electrodes_columns", "x": "probe_horizontal_position", "y": "probe_vertical_position"}
+        ],
+        "quality_metrics": [{"source": "units_columns", "columns": {"isi_violations_ratio": "isi_violations"}}],
+    }
+    analyzer = read_nwb_sorting_analyzer(
+        str(path),
+        electrical_series_path=None,
+        sampling_frequency=30000.0,
+        recording_duration=100.0,
+        group_name="probe1",
+        nwb_to_analyzer_mapping=allen_mapping,
+        rescale_templates_to_uV=False,
+        compute_extra=None,
+    )
+
+    # the probe1 units cover every probe1 electrode, at the positions of the mapped columns
+    assert analyzer.get_num_units() == 3
+    channel_ids = np.asarray(analyzer.channel_ids).astype(int)
+    assert set(channel_ids) == set(range(6, 12))
+    local_index_of_channel = 11 - channel_ids  # ids 6-11 hold local_index 5-0
+    expected_locations = np.stack([(local_index_of_channel % 2) * 16.0, local_index_of_channel * 20.0], axis=1)
+    np.testing.assert_allclose(analyzer.get_channel_locations(), expected_locations)
+
+    # each waveform row lands on the electrode with that local_index
+    average = analyzer.get_extension("templates").data["average"]
+    for selected_index, unit_index in enumerate([2, 3, 4]):
+        expected = blocks[unit_index][local_index_of_channel].T
+        np.testing.assert_allclose(average[selected_index], expected, rtol=1e-6)
+
+    metrics = analyzer.get_extension("quality_metrics").get_data()
+    assert list(metrics.columns) == ["isi_violations_ratio"]
+
+
 def test_read_nwb_sorting_analyzer_electrical_series_path(tmp_path):
     pytest.importorskip("pynwb")
     from spikeinterface.extractors import read_nwb_sorting_analyzer

From 361a37512b316a06cfde13c57c6bea3ee3089aaf Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 14:25:47 -0600
Subject: [PATCH 36/46] Say why the sampling frequency is unknown when the
 Units table has no resolution

---
 src/spikeinterface/extractors/nwbextractors.py            | 5 +++--
 src/spikeinterface/extractors/tests/test_nwbextractors.py | 5 +++--
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index a7120f8499..71dd909745 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1400,8 +1400,9 @@ def __init__(
                 sampling_frequency = self._reader._sampling_frequency_from_units_metadata()
             if sampling_frequency is None:
                 raise ValueError(
-                    "Couldn't determine the sampling frequency. Provide it with the 'sampling_frequency' "
-                    "argument, set the Units table 'resolution' attribute, or pass 'electrical_series_path'."
+                    "The Units table has no 'resolution' attribute, so the sampling frequency of the spike times "
+                    "is unknown. Pass it with 'sampling_frequency' (e.g. sampling_frequency=30000.0), or pass "
+                    "'electrical_series_path' to take it from the recording the sorting was run on."
                 )
 
         unit_ids = self._reader.unit_ids()
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index 87ad7ca044..7f30274718 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -838,8 +838,9 @@ def test_sorting_sampling_frequency_missing_raises(tmp_path, use_pynwb):
             use_pynwb=use_pynwb,
         )
     expected_error = (
-        "Couldn't determine the sampling frequency. Provide it with the 'sampling_frequency' "
-        "argument, set the Units table 'resolution' attribute, or pass 'electrical_series_path'."
+        "The Units table has no 'resolution' attribute, so the sampling frequency of the spike times "
+        "is unknown. Pass it with 'sampling_frequency' (e.g. sampling_frequency=30000.0), or pass "
+        "'electrical_series_path' to take it from the recording the sorting was run on."
     )
     assert str(exc_info.value) == expected_error
 

From 092d752db7547cb4690798fb7de558610531b118 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 14:35:13 -0600
Subject: [PATCH 37/46] Remove the unused get_last_spike_time

---
 src/spikeinterface/core/basesorting.py        | 17 -----------------
 .../extractors/nwbextractors.py               | 19 -------------------
 2 files changed, 36 deletions(-)

diff --git a/src/spikeinterface/core/basesorting.py b/src/spikeinterface/core/basesorting.py
index f2e66d8817..36a53e51dc 100644
--- a/src/spikeinterface/core/basesorting.py
+++ b/src/spikeinterface/core/basesorting.py
@@ -485,23 +485,6 @@ def get_last_spike_frame(self, segment_index: int | None = None) -> int:
             return 0
         return int(np.max(spikes_in_segment["sample_index"]))
 
-    def get_last_spike_time(self, segment_index: int | None = None) -> float:
-        """Get the time of the last spike in a segment across all units.
-
-        Parameters
-        ----------
-        segment_index : int or None, default: None
-            The segment index (required for multi-segment)
-
-        Returns
-        -------
-        float
-            The time of the last spike in seconds, or 0.0 if no spikes exist.
-        """
-        segment_index = self._check_segment_index(segment_index)
-        last_spike_frame = self.get_last_spike_frame(segment_index=segment_index)
-        return self.sample_index_to_time(last_spike_frame, segment_index=segment_index)
-
     def get_times(
         self,
         segment_index: int | None = None,
diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 71dd909745..ecc1894954 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1595,25 +1595,6 @@ def _get_all_spike_samples_and_unit_indices(self):
         sample_indices = self._times_to_samples(spike_times)
         return sample_indices, unit_indices
 
-    def get_last_spike_time(self, segment_index: int | None = None) -> float:
-        """Get the time of the last spike in a segment across all units.
-        Overridden from BaseSorting to use spike_times_data.
-
-        Parameters
-        ----------
-        segment_index : int or None, default: None
-            The segment index (required for multi-segment)
-
-        Returns
-        -------
-        float
-            The time of the last spike in seconds, or 0.0 if no spikes exist.
-        """
-        segment = self.segments[segment_index] if segment_index is not None else self.segments[0]
-        if segment.spike_times_data.size == 0:
-            return 0.0
-        return segment.spike_times_data[-1]
-
     def get_unit_spike_train_in_seconds(
         self,
         unit_id,

From c05e3a6e854c1360f742bd4ff5d30bc4906c3230 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 14:51:31 -0600
Subject: [PATCH 38/46] Import the core analyzer pieces at the top of
 nwbextractors

---
 src/spikeinterface/extractors/nwbextractors.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index ecc1894954..086296ddcd 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -12,7 +12,11 @@
     BaseRecordingSegment,
     BaseSorting,
     BaseSortingSegment,
+    ChannelSparsity,
+    ComputeRandomSpikes,
+    ComputeTemplates,
     SortingAnalyzer,
+    generate_ground_truth_recording,
     get_default_analyzer_extension_params,
 )
 from spikeinterface.core.base import minimum_spike_dtype
@@ -2438,7 +2442,6 @@ def _make_placeholder_recording_from_electrodes(
     ordered channel ids.
     """
     from probeinterface import Probe
-    from spikeinterface.core import generate_ground_truth_recording
 
     # union of electrode rows referenced by any unit, in sorted order
     electrode_indices_all = []
@@ -2482,8 +2485,6 @@ def _make_sparsity_from_electrodes(sorting, electrodes_table, electrodes_indices
     if electrodes_indices is None:
         return None, None
 
-    from spikeinterface.core.sparsity import ChannelSparsity
-
     num_channels = len(analyzer_channel_ids)
     if "channel_name" in electrodes_table.columns:
         electrode_row_to_channel_id = electrodes_table["channel_name"].to_numpy()
@@ -2517,7 +2518,6 @@ def _make_random_spikes(analyzer, total_num_spikes):
     `method="all"` means the indices are just `arange(total_num_spikes)`; the caller supplies the count
     (read cheaply from `spike_times_index`) so this never materializes the spike vector.
     """
-    from spikeinterface.core.analyzer_extension_core import ComputeRandomSpikes
 
     random_spikes_ext = ComputeRandomSpikes(sorting_analyzer=analyzer)
     random_spikes_ext.set_params(method="all")
@@ -2557,7 +2557,6 @@ def _make_templates(
 ):
     """Attach the `templates` extension from the resolved templates column, scattering each unit's
     sparse waveform block onto its true channel positions."""
-    from spikeinterface.core.analyzer_extension_core import ComputeTemplates
 
     waveform_mean = np.array([np.asarray(t, dtype="float") for t in units[templates_column].values])
     # The schema also allows one waveform per unit with no electrode axis (e.g. DANDI 000363). It belongs to

From 842b82b9cb5512284f0543407e1884c1bc99ffd0 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 15:33:32 -0600
Subject: [PATCH 39/46] Read the spike counts and last spike before the group
 selection instead of keeping sorting_tmp

---
 .../extractors/nwbextractors.py               | 29 ++++++++++---------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 086296ddcd..6726d9c6fb 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2196,7 +2196,7 @@ def read_nwb_sorting_analyzer(
         recording = None
         t_start_tmp = 0 if t_start is None else t_start
 
-    sorting_tmp = NwbSortingExtractor(
+    sorting = NwbSortingExtractor(
         file_path=file_path,
         electrical_series_path=electrical_series_path,
         unit_table_path=unit_table_path,
@@ -2210,7 +2210,16 @@ def read_nwb_sorting_analyzer(
         load_unit_properties=False,  # columns are read deliberately below, into their extensions
     )
 
-    sorting = sorting_tmp
+    # Read the spike count of each unit (from spike_times_index, ~KB) and, when the duration has to be estimated,
+    # the last stored spike time (one chunk) from the NWB arrays now. A group selection below replaces `sorting`
+    # with a UnitsSelectionSorting, which could only get them by building the parent's whole spike vector
+    # (~130-480 MB).
+    sorting_segment = sorting._sorting_segments[0]
+    per_unit_spike_counts = np.diff(sorting_segment.spike_times_index_data[:], prepend=0)
+    last_stored_spike_time = None
+    if recording is None and recording_duration is None:
+        spike_times_data = sorting_segment.spike_times_data
+        last_stored_spike_time = 0.0 if spike_times_data.shape[0] == 0 else float(np.asarray(spike_times_data[-1]))
     # Recordingless case: leave t_start at 0 (set when the sorting was constructed). NWB spike times are
     # in seconds from session start, so they are always >= 0 and map to non-negative frames with
     # t_start = 0. Anchoring t_start to the first spike would only remove empty space before it on the
@@ -2307,6 +2316,7 @@ def read_nwb_sorting_analyzer(
                         print(f"Selecting {sum(unit_mask)} / {len(units)} units from {group_name}")
                     sorting = sorting.select_units(unit_ids=sorting.unit_ids[unit_mask])
                     units = units.loc[units.index[unit_mask]]
+                    per_unit_spike_counts = per_unit_spike_counts[unit_mask]
                     electrodes_indices = [region for region, keep in zip(electrodes_indices, unit_mask) if keep]
 
     # Every scalar column not claimed as a quality metric becomes a sorting property (labels like
@@ -2333,13 +2343,10 @@ def read_nwb_sorting_analyzer(
         # the last stored spike (one chunk, ~4 MB) and double it. spike_times is stored per-unit-
         # concatenated (not globally time-sorted), so the last stored value is a lower bound on the true
         # last spike; doubling generously covers the gap while over-estimating only adds a harmless empty
-        # tail. Read from sorting_tmp because a group selection wraps `sorting` in a UnitsSelectionSorting
-        # whose segment has no spike_times_data.
+        # tail.
         if recording_duration is not None:
             placeholder_duration = recording_duration
         else:
-            spike_times_data = sorting_tmp._sorting_segments[0].spike_times_data
-            last_stored_spike_time = 0.0 if spike_times_data.shape[0] == 0 else float(np.asarray(spike_times_data[-1]))
             placeholder_duration = max(last_stored_spike_time * 2.0, 1.0)
             warnings.warn(
                 "No recording attached and no `recording_duration` given, so the recording duration is estimated. "
@@ -2397,14 +2404,8 @@ def read_nwb_sorting_analyzer(
 
     # templates (from the resolved templates column), and the required random_spikes extension
     if templates_column is not None and templates_column in units and unit_local_channels is not None:
-        # random_spikes only needs the total spike count of the (possibly group-selected) sorting. Read it
-        # from the NWB spike_times_index (per-unit cumulative counts, ~KB) on the unwrapped sorting_tmp,
-        # restricted to the selected units, rather than counting through the sorting object: a group
-        # selection wraps `sorting` in a UnitsSelectionSorting whose count materializes the parent's whole
-        # spike vector (~130-480 MB). np.isin keeps the units still present after any group selection.
-        per_unit_spike_counts = np.diff(sorting_tmp._sorting_segments[0].spike_times_index_data[:], prepend=0)
-        selected_units_mask = np.isin(sorting_tmp.unit_ids, sorting.unit_ids)
-        total_num_spikes = int(per_unit_spike_counts[selected_units_mask].sum())
+        # random_spikes only needs the total spike count of the selected units, read above
+        total_num_spikes = int(per_unit_spike_counts.sum())
         _make_random_spikes(analyzer, total_num_spikes)
         _make_templates(
             analyzer,

From 7a02e2f9c56702a7b4753ab90deeef302c960de3 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 15:44:51 -0600
Subject: [PATCH 40/46] Raise when a unit's electrodes are not channels of the
 attached recording

---
 .../extractors/nwbextractors.py               | 22 ++++++++++++++-----
 .../extractors/tests/test_nwbextractors.py    | 14 ++++++++++--
 2 files changed, 29 insertions(+), 7 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 6726d9c6fb..2ad0753b17 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2493,14 +2493,26 @@ def _make_sparsity_from_electrodes(sorting, electrodes_table, electrodes_indices
         electrode_row_to_channel_id = electrodes_table.index.to_numpy()
     position_of_channel_id = {channel_id: pos for pos, channel_id in enumerate(analyzer_channel_ids)}
 
+    # Every electrode of a unit has to be an analyzer channel. Only an attached recording can miss one (without it
+    # the channels are built from the units' electrodes), and then the recording does not match the sorting.
+    missing_electrodes = {}
+    for unit_id, region in zip(sorting.unit_ids, electrodes_indices):
+        missing = [electrode_row_to_channel_id[int(row)] for row in region]
+        missing = [channel_id for channel_id in missing if channel_id not in position_of_channel_id]
+        if len(missing) > 0:
+            missing_electrodes[unit_id] = missing
+    if len(missing_electrodes) > 0:
+        examples = ", ".join(f"{unit_id}: {missing}" for unit_id, missing in list(missing_electrodes.items())[:5])
+        raise ValueError(
+            f"{len(missing_electrodes)} units have electrodes that are not channels of the recording ({examples}), "
+            "so the recording does not match the sorting. Check that `electrical_series_path` is the series the "
+            "sorting was run on, or pass `electrical_series_path=None`."
+        )
+
     unit_local_channels = []
     sparsity_mask = np.zeros((sorting.get_num_units(), num_channels), dtype=bool)
     for unit_index, region in enumerate(electrodes_indices):
-        positions = [
-            position_of_channel_id[electrode_row_to_channel_id[int(electrode_row)]]
-            for electrode_row in region
-            if electrode_row_to_channel_id[int(electrode_row)] in position_of_channel_id
-        ]
+        positions = [position_of_channel_id[electrode_row_to_channel_id[int(row)]] for row in region]
         unit_local_channels.append(positions)
         sparsity_mask[unit_index, positions] = True
     sparsity = ChannelSparsity(
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index 7f30274718..f1721e54ac 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1276,6 +1276,7 @@ def _make_units_nwb(
     with_per_spike=False,
     with_electrical_series=False,
     single_channel=False,
+    series_electrodes=None,
 ):
     """Build a minimal NWB file with a Units table for read_nwb_sorting_analyzer tests.
 
@@ -1283,6 +1284,7 @@ def _make_units_nwb(
     quality metric (`snr`), a non-canonical numeric column (`custom_score`), a string label
     (`ks_label`), and a per-unit `electrodes` region; the electrodes table carries `rel_x`/`rel_y`.
     With `single_channel`, each unit has one electrode (`unit % n_ch`) and a 1D `waveform_mean`.
+    `series_electrodes` are the electrodes of the ElectricalSeries (all of them by default).
     """
     from datetime import datetime
     from pynwb import NWBFile, NWBHDF5IO
@@ -1297,10 +1299,11 @@ def _make_units_nwb(
     if with_electrical_series:
         from pynwb.ecephys import ElectricalSeries
 
-        region = nwbfile.create_electrode_table_region(region=list(range(n_ch)), description="all electrodes")
+        series_electrodes = list(range(n_ch)) if series_electrodes is None else list(series_electrodes)
+        region = nwbfile.create_electrode_table_region(region=series_electrodes, description="series electrodes")
         series = ElectricalSeries(
             name="ElectricalSeries",
-            data=np.zeros((1_000, n_ch), dtype="int16"),
+            data=np.zeros((1_000, len(series_electrodes)), dtype="int16"),
             electrodes=region,
             rate=30_000.0,
             starting_time=5.0,
@@ -1611,6 +1614,13 @@ def test_read_nwb_sorting_analyzer_electrical_series_path(tmp_path):
     )
     assert not analyzer.has_recording()
 
+    # a series that misses some of the units' electrodes does not match the sorting
+    path = _make_units_nwb(
+        tmp_path / "partial_series.nwb", n_ch=8, with_electrical_series=True, series_electrodes=range(6)
+    )
+    with pytest.raises(ValueError, match="not channels of the recording"):
+        read_nwb_sorting_analyzer(path, electrical_series_path="acquisition/ElectricalSeries", compute_extra=None)
+
 
 if __name__ == "__main__":
     tmp_path = Path("tmp")

From 6a3838dda42df2fd7cbf4b8af120c8e990cd1f57 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 16:03:15 -0600
Subject: [PATCH 41/46] List the electrode groups when group_name is missing or
 matches no units

---
 .../extractors/nwbextractors.py               | 30 +++++++++++--------
 .../extractors/tests/test_nwbextractors.py    |  9 ++++++
 2 files changed, 26 insertions(+), 13 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 2ad0753b17..409d7ce862 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2305,19 +2305,23 @@ def read_nwb_sorting_analyzer(
         # here we assume all groups are the same for each unit, so we just check one.
         if "group_name" in electrodes_table.columns:
             group_names = np.array([electrodes_table.iloc[int(ei[0])]["group_name"] for ei in electrodes_indices])
-            if len(np.unique(group_names)) > 1:
-                if group_name is None:
-                    raise Exception(
-                        f"More than one group, use group_name option to select units. Available groups: {np.unique(group_names)}"
-                    )
-                else:
-                    unit_mask = group_names == group_name
-                    if verbose:
-                        print(f"Selecting {sum(unit_mask)} / {len(units)} units from {group_name}")
-                    sorting = sorting.select_units(unit_ids=sorting.unit_ids[unit_mask])
-                    units = units.loc[units.index[unit_mask]]
-                    per_unit_spike_counts = per_unit_spike_counts[unit_mask]
-                    electrodes_indices = [region for region, keep in zip(electrodes_indices, unit_mask) if keep]
+            unique_group_names = np.unique(group_names)
+            available_groups = ", ".join(repr(str(name)) for name in unique_group_names)
+            if group_name is None and len(unique_group_names) > 1:
+                raise ValueError(
+                    f"The units span {len(unique_group_names)} electrode groups. Pass `group_name` to load one: "
+                    f"{available_groups}"
+                )
+            if group_name is not None and group_name not in unique_group_names:
+                raise ValueError(f"No units in group {group_name!r}. The units' groups are: {available_groups}")
+            if group_name is not None and len(unique_group_names) > 1:
+                unit_mask = group_names == group_name
+                if verbose:
+                    print(f"Selecting {sum(unit_mask)} / {len(units)} units from {group_name}")
+                sorting = sorting.select_units(unit_ids=sorting.unit_ids[unit_mask])
+                units = units.loc[units.index[unit_mask]]
+                per_unit_spike_counts = per_unit_spike_counts[unit_mask]
+                electrodes_indices = [region for region, keep in zip(electrodes_indices, unit_mask) if keep]
 
     # Every scalar column not claimed as a quality metric becomes a sorting property (labels like
     # cluster_uuid, and any non-canonical numeric column). Properties show up in the GUI's unit table
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index f1721e54ac..5bc54acae6 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1560,6 +1560,15 @@ def test_read_nwb_sorting_analyzer_peak_channel_group_mapping(tmp_path):
         ],
         "quality_metrics": [{"source": "units_columns", "columns": {"isi_violations_ratio": "isi_violations"}}],
     }
+    # the units span two groups, so one has to be named, and it has to exist
+    read_kwargs = dict(electrical_series_path=None, sampling_frequency=30000.0, recording_duration=100.0)
+    with pytest.raises(ValueError, match="span 2 electrode groups"):
+        read_nwb_sorting_analyzer(str(path), nwb_to_analyzer_mapping=allen_mapping, compute_extra=None, **read_kwargs)
+    with pytest.raises(ValueError, match="No units in group 'probeX'"):
+        read_nwb_sorting_analyzer(
+            str(path), group_name="probeX", nwb_to_analyzer_mapping=allen_mapping, compute_extra=None, **read_kwargs
+        )
+
     analyzer = read_nwb_sorting_analyzer(
         str(path),
         electrical_series_path=None,

From d39b99c0c384456c09c32181e1cd40bb709f4075 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 16:06:29 -0600
Subject: [PATCH 42/46] Remove the typed_container places from the default
 mapping until ndx-spikesorting is read natively

---
 src/spikeinterface/extractors/nwbextractors.py | 12 +-----------
 1 file changed, 1 insertion(+), 11 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index 409d7ce862..c6bd2f1fcf 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1933,8 +1933,6 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
 #   * "electrodes_columns": columns of the electrodes table
 #   * "peak_channel_group": a Units column holding each unit's peak electrode id; the unit covers every
 #     electrode in that electrode's group, ordered by the electrodes column "order_by"
-#   * "typed_container": an ndx-spikesorting object. This is a hook: the typed reader lives in
-#     ndx-spikesorting for now, so it never resolves here and the next place is used.
 DEFAULT_NWB_TO_ANALYZER_MAPPING = {
     # which electrodes each unit's waveform covers (sparsity, and where the templates are placed)
     "sparsity": [
@@ -1945,21 +1943,13 @@ def read_nwb(file_path, load_recording=True, load_sorting=False, electrical_seri
         {"source": "electrodes_columns", "x": "rel_x", "y": "rel_y"},
     ],
     "templates": [
-        {"source": "typed_container", "type": "Templates"},
         {"source": "units_columns", "mean": "waveform_mean", "std": "waveform_sd"},
     ],
     # "columns" is "canonical_quality" (the columns named like a SpikeInterface quality metric), a list of
     # column names, or a dict from SpikeInterface metric name to column name
     "quality_metrics": [
-        {"source": "typed_container", "type": "UnitsMetrics"},
         {"source": "units_columns", "columns": "canonical_quality"},
     ],
-    "template_metrics": [
-        {"source": "typed_container", "type": "UnitsMetrics"},
-    ],
-    "unit_locations": [
-        {"source": "typed_container", "type": "UnitLocations"},
-    ],
 }
 
 
@@ -1979,7 +1969,7 @@ def _resolve_nwb_to_analyzer_mapping(mapping, scalar_colnames, all_colnames):
 
     Returns which place gives each unit's electrodes, which columns hold the templates, which columns are
     quality metrics (as a dict from metric name to column) and which scalar columns become sorting
-    properties. "typed_container" places are skipped (see `DEFAULT_NWB_TO_ANALYZER_MAPPING`).
+    properties.
     `scalar_colnames` are the 1D per-unit columns eligible to be metrics or properties, `all_colnames` is
     every column name.
     """

From 2321e7699e5a71b44a39079e156bcc5af9db5e2a Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 16:16:11 -0600
Subject: [PATCH 43/46] Use the NWB 2.11 waveform unit and
 time_before_peak_in_ms attributes and test template values on both read paths

---
 .../extractors/nwbextractors.py               | 75 +++++++++++++------
 .../extractors/tests/test_nwbextractors.py    | 74 +++++++++++++++---
 2 files changed, 117 insertions(+), 32 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index c6bd2f1fcf..e1f8a9d459 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2122,7 +2122,8 @@ def read_nwb_sorting_analyzer(
         "templates", "quality_metrics", ...) and maps to a list of places tried in order. Setting a key to
         None disables it.
     rescale_templates_to_uV : bool, default: True
-        Infer the scale of `waveform_mean` from its magnitude and convert the templates to microvolts.
+        Convert the templates to microvolts. The scale comes from the `unit` attribute of `waveform_mean` when
+        it is set to something other than the default "volts", and is inferred from the magnitude otherwise.
     verbose : bool, default: False
         If True, print progress messages.
 
@@ -2268,6 +2269,11 @@ def read_nwb_sorting_analyzer(
     if unit_electrodes is not None and unit_electrodes["source"] == "units_columns":
         electrodes_column = unit_electrodes["electrodes"]
 
+    # the unit and alignment attributes of the templates column (read before a group selection replaces `sorting`)
+    templates_attributes = {}
+    if templates_column is not None:
+        templates_attributes = _column_attributes(sorting.units_table, templates_column)
+
     if not use_pynwb:
         # deliberate read: only the resolved template column(s), the electrodes region, and the scalar
         # columns that become metrics or properties are materialized; ragged/per-spike columns are not.
@@ -2409,6 +2415,8 @@ def read_nwb_sorting_analyzer(
             unit_local_channels,
             num_channels,
             rescale_to_uV=rescale_templates_to_uV,
+            waveform_unit=templates_attributes.get("unit"),
+            time_before_peak_ms=templates_attributes.get("time_before_peak_in_ms"),
         )
 
     # the resolved quality-metric columns -> quality_metrics extension
@@ -2533,14 +2541,28 @@ def _make_random_spikes(analyzer, total_num_spikes):
     analyzer.extensions["random_spikes"] = random_spikes_ext
 
 
+# factor from a waveform `unit` attribute to microvolts
+_WAVEFORM_UNIT_TO_UV = {"volts": 1e6, "V": 1e6, "millivolts": 1e3, "mV": 1e3, "microvolts": 1.0, "uV": 1.0, "µV": 1.0}
+
+
+def _column_attributes(units_table, column):
+    """The attributes of a Units column, from a pynwb table or an h5py/zarr group."""
+    column_object = units_table[column]
+    column_object = getattr(column_object, "target", column_object)  # a ragged pynwb column is its VectorIndex
+    dataset = getattr(column_object, "data", column_object)
+    attributes = {}
+    for key, value in getattr(dataset, "attrs", {}).items():
+        attributes[key] = value.decode("utf-8") if isinstance(value, bytes) else value
+    return attributes
+
+
 def _infer_template_scale_to_uV(dense_templates):
     """Infer the factor that brings templates into microvolts, SpikeInterface's convention.
 
-    The NWB waveform unit attribute is schema-fixed to "volts" and never reflects the true scale
-    (pynwb issue #2162), so it cannot drive the conversion. The template magnitude can: an extracellular
-    spike peaks at roughly tens to hundreds of microvolts, so pick the factor in {1, 1e3, 1e6} that lands
-    the median per-unit peak in that physiological window. This reads no extra data (it operates on the
-    already-built templates). Returns 1.0 if the magnitude is unusable (all zero).
+    Used when the waveform `unit` attribute cannot be trusted (see `_make_templates`). An extracellular spike
+    peaks at roughly tens to hundreds of microvolts, so pick the factor in {1, 1e3, 1e6} that lands the median
+    per-unit peak in that physiological window. This reads no extra data (it operates on the already-built
+    templates). Returns 1.0 if the magnitude is unusable (all zero).
     """
     per_unit_peak = np.abs(dense_templates).max(axis=(1, 2))
     per_unit_peak = per_unit_peak[per_unit_peak > 0]
@@ -2561,9 +2583,12 @@ def _make_templates(
     unit_local_channels,
     num_channels,
     rescale_to_uV=False,
+    waveform_unit=None,
+    time_before_peak_ms=None,
 ):
     """Attach the `templates` extension from the resolved templates column, scattering each unit's
-    sparse waveform block onto its true channel positions."""
+    sparse waveform block onto its true channel positions. `waveform_unit` and `time_before_peak_ms` are the
+    `unit` and `time_before_peak_in_ms` attributes of the templates column, when it has them."""
 
     waveform_mean = np.array([np.asarray(t, dtype="float") for t in units[templates_column].values])
     # The schema also allows one waveform per unit with no electrode axis (e.g. DANDI 000363). It belongs to
@@ -2588,29 +2613,35 @@ def _make_templates(
         waveform_mean = waveform_mean.transpose(0, 2, 1)
     num_samples_template = waveform_mean.shape[1]
 
-    # Scatter each unit's waveform onto its channels. This relies on the alignment contract that
-    # waveform_mean[:, i] corresponds to the unit's electrodes-region entry i. The IBL writer guarantees
-    # it by ordering each unit's real channels first (padding last) and building the electrodes region
-    # from exactly those real channels in the same order, so the first k = len(region) columns are the
-    # real channels in region order and any trailing padding columns are correctly dropped by [:, :k].
+    # Scatter each unit's waveform onto its channels: column i is the unit's electrodes-region entry i. In
+    # the schema's single waveform array every unit has the same number of columns, so a writer whose units
+    # cover different numbers of electrodes pads the shorter ones at the end (e.g. the IBL files). The first
+    # k = len(region) columns are the unit's electrodes and the padding after them is dropped by [:, :k].
     dense_templates = np.zeros((analyzer.get_num_units(), num_samples_template, num_channels), dtype="float32")
     for unit_index, positions in enumerate(unit_local_channels):
         k = len(positions)
         dense_templates[unit_index][:, positions] = waveform_mean[unit_index][:, :k]
 
-    # The stored waveform unit attribute is schema-fixed to "volts" and never reflects the real scale
-    # (pynwb #2162), so it cannot drive a conversion. When rescaling is requested, infer the factor from
-    # the template magnitude instead and apply the same factor to the std templates below so both stay in
-    # the same unit. On by default (read_nwb_sorting_analyzer) so templates read in microvolts; it can be
-    # turned off to keep the file's native unit, since the scale does not affect peak channel, shape,
-    # location, or sparsity.
-    template_scale = _infer_template_scale_to_uV(dense_templates) if rescale_to_uV else 1.0
+    # Convert to microvolts when requested, applying the same factor to the std templates below. Before NWB
+    # 2.11 the waveform `unit` attribute was fixed to "volts" whatever the real scale, and since then it can
+    # be set but still defaults to "volts". So a unit other than "volts" is trusted, and otherwise the factor
+    # is inferred from the template magnitude. It can be turned off to keep the file's native unit, since the
+    # scale does not affect peak channel, shape, location, or sparsity.
+    template_scale = 1.0
+    if rescale_to_uV:
+        template_scale = None if waveform_unit == "volts" else _WAVEFORM_UNIT_TO_UV.get(waveform_unit)
+        if template_scale is None:
+            template_scale = _infer_template_scale_to_uV(dense_templates)
     if template_scale != 1.0:
         dense_templates *= template_scale
 
-    # scalar alignment (nbefore): the sample where the average peak-channel amplitude is largest
-    peak_amplitude_per_sample = np.abs(np.nan_to_num(waveform_mean)).max(axis=2).mean(axis=0)
-    nbefore = int(np.argmax(peak_amplitude_per_sample))
+    # scalar alignment (nbefore): the `time_before_peak_in_ms` attribute (NWB 2.11) when the file has it,
+    # otherwise the sample where the average peak-channel amplitude is largest
+    if time_before_peak_ms is not None:
+        nbefore = int(round(float(time_before_peak_ms) * analyzer.sampling_frequency / 1000))
+    else:
+        peak_amplitude_per_sample = np.abs(np.nan_to_num(waveform_mean)).max(axis=2).mean(axis=0)
+        nbefore = int(np.argmax(peak_amplitude_per_sample))
 
     templates_ext = ComputeTemplates(sorting_analyzer=analyzer)
     templates_ext.data["average"] = dense_templates
diff --git a/src/spikeinterface/extractors/tests/test_nwbextractors.py b/src/spikeinterface/extractors/tests/test_nwbextractors.py
index 5bc54acae6..a1e7804042 100644
--- a/src/spikeinterface/extractors/tests/test_nwbextractors.py
+++ b/src/spikeinterface/extractors/tests/test_nwbextractors.py
@@ -1277,6 +1277,10 @@ def _make_units_nwb(
     with_electrical_series=False,
     single_channel=False,
     series_electrodes=None,
+    n_unit_electrodes=None,
+    padding=0,
+    waveform_scale=1.0,
+    waveform_attributes=None,
 ):
     """Build a minimal NWB file with a Units table for read_nwb_sorting_analyzer tests.
 
@@ -1285,6 +1289,9 @@ def _make_units_nwb(
     (`ks_label`), and a per-unit `electrodes` region; the electrodes table carries `rel_x`/`rel_y`.
     With `single_channel`, each unit has one electrode (`unit % n_ch`) and a 1D `waveform_mean`.
     `series_electrodes` are the electrodes of the ElectricalSeries (all of them by default).
+    With `n_unit_electrodes`, unit u covers electrodes `(u + j) % n_ch` for j < n_unit_electrodes, and every
+    waveform has `padding` extra columns after them, as in the IBL files. `waveform_scale` multiplies the
+    waveforms and `waveform_attributes` are written on the waveform datasets.
     """
     from datetime import datetime
     from pynwb import NWBFile, NWBHDF5IO
@@ -1315,33 +1322,55 @@ def _make_units_nwb(
     if with_per_spike:
         nwbfile.add_unit_column(name="spike_amplitudes_uV", description="per-spike amplitudes", index=True)
     rng = np.random.default_rng(0)
-    waveform_shape = (n_samp,) if single_channel else (n_samp, n_ch)
+    if single_channel:
+        waveform_shape = (n_samp,)
+    elif n_unit_electrodes is not None:
+        waveform_shape = (n_samp, n_unit_electrodes + padding)
+    else:
+        waveform_shape = (n_samp, n_ch)
     for u in range(n_units):
+        if single_channel:
+            electrodes = [u % n_ch]
+        elif n_unit_electrodes is not None:
+            electrodes = [(u + j) % n_ch for j in range(n_unit_electrodes)]
+        else:
+            electrodes = list(range(n_ch))
         kwargs = dict(
             spike_times=np.sort(rng.uniform(0, 100, size=int(rng.integers(50, 500)))),
-            electrodes=[u % n_ch] if single_channel else list(range(n_ch)),
-            waveform_mean=rng.standard_normal(waveform_shape),
+            electrodes=electrodes,
+            waveform_mean=rng.standard_normal(waveform_shape) * waveform_scale,
             snr=float(rng.uniform(2, 10)),
             custom_score=float(rng.uniform(0, 1)),
             ks_label="good" if u % 2 == 0 else "mua",
         )
         if with_std:
-            kwargs["waveform_sd"] = np.abs(rng.standard_normal(waveform_shape))
+            kwargs["waveform_sd"] = np.abs(rng.standard_normal(waveform_shape)) * waveform_scale
         if with_per_spike:
             kwargs["spike_amplitudes_uV"] = rng.standard_normal(kwargs["spike_times"].size)
         nwbfile.add_unit(**kwargs)
     with NWBHDF5IO(str(path), "w") as io:
         io.write(nwbfile)
+    if waveform_attributes is not None:
+        import h5py
+
+        with h5py.File(path, "r+") as file:
+            for column in ("waveform_mean", "waveform_sd"):
+                if column in file["units"]:
+                    file["units"][column].attrs.update(waveform_attributes)
     return str(path)
 
 
-def test_read_nwb_sorting_analyzer_default(tmp_path):
+@pytest.mark.parametrize("use_pynwb", [True, False])
+def test_read_nwb_sorting_analyzer_default(tmp_path, use_pynwb):
     pytest.importorskip("pynwb")
+    import h5py
     from spikeinterface.extractors import read_nwb_sorting_analyzer
 
-    path = _make_units_nwb(tmp_path / "units.nwb")
+    # waveforms in volts, each unit on 5 of the 8 electrodes with 3 padding columns, as in the IBL files
+    n_unit_electrodes = 5
+    path = _make_units_nwb(tmp_path / "units.nwb", n_unit_electrodes=n_unit_electrodes, padding=3, waveform_scale=1e-4)
     analyzer = read_nwb_sorting_analyzer(
-        path, electrical_series_path=None, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None
+        path, electrical_series_path=None, use_pynwb=use_pynwb, sampling_frequency=30000.0, compute_extra=None
     )
 
     # recordingless, sparse, with templates for every unit
@@ -1350,6 +1379,18 @@ def test_read_nwb_sorting_analyzer_default(tmp_path):
     assert analyzer.has_extension("templates")
     assert analyzer.get_extension("templates").get_data().shape[0] == analyzer.get_num_units()
 
+    # each waveform column is on its electrode, in microvolts, and the padding is dropped
+    with h5py.File(path, "r") as file:
+        stored_mean = file["units"]["waveform_mean"][:]
+    average = analyzer.get_extension("templates").data["average"]
+    position_of_electrode = {int(channel_id): pos for pos, channel_id in enumerate(analyzer.channel_ids)}
+    for unit_index in range(analyzer.get_num_units()):
+        electrodes = [(unit_index + j) % 8 for j in range(n_unit_electrodes)]
+        positions = [position_of_electrode[electrode] for electrode in electrodes]
+        expected = stored_mean[unit_index][:, :n_unit_electrodes] * 1e6
+        np.testing.assert_allclose(average[unit_index][:, positions], expected, rtol=1e-5)
+        assert not np.any(np.delete(average[unit_index], positions, axis=1))
+
     # by default only canonical quality columns become quality_metrics; the rest become properties, and
     # template_metrics is never synthesized from value columns
     assert list(analyzer.get_extension("quality_metrics").get_data().columns) == ["snr"]
@@ -1395,19 +1436,32 @@ def test_read_nwb_sorting_analyzer_mapping_disable(tmp_path):
     assert {"snr", "custom_score", "ks_label"} <= set(analyzer.sorting.get_property_keys())
 
 
-def test_read_nwb_sorting_analyzer_waveform_sd(tmp_path):
+@pytest.mark.parametrize("use_pynwb", [True, False])
+def test_read_nwb_sorting_analyzer_waveform_sd(tmp_path, use_pynwb):
     pytest.importorskip("pynwb")
+    import h5py
     from spikeinterface.extractors import read_nwb_sorting_analyzer
 
-    path = _make_units_nwb(tmp_path / "units.nwb", with_std=True)
+    # waveforms in millivolts, declared with the NWB 2.11 unit and alignment attributes. The magnitude alone
+    # (peaks of a few units) would not be rescaled, so the conversion has to come from the unit attribute.
+    path = _make_units_nwb(
+        tmp_path / "units.nwb", with_std=True, waveform_attributes={"unit": "millivolts", "time_before_peak_in_ms": 0.5}
+    )
     analyzer = read_nwb_sorting_analyzer(
-        path, electrical_series_path=None, use_pynwb=False, sampling_frequency=30000.0, compute_extra=None
+        path, electrical_series_path=None, use_pynwb=use_pynwb, sampling_frequency=30000.0, compute_extra=None
     )
     templates = analyzer.get_extension("templates")
     # the std operator is populated only because the file stores waveform_sd
     assert "std" in templates.params["operators"]
     assert "std" in templates.data
 
+    with h5py.File(path, "r") as file:
+        stored_mean = file["units"]["waveform_mean"][:]
+        stored_sd = file["units"]["waveform_sd"][:]
+    np.testing.assert_allclose(templates.data["average"], stored_mean * 1e3, rtol=1e-5)
+    np.testing.assert_allclose(templates.data["std"], stored_sd * 1e3, rtol=1e-5)
+    assert templates.params["ms_before"] == pytest.approx(0.5)
+
 
 @pytest.mark.parametrize("use_pynwb", [True, False])
 def test_read_nwb_sorting_analyzer_skips_per_spike_columns(tmp_path, use_pynwb, monkeypatch):

From 268ad19e925db75d10cb3bcd4cb3d19f31b73884 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 16:40:43 -0600
Subject: [PATCH 44/46] Say that t_start is left at 0 only when it is not
 passed

---
 src/spikeinterface/extractors/nwbextractors.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index e1f8a9d459..f247028b28 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -2211,7 +2211,7 @@ def read_nwb_sorting_analyzer(
     if recording is None and recording_duration is None:
         spike_times_data = sorting_segment.spike_times_data
         last_stored_spike_time = 0.0 if spike_times_data.shape[0] == 0 else float(np.asarray(spike_times_data[-1]))
-    # Recordingless case: leave t_start at 0 (set when the sorting was constructed). NWB spike times are
+    # Recordingless case without a `t_start`: leave it at 0 (set when the sorting was constructed). NWB spike times are
     # in seconds from session start, so they are always >= 0 and map to non-negative frames with
     # t_start = 0. Anchoring t_start to the first spike would only remove empty space before it on the
     # timeline (cosmetic) and would cost a streamed read, so we skip it to keep the build cheap.

From 5278aa016c3732a82ac27a18999750210cad5a12 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 12 Oct 2026 18:13:51 -0600
Subject: [PATCH 45/46] Read the waveforms only for the selected units in
 read_nwb_sorting_analyzer

---
 .../extractors/nwbextractors.py               | 96 +++++++++++++++----
 1 file changed, 78 insertions(+), 18 deletions(-)

diff --git a/src/spikeinterface/extractors/nwbextractors.py b/src/spikeinterface/extractors/nwbextractors.py
index f247028b28..fcceaae600 100644
--- a/src/spikeinterface/extractors/nwbextractors.py
+++ b/src/spikeinterface/extractors/nwbextractors.py
@@ -1964,6 +1964,49 @@ def _mapping_structural_colnames(mapping):
     return structural
 
 
+def _mapping_waveform_colnames(mapping):
+    """The Units columns that the mapping reads as waveforms. They are read apart from the rest of the table,
+    only for the selected units (see `_read_unit_rows`)."""
+    waveform_colnames = set()
+    for place in mapping.get("templates") or []:
+        if place.get("source") == "units_columns":
+            waveform_colnames.update(place.get(role) for role in ("mean", "std") if place.get(role))
+    return waveform_colnames
+
+
+def _read_unit_rows(units_table, column, unit_rows):
+    """Read the values of a Units column for the units at the (sorted) positions `unit_rows` only.
+
+    Works on a pynwb table and on an h5py/zarr group, for regular and ragged columns, and returns one array per
+    unit. It reads the range of rows from the first to the last selected unit in one request: over a stream a
+    few large reads are much faster than many small ones, and the selected units are usually stored together
+    (e.g. the units of one probe). At worst the range is the whole column, which is what reading it all costs.
+    """
+    column_object = units_table[column]
+    if hasattr(units_table, "colnames"):  # pynwb: a ragged column is its VectorIndex
+        ragged = hasattr(column_object, "target")
+        data = column_object.target.data if ragged else column_object.data
+        index = column_object.data if ragged else None
+    else:  # h5py or zarr group
+        data = column_object
+        index = units_table[f"{column}_index"] if f"{column}_index" in units_table else None
+
+    unit_rows = np.asarray(unit_rows)
+    if unit_rows.size == 0:
+        return []
+    first, last = int(unit_rows[0]), int(unit_rows[-1])
+    if index is None:
+        values = np.asarray(data[first : last + 1])
+        return [values[row - first] for row in unit_rows]
+
+    # ragged: unit u holds the data rows stops[u - 1]:stops[u]
+    stops = np.asarray(index[: last + 1], dtype="int64")
+    starts = np.concatenate(([0], stops[:-1]))
+    span_start = int(starts[first])
+    values = np.asarray(data[span_start : int(stops[last])])
+    return [values[starts[row] - span_start : stops[row] - span_start] for row in unit_rows]
+
+
 def _resolve_nwb_to_analyzer_mapping(mapping, scalar_colnames, all_colnames):
     """Resolve the Units-table places of the mapping against the columns of the Units table.
 
@@ -2228,10 +2271,16 @@ def read_nwb_sorting_analyzer(
         colnames = list(units_table.colnames)
         # ragged columns (per-spike data) are excluded so to_dataframe never reads them; electrodes is
         # ragged too but is the per-unit region needed for sparsity
-        # the structural columns (electrodes region, waveforms) are kept even when ragged
+        # the structural columns (electrodes region, waveforms) are kept even when ragged, and the waveforms
+        # are read later, only for the selected units
         structural = _mapping_structural_colnames(mapping)
-        ragged_colnames = {c for c in colnames if c not in structural and isinstance(units_table[c], VectorIndex)}
-        units = units_table.to_dataframe(index=True, exclude=ragged_colnames)
+        waveform_colnames = _mapping_waveform_colnames(mapping)
+        excluded_colnames = {
+            c
+            for c in colnames
+            if c in waveform_colnames or (c not in structural and isinstance(units_table[c], VectorIndex))
+        }
+        units = units_table.to_dataframe(index=True, exclude=excluded_colnames)
         metric_colnames = [c for c in units.columns if c not in structural and units[c].dtype.kind in "fiu"]
         label_colnames = [
             c
@@ -2270,17 +2319,22 @@ def read_nwb_sorting_analyzer(
         electrodes_column = unit_electrodes["electrodes"]
 
     # the unit and alignment attributes of the templates column (read before a group selection replaces `sorting`)
+    units_table_object = sorting.units_table
     templates_attributes = {}
     if templates_column is not None:
-        templates_attributes = _column_attributes(sorting.units_table, templates_column)
+        templates_attributes = _column_attributes(units_table_object, templates_column)
 
     if not use_pynwb:
-        # deliberate read: only the resolved template column(s), the electrodes region, and the scalar
-        # columns that become metrics or properties are materialized; ragged/per-spike columns are not.
-        needed = [c for c in (templates_column, templates_std_column, electrodes_column) if c and c in colnames]
+        # deliberate read: only the electrodes region and the scalar columns that become metrics or properties
+        # are materialized here; ragged/per-spike columns are not, and the waveforms are read later, only for
+        # the selected units
+        needed = [electrodes_column] if electrodes_column in colnames else []
         needed += scalar_colnames
         units = _create_df_from_nwb_table(units_group, columns=needed)
 
+    # positions in the Units table of the units that are loaded, narrowed by a group selection below
+    unit_rows = np.arange(len(units))
+
     if use_pynwb:
         electrodes_table = sorting._reader.nwbfile.electrodes.to_dataframe(index=True)
     else:
@@ -2316,6 +2370,7 @@ def read_nwb_sorting_analyzer(
                     print(f"Selecting {sum(unit_mask)} / {len(units)} units from {group_name}")
                 sorting = sorting.select_units(unit_ids=sorting.unit_ids[unit_mask])
                 units = units.loc[units.index[unit_mask]]
+                unit_rows = unit_rows[unit_mask]
                 per_unit_spike_counts = per_unit_spike_counts[unit_mask]
                 electrodes_indices = [region for region, keep in zip(electrodes_indices, unit_mask) if keep]
 
@@ -2403,15 +2458,19 @@ def read_nwb_sorting_analyzer(
     num_channels = len(analyzer_channel_ids)
 
     # templates (from the resolved templates column), and the required random_spikes extension
-    if templates_column is not None and templates_column in units and unit_local_channels is not None:
+    if templates_column is not None and unit_local_channels is not None:
         # random_spikes only needs the total spike count of the selected units, read above
         total_num_spikes = int(per_unit_spike_counts.sum())
         _make_random_spikes(analyzer, total_num_spikes)
+        waveform_mean_rows = _read_unit_rows(units_table_object, templates_column, unit_rows)
+        waveform_sd_rows = None
+        if templates_std_column is not None:
+            waveform_sd_rows = _read_unit_rows(units_table_object, templates_std_column, unit_rows)
         _make_templates(
             analyzer,
-            units,
             templates_column,
-            templates_std_column,
+            waveform_mean_rows,
+            waveform_sd_rows,
             unit_local_channels,
             num_channels,
             rescale_to_uV=rescale_templates_to_uV,
@@ -2577,20 +2636,21 @@ def _infer_template_scale_to_uV(dense_templates):
 
 def _make_templates(
     analyzer,
-    units,
     templates_column,
-    templates_std_column,
+    waveform_mean_rows,
+    waveform_sd_rows,
     unit_local_channels,
     num_channels,
     rescale_to_uV=False,
     waveform_unit=None,
     time_before_peak_ms=None,
 ):
-    """Attach the `templates` extension from the resolved templates column, scattering each unit's
-    sparse waveform block onto its true channel positions. `waveform_unit` and `time_before_peak_ms` are the
-    `unit` and `time_before_peak_in_ms` attributes of the templates column, when it has them."""
+    """Attach the `templates` extension from the waveforms of the resolved templates column (and its std
+    column, when there is one), one array per unit, scattering each unit's sparse waveform block onto its true
+    channel positions. `waveform_unit` and `time_before_peak_ms` are the `unit` and `time_before_peak_in_ms`
+    attributes of the templates column, when it has them."""
 
-    waveform_mean = np.array([np.asarray(t, dtype="float") for t in units[templates_column].values])
+    waveform_mean = np.array([np.asarray(t, dtype="float") for t in waveform_mean_rows])
     # The schema also allows one waveform per unit with no electrode axis (e.g. DANDI 000363). It belongs to
     # the unit's electrode only when the region has exactly one, otherwise we cannot tell which one it is.
     single_channel = waveform_mean.ndim == 2
@@ -2650,8 +2710,8 @@ def _make_templates(
     # Only expose the "std" operator when the file actually stores a std column (resolved). Fabricating a
     # zero std would falsely imply the unit's spikes have no waveform variability, so when it is absent we
     # declare only "average" rather than inventing data.
-    if templates_std_column is not None and templates_std_column in units.columns:
-        waveform_sd = np.array([np.asarray(t, dtype="float") for t in units[templates_std_column].values])
+    if waveform_sd_rows is not None:
+        waveform_sd = np.array([np.asarray(t, dtype="float") for t in waveform_sd_rows])
         if single_channel:
             waveform_sd = waveform_sd[:, :, np.newaxis]
         if electrodes_by_samples:

From bf29f5cd35d64a1cabfa6eb5e1b8b08875fb4034 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 12 Nov 2026 09:22:46 +0000
Subject: [PATCH 46/46] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 .../generation/session_displacement_generator.py         | 9 +++------
 .../tests/test_session_displacement_generator.py         | 6 ------
 2 files changed, 3 insertions(+), 12 deletions(-)

diff --git a/src/spikeinterface/generation/session_displacement_generator.py b/src/spikeinterface/generation/session_displacement_generator.py
index ce81c2d4f9..d6421b9886 100644
--- a/src/spikeinterface/generation/session_displacement_generator.py
+++ b/src/spikeinterface/generation/session_displacement_generator.py
@@ -196,7 +196,6 @@ def generate_session_displacement_recordings(
     amplitude_scaling_order = None
 
     for rec_idx, (shift, duration) in enumerate(zip(recording_shifts, recording_durations)):
-
         displacement_vector, displacement_unit_factor = _get_inter_session_displacements(
             shift,
             non_rigid_gradient[rec_idx] if isinstance(non_rigid_gradient, list) else non_rigid_gradient,
@@ -238,7 +237,6 @@ def generate_session_displacement_recordings(
         )
 
         if recording_amplitude_scalings is not None:
-
             if rec_idx == 0:
                 amplitude_scaling_order = _get_amplitude_scaling_order(
                     templates_array_moved,
@@ -399,13 +397,12 @@ def _check_generate_session_displacement_arguments(
         raise ValueError("Each record entry for `recording_shifts` must have two elements, the x and y shift.")
 
     if recording_amplitude_scalings is not None:
-
         keys = recording_amplitude_scalings.keys()
-        if not "method" in keys or not "scalings" in keys:
+        if "method" not in keys or "scalings" not in keys:
             raise ValueError("`recording_amplitude_scalings` must be a dict with keys `method` and `scalings`.")
 
         allowed_methods = ["by_passed_order", "by_amplitude_and_firing_rate", "by_firing_rate"]
-        if not recording_amplitude_scalings["method"] in allowed_methods:
+        if recording_amplitude_scalings["method"] not in allowed_methods:
             raise ValueError(f"`recording_amplitude_scalings` must be one of {allowed_methods}")
 
         rec_scalings = recording_amplitude_scalings["scalings"]
@@ -420,7 +417,7 @@ def _check_generate_session_displacement_arguments(
 
     if isinstance(non_rigid_gradient, list):
         if not len(non_rigid_gradient) == expected_num_recs:
-            raise ValueError("If `non_rigid_gradient` is a list, it must " "contain one value for each recording.")
+            raise ValueError("If `non_rigid_gradient` is a list, it must contain one value for each recording.")
 
 
 def _update_kwargs_for_extended_units(
diff --git a/src/spikeinterface/generation/tests/test_session_displacement_generator.py b/src/spikeinterface/generation/tests/test_session_displacement_generator.py
index 05d6db4c8e..b1d8cffbfa 100644
--- a/src/spikeinterface/generation/tests/test_session_displacement_generator.py
+++ b/src/spikeinterface/generation/tests/test_session_displacement_generator.py
@@ -76,7 +76,6 @@ def test_x_y_rigid_shifts_are_properly_set(self, options):
         locations_1 = extra_outputs["unit_locations"][0]
 
         for rec_idx in range(1, 4):
-
             shifts = recording_shifts[rec_idx]
 
             assert np.array_equal(
@@ -91,14 +90,12 @@ def test_x_y_rigid_shifts_are_properly_set(self, options):
         min_channel_loc = output_recordings[0].get_channel_locations()[0, 1]
         max_channel_loc = output_recordings[0].get_channel_locations()[-1, 1]
         for unit_idx in range(num_units):
-
             start_pos = self._get_peak_chan_loc_in_um(
                 extra_outputs["templates_array_moved"][0][unit_idx],
                 options["y_bin_um"],
             )
 
             for rec_idx in range(1, options["num_recs"]):
-
                 new_pos = self._get_peak_chan_loc_in_um(
                     extra_outputs["templates_array_moved"][rec_idx][unit_idx], options["y_bin_um"]
                 )
@@ -156,7 +153,6 @@ def test_spike_times_and_firing_rates_across_recordings(self, options):
 
         for unit_idx in range(options["kwargs"]["num_units"]):
             for rec_idx in range(1, options["num_recs"]):
-
                 # Exact spike times are not preserved when seed is None
                 assert np.array_equal(
                     output_sortings_same[0].get_unit_spike_train(str(unit_idx)),
@@ -222,7 +218,6 @@ def test_x_y_shift_non_rigid(self, options, dim_idx):
         # w.r.t the first recording), check the rigid and nonrigid shifts
         # are as expected.
         for rec_idx in range(1, options["num_recs"]):
-
             shift = options["kwargs"]["recording_shifts"][rec_idx][dim_idx]
 
             # Get the rigid shift between the first recording and this shifted recording
@@ -437,7 +432,6 @@ def test_shift_units_outside_probe(self, options):
         num_extended_units = num_units * 3
 
         for ses_idx in range(num_sessions):
-
             # There are 3x the number of units when new units are created
             # (one new set above, and one new set below the probe).
             for key in ["unit_locations", "templates_array_moved", "firing_rates"]: