Delete code/GP_Code.py
This commit is contained in:
parent
53358acbde
commit
337b9d8e0e
329
code/GP_Code.py
329
code/GP_Code.py
@ -1,329 +0,0 @@
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Created on Tue Oct 22 15:21:41 2024
|
||||
|
||||
@author: diana
|
||||
"""
|
||||
|
||||
import glob
|
||||
import os
|
||||
import rlxnix as rlx
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import scipy.signal as sig
|
||||
from scipy.integrate import quad
|
||||
|
||||
|
||||
### FUNCTIONS ###
|
||||
def binary_spikes(spike_times, duration, dt):
|
||||
"""Converts the spike times to a binary representation.
|
||||
Zeros when there is no spike, one when there is.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
spike_times : np.array
|
||||
The spike times.
|
||||
duration : float
|
||||
The trial duration.
|
||||
dt : float
|
||||
The temporal resolution.
|
||||
|
||||
Returns
|
||||
-------
|
||||
binary : np.array
|
||||
The binary representation of the spike times.
|
||||
"""
|
||||
binary = np.zeros(int(np.round(duration / dt))) #Vektor, der genauso lang ist wie die stim time
|
||||
spike_indices = np.asarray(np.round(spike_times / dt), dtype=int)
|
||||
binary[spike_indices] = 1
|
||||
return binary
|
||||
|
||||
|
||||
def firing_rate(binary_spikes, box_width, dt=0.000025):
|
||||
"""Calculate the firing rate from binary spike data.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
binary_spikes : np.array
|
||||
A binary array representing spike occurrences.
|
||||
box_width : float
|
||||
The width of the box filter in seconds.
|
||||
dt : float, optional
|
||||
The temporal resolution (time step) in seconds. Default is 0.000025 seconds.
|
||||
|
||||
Returns
|
||||
-------
|
||||
rate : np.array
|
||||
An array representing the firing rate at each time step.
|
||||
"""
|
||||
box = np.ones(int(box_width // dt))
|
||||
box /= np.sum(box) * dt # Normalization of box kernel to an integral of 1
|
||||
rate = np.convolve(binary_spikes, box, mode="same")
|
||||
return rate
|
||||
|
||||
|
||||
def powerspectrum(rate, dt):
|
||||
"""Compute the power spectrum of a given firing rate.
|
||||
|
||||
This function calculates the power spectrum using the Welch method.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rate : np.array
|
||||
An array of firing rates.
|
||||
dt : float
|
||||
The temporal resolution (time step) in seconds.
|
||||
|
||||
Returns
|
||||
-------
|
||||
frequency : np.array
|
||||
An array of frequencies corresponding to the power values.
|
||||
power : np.array
|
||||
An array of power spectral density values.
|
||||
"""
|
||||
frequency, power = sig.welch(rate, fs=1/dt, nperseg=2**15, noverlap=2**14)
|
||||
return frequency, power
|
||||
|
||||
|
||||
def calculate_integral(frequency, power, point, delta):
|
||||
"""
|
||||
Calculate the integral around a single specified point.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
frequency : np.array
|
||||
An array of frequencies corresponding to the power values.
|
||||
power : np.array
|
||||
An array of power spectral density values.
|
||||
point : float
|
||||
The harmonic frequency at which to calculate the integral.
|
||||
delta : float
|
||||
Half-width of the range for integration around the point.
|
||||
|
||||
Returns
|
||||
-------
|
||||
integral : float
|
||||
The calculated integral around the point.
|
||||
local_mean : float
|
||||
The local mean value (adjacent integrals).
|
||||
"""
|
||||
indices = (frequency >= point - delta) & (frequency <= point + delta)
|
||||
integral = np.trapz(power[indices], frequency[indices])
|
||||
|
||||
left_indices = (frequency >= point - 5 * delta) & (frequency < point - delta)
|
||||
right_indices = (frequency > point + delta) & (frequency <= point + 5 * delta)
|
||||
|
||||
l_integral = np.trapz(power[left_indices], frequency[left_indices])
|
||||
r_integral = np.trapz(power[right_indices], frequency[right_indices])
|
||||
|
||||
local_mean = np.mean([l_integral, r_integral])
|
||||
return integral, local_mean
|
||||
|
||||
|
||||
def valid_integrals(integral, local_mean, threshold, point):
|
||||
"""
|
||||
Check if the integral exceeds the threshold compared to the local mean and
|
||||
provide feedback on whether the given point is valid or not.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
integral : float
|
||||
The calculated integral around the point.
|
||||
local_mean : float
|
||||
The local mean value (adjacent integrals).
|
||||
threshold : float
|
||||
Threshold value to compare integrals with local mean.
|
||||
point : float
|
||||
The harmonic frequency point being evaluated.
|
||||
|
||||
Returns
|
||||
-------
|
||||
valid : bool
|
||||
True if the integral exceeds the local mean by the threshold, otherwise False.
|
||||
message : str
|
||||
A message stating whether the point is valid or not.
|
||||
"""
|
||||
valid = integral > (local_mean * threshold)
|
||||
if valid:
|
||||
message = f"The point {point} is valid, as its integral exceeds the threshold."
|
||||
else:
|
||||
message = f"The point {point} is not valid, as its integral does not exceed the threshold."
|
||||
return valid, message
|
||||
|
||||
|
||||
def prepare_harmonics(frequencies, categories, num_harmonics, colors):
|
||||
"""
|
||||
Prepare harmonic frequencies and assign colors based on categories.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
frequencies : list
|
||||
Base frequencies to generate harmonics.
|
||||
categories : list
|
||||
Corresponding categories for the base frequencies.
|
||||
num_harmonics : list
|
||||
Number of harmonics for each base frequency.
|
||||
colors : list
|
||||
List of colors corresponding to the categories.
|
||||
|
||||
Returns
|
||||
-------
|
||||
points : list
|
||||
A flat list of harmonic frequencies.
|
||||
color_mapping : dict
|
||||
A dictionary mapping each category to its corresponding color.
|
||||
points_categories : dict
|
||||
A mapping of categories to their harmonic frequencies.
|
||||
"""
|
||||
points_categories = {}
|
||||
for idx, (freq, category) in enumerate(zip(frequencies, categories)):
|
||||
points_categories[category] = [freq * (i + 1) for i in range(num_harmonics[idx])]
|
||||
|
||||
points = [p for harmonics in points_categories.values() for p in harmonics]
|
||||
color_mapping = {category: colors[idx] for idx, category in enumerate(categories)}
|
||||
|
||||
return points, color_mapping, points_categories
|
||||
|
||||
|
||||
def find_exceeding_points(frequency, power, points, delta, threshold):
|
||||
"""
|
||||
Find the points where the integral exceeds the local mean by a given threshold.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
frequency : np.array
|
||||
An array of frequencies corresponding to the power values.
|
||||
power : np.array
|
||||
An array of power spectral density values.
|
||||
points : list
|
||||
A list of harmonic frequencies to evaluate.
|
||||
delta : float
|
||||
Half-width of the range for integration around the point.
|
||||
threshold : float
|
||||
Threshold value to compare integrals with local mean.
|
||||
|
||||
Returns
|
||||
-------
|
||||
exceeding_points : list
|
||||
A list of points where the integral exceeds the local mean by the threshold.
|
||||
"""
|
||||
exceeding_points = []
|
||||
|
||||
for point in points:
|
||||
# Calculate the integral and local mean for the current point
|
||||
integral, local_mean = calculate_integral(frequency, power, point, delta)
|
||||
|
||||
# Check if the integral exceeds the threshold
|
||||
valid, message = valid_integrals(integral, local_mean, threshold, point)
|
||||
|
||||
if valid:
|
||||
exceeding_points.append(point)
|
||||
|
||||
return exceeding_points
|
||||
|
||||
|
||||
def plot_highlighted_integrals(frequency, power, exceeding_points, delta, threshold, color_mapping, points_categories):
|
||||
"""
|
||||
Plot the power spectrum and highlight integrals that exceed the threshold.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
frequency : np.array
|
||||
An array of frequencies corresponding to the power values.
|
||||
power : np.array
|
||||
An array of power spectral density values.
|
||||
exceeding_points : list
|
||||
A list of harmonic frequencies that exceed the threshold.
|
||||
delta : float
|
||||
Half-width of the range for integration around each point.
|
||||
threshold : float
|
||||
Threshold value to compare integrals with local mean.
|
||||
color_mapping : dict
|
||||
A dictionary mapping each category to its color.
|
||||
points_categories : dict
|
||||
A mapping of categories to lists of points.
|
||||
|
||||
Returns
|
||||
-------
|
||||
fig : matplotlib.figure.Figure
|
||||
The created figure object with highlighted integrals.
|
||||
"""
|
||||
fig, ax = plt.subplots()
|
||||
ax.plot(frequency, power) # Plot power spectrum
|
||||
|
||||
for point in exceeding_points:
|
||||
integral, local_mean = calculate_integral(frequency, power, point, delta)
|
||||
valid, _ = valid_integrals(integral, local_mean, threshold, point)
|
||||
if valid:
|
||||
# Define color based on the category of the point
|
||||
color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
|
||||
# Shade the region around the point where the integral was calculated
|
||||
ax.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz')
|
||||
print(f"Integral around {point:.2f} Hz: {integral:.5e}")
|
||||
|
||||
|
||||
# Define left and right boundaries of adjacent regions
|
||||
left_boundary = frequency[np.where((frequency >= point - 5 * delta) & (frequency < point - delta))[0][0]]
|
||||
right_boundary = frequency[np.where((frequency > point + delta) & (frequency <= point + 5 * delta))[0][-1]]
|
||||
|
||||
# Add vertical dashed lines at the boundaries of the adjacent regions
|
||||
ax.axvline(x=left_boundary, color="k", linestyle="--")
|
||||
ax.axvline(x=right_boundary, color="k", linestyle="--")
|
||||
|
||||
|
||||
ax.set_xlim([0, 1200])
|
||||
ax.set_xlabel('Frequency (Hz)')
|
||||
ax.set_ylabel('Power')
|
||||
ax.set_title('Power Spectrum with Highlighted Integrals')
|
||||
ax.legend()
|
||||
|
||||
return fig
|
||||
|
||||
|
||||
### Data retrieval ###
|
||||
datafolder = "../data"
|
||||
example_file = os.path.join("..", "data", "2024-10-16-ad-invivo-1.nix")
|
||||
dataset = rlx.Dataset(example_file)
|
||||
sams = dataset.repro_runs("SAM")
|
||||
sam = sams[2]
|
||||
|
||||
## Data for functions
|
||||
df = sam.metadata["RePro-Info"]["settings"]["deltaf"][0][0]
|
||||
stim = sam.stimuli[1]
|
||||
potential, time = stim.trace_data("V-1")
|
||||
spikes, _ = stim.trace_data("Spikes-1")
|
||||
duration = stim.duration
|
||||
dt = stim.trace_info("V-1").sampling_interval
|
||||
|
||||
|
||||
### Apply Functions to calculate data ###
|
||||
b = binary_spikes(spikes, duration, dt)
|
||||
rate = firing_rate(b, box_width=0.05, dt=dt)
|
||||
frequency, power = powerspectrum(b, dt)
|
||||
|
||||
|
||||
### Important stuff ###
|
||||
## Frequencies
|
||||
eodf = stim.metadata[stim.name]["EODf"][0][0]
|
||||
stimulus_frequency = eodf + df
|
||||
AM = 50 # Hz
|
||||
frequencies = [AM, eodf, stimulus_frequency]
|
||||
|
||||
categories = ["AM", "EODf", "Stimulus frequency"]
|
||||
num_harmonics = [4, 2, 2]
|
||||
colors = ["green", "orange", "red"]
|
||||
|
||||
delta = 2.5
|
||||
threshold = 10
|
||||
|
||||
### Apply functions to make powerspectrum ###
|
||||
integral, local = calculate_integral(frequency, power, eodf, delta)
|
||||
valid = valid_integrals(integral, local, threshold, eodf)
|
||||
points, color, categories = prepare_harmonics(frequencies, categories, num_harmonics, colors)
|
||||
print(len(points))
|
||||
exceeding = find_exceeding_points(frequency, power, points, delta, threshold)
|
||||
print(len(exceeding))
|
||||
|
||||
## Plot power spectrum and highlight integrals
|
||||
fig = plot_highlighted_integrals(frequency, power, points, delta, threshold, color, categories)
|
||||
plt.show()
|
||||
Loading…
Reference in New Issue
Block a user