sorting.py

spike_times_to_valid_samples(recording_times, spike_times, n_samples, unit_id)

Convert spike times (seconds) to sample indices within recording bounds.

Spike times are persisted in absolute seconds in NWB. On readback, floating-point rounding in the seconds-to-samples round-trip can cause np.searchsorted to return an index equal to n_samples (one past the last valid sample) for spikes at or near the end of the recording. SpikeInterface rejects such a sorting with ValueError: "The sorting object has spikes exceeding the recording duration". This helper drops those out-of-bounds indices and emits a warning identifying the affected unit and the count removed.

np.searchsorted is called with the default side='left': a spike that exactly equals recording_times[-1] maps to index n_samples - 1 (valid). Switching to side='right' would map the same spike to n_samples and reintroduce the bug.

Parameters:

Name	Type	Description	Default
recording_times	(ndarray, shape(n_samples))	Recording timestamps in seconds, monotonically increasing.	required
spike_times	(ndarray, shape(n_spikes))	Spike times for a single unit in seconds.	required
n_samples	int	Total number of samples in the recording.	required
unit_id	int or str	Identifier of the unit, used only in the warning message.	required

Returns:

Name	Type	Description
spike_samples	(ndarray, shape(n_valid_spikes))	Sample indices in [0, n_samples) corresponding to spike_times, with any out-of-bounds indices removed. n_valid_spikes <= n_spikes.

Source code in src/spyglass/spikesorting/v1/sorting.py

def spike_times_to_valid_samples(
    recording_times: np.ndarray,
    spike_times: np.ndarray,
    n_samples: int,
    unit_id,
) -> np.ndarray:
    """Convert spike times (seconds) to sample indices within recording bounds.

    Spike times are persisted in absolute seconds in NWB. On readback,
    floating-point rounding in the seconds-to-samples round-trip can cause
    ``np.searchsorted`` to return an index equal to ``n_samples`` (one past
    the last valid sample) for spikes at or near the end of the recording.
    SpikeInterface rejects such a sorting with ``ValueError: "The sorting
    object has spikes exceeding the recording duration"``. This helper drops
    those out-of-bounds indices and emits a warning identifying the affected
    unit and the count removed.

    ``np.searchsorted`` is called with the default ``side='left'``: a spike
    that exactly equals ``recording_times[-1]`` maps to index ``n_samples - 1``
    (valid). Switching to ``side='right'`` would map the same spike to
    ``n_samples`` and reintroduce the bug.

    Parameters
    ----------
    recording_times : np.ndarray, shape (n_samples,)
        Recording timestamps in seconds, monotonically increasing.
    spike_times : np.ndarray, shape (n_spikes,)
        Spike times for a single unit in seconds.
    n_samples : int
        Total number of samples in the recording.
    unit_id : int or str
        Identifier of the unit, used only in the warning message.

    Returns
    -------
    spike_samples : np.ndarray, shape (n_valid_spikes,)
        Sample indices in ``[0, n_samples)`` corresponding to ``spike_times``,
        with any out-of-bounds indices removed. ``n_valid_spikes <= n_spikes``.
    """
    spike_samples = np.searchsorted(recording_times, spike_times)
    excess_mask = spike_samples >= n_samples
    n_excess = int(excess_mask.sum())
    if n_excess > 0:
        logger.warning(
            f"Unit {unit_id} has {n_excess} spike(s) exceeding the "
            "recording duration. Removing excess spikes. This may be "
            "caused by floating-point rounding during the seconds-to-"
            "samples conversion."
        )
        spike_samples = spike_samples[~excess_mask]
    return spike_samples

SpikeSorterParameters

Bases: SpyglassMixin, Lookup

Parameters for spike sorting algorithms.

Attributes:

