Refractory Period Distance?

[karmenrai2]

when calculating refractory period violations, I am getting nan for several of the values. Does this mean there are no violations as I can't see any value with 0 or did the calculation fail and it's inconclusive?

import spikeinterface as si
import spikeinterface.extractors as se
import spikeinterface.widgets as sw
import spikeinterface.qualitymetrics as sqm
from spikeinterface.extractors import read_phy
from probeinterface import get_probe
from pathlib import Path
import numpy as np
import csv


# Cluster iterator
i = 0

# File paths 
rec_file_path = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\continuous1.dat") # Recording file
sort_file_path = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4") # Folder with the kilosort data
cluster_label = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\spike_clusters.npy") # spike_clusters npy file
pcs = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\spike_positions.npy") # spike_positions npy file

# Define recording parameters
sampling_frequency = 30_000.0  
num_channels = 384  
dtype = "float64"  

# Load data using SpikeInterface
recording = si.read_binary(file_paths=rec_file_path, sampling_frequency=sampling_frequency,
                           num_channels=num_channels, dtype=dtype)
sorting = read_phy(sort_file_path)

# Load data using numpy
all_labels = np.load(cluster_label)
all_pcs = np.load(pcs)

# Set the probe
other_probe = get_probe('cambridgeneurotech', 'ASSY-37-E-1')
other_probe.set_device_channel_indices(np.arange(32))
recording = recording.set_probe(other_probe, group_mode='by_shank')

# Create sorting analyzer
analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording, format="memory")

# Quality Metrics
presence_ratio = sqm.compute_presence_ratios(sorting_analyzer=analyzer)
refractory_period = sqm.compute_sliding_rp_violations(sorting_analyzer=analyzer, bin_size_ms=0.25, exclude_ref_period_below_ms=0.5)

# CSV file 
with open('spikesorting.csv', 'w') as file:
    writer = csv.writer(file)
    writer.writerow(["Cluster", "Presence Ratio", "Amplitude Cutoff", "Isolation Distance", "Refractory Period Violations"])
    for cluster in presence_ratio.keys():
        isolation_distance, _ = sqm.mahalanobis_metrics(all_pcs=all_pcs, all_labels=all_labels, this_unit_id=cluster)
        i+=1
        print(str(i) + "th cluster done")
        writer.writerow([cluster, presence_ratio[cluster], isolation_distance, refractory_period.get(cluster, "")])

# Print data for verification
print(recording)
print(sorting)
print(analyzer)
print(presence_ratio)
print(refractory_period)
exit()`

Cluster	Presence Ratio	Amplitude Cutoff	Isolation Distance	Refractory Period Violations
0	1	nan	1.89105441986244	0.11
1	1	nan	21.9754870562538	0.015
2	1	nan	2.08639005321337	0.01
3	1	nan	1.42021378063038	nan
4	1	nan	0.614645809462376	0.185
5	1	nan	0.209511255264518	nan
6	1	nan	0.548317966330339	0.265
7	1	nan	0.492024220726647	nan
8	1	nan	0.0672568742271259	nan
9	1	nan	2.09064727504505	0.065
10	1	nan	0.64246598552092	0.2
11	1	nan	1.75604923652422	0.035
12	1	nan	16.469823279801	0.02
13	1	nan	19.148817333648	0.025
14	1	nan	81.7195861151427	nan
15	1	nan	24.3185509096019	0.01
16	1	nan	2.01138983826387	0.02
17	1	nan	0.159849371240687	nan
18	1	nan	2.73207881475302	0.11
19	1	nan	6.44073858039389	0.025

[zm711]

Hey @karmenrai2

So if you look here:

spikeinterface/src/spikeinterface/qualitymetrics/misc_metrics.py

Lines 495 to 497 in d991382

	if unit_n_spikes <= min_spikes:
	contamination[unit_id] = np.nan
	continue

you see the assumption you like failed in using the sliding_rp. You could try using one of our other rp violation metrics which all have their different assumptions (we explain some of the assumptions/math in the docs see:
https://spikeinterface.readthedocs.io/en/stable/modules/qualitymetrics/sliding_rp_violations.html
and
https://spikeinterface.readthedocs.io/en/stable/modules/qualitymetrics/isi_violations.html

If you need more detail then that you could look in the actual code to see what is happening or ask more questions here :)

In general nan would indicate failed assumptions for any of our metrics (for example isolation distance fails if the covariance matrix is singular for a unit).

[zm711]

Out of curiosity is there a reason you want to do all of this at the low level rather than using the SortingAnalyzer object that will already great and save a pandas dataframe of all quality metrics for you?