python script to create firing rates in different ways
This commit is contained in:
parent
cc2515f22b
commit
c842726950
198
programming/code/firing_rates.py
Normal file
198
programming/code/firing_rates.py
Normal file
@ -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)
|
||||
Reference in New Issue
Block a user