Name	Type	Description
sorter	str	Name of the spike sorting algorithm.
sorter_params_name	str	Name of the parameter set for the spike sorting algorithm.
sorter_params	dict	Dictionary of parameters for the spike sorting algorithm. The keys and values depend on the specific algorithm being used. For example, for the "mountainsort4" algorithm, the parameters are... detect_sign: int Sign of the detected spikes. 1 for positive, -1 for negative. adjacency_radius: int Radius for adjacency graph. Determines which channels are considered neighbors. freq_min: int Minimum frequency for bandpass filter. freq_max: int Maximum frequency for bandpass filter. filter: bool Whether to apply bandpass filter. whiten: bool Whether to whiten the data. num_workers: int Number of workers to use for parallel processing. clip_size: int Size of the clips to extract for spike detection. detect_threshold: float Threshold for spike detection. detect_interval: int Minimum interval between detected spikes. For the "clusterless_thresholder" algorithm, the parameters are... detect_threshold: float microvolt detection threshold for spike detection. method: str Method for spike detection. Options are "locally_exclusive" or "global". peak_sign: enum ("neg", "pos") Sign of the detected peaks. exclude_sweep_ms: float Exclusion time in milliseconds for detected spikes. local_radius_um: int Local radius in micrometers for spike detection. noise_levels: np.ndarray Noise levels for spike detection. random_chunk_kwargs: dict Additional arguments for random chunk processing. outputs: str Output type for spike detection. Options are "sorting" or "labels".

Source code in src/spyglass/spikesorting/v1/sorting.py

@schema
class SpikeSorterParameters(SpyglassMixin, dj.Lookup):
    """Parameters for spike sorting algorithms.

    Attributes
    ----------
    sorter: str
        Name of the spike sorting algorithm.
    sorter_params_name: str
        Name of the parameter set for the spike sorting algorithm.
    sorter_params: dict
        Dictionary of parameters for the spike sorting algorithm.
        The keys and values depend on the specific algorithm being used.
        For example, for the "mountainsort4" algorithm, the parameters are...
            detect_sign: int
                Sign of the detected spikes. 1 for positive, -1 for negative.
            adjacency_radius: int
                Radius for adjacency graph. Determines which channels are
                considered neighbors.
            freq_min: int
                Minimum frequency for bandpass filter.
            freq_max: int
                Maximum frequency for bandpass filter.
            filter: bool
                Whether to apply bandpass filter.
            whiten: bool
                Whether to whiten the data.
            num_workers: int
                Number of workers to use for parallel processing.
            clip_size: int
                Size of the clips to extract for spike detection.
            detect_threshold: float
                Threshold for spike detection.
            detect_interval: int
                Minimum interval between detected spikes.
        For the "clusterless_thresholder" algorithm, the parameters are...
            detect_threshold: float
                microvolt detection threshold for spike detection.
            method: str
                Method for spike detection. Options are "locally_exclusive" or
                "global".
            peak_sign: enum ("neg", "pos")
                Sign of the detected peaks.
            exclude_sweep_ms: float
                Exclusion time in milliseconds for detected spikes.
            local_radius_um: int
                Local radius in micrometers for spike detection.
            noise_levels: np.ndarray
                Noise levels for spike detection.
            random_chunk_kwargs: dict
                Additional arguments for random chunk processing.
            outputs: str
                Output type for spike detection. Options are "sorting" or
                "labels".
    """

    definition = """
    # Spike sorting algorithm and associated parameters.
    sorter: varchar(200)
    sorter_param_name: varchar(200)
    ---
    sorter_params: blob
    """
    mountain_default = {
        "detect_sign": -1,
        "adjacency_radius": 100,
        "filter": False,
        "whiten": True,
        "num_workers": 1,
        "clip_size": 40,
        "detect_threshold": 3,
        "detect_interval": 10,
    }
    contents = [
        [
            "mountainsort4",
            "franklab_tetrode_hippocampus_30KHz",
            {**mountain_default, "freq_min": 600, "freq_max": 6000},
        ],
        [
            "mountainsort4",
            "franklab_probe_ctx_30KHz",
            {**mountain_default, "freq_min": 300, "freq_max": 6000},
        ],
        [
            "clusterless_thresholder",
            "default_clusterless",
            {
                "detect_threshold": 100.0,  # uV
                # Locally exclusive means one unit per spike detected
                "method": "locally_exclusive",
                "peak_sign": "neg",
                "exclude_sweep_ms": 0.1,
                "local_radius_um": 100,
                # noise levels needs to be 1.0 so the units are in uV and not MAD
                "noise_levels": np.asarray([1.0]),
                "random_chunk_kwargs": {},
                # output needs to be set to sorting for the rest of the pipeline
                "outputs": "sorting",
            },
        ],
    ]
    contents.extend(
        [
            [sorter, "default", sis.get_default_sorter_params(sorter)]
            for sorter in sis.available_sorters()
        ]
    )

    @classmethod
    def insert_default(cls):
        """Insert default sorter parameters into SpikeSorterParameters table."""
        cls.insert(cls.contents, skip_duplicates=True)

