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"])