Merge branch 'main' of https://whale.am28.uni-tuebingen.de/git/mbergmann/gpgrewe2024

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()

code/analysis_1.py

@@ -1,154 +0,0 @@
-import rlxnix as rlx
-import numpy as np
-import matplotlib.pyplot as plt
-import os
-from scipy.signal import welch
-
-# close all currently open figures
-plt.close('all')
-
-'''FUNCTIONS'''
-def plot_vt_spikes(t, v, spike_t):
-    fig = plt.figure(figsize=(5, 2.5))
-    # alternative to ax = axs[0]
-    ax = fig.add_subplot()
-    # plot vt diagram
-    ax.plot(t[t<0.1], v[t<0.1])
-    # plot spikes into vt diagram, at max V
-    ax.scatter(spike_t[spike_t<0.1], np.ones_like(spike_t[spike_t<0.1]) * np.max(v))
-    plt.show()
-
-def scatter_plot(colormap, stimuli_list, stimulus_count):
-    '''plot scatter plot for one sam with all 3 stims'''
-    fig = plt.figure()
-    ax = fig.add_subplot()
-        
-    ax.eventplot(stimuli_list, colors=colormap)
-    ax.set_xlabel('Spike Times [ms]')
-    ax.set_ylabel('Loop #')
-    ax.set_yticks(range(stimulus_count))
-    ax.set_title('Spikes of SAM 3')
-    plt.show()
-
-# create binary array with ones for spike times
-def binary_spikes(spike_times, duration , dt):
-    '''Converts spike times to binary representation
-    Params
-    ------
-    spike_times: np.array
-        spike times
-    duration: float
-        trial duration
-    dt: float
-        temporal resolution
-    
-    Returns
-    --------
-    binary: np.array
-        The binary representation of the spike times
-    '''
-    binary = np.zeros(int(duration//dt)) # // is truncated division, returns number w/o decimals, same as np.round
-    spike_indices = np.asarray(np.round(spike_times//dt), dtype=int)
-    binary[spike_indices] = 1 
-    return binary
-
-# function to plot psth
-def firing_rates(binary_spikes, box_width=0.01, dt=0.000025):
-    box = np.ones(int(box_width // dt))
-    box /= np.sum(box * dt) # normalize box kernel w interal of 1
-    rate = np.convolve(binary_spikes, box, mode='same')
-    return rate
-
-def power_spectrum(rate, dt):
-    f, p = welch(rate, fs = 1./dt, nperseg=2**16, noverlap=2**15) 
-    # algorithm makes rounding mistakes, we want to calc many spectra and take mean of those 
-    # nperseg: length of segments in # datapoints
-    # noverlap: # datapoints that overlap in segments
-    return f, p
-    
-def power_spectrum_plot(f, p):
-    # plot power spectrum
-    fig = plt.figure()
-    ax = fig.add_subplot()
-    ax.plot(freq, power)
-    ax.set_xlabel('Frequency [Hz]')
-    ax.set_ylabel('Power [1/Hz]')
-    ax.set_xlim(0, 1000)
-    plt.show()
-
-'''IMPORT DATA'''
-datafolder = '../data' #./ wo ich gerade bin; ../ eine ebene höher; ../../ zwei ebenen höher
-
-example_file = os.path.join('..', 'data', '2024-10-16-ac-invivo-1.nix')
-
-'''EXTRACT DATA'''
-dataset = rlx.Dataset(example_file)
-
-# get sams
-sams = dataset.repro_runs('SAM')
-sam = sams[2]
-
-# get potetial over time (vt curve)
-potential, time = sam.trace_data('V-1')
-
-# get spike times
-spike_times, _ = sam.trace_data('Spikes-1')
-
-# get stim count
-stim_count = sam.stimulus_count
-
-# extract spike times of all 3 loops of current sam
-stimuli = []
-for i in range(stim_count):
-    # get stim i from sam
-    stim = sam.stimuli[i]
-    potential_stim, time_stim = stim.trace_data('V-1')
-    # get spike_times
-    spike_times_stim, _ = stim.trace_data('Spikes-1')
-    stimuli.append(spike_times_stim)
-    
-eodf = stim.metadata[stim.name]['EODF'][0][0]
-df = stim.metadata['RePro-Info']['settings']['deltaf'][0][0]
-stimulus_freq = df + eodf
-
-'''PLOT'''
-# create colormap
-colors = plt.cm.prism(np.linspace(0, 1, stim_count))
-
-# timeline of whole rec
-dataset.plot_timeline()
-
-# voltage and spikes of current sam
-plot_vt_spikes(time, potential, spike_times)
-
-# spike times of all loops
-scatter_plot(colors, stimuli, stim_count)
-
-
-'''POWER SPECTRUM'''
-# define variables for binary spikes function
-spikes, _ = stim.trace_data('Spikes-1')
-ti = stim.trace_info('V-1')
-dt = ti.sampling_interval
-duration = stim.duration
-
-### spectrum
-# vector with binary values for wholes length of stim
-binary = binary_spikes(spikes, duration, dt)
-
-# calculate firing rate 
-rate = firing_rates(binary, 0.01, dt) # box width of 10 ms
-
-# plot psth or whatever
-# plt.plot(time_stim, rate)
-# plt.show()
-
-freq, power = power_spectrum(binary, dt)
-
-power_spectrum_plot(freq, power)
-
-
-### TODO:
-    # then loop over sams/dfs, all stims, intensities
-    # when does stim start in eodf/ at which phase and how does that influence our signal --> alignment problem: egal wenn wir spectren haben
-    # we want to see peaks at phase locking to own and stim frequency, and at amp modulation frequency

code/plot_functions.py

@@ -71,13 +71,20 @@ def power_spectrum_plot(f, p):
 functions_path = r"C:\Users\diana\OneDrive - UT Cloud\Master\GPs\GP1_Grewe\Projekt\gpgrewe2024\code"
 sys.path.append(functions_path)
 import useful_functions as u
+import matplotlib.ticker as ticker
 
-def plot_highlighted_integrals(frequency, power, points, color_mapping, points_categories, delta=2.5):
+def float_formatter(x, _):
+    """Format the y-axis values as floats with a specified precision."""
+    return f'{x:.5f}'
+
+def plot_highlighted_integrals(ax, frequency, power, points, color_mapping, points_categories, delta=2.5):
     """
-    Plot the power spectrum and highlight integrals that exceed the threshold.
+    Highlight integrals on the existing axes of the power spectrum.
 
     Parameters
     ----------
+    ax : matplotlib.axes.Axes
+        The axes on which to plot the highlighted integrals.
     frequency : np.array
         An array of frequencies corresponding to the power values.
     power : np.array
@@ -93,50 +100,40 @@ def plot_highlighted_integrals(frequency, power, points, color_mapping, points_c
 
     Returns
     -------
-    fig : matplotlib.figure.Figure
-        The created figure object with highlighted integrals.
+    None
     """
-    fig, ax = plt.subplots()
-    ax.plot(frequency, power)  # Plot power spectrum
-    
+    ax.plot(frequency, power, color = "k")  # Plot power spectrum on the existing axes
+
     for point in points:
-        # Use the imported function to calculate the integral and local mean
-        integral, local_mean, _ = u.calculate_integral(frequency, power, point)
+        # Calculate the integral and local mean
+        integral, local_mean = u.calculate_integral_2(frequency, power, point)
         
-        # Use the imported function to check if the point is valid
+        # Check if the point is valid
         valid = u.valid_integrals(integral, local_mean, 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')
-            
-            # Find the category of the point
             point_category = next((cat for cat, pts in points_categories.items() if point in pts), "Unknown")
+            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 out point, category, and color
-            print(f"{point_category}: Integral: {integral:.5e}, Color: {color}")
-            
-            # Annotate the plot with the point and its color
-            ax.text(point, max(power) * 0.9, f'{point:.2f}', color=color, fontsize=10, ha='center')
-
-            # 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.axvspan(point - delta, point + delta, color=color, alpha=0.2, label=f'{point_category}')
+
+            # Text with categories and colors
+            ax.text(1000, 5.8e-5, "AM", fontsize=10, color="green", alpha=0.2)
+            ax.text(1000, 5.6e-5, "Nyquist", fontsize=10, color="blue", alpha=0.2)
+            ax.text(1000, 5.4e-5, "EODf", fontsize=10, color="red", alpha=0.2)
+            ax.text(1000, 5.2e-5, "Stimulus frequency", fontsize=10, color="orange", alpha=0.2)
+            ax.text(1000, 5.0e-5, "EODf of awake fish", fontsize=10, color="purple", alpha=0.2)
     
     ax.set_xlim([0, 1200])
+    ax.set_ylim([0, 6e-5])
     ax.set_xlabel('Frequency (Hz)')
     ax.set_ylabel('Power')
-    ax.set_title('Power Spectrum with Highlighted Integrals')
-    ax.legend()
-    
-    return fig, ax
+    ax.set_title('Power Spectrum with highlighted Integrals')
+
+    # Apply float formatting to the y-axis
+    ax.yaxis.set_major_formatter(ticker.FuncFormatter(float_formatter))
 
 
 

code/useful_functions.py

@@ -47,7 +47,7 @@ def all_coming_together(freq_array, power_array, points_list, categories, num_ha
         color = colors[i]
         
         # Step 1: Calculate the integral for the point
-        integral, local_mean, _ = calculate_integral(freq_array, power_array, point, delta)
+        integral, local_mean = calculate_integral_2(freq_array, power_array, point, delta)
         
         # Step 2: Check if the point is valid
         valid = valid_integrals(integral, local_mean, point, threshold)
@@ -150,6 +150,42 @@ def calculate_integral(freq, power, point, delta = 2.5):
     local_mean = np.mean([l_integral, r_integral])
     return integral, local_mean, p_power
 
+def calculate_integral_2(freq, power, point, delta = 2.5):   
+    """
+    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, optional
+        Radius of the range for integration around the point. The default is 2.5.
+
+    Returns
+    -------
+    integral : float
+        The calculated integral around the point.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    p_power : float
+        The local maxiumum power.
+    """
+    indices = (freq >= point - delta) & (freq <= point + delta)
+    integral = np.trapz(power[indices], freq[indices])
+            
+    left_indices = (freq >= point - 5 * delta) & (freq < point - delta)
+    right_indices = (freq > point + delta) & (freq <= point + 5 * delta)
+           
+    l_integral = np.trapz(power[left_indices], freq[left_indices])
+    r_integral = np.trapz(power[right_indices], freq[right_indices])
+           
+    local_mean = np.mean([l_integral, r_integral])
+    return integral, local_mean
+
 def contrast_sorting(sams, con_1 = 20, con_2 = 10, con_3 = 5, stim_count = 3, stim_dur = 2):
     '''
     sorts the sams into three contrasts