Merge branch 'develop' into one_globus

Author: mayofaulkner

.github/workflows/ibllib_ci.yml

@@ -17,8 +17,8 @@ jobs:
       fail-fast: false  # Whether to stop execution of other instances
       max-parallel: 4
       matrix:
-        os: ["ubuntu-latest"] #  "windows-latest"]  #  , "macos-latest"
-        python-version: ["3.8"]  # "3.7",
+        os: ["ubuntu-latest"]
+        python-version: ["3.8"]
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}

MANIFEST.in

@@ -1,6 +1,7 @@
 include ibllib/atlas/allen_structure_tree.csv
 include ibllib/atlas/beryl.npy
 include ibllib/atlas/cosmos.npy
+include ibllib/atlas/swanson.npy
 include ibllib/atlas/mappings.pqt
 include ibllib/io/extractors/extractor_types.json
 include brainbox/tests/wheel_test.p

README.md

@@ -12,9 +12,9 @@ The library is currently 2 main modules:
 [Release Notes here](release_notes.md)
 
 ## Requirements
-**OS**: Deployed on Linux and Windows. Minimally tested for Mac.
+**OS**: Only tested on Linux. Windows and Mac may work, but are not supported.
 
-**Python Module**: Python 3.7 or higher, we develop on 3.8.
+**Python Module**: Python 3.8 or higher
 
 ## Installation, documentation and examples
 https://docs.internationalbrainlab.org
@@ -25,7 +25,7 @@ See https://int-brain-lab.github.io/iblenv/07_contribution.html
 
 We use gitflow and Semantic Versioning.
 
-Before commiting to your branch:
+Before committing to your branch:
 -   run tests
 -   flake8
 This is also enforced by continuous integration.

brainbox/__init__.py

@@ -1,4 +1,5 @@
+import logging
 try:
     import one
 except ModuleNotFoundError:
-    logging.getLogger('ibllib').error('Missing dependency, please run `pip install ONE-api`')
+    logging.getLogger(__name__).error('Missing dependency, please run `pip install ONE-api`')

brainbox/behavior/dlc.py

@@ -11,7 +11,7 @@
 import scipy.interpolate as interpolate
 from scipy.stats import zscore
 
-from ibllib.dsp.smooth import smooth_interpolate_savgol
+from neurodsp.smooth import smooth_interpolate_savgol
 from brainbox.processing import bincount2D
 import brainbox.behavior.wheel as bbox_wheel
 
@@ -350,11 +350,12 @@ def plot_wheel_position(wheel_position, wheel_time, trials_df):
 
     plt.axvline(x=0, linestyle='--', c='k', label='stimOn')
     plt.axhline(y=-0.26, linestyle='--', c='g', label='reward')
+    plt.axhline(y=0.26, linestyle='--', c='g', label='reward')
     plt.ylim([-0.27, 0.27])
     plt.xlabel('time [sec]')
-    plt.ylabel('wheel position [rad]')
+    plt.ylabel('wheel position diff to first value [rad]')
     plt.legend(loc='center right')
-    plt.title('Wheel position')
+    plt.title('Wheel position trial avg\n(and individual trials)')
     plt.tight_layout()
 
     return plt.gca()
@@ -406,7 +407,7 @@ def plot_lick_hist(lick_times, trials_df):
     plt.plot(times, pd.DataFrame.from_dict(dict(zip(correct.index, incorrect.values))).mean(axis=1),
              c='gray', label='incorrect trial')
     plt.axvline(x=0, label='feedback', linestyle='--', c='purple')
-    plt.title('Lick events')
+    plt.title('Lick events trial avg')
     plt.xlabel('time [sec]')
     plt.ylabel('lick events [a.u.]')
     plt.legend(loc='lower right')
@@ -457,16 +458,17 @@ def plot_motion_energy_hist(camera_dict, trials_df):
                 motion_energy = zscore(camera_dict[cam]['motion_energy'], nan_policy='omit')
                 try:
                     start_idx = insert_idx(camera_dict[cam]['times'], start_window)
