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 1/5] 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 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 2/5] 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 3/5] 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 4/5] 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 5/5] 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