From 60f1c7559ee6f80af795cb32132b068d190b2e55 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Wed, 28 Oct 2020 18:11:41 -0400 Subject: [PATCH 01/39] Initial code draft --- nwbwidgets/__init__.py | 1 + nwbwidgets/placefield.py | 291 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 292 insertions(+) create mode 100644 nwbwidgets/placefield.py diff --git a/nwbwidgets/__init__.py b/nwbwidgets/__init__.py index 7eb06b8f..8372e02c 100644 --- a/nwbwidgets/__init__.py +++ b/nwbwidgets/__init__.py @@ -1,6 +1,7 @@ import plotly.io as pio from .view import nwb2widget, default_neurodata_vis_spec +import .placefield # from .ephys_viz_interface import ephys_viz_neurodata_vis_spec diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py new file mode 100644 index 00000000..c27fe8c6 --- /dev/null +++ b/nwbwidgets/placefield.py @@ -0,0 +1,291 @@ +import numpy as np +from scipy.ndimage import label +from scipy.ndimage.filters import gaussian_filter, maximum_filter + +import pynwb +from pynwb.misc import Units +from .behavior import plotly_show_spatial_trace +from .base import vis2widget +from ipywidgets import widgets +from .utils.units import get_spike_times + +## To-do +# [] Create PlaceFieldWidget class + # [X] Refactor place field calculation code to deal with nwb data type + # [X] Incorporate place field fxns into class + # [X] Change all internal attributes references + # [X]Change all internal method references + + # [X] Get pos + # [X] Get time + # [X] Get spikes + # [] Get trials / epochs + + # [] Submit draft PR + + # [] Modify plotly_show_spatial_trace to plot 2D heatmap representing place fields or create new figure function? + # [] Dropdown that controls which unit + + # [] Work in buttons / dropdowns / sliders to modify following parameters in place field calculation: + # [] Different epochs + # [] Gaussian SD + # [] Speed threshold + # [] Minimum firing rate + # [] Place field thresh (% of local max) + +# Put widget rendering here +class PlaceFieldWidget(widgets.HBox): + + def __init__(self, unit: Units, node: pynwb.behavior.SpatialSeries): + + super().__init__() + + # Initialize receptive fields + self.receptive_fields = np.zeros(firing_rate.shape, dtype=int) + # Get pos + pos, unit = get_timeseries_in_units(node) + # Get time + pos_tt = get_timeseries_tt(node) + # Get spikes + spikes = get_spike_times(unit) + # Pixel width? + pixel_width = + + # Put widget controls here: + # - Gaussian SD + # - Speed threshold + # - Minimum firing rate + # - Place field thresh (% of local max) + + self.compute_2d_firing_rate(pos, pos_tt, spikes, pixel_width) # speed_thresh=0.03, gaussian_sd=0.0184, + # x_start=None, x_stop=None, y_start=None, y_stop=None) + + self.compute_2d_place_fields() # min_firing_rate=1, thresh=0.2, + # min_size=100): + + + return vis2widget(plotly_show_spatial_trace(self.receptive_fields)) + + + # Put place field code here + def find_nearest(self, arr, tt): + """Used for picking out elements of a TimeSeries based on spike times + + Parameters + ---------- + arr + tt + + Returns + ------- + + """ + arr = arr[arr > tt[0]] + arr = arr[arr < tt[-1]] + return np.searchsorted(tt, arr) + + def smooth(self, y, box_pts): + """Moving average + + Parameters + ---------- + y + box_pts + + Returns + ------- + + """ + box = np.ones(box_pts) / box_pts + return np.convolve(y, box, mode='same') + + def compute_speed(self, pos, pos_tt, smooth_param=40): + """Compute boolean of whether the speed of the animal was above a threshold + for each time point + + Parameters + ---------- + pos: np.ndarray(dtype=float) + in meters + pos_tt: np.ndarray(dtype=float) + in seconds + smooth_param: float, optional + + Returns + ------- + running: np.ndarray(dtype=bool) + + """ + speed = np.hstack((0, np.sqrt(np.sum(np.diff(pos.T) ** 2, axis=0)) / np.diff(pos_tt))) + return self.smooth(speed, smooth_param) + + + def compute_2d_occupancy(self, pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): + """Computes occupancy per bin in seconds + + Parameters + ---------- + pos: np.ndarray(dtype=float) + in meters + pos_tt: np.ndarray(dtype=float) + in seconds + edges_x: array-like + edges of histogram in meters + edges_y: array-like + edges of histogram in meters + speed_thresh: float, optional + in meters. Default = 3.0 cm/s + + Returns + ------- + occupancy: np.ndarray(dtype=float) + in seconds + running: np.ndarray(dtype=bool) + + """ + + sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) + is_running = self.compute_speed(pos, pos_tt) > speed_thresh + run_pos = pos[is_running, :] + occupancy = np.histogram2d(run_pos[:, 0], + run_pos[:, 1], + [edges_x, edges_y])[0] * sampling_period # in seconds + + return occupancy, is_running + + + def compute_2d_n_spikes(self, pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03): + """Returns speed-gated position during spikes + + Parameters + ---------- + pos: np.ndarray(dtype=float) + (time x 2) in meters + pos_tt: np.ndarray(dtype=float) + (time,) in seconds + spikes: np.ndarray(dtype=float) + in seconds + edges_x: np.ndarray(dtype=float) + edges of histogram in meters + edges_y: np.ndarray(dtype=float) + edges of histogram in meters + speed_thresh: float + in meters. Default = 3.0 cm/s + + Returns + ------- + """ + + is_running = self.compute_speed(pos, pos_tt) > speed_thresh + + spike_pos_inds = self.find_nearest(spikes, pos_tt) + spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] + pos_on_spikes = pos[spike_pos_inds, :] + + n_spikes = np.histogram2d(pos_on_spikes[:, 0], + pos_on_spikes[:, 1], + [edges_x, edges_y])[0] + + return n_spikes + + + def compute_2d_firing_rate(self, pos, pos_tt, spikes, + pixel_width, + speed_thresh=0.03, + gaussian_sd=0.0184, + x_start=None, x_stop=None, + y_start=None, y_stop=None): + """Returns speed-gated occupancy and speed-gated and + Gaussian-filtered firing rate + + Parameters + ---------- + pos: np.ndarray(dtype=float) + (time x 2), in meters + pos_tt: np.ndarray(dtype=float) + (time,) in seconds + spikes: np.ndarray(dtype=float) + in seconds + pixel_width: float + speed_thresh: float, optional + in meters. Default = 3.0 cm/s + gaussian_sd: float, optional + in meters. Default = 1.84 cm + x_start: float, optional + x_stop: float, optional + y_start: float, optional + y_stop: float, optional + + + Returns + ------- + + occupancy: np.ndarray + in seconds + filtered_firing_rate: np.ndarray(shape=(cell, x, y), dtype=float) + in Hz + + """ + # pixel_width=0.0092, + # field_len = 0.46 + # edges = np.arange(0, field_len + pixel_width, pixel_width) + + x_start = np.min(pos[:, 0]) if x_start is None else x_start + x_stop = np.max(pos[:, 0]) if x_stop is None else x_stop + + y_start = np.min(pos[:, 1]) if y_start is None else y_start + y_stop = np.max(pos[:, 1]) if y_stop is None else y_stop + + self.edges_x = np.arange(x_start, x_stop, pixel_width) + self.edges_y = np.arange(y_start, y_stop, pixel_width) + + self.occupancy, running = self.compute_2d_occupancy(pos, pos_tt, self.edges_x, self.edges_y, speed_thresh) + + n_spikes = self.compute_2d_n_spikes(pos, pos_tt, spikes, self.edges_x, self.edges_y, speed_thresh) + + firing_rate = n_spikes / self.occupancy # in Hz + firing_rate[np.isnan(firing_rate)] = 0 + + self.firing_rate = gaussian_filter(firing_rate, gaussian_sd / pixel_width) + + # return occupancy, filtered_firing_rate, [edges_x, edges_y] + + + def compute_2d_place_fields(self, min_firing_rate=1, thresh=0.2, + min_size=100): + """Compute place fields + + Parameters + ---------- + firing_rate: np.ndarray(NxN, dtype=float) + min_firing_rate: float + in Hz + thresh: float + % of local max + min_size: float + minimum size of place field in pixels + + Returns + ------- + receptive_fields: np.ndarray(NxN, dtype=int) + Each receptive field is labeled with a unique integer + """ + firing_rate = self.firing_rate + local_maxima_inds = firing_rate == maximum_filter(firing_rate, 3) + n_receptive_fields = 0 + firing_rate = firing_rate.copy() + for local_max in np.flipud(np.sort(firing_rate[local_maxima_inds])): + labeled_image, num_labels = label(firing_rate > max(local_max * thresh, + min_firing_rate)) + if not num_labels: # nothing above min_firing_thresh + return + for i in range(1, num_labels + 1): + image_label = labeled_image == i + if local_max in firing_rate[image_label]: + break + if np.sum(image_label) >= min_size: + n_receptive_fields += 1 + self.receptive_fields[image_label] = n_receptive_fields + firing_rate[image_label] = 0 + + # return receptive_fields \ No newline at end of file From 9e9c7c3448b1e6c79fae5c4dd937d260780a592b Mon Sep 17 00:00:00 2001 From: "!git for-each-ref --format='%(refname:short)' `git symbolic-ref HEAD`" Date: Wed, 28 Oct 2020 19:20:09 -0400 Subject: [PATCH 02/39] place fields almost working --- nwbwidgets/__init__.py | 2 +- nwbwidgets/placefield.py | 530 ++++++++++++++++++++------------------- nwbwidgets/view.py | 5 +- 3 files changed, 281 insertions(+), 256 deletions(-) diff --git a/nwbwidgets/__init__.py b/nwbwidgets/__init__.py index 8372e02c..3e01a335 100644 --- a/nwbwidgets/__init__.py +++ b/nwbwidgets/__init__.py @@ -1,7 +1,7 @@ import plotly.io as pio from .view import nwb2widget, default_neurodata_vis_spec -import .placefield +from .placefield import * # from .ephys_viz_interface import ephys_viz_neurodata_vis_spec diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index c27fe8c6..fa305ea9 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -2,290 +2,314 @@ from scipy.ndimage import label from scipy.ndimage.filters import gaussian_filter, maximum_filter +import matplotlib.pyplot as plt + import pynwb from pynwb.misc import Units from .behavior import plotly_show_spatial_trace from .base import vis2widget from ipywidgets import widgets from .utils.units import get_spike_times +from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt + +import plotly.graph_objects as go + + ## To-do # [] Create PlaceFieldWidget class - # [X] Refactor place field calculation code to deal with nwb data type - # [X] Incorporate place field fxns into class - # [X] Change all internal attributes references - # [X]Change all internal method references +# [X] Refactor place field calculation code to deal with nwb data type +# [X] Incorporate place field fxns into class +# [X] Change all internal attributes references +# [X]Change all internal method references - # [X] Get pos - # [X] Get time - # [X] Get spikes - # [] Get trials / epochs +# [X] Get pos +# [X] Get time +# [X] Get spikes +# [] Get trials / epochs - # [] Submit draft PR +# [] Submit draft PR - # [] Modify plotly_show_spatial_trace to plot 2D heatmap representing place fields or create new figure function? - # [] Dropdown that controls which unit +# [] Modify plotly_show_spatial_trace to plot 2D heatmap representing place fields or create new figure function? +# [] Dropdown that controls which unit - # [] Work in buttons / dropdowns / sliders to modify following parameters in place field calculation: - # [] Different epochs - # [] Gaussian SD - # [] Speed threshold - # [] Minimum firing rate - # [] Place field thresh (% of local max) +# [] Work in buttons / dropdowns / sliders to modify following parameters in place field calculation: +# [] Different epochs +# [] Gaussian SD +# [] Speed threshold +# [] Minimum firing rate +# [] Place field thresh (% of local max) # Put widget rendering here +def find_nearest(arr, tt): + """Used for picking out elements of a TimeSeries based on spike times + + Parameters + ---------- + arr + tt + + Returns + ------- + + """ + arr = arr[arr > tt[0]] + arr = arr[arr < tt[-1]] + return np.searchsorted(tt, arr) + + +def smooth(y, box_pts): + """Moving average + + Parameters + ---------- + y + box_pts + + Returns + ------- + + """ + box = np.ones(box_pts) / box_pts + return np.convolve(y, box, mode='same') + + +def compute_speed(pos, pos_tt, smooth_param=40): + """Compute boolean of whether the speed of the animal was above a threshold + for each time point + + Parameters + ---------- + pos: np.ndarray(dtype=float) + in meters + pos_tt: np.ndarray(dtype=float) + in seconds + smooth_param: float, optional + + Returns + ------- + running: np.ndarray(dtype=bool) + + """ + speed = np.hstack((0, np.sqrt(np.sum(np.diff(pos.T) ** 2, axis=0)) / np.diff(pos_tt))) + return smooth(speed, smooth_param) + + +def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): + """Computes occupancy per bin in seconds + + Parameters + ---------- + pos: np.ndarray(dtype=float) + in meters + pos_tt: np.ndarray(dtype=float) + in seconds + edges_x: array-like + edges of histogram in meters + edges_y: array-like + edges of histogram in meters + speed_thresh: float, optional + in meters. Default = 3.0 cm/s + + Returns + ------- + occupancy: np.ndarray(dtype=float) + in seconds + running: np.ndarray(dtype=bool) + + """ + + sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) + is_running = compute_speed(pos, pos_tt) > speed_thresh + run_pos = pos[is_running, :] + occupancy = np.histogram2d(run_pos[:, 0], + run_pos[:, 1], + [edges_x, edges_y])[0] * sampling_period # in seconds + + return occupancy, is_running + + +def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03): + """Returns speed-gated position during spikes + + Parameters + ---------- + pos: np.ndarray(dtype=float) + (time x 2) in meters + pos_tt: np.ndarray(dtype=float) + (time,) in seconds + spikes: np.ndarray(dtype=float) + in seconds + edges_x: np.ndarray(dtype=float) + edges of histogram in meters + edges_y: np.ndarray(dtype=float) + edges of histogram in meters + speed_thresh: float + in meters. Default = 3.0 cm/s + + Returns + ------- + """ + + is_running = compute_speed(pos, pos_tt) > speed_thresh + + spike_pos_inds = find_nearest(spikes, pos_tt) + spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] + pos_on_spikes = pos[spike_pos_inds, :] + + n_spikes = np.histogram2d(pos_on_spikes[:, 0], + pos_on_spikes[:, 1], + [edges_x, edges_y])[0] + + return n_spikes + + +def compute_2d_firing_rate(pos, pos_tt, spikes, + pixel_width, + speed_thresh=0.03, + gaussian_sd=0.0184, + x_start=None, x_stop=None, + y_start=None, y_stop=None): + """Returns speed-gated occupancy and speed-gated and + Gaussian-filtered firing rate + + Parameters + ---------- + pos: np.ndarray(dtype=float) + (time x 2), in meters + pos_tt: np.ndarray(dtype=float) + (time,) in seconds + spikes: np.ndarray(dtype=float) + in seconds + pixel_width: float + speed_thresh: float, optional + in meters. Default = 3.0 cm/s + gaussian_sd: float, optional + in meters. Default = 1.84 cm + x_start: float, optional + x_stop: float, optional + y_start: float, optional + y_stop: float, optional + + + Returns + ------- + + occupancy: np.ndarray + in seconds + filtered_firing_rate: np.ndarray(shape=(cell, x, y), dtype=float) + in Hz + + """ + # pixel_width=0.0092, + # field_len = 0.46 + # edges = np.arange(0, field_len + pixel_width, pixel_width) + + x_start = np.nanmin(pos[:, 0]) if x_start is None else x_start + x_stop = np.nanmax(pos[:, 0]) if x_stop is None else x_stop + + y_start = np.nanmin(pos[:, 1]) if y_start is None else y_start + y_stop = np.nanmax(pos[:, 1]) if y_stop is None else y_stop + + edges_x = np.arange(x_start, x_stop, pixel_width) + edges_y = np.arange(y_start, y_stop, pixel_width) + + occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh) + + n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh) + + firing_rate = n_spikes / occupancy # in Hz + firing_rate[np.isnan(firing_rate)] = 0 + + filtered_firing_rate = gaussian_filter(firing_rate, gaussian_sd / pixel_width) + + return occupancy, filtered_firing_rate, [edges_x, edges_y] + + +def compute_2d_place_fields(firing_rate, min_firing_rate=1, thresh=0.2, + min_size=100): + """Compute place fields + + Parameters + ---------- + firing_rate: np.ndarray(NxN, dtype=float) + min_firing_rate: float + in Hz + thresh: float + % of local max + min_size: float + minimum size of place field in pixels + + Returns + ------- + receptive_fields: np.ndarray(NxN, dtype=int) + Each receptive field is labeled with a unique integer + """ + local_maxima_inds = firing_rate == maximum_filter(firing_rate, 3) + n_receptive_fields = 0 + firing_rate = firing_rate.copy() + receptive_fields = {} + for local_max in np.flipud(np.sort(firing_rate[local_maxima_inds])): + labeled_image, num_labels = label(firing_rate > max(local_max * thresh, + min_firing_rate)) + if not num_labels: # nothing above min_firing_thresh + return + for i in range(1, num_labels + 1): + image_label = labeled_image == i + if local_max in firing_rate[image_label]: + break + if np.sum(image_label) >= min_size: + n_receptive_fields += 1 + receptive_fields[image_label] = n_receptive_fields + firing_rate[image_label] = 0 + + return receptive_fields + + class PlaceFieldWidget(widgets.HBox): - def __init__(self, unit: Units, node: pynwb.behavior.SpatialSeries): + def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, index=0, **kwargs): super().__init__() + units = spatial_series.get_ancestor('NWBFile').units + # Initialize receptive fields - self.receptive_fields = np.zeros(firing_rate.shape, dtype=int) # Get pos - pos, unit = get_timeseries_in_units(node) + pos, unit = get_timeseries_in_units(spatial_series) # Get time - pos_tt = get_timeseries_tt(node) + pos_tt = get_timeseries_tt(spatial_series) # Get spikes - spikes = get_spike_times(unit) + spikes = get_spike_times(units, index, [min(pos_tt), max(pos_tt)]) # Pixel width? - pixel_width = + pixel_width = (np.nanmax(pos) - np.nanmin(pos)) / 1000 # Put widget controls here: - # - Gaussian SD - # - Speed threshold - # - Minimum firing rate - # - Place field thresh (% of local max) + # - Gaussian SD + # - Speed threshold + # - Minimum firing rate + # - Place field thresh (% of local max) + + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(pos, pos_tt, spikes, pixel_width) # speed_thresh=0.03, gaussian_sd=0.0184, + # x_start=None, x_stop=None, y_start=None, y_stop=None) - self.compute_2d_firing_rate(pos, pos_tt, spikes, pixel_width) # speed_thresh=0.03, gaussian_sd=0.0184, - # x_start=None, x_stop=None, y_start=None, y_stop=None) + # self.compute_2d_place_fields() # min_firing_rate=1, thresh=0.2, + # # min_size=100): - self.compute_2d_place_fields() # min_firing_rate=1, thresh=0.2, - # min_size=100): + #fig, ax = plt.subplots() + #ax.imshow(filtered_firing_rate, + # extent=[edges_x[0], edges_x[-1], edges_y[0], edges_y[-1]], + # aspect='equal') + # + #self.children = [vis2widget(fig)] - return vis2widget(plotly_show_spatial_trace(self.receptive_fields)) + fig = go.FigureWidget() + fig.add_trace(go.Image(z=filtered_firing_rate)) + self.children = [fig] + + # return vis2widget(plotly_show_spatial_trace(self.receptive_fields)) # Put place field code here - def find_nearest(self, arr, tt): - """Used for picking out elements of a TimeSeries based on spike times - - Parameters - ---------- - arr - tt - - Returns - ------- - - """ - arr = arr[arr > tt[0]] - arr = arr[arr < tt[-1]] - return np.searchsorted(tt, arr) - - def smooth(self, y, box_pts): - """Moving average - - Parameters - ---------- - y - box_pts - - Returns - ------- - - """ - box = np.ones(box_pts) / box_pts - return np.convolve(y, box, mode='same') - - def compute_speed(self, pos, pos_tt, smooth_param=40): - """Compute boolean of whether the speed of the animal was above a threshold - for each time point - - Parameters - ---------- - pos: np.ndarray(dtype=float) - in meters - pos_tt: np.ndarray(dtype=float) - in seconds - smooth_param: float, optional - - Returns - ------- - running: np.ndarray(dtype=bool) - - """ - speed = np.hstack((0, np.sqrt(np.sum(np.diff(pos.T) ** 2, axis=0)) / np.diff(pos_tt))) - return self.smooth(speed, smooth_param) - - - def compute_2d_occupancy(self, pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): - """Computes occupancy per bin in seconds - - Parameters - ---------- - pos: np.ndarray(dtype=float) - in meters - pos_tt: np.ndarray(dtype=float) - in seconds - edges_x: array-like - edges of histogram in meters - edges_y: array-like - edges of histogram in meters - speed_thresh: float, optional - in meters. Default = 3.0 cm/s - - Returns - ------- - occupancy: np.ndarray(dtype=float) - in seconds - running: np.ndarray(dtype=bool) - - """ - - sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) - is_running = self.compute_speed(pos, pos_tt) > speed_thresh - run_pos = pos[is_running, :] - occupancy = np.histogram2d(run_pos[:, 0], - run_pos[:, 1], - [edges_x, edges_y])[0] * sampling_period # in seconds - - return occupancy, is_running - - - def compute_2d_n_spikes(self, pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03): - """Returns speed-gated position during spikes - - Parameters - ---------- - pos: np.ndarray(dtype=float) - (time x 2) in meters - pos_tt: np.ndarray(dtype=float) - (time,) in seconds - spikes: np.ndarray(dtype=float) - in seconds - edges_x: np.ndarray(dtype=float) - edges of histogram in meters - edges_y: np.ndarray(dtype=float) - edges of histogram in meters - speed_thresh: float - in meters. Default = 3.0 cm/s - - Returns - ------- - """ - - is_running = self.compute_speed(pos, pos_tt) > speed_thresh - - spike_pos_inds = self.find_nearest(spikes, pos_tt) - spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] - pos_on_spikes = pos[spike_pos_inds, :] - - n_spikes = np.histogram2d(pos_on_spikes[:, 0], - pos_on_spikes[:, 1], - [edges_x, edges_y])[0] - - return n_spikes - - - def compute_2d_firing_rate(self, pos, pos_tt, spikes, - pixel_width, - speed_thresh=0.03, - gaussian_sd=0.0184, - x_start=None, x_stop=None, - y_start=None, y_stop=None): - """Returns speed-gated occupancy and speed-gated and - Gaussian-filtered firing rate - - Parameters - ---------- - pos: np.ndarray(dtype=float) - (time x 2), in meters - pos_tt: np.ndarray(dtype=float) - (time,) in seconds - spikes: np.ndarray(dtype=float) - in seconds - pixel_width: float - speed_thresh: float, optional - in meters. Default = 3.0 cm/s - gaussian_sd: float, optional - in meters. Default = 1.84 cm - x_start: float, optional - x_stop: float, optional - y_start: float, optional - y_stop: float, optional - - - Returns - ------- - - occupancy: np.ndarray - in seconds - filtered_firing_rate: np.ndarray(shape=(cell, x, y), dtype=float) - in Hz - - """ - # pixel_width=0.0092, - # field_len = 0.46 - # edges = np.arange(0, field_len + pixel_width, pixel_width) - - x_start = np.min(pos[:, 0]) if x_start is None else x_start - x_stop = np.max(pos[:, 0]) if x_stop is None else x_stop - - y_start = np.min(pos[:, 1]) if y_start is None else y_start - y_stop = np.max(pos[:, 1]) if y_stop is None else y_stop - - self.edges_x = np.arange(x_start, x_stop, pixel_width) - self.edges_y = np.arange(y_start, y_stop, pixel_width) - - self.occupancy, running = self.compute_2d_occupancy(pos, pos_tt, self.edges_x, self.edges_y, speed_thresh) - - n_spikes = self.compute_2d_n_spikes(pos, pos_tt, spikes, self.edges_x, self.edges_y, speed_thresh) - - firing_rate = n_spikes / self.occupancy # in Hz - firing_rate[np.isnan(firing_rate)] = 0 - - self.firing_rate = gaussian_filter(firing_rate, gaussian_sd / pixel_width) - - # return occupancy, filtered_firing_rate, [edges_x, edges_y] - - - def compute_2d_place_fields(self, min_firing_rate=1, thresh=0.2, - min_size=100): - """Compute place fields - - Parameters - ---------- - firing_rate: np.ndarray(NxN, dtype=float) - min_firing_rate: float - in Hz - thresh: float - % of local max - min_size: float - minimum size of place field in pixels - - Returns - ------- - receptive_fields: np.ndarray(NxN, dtype=int) - Each receptive field is labeled with a unique integer - """ - firing_rate = self.firing_rate - local_maxima_inds = firing_rate == maximum_filter(firing_rate, 3) - n_receptive_fields = 0 - firing_rate = firing_rate.copy() - for local_max in np.flipud(np.sort(firing_rate[local_maxima_inds])): - labeled_image, num_labels = label(firing_rate > max(local_max * thresh, - min_firing_rate)) - if not num_labels: # nothing above min_firing_thresh - return - for i in range(1, num_labels + 1): - image_label = labeled_image == i - if local_max in firing_rate[image_label]: - break - if np.sum(image_label) >= min_size: - n_receptive_fields += 1 - self.receptive_fields[image_label] = n_receptive_fields - firing_rate[image_label] = 0 - - # return receptive_fields \ No newline at end of file + diff --git a/nwbwidgets/view.py b/nwbwidgets/view.py index a3c4cb96..05b88e7a 100644 --- a/nwbwidgets/view.py +++ b/nwbwidgets/view.py @@ -7,7 +7,7 @@ import zarr from ipywidgets import widgets from ndx_icephys_meta.icephys import SweepSequences -from nwbwidgets import behavior, misc, base, ecephys, image, ophys, icephys, timeseries, file +from nwbwidgets import behavior, misc, base, ecephys, image, ophys, icephys, timeseries, file, placefield # def show_dynamic_table(node: DynamicTable, **kwargs): @@ -47,7 +47,8 @@ def show_dynamic_table(node, **kwargs) -> widgets.Widget: pynwb.behavior.Position: behavior.show_position, pynwb.behavior.SpatialSeries: OrderedDict({ 'over time': timeseries.SeparateTracesPlotlyWidget, - 'trace': behavior.plotly_show_spatial_trace}), + 'trace': behavior.plotly_show_spatial_trace, + 'rate map': placefield.PlaceFieldWidget}), pynwb.image.GrayscaleImage: image.show_grayscale_image, pynwb.image.RGBImage: image.show_rbga_image, pynwb.image.RGBAImage: image.show_rbga_image, From b75288285163f86cd5cb4dfe1a7bf8fa65d43295 Mon Sep 17 00:00:00 2001 From: "!git for-each-ref --format='%(refname:short)' `git symbolic-ref HEAD`" Date: Thu, 29 Oct 2020 11:34:29 -0400 Subject: [PATCH 03/39] interactive widget working --- nwbwidgets/placefield.py | 79 +++++++++++++++++++++++----------------- 1 file changed, 46 insertions(+), 33 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index fa305ea9..4b318c15 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -5,15 +5,15 @@ import matplotlib.pyplot as plt import pynwb -from pynwb.misc import Units -from .behavior import plotly_show_spatial_trace -from .base import vis2widget -from ipywidgets import widgets +from ipywidgets import widgets, BoundedFloatText, Dropdown + +from .utils.widgets import interactive_output from .utils.units import get_spike_times from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt +from .base import vis2widget -import plotly.graph_objects as go +import plotly.graph_objects as go ## To-do @@ -268,48 +268,61 @@ def compute_2d_place_fields(firing_rate, min_firing_rate=1, thresh=0.2, class PlaceFieldWidget(widgets.HBox): - def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, index=0, **kwargs): - + def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): super().__init__() - units = spatial_series.get_ancestor('NWBFile').units + self.units = spatial_series.get_ancestor('NWBFile').units - # Initialize receptive fields - # Get pos - pos, unit = get_timeseries_in_units(spatial_series) - # Get time - pos_tt = get_timeseries_tt(spatial_series) - # Get spikes - spikes = get_spike_times(units, index, [min(pos_tt), max(pos_tt)]) - # Pixel width? - pixel_width = (np.nanmax(pos) - np.nanmin(pos)) / 1000 + self.pos, self.unit = get_timeseries_in_units(spatial_series) + self.pos_tt = get_timeseries_tt(spatial_series) + + self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 # Put widget controls here: - # - Gaussian SD - # - Speed threshold # - Minimum firing rate # - Place field thresh (% of local max) - occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(pos, pos_tt, spikes, pixel_width) # speed_thresh=0.03, gaussian_sd=0.0184, - # x_start=None, x_stop=None, y_start=None, y_stop=None) + bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=np.Inf, desciription='gaussian sd (cm)') + bft_speed = BoundedFloatText(value=0.03, min=0, max=np.Inf, description='speed threshold (cm/s)') + dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') - # self.compute_2d_place_fields() # min_firing_rate=1, thresh=0.2, - # # min_size=100): + self.controls = dict( + gaussian_sd=bft_gaussian, + speed_thresh=bft_speed, + index=dd_unit_select + ) - #fig, ax = plt.subplots() + out_fig = interactive_output(self.do_rate_map, self.controls) - #ax.imshow(filtered_firing_rate, - # extent=[edges_x[0], edges_x[-1], edges_y[0], edges_y[-1]], - # aspect='equal') - # - #self.children = [vis2widget(fig)] + self.children = [ + widgets.VBox([ + bft_gaussian, + bft_speed, + dd_unit_select + ]), + vis2widget(out_fig) + ] - fig = go.FigureWidget() - fig.add_trace(go.Image(z=filtered_firing_rate)) + # fig = go.FigureWidget() + # fig.add_trace(go.Image(z=filtered_firing_rate)) - self.children = [fig] + # self.children = [fig] # return vis2widget(plotly_show_spatial_trace(self.receptive_fields)) - # Put place field code here + def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184): + tmin = min(self.pos_tt) + tmax = max(self.pos_tt) + + spikes = get_spike_times(self.units, index, [tmin, tmax]) + + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( + self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd) + + fig, ax = plt.subplots() + + ax.imshow(filtered_firing_rate, + extent=[edges_x[0], edges_x[-1], edges_y[0], edges_y[-1]], + aspect='equal') + return fig From 8d5f754d142431daf7ce646f8d32494d16b74d4a Mon Sep 17 00:00:00 2001 From: "!git for-each-ref --format='%(refname:short)' `git symbolic-ref HEAD`" Date: Thu, 29 Oct 2020 12:19:33 -0400 Subject: [PATCH 04/39] interactive widget working --- nwbwidgets/placefield.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 4b318c15..fd130972 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -282,8 +282,8 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): # - Minimum firing rate # - Place field thresh (% of local max) - bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=np.Inf, desciription='gaussian sd (cm)') - bft_speed = BoundedFloatText(value=0.03, min=0, max=np.Inf, description='speed threshold (cm/s)') + bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd (cm)') + bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)') dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') self.controls = dict( @@ -321,8 +321,13 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184): fig, ax = plt.subplots() - ax.imshow(filtered_firing_rate, - extent=[edges_x[0], edges_x[-1], edges_y[0], edges_y[-1]], - aspect='equal') + im = ax.imshow(filtered_firing_rate, + extent=[edges_x[0], edges_x[-1], edges_y[0], edges_y[-1]], + aspect='equal') + ax.set_xlabel('x ({})'.format(self.unit)) + ax.set_ylabel('y ({})'.format(self.unit)) + + cbar = plt.colorbar(im) + cbar.ax.set_ylabel('firing rate (Hz)') return fig From 50e284ccba542262b1b4c84f1427f04c5da58851 Mon Sep 17 00:00:00 2001 From: "!git for-each-ref --format='%(refname:short)' `git symbolic-ref HEAD`" Date: Thu, 29 Oct 2020 14:42:05 -0400 Subject: [PATCH 05/39] NaN out unexplored territory --- nwbwidgets/placefield.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index fd130972..1ecf3b4c 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -33,10 +33,10 @@ # [] Modify plotly_show_spatial_trace to plot 2D heatmap representing place fields or create new figure function? # [] Dropdown that controls which unit -# [] Work in buttons / dropdowns / sliders to modify following parameters in place field calculation: +# [x] Work in buttons / dropdowns / sliders to modify following parameters in place field calculation: # [] Different epochs -# [] Gaussian SD -# [] Speed threshold +# [x] Gaussian SD +# [x] Speed threshold # [] Minimum firing rate # [] Place field thresh (% of local max) @@ -201,9 +201,6 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, in Hz """ - # pixel_width=0.0092, - # field_len = 0.46 - # edges = np.arange(0, field_len + pixel_width, pixel_width) x_start = np.nanmin(pos[:, 0]) if x_start is None else x_start x_stop = np.nanmax(pos[:, 0]) if x_stop is None else x_stop @@ -219,10 +216,14 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh) firing_rate = n_spikes / occupancy # in Hz - firing_rate[np.isnan(firing_rate)] = 0 + firing_rate[np.isnan(firing_rate)] = 0 # get rid of NaNs so convolution works filtered_firing_rate = gaussian_filter(firing_rate, gaussian_sd / pixel_width) + # filter occupancy to create a mask so non-explored regions are nan'ed + filtered_occupancy = gaussian_filter(occupancy, gaussian_sd / pixel_width / 8) + filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan + return occupancy, filtered_firing_rate, [edges_x, edges_y] From 03f617b40190ef2f26cbbe84ea1ed3eca35946fc Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 29 Oct 2020 19:18:55 -0400 Subject: [PATCH 06/39] 1D place code widget --- nwbwidgets/__init__.py | 1 - nwbwidgets/placefield.py | 140 ++++++++++++++++++++++++++++++++++++--- nwbwidgets/view.py | 3 +- 3 files changed, 132 insertions(+), 12 deletions(-) diff --git a/nwbwidgets/__init__.py b/nwbwidgets/__init__.py index 3e01a335..7eb06b8f 100644 --- a/nwbwidgets/__init__.py +++ b/nwbwidgets/__init__.py @@ -1,7 +1,6 @@ import plotly.io as pio from .view import nwb2widget, default_neurodata_vis_spec -from .placefield import * # from .ephys_viz_interface import ephys_viz_neurodata_vis_spec diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 1ecf3b4c..3c5c55fe 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -12,7 +12,6 @@ from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt from .base import vis2widget - import plotly.graph_objects as go @@ -273,10 +272,12 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): super().__init__() self.units = spatial_series.get_ancestor('NWBFile').units - - self.pos, self.unit = get_timeseries_in_units(spatial_series) self.pos_tt = get_timeseries_tt(spatial_series) + istart = 0 + istop = None + self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) + self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 # Put widget controls here: @@ -304,13 +305,6 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): vis2widget(out_fig) ] - # fig = go.FigureWidget() - # fig.add_trace(go.Image(z=filtered_firing_rate)) - - # self.children = [fig] - - # return vis2widget(plotly_show_spatial_trace(self.receptive_fields)) - def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184): tmin = min(self.pos_tt) tmax = max(self.pos_tt) @@ -332,3 +326,129 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184): cbar.ax.set_ylabel('firing rate (Hz)') return fig + + +def compute_1d_occupancy(pos, spatial_bins, sampling_rate): + finite_lin_pos = pos[np.isfinite(pos)] + + occupancy = np.histogram( + finite_lin_pos, bins=spatial_bins)[0][:-2] / sampling_rate + + return occupancy + + +def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, + spatial_bin_len=0.0168): + """The occupancy and number of spikes, speed-gated, binned, and smoothed + over position + + Parameters + ---------- + pos: np.ndarray + linearized position + pos_tt: np.ndarray + sample times in seconds + spikes: np.ndarray + for a single cell in seconds + gaussian_sd: float (optional) + in meters. Default = 5.57 cm + spatial_bin_len: float (optional) + in meters. Default = 1.68 cm + + + Returns + ------- + xx: np.ndarray + center of position bins in meters + occupancy: np.ndarray + time in each spatial bin in seconds, during appropriate trials and + while running + filtered_n_spikes: np.ndarray + number of spikes in each spatial bin, during appropriate trials, while + running, and processed with a Gaussian filter + + """ + + spatial_bins = np.arange(np.min(pos), np.max(pos) + spatial_bin_len, + spatial_bin_len) + + sampling_rate = len(pos_tt) / (np.max(pos_tt) - np.min(pos_tt)) + + occupancy = compute_1d_occupancy(pos, spatial_bins, sampling_rate) + + # find pos_tt bin associated with each spike + spike_pos_inds = find_nearest(spikes, pos_tt) + + pos_on_spikes = pos[spike_pos_inds] + finite_pos_on_spikes = pos_on_spikes[np.isfinite(pos_on_spikes)] + + n_spikes = np.histogram(finite_pos_on_spikes, bins=spatial_bins)[0][:-2] + + firing_rate = n_spikes / occupancy + + filtered_firing_rate = gaussian_filter( + firing_rate, gaussian_sd / spatial_bin_len) + xx = spatial_bins[:-3] + (spatial_bins[1] - spatial_bins[0]) / 2 + + return xx, occupancy, filtered_firing_rate + + +class PlaceField_1D_Widget(widgets.HBox): + + def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): + super().__init__() + + self.units = spatial_series.get_ancestor('NWBFile').units + self.pos_tt = get_timeseries_tt(spatial_series) + + istart = 0 + istop = None + self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) + + self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 + + # Put widget controls here: + # - Minimum firing rate + # - Place field thresh (% of local max) + + bft_gaussian = BoundedFloatText(value=0.0557, min=0, max=99999, description='gaussian sd (m)') + bft_spatial_bin_len = BoundedFloatText(value=0.0168, min=0, max=99999, description='spatial bin length (m)') + dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') + + self.controls = dict( + gaussian_sd=bft_gaussian, + spatial_bin_len=bft_spatial_bin_len, + index=dd_unit_select + ) + + out_fig = interactive_output(self.do_1d_rate_map, self.controls) + + self.children = [ + widgets.VBox([ + bft_gaussian, + dd_unit_select + ]), + vis2widget(out_fig) + ] + + def do_1d_rate_map(self, index=0, gaussian_sd=0.0557, spatial_bin_len=0.0168): + tmin = min(self.pos_tt) + tmax = max(self.pos_tt) + + spikes = get_spike_times(self.units, index, [tmin, tmax]) + + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_linear_firing_rate( + self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) + + fig, ax = plt.subplots() + + im = ax.imshow(filtered_firing_rate, + extent=[edges_x[0], edges_x[-1], edges_y[0], edges_y[-1]], + aspect='equal') + ax.set_xlabel('x ({})'.format(self.unit)) + ax.set_ylabel('y ({})'.format(self.unit)) + + cbar = plt.colorbar(im) + cbar.ax.set_ylabel('firing rate (Hz)') + + return fig diff --git a/nwbwidgets/view.py b/nwbwidgets/view.py index 0a17bd8d..d6a7a479 100644 --- a/nwbwidgets/view.py +++ b/nwbwidgets/view.py @@ -47,7 +47,8 @@ def show_dynamic_table(node, **kwargs) -> widgets.Widget: pynwb.behavior.SpatialSeries: OrderedDict({ 'over time': timeseries.SeparateTracesPlotlyWidget, 'trace': behavior.plotly_show_spatial_trace, - 'rate map': placefield.PlaceFieldWidget}), + 'rate map': placefield.PlaceFieldWidget, + '1d rate map': placefield.PlaceField_1D_Widget}), pynwb.image.GrayscaleImage: image.show_grayscale_image, pynwb.image.RGBImage: image.show_rbga_image, pynwb.image.RGBAImage: image.show_rbga_image, From 09e066962301c003634e1799eb5fb38dbe025b85 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 29 Oct 2020 20:35:47 -0400 Subject: [PATCH 07/39] 1D place code widget working --- nwbwidgets/placefield.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 3c5c55fe..2c2dea1c 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -369,10 +369,9 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, """ - spatial_bins = np.arange(np.min(pos), np.max(pos) + spatial_bin_len, - spatial_bin_len) + spatial_bins = np.arange(np.nanmin(pos), np.nanmax(pos) + spatial_bin_len, spatial_bin_len) - sampling_rate = len(pos_tt) / (np.max(pos_tt) - np.min(pos_tt)) + sampling_rate = len(pos_tt) / (np.nanmax(pos_tt) - np.nanmin(pos_tt)) occupancy = compute_1d_occupancy(pos, spatial_bins, sampling_rate) @@ -384,7 +383,7 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, n_spikes = np.histogram(finite_pos_on_spikes, bins=spatial_bins)[0][:-2] - firing_rate = n_spikes / occupancy + firing_rate = np.nan_to_num(n_spikes / occupancy) filtered_firing_rate = gaussian_filter( firing_rate, gaussian_sd / spatial_bin_len) @@ -426,6 +425,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): self.children = [ widgets.VBox([ bft_gaussian, + bft_spatial_bin_len, dd_unit_select ]), vis2widget(out_fig) @@ -437,18 +437,14 @@ def do_1d_rate_map(self, index=0, gaussian_sd=0.0557, spatial_bin_len=0.0168): spikes = get_spike_times(self.units, index, [tmin, tmax]) - occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_linear_firing_rate( + xx, occupancy, filtered_firing_rate = compute_linear_firing_rate( self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) fig, ax = plt.subplots() - im = ax.imshow(filtered_firing_rate, - extent=[edges_x[0], edges_x[-1], edges_y[0], edges_y[-1]], - aspect='equal') + fig = ax.plot(xx, filtered_firing_rate, '-') ax.set_xlabel('x ({})'.format(self.unit)) - ax.set_ylabel('y ({})'.format(self.unit)) + ax.set_ylabel('firing rate (Hz)') - cbar = plt.colorbar(im) - cbar.ax.set_ylabel('firing rate (Hz)') return fig From c346a60a633077c3006f11aece405d9e5a9a06b1 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Fri, 30 Oct 2020 20:04:42 -0400 Subject: [PATCH 08/39] nelpy package integrated into 1D rate map widget --- nwbwidgets/placefield.py | 133 ++++++++++++++++++++++++++++++++------- 1 file changed, 110 insertions(+), 23 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 2c2dea1c..f4305252 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -3,33 +3,41 @@ from scipy.ndimage.filters import gaussian_filter, maximum_filter import matplotlib.pyplot as plt +import nelpy.plotting as nlp import pynwb from ipywidgets import widgets, BoundedFloatText, Dropdown + from .utils.widgets import interactive_output from .utils.units import get_spike_times from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt from .base import vis2widget -import plotly.graph_objects as go - ## To-do -# [] Create PlaceFieldWidget class -# [X] Refactor place field calculation code to deal with nwb data type -# [X] Incorporate place field fxns into class -# [X] Change all internal attributes references -# [X]Change all internal method references - -# [X] Get pos -# [X] Get time -# [X] Get spikes -# [] Get trials / epochs +# [X] Create PlaceFieldWidget class + # [X] Refactor place field calculation code to deal with nwb data type + # [X] Incorporate place field fxns into class + # [X] Change all internal attributes references + # [X]Change all internal method references + + # [X] Get pos + # [X] Get time + # [X] Get spikes + # [] Get trials / epochs + +# [X] Submit draft PR + + # [] 1D Place Field Widget + # [X] Incorporate nelpy package into widget + # [] Add foreign group and sort controller to pick unit groups and ranges? + # [] Normalized firing rate figure? + # [] Add collapsed unit vizualization? + # [] Scale bar? + # [] Sort place cell tuning curves by peak firing rate position? + # [] Color palette control? -# [] Submit draft PR - -# [] Modify plotly_show_spatial_trace to plot 2D heatmap representing place fields or create new figure function? # [] Dropdown that controls which unit # [x] Work in buttons / dropdowns / sliders to modify following parameters in place field calculation: @@ -395,9 +403,18 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, class PlaceField_1D_Widget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): + # foreign_group_and_sort_controller: GroupAndSortController = None, + # group_by=None, + super().__init__() + # if foreign_group_and_sort_controller: + # self.gas = foreign_group_and_sort_controller + # else: + # self.gas = self.make_group_and_sort(group_by=group_by, control_order=False) + self.units = spatial_series.get_ancestor('NWBFile').units + self.pos_tt = get_timeseries_tt(spatial_series) istart = 0 @@ -417,7 +434,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): self.controls = dict( gaussian_sd=bft_gaussian, spatial_bin_len=bft_spatial_bin_len, - index=dd_unit_select + # index=dd_unit_select ) out_fig = interactive_output(self.do_1d_rate_map, self.controls) @@ -426,25 +443,95 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): widgets.VBox([ bft_gaussian, bft_spatial_bin_len, - dd_unit_select + # dd_unit_select ]), vis2widget(out_fig) ] - def do_1d_rate_map(self, index=0, gaussian_sd=0.0557, spatial_bin_len=0.0168): + def do_1d_rate_map(self, gaussian_sd=0.0557, spatial_bin_len=0.0168): tmin = min(self.pos_tt) tmax = max(self.pos_tt) - spikes = get_spike_times(self.units, index, [tmin, tmax]) + index = np.arange(len(self.units)) + spikes = get_spike_times(self.units, index[0], [tmin, tmax]) xx, occupancy, filtered_firing_rate = compute_linear_firing_rate( self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) - fig, ax = plt.subplots() + all_unit_firing_rate = np.zeros([len(self.units), len(xx)]) + all_unit_firing_rate[0] = filtered_firing_rate - fig = ax.plot(xx, filtered_firing_rate, '-') - ax.set_xlabel('x ({})'.format(self.unit)) - ax.set_ylabel('firing rate (Hz)') + for ind in index[1:]: + spikes = get_spike_times(self.units, ind, [tmin, tmax]) + _, _, all_unit_firing_rate[ind] = compute_linear_firing_rate( + self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) + + # npl.set_palette(npl.colors.rainbow) + # with npl.FigureManager(show=True, figsize=(8, 8)) as (fig, ax): + # npl.utils.skip_if_no_output(fig) + fig, ax = plt.subplots() + plot_tuning_curves1D(all_unit_firing_rate, xx, ax=ax, unit_labels=index) + # fig = ax.plot(xx, filtered_firing_rate, '-') + # ax.set_xlabel('x ({})'.format(self.unit)) + # ax.set_ylabel('firing rate (Hz)') return fig + +def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, unit_labels=None, fill=True, color=None): + """ + WARNING! This function is not complete, and hence 'private', + and may be moved somewhere else later on. + + If pad=0 then the y-axis is assumed to be firing rate + """ + xmin = bin_pos[0] + xmax = bin_pos[-1] + xvals = bin_pos + + n_units, n_ext = ratemap.shape + + # if normalize: + # peak_firing_rates = ratemap.max(axis=1) + # ratemap = (ratemap.T / peak_firing_rates).T + + # determine max firing rate + max_firing_rate = ratemap.max() + + if xvals is None: + xvals = np.arange(n_ext) + if xmin is None: + xmin = xvals[0] + if xmax is None: + xmax = xvals[-1] + + for unit, curve in enumerate(ratemap): + if color is None: + line = ax.plot(xvals, unit*pad + curve, zorder=int(10+2*n_units-2*unit)) + else: + line = ax.plot(xvals, unit*pad + curve, zorder=int(10+2*n_units-2*unit), color=color) + if fill: + # Get the color from the current curve + fillcolor = line[0].get_color() + ax.fill_between(xvals, unit*pad, unit*pad + curve, alpha=0.3, color=fillcolor, zorder=int(10+2*n_units-2*unit-1)) + + ax.set_xlim(xmin, xmax) + if pad != 0: + yticks = np.arange(n_units)*pad + 0.5*pad + ax.set_yticks(yticks) + ax.set_yticklabels(unit_labels) + ax.set_xlabel('external variable') + ax.set_ylabel('unit') + nlp.utils.no_yticks(ax) + nlp.utils.clear_left(ax) + else: + if normalize: + ax.set_ylabel('normalized firing rate') + else: + ax.set_ylabel('firing rate [Hz]') + ax.set_ylim(0) + + nlp.utils.clear_top(ax) + nlp.utils.clear_right(ax) + + return ax From 2e92ed4ccaefffcb6ea05fbe856d22db7cf2c717 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Mon, 2 Nov 2020 13:28:41 -0500 Subject: [PATCH 09/39] 1D rate map speed threshold implementation --- nwbwidgets/placefield.py | 56 +++++++--------------------------------- 1 file changed, 10 insertions(+), 46 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index f4305252..d7ea702a 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -234,46 +234,6 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, return occupancy, filtered_firing_rate, [edges_x, edges_y] -def compute_2d_place_fields(firing_rate, min_firing_rate=1, thresh=0.2, - min_size=100): - """Compute place fields - - Parameters - ---------- - firing_rate: np.ndarray(NxN, dtype=float) - min_firing_rate: float - in Hz - thresh: float - % of local max - min_size: float - minimum size of place field in pixels - - Returns - ------- - receptive_fields: np.ndarray(NxN, dtype=int) - Each receptive field is labeled with a unique integer - """ - local_maxima_inds = firing_rate == maximum_filter(firing_rate, 3) - n_receptive_fields = 0 - firing_rate = firing_rate.copy() - receptive_fields = {} - for local_max in np.flipud(np.sort(firing_rate[local_maxima_inds])): - labeled_image, num_labels = label(firing_rate > max(local_max * thresh, - min_firing_rate)) - if not num_labels: # nothing above min_firing_thresh - return - for i in range(1, num_labels + 1): - image_label = labeled_image == i - if local_max in firing_rate[image_label]: - break - if np.sum(image_label) >= min_size: - n_receptive_fields += 1 - receptive_fields[image_label] = n_receptive_fields - firing_rate[image_label] = 0 - - return receptive_fields - - class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): @@ -336,8 +296,11 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184): return fig -def compute_1d_occupancy(pos, spatial_bins, sampling_rate): - finite_lin_pos = pos[np.isfinite(pos)] +def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03): + + is_running = compute_speed(pos, pos_tt) > speed_thresh + run_pos = pos[is_running, :] + finite_lin_pos = run_pos[np.isfinite(run_pos)] occupancy = np.histogram( finite_lin_pos, bins=spatial_bins)[0][:-2] / sampling_rate @@ -346,7 +309,7 @@ def compute_1d_occupancy(pos, spatial_bins, sampling_rate): def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, - spatial_bin_len=0.0168): + spatial_bin_len=0.0168, speed_thresh=0.03): """The occupancy and number of spikes, speed-gated, binned, and smoothed over position @@ -376,16 +339,17 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, running, and processed with a Gaussian filter """ - spatial_bins = np.arange(np.nanmin(pos), np.nanmax(pos) + spatial_bin_len, spatial_bin_len) sampling_rate = len(pos_tt) / (np.nanmax(pos_tt) - np.nanmin(pos_tt)) - occupancy = compute_1d_occupancy(pos, spatial_bins, sampling_rate) + occupancy = compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate) + + is_running = compute_speed(pos, pos_tt) > speed_thresh # find pos_tt bin associated with each spike spike_pos_inds = find_nearest(spikes, pos_tt) - + spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] pos_on_spikes = pos[spike_pos_inds] finite_pos_on_spikes = pos_on_spikes[np.isfinite(pos_on_spikes)] From 7b48dc0334b6a3de5a278f82a0d74a3db1731024 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Tue, 3 Nov 2020 13:31:51 -0500 Subject: [PATCH 10/39] 1D rate map normalization and collapsed view feature addition --- nwbwidgets/placefield.py | 75 ++++++++++++++++++++++++++-------------- 1 file changed, 49 insertions(+), 26 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index d7ea702a..6251a1bd 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -6,13 +6,14 @@ import nelpy.plotting as nlp import pynwb -from ipywidgets import widgets, BoundedFloatText, Dropdown +from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox from .utils.widgets import interactive_output from .utils.units import get_spike_times from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt from .base import vis2widget +from .controllers import GroupAndSortController ## To-do @@ -31,6 +32,7 @@ # [] 1D Place Field Widget # [X] Incorporate nelpy package into widget + # [X] Speed threshold implementation # [] Add foreign group and sort controller to pick unit groups and ranges? # [] Normalized firing rate figure? # [] Add collapsed unit vizualization? @@ -366,21 +368,21 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, class PlaceField_1D_Widget(widgets.HBox): - def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): - # foreign_group_and_sort_controller: GroupAndSortController = None, - # group_by=None, + def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, + foreign_group_and_sort_controller: GroupAndSortController = None, + group_by=None, **kwargs): super().__init__() - # if foreign_group_and_sort_controller: - # self.gas = foreign_group_and_sort_controller - # else: - # self.gas = self.make_group_and_sort(group_by=group_by, control_order=False) - self.units = spatial_series.get_ancestor('NWBFile').units self.pos_tt = get_timeseries_tt(spatial_series) + if foreign_group_and_sort_controller: + self.gas = foreign_group_and_sort_controller + else: + self.gas = self.make_group_and_sort(group_by=group_by, control_order=False) + istart = 0 istop = None self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) @@ -393,48 +395,66 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): bft_gaussian = BoundedFloatText(value=0.0557, min=0, max=99999, description='gaussian sd (m)') bft_spatial_bin_len = BoundedFloatText(value=0.0168, min=0, max=99999, description='spatial bin length (m)') - dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') + cb_normalize_select = Checkbox(value=False, description='normalize', indent=False) + cb_collapsed_select = Checkbox(value=False, description='collapsed', indent=False) self.controls = dict( + gas=self.gas, gaussian_sd=bft_gaussian, spatial_bin_len=bft_spatial_bin_len, - # index=dd_unit_select + normalize=cb_normalize_select, + collapsed=cb_collapsed_select ) out_fig = interactive_output(self.do_1d_rate_map, self.controls) - + checkboxes = widgets.HBox([cb_normalize_select, cb_collapsed_select]) self.children = [ widgets.VBox([ bft_gaussian, bft_spatial_bin_len, - # dd_unit_select + checkboxes, + self.gas, ]), vis2widget(out_fig) ] - def do_1d_rate_map(self, gaussian_sd=0.0557, spatial_bin_len=0.0168): + def make_group_and_sort(self, group_by=None, control_order=True): + return GroupAndSortController(self.units, group_by=group_by, control_order=control_order) + + def do_1d_rate_map(self, units_window=None, order=None, group_inds=None, labels=None, normalize=False, + collapsed=False, gaussian_sd=0.0557, spatial_bin_len=0.0168, **kwargs): tmin = min(self.pos_tt) tmax = max(self.pos_tt) - index = np.arange(len(self.units)) + + print('Order is {}'.format(order)) + print('Units_window is {}'.format(units_window)) + print('Group_inds is {}'.format(group_inds)) + print('Labels is {}'.format(labels)) + if order.any() == None: + index = np.arange(0, len(self.units)) + else: + index = order spikes = get_spike_times(self.units, index[0], [tmin, tmax]) xx, occupancy, filtered_firing_rate = compute_linear_firing_rate( self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) - all_unit_firing_rate = np.zeros([len(self.units), len(xx)]) + all_unit_firing_rate = np.zeros([len(index), len(xx)]) all_unit_firing_rate[0] = filtered_firing_rate - + firing_rate_ind = 0 for ind in index[1:]: spikes = get_spike_times(self.units, ind, [tmin, tmax]) - _, _, all_unit_firing_rate[ind] = compute_linear_firing_rate( + _, _, all_unit_firing_rate[firing_rate_ind] = compute_linear_firing_rate( self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) + firing_rate_ind += 1 # npl.set_palette(npl.colors.rainbow) # with npl.FigureManager(show=True, figsize=(8, 8)) as (fig, ax): # npl.utils.skip_if_no_output(fig) fig, ax = plt.subplots() - plot_tuning_curves1D(all_unit_firing_rate, xx, ax=ax, unit_labels=index) + plot_tuning_curves1D(all_unit_firing_rate, xx, ax=ax, unit_labels=index, normalize=normalize, + collapsed=collapsed) # fig = ax.plot(xx, filtered_firing_rate, '-') # ax.set_xlabel('x ({})'.format(self.unit)) @@ -442,7 +462,8 @@ def do_1d_rate_map(self, gaussian_sd=0.0557, spatial_bin_len=0.0168): return fig -def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, unit_labels=None, fill=True, color=None): +def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, unit_labels=None, fill=True, color=None, + collapsed=False): """ WARNING! This function is not complete, and hence 'private', and may be moved somewhere else later on. @@ -455,12 +476,14 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni n_units, n_ext = ratemap.shape - # if normalize: - # peak_firing_rates = ratemap.max(axis=1) - # ratemap = (ratemap.T / peak_firing_rates).T - - # determine max firing rate - max_firing_rate = ratemap.max() + if normalize: + peak_firing_rates = ratemap.max(axis=1) + ratemap = (ratemap.T / peak_firing_rates).T + pad = 1 + # # determine max firing rate + # max_firing_rate = ratemap.max() + if collapsed: + pad = 0 if xvals is None: xvals = np.arange(n_ext) From ece7795e59431b211d8d9c99b2a939a71b2d1a03 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Tue, 3 Nov 2020 18:55:28 -0500 Subject: [PATCH 11/39] - Added Numpy docstring for plot_tuning_curves1D fxn - Checkbox for normalize added in last commit --- nwbwidgets/placefield.py | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 6251a1bd..06657d57 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -426,7 +426,6 @@ def do_1d_rate_map(self, units_window=None, order=None, group_inds=None, labels= tmin = min(self.pos_tt) tmax = max(self.pos_tt) - print('Order is {}'.format(order)) print('Units_window is {}'.format(units_window)) print('Group_inds is {}'.format(group_inds)) @@ -465,10 +464,33 @@ def do_1d_rate_map(self, units_window=None, order=None, group_inds=None, labels= def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, unit_labels=None, fill=True, color=None, collapsed=False): """ - WARNING! This function is not complete, and hence 'private', - and may be moved somewhere else later on. - If pad=0 then the y-axis is assumed to be firing rate + Parameters + ---------- + ratemap: array-like + An array of dim: [number of units, bin positions] with the spike rates for a unit, at every pos, in each row + bin_pos: array-like + An array representing the bin positions of ratemap for each column + ax: matplotlib.pyplot.ax + Axes object for the figure on which the ratemaps will be plotted + normalize: bool + default = False + Input to determine whether or not to normalize firing rates + pad: int + default = 10 + Changes to 0 if 'collapsed' is true + Amount of space to put between each unit (i.e. row) in the figure + unit_labels: array-like + Unit ids for each unit in ratemap + collapsed: bool + default = False + Determines whether to plot the ratemaps with zero padding, i.e. at the same y coordinate, on the ratemap + fill: bool, optional + + Returns + ------- + matplotlib.pyplot.ax + """ xmin = bin_pos[0] xmax = bin_pos[-1] From 07a7ae7091cd372c4d8ff3102723a424a8a8dc8f Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Tue, 3 Nov 2020 19:05:30 -0500 Subject: [PATCH 12/39] Removed nelpy dependency --- nwbwidgets/placefield.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 06657d57..e6203b1c 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -3,7 +3,6 @@ from scipy.ndimage.filters import gaussian_filter, maximum_filter import matplotlib.pyplot as plt -import nelpy.plotting as nlp import pynwb from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox @@ -531,8 +530,9 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni ax.set_yticklabels(unit_labels) ax.set_xlabel('external variable') ax.set_ylabel('unit') - nlp.utils.no_yticks(ax) - nlp.utils.clear_left(ax) + ax.tick_params(axis=u'y', which=u'both', length=0) + ax.spines['left'].set_color('none') + ax.yaxis.set_ticks_position('right') else: if normalize: ax.set_ylabel('normalized firing rate') @@ -540,7 +540,9 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni ax.set_ylabel('firing rate [Hz]') ax.set_ylim(0) - nlp.utils.clear_top(ax) - nlp.utils.clear_right(ax) + ax.spines['top'].set_color('none') + ax.xaxis.set_ticks_position('bottom') + ax.spines['right'].set_color('none') + ax.yaxis.set_ticks_position('left') return ax From 7f8a790d90b444aa3be3c19ef69996741fa3cf62 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Tue, 3 Nov 2020 19:16:08 -0500 Subject: [PATCH 13/39] Moved analysis code to analysis/placefields.py --- nwbwidgets/analysis/placefields.py | 255 +++++++++++++++++++++++++++++ nwbwidgets/placefield.py | 255 +---------------------------- 2 files changed, 256 insertions(+), 254 deletions(-) create mode 100644 nwbwidgets/analysis/placefields.py diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py new file mode 100644 index 00000000..d44897c1 --- /dev/null +++ b/nwbwidgets/analysis/placefields.py @@ -0,0 +1,255 @@ +import numpy as np +from scipy.ndimage.filters import gaussian_filter, maximum_filter + +def find_nearest(arr, tt): + """Used for picking out elements of a TimeSeries based on spike times + + Parameters + ---------- + arr + tt + + Returns + ------- + + """ + arr = arr[arr > tt[0]] + arr = arr[arr < tt[-1]] + return np.searchsorted(tt, arr) + + +def smooth(y, box_pts): + """Moving average + + Parameters + ---------- + y + box_pts + + Returns + ------- + + """ + box = np.ones(box_pts) / box_pts + return np.convolve(y, box, mode='same') + + +def compute_speed(pos, pos_tt, smooth_param=40): + """Compute boolean of whether the speed of the animal was above a threshold + for each time point + + Parameters + ---------- + pos: np.ndarray(dtype=float) + in meters + pos_tt: np.ndarray(dtype=float) + in seconds + smooth_param: float, optional + + Returns + ------- + running: np.ndarray(dtype=bool) + + """ + speed = np.hstack((0, np.sqrt(np.sum(np.diff(pos.T) ** 2, axis=0)) / np.diff(pos_tt))) + return smooth(speed, smooth_param) + + +def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): + """Computes occupancy per bin in seconds + + Parameters + ---------- + pos: np.ndarray(dtype=float) + in meters + pos_tt: np.ndarray(dtype=float) + in seconds + edges_x: array-like + edges of histogram in meters + edges_y: array-like + edges of histogram in meters + speed_thresh: float, optional + in meters. Default = 3.0 cm/s + + Returns + ------- + occupancy: np.ndarray(dtype=float) + in seconds + running: np.ndarray(dtype=bool) + + """ + + sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) + is_running = compute_speed(pos, pos_tt) > speed_thresh + run_pos = pos[is_running, :] + occupancy = np.histogram2d(run_pos[:, 0], + run_pos[:, 1], + [edges_x, edges_y])[0] * sampling_period # in seconds + + return occupancy, is_running + + +def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03): + """Returns speed-gated position during spikes + + Parameters + ---------- + pos: np.ndarray(dtype=float) + (time x 2) in meters + pos_tt: np.ndarray(dtype=float) + (time,) in seconds + spikes: np.ndarray(dtype=float) + in seconds + edges_x: np.ndarray(dtype=float) + edges of histogram in meters + edges_y: np.ndarray(dtype=float) + edges of histogram in meters + speed_thresh: float + in meters. Default = 3.0 cm/s + + Returns + ------- + """ + + is_running = compute_speed(pos, pos_tt) > speed_thresh + + spike_pos_inds = find_nearest(spikes, pos_tt) + spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] + pos_on_spikes = pos[spike_pos_inds, :] + + n_spikes = np.histogram2d(pos_on_spikes[:, 0], + pos_on_spikes[:, 1], + [edges_x, edges_y])[0] + + return n_spikes + + +def compute_2d_firing_rate(pos, pos_tt, spikes, + pixel_width, + speed_thresh=0.03, + gaussian_sd=0.0184, + x_start=None, x_stop=None, + y_start=None, y_stop=None): + """Returns speed-gated occupancy and speed-gated and + Gaussian-filtered firing rate + + Parameters + ---------- + pos: np.ndarray(dtype=float) + (time x 2), in meters + pos_tt: np.ndarray(dtype=float) + (time,) in seconds + spikes: np.ndarray(dtype=float) + in seconds + pixel_width: float + speed_thresh: float, optional + in meters. Default = 3.0 cm/s + gaussian_sd: float, optional + in meters. Default = 1.84 cm + x_start: float, optional + x_stop: float, optional + y_start: float, optional + y_stop: float, optional + + + Returns + ------- + + occupancy: np.ndarray + in seconds + filtered_firing_rate: np.ndarray(shape=(cell, x, y), dtype=float) + in Hz + + """ + + x_start = np.nanmin(pos[:, 0]) if x_start is None else x_start + x_stop = np.nanmax(pos[:, 0]) if x_stop is None else x_stop + + y_start = np.nanmin(pos[:, 1]) if y_start is None else y_start + y_stop = np.nanmax(pos[:, 1]) if y_stop is None else y_stop + + edges_x = np.arange(x_start, x_stop, pixel_width) + edges_y = np.arange(y_start, y_stop, pixel_width) + + occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh) + + n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh) + + firing_rate = n_spikes / occupancy # in Hz + firing_rate[np.isnan(firing_rate)] = 0 # get rid of NaNs so convolution works + + filtered_firing_rate = gaussian_filter(firing_rate, gaussian_sd / pixel_width) + + # filter occupancy to create a mask so non-explored regions are nan'ed + filtered_occupancy = gaussian_filter(occupancy, gaussian_sd / pixel_width / 8) + filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan + + return occupancy, filtered_firing_rate, [edges_x, edges_y] + + +def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03): + + is_running = compute_speed(pos, pos_tt) > speed_thresh + run_pos = pos[is_running, :] + finite_lin_pos = run_pos[np.isfinite(run_pos)] + + occupancy = np.histogram( + finite_lin_pos, bins=spatial_bins)[0][:-2] / sampling_rate + + return occupancy + + +def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, + spatial_bin_len=0.0168, speed_thresh=0.03): + """The occupancy and number of spikes, speed-gated, binned, and smoothed + over position + + Parameters + ---------- + pos: np.ndarray + linearized position + pos_tt: np.ndarray + sample times in seconds + spikes: np.ndarray + for a single cell in seconds + gaussian_sd: float (optional) + in meters. Default = 5.57 cm + spatial_bin_len: float (optional) + in meters. Default = 1.68 cm + + + Returns + ------- + xx: np.ndarray + center of position bins in meters + occupancy: np.ndarray + time in each spatial bin in seconds, during appropriate trials and + while running + filtered_n_spikes: np.ndarray + number of spikes in each spatial bin, during appropriate trials, while + running, and processed with a Gaussian filter + + """ + spatial_bins = np.arange(np.nanmin(pos), np.nanmax(pos) + spatial_bin_len, spatial_bin_len) + + sampling_rate = len(pos_tt) / (np.nanmax(pos_tt) - np.nanmin(pos_tt)) + + occupancy = compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate) + + is_running = compute_speed(pos, pos_tt) > speed_thresh + + # find pos_tt bin associated with each spike + spike_pos_inds = find_nearest(spikes, pos_tt) + spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] + pos_on_spikes = pos[spike_pos_inds] + finite_pos_on_spikes = pos_on_spikes[np.isfinite(pos_on_spikes)] + + n_spikes = np.histogram(finite_pos_on_spikes, bins=spatial_bins)[0][:-2] + + firing_rate = np.nan_to_num(n_spikes / occupancy) + + filtered_firing_rate = gaussian_filter( + firing_rate, gaussian_sd / spatial_bin_len) + xx = spatial_bins[:-3] + (spatial_bins[1] - spatial_bins[0]) / 2 + + return xx, occupancy, filtered_firing_rate \ No newline at end of file diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index e6203b1c..302dd868 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -7,6 +7,7 @@ import pynwb from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox +from .analysis.placefields import compute_2d_firing_rate, compute_linear_firing_rate from .utils.widgets import interactive_output from .utils.units import get_spike_times @@ -49,192 +50,6 @@ # [] Place field thresh (% of local max) # Put widget rendering here -def find_nearest(arr, tt): - """Used for picking out elements of a TimeSeries based on spike times - - Parameters - ---------- - arr - tt - - Returns - ------- - - """ - arr = arr[arr > tt[0]] - arr = arr[arr < tt[-1]] - return np.searchsorted(tt, arr) - - -def smooth(y, box_pts): - """Moving average - - Parameters - ---------- - y - box_pts - - Returns - ------- - - """ - box = np.ones(box_pts) / box_pts - return np.convolve(y, box, mode='same') - - -def compute_speed(pos, pos_tt, smooth_param=40): - """Compute boolean of whether the speed of the animal was above a threshold - for each time point - - Parameters - ---------- - pos: np.ndarray(dtype=float) - in meters - pos_tt: np.ndarray(dtype=float) - in seconds - smooth_param: float, optional - - Returns - ------- - running: np.ndarray(dtype=bool) - - """ - speed = np.hstack((0, np.sqrt(np.sum(np.diff(pos.T) ** 2, axis=0)) / np.diff(pos_tt))) - return smooth(speed, smooth_param) - - -def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): - """Computes occupancy per bin in seconds - - Parameters - ---------- - pos: np.ndarray(dtype=float) - in meters - pos_tt: np.ndarray(dtype=float) - in seconds - edges_x: array-like - edges of histogram in meters - edges_y: array-like - edges of histogram in meters - speed_thresh: float, optional - in meters. Default = 3.0 cm/s - - Returns - ------- - occupancy: np.ndarray(dtype=float) - in seconds - running: np.ndarray(dtype=bool) - - """ - - sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) - is_running = compute_speed(pos, pos_tt) > speed_thresh - run_pos = pos[is_running, :] - occupancy = np.histogram2d(run_pos[:, 0], - run_pos[:, 1], - [edges_x, edges_y])[0] * sampling_period # in seconds - - return occupancy, is_running - - -def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03): - """Returns speed-gated position during spikes - - Parameters - ---------- - pos: np.ndarray(dtype=float) - (time x 2) in meters - pos_tt: np.ndarray(dtype=float) - (time,) in seconds - spikes: np.ndarray(dtype=float) - in seconds - edges_x: np.ndarray(dtype=float) - edges of histogram in meters - edges_y: np.ndarray(dtype=float) - edges of histogram in meters - speed_thresh: float - in meters. Default = 3.0 cm/s - - Returns - ------- - """ - - is_running = compute_speed(pos, pos_tt) > speed_thresh - - spike_pos_inds = find_nearest(spikes, pos_tt) - spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] - pos_on_spikes = pos[spike_pos_inds, :] - - n_spikes = np.histogram2d(pos_on_spikes[:, 0], - pos_on_spikes[:, 1], - [edges_x, edges_y])[0] - - return n_spikes - - -def compute_2d_firing_rate(pos, pos_tt, spikes, - pixel_width, - speed_thresh=0.03, - gaussian_sd=0.0184, - x_start=None, x_stop=None, - y_start=None, y_stop=None): - """Returns speed-gated occupancy and speed-gated and - Gaussian-filtered firing rate - - Parameters - ---------- - pos: np.ndarray(dtype=float) - (time x 2), in meters - pos_tt: np.ndarray(dtype=float) - (time,) in seconds - spikes: np.ndarray(dtype=float) - in seconds - pixel_width: float - speed_thresh: float, optional - in meters. Default = 3.0 cm/s - gaussian_sd: float, optional - in meters. Default = 1.84 cm - x_start: float, optional - x_stop: float, optional - y_start: float, optional - y_stop: float, optional - - - Returns - ------- - - occupancy: np.ndarray - in seconds - filtered_firing_rate: np.ndarray(shape=(cell, x, y), dtype=float) - in Hz - - """ - - x_start = np.nanmin(pos[:, 0]) if x_start is None else x_start - x_stop = np.nanmax(pos[:, 0]) if x_stop is None else x_stop - - y_start = np.nanmin(pos[:, 1]) if y_start is None else y_start - y_stop = np.nanmax(pos[:, 1]) if y_stop is None else y_stop - - edges_x = np.arange(x_start, x_stop, pixel_width) - edges_y = np.arange(y_start, y_stop, pixel_width) - - occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh) - - n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh) - - firing_rate = n_spikes / occupancy # in Hz - firing_rate[np.isnan(firing_rate)] = 0 # get rid of NaNs so convolution works - - filtered_firing_rate = gaussian_filter(firing_rate, gaussian_sd / pixel_width) - - # filter occupancy to create a mask so non-explored regions are nan'ed - filtered_occupancy = gaussian_filter(occupancy, gaussian_sd / pixel_width / 8) - filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan - - return occupancy, filtered_firing_rate, [edges_x, edges_y] - - class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): @@ -297,74 +112,6 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184): return fig -def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03): - - is_running = compute_speed(pos, pos_tt) > speed_thresh - run_pos = pos[is_running, :] - finite_lin_pos = run_pos[np.isfinite(run_pos)] - - occupancy = np.histogram( - finite_lin_pos, bins=spatial_bins)[0][:-2] / sampling_rate - - return occupancy - - -def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, - spatial_bin_len=0.0168, speed_thresh=0.03): - """The occupancy and number of spikes, speed-gated, binned, and smoothed - over position - - Parameters - ---------- - pos: np.ndarray - linearized position - pos_tt: np.ndarray - sample times in seconds - spikes: np.ndarray - for a single cell in seconds - gaussian_sd: float (optional) - in meters. Default = 5.57 cm - spatial_bin_len: float (optional) - in meters. Default = 1.68 cm - - - Returns - ------- - xx: np.ndarray - center of position bins in meters - occupancy: np.ndarray - time in each spatial bin in seconds, during appropriate trials and - while running - filtered_n_spikes: np.ndarray - number of spikes in each spatial bin, during appropriate trials, while - running, and processed with a Gaussian filter - - """ - spatial_bins = np.arange(np.nanmin(pos), np.nanmax(pos) + spatial_bin_len, spatial_bin_len) - - sampling_rate = len(pos_tt) / (np.nanmax(pos_tt) - np.nanmin(pos_tt)) - - occupancy = compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate) - - is_running = compute_speed(pos, pos_tt) > speed_thresh - - # find pos_tt bin associated with each spike - spike_pos_inds = find_nearest(spikes, pos_tt) - spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] - pos_on_spikes = pos[spike_pos_inds] - finite_pos_on_spikes = pos_on_spikes[np.isfinite(pos_on_spikes)] - - n_spikes = np.histogram(finite_pos_on_spikes, bins=spatial_bins)[0][:-2] - - firing_rate = np.nan_to_num(n_spikes / occupancy) - - filtered_firing_rate = gaussian_filter( - firing_rate, gaussian_sd / spatial_bin_len) - xx = spatial_bins[:-3] + (spatial_bins[1] - spatial_bins[0]) / 2 - - return xx, occupancy, filtered_firing_rate - - class PlaceField_1D_Widget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, From 55611dd1e463272b638005c064ae3ca9f712c171 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Wed, 4 Nov 2020 13:34:13 -0500 Subject: [PATCH 14/39] Dynamically route rate maps to correspond to dimensionality of the data --- nwbwidgets/placefield.py | 17 ++++++++++++++--- nwbwidgets/view.py | 3 +-- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 302dd868..523ecfee 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -33,9 +33,9 @@ # [] 1D Place Field Widget # [X] Incorporate nelpy package into widget # [X] Speed threshold implementation - # [] Add foreign group and sort controller to pick unit groups and ranges? - # [] Normalized firing rate figure? - # [] Add collapsed unit vizualization? + # [X] Add foreign group and sort controller to pick unit groups and ranges? + # [X] Normalized firing rate figure? + # [X] Add collapsed unit vizualization? # [] Scale bar? # [] Sort place cell tuning curves by peak firing rate position? # [] Color palette control? @@ -49,6 +49,17 @@ # [] Minimum firing rate # [] Place field thresh (% of local max) +def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): + print(spatial_series.data.shape) + if spatial_series.data.shape[1] == 2: + return PlaceFieldWidget(spatial_series) + elif spatial_series.data.shape[1] == 1: + return PlaceField_1D_Widget(spatial_series) + else: + print('Spatial series exceeds dimensionality for visualization') + return + + # Put widget rendering here class PlaceFieldWidget(widgets.HBox): diff --git a/nwbwidgets/view.py b/nwbwidgets/view.py index d6a7a479..4abedf9c 100644 --- a/nwbwidgets/view.py +++ b/nwbwidgets/view.py @@ -47,8 +47,7 @@ def show_dynamic_table(node, **kwargs) -> widgets.Widget: pynwb.behavior.SpatialSeries: OrderedDict({ 'over time': timeseries.SeparateTracesPlotlyWidget, 'trace': behavior.plotly_show_spatial_trace, - 'rate map': placefield.PlaceFieldWidget, - '1d rate map': placefield.PlaceField_1D_Widget}), + 'rate map': placefield.route_placefield}), pynwb.image.GrayscaleImage: image.show_grayscale_image, pynwb.image.RGBImage: image.show_rbga_image, pynwb.image.RGBAImage: image.show_rbga_image, From 0ff00372e1e92ec8eeb92b9f8ceb3de41b1637d0 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Wed, 4 Nov 2020 16:49:40 -0500 Subject: [PATCH 15/39] Add optional TimeSeries of velocity as input --- nwbwidgets/analysis/placefields.py | 22 +++++++++++++++------- nwbwidgets/placefield.py | 29 ++++++++++++++++++++++------- 2 files changed, 37 insertions(+), 14 deletions(-) diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index d44897c1..ee6a57bb 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -55,7 +55,7 @@ def compute_speed(pos, pos_tt, smooth_param=40): return smooth(speed, smooth_param) -def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): +def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, velocity=[]): """Computes occupancy per bin in seconds Parameters @@ -80,7 +80,11 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): """ sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) - is_running = compute_speed(pos, pos_tt) > speed_thresh + if len(velocity) == 0: + is_running = compute_speed(pos, pos_tt) > speed_thresh + else: + is_running = self.velocity > speed_thresh + run_pos = pos[is_running, :] occupancy = np.histogram2d(run_pos[:, 0], run_pos[:, 1], @@ -89,7 +93,7 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03): return occupancy, is_running -def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03): +def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03, velocity=[]): """Returns speed-gated position during spikes Parameters @@ -111,7 +115,10 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 ------- """ - is_running = compute_speed(pos, pos_tt) > speed_thresh + if len(velocity) == 0: + is_running = compute_speed(pos, pos_tt) > speed_thresh + else: + is_running = velocity > speed_thresh spike_pos_inds = find_nearest(spikes, pos_tt) spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] @@ -129,7 +136,8 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, speed_thresh=0.03, gaussian_sd=0.0184, x_start=None, x_stop=None, - y_start=None, y_stop=None): + y_start=None, y_stop=None, + velocity=[]): """Returns speed-gated occupancy and speed-gated and Gaussian-filtered firing rate @@ -171,9 +179,9 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, edges_x = np.arange(x_start, x_stop, pixel_width) edges_y = np.arange(y_start, y_stop, pixel_width) - occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh) + occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh, velocity) - n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh) + n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh, velocity) firing_rate = n_spikes / occupancy # in Hz firing_rate[np.isnan(firing_rate)] = 0 # get rid of NaNs so convolution works diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 523ecfee..53c20d08 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -5,7 +5,7 @@ import matplotlib.pyplot as plt import pynwb -from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox +from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox, fixed from .analysis.placefields import compute_2d_firing_rate, compute_linear_firing_rate @@ -66,11 +66,18 @@ class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): super().__init__() + + if hasattr(spatial_series.get_ancestor('NWBFile'),'velocity'): + velocity = spatial_series.get_ancestor('NWBFile').velocity + else: + velocity = [] + self.units = spatial_series.get_ancestor('NWBFile').units self.pos_tt = get_timeseries_tt(spatial_series) istart = 0 istop = None + self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 @@ -82,11 +89,14 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd (cm)') bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)') dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') + cb_velocity = Checkbox(value=False, description='use velocity', indent=False) self.controls = dict( gaussian_sd=bft_gaussian, speed_thresh=bft_speed, - index=dd_unit_select + index=dd_unit_select, + use_velocity=cb_velocity, + velocity=fixed(velocity) ) out_fig = interactive_output(self.do_rate_map, self.controls) @@ -95,19 +105,24 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): widgets.VBox([ bft_gaussian, bft_speed, - dd_unit_select + dd_unit_select, + cb_velocity ]), vis2widget(out_fig) ] - def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184): + def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_velocity=False, velocity=[]): tmin = min(self.pos_tt) tmax = max(self.pos_tt) spikes = get_spike_times(self.units, index, [tmin, tmax]) - - occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( - self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd) + if use_velocity == False: + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( + self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd) + else: + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( + self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, + velocity=velocity) fig, ax = plt.subplots() From 81c40b2b995dabcab1f9f1ea5d916863be06ba20 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:37:06 -0500 Subject: [PATCH 16/39] Update nwbwidgets/placefield.py Co-authored-by: Ben Dichter --- nwbwidgets/placefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 53c20d08..74544a80 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -63,7 +63,7 @@ def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): # Put widget rendering here class PlaceFieldWidget(widgets.HBox): - def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, **kwargs): + def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb.TimeSeries = None, **kwargs): super().__init__() From 2601b86bd8b7b4044868621edec6e17fe0ac8ec7 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:38:23 -0500 Subject: [PATCH 17/39] Update nwbwidgets/placefield.py Co-authored-by: Ben Dichter --- nwbwidgets/placefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 74544a80..f0226a35 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -74,7 +74,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb self.units = spatial_series.get_ancestor('NWBFile').units self.pos_tt = get_timeseries_tt(spatial_series) - + self.velocity = velocity istart = 0 istop = None From 7290f008090f72c3146f2662d3396907f95edb54 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:38:42 -0500 Subject: [PATCH 18/39] Update nwbwidgets/placefield.py Co-authored-by: Ben Dichter --- nwbwidgets/placefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index f0226a35..d78fef01 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -118,7 +118,7 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_veloci spikes = get_spike_times(self.units, index, [tmin, tmax]) if use_velocity == False: occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( - self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd) + self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, velocity=self.velocity) else: occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, From 33f3c70d9c78cc70ca03a3dc513ac5e08f5d5bf2 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:39:05 -0500 Subject: [PATCH 19/39] Update nwbwidgets/placefield.py Co-authored-by: Ben Dichter --- nwbwidgets/placefield.py | 1 - 1 file changed, 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index d78fef01..180da8ab 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -96,7 +96,6 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb speed_thresh=bft_speed, index=dd_unit_select, use_velocity=cb_velocity, - velocity=fixed(velocity) ) out_fig = interactive_output(self.do_rate_map, self.controls) From 0f77f2db21ac3f23ca93f45657578835dfc02b39 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:39:22 -0500 Subject: [PATCH 20/39] Update nwbwidgets/placefield.py Co-authored-by: Ben Dichter --- nwbwidgets/placefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 180da8ab..fc8844e2 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -110,7 +110,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb vis2widget(out_fig) ] - def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_velocity=False, velocity=[]): + def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_velocity=False): tmin = min(self.pos_tt) tmax = max(self.pos_tt) From 45ae3ec37bd79891f50719baddf76ac2db423561 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:40:37 -0500 Subject: [PATCH 21/39] Update nwbwidgets/placefield.py Co-authored-by: Ben Dichter --- nwbwidgets/placefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index fc8844e2..41876599 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -121,7 +121,7 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_veloci else: occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, - velocity=velocity) + velocity=self.velocity) fig, ax = plt.subplots() From c7b020816d697e63bc033bf4a049edda0a301995 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:40:52 -0500 Subject: [PATCH 22/39] Update nwbwidgets/placefield.py Co-authored-by: Ben Dichter --- nwbwidgets/placefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 41876599..13c0869f 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -201,7 +201,7 @@ def do_1d_rate_map(self, units_window=None, order=None, group_inds=None, labels= print('Units_window is {}'.format(units_window)) print('Group_inds is {}'.format(group_inds)) print('Labels is {}'.format(labels)) - if order.any() == None: + if order is None: index = np.arange(0, len(self.units)) else: index = order From fdb0d5bde0bc715829f881590a2f78986d004c18 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 12:43:14 -0500 Subject: [PATCH 23/39] Corrected velocity as optional input --- nwbwidgets/analysis/placefields.py | 36 +++++++++++++++++++----------- nwbwidgets/view.py | 3 ++- 2 files changed, 25 insertions(+), 14 deletions(-) diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index ee6a57bb..dda1691c 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -55,7 +55,7 @@ def compute_speed(pos, pos_tt, smooth_param=40): return smooth(speed, smooth_param) -def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, velocity=[]): +def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, velocity=None): """Computes occupancy per bin in seconds Parameters @@ -70,6 +70,8 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc edges of histogram in meters speed_thresh: float, optional in meters. Default = 3.0 cm/s + velocity: np.ndarray(dtype=float) + pre-computed velocity Returns ------- @@ -80,10 +82,10 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc """ sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) - if len(velocity) == 0: + if velocity is None: is_running = compute_speed(pos, pos_tt) > speed_thresh else: - is_running = self.velocity > speed_thresh + is_running = np.linalg.norm(self.velocity) > speed_thresh run_pos = pos[is_running, :] occupancy = np.histogram2d(run_pos[:, 0], @@ -93,7 +95,7 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc return occupancy, is_running -def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03, velocity=[]): +def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03, velocity=None): """Returns speed-gated position during spikes Parameters @@ -110,15 +112,17 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 edges of histogram in meters speed_thresh: float in meters. Default = 3.0 cm/s + velocity: np.ndarray(dtype=float) + pre-computed velocity Returns ------- """ - if len(velocity) == 0: + if velocity is None: is_running = compute_speed(pos, pos_tt) > speed_thresh else: - is_running = velocity > speed_thresh + is_running = np.linalg.norm(velocity) > speed_thresh spike_pos_inds = find_nearest(spikes, pos_tt) spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] @@ -137,7 +141,7 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0184, x_start=None, x_stop=None, y_start=None, y_stop=None, - velocity=[]): + velocity=None): """Returns speed-gated occupancy and speed-gated and Gaussian-filtered firing rate @@ -158,7 +162,8 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, x_stop: float, optional y_start: float, optional y_stop: float, optional - + velocity: np.ndarray(dtype=float) + pre-computed velocity Returns ------- @@ -195,9 +200,13 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, return occupancy, filtered_firing_rate, [edges_x, edges_y] -def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03): +def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03, velocity=None): + + if velocity is None: + is_running = compute_speed(pos, pos_tt) > speed_thresh + else: + is_running = np.linalg.norm(velocity) > speed_thresh - is_running = compute_speed(pos, pos_tt) > speed_thresh run_pos = pos[is_running, :] finite_lin_pos = run_pos[np.isfinite(run_pos)] @@ -208,7 +217,7 @@ def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh= def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, - spatial_bin_len=0.0168, speed_thresh=0.03): + spatial_bin_len=0.0168, speed_thresh=0.03, velocity=None): """The occupancy and number of spikes, speed-gated, binned, and smoothed over position @@ -224,7 +233,8 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, in meters. Default = 5.57 cm spatial_bin_len: float (optional) in meters. Default = 1.68 cm - + velocity: np.ndarray(dtype=float) + pre-computed velocity Returns ------- @@ -242,7 +252,7 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, sampling_rate = len(pos_tt) / (np.nanmax(pos_tt) - np.nanmin(pos_tt)) - occupancy = compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate) + occupancy = compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, velocity) is_running = compute_speed(pos, pos_tt) > speed_thresh diff --git a/nwbwidgets/view.py b/nwbwidgets/view.py index b0b36d5d..4e172be1 100644 --- a/nwbwidgets/view.py +++ b/nwbwidgets/view.py @@ -49,7 +49,8 @@ def show_dynamic_table(node, **kwargs) -> widgets.Widget: pynwb.behavior.SpatialSeries: OrderedDict({ 'over time': timeseries.SeparateTracesPlotlyWidget, 'trace': behavior.plotly_show_spatial_trace, - 'rate map': placefield.route_placefield}), + 'rate map': placefield.route_placefield, + '1D rate map': placefield.PlaceField_1D_Widget}), pynwb.image.GrayscaleImage: image.show_grayscale_image, pynwb.image.RGBImage: image.show_rbga_image, pynwb.image.RGBAImage: image.show_rbga_image, From 69bc054d0ba74703e4cadf3d26ea69298aff23ed Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 5 Nov 2020 20:27:28 -0500 Subject: [PATCH 24/39] Cleaned up code, fixed single unit input bug --- nwbwidgets/analysis/placefields.py | 9 ++- nwbwidgets/placefield.py | 112 ++++++++--------------------- 2 files changed, 38 insertions(+), 83 deletions(-) diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index dda1691c..5235d85e 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -82,6 +82,7 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc """ sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) + np.seterr(invalid='ignore') if velocity is None: is_running = compute_speed(pos, pos_tt) > speed_thresh else: @@ -118,7 +119,7 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 Returns ------- """ - + np.seterr(invalid='ignore') if velocity is None: is_running = compute_speed(pos, pos_tt) > speed_thresh else: @@ -188,6 +189,7 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh, velocity) + np.seterr(divide='ignore') firing_rate = n_spikes / occupancy # in Hz firing_rate[np.isnan(firing_rate)] = 0 # get rid of NaNs so convolution works @@ -202,6 +204,7 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03, velocity=None): + np.seterr(invalid='ignore') if velocity is None: is_running = compute_speed(pos, pos_tt) > speed_thresh else: @@ -252,8 +255,9 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, sampling_rate = len(pos_tt) / (np.nanmax(pos_tt) - np.nanmin(pos_tt)) - occupancy = compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, velocity) + occupancy = compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh, velocity) + np.seterr(invalid='ignore') is_running = compute_speed(pos, pos_tt) > speed_thresh # find pos_tt bin associated with each spike @@ -264,6 +268,7 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, n_spikes = np.histogram(finite_pos_on_spikes, bins=spatial_bins)[0][:-2] + np.seterr(divide='ignore') firing_rate = np.nan_to_num(n_spikes / occupancy) filtered_firing_rate = gaussian_filter( diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 13c0869f..7bd27d74 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -15,42 +15,7 @@ from .base import vis2widget from .controllers import GroupAndSortController - -## To-do -# [X] Create PlaceFieldWidget class - # [X] Refactor place field calculation code to deal with nwb data type - # [X] Incorporate place field fxns into class - # [X] Change all internal attributes references - # [X]Change all internal method references - - # [X] Get pos - # [X] Get time - # [X] Get spikes - # [] Get trials / epochs - -# [X] Submit draft PR - - # [] 1D Place Field Widget - # [X] Incorporate nelpy package into widget - # [X] Speed threshold implementation - # [X] Add foreign group and sort controller to pick unit groups and ranges? - # [X] Normalized firing rate figure? - # [X] Add collapsed unit vizualization? - # [] Scale bar? - # [] Sort place cell tuning curves by peak firing rate position? - # [] Color palette control? - -# [] Dropdown that controls which unit - -# [x] Work in buttons / dropdowns / sliders to modify following parameters in place field calculation: -# [] Different epochs -# [x] Gaussian SD -# [x] Speed threshold -# [] Minimum firing rate -# [] Place field thresh (% of local max) - def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): - print(spatial_series.data.shape) if spatial_series.data.shape[1] == 2: return PlaceFieldWidget(spatial_series) elif spatial_series.data.shape[1] == 1: @@ -59,33 +24,21 @@ def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): print('Spatial series exceeds dimensionality for visualization') return - -# Put widget rendering here class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb.TimeSeries = None, **kwargs): super().__init__() - - - if hasattr(spatial_series.get_ancestor('NWBFile'),'velocity'): - velocity = spatial_series.get_ancestor('NWBFile').velocity - else: - velocity = [] - self.units = spatial_series.get_ancestor('NWBFile').units self.pos_tt = get_timeseries_tt(spatial_series) - self.velocity = velocity + if velocity is not None: + self.velocity = velocity + else: + self.velocity = None istart = 0 istop = None - self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) - self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 - # Put widget controls here: - # - Minimum firing rate - # - Place field thresh (% of local max) - bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd (cm)') bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)') dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') @@ -138,16 +91,20 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_veloci class PlaceField_1D_Widget(widgets.HBox): - def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, foreign_group_and_sort_controller: GroupAndSortController = None, - group_by=None, **kwargs): + group_by=None, + velocity: pynwb.TimeSeries = None, + **kwargs): super().__init__() self.units = spatial_series.get_ancestor('NWBFile').units - self.pos_tt = get_timeseries_tt(spatial_series) + if velocity is not None: + self.velocity = velocity + else: + self.velocity = None if foreign_group_and_sort_controller: self.gas = foreign_group_and_sort_controller @@ -160,10 +117,6 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 - # Put widget controls here: - # - Minimum firing rate - # - Place field thresh (% of local max) - bft_gaussian = BoundedFloatText(value=0.0557, min=0, max=99999, description='gaussian sd (m)') bft_spatial_bin_len = BoundedFloatText(value=0.0168, min=0, max=99999, description='spatial bin length (m)') cb_normalize_select = Checkbox(value=False, description='normalize', indent=False) @@ -197,39 +150,37 @@ def do_1d_rate_map(self, units_window=None, order=None, group_inds=None, labels= tmin = min(self.pos_tt) tmax = max(self.pos_tt) - print('Order is {}'.format(order)) - print('Units_window is {}'.format(units_window)) - print('Group_inds is {}'.format(group_inds)) - print('Labels is {}'.format(labels)) if order is None: index = np.arange(0, len(self.units)) - else: + elif isinstance(order, np.ndarray): index = order + else: + index = [order] - spikes = get_spike_times(self.units, index[0], [tmin, tmax]) - xx, occupancy, filtered_firing_rate = compute_linear_firing_rate( - self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) - - all_unit_firing_rate = np.zeros([len(index), len(xx)]) - all_unit_firing_rate[0] = filtered_firing_rate firing_rate_ind = 0 - for ind in index[1:]: + for ind in index: + if firing_rate_ind == 0: + spikes = get_spike_times(self.units, ind, [tmin, tmax]) + xx, _, all_unit_firing_rate_temp = compute_linear_firing_rate(self.pos, self.pos_tt, spikes, + gaussian_sd=gaussian_sd, + spatial_bin_len=spatial_bin_len, + velocity=self.velocity) + all_unit_firing_rate = np.zeros([len(index), len(xx)]) + all_unit_firing_rate[0] = all_unit_firing_rate_temp + firing_rate_ind += 1 + continue spikes = get_spike_times(self.units, ind, [tmin, tmax]) - _, _, all_unit_firing_rate[firing_rate_ind] = compute_linear_firing_rate( - self.pos, self.pos_tt, spikes, gaussian_sd=gaussian_sd, spatial_bin_len=spatial_bin_len) + xx, _, all_unit_firing_rate[firing_rate_ind] = compute_linear_firing_rate(self.pos, self.pos_tt, spikes, + gaussian_sd=gaussian_sd, + spatial_bin_len=spatial_bin_len, + velocity=self.velocity) firing_rate_ind += 1 - # npl.set_palette(npl.colors.rainbow) - # with npl.FigureManager(show=True, figsize=(8, 8)) as (fig, ax): - # npl.utils.skip_if_no_output(fig) + fig, ax = plt.subplots() plot_tuning_curves1D(all_unit_firing_rate, xx, ax=ax, unit_labels=index, normalize=normalize, collapsed=collapsed) - # fig = ax.plot(xx, filtered_firing_rate, '-') - # ax.set_xlabel('x ({})'.format(self.unit)) - # ax.set_ylabel('firing rate (Hz)') - return fig def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, unit_labels=None, fill=True, color=None, @@ -273,8 +224,7 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni peak_firing_rates = ratemap.max(axis=1) ratemap = (ratemap.T / peak_firing_rates).T pad = 1 - # # determine max firing rate - # max_firing_rate = ratemap.max() + if collapsed: pad = 0 From 6fbd00654044cc9f8702953b9b66c6325912e08f Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Tue, 17 Nov 2020 23:40:51 -0500 Subject: [PATCH 25/39] Removed group and sort controller. Added multi-select. --- nwbwidgets/placefield.py | 59 ++++++++++++++++++++-------------------- 1 file changed, 30 insertions(+), 29 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 7bd27d74..07bdc1c8 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -5,7 +5,7 @@ import matplotlib.pyplot as plt import pynwb -from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox, fixed +from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox, Layout from .analysis.placefields import compute_2d_firing_rate, compute_linear_firing_rate @@ -39,8 +39,9 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 - bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd (cm)') - bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)') + style = {'description_width': 'initial'} + bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd (cm)', style=style) + bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)', style=style) dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') cb_velocity = Checkbox(value=False, description='use velocity', indent=False) @@ -100,62 +101,62 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, super().__init__() self.units = spatial_series.get_ancestor('NWBFile').units + index = np.arange(1, len(self.units)) + self.pos_tt = get_timeseries_tt(spatial_series) if velocity is not None: self.velocity = velocity else: self.velocity = None - if foreign_group_and_sort_controller: - self.gas = foreign_group_and_sort_controller - else: - self.gas = self.make_group_and_sort(group_by=group_by, control_order=False) - istart = 0 istop = None self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 - bft_gaussian = BoundedFloatText(value=0.0557, min=0, max=99999, description='gaussian sd (m)') - bft_spatial_bin_len = BoundedFloatText(value=0.0168, min=0, max=99999, description='spatial bin length (m)') + style = {'description_width': 'initial'} + bft_gaussian = BoundedFloatText(value=0.0557, min=0, max=99999, description='gaussian sd (m)', style=style) + bft_spatial_bin_len = BoundedFloatText(value=0.0168, min=0, max=99999, description='spatial bin length (m)', + style=style) cb_normalize_select = Checkbox(value=False, description='normalize', indent=False) cb_collapsed_select = Checkbox(value=False, description='collapsed', indent=False) + sm_unit_select = widgets.SelectMultiple(options=index, + value=[1, 2, 3, 4, 5], rows= 20, + description='Select units', disabled=False + ) self.controls = dict( - gas=self.gas, gaussian_sd=bft_gaussian, spatial_bin_len=bft_spatial_bin_len, normalize=cb_normalize_select, - collapsed=cb_collapsed_select + collapsed=cb_collapsed_select, + order=sm_unit_select ) out_fig = interactive_output(self.do_1d_rate_map, self.controls) checkboxes = widgets.HBox([cb_normalize_select, cb_collapsed_select]) - self.children = [ + widget_fig = vis2widget(out_fig) + self.children = [widgets.HBox([ widgets.VBox([ bft_gaussian, bft_spatial_bin_len, checkboxes, - self.gas, - ]), - vis2widget(out_fig) + sm_unit_select + ], + layout = Layout(max_width="40%")), + widget_fig], + layout = Layout(width="100%", height="100%")) ] def make_group_and_sort(self, group_by=None, control_order=True): return GroupAndSortController(self.units, group_by=group_by, control_order=control_order) - def do_1d_rate_map(self, units_window=None, order=None, group_inds=None, labels=None, normalize=False, - collapsed=False, gaussian_sd=0.0557, spatial_bin_len=0.0168, **kwargs): + def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian_sd=0.0557, + spatial_bin_len=0.0168, **kwargs): tmin = min(self.pos_tt) tmax = max(self.pos_tt) - - if order is None: - index = np.arange(0, len(self.units)) - elif isinstance(order, np.ndarray): - index = order - else: - index = [order] + index = np.asarray(order) firing_rate_ind = 0 for ind in index: @@ -256,11 +257,11 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni ax.spines['left'].set_color('none') ax.yaxis.set_ticks_position('right') else: - if normalize: - ax.set_ylabel('normalized firing rate') - else: - ax.set_ylabel('firing rate [Hz]') ax.set_ylim(0) + if normalize: + ax.set_ylabel('normalized firing rate') + else: + ax.set_ylabel('firing rate [Hz]') ax.spines['top'].set_color('none') ax.xaxis.set_ticks_position('bottom') From 9dab486bce0682a775f11080e6aa6d62fc6b7e95 Mon Sep 17 00:00:00 2001 From: "!git for-each-ref --format='%(refname:short)' `git symbolic-ref HEAD`" Date: Wed, 18 Nov 2020 16:31:53 -0500 Subject: [PATCH 26/39] * lru_cache * figure size --- nwbwidgets/placefield.py | 47 ++++++++++++++++++++-------------------- 1 file changed, 24 insertions(+), 23 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 07bdc1c8..ffcd8391 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -1,6 +1,7 @@ import numpy as np from scipy.ndimage import label from scipy.ndimage.filters import gaussian_filter, maximum_filter +from functools import lru_cache import matplotlib.pyplot as plt @@ -13,7 +14,8 @@ from .utils.units import get_spike_times from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt from .base import vis2widget -from .controllers import GroupAndSortController +from .controllers import GroupAndSortController + def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): if spatial_series.data.shape[1] == 2: @@ -144,9 +146,9 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, checkboxes, sm_unit_select ], - layout = Layout(max_width="40%")), + layout=Layout(max_width="40%")), widget_fig], - layout = Layout(width="100%", height="100%")) + layout=Layout(width="100%", height="100%")) ] def make_group_and_sort(self, group_by=None, control_order=True): @@ -158,32 +160,31 @@ def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian tmax = max(self.pos_tt) index = np.asarray(order) - firing_rate_ind = 0 - for ind in index: - if firing_rate_ind == 0: - spikes = get_spike_times(self.units, ind, [tmin, tmax]) - xx, _, all_unit_firing_rate_temp = compute_linear_firing_rate(self.pos, self.pos_tt, spikes, - gaussian_sd=gaussian_sd, - spatial_bin_len=spatial_bin_len, - velocity=self.velocity) + for i, ind in enumerate(index): + + all_unit_firing_rate_temp, xx = self.compute_1d_firing_rate( + ind, tmin, tmax, gaussian_sd, spatial_bin_len) + if not i: all_unit_firing_rate = np.zeros([len(index), len(xx)]) - all_unit_firing_rate[0] = all_unit_firing_rate_temp - firing_rate_ind += 1 - continue - spikes = get_spike_times(self.units, ind, [tmin, tmax]) - xx, _, all_unit_firing_rate[firing_rate_ind] = compute_linear_firing_rate(self.pos, self.pos_tt, spikes, - gaussian_sd=gaussian_sd, - spatial_bin_len=spatial_bin_len, - velocity=self.velocity) - firing_rate_ind += 1 + all_unit_firing_rate[ind] = all_unit_firing_rate_temp - fig, ax = plt.subplots() + fig, ax = plt.subplots(figsize=(7, 7)) plot_tuning_curves1D(all_unit_firing_rate, xx, ax=ax, unit_labels=index, normalize=normalize, collapsed=collapsed) return fig + @lru_cache + def compute_1d_firing_rate(self, ind, tmin, tmax, gaussian_sd, spatial_bin_len): + spikes = get_spike_times(self.units, ind, [tmin, tmax]) + xx, _, all_unit_firing_rate_temp = compute_linear_firing_rate(self.pos, self.pos_tt, spikes, + gaussian_sd=gaussian_sd, + spatial_bin_len=spatial_bin_len, + velocity=self.velocity) + return all_unit_firing_rate_temp, xx + + def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, unit_labels=None, fill=True, color=None, collapsed=False): """ @@ -194,7 +195,7 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni An array of dim: [number of units, bin positions] with the spike rates for a unit, at every pos, in each row bin_pos: array-like An array representing the bin positions of ratemap for each column - ax: matplotlib.pyplot.ax + ax: matplotlib.pyplot.Axes Axes object for the figure on which the ratemaps will be plotted normalize: bool default = False @@ -212,7 +213,7 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni Returns ------- - matplotlib.pyplot.ax + matplotlib.pyplot.Axes """ xmin = bin_pos[0] From 52d17062ef1384b1b86918ce2ab4fa52796be987 Mon Sep 17 00:00:00 2001 From: "!git for-each-ref --format='%(refname:short)' `git symbolic-ref HEAD`" Date: Sat, 21 Nov 2020 16:21:05 -0500 Subject: [PATCH 27/39] correct index and cache --- nwbwidgets/placefield.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index ffcd8391..490b29d3 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -167,7 +167,7 @@ def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian if not i: all_unit_firing_rate = np.zeros([len(index), len(xx)]) - all_unit_firing_rate[ind] = all_unit_firing_rate_temp + all_unit_firing_rate[i] = all_unit_firing_rate_temp fig, ax = plt.subplots(figsize=(7, 7)) plot_tuning_curves1D(all_unit_firing_rate, xx, ax=ax, unit_labels=index, normalize=normalize, @@ -175,7 +175,7 @@ def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian return fig - @lru_cache + @lru_cache() def compute_1d_firing_rate(self, ind, tmin, tmax, gaussian_sd, spatial_bin_len): spikes = get_spike_times(self.units, ind, [tmin, tmax]) xx, _, all_unit_firing_rate_temp = compute_linear_firing_rate(self.pos, self.pos_tt, spikes, From 84e87aab0eda95a1d9cbdd84886b30eb3a05b5fd Mon Sep 17 00:00:00 2001 From: "!git for-each-ref --format='%(refname:short)' `git symbolic-ref HEAD`" Date: Tue, 24 Nov 2020 15:24:43 -0500 Subject: [PATCH 28/39] style and caching some analyses --- nwbwidgets/analysis/placefields.py | 27 ++++++++++++++------ nwbwidgets/ecephys.py | 8 +++--- nwbwidgets/placefield.py | 40 ++++++++++++++---------------- nwbwidgets/view.py | 2 +- 4 files changed, 44 insertions(+), 33 deletions(-) diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index 5235d85e..d1d2d82c 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -1,16 +1,20 @@ +from functools import lru_cache + import numpy as np -from scipy.ndimage.filters import gaussian_filter, maximum_filter +from scipy.ndimage.filters import gaussian_filter + def find_nearest(arr, tt): """Used for picking out elements of a TimeSeries based on spike times Parameters ---------- - arr - tt + arr: array-like + tt: array-like Returns ------- + indices: array-like """ arr = arr[arr > tt[0]] @@ -23,17 +27,19 @@ def smooth(y, box_pts): Parameters ---------- - y - box_pts + y: array-like + box_pts: int Returns ------- + output: np.array(dtype=float) """ box = np.ones(box_pts) / box_pts return np.convolve(y, box, mode='same') +@lru_cache() def compute_speed(pos, pos_tt, smooth_param=40): """Compute boolean of whether the speed of the animal was above a threshold for each time point @@ -55,6 +61,7 @@ def compute_speed(pos, pos_tt, smooth_param=40): return smooth(speed, smooth_param) +@lru_cache() def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, velocity=None): """Computes occupancy per bin in seconds @@ -86,7 +93,7 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc if velocity is None: is_running = compute_speed(pos, pos_tt) > speed_thresh else: - is_running = np.linalg.norm(self.velocity) > speed_thresh + is_running = np.linalg.norm(velocity) > speed_thresh run_pos = pos[is_running, :] occupancy = np.histogram2d(run_pos[:, 0], @@ -96,6 +103,7 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc return occupancy, is_running +@lru_cache() def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03, velocity=None): """Returns speed-gated position during spikes @@ -136,6 +144,7 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 return n_spikes +@lru_cache() def compute_2d_firing_rate(pos, pos_tt, spikes, pixel_width, speed_thresh=0.03, @@ -202,6 +211,7 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, return occupancy, filtered_firing_rate, [edges_x, edges_y] +@lru_cache() def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03, velocity=None): np.seterr(invalid='ignore') @@ -219,6 +229,7 @@ def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh= return occupancy +@lru_cache() def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, spatial_bin_len=0.0168, speed_thresh=0.03, velocity=None): """The occupancy and number of spikes, speed-gated, binned, and smoothed @@ -236,6 +247,8 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, in meters. Default = 5.57 cm spatial_bin_len: float (optional) in meters. Default = 1.68 cm + speed_thresh: float (optional) + in m/s. Default = 0.03 velocity: np.ndarray(dtype=float) pre-computed velocity @@ -275,4 +288,4 @@ def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, firing_rate, gaussian_sd / spatial_bin_len) xx = spatial_bins[:-3] + (spatial_bins[1] - spatial_bins[0]) / 2 - return xx, occupancy, filtered_firing_rate \ No newline at end of file + return xx, occupancy, filtered_firing_rate diff --git a/nwbwidgets/ecephys.py b/nwbwidgets/ecephys.py index dd0cfc34..b633e5b5 100644 --- a/nwbwidgets/ecephys.py +++ b/nwbwidgets/ecephys.py @@ -1,11 +1,11 @@ import matplotlib.pyplot as plt import numpy as np import plotly.graph_objects as go -from plotly.colors import DEFAULT_PLOTLY_COLORS +import pynwb from ipywidgets import widgets, ValueWidget +from plotly.colors import DEFAULT_PLOTLY_COLORS from pynwb.ecephys import LFP, SpikeEventSeries, ElectricalSeries from scipy.signal import stft -import pynwb from .base import fig2widget, nwb2widget, lazy_tabs, render_dataframe from .timeseries import BaseGroupedTraceWidget @@ -87,8 +87,8 @@ def show_electrodes(electrodes_table): if np.isnan(electrodes_table.x[0]): # position is not defined in_dict.update(electrode_groups=ElectrodeGroupsWidget) else: - if electrodes_table.get_ancestor('NWBFile').subject.species \ - in ('mouse', 'Mus musculus'): + subject = electrodes_table.get_ancestor('NWBFile').subject + if subject is not None and subject.species in ('mouse', 'Mus musculus'): in_dict.update(CCF=show_ccf) return lazy_tabs(in_dict, electrodes_table) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 490b29d3..12eb08e7 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -1,31 +1,28 @@ -import numpy as np -from scipy.ndimage import label -from scipy.ndimage.filters import gaussian_filter, maximum_filter from functools import lru_cache import matplotlib.pyplot as plt - +import numpy as np import pynwb from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox, Layout from .analysis.placefields import compute_2d_firing_rate, compute_linear_firing_rate - -from .utils.widgets import interactive_output -from .utils.units import get_spike_times -from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt from .base import vis2widget from .controllers import GroupAndSortController +from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt +from .utils.units import get_spike_times +from .utils.widgets import interactive_output def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): if spatial_series.data.shape[1] == 2: return PlaceFieldWidget(spatial_series) elif spatial_series.data.shape[1] == 1: - return PlaceField_1D_Widget(spatial_series) + return PlaceField1DWidget(spatial_series) else: print('Spatial series exceeds dimensionality for visualization') return + class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb.TimeSeries = None, **kwargs): @@ -71,9 +68,10 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_veloci tmax = max(self.pos_tt) spikes = get_spike_times(self.units, index, [tmin, tmax]) - if use_velocity == False: + if not use_velocity: occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( - self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, velocity=self.velocity) + self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, + velocity=self.velocity) else: occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, @@ -93,7 +91,7 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_veloci return fig -class PlaceField_1D_Widget(widgets.HBox): +class PlaceField1DWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, foreign_group_and_sort_controller: GroupAndSortController = None, group_by=None, @@ -111,9 +109,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, else: self.velocity = None - istart = 0 - istop = None - self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) + self.pos, self.unit = get_timeseries_in_units(spatial_series) self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 @@ -124,7 +120,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, cb_normalize_select = Checkbox(value=False, description='normalize', indent=False) cb_collapsed_select = Checkbox(value=False, description='collapsed', indent=False) sm_unit_select = widgets.SelectMultiple(options=index, - value=[1, 2, 3, 4, 5], rows= 20, + value=[1, 2, 3, 4, 5], rows=20, description='Select units', disabled=False ) @@ -154,7 +150,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, def make_group_and_sort(self, group_by=None, control_order=True): return GroupAndSortController(self.units, group_by=group_by, control_order=control_order) - def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian_sd=0.0557, + def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian_sd=0.0557, spatial_bin_len=0.0168, **kwargs): tmin = min(self.pos_tt) tmax = max(self.pos_tt) @@ -206,6 +202,7 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni Amount of space to put between each unit (i.e. row) in the figure unit_labels: array-like Unit ids for each unit in ratemap + color: matplotlib color collapsed: bool default = False Determines whether to plot the ratemaps with zero padding, i.e. at the same y coordinate, on the ratemap @@ -239,17 +236,18 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni for unit, curve in enumerate(ratemap): if color is None: - line = ax.plot(xvals, unit*pad + curve, zorder=int(10+2*n_units-2*unit)) + line = ax.plot(xvals, unit * pad + curve, zorder=int(10 + 2 * n_units - 2 * unit)) else: - line = ax.plot(xvals, unit*pad + curve, zorder=int(10+2*n_units-2*unit), color=color) + line = ax.plot(xvals, unit * pad + curve, zorder=int(10 + 2 * n_units - 2 * unit), color=color) if fill: # Get the color from the current curve fillcolor = line[0].get_color() - ax.fill_between(xvals, unit*pad, unit*pad + curve, alpha=0.3, color=fillcolor, zorder=int(10+2*n_units-2*unit-1)) + ax.fill_between(xvals, unit * pad, unit * pad + curve, alpha=0.3, color=fillcolor, + zorder=int(10 + 2 * n_units - 2 * unit - 1)) ax.set_xlim(xmin, xmax) if pad != 0: - yticks = np.arange(n_units)*pad + 0.5*pad + yticks = np.arange(n_units) * pad + 0.5 * pad ax.set_yticks(yticks) ax.set_yticklabels(unit_labels) ax.set_xlabel('external variable') diff --git a/nwbwidgets/view.py b/nwbwidgets/view.py index 1b05d00c..a35a1317 100644 --- a/nwbwidgets/view.py +++ b/nwbwidgets/view.py @@ -50,7 +50,7 @@ def show_dynamic_table(node, **kwargs) -> widgets.Widget: 'over time': timeseries.SeparateTracesPlotlyWidget, 'trace': behavior.plotly_show_spatial_trace, 'rate map': placefield.route_placefield, - '1D rate map': placefield.PlaceField_1D_Widget}), + '1D rate map': placefield.PlaceField1DWidget}), pynwb.image.GrayscaleImage: image.show_grayscale_image, pynwb.image.RGBImage: image.show_rbga_image, pynwb.image.RGBAImage: image.show_rbga_image, From 6bb219e36ec95717108749eb9aea632c9ab2d4bd Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Wed, 9 Dec 2020 13:22:19 -0500 Subject: [PATCH 29/39] Place field modification for towers task. --- nwbwidgets/placefield.py | 44 +++++++++++++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 10 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index ffcd8391..3ca417df 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -28,7 +28,10 @@ def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): class PlaceFieldWidget(widgets.HBox): - def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb.TimeSeries = None, **kwargs): + def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, + velocity: pynwb.TimeSeries = None, + towers: pynwb.epoch.TimeIntervals = None, + **kwargs): super().__init__() self.units = spatial_series.get_ancestor('NWBFile').units self.pos_tt = get_timeseries_tt(spatial_series) @@ -36,9 +39,26 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb self.velocity = velocity else: self.velocity = None + istart = 0 istop = None self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) + if towers is not None: + left_towers = towers.left cue_onset + right_towers = towers.right_cue_onset + tt = self.pos_tt[:, 0] + ss = np.zeros_like(tt) + ss[np.searchsorted(tt, right_towers)] += 1 + ss[np.searchsorted(tt, left_towers)] -= 1 + starts = np.searchsorted(tt, towers.start_time) + ends = np.searchsorted(tt, towers.stop_time) + ss = np.zeros_like(tt) + for start, end in zip(starts,ends): + ss[start:end] = np.cumsum(ss[start:end]) + + self.pos[:, 0] = self.pos[:, 1] + self.pos[:, 1] = states + self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 style = {'description_width': 'initial'} @@ -51,7 +71,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb gaussian_sd=bft_gaussian, speed_thresh=bft_speed, index=dd_unit_select, - use_velocity=cb_velocity, + use_velocity=cb_velocity ) out_fig = interactive_output(self.do_rate_map, self.controls) @@ -61,7 +81,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb bft_gaussian, bft_speed, dd_unit_select, - cb_velocity + cb_velocity, ]), vis2widget(out_fig) ] @@ -72,12 +92,16 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_veloci spikes = get_spike_times(self.units, index, [tmin, tmax]) if use_velocity == False: - occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( - self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, velocity=self.velocity) + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(self.pos, self.pos_tt, spikes, + self.pixel_width, + speed_thresh=speed_thresh, + gaussian_sd=gaussian_sd) else: - occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate( - self.pos, self.pos_tt, spikes, self.pixel_width, speed_thresh=speed_thresh, gaussian_sd=gaussian_sd, - velocity=self.velocity) + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(self.pos, self.pos_tt, spikes, + self.pixel_width, + speed_thresh=speed_thresh, + gaussian_sd=gaussian_sd, + velocity=self.velocity) fig, ax = plt.subplots() @@ -167,7 +191,7 @@ def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian if not i: all_unit_firing_rate = np.zeros([len(index), len(xx)]) - all_unit_firing_rate[ind] = all_unit_firing_rate_temp + all_unit_firing_rate[i] = all_unit_firing_rate_temp fig, ax = plt.subplots(figsize=(7, 7)) plot_tuning_curves1D(all_unit_firing_rate, xx, ax=ax, unit_labels=index, normalize=normalize, @@ -175,7 +199,7 @@ def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian return fig - @lru_cache + @lru_cache() def compute_1d_firing_rate(self, ind, tmin, tmax, gaussian_sd, spatial_bin_len): spikes = get_spike_times(self.units, ind, [tmin, tmax]) xx, _, all_unit_firing_rate_temp = compute_linear_firing_rate(self.pos, self.pos_tt, spikes, From 43080e8f0d4a03fb4d277e998025a51e62b99091 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Wed, 9 Dec 2020 13:34:37 -0500 Subject: [PATCH 30/39] Place field modification for towers task. --- nwbwidgets/placefield.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 3ca417df..b9675f58 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -1,21 +1,17 @@ -import numpy as np -from scipy.ndimage import label -from scipy.ndimage.filters import gaussian_filter, maximum_filter from functools import lru_cache import matplotlib.pyplot as plt - +import numpy as np import pynwb from ipywidgets import widgets, BoundedFloatText, Dropdown, Checkbox, Layout from .analysis.placefields import compute_2d_firing_rate, compute_linear_firing_rate +from .base import vis2widget +from .controllers import GroupAndSortController from .utils.widgets import interactive_output from .utils.units import get_spike_times from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt -from .base import vis2widget -from .controllers import GroupAndSortController - def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): if spatial_series.data.shape[1] == 2: From c365f2b25bc2e7e5da169e4ed798aaee25fb3706 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Wed, 9 Dec 2020 23:39:50 -0500 Subject: [PATCH 31/39] Fixed the "TypeError: unhashable type: 'numpy.ndarray" issue that arose from caching the calls in placefields.py with numpy arrays im the cache keys. Incorporated an additional input in the widget to specify the width of the gaussian kernel independently for x and y respectively. Built a wrapper function around the 2D place field calculation call so that the place field calculation function calls could be properly cached without the "unhashable key" errors. --- nwbwidgets/analysis/placefields.py | 22 +++---- nwbwidgets/placefield.py | 99 +++++++++++++----------------- 2 files changed, 53 insertions(+), 68 deletions(-) diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index d1d2d82c..b0cb078a 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -39,7 +39,6 @@ def smooth(y, box_pts): return np.convolve(y, box, mode='same') -@lru_cache() def compute_speed(pos, pos_tt, smooth_param=40): """Compute boolean of whether the speed of the animal was above a threshold for each time point @@ -61,7 +60,6 @@ def compute_speed(pos, pos_tt, smooth_param=40): return smooth(speed, smooth_param) -@lru_cache() def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, velocity=None): """Computes occupancy per bin in seconds @@ -103,7 +101,6 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc return occupancy, is_running -@lru_cache() def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03, velocity=None): """Returns speed-gated position during spikes @@ -144,11 +141,11 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 return n_spikes -@lru_cache() def compute_2d_firing_rate(pos, pos_tt, spikes, pixel_width, speed_thresh=0.03, - gaussian_sd=0.0184, + gaussian_sd_x=0.0184, + gaussian_sd_y=0.0184, x_start=None, x_stop=None, y_start=None, y_stop=None, velocity=None): @@ -166,8 +163,10 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, pixel_width: float speed_thresh: float, optional in meters. Default = 3.0 cm/s - gaussian_sd: float, optional - in meters. Default = 1.84 cm + gaussian_sd_x: float, optional + width of gaussian kernel in x-dim, in meters. Default = 1.84 cm + gaussian_sd_y: float, optional + width of gaussian kernel in y-dim, in meters. Default = 1.84 cm x_start: float, optional x_stop: float, optional y_start: float, optional @@ -201,17 +200,17 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, np.seterr(divide='ignore') firing_rate = n_spikes / occupancy # in Hz firing_rate[np.isnan(firing_rate)] = 0 # get rid of NaNs so convolution works - - filtered_firing_rate = gaussian_filter(firing_rate, gaussian_sd / pixel_width) + sigmas = [gaussian_sd_y / pixel_width, gaussian_sd_x / pixel_width] + filtered_firing_rate = gaussian_filter(firing_rate, sigmas) # filter occupancy to create a mask so non-explored regions are nan'ed - filtered_occupancy = gaussian_filter(occupancy, gaussian_sd / pixel_width / 8) + sigmas_occ = [gaussian_sd_y / pixel_width / 8, gaussian_sd_x / pixel_width / 8] + filtered_occupancy = gaussian_filter(occupancy, sigmas_occ) filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan return occupancy, filtered_firing_rate, [edges_x, edges_y] -@lru_cache() def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh=0.03, velocity=None): np.seterr(invalid='ignore') @@ -229,7 +228,6 @@ def compute_1d_occupancy(pos, pos_tt, spatial_bins, sampling_rate, speed_thresh= return occupancy -@lru_cache() def compute_linear_firing_rate(pos, pos_tt, spikes, gaussian_sd=0.0557, spatial_bin_len=0.0168, speed_thresh=0.03, velocity=None): """The occupancy and number of spikes, speed-gated, binned, and smoothed diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index b9675f58..0e6baa4d 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -17,7 +17,7 @@ def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): if spatial_series.data.shape[1] == 2: return PlaceFieldWidget(spatial_series) elif spatial_series.data.shape[1] == 1: - return PlaceField_1D_Widget(spatial_series) + return PlaceField1DWidget(spatial_series) else: print('Spatial series exceeds dimensionality for visualization') return @@ -26,7 +26,6 @@ class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb.TimeSeries = None, - towers: pynwb.epoch.TimeIntervals = None, **kwargs): super().__init__() self.units = spatial_series.get_ancestor('NWBFile').units @@ -36,35 +35,20 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, else: self.velocity = None - istart = 0 - istop = None - self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) - if towers is not None: - left_towers = towers.left cue_onset - right_towers = towers.right_cue_onset - tt = self.pos_tt[:, 0] - ss = np.zeros_like(tt) - ss[np.searchsorted(tt, right_towers)] += 1 - ss[np.searchsorted(tt, left_towers)] -= 1 - starts = np.searchsorted(tt, towers.start_time) - ends = np.searchsorted(tt, towers.stop_time) - ss = np.zeros_like(tt) - for start, end in zip(starts,ends): - ss[start:end] = np.cumsum(ss[start:end]) - - self.pos[:, 0] = self.pos[:, 1] - self.pos[:, 1] = states + self.pos, self.unit = get_timeseries_in_units(spatial_series) self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 style = {'description_width': 'initial'} - bft_gaussian = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd (cm)', style=style) + bft_gaussian_x = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd x (cm)', style=style) + bft_gaussian_y = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd y (cm)', style=style) bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)', style=style) dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') cb_velocity = Checkbox(value=False, description='use velocity', indent=False) self.controls = dict( - gaussian_sd=bft_gaussian, + gaussian_sd_x=bft_gaussian_x, + gaussian_sd_y=bft_gaussian_y, speed_thresh=bft_speed, index=dd_unit_select, use_velocity=cb_velocity @@ -74,7 +58,8 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, self.children = [ widgets.VBox([ - bft_gaussian, + bft_gaussian_x, + bft_gaussian_y, bft_speed, dd_unit_select, cb_velocity, @@ -82,23 +67,12 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, vis2widget(out_fig) ] - def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_velocity=False): - tmin = min(self.pos_tt) - tmax = max(self.pos_tt) - - spikes = get_spike_times(self.units, index, [tmin, tmax]) - if use_velocity == False: - occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(self.pos, self.pos_tt, spikes, - self.pixel_width, + def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian_sd_y=0.0184, use_velocity=False): + occupancy, filtered_firing_rate, [edges_x, edges_y] = self.compute_twodim_firing_rate(index=index, speed_thresh=speed_thresh, - gaussian_sd=gaussian_sd) - else: - occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(self.pos, self.pos_tt, spikes, - self.pixel_width, - speed_thresh=speed_thresh, - gaussian_sd=gaussian_sd, - velocity=self.velocity) - + gaussian_sd_x=gaussian_sd_x, + gaussian_sd_y=gaussian_sd_y, + use_velocity=use_velocity) fig, ax = plt.subplots() im = ax.imshow(filtered_firing_rate, @@ -112,11 +86,29 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd=0.0184, use_veloci return fig + @lru_cache() + def compute_twodim_firing_rate(self, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian_sd_y=0.0184, + use_velocity=False): + tmin = min(self.pos_tt) + tmax = max(self.pos_tt) + spikes = get_spike_times(self.units, index, [tmin, tmax]) + if use_velocity == False: + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(self.pos, self.pos_tt, spikes, + self.pixel_width, + speed_thresh=speed_thresh, + gaussian_sd_x=gaussian_sd_x, + gaussian_sd_y=gaussian_sd_y) + else: + occupancy, filtered_firing_rate, [edges_x, edges_y] = compute_2d_firing_rate(self.pos, self.pos_tt, spikes, + self.pixel_width, + speed_thresh=speed_thresh, + gaussian_sd_x=gaussian_sd_x, + gaussian_sd_y=gaussian_sd_y, + velocity=self.velocity) + return occupancy, filtered_firing_rate, [edges_x, edges_y] -class PlaceField_1D_Widget(widgets.HBox): +class PlaceField1DWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, - foreign_group_and_sort_controller: GroupAndSortController = None, - group_by=None, velocity: pynwb.TimeSeries = None, **kwargs): @@ -131,9 +123,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, else: self.velocity = None - istart = 0 - istop = None - self.pos, self.unit = get_timeseries_in_units(spatial_series, istart, istop) + self.pos, self.unit = get_timeseries_in_units(spatial_series) self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 @@ -144,7 +134,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, cb_normalize_select = Checkbox(value=False, description='normalize', indent=False) cb_collapsed_select = Checkbox(value=False, description='collapsed', indent=False) sm_unit_select = widgets.SelectMultiple(options=index, - value=[1, 2, 3, 4, 5], rows= 20, + value=[1, 2, 3, 4, 5], rows=20, description='Select units', disabled=False ) @@ -171,10 +161,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, layout=Layout(width="100%", height="100%")) ] - def make_group_and_sort(self, group_by=None, control_order=True): - return GroupAndSortController(self.units, group_by=group_by, control_order=control_order) - - def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian_sd=0.0557, + def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian_sd=0.0557, spatial_bin_len=0.0168, **kwargs): tmin = min(self.pos_tt) tmax = max(self.pos_tt) @@ -182,8 +169,7 @@ def do_1d_rate_map(self, order=None, normalize=False, collapsed=False, gaussian for i, ind in enumerate(index): - all_unit_firing_rate_temp, xx = self.compute_1d_firing_rate( - ind, tmin, tmax, gaussian_sd, spatial_bin_len) + all_unit_firing_rate_temp, xx = self.compute_1d_firing_rate(ind, tmin, tmax, gaussian_sd, spatial_bin_len) if not i: all_unit_firing_rate = np.zeros([len(index), len(xx)]) @@ -259,17 +245,18 @@ def plot_tuning_curves1D(ratemap, bin_pos, ax=None, normalize=False, pad=10, uni for unit, curve in enumerate(ratemap): if color is None: - line = ax.plot(xvals, unit*pad + curve, zorder=int(10+2*n_units-2*unit)) + line = ax.plot(xvals, unit * pad + curve, zorder=int(10 + 2 * n_units - 2 * unit)) else: - line = ax.plot(xvals, unit*pad + curve, zorder=int(10+2*n_units-2*unit), color=color) + line = ax.plot(xvals, unit * pad + curve, zorder=int(10 + 2 * n_units - 2 * unit), color=color) if fill: # Get the color from the current curve fillcolor = line[0].get_color() - ax.fill_between(xvals, unit*pad, unit*pad + curve, alpha=0.3, color=fillcolor, zorder=int(10+2*n_units-2*unit-1)) + ax.fill_between(xvals, unit * pad, unit * pad + curve, alpha=0.3, color=fillcolor, + zorder=int(10 + 2 * n_units - 2 * unit - 1)) ax.set_xlim(xmin, xmax) if pad != 0: - yticks = np.arange(n_units)*pad + 0.5*pad + yticks = np.arange(n_units) * pad + 0.5 * pad ax.set_yticks(yticks) ax.set_yticklabels(unit_labels) ax.set_xlabel('external variable') From 194f457e5262422482e452cdb9780b638eb60a96 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 10 Dec 2020 02:08:57 -0500 Subject: [PATCH 32/39] Remove group and sort controller import --- nwbwidgets/placefield.py | 1 - 1 file changed, 1 deletion(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 0e6baa4d..376db54b 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -8,7 +8,6 @@ from .analysis.placefields import compute_2d_firing_rate, compute_linear_firing_rate from .base import vis2widget -from .controllers import GroupAndSortController from .utils.widgets import interactive_output from .utils.units import get_spike_times from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt From 2624add3dce52f8d63fb051a06a0e23020938ee7 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Fri, 11 Dec 2020 21:13:10 -0500 Subject: [PATCH 33/39] Re-factor code Re-factor position and widget controls to make the PlaceFieldWidget class more easily extensible for the towers task place field widget --- nwbwidgets/placefield.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 376db54b..e6770e53 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -34,16 +34,11 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, else: self.velocity = None - self.pos, self.unit = get_timeseries_in_units(spatial_series) + self.get_position(spatial_series) self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 - style = {'description_width': 'initial'} - bft_gaussian_x = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd x (cm)', style=style) - bft_gaussian_y = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd y (cm)', style=style) - bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)', style=style) - dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') - cb_velocity = Checkbox(value=False, description='use velocity', indent=False) + bft_gaussian_x, bft_gaussian_y, bft_speed, dd_unit_select, cb_velocity = self.get_controls() self.controls = dict( gaussian_sd_x=bft_gaussian_x, @@ -66,6 +61,19 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, vis2widget(out_fig) ] + def get_position(self, spatial_series): + self.pos, self.unit = get_timeseries_in_units(spatial_series) + + def get_controls(self): + style = {'description_width': 'initial'} + bft_gaussian_x = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd x (cm)', style=style) + bft_gaussian_y = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd y (cm)', style=style) + bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)', style=style) + dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') + cb_velocity = Checkbox(value=False, description='use velocity', indent=False) + + return bft_gaussian_x, bft_gaussian_y, bft_speed, dd_unit_select, cb_velocity + def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian_sd_y=0.0184, use_velocity=False): occupancy, filtered_firing_rate, [edges_x, edges_y] = self.compute_twodim_firing_rate(index=index, speed_thresh=speed_thresh, From d6905bcba4ec6033028792756b740d42ff494e9c Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Wed, 16 Dec 2020 18:23:33 -0500 Subject: [PATCH 34/39] Allow for different pixel width in x and y dimension Refactored placefield to take separate pixel_widths. for x and y dimension. Modified placefields to reflect this change. --- .gitignore | 3 +++ nwbwidgets/analysis/placefields.py | 10 +++++----- nwbwidgets/placefield.py | 12 +++++++++--- 3 files changed, 17 insertions(+), 8 deletions(-) diff --git a/.gitignore b/.gitignore index fba6f61b..cfcd4378 100644 --- a/.gitignore +++ b/.gitignore @@ -129,3 +129,6 @@ dmypy.json # nwb **.nwb +nwbwidgets/controllers/group_and_sort_controllers_multi_select.py +nwbwidgets/placefield_group_and_sort.py +nwbwidgets/towers_task_placefield.py diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index b0cb078a..dc01dd33 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -190,8 +190,8 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, y_start = np.nanmin(pos[:, 1]) if y_start is None else y_start y_stop = np.nanmax(pos[:, 1]) if y_stop is None else y_stop - edges_x = np.arange(x_start, x_stop, pixel_width) - edges_y = np.arange(y_start, y_stop, pixel_width) + edges_x = np.arange(x_start, x_stop, pixel_width[1]) + edges_y = np.arange(y_start, y_stop, pixel_width[0]) occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh, velocity) @@ -200,13 +200,13 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, np.seterr(divide='ignore') firing_rate = n_spikes / occupancy # in Hz firing_rate[np.isnan(firing_rate)] = 0 # get rid of NaNs so convolution works - sigmas = [gaussian_sd_y / pixel_width, gaussian_sd_x / pixel_width] + sigmas = [gaussian_sd_y / pixel_width[1], gaussian_sd_x / pixel_width[0]] filtered_firing_rate = gaussian_filter(firing_rate, sigmas) # filter occupancy to create a mask so non-explored regions are nan'ed - sigmas_occ = [gaussian_sd_y / pixel_width / 8, gaussian_sd_x / pixel_width / 8] + sigmas_occ = [gaussian_sd_y / pixel_width[1] / 8, gaussian_sd_x / pixel_width[0] / 8] filtered_occupancy = gaussian_filter(occupancy, sigmas_occ) - filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan + # filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan return occupancy, filtered_firing_rate, [edges_x, edges_y] diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index e6770e53..1d9dae41 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -24,10 +24,13 @@ def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, - velocity: pynwb.TimeSeries = None, + velocity: pynwb.TimeSeries = None, units = None, **kwargs): super().__init__() - self.units = spatial_series.get_ancestor('NWBFile').units + if units is None: + self.units = spatial_series.get_ancestor('NWBFile').units + else: + self.units = units self.pos_tt = get_timeseries_tt(spatial_series) if velocity is not None: self.velocity = velocity @@ -36,7 +39,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, self.get_position(spatial_series) - self.pixel_width = (np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000 + self.get_pixel_width() bft_gaussian_x, bft_gaussian_y, bft_speed, dd_unit_select, cb_velocity = self.get_controls() @@ -61,6 +64,9 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, vis2widget(out_fig) ] + def get_pixel_width(self): + self.pixel_width = [(np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000] * 2 + def get_position(self, spatial_series): self.pos, self.unit = get_timeseries_in_units(spatial_series) From 04cca7172086898d0a59a2041c86fbe8b9c649c2 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Thu, 17 Dec 2020 00:41:09 -0500 Subject: [PATCH 35/39] Changes to accomodate towers task place field Modifications: - Disabled velocity button when not in use - Modified both compute_2d_occupancy and compute_2d_n_spikes to recognize when the towers task place field widget is calling it so that speed is computed on just the x-dim. - Modified compute_speed to work with 1-dim input - Fixed the order of inputs to np.histogram2d. Previously the x and y input were reversed from the order expected by the function. - Refactored the pixel_width instance atrribute definition to be separate for x and y. This allows the towers task place field widget extension to overwrite this method. --- nwbwidgets/analysis/placefields.py | 43 +++++++++++++++++++----------- nwbwidgets/placefield.py | 6 +++-- 2 files changed, 31 insertions(+), 18 deletions(-) diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index dc01dd33..ce585b99 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -56,11 +56,14 @@ def compute_speed(pos, pos_tt, smooth_param=40): running: np.ndarray(dtype=bool) """ - speed = np.hstack((0, np.sqrt(np.sum(np.diff(pos.T) ** 2, axis=0)) / np.diff(pos_tt))) + if len(pos.shape) > 1: + speed = np.hstack((0, np.sqrt(np.sum(np.diff(pos.T) ** 2, axis=0)) / np.diff(pos_tt))) + else: + speed = np.hstack((0, np.sqrt(np.diff(pos.T) ** 2) / np.diff(pos_tt))) return smooth(speed, smooth_param) -def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, velocity=None): +def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, pixel_width, speed_thresh=0.03, velocity=None): """Computes occupancy per bin in seconds Parameters @@ -73,6 +76,7 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc edges of histogram in meters edges_y: array-like edges of histogram in meters + pixel_width: array-like speed_thresh: float, optional in meters. Default = 3.0 cm/s velocity: np.ndarray(dtype=float) @@ -89,19 +93,22 @@ def compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh=0.03, veloc sampling_period = (np.max(pos_tt) - np.min(pos_tt)) / len(pos_tt) np.seterr(invalid='ignore') if velocity is None: - is_running = compute_speed(pos, pos_tt) > speed_thresh + if pixel_width[1] is not int(1): + is_running = compute_speed(pos, pos_tt) > speed_thresh + else: + is_running = compute_speed(pos[:, 0], pos_tt) > speed_thresh else: is_running = np.linalg.norm(velocity) > speed_thresh run_pos = pos[is_running, :] - occupancy = np.histogram2d(run_pos[:, 0], - run_pos[:, 1], - [edges_x, edges_y])[0] * sampling_period # in seconds + occupancy = np.histogram2d(run_pos[:, 1], + run_pos[:, 0], + [edges_y, edges_x])[0] * sampling_period # in seconds return occupancy, is_running -def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03, velocity=None): +def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, pixel_width, speed_thresh=0.03, velocity=None): """Returns speed-gated position during spikes Parameters @@ -116,6 +123,7 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 edges of histogram in meters edges_y: np.ndarray(dtype=float) edges of histogram in meters + pixel_width: array speed_thresh: float in meters. Default = 3.0 cm/s velocity: np.ndarray(dtype=float) @@ -126,7 +134,10 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 """ np.seterr(invalid='ignore') if velocity is None: - is_running = compute_speed(pos, pos_tt) > speed_thresh + if pixel_width[1] is not int(1): + is_running = compute_speed(pos, pos_tt) > speed_thresh + else: + is_running = compute_speed(pos[:, 0], pos_tt) > speed_thresh else: is_running = np.linalg.norm(velocity) > speed_thresh @@ -134,9 +145,9 @@ def compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh=0.03 spike_pos_inds = spike_pos_inds[is_running[spike_pos_inds]] pos_on_spikes = pos[spike_pos_inds, :] - n_spikes = np.histogram2d(pos_on_spikes[:, 0], - pos_on_spikes[:, 1], - [edges_x, edges_y])[0] + n_spikes = np.histogram2d(pos_on_spikes[:, 1], + pos_on_spikes[:, 0], + [edges_y, edges_x])[0] return n_spikes @@ -160,7 +171,7 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, (time,) in seconds spikes: np.ndarray(dtype=float) in seconds - pixel_width: float + pixel_width: array-like speed_thresh: float, optional in meters. Default = 3.0 cm/s gaussian_sd_x: float, optional @@ -190,12 +201,12 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, y_start = np.nanmin(pos[:, 1]) if y_start is None else y_start y_stop = np.nanmax(pos[:, 1]) if y_stop is None else y_stop - edges_x = np.arange(x_start, x_stop, pixel_width[1]) - edges_y = np.arange(y_start, y_stop, pixel_width[0]) + edges_x = np.arange(x_start, x_stop, pixel_width[0]) + edges_y = np.arange(y_start, y_stop, pixel_width[1]) - occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, speed_thresh, velocity) + occupancy, running = compute_2d_occupancy(pos, pos_tt, edges_x, edges_y, pixel_width, speed_thresh, velocity) - n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, speed_thresh, velocity) + n_spikes = compute_2d_n_spikes(pos, pos_tt, spikes, edges_x, edges_y, pixel_width, speed_thresh, velocity) np.seterr(divide='ignore') firing_rate = n_spikes / occupancy # in Hz diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 1d9dae41..b77a7d94 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -34,8 +34,10 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, self.pos_tt = get_timeseries_tt(spatial_series) if velocity is not None: self.velocity = velocity + self.disable = False else: self.velocity = None + self.disable = True self.get_position(spatial_series) @@ -74,9 +76,9 @@ def get_controls(self): style = {'description_width': 'initial'} bft_gaussian_x = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd x (cm)', style=style) bft_gaussian_y = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd y (cm)', style=style) - bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (cm/s)', style=style) + bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (m/s)', style=style) dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') - cb_velocity = Checkbox(value=False, description='use velocity', indent=False) + cb_velocity = Checkbox(value=False, description='use velocity', indent=False, disabled= self.disable) return bft_gaussian_x, bft_gaussian_y, bft_speed, dd_unit_select, cb_velocity From 6c68ddfb9a5699bbe08e8816f30b9dc6792d2b7b Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Mon, 25 Jan 2021 18:50:08 -0500 Subject: [PATCH 36/39] Added bin number control Made modifications to: - NaN out any unexplored areas in the place field - Add the ability to control the number of bins in the place field calculation. --- nwbwidgets/analysis/placefields.py | 2 +- nwbwidgets/placefield.py | 32 +++++++++++++++++------------- 2 files changed, 19 insertions(+), 15 deletions(-) diff --git a/nwbwidgets/analysis/placefields.py b/nwbwidgets/analysis/placefields.py index ce585b99..30bb4c51 100644 --- a/nwbwidgets/analysis/placefields.py +++ b/nwbwidgets/analysis/placefields.py @@ -217,7 +217,7 @@ def compute_2d_firing_rate(pos, pos_tt, spikes, # filter occupancy to create a mask so non-explored regions are nan'ed sigmas_occ = [gaussian_sd_y / pixel_width[1] / 8, gaussian_sd_x / pixel_width[0] / 8] filtered_occupancy = gaussian_filter(occupancy, sigmas_occ) - # filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan + filtered_firing_rate[filtered_occupancy.astype('bool') < .00001] = np.nan return occupancy, filtered_firing_rate, [edges_x, edges_y] diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index b77a7d94..0c9e895a 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -41,13 +41,12 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, self.get_position(spatial_series) - self.get_pixel_width() - - bft_gaussian_x, bft_gaussian_y, bft_speed, dd_unit_select, cb_velocity = self.get_controls() + bft_gaussian_x, bft_gaussian_y, bft_bin_num, bft_speed, dd_unit_select, cb_velocity = self.get_controls() self.controls = dict( gaussian_sd_x=bft_gaussian_x, gaussian_sd_y=bft_gaussian_y, + bin_num=bft_bin_num, speed_thresh=bft_speed, index=dd_unit_select, use_velocity=cb_velocity @@ -59,6 +58,7 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, widgets.VBox([ bft_gaussian_x, bft_gaussian_y, + bft_bin_num, bft_speed, dd_unit_select, cb_velocity, @@ -66,8 +66,8 @@ def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, vis2widget(out_fig) ] - def get_pixel_width(self): - self.pixel_width = [(np.nanmax(self.pos) - np.nanmin(self.pos)) / 1000] * 2 + def get_pixel_width(self, bin_num): + self.pixel_width = [(np.nanmax(self.pos) - np.nanmin(self.pos)) / bin_num] * 2 def get_position(self, spatial_series): self.pos, self.unit = get_timeseries_in_units(spatial_series) @@ -76,18 +76,22 @@ def get_controls(self): style = {'description_width': 'initial'} bft_gaussian_x = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd x (cm)', style=style) bft_gaussian_y = BoundedFloatText(value=0.0184, min=0, max=99999, description='gaussian sd y (cm)', style=style) + bft_bin_num = BoundedFloatText(value=1000, min=0, max=99999, description='number of bins', style=style) bft_speed = BoundedFloatText(value=0.03, min=0, max=99999, description='speed threshold (m/s)', style=style) dd_unit_select = Dropdown(options=np.arange(len(self.units)), description='unit') cb_velocity = Checkbox(value=False, description='use velocity', indent=False, disabled= self.disable) - return bft_gaussian_x, bft_gaussian_y, bft_speed, dd_unit_select, cb_velocity - - def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian_sd_y=0.0184, use_velocity=False): - occupancy, filtered_firing_rate, [edges_x, edges_y] = self.compute_twodim_firing_rate(index=index, - speed_thresh=speed_thresh, - gaussian_sd_x=gaussian_sd_x, - gaussian_sd_y=gaussian_sd_y, - use_velocity=use_velocity) + return bft_gaussian_x, bft_gaussian_y, bft_bin_num, bft_speed, dd_unit_select, cb_velocity + + def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian_sd_y=0.0184, bin_num=1000, + use_velocity=False): + self.get_pixel_width(bin_num) + occupancy, filtered_firing_rate, [edges_x, edges_y] = self.compute_twodim_firing_rate(self.pixel_width[0], + index=index, + speed_thresh=speed_thresh, + gaussian_sd_x=gaussian_sd_x, + gaussian_sd_y=gaussian_sd_y, + use_velocity=use_velocity) fig, ax = plt.subplots() im = ax.imshow(filtered_firing_rate, @@ -102,7 +106,7 @@ def do_rate_map(self, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian return fig @lru_cache() - def compute_twodim_firing_rate(self, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian_sd_y=0.0184, + def compute_twodim_firing_rate(self, pixel_width, index=0, speed_thresh=0.03, gaussian_sd_x=0.0184, gaussian_sd_y=0.0184, use_velocity=False): tmin = min(self.pos_tt) tmax = max(self.pos_tt) From 6e3c1b176d89e77934e6a3fa160140a8f913b30b Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Mon, 8 Mar 2021 21:47:20 -0500 Subject: [PATCH 37/39] Simple formatting --- nwbwidgets/placefield.py | 3 +++ nwbwidgets/view.py | 3 +-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/nwbwidgets/placefield.py b/nwbwidgets/placefield.py index 0c9e895a..a7c973d3 100644 --- a/nwbwidgets/placefield.py +++ b/nwbwidgets/placefield.py @@ -12,6 +12,7 @@ from .utils.units import get_spike_times from .utils.timeseries import get_timeseries_in_units, get_timeseries_tt + def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): if spatial_series.data.shape[1] == 2: return PlaceFieldWidget(spatial_series) @@ -21,6 +22,7 @@ def route_placefield(spatial_series: pynwb.behavior.SpatialSeries): print('Spatial series exceeds dimensionality for visualization') return + class PlaceFieldWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, @@ -126,6 +128,7 @@ def compute_twodim_firing_rate(self, pixel_width, index=0, speed_thresh=0.03, ga velocity=self.velocity) return occupancy, filtered_firing_rate, [edges_x, edges_y] + class PlaceField1DWidget(widgets.HBox): def __init__(self, spatial_series: pynwb.behavior.SpatialSeries, velocity: pynwb.TimeSeries = None, diff --git a/nwbwidgets/view.py b/nwbwidgets/view.py index a35a1317..da34846e 100644 --- a/nwbwidgets/view.py +++ b/nwbwidgets/view.py @@ -49,8 +49,7 @@ def show_dynamic_table(node, **kwargs) -> widgets.Widget: pynwb.behavior.SpatialSeries: OrderedDict({ 'over time': timeseries.SeparateTracesPlotlyWidget, 'trace': behavior.plotly_show_spatial_trace, - 'rate map': placefield.route_placefield, - '1D rate map': placefield.PlaceField1DWidget}), + 'rate map': placefield.route_placefield}), pynwb.image.GrayscaleImage: image.show_grayscale_image, pynwb.image.RGBImage: image.show_rbga_image, pynwb.image.RGBAImage: image.show_rbga_image, From c9fa416a2f7177d9a815bfd5302c3f6b57752f8d Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Mon, 8 Mar 2021 22:46:58 -0500 Subject: [PATCH 38/39] Update imports --- nwbwidgets/view.py | 1 + 1 file changed, 1 insertion(+) diff --git a/nwbwidgets/view.py b/nwbwidgets/view.py index a10f33ca..2813882d 100644 --- a/nwbwidgets/view.py +++ b/nwbwidgets/view.py @@ -20,6 +20,7 @@ timeseries, file, spectrum, + placefield ) From 87a24336c0ad18f7640530fc49f559a8a6f1ca29 Mon Sep 17 00:00:00 2001 From: Michael Scheid Date: Mon, 17 May 2021 11:28:38 -0400 Subject: [PATCH 39/39] Removed old ignore files Removed outdated files being ignored --- .gitignore | 3 --- 1 file changed, 3 deletions(-) diff --git a/.gitignore b/.gitignore index cfcd4378..fba6f61b 100644 --- a/.gitignore +++ b/.gitignore @@ -129,6 +129,3 @@ dmypy.json # nwb **.nwb -nwbwidgets/controllers/group_and_sort_controllers_multi_select.py -nwbwidgets/placefield_group_and_sort.py -nwbwidgets/towers_task_placefield.py