# Allen Institute Software License - This software license is the 2-clause BSD
# license plus a third clause that prohibits redistribution for commercial
# purposes without further permission.
#
# Copyright 2017. Allen Institute. All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
# 3. Redistributions for commercial purposes are not permitted without the
# Allen Institute's written permission.
# For purposes of this license, commercial purposes is the incorporation of the
# Allen Institute's software into anything for which you will charge fees or
# other compensation. Contact terms@alleninstitute.org for commercial licensing
# opportunities.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#
from .utilities import smooth
import numpy as np
import scipy.stats as sps
[docs]def detect_events(data, cell_index, stimulus, debug_plots=False):
stimulus_table = data.get_stimulus_table(stimulus)
dff_trace = data.get_dff_traces()[1][cell_index, :]
k_min = 0
k_max = 10
delta = 3
dff_trace = smooth(dff_trace, 5)
var_dict = {}
debug_dict = {}
for ii, fi in enumerate(stimulus_table['start'].values):
if ii > 0 and stimulus_table.iloc[ii].start == stimulus_table.iloc[ii-1].end:
offset = 1
else:
offset = 0
if fi + k_min >= 0 and fi + k_max <= len(dff_trace):
trace = dff_trace[fi + k_min+1+offset:fi + k_max+1+offset]
xx = (trace - trace[0])[delta] - (trace - trace[0])[0]
yy = max((trace - trace[0])[delta + 2] - (trace - trace[0])[0 + 2],
(trace - trace[0])[delta + 3] - (trace - trace[0])[0 + 3],
(trace - trace[0])[delta + 4] - (trace - trace[0])[0 + 4])
var_dict[ii] = (trace[0], trace[-1], xx, yy)
debug_dict[fi + k_min+1+offset] = (ii, trace)
xx_list, yy_list = [], []
for _, _, xx, yy in var_dict.values():
xx_list.append(xx)
yy_list.append(yy)
mu_x = np.median(xx_list)
mu_y = np.median(yy_list)
xx_centered = np.array(xx_list)-mu_x
yy_centered = np.array(yy_list)-mu_y
std_factor = 1
std_x = 1./std_factor*np.percentile(np.abs(xx_centered), [100*(1-2*(1-sps.norm.cdf(std_factor)))])
std_y = 1./std_factor*np.percentile(np.abs(yy_centered), [100*(1-2*(1-sps.norm.cdf(std_factor)))])
curr_inds = []
allowed_sigma = 4
for ii, (xi, yi) in enumerate(zip(xx_centered, yy_centered)):
if np.sqrt(((xi)/std_x)**2+((yi)/std_y)**2) < allowed_sigma:
curr_inds.append(True)
else:
curr_inds.append(False)
curr_inds = np.array(curr_inds)
data_x = xx_centered[curr_inds]
data_y = yy_centered[curr_inds]
Cov = np.cov(data_x, data_y)
Cov_Factor = np.linalg.cholesky(Cov)
Cov_Factor_Inv = np.linalg.inv(Cov_Factor)
#===================================================================================================================
noise_threshold = max(allowed_sigma * std_x + mu_x, allowed_sigma * std_y + mu_y)
mu_array = np.array([mu_x, mu_y])
yes_set, no_set = set(), set()
for ii, (t0, tf, xx, yy) in var_dict.items():
xi_z, yi_z = Cov_Factor_Inv.dot((np.array([xx,yy]) - mu_array))
# Conditions in order:
# 1) Outside noise blob
# 2) Minimum change in df/f
# 3) Change evoked by this trial, not previous
# 4) At end of trace, ended up outside of noise floor
if np.sqrt(xi_z**2 + yi_z**2) > 4 and yy > .05 and xx < yy and tf > noise_threshold/2:
yes_set.add(ii)
else:
no_set.add(ii)
assert len(var_dict) == len(stimulus_table)
b = np.zeros(len(stimulus_table), dtype=np.bool)
for yi in yes_set:
b[yi] = True
if debug_plots == True:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,2)
# ax[0].plot(dff_trace)
for key, val in debug_dict.items():
ti, trace = val
if ti in no_set:
ax[0].plot(np.arange(key, key+len(trace)), trace, 'b')
elif ti in yes_set:
ax[0].plot(np.arange(key, key + len(trace)), trace, 'r', linewidth=2)
else:
raise Exception
for ii in yes_set:
ax[1].plot([var_dict[ii][2]], [var_dict[ii][3]], 'r.')
for ii in no_set:
ax[1].plot([var_dict[ii][2]], [var_dict[ii][3]], 'b.')
print('number_of_events: %d' % b.sum())
plt.show()
return b