Merge branch 'release/2.10.0'

Author: mayofaulkner

brainbox/behavior/training.py

@@ -408,7 +408,7 @@ def compute_performance_easy(trials):
     return np.sum(trials['feedbackType'][easy_trials] == 1) / easy_trials.shape[0]
 
 
-def compute_performance(trials, signed_contrast=None, block=None):
+def compute_performance(trials, signed_contrast=None, block=None, prob_right=False):
     """
     Compute performance on all trials at each contrast level from trials object
 
@@ -429,12 +429,16 @@ def compute_performance(trials, signed_contrast=None, block=None):
         return np.nan * np.zeros(3)
 
     contrasts, n_contrasts = np.unique(signed_contrast[block_idx], return_counts=True)
-    rightward = trials.choice == -1
-    # Calculate the proportion rightward for each contrast type
-    prob_choose_right = np.vectorize(lambda x: np.mean(rightward[(x == signed_contrast) &
-                                                                 block_idx]))(contrasts)
 
-    return prob_choose_right, contrasts, n_contrasts
+    if not prob_right:
+        correct = trials.feedbackType == 1
+        performance = np.vectorize(lambda x: np.mean(correct[(x == signed_contrast) & block_idx]))(contrasts)
+    else:
+        rightward = trials.choice == -1
+        # Calculate the proportion rightward for each contrast type
+        performance = np.vectorize(lambda x: np.mean(rightward[(x == signed_contrast) & block_idx]))(contrasts)
+
+    return performance, contrasts, n_contrasts
 
 
 def compute_n_trials(trials):
@@ -472,7 +476,8 @@ def compute_psychometric(trials, signed_contrast=None, block=None):
     if not np.any(block_idx):
         return np.nan * np.zeros(4)
 
-    prob_choose_right, contrasts, n_contrasts = compute_performance(trials, signed_contrast=signed_contrast, block=block)
+    prob_choose_right, contrasts, n_contrasts = compute_performance(trials, signed_contrast=signed_contrast, block=block,
+                                                                    prob_right=True)
 
     psych, _ = psy.mle_fit_psycho(
         np.vstack([contrasts, n_contrasts, prob_choose_right]),
@@ -584,15 +589,15 @@ def plot_psychometric(trials, ax=None, title=None, **kwargs):
     signed_contrast = get_signed_contrast(trials)
     contrasts_fit = np.arange(-100, 100)
 
-    prob_right_50, contrasts_50, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.5)
+    prob_right_50, contrasts_50, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.5, prob_right=True)
     pars_50 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.5)
     prob_right_fit_50 = psy.erf_psycho_2gammas(pars_50, contrasts_fit)
 
-    prob_right_20, contrasts_20, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.2)
+    prob_right_20, contrasts_20, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.2, prob_right=True)
     pars_20 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.2)
     prob_right_fit_20 = psy.erf_psycho_2gammas(pars_20, contrasts_fit)
 
-    prob_right_80, contrasts_80, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.8)
+    prob_right_80, contrasts_80, _ = compute_performance(trials, signed_contrast=signed_contrast, block=0.8, prob_right=True)
     pars_80 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.8)
     prob_right_fit_80 = psy.erf_psycho_2gammas(pars_80, contrasts_fit)
 

brainbox/ephys_plots.py

@@ -25,8 +25,9 @@ def image_lfp_spectrum_plot(lfp_power, lfp_freq, chn_coords=None, chn_inds=None,
 
     ylabel = 'Channel index' if chn_coords is None else 'Distance from probe tip (um)'
     title = title or 'LFP Power Spectrum'
-    chn_inds = chn_inds or np.arange(lfp_power.shape[1])
+
     y = np.arange(lfp_power.shape[1]) if chn_coords is None else chn_coords[:, 1]
+    chn_inds = np.arange(lfp_power.shape[1]) if chn_inds is None else chn_inds
 
     freq_idx = np.where((lfp_freq >= freq_range[0]) & (lfp_freq < freq_range[1]))[0]
     freqs = lfp_freq[freq_idx]
@@ -80,7 +81,7 @@ def image_rms_plot(rms_amps, rms_times, chn_coords=None, chn_inds=None, avg_acro
 
     ylabel = 'Channel index' if chn_coords is None else 'Distance from probe tip (um)'
     title = title or f'{band} RMS'
-    chn_inds = chn_inds or np.arange(rms_amps.shape[1])
+    chn_inds = np.arange(rms_amps.shape[1]) if chn_inds is None else chn_inds
     y = np.arange(rms_amps.shape[1]) if chn_coords is None else chn_coords[:, 1]
 
     rms = rms_amps[:, chn_inds]

brainbox/examples/docs_load_spike_sorting.py

@@ -23,7 +23,7 @@
 
 insertions = one.alyx.rest('insertions', 'list')
 pid = insertions[0]['id']
-sl = SpikeSortingLoader(pid, one=one, atlas=ba)
+sl = SpikeSortingLoader(pid=pid, one=one, atlas=ba)
 spikes, clusters, channels = sl.load_spike_sorting()
 clusters_labeled = SpikeSortingLoader.merge_clusters(spikes, clusters, channels)
 

brainbox/examples/docs_wheel_moves.ipynb

@@ -577,9 +577,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.9.7"
   }
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}

brainbox/task/closed_loop.py

@@ -11,7 +11,7 @@
 
 
 def responsive_units(spike_times, spike_clusters, event_times,
-                     pre_time=[0.5, 0], post_time=[0, 0.5], alpha=0.05):
+                     pre_time=[0.5, 0], post_time=[0, 0.5], alpha=0.05, use_fr=False):
     """
     Determine responsive neurons by doing a Wilcoxon Signed-Rank test between a baseline period
     before a certain task event (e.g. stimulus onset) and a period after the task event.
