From dde8bd842dc7303a42bda2af4d37529d22ed8693 Mon Sep 17 00:00:00 2001 From: Jan Grewe Date: Thu, 17 Sep 2020 11:48:39 +0200 Subject: [PATCH] extract despine method progress with response plot --- chirp_ams.py | 17 +------- punit_responses.py | 7 +-- response_discriminability.py | 82 ++++++++++++++++++++++++++++++++---- util.py | 16 +++++++ 4 files changed, 94 insertions(+), 28 deletions(-) diff --git a/chirp_ams.py b/chirp_ams.py index cb9f6a1..1ff3170 100644 --- a/chirp_ams.py +++ b/chirp_ams.py @@ -3,22 +3,7 @@ import scipy.signal as sig import matplotlib.pyplot as plt from chirp_stimulation import create_chirp -from IPython import embed - -def despine(axis, spines=None, hide_ticks=True): - - def hide_spine(spine): - spine.set_visible(False) - - for spine in axis.spines.keys(): - if spines is not None: - if spine in spines: - hide_spine(axis.spines[spine]) - else: - hide_spine(axis.spines[spine]) - if hide_ticks: - axis.xaxis.set_ticks([]) - axis.yaxis.set_ticks([]) +from util import despine def get_signals(eodfs, condition, contrast, chirp_size, chirp_duration, chirp_amplitude_dip, diff --git a/punit_responses.py b/punit_responses.py index 93cb7b4..cf1c3c8 100644 --- a/punit_responses.py +++ b/punit_responses.py @@ -39,7 +39,7 @@ def save(filename, name, stimulus_settings, model_settings, self_signal, other_s b.metadata = mdata # save stimulus - stim_da = b.create_data_array("complete_stimulus", "nix.timeseries.sampled", dtype=nix.DataType.Float, + stim_da = b.create_data_array("complete stimulus", "nix.timeseries.sampled.stimulus", dtype=nix.DataType.Float, data=complete_stimulus) stim_da.label = "voltage" stim_da.label = "mV/cm" @@ -49,7 +49,7 @@ def save(filename, name, stimulus_settings, model_settings, self_signal, other_s self_freq_da = None if self_freq is not None: - self_freq_da = b.create_data_array("self frequency", "nix.timeseries.sampled", dtype=nix.DataType.Float, + self_freq_da = b.create_data_array("self frequency", "nix.timeseries.sampled.frequency", dtype=nix.DataType.Float, data=self_freq) self_freq_da.label = "frequency" self_freq_da.label = "Hz" @@ -169,7 +169,8 @@ def simulate_responses(stimulus_params, model_params, repeats=10, deltaf=20): if condition == "self": v_0 = np.random.rand(1)[0] cell_params["v_zero"] = v_0 - no_other_spikes.append(simulate(np.hstack((pre_stim, self_signal)), **cell_params)) + sp = simulate(np.hstack((pre_stim, self_signal)), **cell_params) + no_other_spikes.append(sp[sp > pre_time[-1]] - pre_time[-1]) if condition == "self": name = "contrast_%.3f_condition_no-other_deltaf_%i" %(contrast, deltaf) save(filename, name, params, cell_params, self_signal, None, self_freq, None, self_signal, no_other_spikes) diff --git a/response_discriminability.py b/response_discriminability.py index 069e4eb..96eed60 100644 --- a/response_discriminability.py +++ b/response_discriminability.py @@ -4,7 +4,7 @@ import nixio as nix import numpy as np import matplotlib.pyplot as plt -from util import firing_rate +from util import firing_rate, despine from IPython import embed @@ -42,8 +42,9 @@ def sort_blocks(nix_file): def get_firing_rate(block_map, df, contrast, condition, kernel_width=0.0005): block = block_map[(contrast, df, condition)] - print((contrast, df, condition)) + print(block.name) response_map = {} + spikes = [] for da in block.data_arrays: if "spike_times" in da.type and "response" in da.name: resp_id = int(da.name.split("_")[-1]) @@ -53,18 +54,81 @@ def get_firing_rate(block_map, df, contrast, condition, kernel_width=0.0005): time = np.arange(0.0, duration, dt) rates = np.zeros((len(response_map.keys()), len(time))) for i, k in enumerate(response_map.keys()): - rates[i,:] = firing_rate(response_map[k][:], duration, kernel_width, dt) - return time, rates + spikes.append(response_map[k][:]) + rates[i,:] = firing_rate(spikes[-1], duration, kernel_width, dt) + return time, rates, spikes + + +def get_signals(block): + print(block.name) + self_freq = None + other_freq = None + signal = None + time = None + if "complete stimulus" not in block.data_arrays or "self frequency" not in block.data_arrays: + raise ValueError("Signals not stored in block!") + if "no-other" not in block.name and "other frequency" not in block.data_arrays: + raise ValueError("Signals not stored in block!") + + signal = block.data_arrays["complete stimulus"][:] + time = np.asarray(block.data_arrays["complete stimulus"].dimensions[0].axis(len(signal))) + self_freq = block.data_arrays["self frequency"][:] + if "no-other" not in block.name: + other_freq = block.data_arrays["other frequency"][:] + return signal, self_freq, other_freq, time def create_response_plot(block_map, all_dfs, all_contrasts, all_conditions, current_df): conditions = ["no-other", "self", "other"] condition_labels = ["alone", "self", "other"] + max_time = 0.5 + + fig = plt.figure(figsize=(6.5, 5.5)) + fig_grid = (len(all_contrasts) + 1, len(all_conditions)*3+2) + all_contrasts = sorted(all_contrasts, reverse=True) + + for i, condition in enumerate(conditions): + # plot the signals + block = block_map[(all_contrasts[0], current_df, condition)] + _, self_freq, other_freq, time = get_signals(block) + self_eodf = block.metadata["stimulus parameter"]["eodfs"]["self"] + other_eodf = block.metadata["stimulus parameter"]["eodfs"]["other"] + + ax = plt.subplot2grid(fig_grid, (0, i * 3 + i), rowspan=1, colspan=3, fig=fig) + ax.plot(time[time < max_time], self_freq[time < max_time], color="#ff7f0e", label="%iHz" % self_eodf) + ax.text(-0.05, self_eodf, "%iHz" % self_eodf, color="#ff7f0e", va="center", ha="right", fontsize=9) + if other_freq is not None: + ax.plot(time[time < max_time], other_freq[time < max_time], color="#1f77b4", label="%iHz" % other_eodf) + ax.text(-0.05, other_eodf, "%iHz" % other_eodf, color="#1f77b4", va="center", ha="right", fontsize=9) + ax.set_title(condition_labels[i]) + despine(ax, ["top", "bottom", "left", "right"], True) + + # for the largest contrast plot the raster with psth, only a section of the data (e.g. 1s) + t, rates, spikes = get_firing_rate(block_map, current_df, all_contrasts[0], condition, kernel_width=0.001) + avg_resp = np.mean(rates, axis=0) + error = np.std(rates, axis=0) + + ax = plt.subplot2grid(fig_grid, (1, i * 3 + i), rowspan=1, colspan=3) + ax.plot(t[t < max_time], avg_resp[t < max_time], color="k", lw=0.5) + ax.fill_between(t[t < max_time], (avg_resp - error)[t < max_time], (avg_resp + error)[t < max_time], color="k", lw=None, alpha=0.25) + despine(ax, ["top", "right"], False) + + # for all other contrast plot the firing rate alone + for j in range(1, len(all_contrasts)): + contrast = all_contrasts[j] + t, rates, _ = get_firing_rate(block_map, current_df, contrast, condition) + avg_resp = np.mean(rates, axis=0) + error = np.std(rates, axis=0) + ax = plt.subplot2grid(fig_grid, (j+1, i * 3 + i), rowspan=1, colspan=3) + ax.plot(t[t < max_time], avg_resp[t < max_time], color="k", lw=0.5) + #ax.fill_between(t[t < max_time], (avg_resp - error)[t < max_time], (avg_resp + error)[t < max_time], color="k", lw=None, alpha=0.25) + despine(ax, ["top", "right"], False) + + plt.savefig("chirp_responses.pdf") + plt.close() + return + - fig = plt.figure(figsize=(4.5, 5.5)) - - - embed() def process_cell(filename, dfs=[], contrasts=[], conditions=[]): @@ -75,7 +139,7 @@ def process_cell(filename, dfs=[], contrasts=[], conditions=[]): else: print("ERROR: no baseline data for file %s!" % filename) - + create_response_plot(block_map, all_contrasts, all_dfs, all_conditions, 20) """ if len(dfs) == 0: dfs = all_dfs diff --git a/util.py b/util.py index 12d3e80..2832dc7 100644 --- a/util.py +++ b/util.py @@ -1,5 +1,21 @@ import numpy as np +def despine(axis, spines=None, hide_ticks=True): + + def hide_spine(spine): + spine.set_visible(False) + + for spine in axis.spines.keys(): + if spines is not None: + if spine in spines: + hide_spine(axis.spines[spine]) + else: + hide_spine(axis.spines[spine]) + if hide_ticks: + axis.xaxis.set_ticks([]) + axis.yaxis.set_ticks([]) + + def gaussKernel(sigma, dt): """ Creates a Gaussian kernel with a given standard deviation and an integral of 1.