SPL_single_units_reaching_task/code/python/H5_ISI.py

# The script calculates the Inter Spike Intervals (ISIs) and plots the
# results in a histogram. The proportion of ISIs violating a threshold
# is also calculated and indicated in the histogram.

# authors: Stefano Diomedi
# date: 10/2022

import numpy as np
import matplotlib.pyplot as plt
from supportFunctions import get_all_data_from_level_h5_rec
import tkinter as tk
from tkinter import filedialog
import os


current_dir = os.getcwd()

# Path to open the file dialog window in the data branch
path_folder_data = os.path.abspath(os.path.join(current_dir, "..", ".."))
path_folder_data = os.path.join(path_folder_data, "data")

# Select H5 dataset file
root = tk.Tk()
root.withdraw()
filename = filedialog.askopenfilename(initialdir=path_folder_data)

group_name = '/DATA'  # it can be 'unit_XXX', 'unit_XXX/condition_YY' or 'unit_XXX/condition_YY/trial_ZZ'
threshold = 1  # threshold for violations in ms

markers = []
all_strings_mk = []
spikes = []
all_strings_sp = []
markers, all_strings_mk, spikes, all_strings_sp = get_all_data_from_level_h5_rec.get_all_data_from_level_h5_rec(filename, group_name, markers, all_strings_mk, spikes, all_strings_sp)



ISIs = []

# Iterate over each array of spike data
for spike_data in spikes:
    # Calculate the differences for each spike_data array
    spike_data=np.array(spike_data)
    diff_data = np.diff(spike_data.flatten())
    # Append the differences to the list
    ISIs.append(diff_data)

# Concatenate all the differences into a single array
ISIs = np.concatenate(ISIs)

# Calculate the violations
violations = (ISIs < threshold).sum() / ISIs.size * 100


# Creation of Histogram
fig, ax = plt.subplots()
time = np.arange(0, 100.5, 0.5)
ax.hist(ISIs, bins=time)
ax.axvline(threshold, linestyle='--')
xLimits = ax.get_xlim()
yLimits = ax.get_ylim()
ax.text(xLimits[1] / 4, yLimits[1] * 2 / 3, f'Violations (ISI below {threshold}ms threshold) = {violations:.2f}%', fontsize=12)
ax.set_xlabel('ISI duration (ms)')
ax.set_ylabel('events num.')
ax.set_title('ISI')
plt.show()