@@ -31,6 +31,8 @@ def responsive_units(spike_times, spike_clusters, event_times,
         time (in seconds) to follow the event times
     alpha : float
         alpha to use for statistical significance
+    use_fr : bool
+        whether to use the firing rate instead of total spike count
 
     Returns
     -------
@@ -51,18 +53,12 @@ def responsive_units(spike_times, spike_clusters, event_times,
     times = np.column_stack(((event_times + post_time[0]), (event_times + post_time[1])))
     spike_counts, cluster_ids = get_spike_counts_in_bins(spike_times, spike_clusters, times)
 
-    # Do statistics
-    p_values = np.empty(spike_counts.shape[0])
-    stats = np.empty(spike_counts.shape[0])
-    for i in range(spike_counts.shape[0]):
-        if np.sum(baseline_counts[i, :] - spike_counts[i, :]) == 0:
-            p_values[i] = 1
-            stats[i] = 0
-        else:
-            stats[i], p_values[i] = wilcoxon(baseline_counts[i, :], spike_counts[i, :])
+    if use_fr:
+        baseline_counts = baseline_counts / (pre_time[0] - pre_time[1])
+        spike_counts = spike_counts / (post_time[1] - post_time[0])
 
-    # Perform FDR correction for multiple testing
-    sig_units, p_values, _, _ = multipletests(p_values, alpha, method='fdr_bh')
+    # Do statistics
+    sig_units, stats, p_values = compute_comparison_statistics(baseline_counts, spike_counts, test='signrank', alpha=alpha)
     significant_units = cluster_ids[sig_units]
 
     return significant_units, stats, p_values, cluster_ids
@@ -125,27 +121,66 @@ def differentiate_units(spike_times, spike_clusters, event_times, event_groups,
     counts_2, cluster_ids = get_spike_counts_in_bins(spike_times, spike_clusters, times_2)
 
     # Do statistics
-    p_values = np.empty(len(cluster_ids))
-    stats = np.empty(len(cluster_ids))
-    for i in range(len(cluster_ids)):
-        if (np.sum(counts_1[i, :]) == 0) and (np.sum(counts_2[i, :]) == 0):
-            p_values[i] = 1
-            stats[i] = 0
+    sig_units, stats, p_values = compute_comparison_statistics(counts_1, counts_2, test=test, alpha=alpha)
+    significant_units = cluster_ids[sig_units]
+
+    return significant_units, stats, p_values, cluster_ids
+
+
+def compute_comparison_statistics(value1, value2, test='ranksums', alpha=0.05):
+    """
+    Compute statistical test between two arrays
+
+    Parameters
+    ----------
+    value1 : 1D array
+        first array of values to compare
+    value2 : 1D array
+        second array of values to compare
+    test : string
+        which statistical test to use, options are:
+            'ranksums'      Wilcoxon Rank Sums test
+            'signrank'      Wilcoxon Signed Rank test (for paired observations)
+            'ttest'         independent samples t-test
+            'paired_ttest'  paired t-test
+    alpha : float
+        alpha to use for statistical significance
+
+    Returns
+    -------
+    significant_units : 1D array
+        an array with the indices of values that are significatly modulated
+    stats : 1D array
+        the statistic of the test that was performed
+    p_values : 1D array
+        the p-values of all the values
+    """
+
+    p_values = np.empty(len(value1))
+    stats = np.empty(len(value1))
+    for i in range(len(value1)):
+        if test == 'signrank':
+            if np.sum(value1[i, :] - value2[i, :]) == 0:
+                p_values[i] = 1
+                stats[i] = 0
+            else:
+                stats[i], p_values[i] = wilcoxon(value1[i, :], value2[i, :])
         else:
-            if test == 'ranksums':
-                stats[i], p_values[i] = ranksums(counts_1[i, :], counts_2[i, :])
-            elif test == 'signrank':
-                stats[i], p_values[i] = wilcoxon(counts_1[i, :], counts_2[i, :])
-            elif test == 'ttest':
-                stats[i], p_values[i] = ttest_ind(counts_1[i, :], counts_2[i, :])
-            elif test == 'paired_ttest':
-                stats[i], p_values[i] = ttest_rel(counts_1[i, :], counts_2[i, :])
+            if (np.sum(value1[i, :]) == 0) and (np.sum(value2[i, :]) == 0):
+                p_values[i] = 1
+                stats[i] = 0
+            else:
+                if test == 'ranksums':
+                    stats[i], p_values[i] = ranksums(value1[i, :], value2[i, :])
+                elif test == 'ttest':
+                    stats[i], p_values[i] = ttest_ind(value1[i, :], value2[i, :])
+                elif test == 'paired_ttest':
+                    stats[i], p_values[i] = ttest_rel(value1[i, :], value2[i, :])
 
     # Perform FDR correction for multiple testing
     sig_units, p_values, _, _ = multipletests(p_values, alpha, method='fdr_bh')
-    significant_units = cluster_ids[sig_units]
 
-    return significant_units, stats, p_values, cluster_ids
+    return sig_units, stats, p_values
 
 
 def roc_single_event(spike_times, spike_clusters, event_times,

brainbox/task/passive.py

@@ -17,9 +17,9 @@ def get_on_off_times_and_positions(rf_map):
     Returns
     -------
     rf_map_times: time of each receptive field map frame np.array(len(stim_frames)
-    rf_map_pos: unique position of each pixel on scree np.array(len(x_pos), len(y_pos))
+    rf_map_pos: unique position of each pixel on screen np.array(len(x_pos), len(y_pos))
     rf_stim_frames: for each pixel on screen stores array of stimulus frames where stim onset
-    occured. For both white squares 'on' and black squares 'off'
+    occurred. For both white squares 'on' and black squares 'off'
 
     """
 

brainbox/task/trials.py

@@ -21,6 +21,7 @@ def find_trial_ids(trials, side='all', choice='all', order='trial num', sort='id
     """
     if event:
         idx = ~np.isnan(trials[event])
+        nan_idx = np.where(idx)[0]
     else:
         idx = np.ones_like(trials['feedbackType'], dtype=bool)
 
@@ -118,6 +119,9 @@ def _order_by(_trials, order):
     if side == 'right' and choice == 'incorrect':
         trial_id = _order_by(incor_r, order)
 
+    if event:
+        trial_id = nan_idx[trial_id]
+
     return trial_id, dividers
 
 
@@ -146,8 +150,9 @@ def get_event_aligned_raster(times, events, tbin=0.02, values=None, epoch=[-0.4,
         nbin = t.shape[0]
 
     # remove nan trials
-    events = events[~np.isnan(events)]
-    intervals = np.c_[events + epoch[0], events + epoch[1]]
+    non_nan_events = events[~np.isnan(events)]
+    nan_idx = np.where(~np.isnan(events))
+    intervals = np.c_[non_nan_events + epoch[0], non_nan_events + epoch[1]]
 
     # Remove any trials that are later than the last value in bin_times
     out_intervals = intervals[:, 1] > bin_times[-1]
@@ -170,7 +175,11 @@ def get_event_aligned_raster(times, events, tbin=0.02, values=None, epoch=[-0.4,
                                                     event_raster.shape[1]), np.nan)]
         assert(event_raster.shape[0] == intervals.shape[0])
 
-    return event_raster, t
+    # Reindex if we have removed any nan values
+    all_event_raster = np.full((events.shape[0], event_raster.shape[1]), np.nan)
+    all_event_raster[nan_idx, :] = event_raster
+
+    return all_event_raster, t
 
 
 def get_psth(raster, trial_ids=None):
@@ -298,7 +307,7 @@ def filter_left_right(trials, event_raster, event, contrast, order='trial num'):
     return raster, psth
 
 
-def filter_trials(trials, event_raster, event, contrast, order='trial num', sort='choice'):
+def filter_trials(trials, event_raster, event, contrast=(1, 0.5, 0.25, 0.125, 0.0625, 0), order='trial num', sort='choice'):
     """
     Wrapper to get out psth and raster for trial choice
     :param trials: trials object

brainbox/tests/test_trials.py

@@ -159,13 +159,12 @@ def test_get_event_aligned_rasters(self):
         # Test when nans have trials - these are removed from the raster
         use_trials[10:12] = np.nan
         raster, t = get_event_aligned_raster(spikes, use_trials)
-        assert(raster.shape[0] == len(use_trials) - 2)
-        assert(all(np.isnan(raster[0:5, :]).ravel()))
-        assert(all(np.isnan(raster[-5:, :]).ravel()))
+        assert (raster.shape[0] == len(use_trials))
+        assert (all(np.isnan(raster[10:12, :]).ravel()))
+        assert (all(~np.isnan(raster[12:15, :]).ravel()))
 
         use_trials[0:2] = np.nan
         raster, t = get_event_aligned_raster(spikes, use_trials)
-        assert(raster.shape[0] == len(use_trials) - 4)
-        assert(all(np.isnan(raster[0:3, :]).ravel()))
-        assert(all(~np.isnan(raster[4, :]).ravel()))
-        assert(all(np.isnan(raster[-5:, :]).ravel()))
+        assert (raster.shape[0] == len(use_trials))
+        assert (all(np.isnan(raster[0:2, :]).ravel()))
+        assert (all(np.isnan(raster[-5:, :]).ravel()))