Dataset basics and concepts — PyMVPA 2.4.0 documentation (2023)

Note

This tutorial part is also available for download as an IPython notebook:[ipynb]

Attributes¶

What we have seen so far does not really warrant the use of a dataset over aplain array or a matrix with samples. However, in the MVPA context we oftenneed to know more about each sample than just the value of its features. Forexample, in order to train a supervised-learning algorithm to discriminate twoclasses of samples we need per-sample target values to label eachsample with its respective class. Such information can then be used in orderto, for example, split a dataset into specific groups of samples. For thistype of auxiliary information a dataset can also contain collections of threetypes of attributes: a sample attribute, a featureattribute, and a dataset attribute.

>>> ds.sa['some_attr'] = [ 0., 1, 1, 3 ]>>> ds.sa.keys()['some_attr']

>>> type(ds.sa['some_attr'])<class 'mvpa2.base.collections.ArrayCollectable'>>>> ds.sa['some_attr'].valuearray([ 0., 1., 1., 3.])

>>> ds.sa['some_attr'].uniquearray([ 0., 1., 3.])

>>> ds.sa.some_attrarray([ 0., 1., 1., 3.])

>>> ds.sa['invalid'] = 4Traceback (most recent call last): File "/usr/lib/python2.6/doctest.py", line 1253, in __run compileflags, 1) in test.globs File "<doctest tutorial_datasets.rst[20]>", line 1, in <module> ds.sa['invalid'] = 4 File "/home/test/pymvpa/mvpa2/base/collections.py", line 459, in __setitem__ value = ArrayCollectable(value) File "/home/test/pymvpa/mvpa2/base/collections.py", line 171, in __init__ % self.__class__.__name__)ValueError: ArrayCollectable only takes sequences as value.Traceback (most recent call last): File "/usr/lib/python2.6/doctest.py", line 1253, in __run compileflags, 1) in test.globs File "<doctest tutorial_datasets.rst[20]>", line 1, in <module> ds.sa['invalid'] = 4 File "/home/test/pymvpa/mvpa2/base/collections.py", line 459, in __setitem__ value = ArrayCollectable(value) File "/home/test/pymvpa/mvpa2/base/collections.py", line 171, in __init__ % self.__class__.__name__)ValueError: ArrayCollectable only takes sequences as value.

>>> ds.sa['invalid'] = [ 1, 2, 3, 4, 5, 6 ]Traceback (most recent call last): File "/usr/lib/python2.6/doctest.py", line 1253, in __run compileflags, 1) in test.globs File "<doctest tutorial_datasets.rst[21]>", line 1, in <module> ds.sa['invalid'] = [ 1, 2, 3, 4, 5, 6 ] File "/home/test/pymvpa/mvpa2/base/collections.py", line 468, in __setitem__ str(self)))ValueError: Collectable 'invalid' with length [6] does not match the required length [4] of collection '<SampleAttributesCollection: some_attr>'.Traceback (most recent call last): File "/usr/lib/python2.6/doctest.py", line 1253, in __run compileflags, 1) in test.globs File "<doctest tutorial_datasets.rst[21]>", line 1, in <module> ds.sa['invalid'] = [ 1, 2, 3, 4, 5, 6 ] File "/home/test/pymvpa/mvpa2/base/collections.py", line 468, in __setitem__ str(self)))ValueError: Collectable 'invalid' with length [6] does not match the required length [4] of collection '<SampleAttributesCollection: some_attr>'.

>>> ds.sa['literal'] = ['one', 'two', 'three', 'four']>>> sorted(ds.sa.keys())['literal', 'some_attr']>>> for attr in ds.sa:...  print "%s: %s" % (attr, ds.sa[attr].value.dtype.name)literal: string40some_attr: float64

>>> ds.nfeatures3>>> ds.fa['my_fav'] = [0, 1, 0]>>> ds.fa['responsible'] = ['me', 'you', 'nobody']>>> sorted(ds.fa.keys())['my_fav', 'responsible']

>>> from glob import glob>>> ds.a['pointless'] = glob("*")>>> 'setup.py' in ds.a.pointlessTrue

Slicing, resampling, feature selection¶

At this point we can already construct a dataset from simple arrays and enrichit with an arbitrary number of additional attributes. But just having a datasetisn’t enough. We often need to be able to select subsets of a dataset forfurther processing.

Slicing a dataset (i.e. selecting specific subsets) is very similar toslicing a NumPy array. It actually works almost identically. A datasetsupports Python’s slice syntax, but also selection by boolean masks andindices. The following three slicing operations result in equivalent outputdatasets, by always selecting every other samples in the dataset:

>>> # original>>> ds.samplesarray([[ 1, 1, -1], [ 2, 0, 0], [ 3, 1, 1], [ 4, 0, -1]])>>>>>> # Python-style slicing>>> ds[::2].samplesarray([[ 1, 1, -1], [ 3, 1, 1]])>>>>>> # Boolean mask array>>> mask = np.array([True, False, True, False])>>> ds[mask].samplesarray([[ 1, 1, -1], [ 3, 1, 1]])>>>>>> # Slicing by index -- Python indexing start with 0 !!>>> ds[[0, 2]].samplesarray([[ 1, 1, -1], [ 3, 1, 1]])

All three slicing-styles are equally applicable to the selection of featuresubsets within a dataset. Remember, features are represented on the second axisof a dataset.

>>> ds[:, [1,2]].samplesarray([[ 1, -1], [ 0, 0], [ 1, 1], [ 0, -1]])

By applying a selection by indices to the second axis, we can easily getthe last two features of our example dataset. Please note that the : is suppliedfor the first axis slicing. This is the Python way to indicate take everythingalong this axis, thus including all samples.

As you can guess, it is also possible to select subsets of samples andfeatures at the same time.

>>> subds = ds[[0,1], [0,2]]>>> subds.samplesarray([[ 1, -1], [ 2, 0]])

If you have prior experience with NumPy you might be confused now. What youmight have expected is this:

>>> ds.samples[[0,1], [0,2]]array([1, 0])

The above code applies the same slicing directly to the NumPy array of.samples, and the result is fundamentally different. For NumPy arraysthis style of slicing allows selection of specific elements by their indices oneach axis of an array. For PyMVPA’s datasets this mode is not very useful,instead we typically want to select rows and columns, i.e. samples andfeatures given by their indices.

One last interesting thing to look at, in the context of dataset slicing,are the attributes. What happens to them when a subset of samples and/orfeatures is chosen? Our original dataset had both samples and feature attributes:

>>> print ds.sa.some_attr[ 0. 1. 1. 3.]>>> print ds.fa.responsible['me' 'you' 'nobody']

Now let’s look at what they became in the subset-dataset we previouslycreated:

>>> print subds.sa.some_attr[ 0. 1.]>>> print subds.fa.responsible['me' 'nobody']

We see that both attributes are still there and, moreover, also thecorresponding subsets have been selected. It makes it convenient to selectsubsets of the dataset matching specific values of sample or feature attributes,or both:

>>> subds = ds[ds.sa.some_attr == 1., ds.fa.responsible == 'me']>>> print subds.shape(2, 1)

To simplify such selections based on the values of attributes, it is possibleto specify the desired selection as a dictionary for either samples of featuresdimensions, where each key corresponds to an attribute name, and each valuespecifies a list of desired attribute values. Specifying multiple keys foreither dimension can be used to obtain the intersection of matching elements:

>>> subds = ds[{'some_attr': [1., 0.], 'literal': ['two']}, {'responsible': ['me', 'you']}]>>> print subds.sa.some_attr, subds.sa.literal, subds.fa.responsible[ 1.] ['two'] ['me' 'you']

Load fMRI data¶