+                    end_idx = np.array(start_idx + int(WINDOW_LEN * SAMPLING[cam]), dtype='int64')
+                    me_all = [motion_energy[start_idx[i]:end_idx[i]] for i in range(len(start_idx))]
+                    me_all = [m for m in me_all if len(m) > 0]
+                    times = np.arange(len(me_all[0])) / SAMPLING[cam] + WINDOW_LAG
+                    me_mean = np.mean(me_all, axis=0)
+                    me_std = np.std(me_all, axis=0) / np.sqrt(len(me_all))
+                    plt.plot(times, me_mean, label=f'{cam} cam', color=colors[cam], linewidth=2)
+                    plt.fill_between(times, me_mean + me_std, me_mean - me_std, color=colors[cam], alpha=0.2)
                 except ValueError:
-                    logger.error("Camera.times are outside of the trial windows")
-                    raise
-                end_idx = np.array(start_idx + int(WINDOW_LEN * SAMPLING[cam]), dtype='int64')
-                me_all = [motion_energy[start_idx[i]:end_idx[i]] for i in range(len(start_idx))]
-                times = np.arange(len(me_all[0])) / SAMPLING[cam] + WINDOW_LAG
-                me_mean = np.mean(me_all, axis=0)
-                me_std = np.std(me_all, axis=0) / np.sqrt(len(me_all))
-                plt.plot(times, me_mean, label=f'{cam} cam', color=colors[cam], linewidth=2)
-                plt.fill_between(times, me_mean + me_std, me_mean - me_std, color=colors[cam], alpha=0.2)
+                    logger.error(f"{cam}Camera camera.times are outside of the trial windows")
+                    missing_data.append(cam)
             except AttributeError:
                 logger.warning(f"Cannot load motion energy and/or times data for {cam} camera")
                 missing_data.append(cam)
@@ -479,11 +481,11 @@ def plot_motion_energy_hist(camera_dict, trials_df):
     plt.xlabel('time [sec]')
     plt.axvline(x=0, label='stimOn', linestyle='--', c='k')
     plt.legend(loc='lower right')
-    plt.title('Motion Energy')
+    plt.title('Motion Energy trial avg\n(+/- std)')
     if len(missing_data) > 0:
         ax = plt.gca()
         ax.text(.95, .35, f"Data incomplete for\n{' and '.join(missing_data)} camera", color='r', fontsize=10,
-                horizontalalignment='right', verticalalignment='center', transform=ax.transAxes)
+                fontweight='bold', horizontalalignment='right', verticalalignment='center', transform=ax.transAxes)
     return plt.gca()
 
 
@@ -523,7 +525,7 @@ def plot_speed_hist(dlc_df, cam_times, trials_df, feature='paw_r', cam='left', l
     plt.plot(times, pd.DataFrame.from_dict(dict(zip(incorrect.index, incorrect.values))).mean(axis=1),
              c='gray', label='incorrect trial')
     plt.axvline(x=0, label='stimOn', linestyle='--', c='r')
-    plt.title(f'{feature.split("_")[0].capitalize()} speed ({cam} cam)')
+    plt.title(f'{feature.capitalize()} speed trial avg\n({cam.upper()} cam)')
     plt.xticks([-0.5, 0, 0.5, 1, 1.5])
     plt.xlabel('time [sec]')
     plt.ylabel('speed [px/sec]')
@@ -559,8 +561,8 @@ def plot_pupil_diameter_hist(pupil_diameter, cam_times, trials_df, cam='left'):
         plt.plot(times, pupil_mean, label=align_to.split("_")[0], color=color)
         plt.fill_between(times, pupil_mean + pupil_std, pupil_mean - pupil_std, color=color, alpha=0.5)
     plt.axvline(x=0, linestyle='--', c='k')
-    plt.title(f'Pupil diameter ({cam} cam)')
+    plt.title(f'Pupil diameter trial avg\n({cam.upper()} cam)')
     plt.xlabel('time [sec]')
     plt.xticks([-0.5, 0, 0.5, 1, 1.5])
-    plt.ylabel('pupil diameter [px]')
+    plt.ylabel('z-scored smoothed pupil diameter [px]')
     plt.legend(loc='lower right', title='aligned to')

brainbox/behavior/training.py

@@ -111,14 +111,14 @@ def get_sessions(subj, date=None, one=None):
         latest_sess = date
 
     sessions = one.alyx.rest('sessions', 'list', subject=subj, date_range=[latest_minus_week,
-                             latest_sess], dataset_types='trials.intervals')
+                             latest_sess], dataset_types='trials.goCueTrigger_times')
 
     # If not enough sessions in the last week, then just fetch them all
     if len(sessions) < 3:
         specified_date_plus = (specified_date + datetime.timedelta(days=1)).strftime("%Y-%m-%d")
         django_query = 'start_time__lte,' + specified_date_plus
         sessions = one.alyx.rest('sessions', 'list', subject=subj,
-                                 dataset_types='trials.intervals', django=django_query)
+                                 dataset_types='trials.goCueTrigger_times', django=django_query)
 
         # If still 0 sessions then return with warning
         if len(sessions) == 0:
