Support a SortingAnalyzer with no unit in the core extensions

src/spikeinterface/core/sorting_tools.py

@@ -250,7 +250,11 @@ def random_spikes_selection(
 
             random_spikes_indices.append(selected_unit_indices)
 
-        random_spikes_indices = np.concatenate(random_spikes_indices)
+        if len(random_spikes_indices) > 0:
+            random_spikes_indices = np.concatenate(random_spikes_indices)
+        else:
+            # a sorting with no unit is valid, np.concatenate would raise on the empty list
+            random_spikes_indices = np.zeros(0, dtype="int64")
         random_spikes_indices = np.sort(random_spikes_indices)
 
     else:

src/spikeinterface/core/tests/test_sorting_tools.py

@@ -83,6 +83,25 @@ def test_random_spikes_selection():
     assert random_spikes_indices.size == spikes.size
 
 
+@pytest.mark.parametrize("method", ["uniform", "percentage", "maximum_rate", "all"])
+def test_random_spikes_selection_no_unit(method):
+    # a sorting with no unit is valid and should give an empty selection, not raise
+    recording, sorting = generate_ground_truth_recording(
+        durations=[5.0],
+        sampling_frequency=16000.0,
+        num_channels=4,
+        num_units=3,
+        seed=2205,
+    )
+    empty_sorting = sorting.select_units([])
+    num_samples = [recording.get_num_samples(seg_index) for seg_index in range(recording.get_num_segments())]
+
+    random_spikes_indices = random_spikes_selection(
+        empty_sorting, num_samples, method=method, percentage=0.5, maximum_rate=10.0, seed=2205
+    )
+    assert random_spikes_indices.size == 0
+
+
 def test_apply_merges_to_sorting():
 
     times = np.array([0, 0, 10, 20, 300])

src/spikeinterface/core/tests/test_sortinganalyzer.py

@@ -736,6 +736,49 @@ def test_excess_spikes(dataset):
         create_sorting_analyzer(sorting=sorting, recording=recording.time_slice(0, 1))
 
 
+@pytest.mark.parametrize("sparse", [False, True])
+def test_analyzer_with_no_unit(dataset, sparse):
+    """
+    A sorting with no unit is a valid sorting, so the core extensions should run on it and
+    return empty results rather than raising.
+    """
+    recording, sorting = dataset
+    empty_sorting = sorting.select_units([])
+    assert len(empty_sorting.unit_ids) == 0
+
+    sorting_analyzer = create_sorting_analyzer(empty_sorting, recording, format="memory", sparse=sparse)
+    sorting_analyzer.compute(["random_spikes", "noise_levels", "waveforms", "templates"])
+
+    random_spikes = sorting_analyzer.get_extension("random_spikes").get_data()
+    assert random_spikes.shape == (0,)
+
+    waveforms = sorting_analyzer.get_extension("waveforms").get_data()
+    assert waveforms.shape[0] == 0
+
+    templates = sorting_analyzer.get_extension("templates").get_data()
+    assert templates.shape[0] == 0
+
+
+def test_analyzer_with_only_empty_units(dataset):
+    """
+    Units that exist but have no spike at all should give all-zero templates instead of raising.
+    """
+    from spikeinterface.core import NumpySorting
+
+    recording, _ = dataset
+    no_spikes = np.zeros(0, dtype="int64")
+    sorting = NumpySorting.from_samples_and_labels(
+        [no_spikes], [no_spikes], sampling_frequency=recording.sampling_frequency, unit_ids=np.array([0, 1])
+    )
+
+    sorting_analyzer = create_sorting_analyzer(sorting, recording, format="memory", sparse=False)
+    sorting_analyzer.compute(["random_spikes", "noise_levels", "templates"])
+
+    templates = sorting_analyzer.get_extension("templates").get_data()
+    assert templates.shape[0] == 2
+    assert np.all(templates == 0)
+
+
 def test_extensions_sorting():
 
     # nothing happens if all parents are on the left of the children

src/spikeinterface/core/waveform_tools.py

@@ -523,7 +523,8 @@ def extract_waveforms_to_single_buffer(
     if sparsity_mask is None:
         num_chans = recording.get_num_channels()
     else:
-        num_chans = int(max(np.sum(sparsity_mask, axis=1)))  # This is a numpy scalar, so we cast to int
+        # `initial` keeps this working for a sorting with no unit, where the mask has no row
+        num_chans = int(np.max(np.sum(sparsity_mask, axis=1), initial=0))  # This is a numpy scalar, so we cast to int
     shape = (int(num_spikes), int(n_samples), int(num_chans))
 
     if mode == "memmap":
@@ -907,17 +908,24 @@ def estimate_templates_with_accumulator(
         )
         return_in_uV = return_scaled
 
-    assert spikes.size > 0, "estimate_templates() need non empty sorting"
-
     job_kwargs = fix_job_kwargs(job_kwargs)
     num_worker = job_kwargs["n_jobs"]
 
     if sparsity_mask is None:
         num_chans = int(recording.get_num_channels())
     else:
-        num_chans = int(max(np.sum(sparsity_mask, axis=1)))  # This is a numpy scalar, so we cast to int
+        # `initial` keeps this working for a sorting with no unit, where the mask has no row
+        num_chans = int(np.max(np.sum(sparsity_mask, axis=1), initial=0))  # This is a numpy scalar, so we cast to int
     num_units = len(unit_ids)
 
+    if spikes.size == 0:
+        # A sorting with no unit (or with only empty units) is valid, there is simply nothing to
+        # accumulate. Returning zeros avoids allocating an empty shared memory buffer.
+        template_means = np.zeros((num_units, nbefore + nafter, num_chans), dtype="float32")
+        if return_std:
+            return template_means, np.zeros_like(template_means)
+        return template_means
+
     shape = (num_worker, num_units, nbefore + nafter, num_chans)
 
     dtype = np.dtype("float32")