insert_default() classmethod

Insert default sorter parameters into SpikeSorterParameters table.

Source code in src/spyglass/spikesorting/v1/sorting.py

@classmethod
def insert_default(cls):
    """Insert default sorter parameters into SpikeSorterParameters table."""
    cls.insert(cls.contents, skip_duplicates=True)

SpikeSortingSelection

Bases: SpyglassMixin, Manual

Source code in src/spyglass/spikesorting/v1/sorting.py

@schema
class SpikeSortingSelection(SpyglassMixin, dj.Manual):
    definition = """
    # Processed recording and spike sorting parameters. See `insert_selection`.
    sorting_id: uuid
    ---
    -> SpikeSortingRecording
    -> SpikeSorterParameters
    -> IntervalList
    """

    @classmethod
    def insert_selection(cls, key: dict):
        """Insert a row into SpikeSortingSelection with an
        automatically generated unique sorting ID as the sole primary key.

        Parameters
        ----------
        key : dict
            primary key of SpikeSortingRecording, SpikeSorterParameters, IntervalList tables

        Returns
        -------
        sorting_id : uuid
            the unique sorting ID serving as primary key for SpikeSorting
        """
        query = cls & key
        if query:
            logger.info("Similar row(s) already inserted.")
            return query.fetch(as_dict=True)
        key["sorting_id"] = uuid.uuid4()
        cls.insert1(key, skip_duplicates=True)
        return key

insert_selection(key) classmethod

Insert a row into SpikeSortingSelection with an automatically generated unique sorting ID as the sole primary key.

Parameters:

Name	Type	Description	Default
key	dict	primary key of SpikeSortingRecording, SpikeSorterParameters, IntervalList tables	required

Returns:

Name	Type	Description
sorting_id	uuid	the unique sorting ID serving as primary key for SpikeSorting

Source code in src/spyglass/spikesorting/v1/sorting.py

@classmethod
def insert_selection(cls, key: dict):
    """Insert a row into SpikeSortingSelection with an
    automatically generated unique sorting ID as the sole primary key.

    Parameters
    ----------
    key : dict
        primary key of SpikeSortingRecording, SpikeSorterParameters, IntervalList tables

    Returns
    -------
    sorting_id : uuid
        the unique sorting ID serving as primary key for SpikeSorting
    """
    query = cls & key
    if query:
        logger.info("Similar row(s) already inserted.")
        return query.fetch(as_dict=True)
    key["sorting_id"] = uuid.uuid4()
    cls.insert1(key, skip_duplicates=True)
    return key

SpikeSorting

Bases: SpyglassMixin, Computed

Source code in src/spyglass/spikesorting/v1/sorting.py

