From 55db76ae0c2f9e53eae3427561cfa42836578016 Mon Sep 17 00:00:00 2001
From: Jan Grewe <jan.grewe@g-node.org>
Date: Mon, 2 Nov 2015 18:57:56 +0100
Subject: [PATCH] project fixes, methods figure

---
 .../project_adaptation_fit/adaptation_fit.tex |   3 +-
 projects/project_onset_fi/onset_fi.tex        |   3 +-
 .../stimulus_reconstruction.tex               |   6 +-
 .../vector_strength.tex                       |   5 +-
 spike_trains/code/isi_method.py               | 150 ++++++++++++++++++
 5 files changed, 161 insertions(+), 6 deletions(-)
 create mode 100644 spike_trains/code/isi_method.py

diff --git a/projects/project_adaptation_fit/adaptation_fit.tex b/projects/project_adaptation_fit/adaptation_fit.tex
index c9b6af0..80bca6f 100644
--- a/projects/project_adaptation_fit/adaptation_fit.tex
+++ b/projects/project_adaptation_fit/adaptation_fit.tex
@@ -42,7 +42,8 @@ electroreceptors of the weakly electric fish \textit{Apteronotus
   \question In the accompanying datasets you find the
   \textit{spike\_times} of an P-unit electrorecptor to a stimulus of a
   certain intensity, i.e. the \textit{contrast} which is also stored
-  in the file.
+  in the file. The data is sampled with 20\,kHz sampling frequency and
+  spike times are given in milliseconds relative to the stimulus onset.
   \begin{parts}
     \part Estimate for each stimulus intensity the PSTH and plot
     it. You will see that there are three parts.  (i) The first
diff --git a/projects/project_onset_fi/onset_fi.tex b/projects/project_onset_fi/onset_fi.tex
index b454a25..776c3a7 100644
--- a/projects/project_onset_fi/onset_fi.tex
+++ b/projects/project_onset_fi/onset_fi.tex
@@ -39,7 +39,8 @@ of the stimulus \textbf{I}ntensity.
   \question In the accompanying datasets you find the
   \textit{spike\_times} of an P-unit electroreceptor of the weakly
   electric fish \textit{Apteronotus leptorhynchus} to a stimulus of a
-  certain intensity, i.e. the \textit{contrast}. 
+  certain intensity, i.e. the \textit{contrast}. The spike times are
+  given in milliseconds relative to the stimulus onset.
   \begin{parts}
     \part For each stimulus intensity estimate the average response
     (PSTH) and plot it. You will see that there are three parts.  (i)
diff --git a/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex b/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex
index 1ac3390..3a44cf2 100644
--- a/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex
+++ b/projects/project_stimulus_reconstruction/stimulus_reconstruction.tex
@@ -45,8 +45,10 @@ reconstruct the stimulus a neuron has been stimulated with.
   \question In the accompanying files you find the spike responses of
   P-units and pyramidal neurons of the weakly electric fish
   \textit{Apteronotus leptorhynchus}. The respective stimuli are
-  stored in separate files.  Start with the P-unit and, in the end,
-  apply the same functions to the pyramidal data.
+  stored in separate files. The data is sampled with 20\,kHz temporal
+  resolution and spike times are given in seconds. Start with the
+  P-unit and, in the end, apply the same functions to the pyramidal
+  data.
   \begin{parts}
     \part Estimate the STA and plot it.
     \part Implement a function that estimates the reconstruction 
diff --git a/projects/project_vector_strength/vector_strength.tex b/projects/project_vector_strength/vector_strength.tex
index a98bb46..a5e490d 100755
--- a/projects/project_vector_strength/vector_strength.tex
+++ b/projects/project_vector_strength/vector_strength.tex
@@ -45,11 +45,12 @@ varies between $0$ and $1$ for no phase locking to perfect phase
 locking, respectively.
 
 \begin{questions}
-  \question In the accompanying datasets you find recrordings of the
+  \question In the accompanying datasets you find recordings of the
   ``baseline'' activity of P-unit electroreceptors of different weakly
   electric fish of the species \textit{Apteronotus leptorhynchus}.
   The files further contain respective recordings of the \textit{eod},
-  i.e. the fish's field.
+  i.e. the fish's field The data is sampled with 20\,kHz and the spike
+  times are given in seconds.
   \begin{parts}
     \part Plot an average of the single EOD cylces of each fish 
     together with an respective PSTH.
diff --git a/spike_trains/code/isi_method.py b/spike_trains/code/isi_method.py
new file mode 100644
index 0000000..776c68e
--- /dev/null
+++ b/spike_trains/code/isi_method.py
@@ -0,0 +1,150 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from IPython import embed
+
+def set_rc():
+    plt.rcParams['xtick.labelsize'] = 8
+    plt.rcParams['ytick.labelsize'] = 8
+    plt.rcParams['xtick.major.size'] = 5
+    plt.rcParams['xtick.minor.size'] = 5
+    plt.rcParams['xtick.major.width'] = 2
+    plt.rcParams['xtick.minor.width'] = 2
+    plt.rcParams['ytick.major.size'] = 5
+    plt.rcParams['ytick.minor.size'] = 5
+    plt.rcParams['ytick.major.width'] = 2
+    plt.rcParams['ytick.minor.width'] = 2
+    plt.rcParams['xtick.direction'] = "out"
+    plt.rcParams['ytick.direction'] = "out"
+
+
+def create_spikes(isi=0.08, duration=0.5):
+    times = np.arange(0., duration, isi)
+    times += np.random.randn(len(times)) * (isi / 2.5)
+    times = np.delete(times, np.nonzero(times < 0))
+    times = np.delete(times, np.nonzero(times > duration))
+    times = np.sort(times)
+    return times
+
+
+def gaussian(sigma, dt):
+    x = np.arange(-4*sigma, 4*sigma, dt)
+    y = np.exp(-0.5 * (x / sigma)**2)/np.sqrt(2*np.pi)/sigma; 
+    return x, y
+
+    
+def setup_axis(spikes_ax, rate_ax):
+    spikes_ax.spines["right"].set_visible(False)
+    spikes_ax.spines["top"].set_visible(False)
+    spikes_ax.yaxis.set_ticks_position('left')
+    spikes_ax.xaxis.set_ticks_position('bottom')
+    spikes_ax.set_yticks([0, 1.0])
+    spikes_ax.set_ylim([0, 1.05])
+    spikes_ax.set_ylabel("spikes", fontsize=9)
+    spikes_ax.text(-0.125, 1.2, "A", transform=spikes_ax.transAxes, size=10)
+
+    rate_ax.spines["right"].set_visible(False)
+    rate_ax.spines["top"].set_visible(False)
+    rate_ax.yaxis.set_ticks_position('left')
+    rate_ax.xaxis.set_ticks_position('bottom')
+    rate_ax.set_xlabel('time[s]', fontsize=9)
+    rate_ax.set_ylabel('firing rate [Hz]', fontsize=9)
+    rate_ax.text(-0.125, 1.15, "B", transform=rate_ax.transAxes, size=10)   
+
+
+def plot_bin_method():
+    dt = 1e-5
+    duration = 0.5
+    
+    spike_times = create_spikes(0.025, duration)
+    t = np.arange(0., duration, dt)
+    
+    bins = np.arange(0, 0.5, 0.05)
+    count, _ = np.histogram(spike_times, bins)
+    
+    plt.xkcd()
+    fig = plt.figure()
+    fig.set_size_inches(5., 2.5)
+    fig.set_facecolor('white')
+    spikes = plt.subplot2grid((3,1), (0,0), rowspan=1, colspan=1)
+    rate_ax = plt.subplot2grid((3,1), (1,0), rowspan=2, colspan=1)
+    setup_axis(spikes, rate_ax)
+
+    spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1)    
+    spikes.vlines(bins, 0., 1.05, color="red", lw = 2)
+    for c in count:
+        spikes.text(bins-(bins[1]-bins[0])/2, 0.75, str(c), fontdict={"color":"red"})
+    spikes.set_xlim([0, duration])
+    
+    #rate_ax.step(bins[1:], count/0.05)
+    
+    fig.tight_layout()
+    plt.show()
+    #fig.savefig("bin_method.pdf")
+
+
+def plot_conv_method():
+    dt = 1e-5
+    duration = 0.5
+    spike_times = create_spikes(0.05, duration) 
+    kernel_time, kernel = gaussian(0.02, dt)
+    
+    t = np.arange(0., duration, dt)
+    rate = np.zeros(t.shape)
+    rate[np.asarray(np.round(spike_times/dt), dtype=int)] = 1  
+    rate = np.convolve(rate, kernel, mode='same')
+    rate = np.roll(rate, -1)
+
+    fig = plt.figure()
+    fig.set_size_inches(5., 2.5)
+    fig.set_facecolor('white')
+    spikes = plt.subplot2grid((3,1), (0,0), rowspan=1, colspan=1)
+    rate_ax = plt.subplot2grid((3,1), (1,0), rowspan=2, colspan=1)
+    setup_axis(spikes, rate_ax)
+    
+    spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.5, zorder=2)
+    for i in spike_times:
+        spikes.plot(kernel_time + i, kernel/np.max(kernel), color="orange", lw=0.75, zorder=1)
+    spikes.set_xlim([0, duration])
+
+    rate_ax.plot(t, rate, color="darkblue", lw=1, zorder=2)
+    rate_ax.fill_between(t, rate, np.zeros(len(rate)), color="red", alpha=0.5)
+    rate_ax.set_xlim([0, duration])
+    rate_ax.set_ylim([0, 50])
+    rate_ax.set_yticks(np.arange(0,75,25))
+    fig.tight_layout()
+    fig.savefig("conv_method.pdf")
+
+
+def plot_isi_method():
+    spike_times = create_spikes(0.08, 0.5)
+    
+    plt.xkcd()
+    set_rc()
+    fig = plt.figure()
+    fig.set_size_inches(5., 2.5)
+    fig.set_facecolor('white')
+    
+    spikes = plt.subplot2grid((3,1), (0,0), rowspan=1, colspan=1)
+    rate = plt.subplot2grid((3,1), (1,0), rowspan=2, colspan=1)
+    setup_axis(spikes, rate)
+    
+    spikes.vlines(spike_times, 0., 1., color="darkblue", lw=1.25)
+    spike_times = np.hstack((0, spike_times))
+    for i in range(1, len(spike_times)):
+        t_start = spike_times[i-1]
+        t = spike_times[i]
+        spikes.annotate(s='', xy=(t_start, 0.5), xytext=(t,0.5), arrowprops=dict(arrowstyle='<->'))
+    i_rate = 1./np.diff(spike_times)
+
+    rate.step(spike_times[1:], i_rate,color="darkblue", lw=1.25, where="pre")
+    rate.set_ylim([0, 75])
+    rate.set_yticks(np.arange(0,100,25))
+    
+    fig.tight_layout()
+    fig.savefig("isi_method.pdf")
+    
+    
+if __name__ == '__main__':
+    #plot_isi_method()
+    #plot_conv_method()
+    plot_bin_method()