Analyzing intracranial electrophysiology data with xarray

Over the last few years, it has been exciting to see the xarray project evolve, add new functionality, and mature. This post is an attempt at giving xarray another visit to see how it could integrate into electrophysiology workflows.

A quick background on our data

It is common in neuroscience to ask individuals to perform a task over and over again. You record the activity in the brain each time they perform the task (called an “epoch” or a “trial”). Time is recorded relative to some onset when the task begins. That is t==0. The result is usually a matrix of epochs x channejupyls x time. You can do a lot of stuff with this data, but our task in this paper is to detect changes in neural activity at trial onset (t==0).

In our case, we’ve got a small dataset from an old paper of mine. The repository contains several tutorial notebooks and sample data to describe predictive modeling in cognitive neuroscience. You can find the repository here. The task that individuals were performing was passively listening to spoken sentences through a speaker. While they did this, we recorded electrical activity at the surface of their brain (these were surgical patients, and had implanted electrodes under their scalp).

In the Feature Extraction notebook, I covered how to do some simple data manipulation and feature extraction with timeseries analysis. Let’s try to re-create some of the main steps in that tutorial, but now using xarray as an in-memory structure for our data.

Note: The goal here is to learn a bit about xarray moreso than to discuss ecog modeling, so I’ll spend more time talking about my thoughts on the various functions/methods/etc in Xarray than talking about neuroscience.

In this post, we’ll perform a few common processing and extraction steps. The goal is to do a few munging operations that require manipulating data and visualizing simple statistics.

# Imports we'll use later
import mne
import numpy as np
import matplotlib.pyplot as plt
from download import download
import os
from sklearn.preprocessing import scale
import xarray as xr
xr.set_options(display_style="html")

import warnings
warnings.simplefilter('ignore')
%matplotlib inline

We’ll load the data from my GitHub repository (probably not the most efficient way to store or retrieve the data, but hey, this was 3 years ago :-) ).

url_epochs = "https://github.com/choldgraf/paper-encoding_decoding_electrophysiology/blob/master/raw_data/ecog-epo.fif?raw=true"

path_data = download(url_epochs, './ecog-epo.fif', replace=True)
ecog = mne.read_epochs(path_data, preload=True)
os.remove(path_data)

Here’s what the raw data looks like - each horizontal line is electrical activity in a channel over time. The faint vertical green lines show the onset of each trial (they are concatenated together, but in reality there’s a bit of time between trials). This will be one of the last times we use MNE hopefully.

_ = ecog.plot(scalings='auto', n_epochs=5, n_channels=10)
../../../_images/2019-10-22-xarray-neuro_7_0.png

Converting to xarray

First off, we’ll define a helper function that converts the MNE Epochs object into an xarray DataArray object. DataArrays provide an N-Dimensional representation of data, but with the option to include a lot of extra metadata.

DataArrays are useful because you can include information about each dimension of the data. For example, we can tell our DataArray the name, values, and units of each dimension. In this case, in our case one dimension is “time” so we can label it as such.

def epochs_to_dataarray(epochs):
    """A simple function to convert an Epochs object to DataArray"""
    da = xr.DataArray(
    epochs._data,
    dims=['epoch', 'channel', 'time'],
    coords={
        'time': ecog.times,
        'channel': ecog.ch_names,
        'epoch': range(ecog._data.shape[0])
    },
    name='Sample dataset',
    attrs=dict(ecog.info)
    )
    return da

Just look at all the metadata that we were able to pack into the DataArray. Almost all of MNE’s metadata fit nicely into .attrs.

# There's quite a lot of output, so keep scrolling down!
da = epochs_to_dataarray(ecog)
da
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'Sample dataset'
  • epoch: 29
  • channel: 32
  • time: 2250
  • -118.1 -121.3 -117.6 -111.4 -107.3 ... -108.2 -105.9 -105.4 -102.5
    array([[[-118.10999298, -121.30046844, -117.64116669, ...,
               20.78869057,   11.38018703,    8.12319088],
            [-126.95540619, -120.84265137, -105.75999451, ...,
               25.57149315,   28.6223793 ,   24.82316971],
            [ -28.1061306 ,  -16.39686775,  -21.64873505, ...,
              -31.78049469,  -24.28160286,  -25.21776772],
            ...,
            [  18.55451965,   22.87419701,   25.66038895, ...,
               25.84903336,   19.71798515,   14.92882633],
            [ -44.79590607,  -49.5848465 ,  -52.08906555, ...,
              -48.92469025,  -54.11177826,  -60.23026276],
            [-105.09076691, -108.88887024, -113.57571411, ...,
               -3.28264022,    5.5660615 ,   13.18798733]],
    
           [[  20.84839821,   18.30262566,   20.12442017, ...,
              -19.40711594,  -18.51371574,  -20.21596718],
            [  18.0386982 ,    5.35804272,   -6.45382166, ...,
              -42.1801033 ,  -44.59093857,  -47.80654144],
            [  76.62892914,   73.66176605,   64.21563721, ...,
              -12.31007576,  -21.26686478,  -26.83379555],
            ...,
            [  58.06867599,   53.49727631,   51.89107895, ...,
               42.84626007,   52.36317062,   55.75439453],
            [ -24.79452324,  -18.58691406,   -8.12459183, ...,
              -45.95447159,  -42.66623306,  -41.06335449],
            [-110.1750946 , -101.40830994,  -88.24209595, ...,
              -26.7741375 ,  -32.31249619,  -37.07465744]],
    
           [[  44.85819626,   57.10462952,   68.87900543, ...,
               13.73396015,   28.33650208,   40.76254654],
            [  10.10721684,   19.94246101,   31.81379509, ...,
              120.01114655,  119.50100708,  128.0358429 ],
            [ -96.35434723,  -89.9883194 ,  -72.41400909, ...,
              140.28143311,  131.52635193,  133.30158997],
            ...,
            [   1.81244755,   -1.55214143,   -2.54425907, ...,
              -30.9526062 ,  -31.50832367,  -32.25312042],
            [  44.95462799,   38.7992363 ,   33.4192009 , ...,
               64.5262146 ,   62.74245834,   61.11261368],
            [ -37.66891861,  -37.9514122 ,  -40.97248077, ...,
               61.204319  ,   58.57411194,   56.27007675]],
    
           ...,
    
           [[  88.75709534,   83.42398071,   76.10377502, ...,
              -23.7650528 ,  -23.43914986,  -19.55894852],
            [  63.89505005,   63.79450226,   61.43893433, ...,
              -54.35200119,  -59.907547  ,  -69.36830139],
            [  39.31270599,   46.15828705,   51.12403107, ...,
              -30.19447899,  -27.1382637 ,  -19.75760651],
            ...,
            [  -2.09485602,  -13.05238819,  -19.73264503, ...,
              150.31195068,  129.28895569,  114.28322601],
            [-142.40054321, -142.72599792, -140.20773315, ...,
               90.81705475,   82.73035431,   77.06429291],
            [-161.34713745, -160.04364014, -151.092453  , ...,
              238.87109375,  230.41633606,  222.64045715]],
    
           [[  67.85353088,   72.97834778,   70.99682617, ...,
               99.13801575,   95.14775085,   86.59616089],
            [ 169.86994934,  172.12055969,  173.22143555, ...,
              132.95147705,  135.66346741,  141.64862061],
            [  80.4973526 ,   77.63601685,   73.43898773, ...,
               91.37036896,   84.76886749,   82.48825836],
            ...,
            [ -13.24606895,   -8.65564728,    6.53395557, ...,
             -105.62326813, -104.58720398, -107.99065399],
            [  22.20972443,   49.33536148,   85.349823  , ...,
             -103.54721069, -102.84375   , -104.26371002],
            [ -33.03855133,   -9.5574255 ,   10.97984695, ...,
             -148.20103455, -148.95999146, -154.4241333 ]],
    
           [[ -98.52385712,  -99.22042847,  -99.56749725, ...,
               49.72173309,   41.6667099 ,   31.01161194],
            [ -70.90259552,  -66.83739471,  -68.66448975, ...,
               35.73180008,   33.55708694,   28.1260643 ],
            [ -29.78354645,  -34.5358963 ,  -40.23687363, ...,
              -39.39603043,  -42.75005341,  -50.70994186],
            ...,
            [ 110.57292175,  103.35498047,   95.09171295, ...,
             -162.01602173, -164.45585632, -168.76620483],
            [  38.4355278 ,   39.19929504,   43.87767792, ...,
             -108.41282654, -106.09828949, -103.17415619],
            [  98.9158783 ,  105.16780853,  116.0725174 , ...,
             -105.856987  , -105.39962769, -102.4755249 ]]])
    • time
      (time)
      float64
      -1.5 -1.497 -1.493 ... 5.993 5.997
      array([-1.5     , -1.496667, -1.493333, ...,  5.99    ,  5.993333,  5.996667])
    • channel
      (channel)
      'ch_0' 'ch_1' ... 'ch_30' 'ch_31'
      array(['ch_0', 'ch_1', 'ch_2', 'ch_3', 'ch_4', 'ch_5', 'ch_6', 'ch_7', 'ch_8',
             'ch_9', 'ch_10', 'ch_11', 'ch_12', 'ch_13', 'ch_14', 'ch_15', 'ch_16',
             'ch_17', 'ch_18', 'ch_19', 'ch_20', 'ch_21', 'ch_22', 'ch_23', 'ch_24',
             'ch_25', 'ch_26', 'ch_27', 'ch_28', 'ch_29', 'ch_30', 'ch_31'],
            dtype='<U5')
    • epoch
      (epoch)
      int64
      0 1 2 3 4 5 6 ... 23 24 25 26 27 28
      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28])
  • file_id :
    {'version': 65539, 'machid': array([808661043, 808661043], dtype=int32), 'secs': 1471644431, 'usecs': 91610}
    events :
    []
    hpi_results :
    []
    hpi_meas :
    []
    subject_info :
    None
    device_info :
    None
    helium_info :
    None
    hpi_subsystem :
    None
    proc_history :
    []
    meas_id :
    {'version': 65539, 'machid': array([808661043, 808661043], dtype=int32), 'secs': 1471644431, 'usecs': 91754}
    experimenter :
    None
    description :
    None
    proj_id :
    None
    proj_name :
    None
    meas_date :
    (0, 0)
    utc_offset :
    None
    sfreq :
    300.0
    highpass :
    0.0
    lowpass :
    500.0
    line_freq :
    None
    gantry_angle :
    None
    chs :
    [{'scanno': 1, 'logno': 1, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_0', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 2, 'logno': 2, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_1', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 3, 'logno': 3, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_2', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 4, 'logno': 4, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_3', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 5, 'logno': 5, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_4', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 6, 'logno': 6, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_5', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 7, 'logno': 7, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_6', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 8, 'logno': 8, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_7', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 9, 'logno': 9, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_8', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 10, 'logno': 10, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_9', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 11, 'logno': 11, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_10', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 12, 'logno': 12, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_11', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 13, 'logno': 13, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_12', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 14, 'logno': 14, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_13', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 15, 'logno': 15, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_14', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 16, 'logno': 16, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_15', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 17, 'logno': 17, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_16', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 18, 'logno': 18, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_17', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 19, 'logno': 19, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_18', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 20, 'logno': 20, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_19', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 21, 'logno': 21, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_20', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 22, 'logno': 22, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_21', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 23, 'logno': 23, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_22', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 24, 'logno': 24, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_23', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 25, 'logno': 25, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_24', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 26, 'logno': 26, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_25', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 27, 'logno': 27, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_26', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 28, 'logno': 28, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_27', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 29, 'logno': 29, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_28', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 30, 'logno': 30, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_29', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 31, 'logno': 31, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_30', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 32, 'logno': 32, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_31', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}]
    dev_head_t :
    <Transform | MEG device->head> [[1. 0. 0. 0.] [0. 1. 0. 0.] [0. 0. 1. 0.] [0. 0. 0. 1.]]
    ctf_head_t :
    None
    dev_ctf_t :
    None
    dig :
    []
    bads :
    []
    ch_names :
    ['ch_0', 'ch_1', 'ch_2', 'ch_3', 'ch_4', 'ch_5', 'ch_6', 'ch_7', 'ch_8', 'ch_9', 'ch_10', 'ch_11', 'ch_12', 'ch_13', 'ch_14', 'ch_15', 'ch_16', 'ch_17', 'ch_18', 'ch_19', 'ch_20', 'ch_21', 'ch_22', 'ch_23', 'ch_24', 'ch_25', 'ch_26', 'ch_27', 'ch_28', 'ch_29', 'ch_30', 'ch_31']
    nchan :
    32
    projs :
    []
    comps :
    []
    acq_pars :
    None
    acq_stim :
    None
    custom_ref_applied :
    False
    xplotter_layout :
    None
    kit_system_id :
    None

