actions/stimulus-spike-lfp-response/

actions/stimulus-spike-lfp-response/data/10-calculate-stimulus-spike-lfp-response.html

@@ -13115,7 +13115,7 @@ div#notebook {
 
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[1]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="o">%</span><span class="k">load_ext</span> autoreload
@@ -13129,7 +13129,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[2]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
@@ -13140,6 +13140,8 @@ div#notebook {
 <span class="kn">import</span> <span class="nn">expipe</span>
 <span class="kn">import</span> <span class="nn">os</span>
 <span class="kn">import</span> <span class="nn">pathlib</span>
+<span class="kn">import</span> <span class="nn">scipy</span>
+<span class="kn">import</span> <span class="nn">scipy.signal</span>
 <span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
 <span class="kn">import</span> <span class="nn">exdir</span>
 <span class="kn">import</span> <span class="nn">pandas</span> <span class="k">as</span> <span class="nn">pd</span>
@@ -13162,28 +13164,10 @@ div#notebook {
 </div>
 </div>
 
-<div class="output_wrapper">
-<div class="output">
-
-
-<div class="output_area">
-
-    <div class="prompt"></div>
-
-
-<div class="output_subarea output_stream output_stderr output_text">
-<pre>10:25:26 [I] klustakwik KlustaKwik2 version 0.2.6
-</pre>
-</div>
-</div>
-
-</div>
-</div>
-
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[3]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="n">data_loader</span> <span class="o">=</span> <span class="n">dp</span><span class="o">.</span><span class="n">Data</span><span class="p">()</span>
@@ -13198,7 +13182,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[4]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="n">output</span> <span class="o">=</span> <span class="n">pathlib</span><span class="o">.</span><span class="n">Path</span><span class="p">(</span><span class="s1">&#39;output/stimulus-spike-lfp-response&#39;</span><span class="p">)</span>
@@ -13212,7 +13196,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[5]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="n">identify_neurons</span> <span class="o">=</span> <span class="n">actions</span><span class="p">[</span><span class="s1">&#39;identify-neurons&#39;</span><span class="p">]</span>
@@ -13227,7 +13211,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[6]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">get_lim</span><span class="p">(</span><span class="n">action_id</span><span class="p">):</span>
@@ -13258,7 +13242,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[7]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">signaltonoise</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">ddof</span><span class="o">=</span><span class="mi">0</span><span class="p">):</span>
@@ -13266,16 +13250,6 @@ div#notebook {
     <span class="n">m</span> <span class="o">=</span> <span class="n">a</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">axis</span><span class="p">)</span>
     <span class="n">sd</span> <span class="o">=</span> <span class="n">a</span><span class="o">.</span><span class="n">std</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="n">axis</span><span class="p">,</span> <span class="n">ddof</span><span class="o">=</span><span class="n">ddof</span><span class="p">)</span>
     <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="n">sd</span> <span class="o">==</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">m</span> <span class="o">/</span> <span class="n">sd</span><span class="p">)</span>
-
-
-<span class="k">def</span> <span class="nf">compute_clean_lfp</span><span class="p">(</span><span class="n">anas</span><span class="p">,</span> <span class="n">width</span><span class="o">=</span><span class="mi">500</span><span class="p">,</span> <span class="n">threshold</span><span class="o">=</span><span class="mi">2</span><span class="p">):</span>
-    <span class="n">anas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">anas</span><span class="p">)</span>
-
-    <span class="k">for</span> <span class="n">ch</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">anas</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">1</span><span class="p">]):</span>
-        <span class="n">idxs</span><span class="p">,</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">anas</span><span class="p">[:,</span> <span class="n">ch</span><span class="p">])</span> <span class="o">&gt;</span> <span class="n">threshold</span><span class="p">)</span>
-        <span class="k">for</span> <span class="n">idx</span> <span class="ow">in</span> <span class="n">idxs</span><span class="p">:</span>
-            <span class="n">anas</span><span class="p">[</span><span class="n">idx</span><span class="o">-</span><span class="n">width</span><span class="p">:</span><span class="n">idx</span><span class="o">+</span><span class="n">width</span><span class="p">,</span> <span class="n">ch</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span> <span class="c1"># TODO AR model prediction</span>
-    <span class="k">return</span> <span class="n">anas</span>
 </pre></div>
 
     </div>
@@ -13285,7 +13259,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[8]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">compute_energy</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">f</span><span class="p">,</span> <span class="n">f1</span><span class="p">,</span> <span class="n">f2</span><span class="p">):</span>
@@ -13303,7 +13277,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[9]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">find_theta_peak</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">f</span><span class="p">,</span> <span class="n">f1</span><span class="p">,</span> <span class="n">f2</span><span class="p">):</span>
@@ -13324,7 +13298,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[10]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">compute_half_width</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">f</span><span class="p">,</span> <span class="n">m_p</span><span class="p">,</span> <span class="n">m_f</span><span class="p">):</span>
@@ -13354,7 +13328,36 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[11]:</div>
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># p = np.load(&#39;debug_p.npy&#39;)</span>
+<span class="c1"># f = np.load(&#39;debug_f.npy&#39;)</span>
+<span class="c1"># compute_half_width(p, f, 0.01038941, 30.30187709636872)</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[12]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># plt.plot(f, p)</span>
+<span class="c1"># plt.xlim(29.9,30.6)</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[13]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">compute_stim_peak</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">f</span><span class="p">,</span> <span class="n">s_f</span><span class="p">):</span>
@@ -13367,81 +13370,16 @@ div#notebook {
 </div>
 </div>
 
