From cf4d63966af1031f3ec096f2173f557430a89268 Mon Sep 17 00:00:00 2001
From: Ramona <efish@verliernix.com>
Date: Thu, 15 Nov 2018 15:48:34 +0100
Subject: [PATCH] vector strength

---
 code/analysis_rs.py | 11 ++++++++++-
 code/utility.py     | 10 ++++++----
 2 files changed, 16 insertions(+), 5 deletions(-)

diff --git a/code/analysis_rs.py b/code/analysis_rs.py
index 106fb28..7a50853 100644
--- a/code/analysis_rs.py
+++ b/code/analysis_rs.py
@@ -3,6 +3,7 @@ import matplotlib.pyplot as plt
 from read_baseline_data import *
 from IPython import embed
 from NixFrame import *
+from utility import *
 
 inch_factor = 2.54
 data_dir = '../data'
@@ -53,6 +54,7 @@ max_cut = int(np.max(np.diff(eod_idx)))
 eod_cuts = np.zeros([len(eod_idx)-1, max_cut])
 # eods 15 + 16 are to short
 relative_times = []
+eod_durations = []
 
 for i, idx in enumerate(eod_idx[:-1]):
 	eod_cut = eod[int(idx):int(eod_idx[i+1])]
@@ -63,14 +65,20 @@ for i, idx in enumerate(eod_idx[:-1]):
 	relative_time  = spike_cut - time_cut[0]
 	if len(relative_time) > 0:
 		relative_times.append(relative_time[:][0]*1000)
+		eod_durations.append(len(eod_cut))
 
 
 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
+vs = vector_strength(relative_times, eod_durations)
+embed()
+exit()
+
+#time_axis = np.arange(max_cut)/sampling_rate*1000
 
 #fig = plt.figure(figsize=(12/inch_factor, 8/inch_factor))
+'''
 fig, ax1 = plt.subplots(figsize=(12/inch_factor, 8/inch_factor))
 ax1.hist(relative_times, color='crimson')
 ax1.set_xlabel('time [ms]', fontsize=12)
@@ -90,4 +98,5 @@ plt.xticks(fontsize = 8)
 plt.yticks(fontsize = 8)
 fig.tight_layout()
 plt.show()
+'''
 
diff --git a/code/utility.py b/code/utility.py
index 0c1ce34..a83b5c2 100644
--- a/code/utility.py
+++ b/code/utility.py
@@ -11,8 +11,10 @@ def zero_crossing(eod,time):
 
 
 def vector_strength(spike_times, eod_durations):
-	alphas = spike_times/ eod_durations
-	cs = (1/len(spike_times))*np.sum(np.cos(alphas))^2
-	sn = (1/len(spike_times))*np.sum(np.sin(alphas))^2
-	vs = np.sqrt(cs+sn)
+	n = len(spike_times)
+	phase_times = np.zeros(n)
+	for i, idx in enumerate(spike_times):
+		phase_times[i]= spike_times[i]/eod_durations[i]
+
+	vs = np.sqrt((1/n*sum(np.cos(phase_times)))**2 + (1/n*sum(np.sin(phase_times)))**2)
 	return vs