The data consists of many trials, channels, and timepoints. Let’s start by selecting a time region within each trial that we can visualize more cleanly.

Subsetting out data with da.sel

In xarray, we select items with the sel and isel method. This behaves kind of like the pandas loc and iloc methods, however because we have named dimensions, we can directly specify them in our call.

# We'll drop a subset of timepoints for visualization
da = da.sel(time=slice(-1, 3))

Now let’s calculate the average across all epochs for each electrode/time point. This is a reduction of our data array, in that it reduces the number of dimensions. Xarray has many of the same statistical methods that NumPy does. An interesting twist is that you can specify named dimensions instead of simply an axis=<integer> argument. In addition, we’ll choose the colors that we’ll use for cycling through our channels - because we can quickly reference the channels axis by name, we don’t need to remember which axis corresponds to channels.

fig, ax = plt.subplots(figsize=(15, 5))
n_channels = da['channel'].shape[0]
ax.set_prop_cycle(color=plt.cm.viridis(np.linspace(0, 1, n_channels)))
da.mean(dim='epoch').plot.line(x='time', hue='channel')
ax.get_legend().remove()
../../../_images/2019-10-22-xarray-neuro_15_0.png

It doesn’t look like much is going on…let’s see if we can clean it up a bit.

De-meaning the data with da.where

First off - we’ll subtract the “pre-baseline mean” from each trial. This makes it easier to visualize how each channel’s activity changed at time == 0.

To accomplish this we’ll use da.where. This takes some kind of boolean-style mask, does a bunch of clever projections according to the names of coordinates, and returns the dataarray masked values removed (as NaNs) and other values unchanged. We can use this to calculate the mean of each channel / epoch only for the pre-baseline timepoints.

# This returns a version of the data array with NaNs where the query is False
# The dimensions will intelligently broadcast 
prebaseline_mean = da.where(da.time < 0).mean(dim='time')
da_demeaned = da - prebaseline_mean

Now we can visualize the de-baseline-meaned data

fig, ax = plt.subplots(figsize=(15, 5))
ax.set_prop_cycle(color=plt.cm.viridis(np.linspace(0, 1, da['channel'].shape[0])))
da_demeaned.mean(dim='epoch').plot.line(x='time', hue='channel')
ax.get_legend().remove()
../../../_images/2019-10-22-xarray-neuro_20_0.png

Hmmm, there still doesn’t seem to be much going on (that channel down at the bottom looks noisy to me, rather than having a meaningful signal) so let’s transform this signal into something with a bit more SNR to it.

Extracting a more useful feature with xr.apply_ufunc

Without going into too much details on the neuroscience, iEEG data is particularly useful because there is information about neuronal activity in the higher frequency parts of the signal (AKA, parts of the electrical signal that change very quickly, but have very low amplitude). To pull that out, we’ll do the following:

  • High-pass filter the signal, which will remove all the slow-moving components

  • Calculate the envelope of the signal, which will tell us the power of high-frequency activity over time.

High-pass filtering the signal

MNE has a lot of nice functions for filtering a timeseries. Most of these operate on numpy arrays instead of MNE objects. We’ll use xarray’s apply_ufunc function to simply map that function onto our dataarray. xarray should keep track of the metadata (e.g. coordinates etc) and output a new DataArray with updated values.

flow = 80
fhigh = 140
da_lowpass = xr.apply_ufunc(
    mne.filter.filter_data, da,
   kwargs=dict(
       sfreq=da.sfreq,
       l_freq=flow,
       h_freq=fhigh,
   )
)
Setting up band-pass filter from 80 - 1.4e+02 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 80.00
- Lower transition bandwidth: 20.00 Hz (-6 dB cutoff frequency: 70.00 Hz)
- Upper passband edge: 140.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 145.00 Hz)
- Filter length: 99 samples (0.330 sec)

Visualizing our data, we can see all the slower fluctuations (e.g. long arcs over time) are gone.

fig, ax = plt.subplots(figsize=(15, 5))
da_lowpass.mean(dim='epoch').plot.line(x='time')
ax.get_legend().remove()
../../../_images/2019-10-22-xarray-neuro_25_0.png

Calculate the envelope of this signal with da.groupby

Next, we’ll calculate the envelope of the high-pass-filtered data. This is roughly the power that is present in these high frequencies over time. We do so by using something called a hilbert transform.

MNE also has a function for applying Hilbert transforms to data, but it has a weird quirk that expects the data to be of a particular shape. We can work around this by using our DataArray’s groupby method. This works similar to DataFrame.groupby - we’ll iterate through each channel, which will return a DataArray with shape epochs x timepoints. We can then calculate the Hilbert transform in each and re-combine into the original shape.

Note: This can be an expensive operation depending on the number of channels/epochs and the length of each trial. This might be a good place to insert paralellization via Dask.

def hilbert_2d(array):
    """Perform a Hilbert transforms on an (n_channels, n_times) array."""
    for ii, channel in enumerate(array):
        array[ii] = mne.filter._my_hilbert(channel, envelope=True)
    return array

