From 54d10789b4bee89bbe967c3013760586435aa17b Mon Sep 17 00:00:00 2001 From: mbergmann Date: Tue, 22 Oct 2024 07:21:37 +0000 Subject: [PATCH] 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