@@ -349,7 +349,7 @@ def compute_training_info(trials, trials_all):
     perf_easy = np.array([compute_performance_easy(trials[k]) for k in trials.keys()])
     n_trials = np.array([compute_n_trials(trials[k]) for k in trials.keys()])
     psych = compute_psychometric(trials_all, signed_contrast=signed_contrast)
-    rt = compute_median_reaction_time(trials_all, signed_contrast=signed_contrast)
+    rt = compute_median_reaction_time(trials_all, contrast=0, signed_contrast=signed_contrast)
 
     return perf_easy, n_trials, psych, rt
 
@@ -376,7 +376,7 @@ def compute_bias_info(trials, trials_all):
     n_trials = np.array([compute_n_trials(trials[k]) for k in trials.keys()])
     psych_20 = compute_psychometric(trials_all, signed_contrast=signed_contrast, block=0.2)
     psych_80 = compute_psychometric(trials_all, signed_contrast=signed_contrast, block=0.8)
-    rt = compute_median_reaction_time(trials_all, signed_contrast=signed_contrast)
+    rt = compute_median_reaction_time(trials_all, contrast=0, signed_contrast=signed_contrast)
 
     return perf_easy, n_trials, psych_20, psych_80, rt
 
@@ -452,7 +452,7 @@ def compute_n_trials(trials):
     return trials['choice'].shape[0]
 
 
-def compute_psychometric(trials, signed_contrast=None, block=None):
+def compute_psychometric(trials, signed_contrast=None, block=None, plotting=False):
     """
     Compute psychometric fit parameters for trials object
 
@@ -479,17 +479,27 @@ def compute_psychometric(trials, signed_contrast=None, block=None):
     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]),
-        P_model='erf_psycho_2gammas',
-        parstart=np.array([np.mean(contrasts), 20., 0.05, 0.05]),
-        parmin=np.array([np.min(contrasts), 0., 0., 0.]),
-        parmax=np.array([np.max(contrasts), 100., 1, 1]))
+    if plotting:
+        psych, _ = psy.mle_fit_psycho(
+            np.vstack([contrasts, n_contrasts, prob_choose_right]),
+            P_model='erf_psycho_2gammas',
+            parstart=np.array([0., 40., 0.1, 0.1]),
+            parmin=np.array([-50., 10., 0., 0.]),
+            parmax=np.array([50., 50., 0.2, 0.2]),
+            nfits=10)
+    else:
+
+        psych, _ = psy.mle_fit_psycho(
+            np.vstack([contrasts, n_contrasts, prob_choose_right]),
+            P_model='erf_psycho_2gammas',
+            parstart=np.array([np.mean(contrasts), 20., 0.05, 0.05]),
+            parmin=np.array([np.min(contrasts), 0., 0., 0.]),
+            parmax=np.array([np.max(contrasts), 100., 1, 1]))
 
     return psych
 
 
-def compute_median_reaction_time(trials, stim_on_type='stimOn_times', signed_contrast=None):
+def compute_median_reaction_time(trials, stim_on_type='stimOn_times', contrast=None, signed_contrast=None):
     """
     Compute median reaction time on zero contrast trials from trials object
 
