This commit is contained in:
mbergmann 2024-10-23 10:01:15 +02:00
commit 23db938fdf

View File

@ -1,11 +1,10 @@
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""" """
Created on Tue Oct 22 11:43:41 2024 Created on Tue Oct 22 15:21:41 2024
@author: diana @author: diana
""" """
import glob import glob
import os import os
import rlxnix as rlx import rlxnix as rlx
@ -17,8 +16,8 @@ from scipy.integrate import quad
### FUNCTIONS ### ### FUNCTIONS ###
def binary_spikes(spike_times, duration, dt): def binary_spikes(spike_times, duration, dt):
""" Converts the spike times to a binary representation. """Converts the spike times to a binary representation.
Zeros when there is no spike, One when there is. Zeros when there is no spike, one when there is.
Parameters Parameters
---------- ----------
@ -27,15 +26,14 @@ def binary_spikes(spike_times, duration, dt):
duration : float duration : float
The trial duration. The trial duration.
dt : float dt : float
the temporal resolution. The temporal resolution.
Returns Returns
------- -------
binary : np.array binary : np.array
The binary representation of the spike times. The binary representation of the spike times.
""" """
binary = np.zeros(int(np.round(duration / dt))) #Vektor, der genauso lang ist wie die stim time 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) spike_indices = np.asarray(np.round(spike_times / dt), dtype=int)
binary[spike_indices] = 1 binary[spike_indices] = 1
return binary return binary
@ -44,9 +42,6 @@ def binary_spikes(spike_times, duration, dt):
def firing_rate(binary_spikes, box_width, dt=0.000025): def firing_rate(binary_spikes, box_width, dt=0.000025):
"""Calculate the firing rate from binary spike data. """Calculate the firing rate from binary spike data.
This function computes the firing rate using a boxcar (moving average)
filter of a specified width.
Parameters Parameters
---------- ----------
binary_spikes : np.array binary_spikes : np.array
@ -62,7 +57,7 @@ def firing_rate(binary_spikes, box_width, dt=0.000025):
An array representing the firing rate at each time step. An array representing the firing rate at each time step.
""" """
box = np.ones(int(box_width // dt)) box = np.ones(int(box_width // dt))
box /= np.sum(box) * dt #Normalization of box kernel to an integral of 1 box /= np.sum(box) * dt # Normalization of box kernel to an integral of 1
rate = np.convolve(binary_spikes, box, mode="same") rate = np.convolve(binary_spikes, box, mode="same")
return rate return rate
@ -90,7 +85,96 @@ def powerspectrum(rate, dt):
return frequency, power 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): 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 = {} points_categories = {}
for idx, (freq, category) in enumerate(zip(frequencies, categories)): for idx, (freq, category) in enumerate(zip(frequencies, categories)):
points_categories[category] = [freq * (i + 1) for i in range(num_harmonics[idx])] points_categories[category] = [freq * (i + 1) for i in range(num_harmonics[idx])]
@ -101,12 +185,9 @@ def prepare_harmonics(frequencies, categories, num_harmonics, colors):
return points, color_mapping, points_categories return points, color_mapping, points_categories
def plot_power_spectrum_with_integrals(frequency, power, points, delta): def find_exceeding_points(frequency, power, points, delta, threshold):
""" """
Create a figure of the power spectrum and calculate integrals around specified points. Find the points where the integral exceeds the local mean by a given threshold.
This function generates the plot of the power spectrum and calculates integrals
around specified harmonic points, but it does not color the regions or add vertical lines.
Parameters Parameters
---------- ----------
@ -115,60 +196,35 @@ def plot_power_spectrum_with_integrals(frequency, power, points, delta):
power : np.array power : np.array
An array of power spectral density values. An array of power spectral density values.
points : list points : list
A list of harmonic frequencies to highlight. A list of harmonic frequencies to evaluate.
delta : float delta : float
Half-width of the range for integration around each point. Half-width of the range for integration around the point.
threshold : float
Threshold value to compare integrals with local mean.
Returns Returns
------- -------
integrals : list exceeding_points : list
List of calculated integrals for each point. A list of points where the integral exceeds the local mean by the threshold.
local_means : list
List of local mean values (adjacent integrals).
fig : matplotlib.figure.Figure
The created figure object with the power plot.
ax : matplotlib.axes.Axes
The axes object for further modifications.
""" """
exceeding_points = []
fig, ax = plt.subplots()
ax.plot(frequency, power) # Plot power spectrum
integrals = []
local_means = []
for point in points: for point in points:
# Define indices for the integration window # Calculate the integral and local mean for the current point
indices = (frequency >= point - delta) & (frequency <= point + delta) integral, local_mean = calculate_integral(frequency, power, point, delta)
# Calculate integral around the point
integral = np.trapz(power[indices], frequency[indices])
integrals.append(integral)
# Calculate adjacent region integrals for local mean # Check if the integral exceeds the threshold
left_indices = (frequency >= point - 5 * delta) & (frequency < point - delta) valid, message = valid_integrals(integral, local_mean, threshold, point)
right_indices = (frequency > point + delta) & (frequency <= point + 5 * delta)
l_integral = np.trapz(power[left_indices], frequency[left_indices]) if valid:
r_integral = np.trapz(power[right_indices], frequency[right_indices]) exceeding_points.append(point)
local_mean = np.mean([l_integral, r_integral])
local_means.append(local_mean)
ax.set_xlim([0, 1200]) # Set x-axis limit
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Power')
ax.set_title('Power Spectrum with Integrals (Uncolored)')
return integrals, local_means, fig, ax return exceeding_points
def highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig_orig, ax_orig): def plot_highlighted_integrals(frequency, power, exceeding_points, delta, threshold, color_mapping, points_categories):
""" """
Create a new figure by highlighting integrals that exceed the threshold. Plot the power spectrum and highlight integrals that exceed the threshold.
This function generates a new figure with colored shading around points where the integrals exceed
the local mean by a given threshold and adds vertical lines at the boundaries of adjacent regions.
It leaves the original figure unchanged.
Parameters Parameters
---------- ----------
@ -176,75 +232,62 @@ def highlight_integrals_with_threshold(frequency, power, points, delta, threshol
An array of frequencies corresponding to the power values. An array of frequencies corresponding to the power values.
power : np.array power : np.array
An array of power spectral density values. An array of power spectral density values.
points : list exceeding_points : list
A list of harmonic frequencies to highlight. A list of harmonic frequencies that exceed the threshold.
delta : float delta : float
Half-width of the range for integration around each point. Half-width of the range for integration around each point.
threshold : float threshold : float
Threshold value to compare integrals with local mean. Threshold value to compare integrals with local mean.
integrals : list
List of calculated integrals for each point.
local_means : list
List of local mean values (adjacent integrals).
color_mapping : dict color_mapping : dict
A mapping of point categories to colors. A dictionary mapping each category to its color.
points_categories : dict points_categories : dict
A mapping of categories to lists of points. A mapping of categories to lists of points.
fig_orig : matplotlib.figure.Figure
The original figure object (remains unchanged).
ax_orig : matplotlib.axes.Axes
The original axes object (remains unchanged).
Returns Returns
------- -------
fig_new : matplotlib.figure.Figure fig : matplotlib.figure.Figure
The new figure object with color highlights and vertical lines. The created figure object with highlighted integrals.
""" """
fig, ax = plt.subplots()
ax.plot(frequency, power) # Plot power spectrum
# Create a new figure based on the original power spectrum for point in exceeding_points:
fig_new, ax_new = plt.subplots() integral, local_mean = calculate_integral(frequency, power, point, delta)
ax_new.plot(frequency, power) # Plot the same power spectrum valid, _ = valid_integrals(integral, local_mean, threshold, point)
if valid:
# Loop through each point and check if the integral exceeds the threshold
for i, point in enumerate(points):
exceeds = integrals[i] > (local_means[i] * threshold)
if exceeds:
# Define color based on the category of the point # 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') 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 # Shade the region around the point where the integral was calculated
ax_new.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz') ax.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz')
print(f"Integral around {point:.2f} Hz: {integrals[i]:.5e}") print(f"Integral around {point:.2f} Hz: {integral:.5e}")
# Define left and right boundaries of adjacent regions # Define left and right boundaries of adjacent regions
left_boundary = frequency[np.where((frequency >= point - 5 * delta) & (frequency < point - delta))[0][0]] 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]] 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 # Add vertical dashed lines at the boundaries of the adjacent regions
ax_new.axvline(x=left_boundary, color="k", linestyle="--") ax.axvline(x=left_boundary, color="k", linestyle="--")
ax_new.axvline(x=right_boundary, color="k", linestyle="--") ax.axvline(x=right_boundary, color="k", linestyle="--")
# Update plot legend and return the new figure
ax_new.set_xlim([0, 1200])
ax_new.set_xlabel('Frequency (Hz)')
ax_new.set_ylabel('Power')
ax_new.set_title('Power Spectrum with Highlighted Integrals')
ax_new.legend()
return fig_new
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 ### ### Data retrieval ###
datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner 'data' datafolder = "../data"
example_file = os.path.join("..", "data", "2024-10-16-ad-invivo-1.nix") example_file = os.path.join("..", "data", "2024-10-16-ad-invivo-1.nix")
dataset = rlx.Dataset(example_file) dataset = rlx.Dataset(example_file)
sams = dataset.repro_runs("SAM") sams = dataset.repro_runs("SAM")
sam = sams[2] sam = sams[2]
## Daten für Funktionen ## Data for functions
df = sam.metadata["RePro-Info"]["settings"]["deltaf"][0][0] df = sam.metadata["RePro-Info"]["settings"]["deltaf"][0][0]
stim = sam.stimuli[1] stim = sam.stimuli[1]
potential, time = stim.trace_data("V-1") potential, time = stim.trace_data("V-1")
@ -253,29 +296,34 @@ duration = stim.duration
dt = stim.trace_info("V-1").sampling_interval dt = stim.trace_info("V-1").sampling_interval
### Anwendung Functionen ### ### Apply Functions to calculate data ###
b = binary_spikes(spikes, duration, dt) b = binary_spikes(spikes, duration, dt)
rate = firing_rate(b, box_width=0.05, dt=dt) rate = firing_rate(b, box_width=0.05, dt=dt)
frequency, power = powerspectrum(b, dt) frequency, power = powerspectrum(b, dt)
## Important stuff
### Important stuff ###
## Frequencies
eodf = stim.metadata[stim.name]["EODf"][0][0] eodf = stim.metadata[stim.name]["EODf"][0][0]
stimulus_frequency = eodf + df stimulus_frequency = eodf + df
AM = 50 # Hz AM = 50 # Hz
#print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}")
frequencies = [AM, eodf, stimulus_frequency] frequencies = [AM, eodf, stimulus_frequency]
categories = ["AM", "EODf", "Stimulus frequency"] categories = ["AM", "EODf", "Stimulus frequency"]
num_harmonics = [4, 2, 2] num_harmonics = [4, 2, 2]
colors = ["green", "orange", "red"] colors = ["green", "orange", "red"]
delta = 2.5 delta = 2.5
threshold = 10 threshold = 10
### ### Apply functions to make powerspectrum ###
points, color_mapping, points_categories = prepare_harmonics(frequencies, categories, num_harmonics, colors) integral, local = calculate_integral(frequency, power, eodf, delta)
valid = valid_integrals(integral, local, threshold, eodf)
# First, create the power spectrum plot with integrals (without coloring) points, color, categories = prepare_harmonics(frequencies, categories, num_harmonics, colors)
integrals, local_means, fig1, ax1 = plot_power_spectrum_with_integrals(frequency, power, points, delta) print(len(points))
exceeding = find_exceeding_points(frequency, power, points, delta, threshold)
# Then, create a new separate figure where integrals exceeding the threshold are highlighted print(len(exceeding))
fig2 = highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig1, ax1)
## Plot power spectrum and highlight integrals
fig = plot_highlighted_integrals(frequency, power, points, delta, threshold, color, categories)
plt.show()