From ea8add517c04d393b46ab1f736dfb79b9cb84166 Mon Sep 17 00:00:00 2001 From: Jan Benda Date: Sun, 10 Jan 2021 20:37:40 +0100 Subject: [PATCH] [pointprocesses] python version of sketch --- plotstyle.py | 2 + pointprocesses/lecture/firingrates.py | 26 +++--- pointprocesses/lecture/isihexamples.py | 88 ++++++++++-------- pointprocesses/lecture/pointprocesses.tex | 9 +- pointprocesses/lecture/rasterexamples.py | 85 ++++++++++-------- pointprocesses/lecture/returnmapexamples.py | 94 -------------------- pointprocesses/lecture/serialcorrexamples.py | 94 +++++++++++++------- 7 files changed, 179 insertions(+), 219 deletions(-) delete mode 100644 pointprocesses/lecture/returnmapexamples.py diff --git a/plotstyle.py b/plotstyle.py index 5cbf30e..add3e40 100644 --- a/plotstyle.py +++ b/plotstyle.py @@ -37,6 +37,7 @@ edgewidth = 0.0 if xkcd_style else 1.0 mainline = {'linestyle': '-', 'linewidth': lwthick} minorline = {'linestyle': '-', 'linewidth': lwthin} largemarker = {'marker': 'o', 'markersize': 9, 'markeredgecolor': colors['white'], 'markeredgewidth': edgewidth} +largeopenmarker = {'marker': 'o', 'markersize': 7, 'markerfacecolor': colors['white'], 'markeredgewidth': 2} smallmarker = {'marker': 'o', 'markersize': 6, 'markeredgecolor': colors['white'], 'markeredgewidth': edgewidth} largelinepoints = {'linestyle': '-', 'linewidth': lwthick, 'marker': 'o', 'markersize': 10, 'markeredgecolor': colors['white'], 'markeredgewidth': 1} smalllinepoints = {'linestyle': '-', 'linewidth': 1.4, 'marker': 'o', 'markersize': 7, 'markeredgecolor': colors['white'], 'markeredgewidth': 1} @@ -88,6 +89,7 @@ fsAa = {'facecolor': colors['blue'], 'edgecolor': 'none', 'alpha': fillalpha} lsB = dict({'color': colors['red']}, **mainline) lsBm = dict({'color': colors['red']}, **minorline) psB = dict({'color': colors['red'], 'linestyle': 'none'}, **largemarker) +psBo = dict({'markeredgecolor': colors['red'], 'linestyle': 'none'}, **largeopenmarker) psBm = dict({'color': colors['red'], 'linestyle': 'none'}, **smallmarker) lpsB = dict({'color': colors['red']}, **largelinepoints) lpsBm = dict({'color': colors['red']}, **smalllinepoints) diff --git a/pointprocesses/lecture/firingrates.py b/pointprocesses/lecture/firingrates.py index 087e7ee..1195e3a 100644 --- a/pointprocesses/lecture/firingrates.py +++ b/pointprocesses/lecture/firingrates.py @@ -6,7 +6,7 @@ import matplotlib.pyplot as plt from plotstyle import * -def get_instantaneous_rate(times, max_t=30., dt=1e-4): +def 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))) @@ -19,12 +19,12 @@ def get_instantaneous_rate(times, max_t=30., dt=1e-4): def plot_isi_rate(spike_times, max_t=30, dt=1e-4): times = np.squeeze(spike_times[0][0])[:50000] - time, rate = get_instantaneous_rate(times, max_t=50000*dt) + time, rate = instantaneous_rate(times, max_t=50000*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])[:50000], - max_t=50000*dt) + _, rates[:, i] = instantaneous_rate(np.squeeze(spike_times[i][0])[:50000], + max_t=50000*dt) avg_rate = np.mean(rates, axis=1) rate_std = np.std(rates, axis=1) @@ -48,7 +48,7 @@ def plot_isi_rate(spike_times, max_t=30, dt=1e-4): plt.close() -def get_binned_rate(times, bin_width=0.05, max_t=30., dt=1e-4): +def 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 = np.asarray(bins / dt, np.int) @@ -62,10 +62,10 @@ def get_binned_rate(times, bin_width=0.05, max_t=30., dt=1e-4): 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) + time, rate = 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])) + _, rates[:, i] = binned_rate(np.squeeze(spike_times[i][0])) avg_rate = np.mean(rates, axis=1) rate_std = np.std(rates, axis=1) @@ -93,7 +93,7 @@ def plot_bin_rate(spike_times, bin_width, max_t=30, dt=1e-4): plt.close() -def get_convolved_rate(times, sigma, max_t=30., dt=1.e-4): +def convolved_rate(times, sigma, max_t=30., dt=1.e-4): time = np.arange(0., max_t, dt) kernel = spst.norm.pdf(np.arange(-8*sigma, 8*sigma, dt),loc=0,scale=sigma) indices = np.asarray(times/dt, dtype=int) @@ -105,11 +105,11 @@ def get_convolved_rate(times, sigma, max_t=30., dt=1.e-4): 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) + time, rate = 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) + _, rates[:, i] = convolved_rate(np.squeeze(spike_times[i][0]), sigma) avg_rate = np.mean(rates, axis=1) rate_std = np.std(rates, axis=1) @@ -139,9 +139,9 @@ def plot_conv_rate(spike_times, sigma=0.05, max_t=30, dt=1e-4): def plot_comparison(spike_times, bin_width, sigma, max_t=30., dt=1e-4): times = np.squeeze(spike_times[0][0]) - time, conv_rate = get_convolved_rate(times, bin_width/np.sqrt(12.)) - time, inst_rate = get_instantaneous_rate(times) - time, binn_rate = get_binned_rate(times, bin_width) + time, conv_rate = convolved_rate(times, bin_width/np.sqrt(12.)) + time, inst_rate = instantaneous_rate(times) + time, binn_rate = binned_rate(times, bin_width) fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=cm_size(figure_width, 1.8*figure_height)) fig.subplots_adjust(**adjust_fs(fig, left=6.0, right=1.5, bottom=3.0, top=1.0)) diff --git a/pointprocesses/lecture/isihexamples.py b/pointprocesses/lecture/isihexamples.py index 674bbd7..9d0f7c0 100644 --- a/pointprocesses/lecture/isihexamples.py +++ b/pointprocesses/lecture/isihexamples.py @@ -2,6 +2,15 @@ import numpy as np import matplotlib.pyplot as plt from plotstyle import * + +rate = 20.0 +trials = 10 +duration = 100.0 +dt = 0.001 +drate = 50.0 +tau = 0.1; + + def hompoisson(rate, trials, duration) : spikes = [] for k in range(trials) : @@ -13,6 +22,7 @@ def hompoisson(rate, trials, duration) : spikes.append( times ) return spikes + def inhompoisson(rate, trials, dt) : spikes = [] p = rate*dt @@ -40,12 +50,27 @@ def pifspikes(input, trials, dt, D=0.1) : spikes.append( times ) return spikes + +def oupifspikes(rate, trials, duration, dt, D, drate, tau): + # OU noise: + rng = np.random.RandomState(54637281) + time = np.arange(0.0, duration, dt) + x = np.zeros(time.shape)+rate + n = rng.randn(len(time))*drate*tau/np.sqrt(dt) + rate + for k in range(1,len(x)) : + x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau + x[x<0.0] = 0.0 + spikes = pifspikes(x, trials, dt, D) + return spikes + + def isis( spikes ) : isi = [] for k in range(len(spikes)) : isi.extend(np.diff(spikes[k])) return isi + def plotisih( ax, isis, binwidth=None ) : if binwidth == None : nperbin = 200.0 # average number of isis per bin @@ -61,42 +86,29 @@ def plotisih( ax, isis, binwidth=None ) : ax.set_ylabel('p(ISI)', '1/s') ax.bar( 1000.0*b[:-1], h, bar_fac*1000.0*np.diff(b), **fsA) -# parameter: -rate = 20.0 -drate = 50.0 -trials = 10 -duration = 100.0 -dt = 0.001 -tau = 0.1; + +def plot_hom_isih(ax): + homspikes = hompoisson(rate, trials, duration) + ax.set_xlim(0.0, 150.0) + ax.set_ylim(0.0, 31.0) + ax.set_xticks(np.arange(0.0, 151.0, 50.0)) + ax.set_yticks(np.arange(0.0, 31.0, 10.0)) + plotisih(ax, isis(homspikes), 0.005) + + +def plot_inhom_isih(ax): + inhspikes = oupifspikes(rate, trials, duration, dt, 0.3, drate, tau) + ax.set_xlim(0.0, 150.0) + ax.set_ylim(0.0, 31.0) + ax.set_xticks(np.arange(0.0, 151.0, 50.0)) + ax.set_yticks(np.arange(0.0, 31.0, 10.0)) + plotisih(ax, isis(inhspikes), 0.005) -# homogeneous spike trains: -homspikes = hompoisson(rate, trials, duration) - -# OU noise: -rng = np.random.RandomState(54637281) -time = np.arange(0.0, duration, dt) -x = np.zeros(time.shape)+rate -n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate -for k in range(1,len(x)) : - x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau -x[x<0.0] = 0.0 - -# pif spike trains: -inhspikes = pifspikes(x, trials, dt, D=0.3) - -fig, (ax1, ax2) = plt.subplots(1, 2) -fig.subplots_adjust(**adjust_fs(fig, top=0.5, right=1.5)) -ax1.set_xlim(0.0, 150.0) -ax1.set_ylim(0.0, 31.0) -ax1.set_xticks(np.arange(0.0, 151.0, 50.0)) -ax1.set_yticks(np.arange(0.0, 31.0, 10.0)) -plotisih(ax1, isis(homspikes), 0.005) - -ax2.set_xlim(0.0, 150.0) -ax2.set_ylim(0.0, 31.0) -ax2.set_xticks(np.arange(0.0, 151.0, 50.0)) -ax2.set_yticks(np.arange(0.0, 31.0, 10.0)) -plotisih(ax2, isis(inhspikes), 0.005) - -plt.savefig('isihexamples.pdf') -plt.close() + +if __name__ == "__main__": + fig, (ax1, ax2) = plt.subplots(1, 2) + fig.subplots_adjust(**adjust_fs(fig, top=0.5, right=1.5)) + plot_hom_isih(ax1) + plot_inhom_isih(ax2) + plt.savefig('isihexamples.pdf') + plt.close() diff --git a/pointprocesses/lecture/pointprocesses.tex b/pointprocesses/lecture/pointprocesses.tex index 81bbb9a..da03d13 100644 --- a/pointprocesses/lecture/pointprocesses.tex +++ b/pointprocesses/lecture/pointprocesses.tex @@ -60,8 +60,8 @@ process]{Punktprozess}{point processes}. \texpicture{pointprocessscetch} \titlecaption{\label{pointprocessscetchfig} Statistics of point processes.}{A point process is a sequence of instances in time - $t_i$ that can be characterized through the inter-event-intervals - $T_i=t_{i+1}-t_i$ and the number of events $n_i$. } + $t_i$ that can be also characterized by inter-event intervals + $T_i=t_{i+1}-t_i$ and event counts $n_i$.} \end{figure} \noindent @@ -150,9 +150,8 @@ the interval $T_i$. The parameter $k$ is called the \enterm{lag} maps are distinctly different \figref{returnmapfig}. \begin{figure}[tp] - \includegraphics[width=1\textwidth]{returnmapexamples} \includegraphics[width=1\textwidth]{serialcorrexamples} - \titlecaption{\label{returnmapfig}Interspike interval + \titlecaption{\label{returnmapfig}Interspike-interval correlations}{of the spike trains shown in \figref{rasterexamplesfig}. Upper panels show the return maps and lower panels the serial correlations of successive intervals @@ -191,7 +190,7 @@ with itself and is always 1. % \begin{figure}[t] % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz10ms}\hfill % \includegraphics[width=0.48\textwidth]{poissoncounthist100hz100ms} -% \titlecaption{\label{countstatsfig}Count Statistik.}{} +% \titlecaption{\label{countstatsfig}Count statistic.}{} % \end{figure} Counting the number of events $n_i$ (counts) in time windows $i$ of duration $W$ yields positive integer random numbers that are commonly quantified diff --git a/pointprocesses/lecture/rasterexamples.py b/pointprocesses/lecture/rasterexamples.py index f616c18..07b8628 100644 --- a/pointprocesses/lecture/rasterexamples.py +++ b/pointprocesses/lecture/rasterexamples.py @@ -2,6 +2,15 @@ import numpy as np import matplotlib.pyplot as plt from plotstyle import * + +rate = 20.0 +trials = 10 +duration = 2.0 +dt = 0.001 +drate = 50.0 +tau = 0.1; + + def hompoisson(rate, trials, duration) : spikes = [] for k in range(trials) : @@ -13,6 +22,7 @@ def hompoisson(rate, trials, duration) : spikes.append(times) return spikes + def inhompoisson(rate, trials, dt) : spikes = [] p = rate*dt @@ -40,47 +50,44 @@ def pifspikes(input, trials, dt, D=0.1) : spikes.append(times) return spikes -# parameter: -rate = 20.0 -drate = 50.0 -trials = 10 -duration = 2.0 -dt = 0.001 -tau = 0.1; - -# homogeneous spike trains: -homspikes = hompoisson(rate, trials, duration) -# OU noise: -rng = np.random.RandomState(54637281) -time = np.arange(0.0, duration, dt) -x = np.zeros(time.shape)+rate -n = rng.randn(len(time))*drate*tau/np.sqrt(dt)+rate -for k in range(1,len(x)) : - x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau -x[x<0.0] = 0.0 - -# inhomogeneous spike trains: -#inhspikes = inhompoisson(x, trials, dt) -# pif spike trains: -inhspikes = pifspikes(x, trials, dt, D=0.3) +def oupifspikes(rate, trials, duration, dt, D, drate, tau): + # OU noise: + rng = np.random.RandomState(54637281) + time = np.arange(0.0, duration, dt) + x = np.zeros(time.shape)+rate + n = rng.randn(len(time))*drate*tau/np.sqrt(dt) + rate + for k in range(1,len(x)) : + x[k] = x[k-1] + (n[k]-x[k-1])*dt/tau + x[x<0.0] = 0.0 + spikes = pifspikes(x, trials, dt, D) + return spikes -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=cm_size(figure_width, 0.5*figure_width)) -fig.subplots_adjust(**adjust_fs(fig, left=4.0, right=1.0, top=1.2)) + +def plot_homogeneous_spikes(ax): + homspikes = hompoisson(rate, trials, duration) + ax.set_title('stationary') + ax.set_xlim(0.0, duration) + ax.set_ylim(-0.5, trials-0.5) + ax.set_xlabel('Time [s]') + ax.set_ylabel('Trial') + ax.eventplot(homspikes, colors=[lsA['color']], linelength=0.8, lw=1) -ax1.set_title('stationary') -ax1.set_xlim(0.0, duration) -ax1.set_ylim(-0.5, trials-0.5) -ax1.set_xlabel('Time [s]') -ax1.set_ylabel('Trial') -ax1.eventplot(homspikes, colors=[lsA['color']], linelength=0.8, lw=1) + +def plot_inhomogeneous_spikes(ax): + inhspikes = oupifspikes(rate, trials, duration, dt, 0.3, drate, tau) + ax.set_title('non-stationary') + ax.set_xlim(0.0, duration) + ax.set_ylim(-0.5, trials-0.5) + ax.set_xlabel('Time [s]') + ax.set_ylabel('Trial') + ax.eventplot(inhspikes, colors=[lsA['color']], linelength=0.8, lw=1) -ax2.set_title('non-stationary') -ax2.set_xlim(0.0, duration) -ax2.set_ylim(-0.5, trials-0.5) -ax2.set_xlabel('Time [s]') -ax2.set_ylabel('Trial') -ax2.eventplot(inhspikes, colors=[lsA['color']], linelength=0.8, lw=1) -plt.savefig('rasterexamples.pdf') -plt.close() +if __name__ == "__main__": + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=cm_size(figure_width, 0.5*figure_width)) + fig.subplots_adjust(**adjust_fs(fig, left=4.0, right=1.0, top=1.2)) + plot_homogeneous_spikes(ax1) + plot_inhomogeneous_spikes(ax2) + plt.savefig('rasterexamples.pdf') + plt.close() diff --git a/pointprocesses/lecture/returnmapexamples.py b/pointprocesses/lecture/returnmapexamples.py deleted file mode 100644 index 99bf88d..0000000 --- a/pointprocesses/lecture/returnmapexamples.py +++ /dev/null @@ -1,94 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from plotstyle import * - -def hompoisson(rate, trials, duration) : - spikes = [] - for k in range(trials) : - times = [] - t = 0.0 - while t < duration : - t += np.random.exponential(1/rate) - times.append( t ) - spikes.append( times ) - return spikes - -def inhompoisson(rate, trials, dt) : - spikes = [] - p = rate*dt - for k in range(trials) : - x = np.random.rand(len(rate)) - times = dt*np.nonzero(x= vthresh : - v = vreset - times.append(k*dt) - spikes.append( times ) - return spikes - - -def isis( spikes ) : - isi = [] - for k in range(len(spikes)) : - isi.extend(np.diff(spikes[k])) - return np.array( isi ) - - -def plotreturnmap(ax, isis, lag=1, max=1.0) : - ax.set_xlabel(r'ISI$_i$', 'ms') - ax.set_ylabel(r'ISI$_{i+1}$', 'ms') - ax.set_xlim(0.0, 1000.0*max) - ax.set_ylim(0.0, 1000.0*max) - isiss = isis[isis