{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Working with ragged indices for multivariate connectivity\n\nThis example demonstrates how multivariate connectivity involving different\nnumbers of seeds and targets can be handled in MNE-Connectivity.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Author: Thomas S. Binns \n# License: BSD (3-clause)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n\nfrom mne_connectivity import spectral_connectivity_epochs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Background\n\nWith multivariate connectivity, interactions between multiple signals can be\nconsidered together, and the number of signals designated as seeds and\ntargets does not have to be equal within or across connections. Issues can\narise from this when storing information associated with connectivity in\narrays.\n\nSuch arrays are 'ragged', and support for ragged arrays is limited in NumPy\nto the ``object`` datatype. Not only is working with ragged arrays\ncumbersome, but saving arrays with ``dtype='object'`` is not supported by the\nh5netcdf engine used to save connectivity objects.\n\nThe workaround used in MNE-Connectivity is to pad ragged arrays with some\nknown values according to the largest number of entries in each dimension,\nsuch that there is an equal amount of information across and within\nconnections for each dimension of the arrays.\n\nAs an example, consider we have 5 channels and want to compute 2 connections:\nthe first between channels in indices 0 and 1 with those in indices 2, 3,\nand 4; and the second between channels 0, 1, 2, and 3 with channel 4. The\nseed and target indices can be written as such::\n\n seeds = [[0, 1 ], [0, 1, 2, 3]]\n targets = [[2, 3, 4], [4 ]]\n\nThe ``indices`` parameter passed to\n:func:`~mne_connectivity.spectral_connectivity_epochs` and\n:func:`~mne_connectivity.spectral_connectivity_time` must be a tuple of\narray-likes, meaning\nthat the indices can be passed as a tuple of: lists; tuples; or NumPy arrays.\nExamples of how ``indices`` can be formed are shown below::\n\n # tuple of lists\n ragged_indices = ([[0, 1 ], [0, 1, 2, 3]],\n [[2, 3, 4], [4 ]])\n\n # tuple of tuples\n ragged_indices = (((0, 1 ), (0, 1, 2, 3)),\n ((2, 3, 4), (4, )))\n\n # tuple of arrays\n ragged_indices = (np.array([[0, 1 ], [0, 1, 2, 3]], dtype='object'),\n np.array([[2, 3, 4], [4 ]], dtype='object'))\n\n**N.B. Note that since NumPy v1.19.0, dtype='object' must be specified when\nforming ragged arrays.**\n\nJust as for bivariate connectivity, the length of ``indices[0]`` and\n``indices[1]`` is equal (i.e. the number of connections), however information\nabout the multiple channel indices for each connection is stored in a nested\narray.\n\nImportantly, these indices are ragged, as the first connection will be\ncomputed between 2 seed and 3 target channels, and the second connection\nbetween 4 seed and 1 target channel(s). The connectivity functions will\nrecognise the indices as being ragged, and pad them to a 'full' array by\nadding placeholder values which are masked accordingly. This makes the\nindices easier to work with, and also compatible with the engine used to save\nconnectivity objects. For example, the above indices would become::\n\n padded_indices = (np.array([[0, 1, --, --], [0, 1, 2, 3]]),\n np.array([[2, 3, 4, --], [4, --, --, --]]))\n\nwhere ``--`` are masked entries. These indices are what is stored in the\nreturned connectivity objects.\n\nFor the connectivity results themselves, the methods available in\nMNE-Connectivity combine information across the different channels into a\nsingle (time-)frequency-resolved connectivity spectrum, regardless of the\nnumber of seed and target channels, so ragged arrays are not a concern here.\n\nHowever, the maximised imaginary part of coherency (MIC) method also returns\nspatial patterns of connectivity, which show the contribution of each channel\nto the dimensionality-reduced connectivity estimate (explained in more detail\nin :doc:`mic_mim`). Because these patterns are returned for each channel,\ntheir shape can vary depending on the number of seeds and targets in each\nconnection, making them ragged.\n\nTo avoid this, the patterns are padded along the channel axis with the known\nand invalid entry ``np.nan``, in line with that applied to ``indices``.\nExtracting only the valid spatial patterns from the connectivity object is\ntrivial, as shown below:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# create random data\ndata = np.random.randn(10, 5, 200) # epochs x channels x times\nsfreq = 50\nragged_indices = ([[0, 1], [0, 1, 2, 3]], [[2, 3, 4], [4]]) # seeds # targets\n\n# compute connectivity\ncon = spectral_connectivity_epochs(\n data,\n method=\"mic\",\n indices=ragged_indices,\n sfreq=sfreq,\n fmin=10,\n fmax=30,\n verbose=False,\n)\npatterns = np.array(con.attrs[\"patterns\"])\npadded_indices = con.indices\nn_freqs = con.get_data().shape[-1]\nn_cons = len(ragged_indices[0])\nmax_n_chans = max(len(inds) for inds in ([*ragged_indices[0], *ragged_indices[1]]))\n\n# show that the padded indices entries are masked\nassert np.sum(padded_indices[0][0].mask) == 2 # 2 padded channels\nassert np.sum(padded_indices[1][0].mask) == 1 # 1 padded channels\nassert np.sum(padded_indices[0][1].mask) == 0 # 0 padded channels\nassert np.sum(padded_indices[1][1].mask) == 3 # 3 padded channels\n\n# patterns have shape [seeds/targets x cons x max channels x freqs (x times)]\nassert patterns.shape == (2, n_cons, max_n_chans, n_freqs)\n\n# show that the padded patterns entries are all np.nan\nassert np.all(np.isnan(patterns[0, 0, 2:])) # 2 padded channels\nassert np.all(np.isnan(patterns[1, 0, 3:])) # 1 padded channels\nassert not np.any(np.isnan(patterns[0, 1])) # 0 padded channels\nassert np.all(np.isnan(patterns[1, 1, 1:])) # 3 padded channels\n\n# extract patterns for first connection using the ragged indices\nseed_patterns_con1 = patterns[0, 0, : len(ragged_indices[0][0])]\ntarget_patterns_con1 = patterns[1, 0, : len(ragged_indices[1][0])]\n\n# extract patterns for second connection using the padded, masked indices\nseed_patterns_con2 = patterns[0, 1, : padded_indices[0][1].count()]\ntarget_patterns_con2 = patterns[1, 1, : padded_indices[1][1].count()]\n\n# show that shapes of patterns are correct\nassert seed_patterns_con1.shape == (2, n_freqs) # channels (0, 1)\nassert target_patterns_con1.shape == (3, n_freqs) # channels (2, 3, 4)\nassert seed_patterns_con2.shape == (4, n_freqs) # channels (0, 1, 2, 3)\nassert target_patterns_con2.shape == (1, n_freqs) # channels (4)\n\nprint(\"Assertions completed successfully!\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}