@schema
class SpikeSorting(SpyglassMixin, dj.Computed):
    definition = """
    -> SpikeSortingSelection
    ---
    -> AnalysisNwbfile
    object_id: varchar(40)          # Object ID for the sorting in NWB file
    time_of_sort: int               # in Unix time, to the nearest second
    """

    _parallel_make = True  # True if n_workers > 1

    def make_fetch(self, key: dict) -> list:
        """Runs spike sorting on the data and parameters specified by the
        SpikeSortingSelection table and inserts a new entry to SpikeSorting table.
        """
        # FETCH
        # - information about the recording
        # - artifact free intervals
        # - spike sorter and sorter params

        recording_key = (
            SpikeSortingRecording * SpikeSortingSelection & key
        ).fetch1()

        nwb_file_name = recording_key["nwb_file_name"]

        artifact_removed_intervals = (
            IntervalList
            & {
                "nwb_file_name": nwb_file_name,
                "interval_list_name": recording_key["interval_list_name"],
            }
        ).fetch1("valid_times")

        sorter, sorter_params = (
            SpikeSorterParameters * SpikeSortingSelection & key
        ).fetch1("sorter", "sorter_params")

        return [
            nwb_file_name,
            artifact_removed_intervals,
            sorter,
            sorter_params,
            recording_key,
        ]

    def make_compute(
        self,
        key: dict,
        nwb_file_name: str,
        artifact_removed_intervals: IntervalLike,
        sorter: str,
        sorter_params: dict,
        recording_key: dict,
    ):
        sorting, timestamps = self._run_spike_sorter(
            recording_key=recording_key,
            artifact_removed_intervals=artifact_removed_intervals,
            sorter=sorter,
            sorter_params=sorter_params,
        )

        time_of_sort = int(time.time())
        analysis_file_name, object_id = self._save_sorting_results(
            sorting=sorting,
            timestamps=timestamps,
            artifact_removed_intervals=artifact_removed_intervals,
            nwb_file_name=nwb_file_name,
        )

        return [nwb_file_name, time_of_sort, analysis_file_name, object_id]

    def make_insert(
        self,
        key: dict,
        nwb_file_name: str,
        time_of_sort: int,
        analysis_file_name: str,
        object_id: str,
    ):
        AnalysisNwbfile().add(nwb_file_name, analysis_file_name)
        self.insert1(
            dict(
                key,
                time_of_sort=time_of_sort,
                analysis_file_name=analysis_file_name,
                object_id=object_id,
            ),
            skip_duplicates=True,
        )

    def _run_spike_sorter(
        self,
        recording_key,
        artifact_removed_intervals,
        sorter,
        sorter_params,
    ):
        """Run spike sorting algorithm (external dependency).

        This method wraps all calls to spikeinterface for spike sorting,
        making it easy to mock in tests for faster execution.

        Parameters
        ----------
        recording_key : dict
            Key for the recording
        artifact_removed_intervals : np.ndarray
            Artifact-free time intervals
        sorter : str
            Name of spike sorter algorithm
        sorter_params : dict
            Parameters for spike sorter

        Returns
        -------
        sorting : si.BaseSorting
            Sorted spike times
        timestamps : np.ndarray
            Recording timestamps
        """
        # Load recording (spikeinterface)
        recording = SpikeSortingRecording().get_recording(recording_key)

        timestamps = recording.get_times()

        artifact_removed_intervals_ind = _consolidate_intervals(
            artifact_removed_intervals, timestamps
        )

        # Remove artifacts if needed (spikeinterface)
        if (
            (len(artifact_removed_intervals_ind) > 1)
            or (artifact_removed_intervals_ind[0][0] > 0)
            or (artifact_removed_intervals_ind[-1][1] < len(timestamps))
        ):
            # set the artifact intervals to zero
            list_triggers = []
            if artifact_removed_intervals_ind[0][0] > 0:
                list_triggers.append(
                    np.arange(0, artifact_removed_intervals_ind[0][0])
                )
            for interval_ind in range(len(artifact_removed_intervals_ind) - 1):
                list_triggers.append(
                    np.arange(
                        (artifact_removed_intervals_ind[interval_ind][1] + 1),
                        artifact_removed_intervals_ind[interval_ind + 1][0],
                    )
                )
            if artifact_removed_intervals_ind[-1][1] < len(timestamps):
                list_triggers.append(
                    np.arange(
                        artifact_removed_intervals_ind[-1][1],
                        len(timestamps) - 1,
                    )
                )

            list_triggers = [list(np.concatenate(list_triggers))]
            recording = sip.remove_artifacts(
                recording=recording,
                list_triggers=list_triggers,
                ms_before=None,
                ms_after=None,
                mode="zeros",
            )

        # Run spike sorting (spikeinterface)
        if sorter == "clusterless_thresholder":
            # need to remove tempdir and whiten from sorter_params
            sorter_params.pop("tempdir", None)
            sorter_params.pop("whiten", None)
            sorter_params.pop("outputs", None)
            if "local_radius_um" in sorter_params:
                sorter_params["radius_um"] = sorter_params.pop(
                    "local_radius_um"
                )  # correct existing parameter sets for spikeinterface>=0.99.1

            # Detect peaks for clusterless decoding
            detected_spikes = detect_peaks(recording, **sorter_params)
            sorting = si.NumpySorting.from_times_labels(
                times_list=detected_spikes["sample_index"],
                labels_list=np.zeros(len(detected_spikes), dtype=np.int32),
                sampling_frequency=recording.get_sampling_frequency(),
            )
        else:
            sorter_temp_dir = tempfile.TemporaryDirectory(dir=temp_dir)
            os.chmod(sorter_temp_dir.name, 0o777)

            # Only mountainsort4 declares a `tempdir` scratch-dir param. Passing
            # `tempdir` to any other sorter makes its parameter validation raise
            # `AttributeError: Bad parameters: ['tempdir']`, so inject it solely
            # for sorters that actually declare it. The temp dir is still handed
            # to every sorter below as `output_folder`.
            if "tempdir" in sis.get_default_sorter_params(sorter):
                sorter_params["tempdir"] = sorter_temp_dir.name

            # if whitening is specified in sorter params, apply whitening separately
            # prior to sorting and turn off "sorter whitening"
            if sorter_params.get("whiten", False):
                recording = sip.whiten(recording, dtype=np.float64)
                sorter_params["whiten"] = False

            common_sorter_items = {
                "sorter_name": sorter,
                "recording": recording,
                "output_folder": sorter_temp_dir.name,
                "remove_existing_folder": True,
            }

            if sorter.lower() in ["kilosort2_5", "kilosort3", "ironclust"]:
                sorter_params = {
                    k: v
                    for k, v in sorter_params.items()
                    if k not in ["mp_context", "max_threads_per_process"]
                }
                sorting = sis.run_sorter(
                    **common_sorter_items,
                    singularity_image=True,
                    **sorter_params,
                )
            else:
                sorting = sis.run_sorter(
                    **common_sorter_items,
                    **sorter_params,
                )

        sorting = sic.remove_excess_spikes(sorting, recording)

        return sorting, timestamps

    def _save_sorting_results(
        self,
        sorting,
        timestamps,
        artifact_removed_intervals,
        nwb_file_name,
    ):
        """Save sorting results to NWB file (external I/O).

        This method wraps file I/O operations, making it easy to
        mock in tests to avoid filesystem dependencies.

        Parameters
        ----------
        sorting : si.BaseSorting
            Sorted spike times
        timestamps : np.ndarray
            Recording timestamps
        artifact_removed_intervals : np.ndarray
            Artifact-free time intervals
        nwb_file_name : str
            Name of source NWB file

        Returns
        -------
        analysis_file_name : str
            Name of analysis NWB file
        object_id : str
            Object ID in NWB file
        """
        return _write_sorting_to_nwb(
            sorting,
            timestamps,
            artifact_removed_intervals,
            nwb_file_name,
        )

    @classmethod
    def get_sorting(cls, key: dict) -> si.BaseSorting:
        """Get sorting in the analysis NWB file as spikeinterface BaseSorting

        Parameters
        ----------
        key : dict
            primary key of SpikeSorting

        Returns
        -------
        sorting : si.BaseSorting

        """

        recording_id = (
            SpikeSortingRecording * SpikeSortingSelection & key
        ).fetch1("recording_id")
        recording = SpikeSortingRecording.get_recording(
            {"recording_id": recording_id}
        )
        sampling_frequency = recording.get_sampling_frequency()
        analysis_file_name = (cls & key).fetch1("analysis_file_name")
        analysis_file_abs_path = AnalysisNwbfile.get_abs_path(
            analysis_file_name
        )
        with pynwb.NWBHDF5IO(
            analysis_file_abs_path, "r", load_namespaces=True
        ) as io:
            nwbf = io.read()
            units = nwbf.units.to_dataframe()

        recording_times = recording.get_times()
        n_samples = recording.get_num_samples()
        units_dict = {
            unit_id: spike_times_to_valid_samples(
                recording_times, spike_times, n_samples, unit_id
            )
            for unit_id, spike_times in zip(units.index, units["spike_times"])
        }

        sorting = si.NumpySorting.from_unit_dict(
            [units_dict], sampling_frequency=sampling_frequency
        )

        return sorting

