Source code for allensdk.brain_observatory.receptive_field_analysis.eventdetection

# 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