da_hf_power = da_lowpass.groupby(da.coords['epoch']).apply(hilbert_2d)

The output dataarray should be the exact same shape, because we haven’t done any dimensional reductions. If we take a look at the resulting data, we can see what seems to be more structure in there:

fig, ax = plt.subplots(figsize=(15, 5))
da_hf_power.mean(dim='epoch').plot.line(x='time', hue='channel')
ax.get_legend().remove()
../../../_images/2019-10-22-xarray-neuro_30_0.png

Cleaning up our HFA data

Next let’s clean up this high-frequency activity (HFA) data.

Z-scoring our array

Instead of simple de-meaning the data like before, we’ll re-scale our data using the same baseline timepoints. What we’d like to do is the following:

  • Calculate the mean and standard deviation across trials of all pre-baseline data values, per channel

  • Z-score each channel using this mean and standard deviation

Once again we’ll use the groupby / apply combination to apply our function to subsets of the data.

# For each channel, apply a z-score that uses the mean/std of pre-baseline activity for all trials
def z_score(activity):
    """Take a DataArray and apply a z-score using the baseline"""
    baseline = activity.where(activity.time < -.1 )
    return (activity - np.nanmean(baseline)) / np.nanstd(baseline)

da_hf_zscored = da_hf_power.groupby('channel').apply(z_score)

Taking a look at the result, we can see a much cleaner separation of activity for some of the channels after time==0.

fig, ax = plt.subplots(figsize=(15, 5))
da_hf_zscored.mean(dim='epoch').plot.line(x='time', hue='channel')
ax.get_legend().remove()
../../../_images/2019-10-22-xarray-neuro_35_0.png

Smoothing our HFA data

Finally, let’s smooth this HFA so it has less jitter to it, and pick a smaller window that removes some of the filtering artifacts at the edges.

We’ll use the same filter_data function as before, but this time applied with the .groupby and .apply combination to show two ways of accomplishing the same thing. We’ll also use .sel to pick a subset of time for visualization

da_hf_zscored_lowpass = da_hf_zscored.groupby('epoch').apply(
    mne.filter.filter_data,
    sfreq=da.sfreq,
    l_freq=None,
    h_freq=10,
    verbose=False
)

Note that quickly selecting a subset of timepoints if we used numpy is much more verbose. Here’s a quick comparison:

# Numpy alone
mask_time = (times > -.8) * (times < 2.8)
epoch_dim = 0
da_hf_zscored_lowpass[..., mask_time].mean(epoch_dim)

# xarray
da_hf_zscored_lowpass.sel(time=slice(-.8, 2.8)).mean(dim='epoch')
fig, ax = plt.subplots(figsize=(15, 5))
da_hf_zscored_lowpass.mean(dim='epoch').sel(time=slice(-.8, 2.8)).plot.line(x='time', hue='channel')
ax.get_legend().remove()
../../../_images/2019-10-22-xarray-neuro_39_0.png

Now we can see there are clearly some channels that become active just after t==0. We can reduce our dataarray to a single dimension of “mean post-baseline activity in each channel” and convert it to a DataFrame for further processing:

# Find the channel with the most activity by first converting to a dataframe
total_activity = da_hf_zscored_lowpass.sel(time=slice(0, 2)).mean(dim=['epoch', 'time'])
total_activity = total_activity.to_dataframe()
total_activity.head()
Sample dataset
channel
ch_0 0.011101
ch_1 0.012152
ch_2 0.051911
ch_3 0.176188
ch_4 0.246021

Let’s grab the channel with maximal activation to look into a bit further.

max_chan = total_activity.squeeze().sort_values(ascending=False).index[0]

Time frequency analysis

As a final step, let’s expand our DataArray and add another dimension. In the above steps we specifically focused on high-frequency activity. A more common approach is to first create a spectrogram of your data to see activity across many frequencies.

To do this, we’ll use another MNE function for creating a Time-Frequency Representation or TFR. We’ll define a range of frequencies, and apply MNE’s function directly on our DataArray. This will return a NumPy array with the filtered values.

frequencies = [2**ii for ii in np.arange(2, 9, .5)]
tfr = mne.time_frequency.tfr_array_morlet(
    da,
    sfreq=da.sfreq,
    freqs=frequencies,
    n_cycles=4,
)

# Take the absolute value to throw out the non-real parts of the numbers
tfr = np.abs(tfr)
tfr[:2, :2, :2, :2]
array([[[[160.04103045, 160.38413909],
         [171.73704543, 175.09249553]],

        [[283.28699104, 285.21855726],
         [241.65630528, 245.77098295]]],


       [[[ 93.99546124,  94.4406537 ],
         [ 78.02050045,  79.15324341]],

        [[ 47.85148993,  49.90845151],
         [ 53.97221461,  54.55793674]]]])

Convert this data into a DataArray with .expand_dims

Next, we’ll convert this into a DataArray by using the metadata from our original DataArray. We can use the expand_dims method to create a new dimension for our DataArray. We’ll use this to store frequency information.

We’ll then reshape our new DataArray so that it matches the output of the MNE function, and use the copy method to create a new DataArray. By supplying the data= argument to copy, we directly insert the new data inside the generated DataArray.

da_tfr = (da
    .expand_dims(frequency=frequencies)
    .transpose('epoch', 'channel', 'frequency', 'time')
    .copy(data=np.log(tfr))
)

We can now visualize this time-frequency representation over time

fig, ax = plt.subplots(figsize=(15, 5))
(da_tfr
    .sel({'frequency': slice(None, 180), 'channel': max_chan})
    .mean('epoch')
    .plot.imshow(x='time', y='frequency')
)
<matplotlib.image.AxesImage at 0x7feb35626c88>
../../../_images/2019-10-22-xarray-neuro_50_1.png

Similar to our one-dimensional visualizations above, it can be hard to visualize relative changes in activity over a baseline (particularly because the amplitude scales inversely with the frequency).

Let’s apply a re-scaling function to our data so that we can see things more clearly. This time we’ll use MNE’s rescale function, which acts similarly to our zscore function above.

da_tfr_baselined = xr.apply_ufunc(
    mne.baseline.rescale,
    da_tfr,
    kwargs={'times': da_tfr.coords['time'], 'baseline': (None, -.1), "mode": 'zscore'}
)
Applying baseline correction (mode: zscore)

again, the result should be a DataArray, so we can directly visualize it:

(da_tfr_baselined
    .sel({'frequency': slice(None, 180), 'channel': max_chan, 'time': slice(-.8, 2.5)})
    .mean('epoch')
    .plot.imshow(x='time', y='frequency')
)
<matplotlib.image.AxesImage at 0x7feb06faaa58>
../../../_images/2019-10-22-xarray-neuro_54_1.png

Now we can see a clear increase in activity in the higher frequencies at t==0.

Combining the two with xr.merge

Finally, let’s combine these two DataArrays into one. We know that they share much of the same metadata - the first is “Amplitude of High-Frequency Activity” and the second is “Time-frequency power”. We should be able to merge these into a single xarray DataSet, which will allow us to perform operations across both by using their shared dimensions. DataSets are kind of like collections of DataArrays, with assumptions that the DataArrays share some metadata or coordinates.

First, we’ll rename each DataArray so that we can merge them nicely. Then, we’ll simply use the xr.merge function, which tries to automatically figure out which dimensions are shared based on their names and coordinate values.