@@ -505,10 +515,15 @@ def compute_median_reaction_time(trials, stim_on_type='stimOn_times', signed_con
     """
     if signed_contrast is None:
         signed_contrast = get_signed_contrast(trials)
-    zero_trials = (trials.response_times - trials[stim_on_type])[signed_contrast == 0]
-    if np.any(zero_trials):
+
+    if contrast is None:
+        contrast_idx = np.full(trials.probabilityLeft.shape, True, dtype=bool)
+    else:
+        contrast_idx = signed_contrast == contrast
+
+    if np.any(contrast_idx):
         reaction_time = np.nanmedian((trials.response_times - trials[stim_on_type])
-                                     [signed_contrast == 0])
+                                     [contrast_idx])
     else:
         reaction_time = np.nan
 
@@ -590,15 +605,15 @@ def plot_psychometric(trials, ax=None, title=None, **kwargs):
     contrasts_fit = np.arange(-100, 100)
 
     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)
+    pars_50 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.5, plotting=True)
     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=True)
-    pars_20 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.2)
+    pars_20 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.2, plotting=True)
     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=True)
-    pars_80 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.8)
+    pars_80 = compute_psychometric(trials, signed_contrast=signed_contrast, block=0.8, plotting=True)
     prob_right_fit_80 = psy.erf_psycho_2gammas(pars_80, contrasts_fit)
 
     cmap = sns.diverging_palette(20, 220, n=3, center="dark")

brainbox/behavior/wheel.py

@@ -213,10 +213,8 @@ def movements(t, pos, freq=1000, pos_thresh=8, t_thresh=.2, min_gap=.1, pos_thre
         if i2proc[-1] == t.size - 1:
             break
 
-    moving = max_disp > pos_thresh  # for each window is the change in position greater than
-    # our threshold?
-    moving = np.insert(moving, 0, False)  # First sample should always be not moving to ensure
-    # we have an onset
+    moving = max_disp > pos_thresh  # for each window is the change in position greater than our threshold?
+    moving = np.insert(moving, 0, False)  # First sample should always be not moving to ensure we have an onset
     moving[-1] = False  # Likewise, ensure we always end on an offset
 
     onset_samps = np.where(~moving[:-1] & moving[1:])[0]

brainbox/ephys_plots.py

@@ -496,3 +496,30 @@ def histc(x, bins):
         return data.convert2dict(), fig, ax
 
     return data
+
+
+def image_raw_data(raw, fs, chn_coords=None, cmap='bone', title=None, display=False, gain=-90, **kwargs):
+
+    def gain2level(gain):
+        return 10 ** (gain / 20) * 4 * np.array([-1, 1])
+
+    ylabel = 'Channel index' if chn_coords is None else 'Distance from probe tip (um)'
+    title = title or 'Raw data'
+
+    y = np.arange(raw.shape[1]) if chn_coords is None else chn_coords[:, 1]
+
+    x = np.array([0, raw.shape[0] - 1]) / fs * 1e3
+
+    data = ImagePlot(raw, y=y, cmap=cmap)
+    data.set_labels(title=title, xlabel='Time (ms)',
+                    ylabel=ylabel, clabel='Power (uV)')
+    clim = gain2level(gain)
+    data.set_clim(clim=clim)
+    data.set_xlim(xlim=x)
+    data.set_ylim()
+
+    if display:
+        ax, fig = plot_image(data.convert2dict(), **kwargs)
+        return data.convert2dict(), fig, ax
+
+    return data