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

Experimenteller Aufbau.py

@@ -0,0 +1,42 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Oct 18 11:07:01 2024
+
+@author: diana
+"""
+### df relative to EODf outside of anesthesia
+    # relative range: -400 zu 2000 Hz
+    # Stim.zeit: 2s
+    # 25 Hz Schritte
+    # Intensität: 20%/ 10%
+# Daten die wir analysieren
+    # relative firing rate über stim. frequenz mit baseline als horizontale linie
+        # Firing rate
+
+
+
+### To Dos:
+    # Power spektrum anotieren
+    # AM plot machen 
+    # 
+    
+    
+
+
+
+
+# Verworfen
+### Stimulus Intensity 
+    # relative range zur outside EODf: -600 bis 1000Hz
+    # größere Schritte (50 Hz)
+    # kürzere Stimulationszeiten
+    # Stimulationsintensitäten als log Skala (20, 10, 5, 1, 0.5, 0.1 %)
+# Daten 
+    ## extract power spectrum for phase locking
+    # wo muss phase locking liegen
+    # search window mit normierten Daten 
+    # signal-to-noise ratio: wo ist peak (z-score)
+    # herausfinden, ob organ auf jeden peak einer bestimmten frequenz feuert
+        # bei welcher stim. frequenz gibt es ein 1:1 firing
+    # für jede zelle schauen, ab welcher Intensität eine Antwort bei einer 
+    # bestimmten frequenz zu sehen ist

analysis_1.py

@@ -0,0 +1,161 @@
+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 vt_spikes(dataset):
+    # 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')
+
+    # plot
+    fig = plt.figure(figsize=(5, 2.5))
+    # alternative to ax = axs[0]
+    ax = fig.add_subplot()
+    # plot vt diagram
+    ax.plot(time[time<0.1], potential[time<0.1])
+    # plot spikes into vt diagram, at max V
+    ax.scatter(spike_times[spike_times<0.1], np.ones_like(spike_times[spike_times<0.1]) * np.max(potential))
+    plt.show()
+    
+    return sam
+
+def scatter_plot(sam1):
+
+    ### plot scatter plot for one sam with all 3 stims
+    # get stim count
+    stim_count = sam1.stimulus_count
+    
+    # create colormap
+    colors = plt.cm.prism(np.linspace(0, 1, stim_count))
+    
+    # plot
+    fig = plt.figure()
+    ax = fig.add_subplot()
+    
+    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)
+        
+    ax.eventplot(stimuli, colors=colors)
+    ax.set_xlabel('Spike Times [ms]')
+    ax.set_ylabel('Loop #')
+    ax.set_yticks(range(stim_count))
+    ax.set_title('Spikes of SAM 3')
+    plt.show()
+    return stim, stim_count, time_stim
+
+# 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
+    
+
+
+'''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 metadata
+dataset = rlx.Dataset(example_file)
+
+### plot
+# timeline of whole rec
+dataset.plot_timeline()
+
+# voltage and spikes of current sam
+sam = vt_spikes(dataset)
+
+# spike times of all loops
+stim, stim_count, time_stim = scatter_plot(sam)
+
+
+'''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)
+
+# 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()
+
+eodf = stim.metadata[stim.name]['EODF'][0][0]
+df = stim.metadata['RePro-Info']['settings']['deltaf'][0][0]
+stimulus_freq = df + eodf
+
+
+### 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
+    # clean up current code (define variables outside of functions, plot spectrum in function)
+    # git
+    # tuning curve over stim intensities or over delta f?

code/GP_Code.py

@@ -1,10 +1,10 @@
 # -*- coding: utf-8 -*-
 """
