From 6fb62aadfc24081ff2e98951cf8ce3ab96c377ae Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Fri, 18 Oct 2024 14:47:51 +0200
Subject: [PATCH 01/13] code adding test.py
---
code/test.py | 101 +++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 101 insertions(+)
create mode 100644 code/test.py
diff --git a/code/test.py b/code/test.py
new file mode 100644
index 0000000..8f818b7
--- /dev/null
+++ b/code/test.py
@@ -0,0 +1,101 @@
+import glob
+import pathlib
+import numpy as np
+import matplotlib.pyplot as plt
+import rlxnix as rlx
+from IPython import embed
+from scipy.signal import welch
+
+def binary_spikes(spike_times, duration, dt):
+ """
+ Converts the spike times to a binary representations
+
+ 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 train.
+
+ """
+ binary = np.zeros(int(np.round(duration / dt))) #create the binary array with the same length as potential
+
+ spike_indices = np.asarray(np.round(spike_times / dt), dtype = int) # get the indices
+ binary[spike_indices] = 1 # put the indices into binary
+ return binary
+
+def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
+ box = np.ones(int(box_width // dt))
+ box /= np.sum(box) * dt # normalisierung des box kernels to an integral of one
+ rate = np.convolve(binary_spikes, box, mode = 'same')
+ return rate
+
+def power_spectrum(rate, dt):
+ freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
+ return freq, power
+
+#find example data
+datafolder = "../data"
+
+example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
+
+#load dataset
+dataset = rlx.Dataset(example_file)
+# find all sams
+sams = dataset.repro_runs('SAM')
+sam = sams[2] # our example sam
+potential,time = sam.trace_data("V-1") #membrane potential
+spike_times, _ = sam.trace_data('Spikes-1') #spike times
+df = sam.metadata['RePro-Info']['settings']['deltaf'][0][0] #find df in metadata
+amp = sam.metadata['RePro-Info']['settings']['contrast'][0][0] * 100 #find amplitude in metadata
+
+#figure for a quick plot
+fig = plt.figure(figsize = (5, 2.5))
+ax = fig.add_subplot()
+ax.plot(time[time < 0.1], potential[time < 0.1]) # plot the membrane potential in 0.1s
+ax.scatter(spike_times[spike_times < 0.1],
+ np.ones_like(spike_times[spike_times < 0.1]) * np.max(potential)) #plot teh spike times on top
+plt.show()
+plt.close()
+# get all the stimuli
+stims = sam.stimuli
+# empty list for the spike times
+spikes = []
+#spikes2 = np.array(range(len(stims)))
+# loop over the stimuli
+for stim in stims:
+ # get the spike times
+ spike, _ = stim.trace_data('Spikes-1')
+ # append the first 100ms to spikes
+ spikes.append(spike[spike < 0.1])
+ # get stimulus duration
+ duration = stim.duration
+ ti = stim.trace_info("V-1")
+ dt = ti.sampling_interval # get the stimulus interval
+ bin_spikes = binary_spikes(spike, duration, dt) #binarize the spike_times
+ print(len(bin_spikes))
+ pot,tim= stim.trace_data("V-1") #membrane potential
+ rate = firing_rate(bin_spikes, dt = dt)
+ print(np.mean(rate))
+ fig, [ax1, ax2] = plt.subplots(1, 2,layout = 'constrained')
+ ax1.plot(tim,rate)
+ ax1.set_ylim(0,600)
+ ax1.set_xlim(0, 0.04)
+ freq, power = power_spectrum(rate, dt)
+ ax2.plot(freq,power)
+ ax2.set_xlim(0,1000)
+
+# make an eventplot
+fig = plt.figure(figsize = (5, 3), layout = 'constrained')
+ax = fig.add_subplot()
+ax.eventplot(spikes, linelength = 0.8)
+ax.set_xlabel('time [ms]')
+ax.set_ylabel('loop no.')
+plt.show()
\ No newline at end of file
From 342914a71944b4f19cc1a3ad2921f9fece68de81 Mon Sep 17 00:00:00 2001
From: "sarah.eisele" <sarah.eisele@student.uni-tuebingen.de>
Date: Fri, 18 Oct 2024 14:59:41 +0000
Subject: [PATCH 02/13] Add analysis_1
---
analysis_1 | 0
1 file changed, 0 insertions(+), 0 deletions(-)
create mode 100644 analysis_1
diff --git a/analysis_1 b/analysis_1
new file mode 100644
index 0000000..e69de29
From 7c32c6e4b86bad2fc4c8177d78902acd72f818c1 Mon Sep 17 00:00:00 2001
From: "sarah.eisele" <sarah.eisele@student.uni-tuebingen.de>
Date: Fri, 18 Oct 2024 15:03:44 +0000
Subject: [PATCH 03/13] Upload files to "/"
---
analysis_1.py | 154 ++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 154 insertions(+)
create mode 100644 analysis_1.py
diff --git a/analysis_1.py b/analysis_1.py
new file mode 100644
index 0000000..b6d6497
--- /dev/null
+++ b/analysis_1.py
@@ -0,0 +1,154 @@
+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
\ No newline at end of file
From 72a97f18ce31d9e7a5b41d7d50a3e0ca8cc8930d Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Fri, 18 Oct 2024 15:04:26 +0000
Subject: [PATCH 04/13] =?UTF-8?q?analysis=5F1=20gel=C3=B6scht?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
was empty and not needed
---
analysis_1 | 0
1 file changed, 0 insertions(+), 0 deletions(-)
delete mode 100644 analysis_1
diff --git a/analysis_1 b/analysis_1
deleted file mode 100644
index e69de29..0000000
From 052210c1eadff598467ef034b447b50b512d2f61 Mon Sep 17 00:00:00 2001
From: Diana <diana.stoll@student.uni-tuebingen.de>
Date: Mon, 21 Oct 2024 12:31:36 +0000
Subject: [PATCH 05/13] Dateien nach "/" hochladen
Integrale unter Peaks berechnen
---
GP_Code.py | 192 +++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 192 insertions(+)
create mode 100644 GP_Code.py
diff --git a/GP_Code.py b/GP_Code.py
new file mode 100644
index 0000000..6ff4428
--- /dev/null
+++ b/GP_Code.py
@@ -0,0 +1,192 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct 17 09:23:10 2024
+
+@author: diana
+"""
+# -*- coding: utf-8 -*-
+
+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.
+
+ This function computes the firing rate using a boxcar (moving average)
+ filter of a specified width.
+
+ 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 prepare_harmonics(frequencies, categories, num_harmonics, colors):
+ 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 plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories):
+ """Create a figure of the power spectrum with integrals highlighted around specified points.
+
+ This function creates a plot of the power spectrum and shades areas around
+ specified harmonic points to indicate the calculated integrals.
+
+ 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.
+ color_mapping : dict
+ A mapping of point categories to colors.
+ points_categories : dict
+ A mapping of categories to lists of points.
+
+ Returns
+ -------
+ fig : matplotlib.figure.Figure
+ The created figure object.
+ """
+ fig, ax = plt.subplots()
+ ax.plot(frequency, power)
+ integrals = []
+
+ for point in points:
+ indices = (frequency >= point - delta) & (frequency <= point + delta)
+ integral = np.trapz(power[indices], frequency[indices])
+ integrals.append(integral)
+
+ # Get color based on point category
+ color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
+ 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}")
+
+ ax.set_xlim([0, 1200])
+ ax.set_xlabel('Frequency (Hz)')
+ ax.set_ylabel('Power')
+ ax.set_title('Power Spectrum with marked Integrals')
+ ax.legend()
+
+ return fig
+
+
+
+
+### Data retrieval ###
+datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner '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]
+
+## Daten für Funktionen
+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
+
+
+### Anwendung Functionen ###
+b = binary_spikes(spikes, duration, dt)
+rate = firing_rate(b, box_width=0.05, dt=dt)
+frequency, power = powerspectrum(b, dt)
+
+## Important stuff
+eodf = stim.metadata[stim.name]["EODf"][0][0]
+stimulus_frequency = eodf + df
+AM = 50 # Hz
+#print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}")
+
+frequencies = [AM, eodf, stimulus_frequency]
+categories = ["AM", "EODf", "Stimulus frequency"]
+num_harmonics = [4, 2, 2]
+colors = ["green", "orange", "red"]
+delta = 2.5
+
+
+### 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()
From b0adb74ae33d9b2ad9ba5eed6af84977d9d1a35b Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 07:14:44 +0000
Subject: [PATCH 06/13] Datei test.py nach "code" hochladen
---
code/test.py | 101 +++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 101 insertions(+)
create mode 100644 code/test.py
diff --git a/code/test.py b/code/test.py
new file mode 100644
index 0000000..f2b9a4e
--- /dev/null
+++ b/code/test.py
@@ -0,0 +1,101 @@
+import glob
+import pathlib
+import numpy as np
+import matplotlib.pyplot as plt
+import rlxnix as rlx
+from IPython import embed
+from scipy.signal import welch
+
+def binary_spikes(spike_times, duration, dt):
+ """
+ Converts the spike times to a binary representations
+
+ 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 train.
+
+ """
+ binary = np.zeros(int(np.round(duration / dt))) #create the binary array with the same length as potential
+
+ spike_indices = np.asarray(np.round(spike_times / dt), dtype = int) # get the indices
+ binary[spike_indices] = 1 # put the indices into binary
+ return binary
+
+def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
+ box = np.ones(int(box_width // dt))
+ box /= np.sum(box) * dt # normalisierung des box kernels to an integral of one
+ rate = np.convolve(binary_spikes, box, mode = 'same')
+ return rate
+
+def power_spectrum(rate, dt):
+ freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
+ return freq, power
+
+#find example data
+datafolder = "../data"
+
+example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
+
+#load dataset
+dataset = rlx.Dataset(example_file)
+# find all sams
+sams = dataset.repro_runs('SAM')
+sam = sams[2] # our example sam
+potential,time = sam.trace_data("V-1") #membrane potential
+spike_times, _ = sam.trace_data('Spikes-1') #spike times
+df = sam.metadata['RePro-Info']['settings']['deltaf'][0][0] #find df in metadata
+amp = sam.metadata['RePro-Info']['settings']['contrast'][0][0] * 100 #find amplitude in metadata
+
+#figure for a quick plot
+fig = plt.figure(figsize = (5, 2.5))
+ax = fig.add_subplot()
+ax.plot(time[time < 0.1], potential[time < 0.1]) # plot the membrane potential in 0.1s
+ax.scatter(spike_times[spike_times < 0.1],
+ np.ones_like(spike_times[spike_times < 0.1]) * np.max(potential)) #plot teh spike times on top
+plt.show()
+plt.close()
+# get all the stimuli
+stims = sam.stimuli
+# empty list for the spike times
+spikes = []
+#spikes2 = np.array(range(len(stims)))
+# loop over the stimuli
+for stim in stims:
+ # get the spike times
+ spike, _ = stim.trace_data('Spikes-1')
+ # append the first 100ms to spikes
+ spikes.append(spike[spike < 0.1])
+ # get stimulus duration
+ duration = stim.duration
+ ti = stim.trace_info("V-1")
+ dt = ti.sampling_interval # get the stimulus interval
+ bin_spikes = binary_spikes(spike, duration, dt) #binarize the spike_times
+ print(len(bin_spikes))
+ pot,tim= stim.trace_data("V-1") #membrane potential
+ rate = firing_rate(bin_spikes, dt = dt)
+ print(np.mean(rate))
+ fig, [ax1, ax2] = plt.subplots(1, 2,layout = 'constrained')
+ ax1.plot(tim,rate)
+ ax1.set_ylim(0,600)
+ ax1.set_xlim(0, 0.04)
+ freq, power = power_spectrum(rate, dt)
+ ax2.plot(freq,power)
+ ax2.set_xlim(0,1000)
+
+# make an eventplot
+fig = plt.figure(figsize = (5, 3), layout = 'constrained')
+ax = fig.add_subplot()
+ax.eventplot(spikes, linelength = 0.8)
+ax.set_xlabel('time [ms]')
+ax.set_ylabel('loop no.')
+plt.show()
\ No newline at end of file
From e8d4ca90fe7c13fed6fab75c2042d64801b48994 Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 07:16:25 +0000
Subject: [PATCH 07/13] GP_code nach code verschoben
---
GP_Code.py => code/GP_Code.py | 384 +++++++++++++++++-----------------
1 file changed, 192 insertions(+), 192 deletions(-)
rename GP_Code.py => code/GP_Code.py (96%)
diff --git a/GP_Code.py b/code/GP_Code.py
similarity index 96%
rename from GP_Code.py
rename to code/GP_Code.py
index 6ff4428..754715e 100644
--- a/GP_Code.py
+++ b/code/GP_Code.py
@@ -1,192 +1,192 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Oct 17 09:23:10 2024
-
-@author: diana
-"""
-# -*- coding: utf-8 -*-
-
-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.
-
- This function computes the firing rate using a boxcar (moving average)
- filter of a specified width.
-
- 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 prepare_harmonics(frequencies, categories, num_harmonics, colors):
- 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 plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories):
- """Create a figure of the power spectrum with integrals highlighted around specified points.
-
- This function creates a plot of the power spectrum and shades areas around
- specified harmonic points to indicate the calculated integrals.
-
- 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.
- color_mapping : dict
- A mapping of point categories to colors.
- points_categories : dict
- A mapping of categories to lists of points.
-
- Returns
- -------
- fig : matplotlib.figure.Figure
- The created figure object.
- """
- fig, ax = plt.subplots()
- ax.plot(frequency, power)
- integrals = []
-
- for point in points:
- indices = (frequency >= point - delta) & (frequency <= point + delta)
- integral = np.trapz(power[indices], frequency[indices])
- integrals.append(integral)
-
- # Get color based on point category
- color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
- 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}")
-
- ax.set_xlim([0, 1200])
- ax.set_xlabel('Frequency (Hz)')
- ax.set_ylabel('Power')
- ax.set_title('Power Spectrum with marked Integrals')
- ax.legend()
-
- return fig
-
-
-
-
-### Data retrieval ###
-datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner '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]
-
-## Daten für Funktionen
-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
-
-
-### Anwendung Functionen ###
-b = binary_spikes(spikes, duration, dt)
-rate = firing_rate(b, box_width=0.05, dt=dt)
-frequency, power = powerspectrum(b, dt)
-
-## Important stuff
-eodf = stim.metadata[stim.name]["EODf"][0][0]
-stimulus_frequency = eodf + df
-AM = 50 # Hz
-#print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}")
-
-frequencies = [AM, eodf, stimulus_frequency]
-categories = ["AM", "EODf", "Stimulus frequency"]
-num_harmonics = [4, 2, 2]
-colors = ["green", "orange", "red"]
-delta = 2.5
-
-
-### 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()
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct 17 09:23:10 2024
+
+@author: diana
+"""
+# -*- coding: utf-8 -*-
+
+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.
+
+ This function computes the firing rate using a boxcar (moving average)
+ filter of a specified width.
+
+ 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 prepare_harmonics(frequencies, categories, num_harmonics, colors):
+ 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 plot_power_spectrum_with_integrals(frequency, power, points, delta, color_mapping, points_categories):
+ """Create a figure of the power spectrum with integrals highlighted around specified points.
+
+ This function creates a plot of the power spectrum and shades areas around
+ specified harmonic points to indicate the calculated integrals.
+
+ 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.
+ color_mapping : dict
+ A mapping of point categories to colors.
+ points_categories : dict
+ A mapping of categories to lists of points.
+
+ Returns
+ -------
+ fig : matplotlib.figure.Figure
+ The created figure object.
+ """
+ fig, ax = plt.subplots()
+ ax.plot(frequency, power)
+ integrals = []
+
+ for point in points:
+ indices = (frequency >= point - delta) & (frequency <= point + delta)
+ integral = np.trapz(power[indices], frequency[indices])
+ integrals.append(integral)
+
+ # Get color based on point category
+ color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
+ 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}")
+
+ ax.set_xlim([0, 1200])
+ ax.set_xlabel('Frequency (Hz)')
+ ax.set_ylabel('Power')
+ ax.set_title('Power Spectrum with marked Integrals')
+ ax.legend()
+
+ return fig
+
+
+
+
+### Data retrieval ###
+datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner '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]
+
+## Daten für Funktionen
+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
+
+
+### Anwendung Functionen ###
+b = binary_spikes(spikes, duration, dt)
+rate = firing_rate(b, box_width=0.05, dt=dt)
+frequency, power = powerspectrum(b, dt)
+
+## Important stuff
+eodf = stim.metadata[stim.name]["EODf"][0][0]
+stimulus_frequency = eodf + df
+AM = 50 # Hz
+#print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}")
+
+frequencies = [AM, eodf, stimulus_frequency]
+categories = ["AM", "EODf", "Stimulus frequency"]
+num_harmonics = [4, 2, 2]
+colors = ["green", "orange", "red"]
+delta = 2.5
+
+
+### 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()
From 7e02490a89d17a96689c3f0c0ad4919df2e09b93 Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 07:17:10 +0000
Subject: [PATCH 08/13] analysis_1.py nach code verschoben
---
analysis_1.py => code/analysis_1.py | 306 ++++++++++++++--------------
1 file changed, 153 insertions(+), 153 deletions(-)
rename analysis_1.py => code/analysis_1.py (95%)
diff --git a/analysis_1.py b/code/analysis_1.py
similarity index 95%
rename from analysis_1.py
rename to code/analysis_1.py
index b6d6497..17fb46b 100644
--- a/analysis_1.py
+++ b/code/analysis_1.py
@@ -1,154 +1,154 @@
-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
+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
\ No newline at end of file
From 54d10789b4bee89bbe967c3013760586435aa17b Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 07:21:37 +0000
Subject: [PATCH 09/13] test.py aktualisiert
---
code/test.py | 230 +++++++++++++++++++++++++++++----------------------
1 file changed, 130 insertions(+), 100 deletions(-)
diff --git a/code/test.py b/code/test.py
index f2b9a4e..9274e87 100644
--- a/code/test.py
+++ b/code/test.py
@@ -1,101 +1,131 @@
-import glob
-import pathlib
-import numpy as np
-import matplotlib.pyplot as plt
-import rlxnix as rlx
-from IPython import embed
-from scipy.signal import welch
-
-def binary_spikes(spike_times, duration, dt):
- """
- Converts the spike times to a binary representations
-
- 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 train.
-
- """
- binary = np.zeros(int(np.round(duration / dt))) #create the binary array with the same length as potential
-
- spike_indices = np.asarray(np.round(spike_times / dt), dtype = int) # get the indices
- binary[spike_indices] = 1 # put the indices into binary
- return binary
-
-def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
- box = np.ones(int(box_width // dt))
- box /= np.sum(box) * dt # normalisierung des box kernels to an integral of one
- rate = np.convolve(binary_spikes, box, mode = 'same')
- return rate
-
-def power_spectrum(rate, dt):
- freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
- return freq, power
-
-#find example data
-datafolder = "../data"
-
-example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
-
-#load dataset
-dataset = rlx.Dataset(example_file)
-# find all sams
-sams = dataset.repro_runs('SAM')
-sam = sams[2] # our example sam
-potential,time = sam.trace_data("V-1") #membrane potential
-spike_times, _ = sam.trace_data('Spikes-1') #spike times
-df = sam.metadata['RePro-Info']['settings']['deltaf'][0][0] #find df in metadata
-amp = sam.metadata['RePro-Info']['settings']['contrast'][0][0] * 100 #find amplitude in metadata
-
-#figure for a quick plot
-fig = plt.figure(figsize = (5, 2.5))
-ax = fig.add_subplot()
-ax.plot(time[time < 0.1], potential[time < 0.1]) # plot the membrane potential in 0.1s
-ax.scatter(spike_times[spike_times < 0.1],
- np.ones_like(spike_times[spike_times < 0.1]) * np.max(potential)) #plot teh spike times on top
-plt.show()
-plt.close()
-# get all the stimuli
-stims = sam.stimuli
-# empty list for the spike times
-spikes = []
-#spikes2 = np.array(range(len(stims)))
-# loop over the stimuli
-for stim in stims:
- # get the spike times
- spike, _ = stim.trace_data('Spikes-1')
- # append the first 100ms to spikes
- spikes.append(spike[spike < 0.1])
- # get stimulus duration
- duration = stim.duration
- ti = stim.trace_info("V-1")
- dt = ti.sampling_interval # get the stimulus interval
- bin_spikes = binary_spikes(spike, duration, dt) #binarize the spike_times
- print(len(bin_spikes))
- pot,tim= stim.trace_data("V-1") #membrane potential
- rate = firing_rate(bin_spikes, dt = dt)
- print(np.mean(rate))
- fig, [ax1, ax2] = plt.subplots(1, 2,layout = 'constrained')
- ax1.plot(tim,rate)
- ax1.set_ylim(0,600)
- ax1.set_xlim(0, 0.04)
- freq, power = power_spectrum(rate, dt)
- ax2.plot(freq,power)
- ax2.set_xlim(0,1000)
-
-# make an eventplot
-fig = plt.figure(figsize = (5, 3), layout = 'constrained')
-ax = fig.add_subplot()
-ax.eventplot(spikes, linelength = 0.8)
-ax.set_xlabel('time [ms]')
-ax.set_ylabel('loop no.')
+import glob
+import pathlib
+import numpy as np
+import matplotlib.pyplot as plt
+import rlxnix as rlx
+from IPython import embed
+from scipy.signal import welch
+
+def binary_spikes(spike_times, duration, dt):
+ """
+ Converts the spike times to a binary representations
+
+ 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 train.
+
+ """
+ binary = np.zeros(int(np.round(duration / dt))) #create the binary array with the same length as potential
+
+ spike_indices = np.asarray(np.round(spike_times / dt), dtype = int) # get the indices
+ binary[spike_indices] = 1 # put the indices into binary
+ return binary
+
+def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
+ box = np.ones(int(box_width // dt))
+ box /= np.sum(box) * dt # normalisierung des box kernels to an integral of one
+ rate = np.convolve(binary_spikes, box, mode = 'same')
+ return rate
+
+def power_spectrum(rate, dt):
+ freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
+ return freq, power
+
+def extract_stim_data(stimulus):
+ '''
+ extracts all necessary metadata for each stimulus
+
+ Parameters
+ ----------
+ stimulus : Stimulus object or rlxnix.base.repro module
+ The stimulus from which the data is needed.
+
+ Returns
+ -------
+ amplitude : float
+ The relative signal amplitude in percent.
+ df : float
+ Distance of the stimulus to the current EODf.
+ eodf : float
+ Current EODf.
+ stim_freq : float
+ The total stimulus frequency (EODF+df).
+
+ '''
+ # extract metadata
+ # the stim.name adjusts the first key as it changes with every stimulus
+ amplitude = stim.metadata[stim.name]['Contrast'][0][0]
+ df = stim.metadata[stim.name]['DeltaF'][0][0]
+ eodf = stim.metadata[stim.name]['EODf'][0][0]
+ stim_freq = stim.metadata[stim.name]['Frequency'][0][0]
+ return amplitude, df, eodf, stim_freq
+
+
+#find example data
+datafolder = "../data"
+
+example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
+
+#load dataset
+dataset = rlx.Dataset(example_file)
+# find all sams
+sams = dataset.repro_runs('SAM')
+sam = sams[2] # our example sam
+potential,time = sam.trace_data("V-1") #membrane potential
+spike_times, _ = sam.trace_data('Spikes-1') #spike times
+df = sam.metadata['RePro-Info']['settings']['deltaf'][0][0] #find df in metadata
+amp = sam.metadata['RePro-Info']['settings']['contrast'][0][0] * 100 #find amplitude in metadata
+
+#figure for a quick plot
+fig = plt.figure(figsize = (5, 2.5))
+ax = fig.add_subplot()
+ax.plot(time[time < 0.1], potential[time < 0.1]) # plot the membrane potential in 0.1s
+ax.scatter(spike_times[spike_times < 0.1],
+ np.ones_like(spike_times[spike_times < 0.1]) * np.max(potential)) #plot teh spike times on top
+plt.show()
+plt.close()
+# get all the stimuli
+stims = sam.stimuli
+# empty list for the spike times
+spikes = []
+#spikes2 = np.array(range(len(stims)))
+# loop over the stimuli
+for stim in stims:
+ # get the spike times
+ spike, _ = stim.trace_data('Spikes-1')
+ # append the first 100ms to spikes
+ spikes.append(spike[spike < 0.1])
+ # get stimulus duration
+ duration = stim.duration
+ ti = stim.trace_info("V-1")
+ dt = ti.sampling_interval # get the stimulus interval
+ bin_spikes = binary_spikes(spike, duration, dt) #binarize the spike_times
+ print(len(bin_spikes))
+ pot,tim= stim.trace_data("V-1") #membrane potential
+ rate = firing_rate(bin_spikes, dt = dt)
+ print(np.mean(rate))
+ fig, [ax1, ax2] = plt.subplots(1, 2,layout = 'constrained')
+ ax1.plot(tim,rate)
+ ax1.set_ylim(0,600)
+ ax1.set_xlim(0, 0.04)
+ freq, power = power_spectrum(rate, dt)
+ ax2.plot(freq,power)
+ ax2.set_xlim(0,1000)
+
+# make an eventplot
+fig = plt.figure(figsize = (5, 3), layout = 'constrained')
+ax = fig.add_subplot()
+ax.eventplot(spikes, linelength = 0.8)
+ax.set_xlabel('time [ms]')
+ax.set_ylabel('loop no.')
plt.show()
\ No newline at end of file
From 51ce9d1178ca0c98ded68176b6ad631f1011a692 Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 09:45:27 +0200
Subject: [PATCH 10/13] updated test.py
---
code/test.py | 107 ---------------------------------------------------
1 file changed, 107 deletions(-)
diff --git a/code/test.py b/code/test.py
index a0eddc9..32b704e 100644
--- a/code/test.py
+++ b/code/test.py
@@ -1,4 +1,3 @@
-
import glob
import pathlib
import numpy as np
@@ -42,7 +41,6 @@ def power_spectrum(rate, dt):
freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
return freq, power
-
def extract_stim_data(stimulus):
'''
extracts all necessary metadata for each stimulus
@@ -72,110 +70,6 @@ def extract_stim_data(stimulus):
stim_freq = stim.metadata[stim.name]['Frequency'][0][0]
return amplitude, df, eodf, stim_freq
-
-
-#find example data
-datafolder = "../data"
-
-example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
-
-#load dataset
-dataset = rlx.Dataset(example_file)
-# find all sams
-sams = dataset.repro_runs('SAM')
-sam = sams[2] # our example sam
-potential,time = sam.trace_data("V-1") #membrane potential
-spike_times, _ = sam.trace_data('Spikes-1') #spike times
-df = sam.metadata['RePro-Info']['settings']['deltaf'][0][0] #find df in metadata
-amp = sam.metadata['RePro-Info']['settings']['contrast'][0][0] * 100 #find amplitude in metadata
-
-#figure for a quick plot
-fig = plt.figure(figsize = (5, 2.5))
-ax = fig.add_subplot()
-ax.plot(time[time < 0.1], potential[time < 0.1]) # plot the membrane potential in 0.1s
-ax.scatter(spike_times[spike_times < 0.1],
- np.ones_like(spike_times[spike_times < 0.1]) * np.max(potential)) #plot teh spike times on top
-plt.show()
-plt.close()
-# get all the stimuli
-stims = sam.stimuli
-# empty list for the spike times
-spikes = []
-#spikes2 = np.array(range(len(stims)))
-# loop over the stimuli
-for stim in stims:
- # get the spike times
- spike, _ = stim.trace_data('Spikes-1')
- # append the first 100ms to spikes
- spikes.append(spike[spike < 0.1])
- # get stimulus duration
- duration = stim.duration
- ti = stim.trace_info("V-1")
- dt = ti.sampling_interval # get the stimulus interval
- bin_spikes = binary_spikes(spike, duration, dt) #binarize the spike_times
- print(len(bin_spikes))
- pot,tim= stim.trace_data("V-1") #membrane potential
- rate = firing_rate(bin_spikes, dt = dt)
- print(np.mean(rate))
- fig, [ax1, ax2] = plt.subplots(1, 2,layout = 'constrained')
- ax1.plot(tim,rate)
- ax1.set_ylim(0,600)
- ax1.set_xlim(0, 0.04)
- freq, power = power_spectrum(rate, dt)
- ax2.plot(freq,power)
- ax2.set_xlim(0,1000)
-
-# make an eventplot
-fig = plt.figure(figsize = (5, 3), layout = 'constrained')
-ax = fig.add_subplot()
-ax.eventplot(spikes, linelength = 0.8)
-ax.set_xlabel('time [ms]')
-ax.set_ylabel('loop no.')
-<<<<<<< HEAD
-=======
-import glob
-import pathlib
-import numpy as np
-import matplotlib.pyplot as plt
-import rlxnix as rlx
-from IPython import embed
-from scipy.signal import welch
-
-def binary_spikes(spike_times, duration, dt):
- """
- Converts the spike times to a binary representations
-
- 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 train.
-
- """
- binary = np.zeros(int(np.round(duration / dt))) #create the binary array with the same length as potential
-
- spike_indices = np.asarray(np.round(spike_times / dt), dtype = int) # get the indices
- binary[spike_indices] = 1 # put the indices into binary
- return binary
-
-def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
- box = np.ones(int(box_width // dt))
- box /= np.sum(box) * dt # normalisierung des box kernels to an integral of one
- rate = np.convolve(binary_spikes, box, mode = 'same')
- return rate
-
-def power_spectrum(rate, dt):
- freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
- return freq, power
-
#find example data
datafolder = "../data"
@@ -233,4 +127,3 @@ ax = fig.add_subplot()
ax.eventplot(spikes, linelength = 0.8)
ax.set_xlabel('time [ms]')
ax.set_ylabel('loop no.')
-plt.show()
\ No newline at end of file
From 2401bcbc16132c252e923fee00a6d696afa7d12a Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 10:15:32 +0200
Subject: [PATCH 11/13] updated test.py with new functions
---
code/test.py | 46 ++++++++++++++++++++++++++++++++++++++++++----
1 file changed, 42 insertions(+), 4 deletions(-)
diff --git a/code/test.py b/code/test.py
index 32b704e..b407995 100644
--- a/code/test.py
+++ b/code/test.py
@@ -60,18 +60,52 @@ def extract_stim_data(stimulus):
Current EODf.
stim_freq : float
The total stimulus frequency (EODF+df).
+ amp_mod : float
+ The current amplitude modulation.
+ ny_freq : float
+ The current nyquist frequency.
'''
# extract metadata
# the stim.name adjusts the first key as it changes with every stimulus
amplitude = stim.metadata[stim.name]['Contrast'][0][0]
df = stim.metadata[stim.name]['DeltaF'][0][0]
- eodf = stim.metadata[stim.name]['EODf'][0][0]
- stim_freq = stim.metadata[stim.name]['Frequency'][0][0]
- return amplitude, df, eodf, stim_freq
+ eodf = round(stim.metadata[stim.name]['EODf'][0][0])
+ stim_freq = round(stim.metadata[stim.name]['Frequency'][0][0])
+ # calculates the amplitude modulation
+ amp_mod, ny_freq = AM(eodf, stim_freq)
+ return amplitude, df, eodf, stim_freq, amp_mod, ny_freq
+
+def AM(EODf, stimulus):
+ """
+ Calculates the Amplitude Modulation and Nyquist frequency
+
+ Parameters
+ ----------
+ EODf : float or int
+ The current EODf.
+ stimulus : float or int
+ The absolute frequency of the stimulus.
+
+ Returns
+ -------
+ AM : float
+ The amplitude modulation resulting from the stimulus.
+ nyquist : float
+ The maximum frequency possible to resolve with the EODf.
+
+ """
+ nyquist = EODf * 0.5
+ AM = np.mod(stimulus, nyquist)
+ return AM, nyquist
+
+def remove_poor(files):
+ good_files =files
+ print('x')
+ return good_files
#find example data
-datafolder = "../data"
+datafolder = "../../data"
example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
@@ -120,6 +154,10 @@ for stim in stims:
freq, power = power_spectrum(rate, dt)
ax2.plot(freq,power)
ax2.set_xlim(0,1000)
+ plt.close()
+ if stim == stims[-1]:
+ amplitude, df, eodf, stim_freq = extract_stim_data(stim)
+ print(amplitude, df, eodf, stim_freq)
# make an eventplot
fig = plt.figure(figsize = (5, 3), layout = 'constrained')
From 8ccd633f859fe6164952c20ffbc5b04bfd8b272b Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 10:46:46 +0200
Subject: [PATCH 12/13] updated test.py with new functions
---
code/test.py | 33 ++++++++++++++++++++++++++++++---
1 file changed, 30 insertions(+), 3 deletions(-)
diff --git a/code/test.py b/code/test.py
index b407995..fe621f7 100644
--- a/code/test.py
+++ b/code/test.py
@@ -100,14 +100,41 @@ def AM(EODf, stimulus):
return AM, nyquist
def remove_poor(files):
- good_files =files
- print('x')
+ """
+ Removes poor datasets from the set of files for analysis
+
+ Parameters
+ ----------
+ files : list
+ list of files.
+
+ Returns
+ -------
+ good_files : list
+ list of files without the ones with the label poor.
+
+ """
+ # create list for good files
+ good_files = []
+ # loop over files
+ for i in range(len(files)):
+ # print(files[i])
+ # load the file (takes some time)
+ data = rlx.Dataset(files[i])
+ # get the quality
+ quality = str.lower(data.metadata["Recording"]["Recording quality"][0][0])
+ # check the quality
+ if quality != "poor":
+ # if its good or fair add it to the good files
+ good_files.append(files[i])
return good_files
#find example data
datafolder = "../../data"
-example_file = datafolder + "/" + "2024-10-16-ad-invivo-1.nix"
+example_file = datafolder + "/" + "2024-10-16-ah-invivo-1.nix"
+
+data_files = glob.glob("../../data/*.nix")
#load dataset
dataset = rlx.Dataset(example_file)
From 2f5a1d27546d89fce6ca1e03ca4d1ed9be17b6cb Mon Sep 17 00:00:00 2001
From: mbergmann <maximilian.bergmann@student.uni-tuebingen.de>
Date: Tue, 22 Oct 2024 11:28:12 +0200
Subject: [PATCH 13/13] update test.py & new: useful_functions.py
---
code/useful_functions.py | 173 +++++++++++++++++++++++++++++++++++++++
1 file changed, 173 insertions(+)
create mode 100644 code/useful_functions.py
diff --git a/code/useful_functions.py b/code/useful_functions.py
new file mode 100644
index 0000000..1115ef5
--- /dev/null
+++ b/code/useful_functions.py
@@ -0,0 +1,173 @@
+import glob
+import pathlib
+import numpy as np
+import matplotlib.pyplot as plt
+import rlxnix as rlx
+from IPython import embed
+from scipy.signal import welch
+
+def AM(EODf, stimulus):
+ """
+ Calculates the Amplitude Modulation and Nyquist frequency
+
+ Parameters
+ ----------
+ EODf : float or int
+ The current EODf.
+ stimulus : float or int
+ The absolute frequency of the stimulus.
+
+ Returns
+ -------
+ AM : float
+ The amplitude modulation resulting from the stimulus.
+ nyquist : float
+ The maximum frequency possible to resolve with the EODf.
+
+ """
+ nyquist = EODf * 0.5
+ AM = np.mod(stimulus, nyquist)
+ return AM, nyquist
+
+def binary_spikes(spike_times, duration, dt):
+ """
+ Converts the spike times to a binary representations
+
+ 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 train.
+
+ """
+ binary = np.zeros(int(np.round(duration / dt))) #create the binary array with the same length as potential
+
+ spike_indices = np.asarray(np.round(spike_times / dt), dtype = int) # get the indices
+ binary[spike_indices] = 1 # put the indices into binary
+ return binary
+
+def extract_stim_data(stimulus):
+ '''
+ extracts all necessary metadata for each stimulus
+
+ Parameters
+ ----------
+ stimulus : Stimulus object or rlxnix.base.repro module
+ The stimulus from which the data is needed.
+
+ Returns
+ -------
+ amplitude : float
+ The relative signal amplitude in percent.
+ df : float
+ Distance of the stimulus to the current EODf.
+ eodf : float
+ Current EODf.
+ stim_freq : float
+ The total stimulus frequency (EODF+df).
+ amp_mod : float
+ The current amplitude modulation.
+ ny_freq : float
+ The current nyquist frequency.
+
+ '''
+ # extract metadata
+ # the stim.name adjusts the first key as it changes with every stimulus
+ amplitude = stimulus.metadata[stimulus.name]['Contrast'][0][0]
+ df = stimulus.metadata[stimulus.name]['DeltaF'][0][0]
+ eodf = round(stimulus.metadata[stimulus.name]['EODf'][0][0])
+ stim_freq = round(stimulus.metadata[stimulus.name]['Frequency'][0][0])
+ # calculates the amplitude modulation
+ amp_mod, ny_freq = AM(eodf, stim_freq)
+ return amplitude, df, eodf, stim_freq, amp_mod, ny_freq
+
+def firing_rate(binary_spikes, dt = 0.000025, box_width = 0.01):
+ '''
+ Calculates the firing rate from binary spikes
+
+ Parameters
+ ----------
+ binary_spikes : np.array
+ The binary representation of the spike train.
+ dt : float, optional
+ Time difference between two datapoints. The default is 0.000025.
+ box_width : float, optional
+ Time window on which the rate should be computed on. The default is 0.01.
+
+ Returns
+ -------
+ rate : np.array
+ Array of firing rates.
+
+ '''
+ box = np.ones(int(box_width // dt))
+ box /= np.sum(box) * dt # normalisierung des box kernels to an integral of one
+ rate = np.convolve(binary_spikes, box, mode = 'same')
+ return rate
+
+def power_spectrum(spike_times, duration, dt):
+ '''
+ Computes a power spectrum based on the spike times
+
+ Parameters
+ ----------
+ spike_times : np.array
+ The spike times.
+ duration : float
+ The trial duration:
+ dt : float
+ The temporal resolution.
+
+ Returns
+ -------
+ freq : np.array
+ All the frequencies of the power spectrum.
+ power : np.array
+ Power of the frequencies calculated.
+
+ '''
+ # binarizes spikes
+ binary = binary_spikes(spike_times, duration, dt)
+ # computes firing rates
+ rate = firing_rate(binary, dt = dt)
+ # creates power spectrum
+ freq, power = welch(rate, fs = 1/dt, nperseg = 2**16, noverlap = 2**15)
+ return freq, power
+
+def remove_poor(files):
+ """
+ Removes poor datasets from the set of files for analysis
+
+ Parameters
+ ----------
+ files : list
+ list of files.
+
+ Returns
+ -------
+ good_files : list
+ list of files without the ones with the label poor.
+
+ """
+ # create list for good files
+ good_files = []
+ # loop over files
+ for i in range(len(files)):
+ # print(files[i])
+ # load the file (takes some time)
+ data = rlx.Dataset(files[i])
+ # get the quality
+ quality = str.lower(data.metadata["Recording"]["Recording quality"][0][0])
+ # check the quality
+ if quality != "poor":
+ # if its good or fair add it to the good files
+ good_files.append(files[i])
+ return good_files
\ No newline at end of file