make_fetch(key)

Runs spike sorting on the data and parameters specified by the SpikeSortingSelection table and inserts a new entry to SpikeSorting table.

Source code in src/spyglass/spikesorting/v1/sorting.py

def make_fetch(self, key: dict) -> list:
    """Runs spike sorting on the data and parameters specified by the
    SpikeSortingSelection table and inserts a new entry to SpikeSorting table.
    """
    # FETCH
    # - information about the recording
    # - artifact free intervals
    # - spike sorter and sorter params

    recording_key = (
        SpikeSortingRecording * SpikeSortingSelection & key
    ).fetch1()

    nwb_file_name = recording_key["nwb_file_name"]

    artifact_removed_intervals = (
        IntervalList
        & {
            "nwb_file_name": nwb_file_name,
            "interval_list_name": recording_key["interval_list_name"],
        }
    ).fetch1("valid_times")

    sorter, sorter_params = (
        SpikeSorterParameters * SpikeSortingSelection & key
    ).fetch1("sorter", "sorter_params")

    return [
        nwb_file_name,
        artifact_removed_intervals,
        sorter,
        sorter_params,
        recording_key,
    ]

get_sorting(key) classmethod

Get sorting in the analysis NWB file as spikeinterface BaseSorting

Parameters:

Name	Type	Description	Default
key	dict	primary key of SpikeSorting	required