Enough theoretical foreplay – let’s look at a concrete example of loading anfMRI dataset. PyMVPA has several helper functions to load data from specializedformats, and the one for fMRI data is fmri_dataset(). Theexample dataset we are going to look at is a single subject from Haxby et al.(2001). For more convenience and less typing, we have a short cut for thepath of the directory with the fMRI data: tutorial_data_path`.

In the simplest case, we now let fmri_dataset do its job,by just pointing it to the fMRI data file. The data is stored as a NIfTI filethat has all volumes of one experiment concatenated into a single file.

>>> bold_fname = os.path.join(tutorial_data_path, 'haxby2001', 'sub001',...  'BOLD', 'task001_run001', 'bold.nii.gz')>>> ds = fmri_dataset(bold_fname)>>> len(ds)121>>> ds.nfeatures163840>>> ds.shape(121, 163840)

We can notice two things. First – it worked! Second, we obtained atwo-dimensional dataset with 121 samples (these are volumes in the NIfTI file),and over 160k features (these are voxels in the volume). The voxels arerepresented as a one-dimensional vector, and it seems that they have lost theirassociation with the 3D-voxel-space. However, this is not the case, as we willsee later. PyMVPA represents data in this simple format to make it compatiblewith a vast range of generic algorithms that expect data to be a simple matrix.

We loaded all data from that NIfTI file, but usually we would be interested ina subset only, i.e. “brain voxels”. fmri_dataset iscapable of performing data masking. We just need to specify a mask image. Sucha mask image is generated in pretty much any fMRI analysis pipeline – may itbe a full-brain mask computed during skull-stripping, or an activation map froma functional localizer. We are going to use the original GLM-based localizermask of ventral temporal cortex from Haxby et al. (2001). Let’s reload thedataset:

>>> mask_fname = os.path.join(tutorial_data_path, 'haxby2001', 'sub001',...  'masks', 'orig', 'vt.nii.gz')>>> ds = fmri_dataset(bold_fname, mask=mask_fname)>>> len(ds)121>>> ds.nfeatures577

As expected, we get the same number of samples, but now only 577 features– voxels corresponding to non-zero elements in the mask image. Now, let’sexplore this dataset a little further.

Besides samples, the dataset offers a number of attributes that enhance thedata with information that is present in the NIfTI image file header.Each sample has information about its volume index in the time series and theactual acquisition time (relative to the beginning of the file). Moreover, theoriginal voxel index (sometimes referred to as ijk) for each feature isavailable too. Finally, the dataset also contains information about thedimensionality of the input volumes, voxel size, and any other NIfTI-specificinformation since it also includes a dump of the full NIfTI image header.

>>> ds.sa.time_indices[:5]array([0, 1, 2, 3, 4])>>> ds.sa.time_coords[:5]array([ 0. , 2.5, 5. , 7.5, 10. ])>>> ds.fa.voxel_indices[:5]array([[ 6, 23, 24], [ 7, 18, 25], [ 7, 18, 26], [ 7, 18, 27], [ 7, 19, 25]])>>> ds.a.voxel_eldim(3.5, 3.75, 3.75)>>> ds.a.voxel_dim(40, 64, 64)>>> 'imghdr' in ds.aTrue

In addition to all this information, the dataset also carries a key additionalattribute: the mapper. A mapper is an important concept in PyMVPA, and hencehas its own tutorial chapter.

>>> print ds.a.mapper<Chain: <Flatten>-<StaticFeatureSelection>>

Having all these attributes being part of a dataset is often a useful thing tohave, but in some cases (e.g. when it comes to efficiency, and/or very largedatasets) one might want to have a leaner dataset with just the informationthat is really necessary. One way to achieve this, is to strip all unwantedattributes. The Dataset class’ copy()method can help with that.

>>> stripped = ds.copy(deep=False, sa=['time_coords'], fa=[], a=[])>>> print stripped<Dataset: 121x577@int16, <sa: time_coords>>

We can see that all attributes besides time_coords have been filtered out.Setting the deep arguments to False causes the copy function to reusethe data from the source dataset to generate the new stripped one, withoutduplicating all data in memory – meaning both datasets now share the sampledata and any change done to ds will also affect stripped.

Intermediate storage¶

Some data preprocessing can take a long time. One would rather prevent havingto do it over and over again, and instead just store the preprocessed data intoa file for subsequent analyses. PyMVPA offers functionality to store a largevariety of objects, including datasets, into HDF5 files. A variant of thisformat is also used by recent versions of Matlab to store data.

For HDF5 support, PyMVPA depends on the h5py package. If it is available, anydataset can be saved to a file by simply callingsave() with the desired filename.

>>> import tempfile, shutil>>> # create a temporary directory>>> tempdir = tempfile.mkdtemp()>>> ds.save(os.path.join(tempdir, 'mydataset.hdf5'))

HDF5 is a flexible format that also supports, for example, data compression. Toenable it, you can pass additional arguments tosave() that are supported by h5py’sGroup.create_dataset(). Instead of usingsave() one can also use theh5save() function in a similar way. Saving the same datasetwith maximum gzip-compression looks like this:

>>> ds.save(os.path.join(tempdir, 'mydataset.gzipped.hdf5'), compression=9)>>> h5save(os.path.join(tempdir, 'mydataset.gzipped.hdf5'), ds, compression=9)

Loading datasets from a file is easy too. h5load() takes afilename as an argument and returns the stored dataset. Compressed data will behandled transparently.

>>> loaded = h5load(os.path.join(tempdir, 'mydataset.hdf5'))>>> np.all(ds.samples == loaded.samples)True>>> # cleanup the temporary directory, and everything it includes>>> shutil.rmtree(tempdir, ignore_errors=True)

Note that this type of dataset storage is not appropriate from long-termarchival of data, as it relies on a stable software environment. For long-termstorage, use other formats.