da_tfr_baselined.name = "Time Frequency Representation"
da_hf_zscored_lowpass.name = "Low-pass filtered HFA"
ds = xr.merge([da_tfr_baselined, da_hf_zscored_lowpass])
ds
Show/Hide data repr Show/Hide attributes
xarray.Dataset
    • channel: 32
    • epoch: 29
    • frequency: 14
    • time: 1201
    • frequency
      (frequency)
      float64
      4.0 5.657 8.0 ... 181.0 256.0 362.0
      array([  4.      ,   5.656854,   8.      ,  11.313708,  16.      ,  22.627417,
              32.      ,  45.254834,  64.      ,  90.509668, 128.      , 181.019336,
             256.      , 362.038672])
    • time
      (time)
      float64
      -1.0 -0.9967 -0.9933 ... 2.997 3.0
      array([-1.      , -0.996667, -0.993333, ...,  2.993333,  2.996667,  3.      ])
    • channel
      (channel)
      'ch_0' 'ch_1' ... 'ch_30' 'ch_31'
      array(['ch_0', 'ch_1', 'ch_2', 'ch_3', 'ch_4', 'ch_5', 'ch_6', 'ch_7', 'ch_8',
             'ch_9', 'ch_10', 'ch_11', 'ch_12', 'ch_13', 'ch_14', 'ch_15', 'ch_16',
             'ch_17', 'ch_18', 'ch_19', 'ch_20', 'ch_21', 'ch_22', 'ch_23', 'ch_24',
             'ch_25', 'ch_26', 'ch_27', 'ch_28', 'ch_29', 'ch_30', 'ch_31'],
            dtype='<U5')
    • epoch
      (epoch)
      int64
      0 1 2 3 4 5 6 ... 23 24 25 26 27 28
      array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
             18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28])
    • Time Frequency Representation
      (epoch, channel, frequency, time)
      float64
      -1.068 -1.065 ... 1.245 1.231
      array([[[[-1.06825801e+00, -1.06490811e+00, -1.06269292e+00, ...,
                -2.44700037e-01, -2.42062452e-01, -2.40373002e-01],
               [-1.55585568e+00, -1.51452161e+00, -1.47437532e+00, ...,
                -5.35209555e-01, -5.40963440e-01, -5.48645559e-01],
               [-1.12040649e+00, -1.10631680e+00, -1.09442692e+00, ...,
                 6.84467283e-01,  6.80997403e-01,  6.74142959e-01],
               ...,
               [ 2.39450005e+00,  1.03380142e+00, -6.02721680e-01, ...,
                 2.49169178e-01,  1.43344119e+00,  3.01566173e+00],
               [ 4.65936619e-02,  9.72573317e-02, -9.93517803e-02, ...,
                 8.63789386e-01,  8.61999875e-01,  6.94728495e-01],
               [ 1.85858368e-01,  1.12174242e-01, -8.93810113e-02, ...,
                 8.28417458e-01,  8.34186151e-01,  7.73164525e-01]],
      
              [[-2.58676060e-01, -2.50290063e-01, -2.42548022e-01, ...,
                -9.04835872e-01, -9.34368529e-01, -9.63759221e-01],
               [-3.67942098e-01, -3.31923584e-01, -2.97381598e-01, ...,
                -2.12598677e+00, -2.13302253e+00, -2.13734082e+00],
               [ 2.60993077e-01,  2.89072254e-01,  3.14996787e-01, ...,
                 1.15761137e+00,  1.12048095e+00,  1.08173863e+00],
               ...,
               [ 2.77716710e+00,  6.04084887e-01, -9.68005734e-01, ...,
                 3.35246062e-01, -3.24181266e-02,  2.52724933e+00],
               [-1.90454378e-01, -7.01874341e-02, -1.15811020e-01, ...,
                 2.79311720e-01,  1.90165444e-01, -1.29436138e-01],
               [-1.88453765e-03, -4.46966730e-02, -7.80730531e-02, ...,
                 2.91926403e-01,  2.04639767e-01, -1.55335782e-04]],
      
              [[ 5.94408483e-01,  6.08399504e-01,  6.21475959e-01, ...,
                -4.03466096e+00, -4.04455269e+00, -4.05441933e+00],
               [ 5.66428039e-01,  5.75586680e-01,  5.83911402e-01, ...,
                -7.52782485e-01, -7.69198939e-01, -7.86330725e-01],
               [ 1.32613359e+00,  1.33765946e+00,  1.34502081e+00, ...,
                -7.96730973e-01, -8.60125890e-01, -9.25721374e-01],
               ...,
               [ 3.49323363e+00,  1.41513761e+00,  1.35878987e+00, ...,
                -1.16787577e+00,  5.34834275e-01,  2.08345745e+00],
               [ 9.31422956e-01,  1.10152722e+00,  1.18786301e+00, ...,
                -3.98972626e-01, -2.30411251e-01, -2.36744936e-01],
               [ 1.00765013e+00,  1.06353654e+00,  1.16189670e+00, ...,
                -3.65679246e-01, -2.15655383e-01, -7.86510000e-02]],
      
              ...,
      
              [[-2.53367490e+00, -2.60408304e+00, -2.68439840e+00, ...,
                -6.21454855e-01, -6.35155228e-01, -6.49024843e-01],
               [-1.75374717e+00, -1.77486397e+00, -1.79693935e+00, ...,
                -8.09528205e-01, -8.09433350e-01, -8.10154148e-01],
               [-1.81181076e+00, -1.78782909e+00, -1.76551910e+00, ...,
                -1.46261607e+00, -1.46822221e+00, -1.47758649e+00],
               ...,
               [ 2.07789036e+00,  7.18824827e-02,  1.06197216e+00, ...,
                -6.40794643e-01,  8.93274712e-01,  2.12543506e+00],
               [ 3.05552282e-01,  4.30699015e-01,  4.45286374e-01, ...,
                 2.82451076e-01,  4.04224777e-01,  3.22452827e-01],
               [ 3.99370920e-01,  4.21660940e-01,  4.78979168e-01, ...,
                 2.86166243e-01,  4.19348615e-01,  4.18523459e-01]],
      
              [[-3.80630126e-01, -3.73097727e-01, -3.66552271e-01, ...,
                -1.35468632e+00, -1.36291328e+00, -1.37195167e+00],
               [ 7.21233258e-01,  7.54431025e-01,  7.85343696e-01, ...,
                -1.09777363e+00, -1.11884845e+00, -1.14288773e+00],
               [ 7.22857120e-01,  7.34295663e-01,  7.41211983e-01, ...,
                -3.36577661e-01, -3.69971448e-01, -4.07883701e-01],
               ...,
               [ 3.54135397e+00,  1.61265218e+00,  8.98771010e-01, ...,
                -1.33149657e+00, -8.61629037e-01, -1.01506852e+00],
               [ 8.21206267e-01,  1.01975532e+00,  1.02555562e+00, ...,
                -1.10798717e+00, -1.67655694e+00, -2.66035015e+00],
               [ 9.16905323e-01,  9.82435051e-01,  9.91414455e-01, ...,
                -1.06591457e+00, -1.60015384e+00, -2.71945197e+00]],
      
              [[ 5.74308919e-01,  5.75633495e-01,  5.75495195e-01, ...,
                -3.28954363e-01, -3.40246647e-01, -3.52736858e-01],
               [ 4.34821963e-01,  4.51781613e-01,  4.67282265e-01, ...,
                -8.11740993e-01, -8.18559705e-01, -8.27546589e-01],
               [ 4.99543539e-01,  5.16121468e-01,  5.30156720e-01, ...,
                -2.14349713e-01, -2.51159176e-01, -2.90786618e-01],
               ...,
               [ 3.14354389e+00,  1.50349449e+00,  4.44072762e-01, ...,
                -1.42600314e+00, -8.67657124e-01, -2.00784876e-01],
               [ 9.42681585e-01,  1.06873600e+00,  1.07691673e+00, ...,
                -1.39523100e+00, -2.39979759e+00, -2.52772837e+00],
               [ 1.03394270e+00,  1.04123605e+00,  1.05185201e+00, ...,
                -1.40157882e+00, -2.69561869e+00, -2.10382962e+00]]],
      
      
             [[[-4.04899393e-01, -3.95691270e-01, -3.87413113e-01, ...,
                 1.35327771e+00,  1.31768441e+00,  1.28145280e+00],
               [-9.01461733e-01, -8.72102541e-01, -8.44249288e-01, ...,
                 5.61288275e-01,  5.16915364e-01,  4.71177072e-01],
               [-8.40763247e-01, -8.68534280e-01, -9.02714234e-01, ...,
                 6.52948281e-01,  6.16754552e-01,  5.78028736e-01],
               ...,
               [ 2.30438231e+00,  1.02358803e+00, -5.01482242e-01, ...,
                -1.45905587e+00,  4.19948990e-01,  1.58783492e+00],
               [-8.62313541e-01, -1.01489738e+00, -1.82117188e+00, ...,
                -3.20057966e+00, -2.21350440e+00, -1.79311988e+00],
               [-6.03629679e-01, -9.68351431e-01, -1.88888571e+00, ...,
                -3.28237325e+00, -2.23061491e+00, -1.44224544e+00]],
      
              [[-2.42497202e+00, -2.36642884e+00, -2.30813927e+00, ...,
                 7.13522894e-01,  6.87651119e-01,  6.61306991e-01],
               [-1.61389892e+00, -1.60113166e+00, -1.58760593e+00, ...,
                -2.11488435e-01, -2.43895836e-01, -2.77118695e-01],
               [-9.48637897e-01, -9.42168242e-01, -9.37381235e-01, ...,
                -1.00783090e+00, -1.07000890e+00, -1.13570400e+00],
               ...,
               [ 2.04690681e+00,  8.26658267e-01, -2.58985631e+00, ...,
                -8.20923674e-01, -1.28089044e+00, -1.34458379e+00],
               [-3.00265635e-01, -2.81324318e-01, -5.35245673e-01, ...,
                -2.00535062e+00, -2.77959030e+00, -3.47477990e+00],
               [-1.18554479e-01, -2.69175094e-01, -5.40118895e-01, ...,
                -1.96219684e+00, -2.87443956e+00, -3.34608953e+00]],
      
              [[-7.56859409e-01, -7.55186342e-01, -7.54150016e-01, ...,
                 1.80111302e-01,  1.62432491e-01,  1.44615986e-01],
               [-9.93338277e-01, -9.80677086e-01, -9.68946723e-01, ...,
                -5.66550486e-01, -5.63818933e-01, -5.63154950e-01],
               [-7.74875463e-01, -7.54171421e-01, -7.35150971e-01, ...,
                -2.52713394e-02, -5.37468287e-02, -8.35269325e-02],
               ...,
               [-6.23805966e-01, -1.05898761e+00, -1.47165407e-01, ...,
                -1.98760732e-01,  4.71743148e-01,  1.93527334e+00],
               [-1.62664974e+00, -1.09409273e+00, -8.06065742e-01, ...,
                -3.98544793e-01, -2.10299265e-01, -2.69129751e-01],
               [-1.60505429e+00, -1.07276822e+00, -7.41386145e-01, ...,
                -3.72584737e-01, -1.90235804e-01, -1.37385445e-01]],
      
              ...,
      
              [[-1.09258841e+00, -1.09255268e+00, -1.09324192e+00, ...,
                 1.63927751e-01,  1.38948305e-01,  1.13586373e-01],
               [-1.29821000e+00, -1.28808497e+00, -1.27885989e+00, ...,
                -9.25671704e-01, -9.83648334e-01, -1.04248370e+00],
               [-9.47032587e-01, -9.39314544e-01, -9.33676450e-01, ...,
                -1.35806667e+00, -1.38429806e+00, -1.41309924e+00],
               ...,
               [ 1.37617293e+00, -2.58149345e-01, -4.06315920e-01, ...,
                -2.19198008e-01,  7.61713551e-01,  1.86446923e+00],
               [-2.10935242e-01, -7.29780387e-02, -3.29670895e-02, ...,
                -4.63308853e-01, -1.17989823e-02,  3.29306567e-02],
               [-9.42372123e-02, -6.00808908e-02, -8.07446142e-03, ...,
                -4.73259795e-01,  1.74576890e-02,  1.77407193e-01]],
      
              [[-2.75341623e-01, -2.59113205e-01, -2.43883351e-01, ...,
                 1.26670422e+00,  1.23597671e+00,  1.20471485e+00],
               [-1.47693145e-01, -1.32036875e-01, -1.17427816e-01, ...,
                 7.89866923e-01,  7.80897492e-01,  7.70718152e-01],
               [-7.43276974e-01, -7.13754609e-01, -6.86061947e-01, ...,
                 1.98092423e+00,  1.94527804e+00,  1.90698660e+00],
               ...,
               [-2.09579206e+00, -2.91582027e+00, -1.32048206e+00, ...,
                 4.83749129e-01,  1.61605063e+00,  3.50200142e+00],
               [-2.93105467e+00, -2.21623749e+00, -1.62649697e+00, ...,
                 6.53002780e-01,  7.07502821e-01,  5.74271021e-01],
               [-2.88915091e+00, -2.21537964e+00, -1.60243628e+00, ...,
                 6.46074426e-01,  6.97008195e-01,  6.90917464e-01]],
      
              [[ 9.37896374e-02,  1.15599964e-01,  1.35909427e-01, ...,
                -6.08552783e-01, -6.53630489e-01, -6.99579003e-01],
               [ 1.17046531e+00,  1.22323098e+00,  1.27302607e+00, ...,
                 2.14994629e-01,  1.55372772e-01,  9.29064933e-02],
               [ 7.44233206e-01,  7.70336913e-01,  7.94178828e-01, ...,
                 4.34394855e-01,  4.01451571e-01,  3.66402272e-01],
               ...,
               [ 3.93677792e+00,  1.54015658e+00, -7.56258750e-02, ...,
                -6.41937683e-01,  8.24722248e-01,  3.25418288e+00],
               [ 2.64769665e-01,  4.44586061e-01,  3.39150385e-01, ...,
                -3.04163049e-01, -2.08308626e-01, -3.58056681e-01],
               [ 4.35588787e-01,  4.36937818e-01,  3.36376778e-01, ...,
                -2.74425085e-01, -1.96012045e-01, -1.48244191e-01]]],
      
      
             [[[ 9.40293628e-01,  9.42178452e-01,  9.43509028e-01, ...,
                -1.74752778e+00, -1.76930970e+00, -1.79123152e+00],
               [ 1.35548257e+00,  1.36497398e+00,  1.37385063e+00, ...,
                -1.30976739e+00, -1.33571183e+00, -1.36269150e+00],
               [ 1.75536584e+00,  1.76306862e+00,  1.76878200e+00, ...,
                -1.61556642e+00, -1.69610410e+00, -1.77850137e+00],
               ...,
               [ 3.33898149e+00,  1.42526661e+00,  1.33135102e+00, ...,
                -4.65665975e-01, -7.86327076e-01, -8.28409706e-01],
               [ 1.03087907e+00,  1.30772070e+00,  1.41761954e+00, ...,
                -2.83163882e+00, -3.78134319e+00, -4.84026420e+00],
               [ 1.12348661e+00,  1.27996919e+00,  1.39373302e+00, ...,
                -2.64299846e+00, -3.99663750e+00, -4.24665583e+00]],
      
              [[ 1.58106036e+00,  1.58887331e+00,  1.59621632e+00, ...,
                -8.17278034e-01, -7.94927632e-01, -7.73962011e-01],
               [ 1.66329760e+00,  1.67836808e+00,  1.69239449e+00, ...,
                -5.63162267e-01, -5.58362885e-01, -5.54746062e-01],
               [ 1.48164972e+00,  1.49608867e+00,  1.50809559e+00, ...,
                 4.21535100e-02,  3.34795387e-02,  2.24498264e-02],
               ...,
               [ 3.56895131e+00,  1.63695141e+00,  1.38799070e+00, ...,
                 1.13895322e-01, -6.19546308e-01,  1.79983028e+00],
               [ 1.12192676e+00,  1.32310053e+00,  1.37926076e+00, ...,
                 8.83071865e-02, -8.02014740e-02, -3.67938553e-01],
               [ 1.20434952e+00,  1.27965381e+00,  1.34491852e+00, ...,
                 1.02972232e-01, -7.66114927e-02, -2.50605663e-01]],
      
              [[ 1.39599929e+00,  1.40862370e+00,  1.42076325e+00, ...,
                -1.67586520e-01, -1.73665897e-01, -1.80038386e-01],
               [ 1.32253719e+00,  1.34384444e+00,  1.36387166e+00, ...,
                -1.57573737e+00, -1.54024795e+00, -1.50807094e+00],
               [ 9.58075018e-01,  9.74641801e-01,  9.88807879e-01, ...,
                -6.96927663e-01, -6.89851843e-01, -6.85165125e-01],
               ...,
               [ 3.59637117e+00,  2.01853809e+00,  6.25036437e-01, ...,
                 6.67947323e-01,  6.52591337e-01,  2.59941463e+00],
               [ 1.04562263e+00,  1.12069870e+00,  9.93794992e-01, ...,
                 4.98037027e-01,  4.95936472e-01,  3.19154869e-01],
               [ 1.14050851e+00,  1.08169194e+00,  9.62888419e-01, ...,
                 5.04336014e-01,  4.87487126e-01,  4.25098987e-01]],
      
              ...,
      
              [[ 1.95483056e-01,  2.01455987e-01,  2.06770047e-01, ...,
                -2.95885395e-01, -3.14194143e-01, -3.32773608e-01],
               [ 5.37089624e-01,  5.45194645e-01,  5.52419182e-01, ...,
                 1.80516277e-01,  1.64910528e-01,  1.48673354e-01],
               [ 8.71209480e-01,  8.89607737e-01,  9.06121442e-01, ...,
                 2.05449378e-02, -1.02256492e-02, -4.17304993e-02],
               ...,
               [ 2.88968819e+00,  1.16499073e+00, -1.00236320e+00, ...,
                -1.08990051e+00, -4.76584528e-01,  1.35160672e+00],
               [-9.88485951e-02, -3.74073896e-02, -2.26745670e-01, ...,
                -9.28263859e-01, -8.96539461e-01, -1.13983463e+00],
               [ 6.05317570e-02, -3.23972319e-02, -2.23188673e-01, ...,
                -8.93963755e-01, -8.54545139e-01, -1.00813872e+00]],
      
              [[ 1.64859177e-01,  1.72338638e-01,  1.79177398e-01, ...,
                -1.97759152e+00, -2.00187182e+00, -2.02640457e+00],
               [ 1.08112540e+00,  1.09105001e+00,  1.09991246e+00, ...,
                -3.86797049e+00, -4.00869203e+00, -4.06643733e+00],
               [ 1.79816868e+00,  1.80905175e+00,  1.81831680e+00, ...,
                -5.42584218e-01, -5.37243854e-01, -5.32308112e-01],
               ...,
               [ 4.71954963e+00,  2.75359110e+00,  6.12894250e-01, ...,
                -5.95455957e-01, -6.28969519e-01,  1.49910859e+00],
               [ 1.08454967e+00,  1.15733030e+00,  1.00518011e+00, ...,
                -4.94074351e-01, -6.17681380e-01, -9.92059905e-01],
               [ 1.21278822e+00,  1.12362043e+00,  9.77000989e-01, ...,
                -4.61788668e-01, -5.80241792e-01, -8.85285890e-01]],
      
              [[ 5.48948056e-01,  5.57660978e-01,  5.65621430e-01, ...,
                -4.81122847e-01, -4.81631140e-01, -4.82836899e-01],
               [ 1.09145917e+00,  1.10808801e+00,  1.12338111e+00, ...,
                -2.04917046e-01, -2.06945317e-01, -2.10274181e-01],
               [ 1.20390709e+00,  1.21774890e+00,  1.22980521e+00, ...,
                 1.93006638e-02,  4.19240932e-02,  6.14164214e-02],
               ...,
               [ 4.00792639e+00,  2.14197528e+00,  3.92633520e-01, ...,
                 5.18565236e-01,  1.15142696e+00,  3.32396568e+00],
               [ 7.97097959e-01,  8.64589293e-01,  6.40287194e-01, ...,
                 4.08759777e-01,  4.28817742e-01,  2.25847498e-01],
               [ 9.81963407e-01,  8.32319487e-01,  6.26185878e-01, ...,
                 4.12239785e-01,  4.11676406e-01,  3.98953891e-01]]],
      
      
             ...,
      
      
             [[[-1.52563214e+00, -1.52730731e+00, -1.53021655e+00, ...,
                -4.28417877e+00, -4.29282226e+00, -4.30281977e+00],
               [-1.52692799e+00, -1.48797385e+00, -1.45238784e+00, ...,
                -6.00174248e+00, -6.04704912e+00, -6.09601406e+00],
               [ 5.70433718e-01,  5.91651659e-01,  6.10391032e-01, ...,
                -5.60112685e-01, -5.89342204e-01, -6.21479783e-01],
               ...,
               [ 2.87602274e+00,  8.16411428e-01,  7.01284910e-01, ...,
                 2.43570053e-01,  3.71326993e-01,  2.33412913e+00],
               [ 4.91710530e-01,  7.18220351e-01,  8.07372432e-01, ...,
                 1.20417920e-01,  1.37509966e-01, -1.69702007e-02],
               [ 5.72769393e-01,  6.95907675e-01,  7.85014735e-01, ...,
                 1.48225680e-01,  1.47185829e-01,  1.10255843e-01]],
      
              [[-2.10810588e+00, -2.00240774e+00, -1.89948316e+00, ...,
                -6.64537476e+00, -6.72329757e+00, -6.80428157e+00],
               [-1.57206425e-01, -1.22270720e-01, -8.91170494e-02, ...,
                -1.95770765e+00, -1.96682205e+00, -1.97764568e+00],
               [ 6.76033321e-01,  7.06599704e-01,  7.34636992e-01, ...,
                -7.78958672e-02, -8.54461185e-02, -9.72182688e-02],
               ...,
               [ 2.89666383e+00,  1.38819054e+00,  1.79565827e-01, ...,
                -1.97082225e+00, -1.03177348e+00, -4.93003517e-01],
               [ 4.69972269e-01,  5.26564009e-01,  4.01128665e-01, ...,
                -2.68235465e+00, -3.45535289e+00, -3.05009916e+00],
               [ 5.88157757e-01,  5.17114187e-01,  4.02614928e-01, ...,
                -2.79086346e+00, -3.49076259e+00, -2.59604487e+00]],
      
              [[-3.78406495e-01, -3.62489628e-01, -3.47627457e-01, ...,
                -2.02469923e+00, -2.04840634e+00, -2.07241774e+00],
               [ 1.33232665e-01,  1.72693250e-01,  2.10568338e-01, ...,
                -3.70534718e+00, -3.63998427e+00, -3.58096687e+00],
               [ 6.96990521e-01,  7.31977627e-01,  7.64995473e-01, ...,
                -5.44724689e-01, -5.37134969e-01, -5.32540057e-01],
               ...,
               [ 3.00357479e+00,  1.45365798e+00,  2.00497517e-01, ...,
                 8.61693479e-01,  1.19740097e+00,  2.98849800e+00],
               [ 5.66590964e-01,  6.50473657e-01,  5.48229121e-01, ...,
                 7.95168991e-01,  7.97456681e-01,  6.36295094e-01],
               [ 6.80455943e-01,  6.44777964e-01,  5.45464208e-01, ...,
                 7.94435592e-01,  7.88864013e-01,  7.29683061e-01]],
      
              ...,
      
              [[-1.16536520e+00, -1.12294686e+00, -1.08147454e+00, ...,
                -5.26242112e-01, -5.37738216e-01, -5.50280906e-01],
               [-9.45649911e-01, -8.81376604e-01, -8.18939735e-01, ...,
                 1.54773483e-01,  1.34216653e-01,  1.12442334e-01],
               [-4.15972338e-01, -3.52334215e-01, -2.90655201e-01, ...,
                 3.60201203e-01,  3.26579155e-01,  2.89624800e-01],
               ...,
               [ 2.13139530e+00,  2.10058731e-01, -3.03750852e-01, ...,
                 1.52649368e-01,  3.39380291e-02,  1.96441353e+00],
               [-8.70767317e-02,  1.37259305e-01,  1.92521241e-01, ...,
                -1.30318598e-01, -6.65542359e-02, -2.49248599e-01],
               [ 1.55702889e-02,  1.49771544e-01,  2.00100103e-01, ...,
                -9.07729540e-02, -4.02212085e-02, -1.32073400e-01]],
      
              [[ 1.00409170e+00,  1.02755963e+00,  1.05038806e+00, ...,
                -1.80298758e+00, -1.78158172e+00, -1.76235800e+00],
               [ 1.06554383e+00,  1.08742881e+00,  1.10803481e+00, ...,
                -6.16587921e-01, -6.13392034e-01, -6.12184718e-01],
               [ 1.42348515e+00,  1.39086106e+00,  1.35184508e+00, ...,
                 1.04025220e+00,  1.04141805e+00,  1.03883446e+00],
               ...,
               [ 2.92501604e+00,  6.21244994e-01,  4.71588671e-01, ...,
                 2.50904972e-01, -1.05154979e+00,  1.88515410e+00],
               [ 1.34014321e-02,  2.81850326e-01,  3.62471090e-01, ...,
                -4.51596323e-01, -6.18002426e-01, -9.29349936e-01],
               [ 1.49898658e-01,  2.73009916e-01,  3.62738524e-01, ...,
                -4.03990647e-01, -6.08246582e-01, -7.63262540e-01]],
      
              [[ 4.16655190e-01,  4.31390394e-01,  4.45680603e-01, ...,
                -2.13714408e-01, -2.21400586e-01, -2.29679313e-01],
               [ 4.08799661e-01,  4.35294601e-01,  4.60747861e-01, ...,
                 6.14005784e-01,  5.98590520e-01,  5.81737216e-01],
               [ 2.68031006e-01,  2.81350320e-01,  2.93149208e-01, ...,
                 1.79852755e+00,  1.76909108e+00,  1.73776942e+00],
               ...,
               [-5.63676185e-01, -4.16378904e-01, -7.97361433e-01, ...,
                 7.28801112e-01,  7.33878165e-01,  2.74939552e+00],
               [-2.03707555e+00, -1.24146949e+00, -8.30854296e-01, ...,
                 4.31678238e-01,  3.99617924e-01,  2.08674915e-01],
               [-2.12273774e+00, -1.21052800e+00, -8.07879269e-01, ...,
                 4.40463583e-01,  3.90542616e-01,  3.25687301e-01]]],
      
      
             [[[ 1.57755132e+00,  1.59468807e+00,  1.61072052e+00, ...,
                -1.50734227e+00, -1.53837433e+00, -1.56991161e+00],
               [ 9.26168334e-01,  9.52899821e-01,  9.78389183e-01, ...,
                -1.85013596e+00, -1.86456367e+00, -1.88105053e+00],
               [ 8.32762041e-01,  8.41605966e-01,  8.47636737e-01, ...,
                -9.79470998e-01, -9.80247964e-01, -9.84740723e-01],
               ...,
               [ 2.28661151e+00,  7.95328426e-02,  8.46288206e-01, ...,
                -5.50613123e-01,  3.69386970e-01,  1.97840883e+00],
               [ 2.03213787e-01,  5.07349941e-01,  6.78492462e-01, ...,
                -3.41455937e-01, -2.56000391e-01, -3.18660029e-01],
               [ 2.60098025e-01,  5.00093387e-01,  6.73170027e-01, ...,
                -3.14765440e-01, -2.46838659e-01, -1.81191661e-01]],
      
              [[ 1.28174851e+00,  1.29208814e+00,  1.30192270e+00, ...,
                -6.06831024e-01, -6.08269893e-01, -6.10285919e-01],
               [ 1.17126823e+00,  1.20265272e+00,  1.23245689e+00, ...,
                -1.28554028e+00, -1.32425630e+00, -1.36492883e+00],
               [ 1.34378123e+00,  1.36437185e+00,  1.38058147e+00, ...,
                -8.95742142e-01, -9.52507519e-01, -1.01282565e+00],
               ...,
               [ 2.87435979e+00,  1.18201789e+00,  7.12622858e-01, ...,
                -1.00501004e+00, -7.89351367e-02,  1.17164139e+00],
               [ 8.04359758e-01,  9.93498178e-01,  1.04108831e+00, ...,
                -5.42316123e-01, -4.13511248e-01, -5.60298931e-01],
               [ 8.68940497e-01,  9.69795626e-01,  1.01406285e+00, ...,
                -4.93564990e-01, -3.51899688e-01, -4.28443512e-01]],
      
              [[ 7.67757581e-01,  7.98308994e-01,  8.27728850e-01, ...,
                -2.42094354e+00, -2.40685168e+00, -2.39537720e+00],
               [ 8.91788862e-01,  9.20617701e-01,  9.48121640e-01, ...,
                -4.26287403e-01, -4.51829029e-01, -4.78754520e-01],
               [ 1.63064381e+00,  1.65305814e+00,  1.67249287e+00, ...,
                 4.39571645e-01,  4.08194239e-01,  3.74027646e-01],
               ...,
               [ 2.91344495e+00,  1.19707852e+00,  1.06439859e+00, ...,
                -6.53158360e-01,  6.54483844e-01,  2.06264870e+00],
               [ 7.67589220e-01,  9.49722740e-01,  1.00933641e+00, ...,
                 2.31716108e-01,  2.52540666e-01,  1.24950681e-01],
               [ 8.22274073e-01,  9.37407145e-01,  9.96273242e-01, ...,
                 2.34813550e-01,  2.65973343e-01,  2.11321726e-01]],
      
              ...,
      
              [[ 9.15252914e-01,  9.26335536e-01,  9.36453182e-01, ...,
                 5.94340496e-01,  5.46423498e-01,  4.97816515e-01],
               [ 2.10456661e+00,  2.11821529e+00,  2.13040730e+00, ...,
                 1.65883799e-01,  1.17778684e-01,  6.89634802e-02],
               [ 2.11940366e+00,  2.12928510e+00,  2.13724339e+00, ...,
                -7.06930241e-01, -6.83296883e-01, -6.59632742e-01],
               ...,
               [ 3.47123261e+00,  1.34207044e+00,  1.78140718e+00, ...,
                -1.40069885e-01,  1.06023083e+00,  2.15567505e+00],
               [ 1.23956108e+00,  1.58128079e+00,  1.76614725e+00, ...,
                 1.87669777e-01,  2.63729116e-01, -1.13478408e-02],
               [ 1.28008510e+00,  1.53128806e+00,  1.71362703e+00, ...,
                 1.96006773e-01,  3.01486942e-01,  7.98313942e-02]],
      
              [[ 1.43915669e+00,  1.44731666e+00,  1.45490393e+00, ...,
                 1.04628046e+00,  1.02232380e+00,  9.97895568e-01],
               [ 1.61035451e+00,  1.61981921e+00,  1.62843733e+00, ...,
                 6.24705938e-01,  6.00904405e-01,  5.76322905e-01],
               [ 2.24063829e+00,  2.25311777e+00,  2.26354841e+00, ...,
                 5.61937968e-01,  5.26890950e-01,  4.89437031e-01],
               ...,
               [ 3.63259185e+00,  1.80243317e+00,  1.64489574e+00, ...,
                 1.90693826e-01,  1.39050420e+00,  2.62133085e+00],
               [ 1.23754931e+00,  1.45404970e+00,  1.53981939e+00, ...,
                -5.13626394e-01,  1.15506199e-02,  1.09980123e-01],
               [ 1.27919683e+00,  1.40334738e+00,  1.48814943e+00, ...,
                -5.22085498e-01,  1.62671168e-02,  2.68459039e-01]],
      
              [[ 4.18319048e-01,  4.36920036e-01,  4.54870305e-01, ...,
                 6.89887033e-01,  6.34788601e-01,  5.79161341e-01],
               [ 7.19790781e-01,  7.45435099e-01,  7.70062395e-01, ...,
                 4.45726861e-01,  4.32110473e-01,  4.19831363e-01],
               [ 1.42494597e+00,  1.46565940e+00,  1.50396874e+00, ...,
                 2.12708998e+00,  2.10800847e+00,  2.08715108e+00],
               ...,
               [ 4.10470708e+00,  2.22786753e+00,  2.81484625e-01, ...,
                 2.32900509e+00,  1.29457295e+00,  4.33247656e+00],
               [ 6.41157726e-01,  6.59836605e-01,  3.33337210e-01, ...,
                 1.69622802e+00,  1.52132068e+00,  1.17537167e+00],
               [ 8.24203461e-01,  6.40035111e-01,  3.26896228e-01, ...,
                 1.66668943e+00,  1.46800026e+00,  1.26693894e+00]]],
      
      
             [[[-1.63657916e+00, -1.58410547e+00, -1.53236309e+00, ...,
                 2.15569295e+00,  2.11859040e+00,  2.08056710e+00],
               [-1.55700938e+00, -1.47462991e+00, -1.39550347e+00, ...,
                 2.01259690e+00,  1.97091752e+00,  1.92820392e+00],
               [-1.06284921e+00, -1.00620788e+00, -9.49007811e-01, ...,
                 7.42397410e-01,  6.99437185e-01,  6.55541900e-01],
               ...,
               [ 3.65706131e-01,  7.82396067e-01, -2.38733996e+00, ...,
                 3.71731568e-01,  1.59149131e+00,  3.07004838e+00],
               [-1.73283273e+00, -8.29979364e-01, -5.82783121e-01, ...,
                -4.18810795e-01,  1.54487693e-01,  2.48719879e-01],
               [-1.85176944e+00, -7.05620960e-01, -5.47547538e-01, ...,
                -4.08443957e-01,  1.55826767e-01,  4.41785015e-01]],
      
              [[ 1.07968491e+00,  1.07751849e+00,  1.07417423e+00, ...,
                 2.80855507e+00,  2.76514355e+00,  2.72098253e+00],
               [ 1.07426739e+00,  1.07180616e+00,  1.06722482e+00, ...,
                 1.44608688e+00,  1.38790605e+00,  1.32937172e+00],
               [ 1.11151038e+00,  1.11544067e+00,  1.11682140e+00, ...,
                -4.82949216e-01, -4.48125210e-01, -4.20080706e-01],
               ...,
               [ 3.43740412e+00,  1.32478630e+00,  1.80658666e+00, ...,
                -4.18487471e-02,  1.33946617e+00,  2.89497076e+00],
               [ 1.26153567e+00,  1.46295163e+00,  1.54749684e+00, ...,
                 8.65661265e-01,  8.89714693e-01,  7.78560167e-01],
               [ 1.29262007e+00,  1.41120737e+00,  1.50292649e+00, ...,
                 8.43375228e-01,  8.70604803e-01,  8.55271824e-01]],
      
              [[-1.61814118e+00, -1.55163933e+00, -1.48601455e+00, ...,
                 2.98041695e+00,  2.93983419e+00,  2.89851200e+00],
               [-2.03554714e+00, -1.95398808e+00, -1.87138678e+00, ...,
                 1.46228732e+00,  1.42366269e+00,  1.38517905e+00],
               [-2.64833564e+00, -2.60434398e+00, -2.56238083e+00, ...,
                 1.35493304e+00,  1.36731273e+00,  1.37148300e+00],
               ...,
               [ 4.92982509e-01,  7.35894811e-01,  9.05699846e-01, ...,
                 1.42472860e+00,  2.41823800e+00,  4.38602641e+00],
               [-9.99214258e-01, -2.61635513e-01,  5.46600403e-03, ...,
                 1.43544588e+00,  1.49015066e+00,  1.37030473e+00],
               [-1.06728771e+00, -2.24630018e-01,  5.86759961e-02, ...,
                 1.41225414e+00,  1.45338060e+00,  1.48155691e+00]],
      
              ...,
      
              [[ 2.75594341e-01,  2.82115436e-01,  2.87988355e-01, ...,
                 1.22338508e+00,  1.19722147e+00,  1.17047157e+00],
               [ 3.19367429e-01,  3.46142206e-01,  3.70329884e-01, ...,
                 2.04315352e+00,  1.94863798e+00,  1.85195677e+00],
               [ 5.02235890e-01,  5.08674417e-01,  5.12345305e-01, ...,
                 8.17523109e-01,  7.73060712e-01,  7.26607934e-01],
               ...,
               [ 2.71818771e+00,  1.18613074e+00, -1.43370308e-01, ...,
                -1.60539066e-01,  5.97446182e-01,  2.06865155e+00],
               [ 5.63072684e-01,  6.62232442e-01,  6.80109620e-01, ...,
                -1.86215591e-01,  4.87331830e-02,  3.62484796e-02],
               [ 6.62841375e-01,  6.51329537e-01,  6.65326970e-01, ...,
                -1.55838563e-01,  6.46226566e-02,  1.68399524e-01]],
      
              [[-1.22737014e+00, -1.18883647e+00, -1.15087073e+00, ...,
                 1.23394915e+00,  1.16986001e+00,  1.10466912e+00],
               [-2.00318494e+00, -1.91176706e+00, -1.82212604e+00, ...,
                 3.27647804e-01,  1.71291068e-01,  1.21784785e-02],
               [-2.38394137e+00, -2.29133110e+00, -2.19901656e+00, ...,
                 6.88779362e-01,  6.42517357e-01,  5.94923561e-01],
               ...,
               [ 1.40403486e+00, -3.93977704e+00,  4.74370010e-01, ...,
                 4.28614500e-01,  1.58514132e+00,  3.42973651e+00],
               [-7.33500467e-01, -3.51935295e-01, -9.60072789e-02, ...,
                 7.17555452e-01,  7.85889675e-01,  6.46641420e-01],
               [-6.26451654e-01, -3.48627750e-01, -6.08625383e-02, ...,
                 7.07719890e-01,  7.73596070e-01,  7.62896284e-01]],
      
              [[-1.77083043e+00, -1.73604970e+00, -1.70150793e+00, ...,
                 1.09062712e+00,  1.02517939e+00,  9.59279976e-01],
               [-1.10369910e+00, -1.10286004e+00, -1.10428199e+00, ...,
                -1.54013290e+00, -1.51831853e+00, -1.48521130e+00],
               [-1.18254422e-01, -9.53441185e-02, -7.46309479e-02, ...,
                 2.48581223e-01,  2.43972400e-01,  2.38042422e-01],
               ...,
               [ 2.28713369e+00,  9.84919358e-01, -1.01329497e+00, ...,
                 9.61620967e-01,  2.05887862e+00,  3.86512765e+00],
               [-4.17881594e-01, -5.55130719e-01, -1.21157965e+00, ...,
                 1.25964575e+00,  1.28305741e+00,  1.13842255e+00],
               [-2.06602502e-01, -5.27021695e-01, -1.23039970e+00, ...,
                 1.22521418e+00,  1.24536401e+00,  1.23116620e+00]]]])
    • Low-pass filtered HFA
      (epoch, channel, time)
      float64
      -0.5572 -0.3505 ... -0.7235 -0.7839
      array([[[-0.55722445, -0.35045045, -0.15205402, ...,  0.32071735,
                0.2653267 ,  0.20620309],
              [-0.78504933, -0.71175459, -0.63951729, ..., -0.48619012,
               -0.58628096, -0.68681096],
              [ 0.08141248, -0.02148098, -0.12156504, ..., -0.79634564,
               -0.8382418 , -0.87936113],
              ...,
              [ 1.27220764,  1.04497205,  0.8213403 , ..., -0.31990238,
               -0.34853621, -0.37733814],
              [-0.41907094, -0.41268774, -0.40452625, ..., -0.11338628,
               -0.1352598 , -0.15755266],
              [-0.92384011, -0.86125627, -0.79696591, ..., -0.56037478,
               -0.69240248, -0.82728707]],
      
             [[-0.31745895, -0.35113427, -0.38355809, ..., -0.97980386,
               -1.18953749, -1.4044814 ],
              [-0.52255477, -0.4747038 , -0.42760367, ..., -0.17815445,
               -0.29884801, -0.42237863],
              [-0.5143446 , -0.49939427, -0.48358077, ..., -0.3344933 ,
               -0.56288699, -0.79734178],
              ...,
              [-1.12970668, -0.98675379, -0.84590074, ..., -0.66108719,
               -0.69541692, -0.72963154],
              [-0.76118298, -0.69168421, -0.62380438, ..., -0.44930201,
               -0.53370862, -0.61952314],
              [-0.79698496, -0.69017411, -0.58591732, ..., -0.1928013 ,
               -0.26283561, -0.33512544]],
      
             [[-1.04352089, -0.92781334, -0.81636026, ..., -0.63552018,
               -0.84232751, -1.05556983],
              [-0.70175051, -0.62421523, -0.54808752, ..., -0.70939563,
               -0.70622902, -0.70183793],
              [-0.31527505, -0.30666109, -0.29802473, ..., -0.47679593,
               -0.46786845, -0.45817719],
              ...,
              [-0.75496034, -0.7750953 , -0.79411767, ..., -0.38863549,
               -0.48162525, -0.57616487],
              [-0.35194069, -0.34217873, -0.3328539 , ..., -0.21180356,
               -0.29411081, -0.3787403 ],
              [-0.55860719, -0.52110093, -0.48451795, ..., -0.37747255,
               -0.46366781, -0.55188862]],
      
             ...,
      
             [[-1.53592456, -1.38733416, -1.24067899, ..., -0.73779431,
               -0.90933315, -1.0798213 ],
              [-0.55448721, -0.53254897, -0.51068622, ..., -0.32623476,
               -0.45909002, -0.5949611 ],
              [-0.77226291, -0.70987597, -0.64935919, ..., -0.72200745,
               -0.79295072, -0.86522509],
              ...,
              [-1.00929101, -0.94942079, -0.8907659 , ..., -0.71859234,
               -0.86394619, -1.01397194],
              [-0.50642446, -0.45218583, -0.39939951, ...,  0.07863249,
                0.12609372,  0.17453339],
              [-0.85923643, -0.79593505, -0.73382414, ..., -0.27641309,
               -0.50481494, -0.7391722 ]],
      
             [[-1.46952946, -1.39135579, -1.31447809, ..., -0.91041666,
               -1.08504229, -1.26339741],
              [-0.4635759 , -0.47777889, -0.49091248, ..., -0.00884495,
               -0.03918553, -0.07153994],
              [-0.72275354, -0.66562973, -0.60925247, ..., -0.45559207,
               -0.57205713, -0.69159195],
              ...,
              [-0.81110225, -0.76979449, -0.72913841, ..., -0.01297166,
               -0.17204604, -0.33373046],
              [-0.67504801, -0.64950589, -0.62407296, ..., -0.44131178,
               -0.53747513, -0.63536607],
              [-0.61259686, -0.60635283, -0.59986444, ...,  0.39894872,
                0.39382701,  0.39072565]],
      
             [[-0.50172746, -0.50007591, -0.49844615, ..., -0.91406926,
               -1.12188931, -1.33385636],
              [-0.48818892, -0.44053809, -0.39169382, ..., -0.8046927 ,
               -0.82182853, -0.83531051],
              [-0.68459361, -0.57315787, -0.46320266, ..., -0.51630811,
               -0.52405482, -0.5300185 ],
              ...,
              [-1.02575486, -0.74979224, -0.47995183, ..., -0.30871513,
               -0.44949317, -0.59194219],
              [ 0.04431285,  0.0419622 ,  0.03912131, ..., -0.530455  ,
               -0.58464065, -0.63994083],
              [-0.52155074, -0.51028736, -0.49879508, ..., -0.66415828,
               -0.72354198, -0.78394542]]])