-Created on Thu Oct 17 09:23:10 2024
+Created on Tue Oct 22 11:43:41 2024
 
 @author: diana
 """
-# -*- coding: utf-8 -*-
+
 
 import glob
 import os
@@ -101,11 +101,12 @@ def prepare_harmonics(frequencies, categories, num_harmonics, colors):
     return points, color_mapping, points_categories
 
 
-def plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories):
+def plot_power_spectrum_with_integrals(frequency, power, points, delta):
-    """Create a figure of the power spectrum with integrals highlighted around specified points.
+    """
+    Create a figure of the power spectrum and calculate integrals around specified points.
 
-    This function creates a plot of the power spectrum and shades areas around
+    This function generates the plot of the power spectrum and calculates integrals
-    specified harmonic points to indicate the calculated integrals.
+    around specified harmonic points, but it does not color the regions or add vertical lines.
 
     Parameters
     ----------
@@ -117,37 +118,121 @@ def plot_power_spectrum_with_integrals(frequency, power, points, delta, color_ma
         A list of harmonic frequencies to highlight.
     delta : float
         Half-width of the range for integration around each point.
-    color_mapping : dict
-        A mapping of point categories to colors.
-    points_categories : dict
-        A mapping of categories to lists of points.
 
     Returns
     -------
+    integrals : list
+        List of calculated integrals for each point.
+    local_means : list
+        List of local mean values (adjacent integrals).
     fig : matplotlib.figure.Figure
-        The created figure object.
+        The created figure object with the power plot.
+    ax : matplotlib.axes.Axes
+        The axes object for further modifications.
     """
+    
     fig, ax = plt.subplots()
-    ax.plot(frequency, power)
+    ax.plot(frequency, power)  # Plot power spectrum
+    
     integrals = []
+    local_means = []
 
     for point in points:
+        # Define indices for the integration window
         indices = (frequency >= point - delta) & (frequency <= point + delta)
+        # Calculate integral around the point
         integral = np.trapz(power[indices], frequency[indices])
         integrals.append(integral)
-
+        
-        # Get color based on point category
+        # Calculate adjacent region integrals for local mean
-        color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
+        left_indices = (frequency >= point - 5 * delta) & (frequency < point - delta)
-        ax.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz')
+        right_indices = (frequency > point + delta) & (frequency <= point + 5 * delta)
-        print(f"Integral around {point:.2f} Hz: {integral:.5e}")
+        
-
+        l_integral = np.trapz(power[left_indices], frequency[left_indices])
-    ax.set_xlim([0, 1200])
+        r_integral = np.trapz(power[right_indices], frequency[right_indices])
+        
+        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 marked Integrals')
+    ax.set_title('Power Spectrum with Integrals (Uncolored)')
-    ax.legend()
+
+    return integrals, local_means, fig, ax
+
+
+def highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig_orig, ax_orig):
+    """
+    Create a new figure by highlighting 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
+    ----------
+    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 highlight.
+    delta : float
+        Half-width of the range for integration around each point.
+    threshold : float
+        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
+        A mapping of point categories to colors.
+    points_categories : dict
+        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
+    -------
+    fig_new : matplotlib.figure.Figure
+        The new figure object with color highlights and vertical lines.
+    """
+    
+    # Create a new figure based on the original power spectrum
+    fig_new, ax_new = plt.subplots()
+    ax_new.plot(frequency, power)  # Plot the same power spectrum
+    
+    # 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
+            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_new.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}")
+            
+            # 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_new.axvline(x=left_boundary, color="k", linestyle="--")
+            ax_new.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
 
-    return fig
 
 
 
@@ -184,9 +269,13 @@ categories = ["AM", "EODf", "Stimulus frequency"]
 num_harmonics = [4, 2, 2]
 colors = ["green", "orange", "red"]
 delta = 2.5
+threshold = 10
 
-
+### 
-### Peaks im Powerspektrum finden ###
 points, color_mapping, points_categories = prepare_harmonics(frequencies, categories, num_harmonics, colors)
-fig = plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories)
+
-plt.show()
+# First, create the power spectrum plot with integrals (without coloring)
+integrals, local_means, fig1, ax1 = plot_power_spectrum_with_integrals(frequency, power, points, delta)
+
+# Then, create a new separate figure where integrals exceeding the threshold are highlighted
+fig2 = highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig1, ax1)

code/analysis_1.py

@@ -151,4 +151,8 @@ 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
+    # we want to see peaks at phase locking to own and stim frequency, and at amp modulation frequency
+    
+'''TEST- HI'''
+
+