From c8427269507f677eae528687d2ddfd7fbe6ba4ea Mon Sep 17 00:00:00 2001 From: Jan Grewe Date: Tue, 27 Oct 2015 12:24:33 +0100 Subject: [PATCH] python script to create firing rates in different ways --- programming/code/firing_rates.py | 198 +++++++++++++++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100644 programming/code/firing_rates.py diff --git a/programming/code/firing_rates.py b/programming/code/firing_rates.py new file mode 100644 index 0000000..5d5d8f9 --- /dev/null +++ b/programming/code/firing_rates.py @@ -0,0 +1,198 @@ +import numpy as np +import matplotlib.pyplot as plt +import scipy.io as spio +import scipy as sp +import seaborn as sb +from IPython import embed +sb.set_context("paper") + + +def set_axis_fontsize(axis, label_size, tick_label_size=None, legend_size=None): + """ + Sets axis, tick label and legend font sizes to the desired size. + + :param axis: the axes object + :param label_size: the size of the axis label + :param tick_label_size: the size of the tick labels. If None, lable_size is used + :param legend_size: the size of the font used in the legend.If None, the label_size is used. + + """ + if not tick_label_size: + tick_label_size = label_size + if not legend_size: + legend_size = label_size + axis.xaxis.get_label().set_fontsize(label_size) + axis.yaxis.get_label().set_fontsize(label_size) + for tick in axis.xaxis.get_major_ticks() + axis.yaxis.get_major_ticks(): + tick.label.set_fontsize(tick_label_size) + + l = axis.get_legend() + if l: + for t in l.get_texts(): + t.set_fontsize(legend_size) + + +def get_instantaneous_rate(times, max_t=30., dt=1e-4): + time = np.arange(0., max_t, dt) + indices = np.asarray(times / dt, dtype=int) + intervals = np.diff(np.hstack(([0], times))) + inst_rate = np.zeros(time.shape) + + for i, index in enumerate(indices[1:]): + inst_rate[indices[i-1]:indices[i]] = 1/intervals[i] + return time, inst_rate + + +def plot_isi_rate(spike_times, max_t=30, dt=1e-4): + times = np.squeeze(spike_times[0][0])[:100000] + time, rate = get_instantaneous_rate(times, max_t=100000*dt) + + rates = np.zeros((len(rate), len(spike_times))) + for i in range(len(spike_times)): + _, rates[:, i] = get_instantaneous_rate(np.squeeze(spike_times[i][0])[:100000], + max_t=100000*dt) + avg_rate = np.mean(rates, axis=1) + rate_std = np.std(rates, axis=1) + + fig = plt.figure() + ax1 = fig.add_subplot(311) + ax2 = fig.add_subplot(312) + ax3 = fig.add_subplot(313) + + ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5) + ax1.set_ylabel("skpikes", fontsize=12) + set_axis_fontsize(ax1, 12) + + ax2.plot(time, rate, label="instantaneous rate, trial 1") + ax2.set_ylabel("firing rate [Hz]", fontsize=12) + ax2.legend(fontsize=12) + set_axis_fontsize(ax2, 12) + + ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue", + alpha=0.5, label="standard deviation") + ax3.plot(time, avg_rate, label="average rate") + ax3.set_xlabel("times [s]", fontsize=12) + ax3.set_ylabel("firing rate [Hz]", fontsize=12) + ax3.legend(fontsize=12) + ax3.set_ylim([0, 450]) + set_axis_fontsize(ax3, 12) + + fig.set_size_inches(15, 10) + fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95) + fig.set_facecolor("white") + fig.savefig("../lectures/images/instantaneous_rate.png") + plt.close() + + +def get_binned_rate(times, bin_width=0.05, max_t=30., dt=1e-4): + time = np.arange(0., max_t, dt) + bins = np.arange(0., max_t, bin_width) + bin_indices = bins / dt + hist, _ = sp.histogram(times, bins) + rate = np.zeros(time.shape) + + for i, b in enumerate(bin_indices[1:]): + rate[bin_indices[i-1]:b] = hist[i-1]/bin_width + return time, rate + + +def plot_bin_rate(spike_times, bin_width, max_t=30, dt=1e-4): + times = np.squeeze(spike_times[0][0]) + time, rate = get_binned_rate(times) + rates = np.zeros((len(rate), len(spike_times))) + for i in range(len(spike_times)): + _, rates[:, i] = get_binned_rate(np.squeeze(spike_times[i][0])) + avg_rate = np.mean(rates, axis=1) + rate_std = np.std(rates, axis=1) + + fig = plt.figure() + ax1 = fig.add_subplot(311) + ax2 = fig.add_subplot(312) + ax3 = fig.add_subplot(313) + + ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5) + ax1.set_ylabel("skpikes", fontsize=12) + ax1.set_xlim([0, 10]) + set_axis_fontsize(ax1, 12) + + ax2.plot(time, rate, label="binned rate, trial 1") + ax2.set_ylabel("firing rate [Hz]", fontsize=12) + ax2.legend(fontsize=12) + ax2.set_xlim([0, 10]) + set_axis_fontsize(ax2, 12) + + ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue", + alpha=0.5, label="standard deviation") + ax3.plot(time, avg_rate, label="average rate") + ax3.set_xlabel("times [s]", fontsize=12) + ax3.set_ylabel("firing rate [Hz]", fontsize=12) + ax3.legend(fontsize=12) + ax3.set_xlim([0, 10]) + ax3.set_ylim([0, 450]) + set_axis_fontsize(ax3, 12) + + fig.set_size_inches(15, 10) + fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95) + fig.set_facecolor("white") + fig.savefig("../lectures/images/binned_rate.png") + plt.close() + + +def get_convolved_rate(times, sigma, max_t=30., dt=1.e-4): + time = np.arange(0., max_t, dt) + kernel = sp.stats.norm.pdf(np.arange(-8*sigma, 8*sigma, dt),loc=0,scale=sigma) + indices = np.asarray(times/dt, dtype=int) + rate = np.zeros(time.shape) + rate[indices] = 1.; + conv_rate = np.convolve(rate, kernel, mode="same") + return time, conv_rate + + +def plot_conv_rate(spike_times, sigma=0.05, max_t=30, dt=1e-4): + times = np.squeeze(spike_times[0][0]) + time, rate = get_convolved_rate(times, sigma) + + rates = np.zeros((len(rate), len(spike_times))) + for i in range(len(spike_times)): + _, rates[:, i] = get_convolved_rate(np.squeeze(spike_times[i][0]), sigma) + avg_rate = np.mean(rates, axis=1) + rate_std = np.std(rates, axis=1) + + fig = plt.figure() + ax1 = fig.add_subplot(311) + ax2 = fig.add_subplot(312) + ax3 = fig.add_subplot(313) + + ax1.vlines(times[times < (100000*dt)], ymin=0, ymax=1, color="dodgerblue", lw=1.5) + ax1.set_ylabel("skpikes", fontsize=12) + ax1.set_xlim([0, 10]) + set_axis_fontsize(ax1, 12) + + ax2.plot(time, rate, label="convolved rate, trial 1") + ax2.set_ylabel("firing rate [Hz]", fontsize=12) + ax2.legend(fontsize=12) + ax2.set_xlim([0, 10]) + set_axis_fontsize(ax2, 12) + + ax3.fill_between(time, avg_rate+rate_std, avg_rate-rate_std, color="dodgerblue", + alpha=0.5, label="standard deviation") + ax3.plot(time, avg_rate, label="average rate") + ax3.set_xlabel("times [s]", fontsize=12) + ax3.set_ylabel("firing rate [Hz]", fontsize=12) + ax3.legend(fontsize=12) + ax3.set_xlim([0, 10]) + ax3.set_ylim([0, 450]) + set_axis_fontsize(ax3, 12) + + fig.set_size_inches(15, 10) + fig.subplots_adjust(left=0.1, bottom=0.125, top=0.95, right=0.95) + fig.set_facecolor("white") + fig.savefig("../lectures/images/concolved_rate.png") + plt.close() + + +if __name__ == "__main__": + spike_times = spio.loadmat('lifoustim.mat')["spikes"] + # plot_isi_rate(spike_times) + # plot_bin_rate(spike_times, 0.05) + plot_conv_rate(spike_times, 0.05)