Since we’ve got a single dataset, we can grab subsets along each axis across both DataArrays at the same time. We’ll select a subset of channels, time, and frequency bands to visualize.

ds_plt = ds.sel({'channel': max_chan, 'frequency': slice(10, 150), 'time': slice(-.5, 2)})

Now, we’ll plot both the spectrogram and the HFA in the same Matplotlib figure. As you can see, these plots contain somewhat redundant information. The top plot tells us that there is a general increase in power for high-frequencies. The bottom plot gives us the average increase in power across the higher frequencies.

fig, (ax_tfr, ax_hfa) = plt.subplots(2, 1, figsize=(15, 10))
im = ds_plt['Time Frequency Representation'].mean('epoch').plot.imshow(x='time', y='frequency',
                                                                       ax=ax_tfr)

ds_plt['Low-pass filtered HFA'].mean('epoch').plot.line(x='time', ax=ax_hfa)
[<matplotlib.lines.Line2D at 0x7feb06f5a470>]
../../../_images/2019-10-22-xarray-neuro_61_1.png

Wrapping up

In all, I was pretty happy with what you can do using xarray’s DataArray structure. It’s pretty nice to be able to refer to axes by their names, and to make more intelligent selection / slicing operations using their coordinate values. Moreover, this post is just scratching the surface for how to use this information in a way that speeds up the exploration and analysis post.

For example, we might have sped-up some feature extraction steps by using a distributed processing framework like Dask in the operations above. Dask integrates nicely with xarray, and offers a lot of interesting opportunities to parallelize interactive computation. I’ll explore that in another blog post.

Finally - the goal of this post has largely been to learn a bit more about xarray. This means I might be totally mis-using functionality, or missing something that would have made the above process much easier. If anybody has tips or thoughts on the code above, please do reach out!

Comments

comments powered by Disqus