From 54c3c4896adc2ca2ec57883470164481ecf4f4ba Mon Sep 17 00:00:00 2001
From: efish <efish@localhost.localdomain>
Date: Thu, 22 Nov 2018 15:39:26 +0100
Subject: [PATCH] plot

---
 code/plot_eodform_spikehist.py | 84 ++++++++++++++++++++++++++++++++++
 1 file changed, 84 insertions(+)
 create mode 100644 code/plot_eodform_spikehist.py

diff --git a/code/plot_eodform_spikehist.py b/code/plot_eodform_spikehist.py
new file mode 100644
index 0000000..a715adf
--- /dev/null
+++ b/code/plot_eodform_spikehist.py
@@ -0,0 +1,84 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from read_baseline_data import *
+from NixFrame import *
+from utility import *
+from IPython import embed
+
+# plot and data values
+inch_factor = 2.54
+data_dir = '../data'
+dataset = '2018-11-09-ad-invivo-1'
+
+# read eod and time of baseline
+time, eod = read_baseline_eod(os.path.join(data_dir, dataset))
+
+
+
+# read spikes during baseline activity
+spikes = read_baseline_spikes(os.path.join(data_dir, dataset))
+# calculate interpike intervals and plot them
+interspikeintervals = np.diff(spikes)
+
+fig, ax = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
+plt.hist(interspikeintervals, bins=np.arange(0, np.max(interspikeintervals), 0.0001))
+plt.show()
+
+# calculate coefficient of variation
+mu = np.mean(interspikeintervals)
+sigma = np.std(interspikeintervals)
+cv = sigma/mu
+print(cv)
+
+# calculate eod times and indices by zero crossings
+threshold = 0
+shift_eod = np.roll(eod, 1)
+eod_times = time[(eod >= threshold) & (shift_eod < threshold)]
+sampling_rate = 40000.0
+eod_idx = eod_times*sampling_rate
+
+# align eods and spikes to eods
+max_cut = int(np.max(np.diff(eod_idx)))
+eod_cuts = np.zeros([len(eod_idx)-1, max_cut])
+spike_times = []
+eod_durations = []
+
+for i, idx in enumerate(eod_idx[:-1]):
+	eod_cut = eod[int(idx):int(eod_idx[i+1])]
+	eod_cuts[i, :len(eod_cut)] = eod_cut
+	eod_cuts[i, len(eod_cut):] = np.nan
+	time_cut = time[int(idx):int(eod_idx[i+1])]
+	spike_cut = spikes[(spikes > time_cut[0]) & (spikes < time_cut[-1])]
+	spike_time = spike_cut - time_cut[0]
+	if len(spike_time) > 0:
+		spike_times.append(spike_time[:][0]*1000)
+		eod_durations.append(len(eod_cut)/sampling_rate*1000)
+
+# calculate vector strength
+vs = vector_strength(spike_times, eod_durations)
+
+# determine means and stds of eod for plot
+# determine time axis
+mu_eod = np.nanmean(eod_cuts, axis=0)
+std_eod = np.nanstd(eod_cuts, axis=0)*3
+time_axis = np.arange(max_cut)/sampling_rate*1000
+
+# plot eod form and spike histogram
+fig, ax1 = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
+ax1.hist(spike_times, color='crimson')
+ax1.set_xlabel('time [ms]', fontsize=12)
+ax1.set_ylabel('number', fontsize=12)
+ax1.tick_params(axis='y', labelcolor='crimson')
+plt.yticks(fontsize=8)
+ax1.spines['top'].set_visible(False)
+
+ax2 = ax1.twinx()
+ax2.fill_between(time_axis, mu_eod+std_eod, mu_eod-std_eod, color='dodgerblue', alpha=0.5)
+ax2.plot(time_axis, mu_eod, color='black', lw=2)
+ax2.set_ylabel('voltage [mV]', fontsize=12)
+ax2.tick_params(axis='y', labelcolor='dodgerblue')
+
+plt.xticks(fontsize=8)
+plt.yticks(fontsize=8)
+fig.tight_layout()
+plt.show()
\ No newline at end of file