# 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.
#
import numpy as np
import scipy.stats as sps
from .utilities import smooth
[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.0
/ std_factor
* np.percentile(
np.abs(xx_centered),
[100 * (1 - 2 * (1 - sps.norm.cdf(std_factor)))],
)
)
std_y = (
1.0
/ 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 > 0.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=bool)
for yi in yes_set:
b[yi] = True
if debug_plots:
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