Returns:

Name	Type	Description
sorting	BaseSorting

Source code in src/spyglass/spikesorting/v1/sorting.py

@classmethod
def get_sorting(cls, key: dict) -> si.BaseSorting:
    """Get sorting in the analysis NWB file as spikeinterface BaseSorting

    Parameters
    ----------
    key : dict
        primary key of SpikeSorting

    Returns
    -------
    sorting : si.BaseSorting

    """

    recording_id = (
        SpikeSortingRecording * SpikeSortingSelection & key
    ).fetch1("recording_id")
    recording = SpikeSortingRecording.get_recording(
        {"recording_id": recording_id}
    )
    sampling_frequency = recording.get_sampling_frequency()
    analysis_file_name = (cls & key).fetch1("analysis_file_name")
    analysis_file_abs_path = AnalysisNwbfile.get_abs_path(
        analysis_file_name
    )
    with pynwb.NWBHDF5IO(
        analysis_file_abs_path, "r", load_namespaces=True
    ) as io:
        nwbf = io.read()
        units = nwbf.units.to_dataframe()

    recording_times = recording.get_times()
    n_samples = recording.get_num_samples()
    units_dict = {
        unit_id: spike_times_to_valid_samples(
            recording_times, spike_times, n_samples, unit_id
        )
        for unit_id, spike_times in zip(units.index, units["spike_times"])
    }

    sorting = si.NumpySorting.from_unit_dict(
        [units_dict], sampling_frequency=sampling_frequency
    )

    return sorting