From 7e673267e5c25095ea08ee55896ec01c8ed24df2 Mon Sep 17 00:00:00 2001 From: Jan Grewe Date: Fri, 18 Sep 2020 16:18:33 +0200 Subject: [PATCH] finished response figure ... move data and figs to folders, ignore them in git --- .gitignore | 2 + chirp_ams.py | 5 +- response_discriminability.py | 131 ++++++++++++++++++++++++----------- 3 files changed, 96 insertions(+), 42 deletions(-) diff --git a/.gitignore b/.gitignore index f1bdbf6..bcc47ed 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ .vscode/.ropeproject/config.py .vscode/.ropeproject/objectdb +figures +data diff --git a/chirp_ams.py b/chirp_ams.py index 1ff3170..dbd6c94 100644 --- a/chirp_ams.py +++ b/chirp_ams.py @@ -1,3 +1,4 @@ +import os import numpy as np import scipy.signal as sig import matplotlib.pyplot as plt @@ -5,6 +6,8 @@ import matplotlib.pyplot as plt from chirp_stimulation import create_chirp from util import despine +figure_folder = "figures" + def get_signals(eodfs, condition, contrast, chirp_size, chirp_duration, chirp_amplitude_dip, chirp_times, duration, dt): @@ -107,5 +110,5 @@ if __name__ == "__main__": fig.subplots_adjust(left=0.1, bottom=0.1, top=0.99, right=0.99) - plt.savefig("Chirp_induced_ams.pdf") + plt.savefig(os.path.join(figure_folder, "Chirp_induced_AMs.pdf")) plt.close() \ No newline at end of file diff --git a/response_discriminability.py b/response_discriminability.py index 96eed60..3a4fb66 100644 --- a/response_discriminability.py +++ b/response_discriminability.py @@ -2,11 +2,14 @@ import os import glob import nixio as nix import numpy as np +import scipy.signal as sig import matplotlib.pyplot as plt +from IPython import embed from util import firing_rate, despine +figure_folder = "figures" +data_folder = "data" -from IPython import embed def read_baseline(block): spikes = [] @@ -42,7 +45,6 @@ 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(block.name) response_map = {} spikes = [] for da in block.data_arrays: @@ -60,7 +62,6 @@ def get_firing_rate(block_map, df, contrast, condition, kernel_width=0.0005): def get_signals(block): - print(block.name) self_freq = None other_freq = None signal = None @@ -78,39 +79,79 @@ def get_signals(block): return signal, self_freq, other_freq, time -def create_response_plot(block_map, all_dfs, all_contrasts, all_conditions, current_df): +def extract_am(signal): + # first add some padding + front_pad = np.flip(signal[:int(len(signal)/100)]) + back_pad = np.flip(signal[-int(len(signal)/100):]) + padded = np.hstack((front_pad, signal, back_pad)) + # do the hilbert and take abs + am = np.abs(sig.hilbert(padded)) + am = am[len(front_pad):-len(back_pad)] + return am + + +def create_response_plot(block_map, all_dfs, all_contrasts, all_conditions, current_df, figure_name=None): conditions = ["no-other", "self", "other"] - condition_labels = ["alone", "self", "other"] - max_time = 0.5 + condition_labels = ["soliloquy", "self chirping", "other chirping"] + min_time = 0.5 + max_time = min_time + 0.5 fig = plt.figure(figsize=(6.5, 5.5)) - fig_grid = (len(all_contrasts) + 1, len(all_conditions)*3+2) + fig_grid = (len(all_contrasts)*2 + 6, 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) + signal, self_freq, other_freq, time = get_signals(block) + am = extract_am(signal) + 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) + + # plot frequency traces + ax = plt.subplot2grid(fig_grid, (0, i * 3 + i), rowspan=2, colspan=3, fig=fig) + ax.plot(time[(time > min_time) & (time < max_time)], self_freq[(time > min_time) & (time < max_time)], + color="#ff7f0e", label="%iHz" % self_eodf) + ax.text(min_time-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.plot(time[(time > min_time) & (time < max_time)], other_freq[(time > min_time) & (time < max_time)], + color="#1f77b4", label="%iHz" % other_eodf) + ax.text(min_time-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) + # plot the am + ax = plt.subplot2grid(fig_grid, (3, i * 3 + i), rowspan=2, colspan=3, fig=fig) + ax.plot(time[(time > min_time) & (time < max_time)], signal[(time > min_time) & (time < max_time)], + color="#2ca02c", label="signal") + ax.plot(time[(time > min_time) & (time < max_time)], am[(time > min_time) & (time < max_time)], + color="#d62728", label="am") + despine(ax, ["top", "bottom", "left", "right"], True) + ax.set_ylim([-1.25, 1.25]) + ax.legend(ncol=2, loc=(0.01, -0.5), fontsize=7, markerscale=0.5, frameon=False) + # 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) + ax = plt.subplot2grid(fig_grid, (6, i * 3 + i), rowspan=2, colspan=3) + ax.plot(t[(t > min_time) & (t < max_time)], avg_resp[(t > min_time) & (t < max_time)], + color="k", lw=0.5) + ax.fill_between(t[(t > min_time) & (t < max_time)], (avg_resp - error)[(t > min_time) & (t < max_time)], + (avg_resp + error)[(t > min_time) & (t < max_time)], color="k", lw=0.0, alpha=0.25) + ax.set_ylim([0, 750]) + ax.set_xlabel("") + ax.set_ylabel("") + ax.set_xticks(np.arange(min_time, max_time+.01, 0.250)) + ax.set_xticklabels(map(int, (np.arange(min_time, max_time + .01, 0.250) - min_time) * 1000)) + ax.set_xticks(np.arange(min_time, max_time+.01, 0.0625), minor=True) + ax.set_xticklabels([]) + ax.set_yticks(np.arange(0.0, 751., 500)) + ax.set_yticks(np.arange(0.0, 751., 125), minor=True) + if i > 0: + ax.set_yticklabels([]) despine(ax, ["top", "right"], False) # for all other contrast plot the firing rate alone @@ -119,46 +160,54 @@ def create_response_plot(block_map, all_dfs, all_contrasts, all_conditions, curr 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) + ax = plt.subplot2grid(fig_grid, (j*2 + 6, i * 3 + i), rowspan=2, colspan=3) + ax.plot(t[(t > min_time) & (t < max_time)], avg_resp[(t > min_time) & (t < max_time)], color="k", lw=0.5) + ax.fill_between(t[(t > min_time) & (t < max_time)], (avg_resp - error)[(t > min_time) & (t < max_time)], + (avg_resp + error)[(t > min_time) & (t < max_time)], color="k", lw=0.0, alpha=0.25) + ax.set_ylim([0, 750]) + ax.set_xlabel("") + ax.set_ylabel("") + ax.set_xticks(np.arange(min_time, max_time+.01, 0.250)) + ax.set_xticklabels(map(int, (np.arange(min_time, max_time + .01, 0.250) - min_time) * 1000)) + ax.set_xticks(np.arange(min_time, max_time+.01, 0.125), minor=True) + if j < len(all_contrasts) -1: + ax.set_xticklabels([]) + ax.set_yticks(np.arange(0.0, 751., 500)) + ax.set_yticks(np.arange(0.0, 751., 125), minor=True) + if i > 0: + ax.set_yticklabels([]) despine(ax, ["top", "right"], False) - - plt.savefig("chirp_responses.pdf") + if i == 1: + ax.set_xlabel("time [ms]") + if i == 0: + ax.set_ylabel("frequency [Hz]", va="center") + ax.yaxis.set_label_coords(-0.45, 3.5) + + name = figure_name if figure_name is not None else "chirp_responses.pdf" + name = (name + ".pdf") if ".pdf" not in name else name + plt.savefig(os.path.join(figure_folder, name)) plt.close() - return - - def process_cell(filename, dfs=[], contrasts=[], conditions=[]): nf = nix.File.open(filename, nix.FileMode.ReadOnly) - block_map, all_dfs, all_contrasts, all_conditions = sort_blocks(nf) + block_map, all_contrasts, all_dfs, all_conditions = sort_blocks(nf) if "baseline" in block_map.keys(): baseline_spikes = read_baseline(block_map["baseline"]) 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 - if len(contrasts) == 0: - contrasts = all_contrasts - if len(conditions) == 0: - conditions = all_conditions + fig_name = filename.split(".nix")[0] + "_df_20Hz.pdf" + create_response_plot(block_map, all_dfs, all_contrasts, all_conditions, 20, figure_name=fig_name) + fig_name = filename.split(".nix")[0] + "_df_-100Hz.pdf" + create_response_plot(block_map, all_dfs, all_contrasts, all_conditions, -100, figure_name=fig_name) + - for df in dfs: - for condition in conditions: - for contrast in contrasts: - time, rates = get_firing_rate(block_map, df, contrast, condition, kernel_width=0.0025) - """ nf.close() def main(): - nix_files = sorted(glob.glob("cell*.nix")) + nix_files = sorted(glob.glob(os.path.join(data_folder, "cell*.nix"))) for nix_file in nix_files: process_cell(nix_file, dfs=[20], contrasts=[20], conditions=["self"])