-</div>
-<div class="cell border-box-sizing code_cell rendered">
-<div class="input">
-<div class="prompt input_prompt">In&nbsp;[12]:</div>
-<div class="inner_cell">
-    <div class="input_area">
-<div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">compute_spike_lfp</span><span class="p">(</span><span class="n">anas</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">t_start</span><span class="p">,</span> <span class="n">t_stop</span><span class="p">,</span> <span class="n">NFFT</span><span class="p">):</span>
-
-    <span class="n">t_start</span> <span class="o">=</span> <span class="n">t_start</span> <span class="o">*</span> <span class="n">pq</span><span class="o">.</span><span class="n">s</span> <span class="k">if</span> <span class="n">t_start</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span> <span class="k">else</span> <span class="mi">0</span> <span class="o">*</span> <span class="n">pq</span><span class="o">.</span><span class="n">s</span>
-    <span class="n">sampling_rate</span> <span class="o">=</span> <span class="n">anas</span><span class="o">.</span><span class="n">sampling_rate</span>
-    <span class="n">units</span> <span class="o">=</span> <span class="n">anas</span><span class="o">.</span><span class="n">units</span>
-    <span class="k">if</span> <span class="n">t_start</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span> <span class="ow">and</span> <span class="n">t_stop</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
-        <span class="n">t_stop</span> <span class="o">=</span> <span class="n">t_stop</span> <span class="o">*</span> <span class="n">pq</span><span class="o">.</span><span class="n">s</span>
-        <span class="n">mask</span> <span class="o">=</span> <span class="p">(</span><span class="n">anas</span><span class="o">.</span><span class="n">times</span> <span class="o">&gt;</span> <span class="n">t_start</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">anas</span><span class="o">.</span><span class="n">times</span> <span class="o">&lt;</span> <span class="n">t_stop</span><span class="p">)</span>
-        <span class="n">anas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">anas</span><span class="p">)[</span><span class="n">mask</span><span class="p">,:]</span>
-
-    <span class="n">cleaned_anas</span> <span class="o">=</span> <span class="n">zscore</span><span class="p">(</span><span class="n">compute_clean_lfp</span><span class="p">(</span><span class="n">anas</span><span class="p">))</span>
-    
-    <span class="n">cleaned_anas</span> <span class="o">=</span> <span class="n">neo</span><span class="o">.</span><span class="n">AnalogSignal</span><span class="p">(</span>
-        <span class="n">signal</span><span class="o">=</span><span class="n">cleaned_anas</span> <span class="o">*</span> <span class="n">units</span><span class="p">,</span> <span class="n">sampling_rate</span><span class="o">=</span><span class="n">sampling_rate</span><span class="p">,</span> <span class="n">t_start</span><span class="o">=</span><span class="n">t_start</span>
-    <span class="p">)</span>
-
-    <span class="n">sptr</span> <span class="o">=</span> <span class="n">neo</span><span class="o">.</span><span class="n">SpikeTrain</span><span class="p">(</span>
-        <span class="n">sptr</span><span class="o">.</span><span class="n">times</span><span class="p">[(</span><span class="n">sptr</span><span class="o">.</span><span class="n">times</span> <span class="o">&gt;</span> <span class="n">t_start</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">sptr</span><span class="o">.</span><span class="n">times</span> <span class="o">&lt;</span> <span class="n">cleaned_anas</span><span class="o">.</span><span class="n">times</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">])],</span>
-        <span class="n">t_start</span><span class="o">=</span><span class="n">t_start</span><span class="p">,</span> <span class="n">t_stop</span><span class="o">=</span><span class="n">cleaned_anas</span><span class="o">.</span><span class="n">times</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
-    <span class="p">)</span>
-
-    <span class="n">sigs</span><span class="p">,</span> <span class="n">freqs</span> <span class="o">=</span> <span class="n">el</span><span class="o">.</span><span class="n">sta</span><span class="o">.</span><span class="n">spike_field_coherence</span><span class="p">(</span><span class="n">cleaned_anas</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="o">**</span><span class="p">{</span><span class="s1">&#39;nperseg&#39;</span><span class="p">:</span> <span class="n">NFFT</span><span class="p">})</span>
-    <span class="k">return</span> <span class="n">sigs</span><span class="p">,</span> <span class="n">freqs</span><span class="p">,</span> <span class="n">cleaned_anas</span>
-</pre></div>
-
-    </div>
-</div>
-</div>
-
-</div>
-<div class="cell border-box-sizing code_cell rendered">
-<div class="input">
-<div class="prompt input_prompt">In&nbsp;[13]:</div>
-<div class="inner_cell">
-    <div class="input_area">
-<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># action_id_0 = &#39;1833-200619-1&#39;</span>
-<span class="c1"># action_id_1 = &#39;1833-200619-2&#39;</span>
-<span class="c1"># lfp_0 = data_loader.lfp(action_id_0, 6)</span>
-<span class="c1"># lfp_1 = data_loader.lfp(action_id_1, 6)</span>
-
-<span class="c1"># lim_0 = get_lim(action_id_0)</span>
-<span class="c1"># lim_1 = get_lim(action_id_1)</span>
-
-<span class="c1"># sptrs_0 = data_loader.spike_trains(action_id_0, 6)</span>
-<span class="c1"># sptrs_1 = data_loader.spike_trains(action_id_1, 6)</span>
-
-<span class="c1"># coher_0, freq_0, clean_lfp_0 = compute_spike_lfp(lfp_0, sptrs_0[163], *lim_0, 4096)</span>
-<span class="c1"># coher_1, freq_1, clean_lfp_1 = compute_spike_lfp(lfp_1, sptrs_1[28], *lim_1, 4096)</span>
-
-<span class="c1"># best_channel_0 = np.argmax(signaltonoise(clean_lfp_0))</span>
-<span class="c1"># best_channel_1 = np.argmax(signaltonoise(clean_lfp_1))</span>
-
-<span class="c1"># plt.plot(freq_0, coher_0[:,best_channel_0])</span>
-<span class="c1"># plt.plot(freq_1, coher_1[:,best_channel_1])</span>
-<span class="c1"># plt.xlim(0,20)</span>
-</pre></div>
-
-    </div>
-</div>
-</div>
-
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
 <div class="prompt input_prompt">In&nbsp;[14]:</div>
 <div class="inner_cell">
     <div class="input_area">
-<div class=" highlight hl-ipython3"><pre><span></span><span class="n">NFFT</span> <span class="o">=</span> <span class="mi">4096</span><span class="o">*</span><span class="mi">2</span>
+<div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">compute_spike_lfp_coherence</span><span class="p">(</span><span class="n">anas</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">NFFT</span><span class="p">):</span>
-<span class="n">theta_band_f1</span><span class="p">,</span> <span class="n">theta_band_f2</span> <span class="o">=</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">10</span> 
+
+    <span class="n">sigs</span><span class="p">,</span> <span class="n">freqs</span> <span class="o">=</span> <span class="n">el</span><span class="o">.</span><span class="n">sta</span><span class="o">.</span><span class="n">spike_field_coherence</span><span class="p">(</span><span class="n">anas</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="o">**</span><span class="p">{</span><span class="s1">&#39;nperseg&#39;</span><span class="p">:</span> <span class="n">NFFT</span><span class="p">})</span>
+    <span class="k">return</span> <span class="n">sigs</span><span class="p">,</span> <span class="n">freqs</span>
 </pre></div>
 
     </div>
@@ -13454,23 +13392,327 @@ div#notebook {
 <div class="prompt input_prompt">In&nbsp;[15]:</div>
 <div class="inner_cell">
     <div class="input_area">
-<div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">process</span><span class="p">(</span><span class="n">row</span><span class="p">):</span>
+<div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">butter_bandpass</span><span class="p">(</span><span class="n">lowcut</span><span class="p">,</span> <span class="n">highcut</span><span class="p">,</span> <span class="n">fs</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="mi">5</span><span class="p">):</span>
+    <span class="n">nyq</span> <span class="o">=</span> <span class="mf">0.5</span> <span class="o">*</span> <span class="n">fs</span>
+    <span class="n">low</span> <span class="o">=</span> <span class="n">lowcut</span> <span class="o">/</span> <span class="n">nyq</span>
+    <span class="n">high</span> <span class="o">=</span> <span class="n">highcut</span> <span class="o">/</span> <span class="n">nyq</span>
+    <span class="n">b</span><span class="p">,</span> <span class="n">a</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">signal</span><span class="o">.</span><span class="n">butter</span><span class="p">(</span><span class="n">order</span><span class="p">,</span> <span class="p">[</span><span class="n">low</span><span class="p">,</span> <span class="n">high</span><span class="p">],</span> <span class="n">btype</span><span class="o">=</span><span class="s1">&#39;band&#39;</span><span class="p">)</span>
+    <span class="k">return</span> <span class="n">b</span><span class="p">,</span> <span class="n">a</span>
+
+
+<span class="k">def</span> <span class="nf">butter_bandpass_filter</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">lowcut</span><span class="p">,</span> <span class="n">highcut</span><span class="p">,</span> <span class="n">fs</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="mi">5</span><span class="p">):</span>
+    <span class="n">b</span><span class="p">,</span> <span class="n">a</span> <span class="o">=</span> <span class="n">butter_bandpass</span><span class="p">(</span><span class="n">lowcut</span><span class="p">,</span> <span class="n">highcut</span><span class="p">,</span> <span class="n">fs</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="n">order</span><span class="p">)</span>
+    <span class="n">y</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">signal</span><span class="o">.</span><span class="n">filtfilt</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">data</span><span class="p">)</span>
+    <span class="k">return</span> <span class="n">y</span>
+
+<span class="k">def</span> <span class="nf">compute_spike_phase_func</span><span class="p">(</span><span class="n">lfp</span><span class="p">,</span> <span class="n">times</span><span class="p">,</span> <span class="n">return_degrees</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
+    <span class="n">x_a</span> <span class="o">=</span> <span class="n">hilbert</span><span class="p">(</span><span class="n">lfp</span><span class="p">)</span>
+    <span class="n">x_phase</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">angle</span><span class="p">(</span><span class="n">x_a</span><span class="p">)</span>
+    <span class="k">if</span> <span class="n">return_degrees</span><span class="p">:</span>
+        <span class="n">x_phase</span> <span class="o">=</span> <span class="n">x_phase</span> <span class="o">*</span> <span class="mi">180</span> <span class="o">/</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span>
+    <span class="k">return</span> <span class="n">interp1d</span><span class="p">(</span><span class="n">times</span><span class="p">,</span> <span class="n">x_phase</span><span class="p">)</span>
+
+
+<span class="k">def</span> <span class="nf">vonmises_kde</span><span class="p">(</span><span class="n">data</span><span class="p">,</span> <span class="n">kappa</span><span class="o">=</span><span class="mi">100</span><span class="p">,</span> <span class="n">n_bins</span><span class="o">=</span><span class="mi">100</span><span class="p">):</span>
+    <span class="kn">from</span> <span class="nn">scipy.special</span> <span class="k">import</span> <span class="n">i0</span>
+    <span class="n">bins</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="n">n_bins</span><span class="p">)</span>
+    <span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">,</span> <span class="n">n_bins</span><span class="p">)</span>
+    <span class="c1"># integrate vonmises kernels</span>
+    <span class="n">kde</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">kappa</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">x</span><span class="p">[:,</span> <span class="kc">None</span><span class="p">]</span> <span class="o">-</span> <span class="n">data</span><span class="p">[</span><span class="kc">None</span><span class="p">,</span> <span class="p">:]))</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="mi">1</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="mi">2</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">i0</span><span class="p">(</span><span class="n">kappa</span><span class="p">))</span>
+    <span class="n">kde</span> <span class="o">/=</span> <span class="n">np</span><span class="o">.</span><span class="n">trapz</span><span class="p">(</span><span class="n">kde</span><span class="p">,</span> <span class="n">x</span><span class="o">=</span><span class="n">bins</span><span class="p">)</span>
+    <span class="k">return</span> <span class="n">bins</span><span class="p">,</span> <span class="n">kde</span>
+
+
+<span class="k">def</span> <span class="nf">spike_phase_score</span><span class="p">(</span><span class="n">phase_bins</span><span class="p">,</span> <span class="n">density</span><span class="p">):</span>
+    <span class="kn">import</span> <span class="nn">math</span>
+    <span class="kn">import</span> <span class="nn">pycircstat</span> <span class="k">as</span> <span class="nn">pc</span>
+    <span class="n">ang</span> <span class="o">=</span> <span class="n">pc</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">phase_bins</span><span class="p">,</span> <span class="n">w</span><span class="o">=</span><span class="n">density</span><span class="p">)</span>
+    <span class="n">vec_len</span> <span class="o">=</span> <span class="n">pc</span><span class="o">.</span><span class="n">resultant_vector_length</span><span class="p">(</span><span class="n">phase_bins</span><span class="p">,</span> <span class="n">w</span><span class="o">=</span><span class="n">density</span><span class="p">)</span>
+    <span class="c1"># ci_lim = pc.mean_ci_limits(head_angle_bins, w=rate)</span>
+    <span class="k">return</span> <span class="n">ang</span><span class="p">,</span> <span class="n">vec_len</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[16]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">compute_clean_lfp</span><span class="p">(</span><span class="n">anas</span><span class="p">,</span> <span class="n">width</span><span class="o">=</span><span class="mi">500</span><span class="p">,</span> <span class="n">threshold</span><span class="o">=</span><span class="mi">2</span><span class="p">):</span>
+    <span class="n">anas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">anas</span><span class="p">)</span>
+    <span class="n">idxs</span><span class="p">,</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">where</span><span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">anas</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">threshold</span><span class="p">)</span>
+    <span class="k">for</span> <span class="n">idx</span> <span class="ow">in</span> <span class="n">idxs</span><span class="p">:</span>
+        <span class="n">anas</span><span class="p">[</span><span class="n">idx</span><span class="o">-</span><span class="n">width</span><span class="p">:</span><span class="n">idx</span><span class="o">+</span><span class="n">width</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span> <span class="c1"># TODO AR model prediction</span>
+    <span class="k">return</span> <span class="n">anas</span><span class="p">,</span> <span class="n">idxs</span>
+
+
+<span class="k">def</span> <span class="nf">compute_clean_spikes</span><span class="p">(</span><span class="n">spikes</span><span class="p">,</span> <span class="n">idxs</span><span class="p">,</span> <span class="n">times</span><span class="p">,</span> <span class="n">width</span><span class="o">=</span><span class="mi">500</span><span class="p">):</span>
+
+    <span class="k">for</span> <span class="n">idx</span> <span class="ow">in</span> <span class="n">idxs</span><span class="p">:</span>
+        <span class="n">t0</span> <span class="o">=</span> <span class="n">times</span><span class="p">[</span><span class="n">idx</span><span class="o">-</span><span class="n">width</span><span class="p">]</span>
+        <span class="n">stop</span> <span class="o">=</span> <span class="n">idx</span> <span class="o">+</span> <span class="n">width</span>
+        <span class="k">if</span> <span class="n">stop</span> <span class="o">&gt;</span> <span class="nb">len</span><span class="p">(</span><span class="n">times</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">:</span>
+            <span class="n">stop</span> <span class="o">=</span> <span class="nb">len</span><span class="p">(</span><span class="n">times</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span> 
+        <span class="n">t1</span> <span class="o">=</span> <span class="n">times</span><span class="p">[</span><span class="n">stop</span><span class="p">]</span>
+        <span class="n">mask</span> <span class="o">=</span> <span class="p">(</span><span class="n">spikes</span> <span class="o">&gt;</span> <span class="n">t0</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">spikes</span> <span class="o">&lt;</span> <span class="n">t1</span><span class="p">)</span>
+        <span class="n">spikes</span> <span class="o">=</span> <span class="n">spikes</span><span class="p">[</span><span class="o">~</span><span class="n">mask</span><span class="p">]</span>
+    <span class="n">spikes</span> <span class="o">=</span> <span class="n">spikes</span><span class="p">[</span><span class="n">spikes</span> <span class="o">&lt;=</span> <span class="n">times</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]]</span>
+    <span class="k">return</span> <span class="n">spikes</span>
+
+
+<span class="k">def</span> <span class="nf">prepare_spike_lfp</span><span class="p">(</span><span class="n">anas</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">t_start</span><span class="p">,</span> <span class="n">t_stop</span><span class="p">):</span>
+
+    <span class="n">t_start</span> <span class="o">=</span> <span class="n">t_start</span> <span class="o">*</span> <span class="n">pq</span><span class="o">.</span><span class="n">s</span> <span class="k">if</span> <span class="n">t_start</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span> <span class="k">else</span> <span class="mi">0</span> <span class="o">*</span> <span class="n">pq</span><span class="o">.</span><span class="n">s</span>
+    <span class="n">sampling_rate</span> <span class="o">=</span> <span class="n">anas</span><span class="o">.</span><span class="n">sampling_rate</span>
+    <span class="n">units</span> <span class="o">=</span> <span class="n">anas</span><span class="o">.</span><span class="n">units</span>
+    <span class="n">times</span> <span class="o">=</span> <span class="n">anas</span><span class="o">.</span><span class="n">times</span>
+    <span class="k">if</span> <span class="n">t_start</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span> <span class="ow">and</span> <span class="n">t_stop</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
+        <span class="n">t_stop</span> <span class="o">=</span> <span class="n">t_stop</span> <span class="o">*</span> <span class="n">pq</span><span class="o">.</span><span class="n">s</span>
+        <span class="n">mask</span> <span class="o">=</span> <span class="p">(</span><span class="n">times</span> <span class="o">&gt;</span> <span class="n">t_start</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">times</span> <span class="o">&lt;</span> <span class="n">t_stop</span><span class="p">)</span>
+        <span class="n">anas</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">anas</span><span class="p">)[</span><span class="n">mask</span><span class="p">,:]</span>
+        <span class="n">times</span> <span class="o">=</span> <span class="n">times</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span>
+
+    <span class="n">best_channel</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">argmax</span><span class="p">(</span><span class="n">signaltonoise</span><span class="p">(</span><span class="n">anas</span><span class="p">))</span>
+<span class="c1">#     best_channel = np.random.choice(anas.shape[1])</span>
+    
+    <span class="n">cleaned_anas</span><span class="p">,</span> <span class="n">idxs</span> <span class="o">=</span> <span class="n">compute_clean_lfp</span><span class="p">(</span><span class="n">anas</span><span class="p">[:,</span> <span class="n">best_channel</span><span class="p">])</span>
+    
+    <span class="n">cleaned_anas</span> <span class="o">=</span> <span class="n">neo</span><span class="o">.</span><span class="n">AnalogSignal</span><span class="p">(</span>
+        <span class="n">signal</span><span class="o">=</span><span class="n">cleaned_anas</span> <span class="o">*</span> <span class="n">units</span><span class="p">,</span> <span class="n">sampling_rate</span><span class="o">=</span><span class="n">sampling_rate</span><span class="p">,</span> <span class="n">t_start</span><span class="o">=</span><span class="n">t_start</span>
+    <span class="p">)</span>
+    
+    <span class="n">spike_units</span> <span class="o">=</span> <span class="n">sptr</span><span class="o">.</span><span class="n">units</span>
+    <span class="n">spike_times</span> <span class="o">=</span> <span class="n">sptr</span><span class="o">.</span><span class="n">times</span>
+    <span class="n">spike_times</span> <span class="o">=</span> <span class="n">compute_clean_spikes</span><span class="p">(</span><span class="n">spike_times</span><span class="p">,</span> <span class="n">idxs</span><span class="p">,</span> <span class="n">times</span><span class="p">)</span>
+
+    <span class="n">sptr</span> <span class="o">=</span> <span class="n">neo</span><span class="o">.</span><span class="n">SpikeTrain</span><span class="p">(</span>
+        <span class="n">spike_times</span><span class="p">[(</span><span class="n">spike_times</span> <span class="o">&gt;</span> <span class="n">t_start</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">spike_times</span> <span class="o">&lt;</span> <span class="n">times</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">])],</span> <span class="n">units</span><span class="o">=</span><span class="n">spike_units</span><span class="p">,</span>
+        <span class="n">t_start</span><span class="o">=</span><span class="n">t_start</span><span class="p">,</span> <span class="n">t_stop</span><span class="o">=</span><span class="n">times</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
+    <span class="p">)</span>
+
+    <span class="k">return</span> <span class="n">cleaned_anas</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">best_channel</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[17]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="k">def</span> <span class="nf">compute_spike_phase_func</span><span class="p">(</span><span class="n">lfp</span><span class="p">,</span> <span class="n">times</span><span class="p">,</span> <span class="n">return_degrees</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
+    <span class="kn">from</span> <span class="nn">scipy.fftpack</span> <span class="k">import</span> <span class="n">next_fast_len</span>
+    <span class="n">x_a</span> <span class="o">=</span> <span class="n">scipy</span><span class="o">.</span><span class="n">signal</span><span class="o">.</span><span class="n">hilbert</span><span class="p">(</span>
+        <span class="n">lfp</span><span class="p">,</span> <span class="n">next_fast_len</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">lfp</span><span class="p">)))[:</span><span class="nb">len</span><span class="p">(</span><span class="n">lfp</span><span class="p">)]</span>
+<span class="c1">#     x_a = hilbert(lfp)</span>
+    <span class="n">x_phase</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">angle</span><span class="p">(</span><span class="n">x_a</span><span class="p">,</span> <span class="n">deg</span><span class="o">=</span><span class="n">return_degrees</span><span class="p">)</span>
+    <span class="k">return</span> <span class="n">interp1d</span><span class="p">(</span><span class="n">times</span><span class="p">,</span> <span class="n">x_phase</span><span class="p">)</span>
+
+
+<span class="k">def</span> <span class="nf">compute_spike_phase</span><span class="p">(</span><span class="n">lfp</span><span class="p">,</span> <span class="n">spikes</span><span class="p">,</span> <span class="n">flim</span><span class="o">=</span><span class="p">[</span><span class="mi">6</span><span class="p">,</span><span class="mi">10</span><span class="p">]):</span>
+    
+    <span class="n">sample_rate</span> <span class="o">=</span> <span class="n">lfp</span><span class="o">.</span><span class="n">sampling_rate</span><span class="o">.</span><span class="n">magnitude</span>
+    
+    <span class="c1"># sometimes the position is recorded after LFP recording is ended</span>
+    <span class="n">times</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arange</span><span class="p">(</span><span class="n">lfp</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="o">/</span> <span class="n">sample_rate</span>
+     
+    <span class="n">spikes</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">spikes</span><span class="p">)</span>
+    <span class="n">spikes</span> <span class="o">=</span> <span class="n">spikes</span><span class="p">[(</span><span class="n">spikes</span> <span class="o">&gt;</span> <span class="n">times</span><span class="o">.</span><span class="n">min</span><span class="p">())</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">spikes</span> <span class="o">&lt;</span> <span class="n">times</span><span class="o">.</span><span class="n">max</span><span class="p">())]</span>
+    
+    <span class="n">filtered_lfp</span> <span class="o">=</span> <span class="n">butter_bandpass_filter</span><span class="p">(</span>
+        <span class="n">lfp</span><span class="o">.</span><span class="n">magnitude</span><span class="o">.</span><span class="n">ravel</span><span class="p">(),</span> <span class="o">*</span><span class="n">flim</span><span class="p">,</span> <span class="n">fs</span><span class="o">=</span><span class="n">sample_rate</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
+
+    <span class="n">spike_phase_func</span> <span class="o">=</span> <span class="n">compute_spike_phase_func</span><span class="p">(</span><span class="n">filtered_lfp</span><span class="p">,</span> <span class="n">times</span><span class="p">)</span>
+    
+    <span class="k">return</span> <span class="n">spike_phase_func</span><span class="p">(</span><span class="n">spikes</span><span class="p">),</span> <span class="n">filtered_lfp</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[18]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># plt.figure(figsize=(16,9))</span>
+<span class="c1"># # lfp = data_loader.lfp(&#39;1833-200619-2&#39;, 6)</span>
+<span class="c1"># lfp = data_loader.lfp(&#39;1834-220319-3&#39;, 6)</span>
+<span class="c1"># # lfp = data_loader.lfp(&#39;1849-010319-4&#39;, 6)</span>
+<span class="c1"># times = np.arange(lfp.shape[0]) / lfp.sampling_rate.magnitude</span>
+<span class="c1"># clean_lfp, _ = compute_clean_lfp(lfp.magnitude[:, 0], threshold=2)</span>
+<span class="c1"># plt.plot(times,lfp[:,0])</span>
+<span class="c1"># plt.plot(times,clean_lfp)</span>
+
+<span class="c1"># plt.figure(figsize=(16,9))</span>
+<span class="c1"># plt.psd(lfp[:,0].ravel(), Fs=1000, NFFT=10000)</span>
+<span class="c1"># plt.psd(clean_lfp, Fs=1000, NFFT=10000)</span>
+<span class="c1"># plt.xlim(0,100)</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[19]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># plt.figure(figsize=(16,9))</span>
+
+<span class="c1"># plt.plot(times,lfp[:,0])</span>
+<span class="c1"># # plt.plot(clean_lfp*100)</span>
+<span class="c1"># plt.plot(times[:-1], np.diff(lfp[:,0].magnitude.ravel()))</span>
+<span class="c1"># plt.xlim(64.5,65.5)</span>
+<span class="c1"># # plt.ylim(-250,250)</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[20]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># # action_id_0, channel_0, unit_0 = &#39;1833-200619-1&#39;, 6, 163</span>
+<span class="c1"># # action_id_1, channel_1, unit_1 = &#39;1833-200619-2&#39;, 6, 28</span>
+<span class="c1"># action_id_0, channel_0, unit_0 = &#39;1834-220319-3&#39;, 2, 46</span>
+<span class="c1"># action_id_1, channel_1, unit_1 = &#39;1834-220319-4&#39;, 2, 60</span>
+<span class="c1"># lfp_0 = data_loader.lfp(action_id_0, channel_0)</span>
+<span class="c1"># lfp_1 = data_loader.lfp(action_id_1, channel_1)</span>
+
+<span class="c1"># sample_rate_0 = lfp_0.sampling_rate</span>
+<span class="c1"># sample_rate_1 = lfp_1.sampling_rate</span>
+
+<span class="c1"># lim_0 = get_lim(action_id_0)</span>
+<span class="c1"># lim_1 = get_lim(action_id_1)</span>
+
+<span class="c1"># sptrs_0 = data_loader.spike_trains(action_id_0, channel_0)</span>
+
+<span class="c1"># sptrs_1 = data_loader.spike_trains(action_id_1, channel_1)</span>
+
+<span class="c1"># cleaned_lfps_0, sptr_0, best_channel_0 = prepare_spike_lfp(lfp_0, sptrs_0[unit_0], *lim_0)</span>
+
+<span class="c1"># cleaned_lfps_1, sptr_1, best_channel_1 = prepare_spike_lfp(lfp_1, sptrs_1[unit_1], *lim_1)</span>
+
+<span class="c1"># coher_0, freq_0 = compute_spike_lfp_coherence(cleaned_lfps_0, sptr_0, 4096)</span>
+
+<span class="c1"># coher_1, freq_1 = compute_spike_lfp_coherence(cleaned_lfps_1, sptr_1, 4096)</span>
+
+<span class="c1"># spike_phase_0, filtered_lfp_0 = compute_spike_phase(cleaned_lfps_0, sptrs_0[unit_0], flim=[6,10])</span>
+
+<span class="c1"># spike_phase_1, filtered_lfp_1 = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[6,10])</span>
+
+<span class="c1"># spike_phase_1_stim, filtered_lfp_1_stim = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[10.5,11.5])</span>
+
+<span class="c1"># plt.figure()</span>
+<span class="c1"># plt.plot(freq_0, coher_0.ravel())</span>
+<span class="c1"># plt.plot(freq_1, coher_1.ravel())</span>
+<span class="c1"># plt.xlim(0,20)</span>
+
+<span class="c1"># plt.figure()</span>
+<span class="c1"># bins_0, kde_0 = vonmises_kde(spike_phase_0, 100)</span>
+<span class="c1"># ang_0, vec_len_0 = spike_phase_score(bins_0, kde_0)</span>
+<span class="c1"># plt.polar(bins_0, kde_0, color=&#39;b&#39;)</span>
+<span class="c1"># plt.polar([ang_0, ang_0], [0, vec_len_0], color=&#39;b&#39;)</span>
+
+<span class="c1"># bins_1, kde_1 = vonmises_kde(spike_phase_1, 100)</span>
+<span class="c1"># ang_1, vec_len_1 = spike_phase_score(bins_1, kde_1)</span>
+<span class="c1"># plt.polar(bins_1, kde_1, color=&#39;r&#39;)</span>
+<span class="c1"># plt.polar([ang_1, ang_1], [0, vec_len_1], color=&#39;r&#39;)</span>
+
+<span class="c1"># bins_1_stim, kde_1_stim = vonmises_kde(spike_phase_1_stim, 100)</span>
+<span class="c1"># ang_1_stim, vec_len_1_stim = spike_phase_score(bins_1_stim, kde_1_stim)</span>
+<span class="c1"># plt.polar(bins_1_stim, kde_1_stim, color=&#39;k&#39;)</span>
+<span class="c1"># plt.polar([ang_1_stim, ang_1_stim], [0, vec_len_1_stim], color=&#39;k&#39;)</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[&nbsp;]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span> 
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[21]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># TODO fix artefact stuff from phase precession</span>
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[22]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="n">NFFT</span> <span class="o">=</span> <span class="mi">8192</span>
+<span class="n">theta_band_f1</span><span class="p">,</span> <span class="n">theta_band_f2</span> <span class="o">=</span> <span class="mi">6</span><span class="p">,</span> <span class="mi">10</span> 
+</pre></div>
+
+    </div>
+</div>
+</div>
+
+</div>
+<div class="cell border-box-sizing code_cell rendered">
+<div class="input">
+<div class="prompt input_prompt">In&nbsp;[23]:</div>
+<div class="inner_cell">
+    <div class="input_area">
+<div class=" highlight hl-ipython3"><pre><span></span><span class="n">coherence_data</span><span class="p">,</span> <span class="n">freqency_data</span> <span class="o">=</span> <span class="p">{},</span> <span class="p">{}</span>
+<span class="n">theta_kde_data</span><span class="p">,</span> <span class="n">theta_bins_data</span> <span class="o">=</span> <span class="p">{},</span> <span class="p">{}</span>
+<span class="n">stim_kde_data</span><span class="p">,</span> <span class="n">stim_bins_data</span> <span class="o">=</span> <span class="p">{},</span> <span class="p">{}</span>
+
+<span class="k">def</span> <span class="nf">process</span><span class="p">(</span><span class="n">row</span><span class="p">):</span>
     <span class="n">action_id</span> <span class="o">=</span> <span class="n">row</span><span class="p">[</span><span class="s1">&#39;action&#39;</span><span class="p">]</span>
     <span class="n">channel_group</span> <span class="o">=</span> <span class="n">row</span><span class="p">[</span><span class="s1">&#39;channel_group&#39;</span><span class="p">]</span>
     <span class="n">unit_name</span> <span class="o">=</span> <span class="n">row</span><span class="p">[</span><span class="s1">&#39;unit_name&#39;</span><span class="p">]</span>
     
+    <span class="n">name</span> <span class="o">=</span> <span class="n">f</span><span class="s1">&#39;</span><span class="si">{action_id}</span><span class="s1">_</span><span class="si">{channel_group}</span><span class="s1">_</span><span class="si">{unit_name}</span><span class="s1">&#39;</span>
+    
     <span class="n">lfp</span> <span class="o">=</span> <span class="n">data_loader</span><span class="o">.</span><span class="n">lfp</span><span class="p">(</span><span class="n">action_id</span><span class="p">,</span> <span class="n">channel_group</span><span class="p">)</span> <span class="c1"># TODO consider choosing strongest stim response</span>
     
     <span class="n">sptr</span> <span class="o">=</span> <span class="n">data_loader</span><span class="o">.</span><span class="n">spike_train</span><span class="p">(</span><span class="n">action_id</span><span class="p">,</span> <span class="n">channel_group</span><span class="p">,</span> <span class="n">unit_name</span><span class="p">)</span>
        
     <span class="n">lim</span> <span class="o">=</span> <span class="n">get_lim</span><span class="p">(</span><span class="n">action_id</span><span class="p">)</span>
-
-    <span class="n">p_xys</span><span class="p">,</span> <span class="n">freq</span><span class="p">,</span> <span class="n">clean_lfp</span> <span class="o">=</span> <span class="n">compute_spike_lfp</span><span class="p">(</span><span class="n">lfp</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="o">*</span><span class="n">lim</span><span class="p">,</span> <span class="n">NFFT</span><span class="o">=</span><span class="n">NFFT</span><span class="p">)</span>
     
-    <span class="n">snls</span> <span class="o">=</span> <span class="n">signaltonoise</span><span class="p">(</span><span class="n">clean_lfp</span><span class="p">)</span>
+    <span class="n">cleaned_lfp</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">best_channel</span> <span class="o">=</span> <span class="n">prepare_spike_lfp</span><span class="p">(</span><span class="n">lfp</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="o">*</span><span class="n">lim</span><span class="p">)</span>
-    <span class="n">best_channel</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">argmax</span><span class="p">(</span><span class="n">snls</span><span class="p">)</span>
+       
-    <span class="n">snl</span> <span class="o">=</span> <span class="n">snls</span><span class="p">[</span><span class="n">best_channel</span><span class="p">]</span>
+    <span class="n">p_xys</span><span class="p">,</span> <span class="n">freq</span> <span class="o">=</span> <span class="n">compute_spike_lfp_coherence</span><span class="p">(</span><span class="n">cleaned_lfp</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">NFFT</span><span class="o">=</span><span class="n">NFFT</span><span class="p">)</span>
-    <span class="n">p_xy</span> <span class="o">=</span> <span class="n">p_xys</span><span class="p">[:,</span><span class="n">best_channel</span><span class="p">]</span><span class="o">.</span><span class="n">magnitude</span>
+    
+    <span class="n">p_xy</span> <span class="o">=</span> <span class="n">p_xys</span><span class="o">.</span><span class="n">magnitude</span><span class="o">.</span><span class="n">ravel</span><span class="p">()</span>
     <span class="n">freq</span> <span class="o">=</span> <span class="n">freq</span><span class="o">.</span><span class="n">magnitude</span>
     
     <span class="n">theta_f</span><span class="p">,</span> <span class="n">theta_p_max</span> <span class="o">=</span> <span class="n">find_theta_peak</span><span class="p">(</span><span class="n">p_xy</span><span class="p">,</span> <span class="n">freq</span><span class="p">,</span> <span class="n">theta_band_f1</span><span class="p">,</span> <span class="n">theta_band_f2</span><span class="p">)</span>
@@ -13483,6 +13725,12 @@ div#notebook {
     
     <span class="n">theta_half_energy</span> <span class="o">=</span> <span class="n">compute_energy</span><span class="p">(</span><span class="n">p_xy</span><span class="p">,</span> <span class="n">freq</span><span class="p">,</span> <span class="n">theta_half_f1</span><span class="p">,</span> <span class="n">theta_half_f2</span><span class="p">)</span> <span class="c1"># theta band 6 - 10 Hz</span>
     
+    <span class="n">theta_spike_phase</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">compute_spike_phase</span><span class="p">(</span><span class="n">cleaned_lfp</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">flim</span><span class="o">=</span><span class="p">[</span><span class="n">theta_band_f1</span><span class="p">,</span> <span class="n">theta_band_f2</span><span class="p">])</span>
+    <span class="n">theta_bins</span><span class="p">,</span> <span class="n">theta_kde</span> <span class="o">=</span> <span class="n">vonmises_kde</span><span class="p">(</span><span class="n">theta_spike_phase</span><span class="p">)</span>
+    <span class="n">theta_ang</span><span class="p">,</span> <span class="n">theta_vec_len</span> <span class="o">=</span> <span class="n">spike_phase_score</span><span class="p">(</span><span class="n">theta_bins</span><span class="p">,</span> <span class="n">theta_kde</span><span class="p">)</span>
+    <span class="n">theta_kde_data</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">theta_kde</span><span class="p">})</span>
+    <span class="n">theta_bins_data</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">theta_bins</span><span class="p">})</span>
+
     <span class="c1"># stim</span>
     
     <span class="n">stim_freq</span> <span class="o">=</span> <span class="n">compute_stim_freq</span><span class="p">(</span><span class="n">action_id</span><span class="p">)</span>
@@ -13494,8 +13742,19 @@ div#notebook {
     
     <span class="n">stim_energy</span> <span class="o">=</span> <span class="n">compute_energy</span><span class="p">(</span><span class="n">p_xy</span><span class="p">,</span> <span class="n">freq</span><span class="p">,</span> <span class="n">stim_half_f1</span><span class="p">,</span> <span class="n">stim_half_f2</span><span class="p">)</span>
     
+    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">stim_freq</span><span class="p">):</span>
+        <span class="n">stim_spike_phase</span><span class="p">,</span> <span class="n">stim_bins</span><span class="p">,</span> <span class="n">stim_kde</span><span class="p">,</span> <span class="n">stim_ang</span><span class="p">,</span> <span class="n">stim_vec_len</span> <span class="o">=</span> <span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">nan</span><span class="p">]</span> <span class="o">*</span> <span class="mi">5</span>
+    <span class="k">else</span><span class="p">:</span>
+        <span class="n">stim_spike_phase</span><span class="p">,</span> <span class="n">_</span> <span class="o">=</span> <span class="n">compute_spike_phase</span><span class="p">(</span><span class="n">cleaned_lfp</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="n">flim</span><span class="o">=</span><span class="p">[</span><span class="n">stim_freq</span> <span class="o">-</span> <span class="o">.</span><span class="mi">5</span><span class="p">,</span> <span class="n">stim_freq</span> <span class="o">+</span> <span class="o">.</span><span class="mi">5</span><span class="p">])</span>
+        <span class="n">stim_bins</span><span class="p">,</span> <span class="n">stim_kde</span> <span class="o">=</span> <span class="n">vonmises_kde</span><span class="p">(</span><span class="n">stim_spike_phase</span><span class="p">)</span>
+        <span class="n">stim_ang</span><span class="p">,</span> <span class="n">stim_vec_len</span> <span class="o">=</span> <span class="n">spike_phase_score</span><span class="p">(</span><span class="n">stim_bins</span><span class="p">,</span> <span class="n">stim_kde</span><span class="p">)</span>
+        <span class="n">stim_kde_data</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">stim_kde</span><span class="p">})</span>
+        <span class="n">stim_bins_data</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">stim_bins</span><span class="p">})</span>
+        
+    <span class="n">coherence_data</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">p_xy</span><span class="p">})</span>
+    <span class="n">freqency_data</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">freq</span><span class="p">})</span>
+    
     <span class="n">result</span> <span class="o">=</span> <span class="n">pd</span><span class="o">.</span><span class="n">Series</span><span class="p">({</span>
-        <span class="s1">&#39;signal_to_noise&#39;</span><span class="p">:</span> <span class="n">snl</span><span class="p">,</span>
         <span class="s1">&#39;best_channel&#39;</span><span class="p">:</span> <span class="n">best_channel</span><span class="p">,</span>
         <span class="s1">&#39;theta_freq&#39;</span><span class="p">:</span> <span class="n">theta_f</span><span class="p">,</span>
         <span class="s1">&#39;theta_peak&#39;</span><span class="p">:</span> <span class="n">theta_p_max</span><span class="p">,</span>
@@ -13504,12 +13763,16 @@ div#notebook {
         <span class="s1">&#39;theta_half_f2&#39;</span><span class="p">:</span> <span class="n">theta_half_f2</span><span class="p">,</span>
         <span class="s1">&#39;theta_half_width&#39;</span><span class="p">:</span> <span class="n">theta_half_width</span><span class="p">,</span>
         <span class="s1">&#39;theta_half_energy&#39;</span><span class="p">:</span> <span class="n">theta_half_energy</span><span class="p">,</span>
+        <span class="s1">&#39;theta_ang&#39;</span><span class="p">:</span> <span class="n">theta_ang</span><span class="p">,</span> 
+        <span class="s1">&#39;theta_vec_len&#39;</span><span class="p">:</span> <span class="n">theta_vec_len</span><span class="p">,</span>
         <span class="s1">&#39;stim_freq&#39;</span><span class="p">:</span> <span class="n">stim_freq</span><span class="p">,</span>
         <span class="s1">&#39;stim_p_max&#39;</span><span class="p">:</span> <span class="n">stim_p_max</span><span class="p">,</span>
         <span class="s1">&#39;stim_half_f1&#39;</span><span class="p">:</span> <span class="n">stim_half_f1</span><span class="p">,</span> 
         <span class="s1">&#39;stim_half_f2&#39;</span><span class="p">:</span> <span class="n">stim_half_f2</span><span class="p">,</span>
         <span class="s1">&#39;stim_half_width&#39;</span><span class="p">:</span> <span class="n">stim_half_width</span><span class="p">,</span>
-        <span class="s1">&#39;stim_energy&#39;</span><span class="p">:</span> <span class="n">stim_energy</span>
+        <span class="s1">&#39;stim_energy&#39;</span><span class="p">:</span> <span class="n">stim_energy</span><span class="p">,</span>
+        <span class="s1">&#39;stim_ang&#39;</span><span class="p">:</span> <span class="n">stim_ang</span><span class="p">,</span> 
+        <span class="s1">&#39;stim_vec_len&#39;</span><span class="p">:</span> <span class="n">stim_vec_len</span>
     <span class="p">})</span>
     <span class="k">return</span> <span class="n">result</span>
 </pre></div>
@@ -13521,7 +13784,7 @@ div#notebook {
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[16]:</div>
+<div class="prompt input_prompt">In&nbsp;[24]:</div>
 <div class="inner_cell">
     <div class="input_area">
 <div class=" highlight hl-ipython3"><pre><span></span><span class="n">results</span> <span class="o">=</span> <span class="n">units</span><span class="o">.</span><span class="n">merge</span><span class="p">(</span>
@@ -13547,13 +13810,13 @@ div#notebook {
 
  
  
-<div id="fdf01411-110c-4f09-8a74-86966ed09bad"></div>
+<div id="a0803bf8-bad9-4783-b06a-316f09ff8583"></div>
 <div class="output_subarea output_widget_view ">
 <script type="text/javascript">
-var element = $('#fdf01411-110c-4f09-8a74-86966ed09bad');
+var element = $('#a0803bf8-bad9-4783-b06a-316f09ff8583');
 </script>
 <script type="application/vnd.jupyter.widget-view+json">
-{"model_id": "564d1ef5339f46eebee42ec41cfcfe62", "version_major": 2, "version_minor": 0}
+{"model_id": "199c002aa9604451a318e3f471c9d892", "version_major": 2, "version_minor": 0}
 </script>
 </div>
 
@@ -13565,8 +13828,9 @@ var element = $('#fdf01411-110c-4f09-8a74-86966ed09bad');
 
 
 <div class="output_subarea output_stream output_stderr output_text">
-<pre>/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/signal/spectral.py:1578: RuntimeWarning: invalid value encountered in true_divide
+<pre>/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/signal/spectral.py:1577: RuntimeWarning: invalid value encountered in true_divide
   Cxy = np.abs(Pxy)**2 / Pxx / Pyy
+/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in true_divide
 </pre>
 </div>
 </div>
@@ -13585,97 +13849,54 @@ var element = $('#fdf01411-110c-4f09-8a74-86966ed09bad');
 </div>
 </div>
 
-</div>
-<div class="cell border-box-sizing text_cell rendered"><div class="prompt input_prompt">
-</div><div class="inner_cell">
-<div class="text_cell_render border-box-sizing rendered_html">
-<h1 id="plot">plot<a class="anchor-link" href="#plot">&#182;</a></h1>
-</div>
-</div>
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[17]:</div>
+<div class="prompt input_prompt">In&nbsp;[25]:</div>
 <div class="inner_cell">
     <div class="input_area">
-<div class=" highlight hl-ipython3"><pre><span></span><span class="n">coher</span><span class="p">,</span> <span class="n">freqs</span> <span class="o">=</span> <span class="p">{},</span> <span class="p">{}</span>
+<div class=" highlight hl-ipython3"><pre><span></span><span class="c1"># coher, freqs = {}, {}</span>
-<span class="k">for</span> <span class="n">i</span><span class="p">,</span> <span class="n">row</span> <span class="ow">in</span> <span class="n">tqdm</span><span class="p">(</span><span class="n">units</span><span class="o">.</span><span class="n">iterrows</span><span class="p">(),</span> <span class="n">total</span><span class="o">=</span><span class="nb">len</span><span class="p">(</span><span class="n">units</span><span class="p">)):</span>
+<span class="c1"># for i, row in tqdm(units.iterrows(), total=len(units)):</span>
-    <span class="n">action_id</span> <span class="o">=</span> <span class="n">row</span><span class="p">[</span><span class="s1">&#39;action&#39;</span><span class="p">]</span>
+<span class="c1">#     action_id = row[&#39;action&#39;]</span>
-    <span class="n">channel_group</span> <span class="o">=</span> <span class="n">row</span><span class="p">[</span><span class="s1">&#39;channel_group&#39;</span><span class="p">]</span>
+<span class="c1">#     channel_group = row[&#39;channel_group&#39;]</span>
-    <span class="n">unit_name</span> <span class="o">=</span> <span class="n">row</span><span class="p">[</span><span class="s1">&#39;unit_name&#39;</span><span class="p">]</span>
+<span class="c1">#     unit_name = row[&#39;unit_name&#39;]</span>
     
-    <span class="n">name</span> <span class="o">=</span> <span class="n">f</span><span class="s1">&#39;</span><span class="si">{action_id}</span><span class="s1">_</span><span class="si">{channel_group}</span><span class="s1">_</span><span class="si">{unit_name}</span><span class="s1">&#39;</span>
+<span class="c1">#     name = f&#39;{action_id}_{channel_group}_{unit_name}&#39;</span>
     
-    <span class="n">lfp</span> <span class="o">=</span> <span class="n">data_loader</span><span class="o">.</span><span class="n">lfp</span><span class="p">(</span><span class="n">action_id</span><span class="p">,</span> <span class="n">channel_group</span><span class="p">)</span> <span class="c1"># TODO consider choosing strongest stim response</span>
+<span class="c1">#     lfp = data_loader.lfp(action_id, channel_group) # TODO consider choosing strongest stim response</span>
     
-    <span class="n">sptr</span> <span class="o">=</span> <span class="n">data_loader</span><span class="o">.</span><span class="n">spike_train</span><span class="p">(</span><span class="n">action_id</span><span class="p">,</span> <span class="n">channel_group</span><span class="p">,</span> <span class="n">unit_name</span><span class="p">)</span>
+<span class="c1">#     sptr = data_loader.spike_train(action_id, channel_group, unit_name)</span>
        
-    <span class="n">lim</span> <span class="o">=</span> <span class="n">get_lim</span><span class="p">(</span><span class="n">action_id</span><span class="p">)</span>
+<span class="c1">#     lim = get_lim(action_id)</span>
 
-    <span class="n">p_xys</span><span class="p">,</span> <span class="n">freq</span><span class="p">,</span> <span class="n">clean_lfp</span> <span class="o">=</span> <span class="n">compute_spike_lfp</span><span class="p">(</span><span class="n">lfp</span><span class="p">,</span> <span class="n">sptr</span><span class="p">,</span> <span class="o">*</span><span class="n">lim</span><span class="p">,</span> <span class="n">NFFT</span><span class="o">=</span><span class="n">NFFT</span><span class="p">)</span>
+<span class="c1">#     p_xys, freq, clean_lfp = compute_spike_lfp(lfp, sptr, *lim, NFFT=NFFT)</span>
     
-    <span class="n">snls</span> <span class="o">=</span> <span class="n">signaltonoise</span><span class="p">(</span><span class="n">clean_lfp</span><span class="p">)</span>
+<span class="c1">#     snls = signaltonoise(clean_lfp)</span>
-    <span class="n">best_channel</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">argmax</span><span class="p">(</span><span class="n">snls</span><span class="p">)</span>
+<span class="c1">#     best_channel = np.argmax(snls)</span>
-    <span class="n">snl</span> <span class="o">=</span> <span class="n">snls</span><span class="p">[</span><span class="n">best_channel</span><span class="p">]</span>
+<span class="c1">#     snl = snls[best_channel]</span>
-    <span class="n">p_xy</span> <span class="o">=</span> <span class="n">p_xys</span><span class="p">[:,</span><span class="n">best_channel</span><span class="p">]</span><span class="o">.</span><span class="n">magnitude</span>
+<span class="c1">#     p_xy = p_xys[:,best_channel].magnitude</span>
-    <span class="n">freq</span> <span class="o">=</span> <span class="n">freq</span><span class="o">.</span><span class="n">magnitude</span>
+<span class="c1">#     freq = freq.magnitude</span>
     
-    <span class="n">coher</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">p_xy</span><span class="p">})</span>
+<span class="c1">#     coher.update({name: p_xy})</span>
-    <span class="n">freqs</span><span class="o">.</span><span class="n">update</span><span class="p">({</span><span class="n">name</span><span class="p">:</span> <span class="n">freq</span><span class="p">})</span>
+<span class="c1">#     freqs.update({name: freq})</span>
 </pre></div>
 
     </div>
 </div>
 </div>
 
-<div class="output_wrapper">
-<div class="output">
-
-
-<div class="output_area">
-
-    <div class="prompt"></div>
-
-
-
-
-
- 
- 
-<div id="0aa9c06f-773a-4f88-a2ad-e7de5c5ebe40"></div>
-<div class="output_subarea output_widget_view ">
-<script type="text/javascript">
-var element = $('#0aa9c06f-773a-4f88-a2ad-e7de5c5ebe40');
-</script>
-<script type="application/vnd.jupyter.widget-view+json">
-{"model_id": "b20a9a26e0864f578e1a9aa0c021999b", "version_major": 2, "version_minor": 0}
-</script>
-</div>
-
-</div>
-
-<div class="output_area">
-
-    <div class="prompt"></div>
-
-
-<div class="output_subarea output_stream output_stdout output_text">
-<pre>
-</pre>
-</div>
-</div>
-
-</div>
-</div>
-
 </div>
 <div class="cell border-box-sizing code_cell rendered">
 <div class="input">
-<div class="prompt input_prompt">In&nbsp;[18]:</div>
+<div class="prompt input_prompt">In&nbsp;[26]:</div>
 <div class="inner_cell">
     <div class="input_area">
-<div class=" highlight hl-ipython3"><pre><span></span><span class="n">coher</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;coherence.feather&#39;</span><span class="p">)</span>
+<div class=" highlight hl-ipython3"><pre><span></span><span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">coherence_data</span><span class="p">)</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;coherence.feather&#39;</span><span class="p">)</span>
-<span class="n">freqs</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;freqs.feather&#39;</span><span class="p">)</span>
+<span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">frequency_data</span><span class="p">)</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;freqs.feather&#39;</span><span class="p">)</span>
+<span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">theta_kde_data</span><span class="p">)</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;theta_kde.feather&#39;</span><span class="p">)</span>
+<span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">theta_bins_data</span><span class="p">)</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;theta_bins.feather&#39;</span><span class="p">)</span>
+<span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">stim_kde_data</span><span class="p">)</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;stim_kde.feather&#39;</span><span class="p">)</span>
+<span class="n">pd</span><span class="o">.</span><span class="n">DataFrame</span><span class="p">(</span><span class="n">stim_bins_data</span><span class="p">)</span><span class="o">.</span><span class="n">to_feather</span><span class="p">(</span><span class="n">output</span> <span class="o">/</span> <span class="s1">&#39;data&#39;</span> <span class="o">/</span> <span class="s1">&#39;stim_bins.feather&#39;</span><span class="p">)</span>
 </pre></div>
 
     </div>
@@ -13694,12 +13915,15 @@ var element = $('#0aa9c06f-773a-4f88-a2ad-e7de5c5ebe40');
 <div class="output_subarea output_text output_error">
 <pre>
 <span class="ansi-red-fg">---------------------------------------------------------------------------</span>
-<span class="ansi-red-fg">AttributeError</span>                            Traceback (most recent call last)
+<span class="ansi-red-fg">NameError</span>                                 Traceback (most recent call last)
-<span class="ansi-green-fg">&lt;ipython-input-18-bf8a91a90c5a&gt;</span> in <span class="ansi-cyan-fg">&lt;module&gt;</span>
+<span class="ansi-green-fg">&lt;ipython-input-26-64114f1ea8bd&gt;</span> in <span class="ansi-cyan-fg">&lt;module&gt;</span>
-<span class="ansi-green-fg">----&gt; 1</span><span class="ansi-red-fg"> </span>coher<span class="ansi-blue-fg">.</span>to_feather<span class="ansi-blue-fg">(</span>output <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;data&#39;</span> <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;coherence.feather&#39;</span><span class="ansi-blue-fg">)</span>
+<span class="ansi-green-intense-fg ansi-bold">      1</span> pd<span class="ansi-blue-fg">.</span>DataFrame<span class="ansi-blue-fg">(</span>coherence_data<span class="ansi-blue-fg">)</span><span class="ansi-blue-fg">.</span>to_feather<span class="ansi-blue-fg">(</span>output <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;data&#39;</span> <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;coherence.feather&#39;</span><span class="ansi-blue-fg">)</span>
-<span class="ansi-green-intense-fg ansi-bold">      2</span> freqs<span class="ansi-blue-fg">.</span>to_feather<span class="ansi-blue-fg">(</span>output <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;data&#39;</span> <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;freqs.feather&#39;</span><span class="ansi-blue-fg">)</span>
+<span class="ansi-green-fg">----&gt; 2</span><span class="ansi-red-fg"> </span>pd<span class="ansi-blue-fg">.</span>DataFrame<span class="ansi-blue-fg">(</span>frequency_data<span class="ansi-blue-fg">)</span><span class="ansi-blue-fg">.</span>to_feather<span class="ansi-blue-fg">(</span>output <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;data&#39;</span> <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;freqs.feather&#39;</span><span class="ansi-blue-fg">)</span>
+<span class="ansi-green-intense-fg ansi-bold">      3</span> pd<span class="ansi-blue-fg">.</span>DataFrame<span class="ansi-blue-fg">(</span>theta_kde_data<span class="ansi-blue-fg">)</span><span class="ansi-blue-fg">.</span>to_feather<span class="ansi-blue-fg">(</span>output <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;data&#39;</span> <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;theta_kde.feather&#39;</span><span class="ansi-blue-fg">)</span>
+<span class="ansi-green-intense-fg ansi-bold">      4</span> pd<span class="ansi-blue-fg">.</span>DataFrame<span class="ansi-blue-fg">(</span>theta_bins_data<span class="ansi-blue-fg">)</span><span class="ansi-blue-fg">.</span>to_feather<span class="ansi-blue-fg">(</span>output <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;data&#39;</span> <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;theta_bins.feather&#39;</span><span class="ansi-blue-fg">)</span>
+<span class="ansi-green-intense-fg ansi-bold">      5</span> pd<span class="ansi-blue-fg">.</span>DataFrame<span class="ansi-blue-fg">(</span>stim_kde_data<span class="ansi-blue-fg">)</span><span class="ansi-blue-fg">.</span>to_feather<span class="ansi-blue-fg">(</span>output <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;data&#39;</span> <span class="ansi-blue-fg">/</span> <span class="ansi-blue-fg">&#39;stim_kde.feather&#39;</span><span class="ansi-blue-fg">)</span>
 
-<span class="ansi-red-fg">AttributeError</span>: &#39;dict&#39; object has no attribute &#39;to_feather&#39;</pre>
+<span class="ansi-red-fg">NameError</span>: name &#39;frequency_data&#39; is not defined</pre>
 </div>
 </div>
 

actions/stimulus-spike-lfp-response/data/10-calculate-stimulus-spike-lfp-response.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -12,17 +12,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
+   "outputs": [],
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "10:25:26 [I] klustakwik KlustaKwik2 version 0.2.6\n"
-     ]
-    }
-   ],
    "source": [
     "import matplotlib.pyplot as plt\n",
     "%matplotlib inline\n",
@@ -32,6 +24,8 @@
     "import expipe\n",
     "import os\n",
     "import pathlib\n",
+    "import scipy\n",
+    "import scipy.signal\n",
     "import numpy as np\n",
     "import exdir\n",
     "import pandas as pd\n",
@@ -52,7 +46,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -63,7 +57,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -73,7 +67,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -84,7 +78,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -111,7 +105,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -119,22 +113,12 @@
     "    a = np.asanyarray(a)\n",
     "    m = a.mean(axis)\n",
     "    sd = a.std(axis=axis, ddof=ddof)\n",
-    "    return np.where(sd == 0, 0, m / sd)\n",
+    "    return np.where(sd == 0, 0, m / sd)"
-    "\n",
-    "\n",
-    "def compute_clean_lfp(anas, width=500, threshold=2):\n",
-    "    anas = np.array(anas)\n",
-    "\n",
-    "    for ch in range(anas.shape[1]):\n",
-    "        idxs, = np.where(abs(anas[:, ch]) > threshold)\n",
-    "        for idx in idxs:\n",
-    "            anas[idx-width:idx+width, ch] = 0 # TODO AR model prediction\n",
-    "    return anas"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -148,7 +132,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -165,7 +149,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -191,7 +175,28 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# p = np.load('debug_p.npy')\n",
+    "# f = np.load('debug_f.npy')\n",
+    "# compute_half_width(p, f, 0.01038941, 30.30187709636872)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plt.plot(f, p)\n",
+    "# plt.xlim(29.9,30.6)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -201,73 +206,16 @@
     "    return interp1d(f, p)(s_f)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 12,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def compute_spike_lfp(anas, sptr, t_start, t_stop, NFFT):\n",
-    "\n",
-    "    t_start = t_start * pq.s if t_start is not None else 0 * pq.s\n",
-    "    sampling_rate = anas.sampling_rate\n",
-    "    units = anas.units\n",
-    "    if t_start is not None and t_stop is not None:\n",
-    "        t_stop = t_stop * pq.s\n",
-    "        mask = (anas.times > t_start) & (anas.times < t_stop)\n",
-    "        anas = np.array(anas)[mask,:]\n",
-    "\n",
-    "    cleaned_anas = zscore(compute_clean_lfp(anas))\n",
-    "    \n",
-    "    cleaned_anas = neo.AnalogSignal(\n",
-    "        signal=cleaned_anas * units, sampling_rate=sampling_rate, t_start=t_start\n",
-    "    )\n",
-    "\n",
-    "    sptr = neo.SpikeTrain(\n",
-    "        sptr.times[(sptr.times > t_start) & (sptr.times < cleaned_anas.times[-1])],\n",
-    "        t_start=t_start, t_stop=cleaned_anas.times[-1]\n",
-    "    )\n",
-    "\n",
-    "    sigs, freqs = el.sta.spike_field_coherence(cleaned_anas, sptr, **{'nperseg': NFFT})\n",
-    "    return sigs, freqs, cleaned_anas"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 13,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# action_id_0 = '1833-200619-1'\n",
-    "# action_id_1 = '1833-200619-2'\n",
-    "# lfp_0 = data_loader.lfp(action_id_0, 6)\n",
-    "# lfp_1 = data_loader.lfp(action_id_1, 6)\n",
-    "\n",
-    "# lim_0 = get_lim(action_id_0)\n",
-    "# lim_1 = get_lim(action_id_1)\n",
-    "\n",
-    "# sptrs_0 = data_loader.spike_trains(action_id_0, 6)\n",
-    "# sptrs_1 = data_loader.spike_trains(action_id_1, 6)\n",
-    "\n",
-    "# coher_0, freq_0, clean_lfp_0 = compute_spike_lfp(lfp_0, sptrs_0[163], *lim_0, 4096)\n",
-    "# coher_1, freq_1, clean_lfp_1 = compute_spike_lfp(lfp_1, sptrs_1[28], *lim_1, 4096)\n",
-    "\n",
-    "# best_channel_0 = np.argmax(signaltonoise(clean_lfp_0))\n",
-    "# best_channel_1 = np.argmax(signaltonoise(clean_lfp_1))\n",
-    "\n",
-    "# plt.plot(freq_0, coher_0[:,best_channel_0])\n",
-    "# plt.plot(freq_1, coher_1[:,best_channel_1])\n",
-    "# plt.xlim(0,20)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 14,
    "metadata": {},
    "outputs": [],
    "source": [
-    "NFFT = 4096*2\n",
+    "def compute_spike_lfp_coherence(anas, sptr, NFFT):\n",
-    "theta_band_f1, theta_band_f2 = 6, 10 "
+    "\n",
+    "    sigs, freqs = el.sta.spike_field_coherence(anas, sptr, **{'nperseg': NFFT})\n",
+    "    return sigs, freqs"
    ]
   },
   {
@@ -276,23 +224,289 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "def butter_bandpass(lowcut, highcut, fs, order=5):\n",
+    "    nyq = 0.5 * fs\n",
+    "    low = lowcut / nyq\n",
+    "    high = highcut / nyq\n",
+    "    b, a = scipy.signal.butter(order, [low, high], btype='band')\n",
+    "    return b, a\n",
+    "\n",
+    "\n",
+    "def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):\n",
+    "    b, a = butter_bandpass(lowcut, highcut, fs, order=order)\n",
+    "    y = scipy.signal.filtfilt(b, a, data)\n",
+    "    return y\n",
+    "\n",
+    "def compute_spike_phase_func(lfp, times, return_degrees=False):\n",
+    "    x_a = hilbert(lfp)\n",
+    "    x_phase = np.angle(x_a)\n",
+    "    if return_degrees:\n",
+    "        x_phase = x_phase * 180 / np.pi\n",
+    "    return interp1d(times, x_phase)\n",
+    "\n",
+    "\n",
+    "def vonmises_kde(data, kappa=100, n_bins=100):\n",
+    "    from scipy.special import i0\n",
+    "    bins = np.linspace(-np.pi, np.pi, n_bins)\n",
+    "    x = np.linspace(-np.pi, np.pi, n_bins)\n",
+    "    # integrate vonmises kernels\n",
+    "    kde = np.exp(kappa * np.cos(x[:, None] - data[None, :])).sum(1) / (2 * np.pi * i0(kappa))\n",
+    "    kde /= np.trapz(kde, x=bins)\n",
+    "    return bins, kde\n",
+    "\n",
+    "\n",
+    "def spike_phase_score(phase_bins, density):\n",
+    "    import math\n",
+    "    import pycircstat as pc\n",
+    "    ang = pc.mean(phase_bins, w=density)\n",
+    "    vec_len = pc.resultant_vector_length(phase_bins, w=density)\n",
+    "    # ci_lim = pc.mean_ci_limits(head_angle_bins, w=rate)\n",
+    "    return ang, vec_len"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def compute_clean_lfp(anas, width=500, threshold=2):\n",
+    "    anas = np.array(anas)\n",
+    "    idxs, = np.where(abs(anas) > threshold)\n",
+    "    for idx in idxs:\n",
+    "        anas[idx-width:idx+width] = 0 # TODO AR model prediction\n",
+    "    return anas, idxs\n",
+    "\n",
+    "\n",
+    "def compute_clean_spikes(spikes, idxs, times, width=500):\n",
+    "\n",
+    "    for idx in idxs:\n",
+    "        t0 = times[idx-width]\n",
+    "        stop = idx + width\n",
+    "        if stop > len(times) - 1:\n",
+    "            stop = len(times) - 1 \n",
+    "        t1 = times[stop]\n",
+    "        mask = (spikes > t0) & (spikes < t1)\n",
+    "        spikes = spikes[~mask]\n",
+    "    spikes = spikes[spikes <= times[-1]]\n",
+    "    return spikes\n",
+    "\n",
+    "\n",
+    "def prepare_spike_lfp(anas, sptr, t_start, t_stop):\n",
+    "\n",
+    "    t_start = t_start * pq.s if t_start is not None else 0 * pq.s\n",
+    "    sampling_rate = anas.sampling_rate\n",
+    "    units = anas.units\n",
+    "    times = anas.times\n",
+    "    if t_start is not None and t_stop is not None:\n",
+    "        t_stop = t_stop * pq.s\n",
+    "        mask = (times > t_start) & (times < t_stop)\n",
+    "        anas = np.array(anas)[mask,:]\n",
+    "        times = times[mask]\n",
+    "\n",
+    "    best_channel = np.argmax(signaltonoise(anas))\n",
+    "#     best_channel = np.random.choice(anas.shape[1])\n",
+    "    \n",
+    "    cleaned_anas, idxs = compute_clean_lfp(anas[:, best_channel])\n",
+    "    \n",
+    "    cleaned_anas = neo.AnalogSignal(\n",
+    "        signal=cleaned_anas * units, sampling_rate=sampling_rate, t_start=t_start\n",
+    "    )\n",
+    "    \n",
+    "    spike_units = sptr.units\n",
+    "    spike_times = sptr.times\n",
+    "    spike_times = compute_clean_spikes(spike_times, idxs, times)\n",
+    "\n",
+    "    sptr = neo.SpikeTrain(\n",
+    "        spike_times[(spike_times > t_start) & (spike_times < times[-1])], units=spike_units,\n",
+    "        t_start=t_start, t_stop=times[-1]\n",
+    "    )\n",
+    "\n",
+    "    return cleaned_anas, sptr, best_channel"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def compute_spike_phase_func(lfp, times, return_degrees=False):\n",
+    "    from scipy.fftpack import next_fast_len\n",
+    "    x_a = scipy.signal.hilbert(\n",
+    "        lfp, next_fast_len(len(lfp)))[:len(lfp)]\n",
+    "#     x_a = hilbert(lfp)\n",
+    "    x_phase = np.angle(x_a, deg=return_degrees)\n",
+    "    return interp1d(times, x_phase)\n",
+    "\n",
+    "\n",
+    "def compute_spike_phase(lfp, spikes, flim=[6,10]):\n",
+    "    \n",
+    "    sample_rate = lfp.sampling_rate.magnitude\n",
+    "    \n",
+    "    # sometimes the position is recorded after LFP recording is ended\n",
+    "    times = np.arange(lfp.shape[0]) / sample_rate\n",
+    "     \n",
+    "    spikes = np.array(spikes)\n",
+    "    spikes = spikes[(spikes > times.min()) & (spikes < times.max())]\n",
+    "    \n",
+    "    filtered_lfp = butter_bandpass_filter(\n",
+    "        lfp.magnitude.ravel(), *flim, fs=sample_rate, order=3)\n",
+    "\n",
+    "    spike_phase_func = compute_spike_phase_func(filtered_lfp, times)\n",
+    "    \n",
+    "    return spike_phase_func(spikes), filtered_lfp"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plt.figure(figsize=(16,9))\n",
+    "# # lfp = data_loader.lfp('1833-200619-2', 6)\n",
+    "# lfp = data_loader.lfp('1834-220319-3', 6)\n",
+    "# # lfp = data_loader.lfp('1849-010319-4', 6)\n",
+    "# times = np.arange(lfp.shape[0]) / lfp.sampling_rate.magnitude\n",
+    "# clean_lfp, _ = compute_clean_lfp(lfp.magnitude[:, 0], threshold=2)\n",
+    "# plt.plot(times,lfp[:,0])\n",
+    "# plt.plot(times,clean_lfp)\n",
+    "\n",
+    "# plt.figure(figsize=(16,9))\n",
+    "# plt.psd(lfp[:,0].ravel(), Fs=1000, NFFT=10000)\n",
+    "# plt.psd(clean_lfp, Fs=1000, NFFT=10000)\n",
+    "# plt.xlim(0,100)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plt.figure(figsize=(16,9))\n",
+    "\n",
+    "# plt.plot(times,lfp[:,0])\n",
+    "# # plt.plot(clean_lfp*100)\n",
+    "# plt.plot(times[:-1], np.diff(lfp[:,0].magnitude.ravel()))\n",
+    "# plt.xlim(64.5,65.5)\n",
+    "# # plt.ylim(-250,250)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# # action_id_0, channel_0, unit_0 = '1833-200619-1', 6, 163\n",
+    "# # action_id_1, channel_1, unit_1 = '1833-200619-2', 6, 28\n",
+    "# action_id_0, channel_0, unit_0 = '1834-220319-3', 2, 46\n",
+    "# action_id_1, channel_1, unit_1 = '1834-220319-4', 2, 60\n",
+    "# lfp_0 = data_loader.lfp(action_id_0, channel_0)\n",
+    "# lfp_1 = data_loader.lfp(action_id_1, channel_1)\n",
+    "\n",
+    "# sample_rate_0 = lfp_0.sampling_rate\n",
+    "# sample_rate_1 = lfp_1.sampling_rate\n",
+    "\n",
+    "# lim_0 = get_lim(action_id_0)\n",
+    "# lim_1 = get_lim(action_id_1)\n",
+    "\n",
+    "# sptrs_0 = data_loader.spike_trains(action_id_0, channel_0)\n",
+    "\n",
+    "# sptrs_1 = data_loader.spike_trains(action_id_1, channel_1)\n",
+    "\n",
+    "# cleaned_lfps_0, sptr_0, best_channel_0 = prepare_spike_lfp(lfp_0, sptrs_0[unit_0], *lim_0)\n",
+    "\n",
+    "# cleaned_lfps_1, sptr_1, best_channel_1 = prepare_spike_lfp(lfp_1, sptrs_1[unit_1], *lim_1)\n",
+    "\n",
+    "# coher_0, freq_0 = compute_spike_lfp_coherence(cleaned_lfps_0, sptr_0, 4096)\n",
+    "\n",
+    "# coher_1, freq_1 = compute_spike_lfp_coherence(cleaned_lfps_1, sptr_1, 4096)\n",
+    "\n",
+    "# spike_phase_0, filtered_lfp_0 = compute_spike_phase(cleaned_lfps_0, sptrs_0[unit_0], flim=[6,10])\n",
+    "\n",
+    "# spike_phase_1, filtered_lfp_1 = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[6,10])\n",
+    "\n",
+    "# spike_phase_1_stim, filtered_lfp_1_stim = compute_spike_phase(cleaned_lfps_1, sptrs_1[unit_1], flim=[10.5,11.5])\n",
+    "\n",
+    "# plt.figure()\n",
+    "# plt.plot(freq_0, coher_0.ravel())\n",
+    "# plt.plot(freq_1, coher_1.ravel())\n",
+    "# plt.xlim(0,20)\n",
+    "\n",
+    "# plt.figure()\n",
+    "# bins_0, kde_0 = vonmises_kde(spike_phase_0, 100)\n",
+    "# ang_0, vec_len_0 = spike_phase_score(bins_0, kde_0)\n",
+    "# plt.polar(bins_0, kde_0, color='b')\n",
+    "# plt.polar([ang_0, ang_0], [0, vec_len_0], color='b')\n",
+    "\n",
+    "# bins_1, kde_1 = vonmises_kde(spike_phase_1, 100)\n",
+    "# ang_1, vec_len_1 = spike_phase_score(bins_1, kde_1)\n",
+    "# plt.polar(bins_1, kde_1, color='r')\n",
+    "# plt.polar([ang_1, ang_1], [0, vec_len_1], color='r')\n",
+    "\n",
+    "# bins_1_stim, kde_1_stim = vonmises_kde(spike_phase_1_stim, 100)\n",
+    "# ang_1_stim, vec_len_1_stim = spike_phase_score(bins_1_stim, kde_1_stim)\n",
+    "# plt.polar(bins_1_stim, kde_1_stim, color='k')\n",
+    "# plt.polar([ang_1_stim, ang_1_stim], [0, vec_len_1_stim], color='k')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# TODO fix artefact stuff from phase precession"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "NFFT = 8192\n",
+    "theta_band_f1, theta_band_f2 = 6, 10 "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coherence_data, freqency_data = {}, {}\n",
+    "theta_kde_data, theta_bins_data = {}, {}\n",
+    "stim_kde_data, stim_bins_data = {}, {}\n",
+    "\n",
     "def process(row):\n",
     "    action_id = row['action']\n",
     "    channel_group = row['channel_group']\n",
     "    unit_name = row['unit_name']\n",
     "    \n",
+    "    name = f'{action_id}_{channel_group}_{unit_name}'\n",
+    "    \n",
     "    lfp = data_loader.lfp(action_id, channel_group) # TODO consider choosing strongest stim response\n",
     "    \n",
     "    sptr = data_loader.spike_train(action_id, channel_group, unit_name)\n",
     "       \n",
     "    lim = get_lim(action_id)\n",
-    "\n",
-    "    p_xys, freq, clean_lfp = compute_spike_lfp(lfp, sptr, *lim, NFFT=NFFT)\n",
     "    \n",
-    "    snls = signaltonoise(clean_lfp)\n",
+    "    cleaned_lfp, sptr, best_channel = prepare_spike_lfp(lfp, sptr, *lim)\n",
-    "    best_channel = np.argmax(snls)\n",
+    "       \n",
-    "    snl = snls[best_channel]\n",
+    "    p_xys, freq = compute_spike_lfp_coherence(cleaned_lfp, sptr, NFFT=NFFT)\n",
-    "    p_xy = p_xys[:,best_channel].magnitude\n",
+    "    \n",
+    "    p_xy = p_xys.magnitude.ravel()\n",
     "    freq = freq.magnitude\n",
     "    \n",
     "    theta_f, theta_p_max = find_theta_peak(p_xy, freq, theta_band_f1, theta_band_f2)\n",
@@ -305,6 +519,12 @@
     "    \n",
     "    theta_half_energy = compute_energy(p_xy, freq, theta_half_f1, theta_half_f2) # theta band 6 - 10 Hz\n",
     "    \n",
+    "    theta_spike_phase, _ = compute_spike_phase(cleaned_lfp, sptr, flim=[theta_band_f1, theta_band_f2])\n",
+    "    theta_bins, theta_kde = vonmises_kde(theta_spike_phase)\n",
+    "    theta_ang, theta_vec_len = spike_phase_score(theta_bins, theta_kde)\n",
+    "    theta_kde_data.update({name: theta_kde})\n",
+    "    theta_bins_data.update({name: theta_bins})\n",
+    "\n",
     "    # stim\n",
     "    \n",
     "    stim_freq = compute_stim_freq(action_id)\n",
@@ -316,8 +536,19 @@
     "    \n",
     "    stim_energy = compute_energy(p_xy, freq, stim_half_f1, stim_half_f2)\n",
     "    \n",
+    "    if np.isnan(stim_freq):\n",
+    "        stim_spike_phase, stim_bins, stim_kde, stim_ang, stim_vec_len = [np.nan] * 5\n",
+    "    else:\n",
+    "        stim_spike_phase, _ = compute_spike_phase(cleaned_lfp, sptr, flim=[stim_freq - .5, stim_freq + .5])\n",
+    "        stim_bins, stim_kde = vonmises_kde(stim_spike_phase)\n",
+    "        stim_ang, stim_vec_len = spike_phase_score(stim_bins, stim_kde)\n",
+    "        stim_kde_data.update({name: stim_kde})\n",
+    "        stim_bins_data.update({name: stim_bins})\n",
+    "        \n",
+    "    coherence_data.update({name: p_xy})\n",
+    "    freqency_data.update({name: freq})\n",
+    "    \n",
     "    result = pd.Series({\n",
-    "        'signal_to_noise': snl,\n",
     "        'best_channel': best_channel,\n",
     "        'theta_freq': theta_f,\n",
     "        'theta_peak': theta_p_max,\n",
@@ -326,30 +557,34 @@
     "        'theta_half_f2': theta_half_f2,\n",
     "        'theta_half_width': theta_half_width,\n",
     "        'theta_half_energy': theta_half_energy,\n",
+    "        'theta_ang': theta_ang, \n",
+    "        'theta_vec_len': theta_vec_len,\n",
     "        'stim_freq': stim_freq,\n",
     "        'stim_p_max': stim_p_max,\n",
     "        'stim_half_f1': stim_half_f1, \n",
     "        'stim_half_f2': stim_half_f2,\n",
     "        'stim_half_width': stim_half_width,\n",
-    "        'stim_energy': stim_energy\n",
+    "        'stim_energy': stim_energy,\n",
+    "        'stim_ang': stim_ang, \n",
+    "        'stim_vec_len': stim_vec_len\n",
     "    })\n",
     "    return result"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 24,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "564d1ef5339f46eebee42ec41cfcfe62",
+       "model_id": "199c002aa9604451a318e3f471c9d892",
        "version_major": 2,
        "version_minor": 0
       },
       "text/plain": [
-       "HBox(children=(IntProgress(value=0, max=1298), HTML(value='')))"
+       "HBox(children=(IntProgress(value=0, max=1284), HTML(value='')))"
       ]
      },
      "metadata": {},
@@ -359,8 +594,9 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/signal/spectral.py:1578: RuntimeWarning: invalid value encountered in true_divide\n",
+      "/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/signal/spectral.py:1577: RuntimeWarning: invalid value encountered in true_divide\n",
-      "  Cxy = np.abs(Pxy)**2 / Pxx / Pyy\n"
+      "  Cxy = np.abs(Pxy)**2 / Pxx / Pyy\n",
+      "/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in true_divide\n"
      ]
     },
     {
@@ -377,87 +613,62 @@
     "    left_index=True, right_index=True)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# plot"
-   ]
-  },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 25,
    "metadata": {},
-   "outputs": [
+   "outputs": [],
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "b20a9a26e0864f578e1a9aa0c021999b",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "HBox(children=(IntProgress(value=0, max=1298), HTML(value='')))"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\n"
-     ]
-    }
-   ],
    "source": [
-    "coher, freqs = {}, {}\n",
+    "# coher, freqs = {}, {}\n",
-    "for i, row in tqdm(units.iterrows(), total=len(units)):\n",
+    "# for i, row in tqdm(units.iterrows(), total=len(units)):\n",
-    "    action_id = row['action']\n",
+    "#     action_id = row['action']\n",
-    "    channel_group = row['channel_group']\n",
+    "#     channel_group = row['channel_group']\n",
-    "    unit_name = row['unit_name']\n",
+    "#     unit_name = row['unit_name']\n",
     "    \n",
-    "    name = f'{action_id}_{channel_group}_{unit_name}'\n",
+    "#     name = f'{action_id}_{channel_group}_{unit_name}'\n",
     "    \n",
-    "    lfp = data_loader.lfp(action_id, channel_group) # TODO consider choosing strongest stim response\n",
+    "#     lfp = data_loader.lfp(action_id, channel_group) # TODO consider choosing strongest stim response\n",
     "    \n",
-    "    sptr = data_loader.spike_train(action_id, channel_group, unit_name)\n",
+    "#     sptr = data_loader.spike_train(action_id, channel_group, unit_name)\n",
     "       \n",
-    "    lim = get_lim(action_id)\n",
+    "#     lim = get_lim(action_id)\n",
     "\n",
-    "    p_xys, freq, clean_lfp = compute_spike_lfp(lfp, sptr, *lim, NFFT=NFFT)\n",
+    "#     p_xys, freq, clean_lfp = compute_spike_lfp(lfp, sptr, *lim, NFFT=NFFT)\n",
     "    \n",
-    "    snls = signaltonoise(clean_lfp)\n",
+    "#     snls = signaltonoise(clean_lfp)\n",
-    "    best_channel = np.argmax(snls)\n",
+    "#     best_channel = np.argmax(snls)\n",
-    "    snl = snls[best_channel]\n",
+    "#     snl = snls[best_channel]\n",
-    "    p_xy = p_xys[:,best_channel].magnitude\n",
+    "#     p_xy = p_xys[:,best_channel].magnitude\n",
-    "    freq = freq.magnitude\n",
+    "#     freq = freq.magnitude\n",
     "    \n",
-    "    coher.update({name: p_xy})\n",
+    "#     coher.update({name: p_xy})\n",
-    "    freqs.update({name: freq})"
+    "#     freqs.update({name: freq})"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 26,
    "metadata": {},
    "outputs": [
     {
-     "ename": "AttributeError",
+     "ename": "NameError",
-     "evalue": "'dict' object has no attribute 'to_feather'",
+     "evalue": "name 'frequency_data' is not defined",
      "output_type": "error",
      "traceback": [
       "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
-      "\u001b[0;32m<ipython-input-18-bf8a91a90c5a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mcoher\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_feather\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'data'\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'coherence.feather'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mfreqs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_feather\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'data'\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'freqs.feather'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+      "\u001b[0;32m<ipython-input-26-64114f1ea8bd>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoherence_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_feather\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'data'\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'coherence.feather'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfrequency_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_feather\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'data'\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'freqs.feather'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta_kde_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_feather\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'data'\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'theta_kde.feather'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta_bins_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_feather\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'data'\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'theta_bins.feather'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDataFrame\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstim_kde_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_feather\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'data'\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m'stim_kde.feather'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
-      "\u001b[0;31mAttributeError\u001b[0m: 'dict' object has no attribute 'to_feather'"
+      "\u001b[0;31mNameError\u001b[0m: name 'frequency_data' is not defined"
      ]
     }
    ],
    "source": [
-    "coher.to_feather(output / 'data' / 'coherence.feather')\n",
+    "pd.DataFrame(coherence_data).to_feather(output / 'data' / 'coherence.feather')\n",
-    "freqs.to_feather(output / 'data' / 'freqs.feather')"
+    "pd.DataFrame(frequency_data).to_feather(output / 'data' / 'freqs.feather')\n",
+    "pd.DataFrame(theta_kde_data).to_feather(output / 'data' / 'theta_kde.feather')\n",
+    "pd.DataFrame(theta_bins_data).to_feather(output / 'data' / 'theta_bins.feather')\n",
+    "pd.DataFrame(stim_kde_data).to_feather(output / 'data' / 'stim_kde.feather')\n",
+    "pd.DataFrame(stim_bins_data).to_feather(output / 'data' / 'stim_bins.feather')"
    ]
   },
   {
@@ -545,5 +756,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }