From 337b9d8e0e23404b27fa21581934e5aded144460 Mon Sep 17 00:00:00 2001 From: "sarah.eisele" Date: Fri, 25 Oct 2024 13:39:39 +0000 Subject: [PATCH] Delete code/GP_Code.py --- code/GP_Code.py | 329 ------------------------------------------------ 1 file changed, 329 deletions(-) delete mode 100644 code/GP_Code.py diff --git a/code/GP_Code.py b/code/GP_Code.py deleted file mode 100644 index 13d87b3..0000000 --- a/code/GP_Code.py +++ /dev/null @@ -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()