# 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.interpolate as si
import scipy.ndimage.filters as filt
import scipy.stats as stats
ON_LUMINANCE = 255
OFF_LUMINANCE = 0
[docs]def chi_square_binary(events, LSN_template):
# note: can only be applied to binary events for trial responses
#
# *****INPUT*****
# events: 2D numpy bool with shape (num_trials,num_cells) for presence
# or absence of response on a given trial
# LSN_template: 3D numpy int8 with shape (num_trials,num_y_pixels,
# num_x_pixels) for luminance at each pixel location
#
# *****OUTPUT*****
# chi_square_grid_NLL: 3D numpy float with shape (num_cells,num_y_pixels,
# num_x_pixels) that gives the p value for the hypothesis that a
# receptive field is contained within a 7x7 pixel mask centered on a
# given pixel location as measured by a chi-square test for the responses
# among the pixels (both on and off) that fall within the mask.
num_trials = np.shape(events)[0]
num_cells = np.shape(events)[1]
num_y = np.shape(LSN_template)[1]
num_x = np.shape(LSN_template)[2]
# for each pixel location, get a mask that is centered on that location
# disc_masks has shape (num_y,num_x,num_y,num_x)
disc_masks = get_disc_masks(LSN_template)
# determine which trials each pixel is active (i.e not gray),
# broken up by ON and OFF pixels.
# trial_matrix has shape (num_y,num_x,2,num_trials)
trial_matrix = build_trial_matrix(LSN_template, num_trials)
# get the total number of trials each pixel is active (i.e. not gray)
# trials_per_pixel has shape (num_y,num_x,2)
trials_per_pixel = np.sum(trial_matrix, axis=3)
# get the sum of the number of events across all trials that each pixel is
# active events_per_pixel has shape (num_cells,num_y,num_x,2)
events_per_pixel = get_events_per_pixel(events, trial_matrix)
# smooth stimulus-triggered average spatially with a gaussian
for n in range(num_cells):
for on_off in range(2):
events_per_pixel[n, :, :, on_off] = smooth_STA(
events_per_pixel[n, :, :, on_off]
)
# calculate the p_value for each exclusion region
chi_square_grid = np.zeros((num_cells, num_y, num_x))
for y in range(num_y):
for x in range(num_x):
exclusion_mask = np.ones((num_y, num_x, 2)) * disc_masks[
y, x, :, :
].reshape(num_y, num_x, 1)
p_vals, __ = chi_square_within_mask(
exclusion_mask, events_per_pixel, trials_per_pixel
)
chi_square_grid[:, y, x] = p_vals
return chi_square_grid
[docs]def get_peak_significance(chi_square_grid_NLL, LSN_template, alpha=0.05):
# ****INPUT*****
# chi_square_grid_NLL: result of chi_square_binary(events,LSN_template)
# LSN_template: 3D numpy int8 with shape (num_trials,num_y_pixels,
# num_x_pixels) for luminance at each pixel location
#
# *****OUTPUT*****
# significant_cells: 1D numpy bool with shape (num_cells,) that indicates
# whether or not each cell has a location with a significant chance of a
# true receptive field
# best_exclusion_region_list: list of 2D numpy bool with len (num_cells).
# For each cell, the array is a mask for the best pixel, or all false
num_cells = np.shape(chi_square_grid_NLL)[0]
num_y = np.shape(chi_square_grid_NLL)[1]
num_x = np.shape(chi_square_grid_NLL)[2]
chi_square_grid = NLL_to_pvalue(chi_square_grid_NLL)
# get the average size of all masks in units of number of pixels
disc_masks = get_disc_masks(LSN_template)
pixels_per_mask_per_pixel = np.sum(disc_masks, axis=(2, 3)).astype(float)
# find the smallest p-value and determine if it's significant
significant_cells = np.zeros((num_cells)).astype(bool)
best_p = np.zeros((num_cells))
p_value_correction_factor_per_pixel = (
1.0 * num_y * num_x / pixels_per_mask_per_pixel
)
best_exclusion_region_list = []
corrected_p_value_array_list = []
for n in range(num_cells):
# Sidak correction:
p_value_corrected_per_pixel = 1 - np.power(
(1 - chi_square_grid[n, :, :]), p_value_correction_factor_per_pixel
)
corrected_p_value_array_list.append(p_value_corrected_per_pixel)
y, x = np.unravel_index(
p_value_corrected_per_pixel.argmin(), (num_y, num_x)
)
# if more than one p-value that maxes out, use the median location
if np.sum(p_value_corrected_per_pixel == 0.0) > 1:
y, x = np.unravel_index(
np.argwhere(p_value_corrected_per_pixel.flatten() == 0.0)[
:, 0
],
(num_y, num_x),
)
center_y, center_x = locate_median(y, x)
best_p[n] = p_value_corrected_per_pixel[y, x]
if best_p[n] < alpha:
significant_cells[n] = True
best_exclusion_region_list.append(
disc_masks[y, x, :, :].astype(bool)
)
else:
best_exclusion_region_list.append(
np.zeros(
(disc_masks.shape[0], disc_masks.shape[1]), dtype=bool
)
)
return (
significant_cells,
best_p,
corrected_p_value_array_list,
best_exclusion_region_list,
)
[docs]def pvalue_to_NLL(p_values, max_NLL=10.0):
return np.where(p_values == 0.0, max_NLL, -np.log10(p_values))
[docs]def NLL_to_pvalue(NLLs, log_base=10.0):
return log_base ** (-NLLs)
[docs]def get_events_per_pixel(responses_np, trial_matrix):
"""Obtain a matrix linking cellular responses to pixel activity.
Parameters
----------
responses_np : np.ndarray
Dimensions are (nTrials, nCells). Boolean values indicate
presence/absence of a response on a given trial.
trial_matrix : np.ndarray
Dimensions are (nYPixels, nXPixels, {on, off}, nTrials). Boolean values
indicate that a pixel was on/off on a particular trial.
Returns
-------
events_per_pixel : np.ndarray
Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Values for each
cell, pixel, and on/off state are the sum of events for that cell
across all trials where the pixel was in the on/off state.
"""
num_cells = np.shape(responses_np)[1]
num_y = np.shape(trial_matrix)[0]
num_x = np.shape(trial_matrix)[1]
events_per_pixel = np.zeros((num_cells, num_y, num_x, 2))
for y in range(num_y):
for x in range(num_x):
for on_off in range(2):
frames = np.argwhere(trial_matrix[y, x, on_off, :])[:, 0]
events_per_pixel[:, y, x, on_off] = np.sum(
responses_np[frames, :], axis=0
)
return events_per_pixel
[docs]def smooth_STA(STA, gauss_std=0.75, total_degrees=64):
"""Smooth an image by convolution with a gaussian kernel
Parameters
----------
STA : np.ndarray
Input image
gauss_std : numeric, optional
Standard deviation of the gaussian kernel. Will be applied to the
upsampled image, so units are visual degrees. Default is 0.75
total_degrees : int, optional
Size in visual degrees of the input image along its zeroth (row) axis.
Used to set the scale factor for up/downsampling.
Returns
-------
STA_smoothed : np.ndarray
Smoothed image
"""
deg_per_pnt = total_degrees // STA.shape[0]
STA_interpolated = interpolate_RF(STA, deg_per_pnt)
STA_interpolated_smoothed = filt.gaussian_filter(
STA_interpolated, gauss_std
)
STA_smoothed = deinterpolate_RF(
STA_interpolated_smoothed, STA.shape[1], STA.shape[0], deg_per_pnt
)
return STA_smoothed
[docs]def interpolate_RF(rf_map, deg_per_pnt):
"""Upsample an image
Parameters
----------
rf_map : np.ndarray
Input image
deg_per_pnt : numeric
scale factor
Returns
-------
interpolated : np.ndarray
Upsampled image
"""
x_pnts = np.shape(rf_map)[1]
y_pnts = np.shape(rf_map)[0]
x_coor = np.arange(
-(x_pnts - 1) * deg_per_pnt / 2,
(x_pnts + 1) * deg_per_pnt / 2,
deg_per_pnt,
)
y_coor = np.arange(
-(y_pnts - 1) * deg_per_pnt / 2,
(y_pnts + 1) * deg_per_pnt / 2,
deg_per_pnt,
)
x_interpolated = np.arange(
-(x_pnts - 1) * deg_per_pnt / 2,
deg_per_pnt / 2 + (x_pnts / 2 - 1) * deg_per_pnt + 1,
1,
)
y_interpolated = np.arange(
-(y_pnts - 1) * deg_per_pnt / 2,
deg_per_pnt / 2 + (y_pnts / 2 - 1) * deg_per_pnt + 1,
1,
)
interpolated = si.interp2d(x_coor, y_coor, rf_map)
interpolated = interpolated(x_interpolated, y_interpolated)
return interpolated
[docs]def deinterpolate_RF(rf_map, x_pnts, y_pnts, deg_per_pnt):
"""Downsample an image
Parameters
----------
rf_map : np.ndarray
Input image
x_pnts : np.ndarray
Count of sample points along the first (column) axis
y_pnts : np.ndarray
Count of sample points along the zeroth (row) axis
deg_per_pnt : numeric
scale factor
Returns
-------
sampled_yx : np.ndarray
Downsampled image
"""
# x_pnts = 28
# y_pnts = 16
x_interpolated = np.arange(
-(x_pnts - 1) * deg_per_pnt / 2,
deg_per_pnt / 2 + (x_pnts / 2 - 1) * deg_per_pnt + 1,
1,
)
y_interpolated = np.arange(
-(y_pnts - 1) * deg_per_pnt / 2,
deg_per_pnt / 2 + (y_pnts / 2 - 1) * deg_per_pnt + 1,
1,
)
x_deinterpolate = np.arange(0, len(x_interpolated), deg_per_pnt)
y_deinterpolate = np.arange(0, len(y_interpolated), deg_per_pnt)
sampled_y = rf_map[y_deinterpolate, :]
sampled_yx = sampled_y[:, x_deinterpolate]
return sampled_yx
[docs]def chi_square_within_mask(exclusion_mask, events_per_pixel, trials_per_pixel):
"""Determine if cells respond preferentially to on/off pixels in a mask
using a chi2 test.
Parameters
----------
exclusion_mask : np.ndarray
Dimensions are (nYPixels, nXPixels, {on, off}). Integer indicator for
INCLUSION (!) of a pixel within the testing region.
events_per_pixel : np.ndarray
Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Integer values
are response counts by cell to on/off luminance at each pixel.
trials_per_pixel : np.ndarray
Dimensions are (nYPixels, nXPixels, {on, off}). Integer values are
counts of trials where a pixel is on/off.
Returns
-------
p_vals : np.ndarray
One-dimensional, of length nCells. Float values are p-values
for the hypothesis that a given cell has a receptive field within the
exclusion mask.
chi : np.ndarray
Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Values (float)
are squared residual event counts divided by expected event counts.
"""
num_y = np.shape(exclusion_mask)[0]
num_x = np.shape(exclusion_mask)[1]
# d.f. is number of pixels in mask minus one
degrees_of_freedom = int(np.sum(exclusion_mask)) - 1
# observed_by_pixel has shape (num_cells,num_y,num_x,2)
expected_by_pixel = get_expected_events_by_pixel(
exclusion_mask, events_per_pixel, trials_per_pixel
)
observed_by_pixel = (
events_per_pixel * exclusion_mask.reshape(1, num_y, num_x, 2)
).astype(float)
# calculate test statistic given observed and expected
residual_by_pixel = observed_by_pixel - expected_by_pixel
chi = (residual_by_pixel**2) / expected_by_pixel
chi_sum = np.nansum(chi, axis=(1, 2, 3))
# get p-value given test statistic and degrees of freedom
p_vals = 1.0 - stats.chi2.cdf(chi_sum, degrees_of_freedom)
return p_vals, chi
[docs]def get_expected_events_by_pixel(
exclusion_mask, events_per_pixel, trials_per_pixel
):
"""Calculate expected number of events per pixel
Parameters
----------
exclusion_mask : np.ndarray
Dimensions are (nYPixels, nXPixels, {on, off}). Integer indicator for
INCLUSION (!) of a pixel within the testing region.
events_per_pixel : np.ndarray
Dimensions are (nCells, nYPixels, nXPixels, {on, off}). Integer values
are response counts by cell to on/off luminance at each pixel.
trials_per_pixel : np.ndarray
Dimensions are (nYPixels, nXPixels, {on, off}). Integer values are
counts of trials where a pixel is on/off.
Returns
-------
np.ndarray :
Dimensions (nCells, nYPixels, nXPixels, {on, off}). Float values are
pixelwise counts of events expected if events are evenly distributed
in mask across trials.
"""
num_y = np.shape(exclusion_mask)[0]
num_x = np.shape(exclusion_mask)[1]
num_cells = np.shape(events_per_pixel)[0]
exclusion_mask = exclusion_mask.reshape(1, num_y, num_x, 2)
trials_per_pixel = trials_per_pixel.reshape(1, num_y, num_x, 2)
masked_trials = exclusion_mask * trials_per_pixel
masked_events = exclusion_mask * events_per_pixel
total_trials = np.sum(masked_trials).astype(float)
total_events_by_cell = np.sum(masked_events, axis=(1, 2, 3)).astype(float)
expected_by_cell_per_trial = total_events_by_cell / total_trials
return masked_trials * expected_by_cell_per_trial.reshape(
num_cells, 1, 1, 1
)
[docs]def build_trial_matrix(
LSN_template, num_trials, on_off_luminance=(ON_LUMINANCE, OFF_LUMINANCE)
):
"""Construct indicator arrays for on/off pixels across trials.
Parameters
----------
LSN_template : np.ndarray
Dimensions are (nTrials, nYPixels, nXPixels). Luminance values per
pixel and trial. The size of the first dimension may be larger than the
num_trials argument (in which case only the first num_trials slices
will be used) but may not be smaller.
num_trials : int
The number of trials (left-justified) to build indicators for.
on_off_luminance : array-like, optional
The zeroth element is the luminance value of a pixel when on, the first
when off. Defaults are [255, 0].
Returns
-------
trial_mat : np.ndarray
Dimensions are (nYPixels, nXPixels, {on, off}, nTrials). Boolean values
indicate that a pixel was on/off on a particular trial.
"""
_, num_y, num_x = np.shape(LSN_template)
trial_mat = np.zeros((num_y, num_x, 2, num_trials), dtype=bool)
for y in range(num_y):
for x in range(num_x):
for oo, on_off in enumerate(on_off_luminance):
frame = np.argwhere(LSN_template[:num_trials, y, x] == on_off)[
:, 0
]
trial_mat[y, x, oo, frame] = True
return trial_mat
[docs]def get_disc_masks(
LSN_template,
radius=3,
on_luminance=ON_LUMINANCE,
off_luminance=OFF_LUMINANCE,
):
"""Obtain an indicator mask surrounding each pixel. The mask is a square,
excluding pixels which are coactive on any trial with the main pixel.
Parameters
----------
LSN_template : np.ndarray
Dimensions are (nTrials, nYPixels, nXPixels). Luminance values per
pixel and trial.
radius : int
The base mask will be a box whose sides are 2 * radius + 1 in length.
on_luminance : int, optional
The value of the luminance for on trials. Default is 255
off_luminance : int, optional
The value of the luminance for off trials. Default is 0
Returns
-------
masks : np.ndarray
Dimensions are (nYPixels, nXPixels, nYPixels, nXPixels). The first 2
dimensions describe the pixel from which the mask was computed. The
last 2 serve as the dimensions of the mask images themselves. Masks are
binary arrays of type float, with 1 indicating inside, 0 outside.
"""
num_y = np.shape(LSN_template)[1]
num_x = np.shape(LSN_template)[2]
# convert template to true on trials a pixel is not gray and false when
# gray
LSN_binary = np.where(LSN_template == off_luminance, 1, LSN_template)
LSN_binary = np.where(LSN_binary == on_luminance, 1, LSN_binary)
LSN_binary = np.where(LSN_binary == 1, 1.0, 0.0)
# get number of trials each pixel is not gray
on_trials = LSN_binary.sum(axis=0).astype(float) # shape is (num_y,num_x)
masks = np.zeros((num_y, num_x, num_y, num_x))
for y in range(num_y):
for x in range(num_x):
trials_not_gray = np.argwhere(LSN_binary[:, y, x] > 0)[:, 0]
raw_mask = np.divide(
LSN_binary[trials_not_gray, :, :].sum(axis=0), on_trials
)
center_y, center_x = np.unravel_index(
raw_mask.argmax(), (num_y, num_x)
)
# include center pixel in mask
raw_mask[center_y, center_x] = 0.0
x_max = center_x + radius + 1
if x_max > num_x:
x_max = num_x
x_min = center_x - radius
if x_min < 0:
x_min = 0
y_max = center_y + radius + 1
if y_max > num_y:
y_max = num_y
y_min = center_y - radius
if y_min < 0:
y_min = 0
# don't include far away pixels that just happen
# to not have any trials in common with center pixel
clean_mask = np.ones(np.shape(raw_mask))
clean_mask[y_min:y_max, x_min:x_max] = raw_mask[
y_min:y_max, x_min:x_max
]
masks[y, x, :, :] = clean_mask
masks = np.where(masks > 0, 0.0, 1.0)
return masks