tutorial_datasets.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:89ef4f34436395fcdc91cdae31008717f879ee9ca42cf1ca17726c0bfdf2d5f3"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Dataset basics and concepts"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "A [Dataset](http://pymvpa.org/generated/mvpa2.datasets.base.Dataset.html#mvpa2-datasets-base-dataset) is the basic data container in PyMVPA. It\n",
      "serves as the primary form of data storage, but also as a common container for\n",
      "results returned by most algorithms. In this tutorial part we will take a look\n",
      "at what a dataset consists of, and how it works.\n",
      "\n",
      "Most datasets in PyMVPA are represented as a two-dimensional array, where the\n",
      "first axis is the [sample](http://pymvpa.org/glossary.html#term-sample)s axis, and the second axis represents the\n",
      "[feature](http://pymvpa.org/glossary.html#term-feature)s of the samples.  In the simplest case, a dataset only\n",
      "contains *data* that is a matrix of numerical values."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from mvpa2.tutorial_suite import *\n",
      "data = [[  1,  1, -1],\n",
      "        [  2,  0,  0],\n",
      "        [  3,  1,  1],\n",
      "        [  4,  0, -1]]\n",
      "ds = Dataset(data)\n",
      "ds.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/python2.7/dist-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
        "  \"This module will be removed in 0.20.\", DeprecationWarning)\n",
        "/usr/lib/python2.7/dist-packages/sklearn/grid_search.py:43: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. This module will be removed in 0.20.\n",
        "  DeprecationWarning)\n",
        "/usr/lib/python2.7/dist-packages/sklearn/learning_curve.py:23: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the functions are moved. This module will be removed in 0.20\n",
        "  DeprecationWarning)\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "(4, 3)"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "len(ds)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.nfeatures"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "In the above example, every row vector in the `data` matrix becomes an\n",
      "observation, a [sample](http://pymvpa.org/glossary.html#term-sample), in the dataset, and every column vector\n",
      "represents an individual variable, a [feature](http://pymvpa.org/glossary.html#term-feature). The concepts of samples\n",
      "and features are essential for a dataset, hence we take a closer look.\n",
      "\n",
      "The dataset assumes that the first axis of the data is to be used to define\n",
      "individual samples. If the dataset is created using a one-dimensional vector it will\n",
      "therefore have as many samples as elements in the vector, and only one feature."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "one_d = [ 0, 1, 2, 3 ]\n",
      "one_ds = Dataset(one_d)\n",
      "one_ds.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "On the other hand, if a dataset is created from multi-dimensional data, only its\n",
      "second axis represents the features"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "m_ds = Dataset(np.random.random((3, 4, 2, 3)))\n",
      "m_ds.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "m_ds.nfeatures"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "In this case we have a dataset with three samples and four features, where each\n",
      "feature is a 2x3 matrix. In case somebody is wondering now why not simply\n",
      "treat each value in the data array as its own feature (yielding 24 features) --\n",
      "stay tuned, as this is going to be of importance later on."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Attributes"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "What we have seen so far does not really warrant the use of a dataset over a\n",
      "plain array or a matrix with samples. However, in the MVPA context we often\n",
      "need to know more about each sample than just the value of its features.  For\n",
      "example, in order to train a supervised-learning algorithm to discriminate two\n",
      "classes of samples we need per-sample [target](http://pymvpa.org/glossary.html#term-target) values to label each\n",
      "sample with its respective class.  Such information can then be used in order\n",
      "to, for example, split a dataset into specific groups of samples.  For this\n",
      "type of auxiliary information a dataset can also contain collections of three\n",
      "types of [attribute](http://pymvpa.org/glossary.html#term-attribute)s: a [sample attribute](http://pymvpa.org/glossary.html#term-sample-attribute), a [feature](http://pymvpa.org/glossary.html#term-feature), and a [dataset attribute](http://pymvpa.org/glossary.html#term-dataset-attribute)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "For samples"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Each [sample](http://pymvpa.org/glossary.html#term-sample) in a dataset can have an arbitrary number of additional\n",
      "attributes. They are stored as vectors of the same length as the number of\n",
      "samples in a collection, and are accessible via the `sa` attribute. A\n",
      "collection is similar to a standard Python `dict`, and hence adding sample\n",
      "attributes works just like adding elements to a dictionary:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa['some_attr'] = [ 0., 1, 1, 3 ]\n",
      "ds.sa.keys()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "However, sample attributes are not directly stored as plain data, but for\n",
      "various reasons as a so-called [Collectable](http://pymvpa.org/generated/mvpa2.base.collections.Collectable.html#mvpa2-base-collections-collectable) that in\n",
      "turn embeds a NumPy array with the actual attribute:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "type(ds.sa['some_attr'])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa['some_attr'].value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "This \"complication\" is done to be able to extend attributes with additional\n",
      "functionality that is often needed and can offer a significant speed-up of\n",
      "processing. For example, sample attributes carry a list of their unique values.\n",
      "This list is only computed once (upon first request) and can subsequently be\n",
      "accessed directly without repeated and expensive searches:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa['some_attr'].unique"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "However, for most interactive uses of PyMVPA this type of access to attributes'\n",
      "`.value` is relatively cumbersome (too much typing), therefore collections\n",
      "support direct access by name:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa.some_attr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Another purpose of the sample attribute collection is to preserve data\n",
      "integrity, by disallowing improper attributes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa['invalid'] = 4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa['invalid'] = [ 1, 2, 3, 4, 5, 6 ]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "But other than basic plausibility checks, no further constraints on values of\n",
      "samples attributes exist. As long as the length of the attribute vector matches\n",
      "the number of samples in the dataset, and the attributes values can be stored\n",
      "in a NumPy array, any value is allowed. Consequently, it is even possible to\n",
      "have n-dimensional arrays, not just vectors, as attributes -- as long as their\n",
      "first axis matched the number of samples in a dataset. Moreover, it is\n",
      "perfectly possible and supported to store literal (non-numerical) attributes.\n",
      "It should also be noted that each attribute may have its own individual data\n",
      "type, hence it is possible to have literal and numeric attributes in the same\n",
      "dataset."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa['literal'] = ['one', 'two', 'three', 'four']\n",
      "sorted(ds.sa.keys())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for attr in ds.sa:\n",
      "   print \"%s: %s\" % (attr, ds.sa[attr].value.dtype.name)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "For features"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "[Feature attribute](http://pymvpa.org/glossary.html#term-feature-attribute)s are almost identical to [sample attribute](http://pymvpa.org/glossary.html#term-sample-attribute)s, the *only* difference is that instead of having one attribute value per\n",
      "sample, feature attributes have one value per (guess what? ...) *feature*.\n",
      "Moreover, they are stored in a separate collection in the dataset that is\n",
      "called `fa`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.nfeatures"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.fa['my_fav'] = [0, 1, 0]\n",
      "ds.fa['responsible'] = ['me', 'you', 'nobody']\n",
      "sorted(ds.fa.keys())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "For the entire dataset"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Lastly, there can be also attributes, not per-sample, or per-feature, but for\n",
      "the dataset as a whole: so called [dataset attribute](http://pymvpa.org/glossary.html#term-dataset-attribute)s.  Both assigning\n",
      "such attributes and accessing them later on work in exactly the same way as for\n",
      "the other two types of attributes, except that dataset attributes are stored in\n",
      "their own collection which is accessible via the `a` property of the dataset.\n",
      "However, in contrast to sample and feature attribute, no constraints on the\n",
      "type or size are imposed -- anything can be stored. Let's store a list with the\n",
      "names of all files in the current directory, just because we can:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from glob import glob\n",
      "ds.a['pointless'] = glob(\"*\")\n",
      "'setup.py' in ds.a.pointless"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Slicing, resampling, feature selection"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "At this point we can already construct a dataset from simple arrays and enrich\n",
      "it with an arbitrary number of additional attributes. But just having a dataset\n",
      "isn't enough. We often need to be able to select subsets of a dataset for\n",
      "further processing.\n",
      "\n",
      "Slicing a dataset (i.e. selecting specific subsets) is very similar to\n",
      "slicing a NumPy array. It actually works *almost* identically. A dataset\n",
      "supports Python's `slice` syntax, but also selection by boolean masks and\n",
      "indices. The following three slicing operations result in equivalent output\n",
      "datasets, by always selecting every other samples in the dataset:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds[::2].samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mask = np.array([True, False, True, False])\n",
      "ds[mask].samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds[[0, 2]].samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n",
      "**Exercise**\n",
      "\n",
      "\n",
      "\n",
      "Search the \n",
      "[NumPy documentation](http://docs.scipy.org/doc/) for the\n",
      "difference between \"basic slicing\" and \"advanced indexing\". The aspect of\n",
      "memory consumption, especially, applies to dataset slicing as well, and being\n",
      "aware of this fact might help to write more efficient analysis scripts. Which\n",
      "of the three slicing approaches above is the most memory-efficient?  Which of\n",
      "the three slicing approaches above might lead to unexpected side-effects if\n",
      "the output dataset gets modified?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# you can use this cell for this exercise"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n",
      "\n",
      "\n",
      "All three slicing-styles are equally applicable to the selection of feature\n",
      "subsets within a dataset. Remember, features are represented on the second axis\n",
      "of a dataset."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds[:, [1,2]].samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "By applying a selection by indices to the second axis, we can easily get\n",
      "the last two features of our example dataset. Please note that the `:` is supplied\n",
      "for the first axis slicing. This is the Python way to indicate *take everything\n",
      "along this axis*, thus including all samples.\n",
      "\n",
      "As you can guess, it is also possible to select subsets of samples and\n",
      "features at the same time."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subds = ds[[0,1], [0,2]]\n",
      "subds.samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "If you have prior experience with NumPy you might be confused now. What you\n",
      "might have expected is this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.samples[[0,1], [0,2]]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "The above code applies the same slicing directly to the NumPy array of\n",
      "`.samples`, and the result is fundamentally different. For NumPy arrays\n",
      "this style of slicing allows selection of specific elements by their indices on\n",
      "each axis of an array. For PyMVPA's datasets this mode is not very useful,\n",
      "instead we typically want to select rows and columns, i.e. samples and\n",
      "features given by their indices."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n",
      "**Exercise**\n",
      "\n",
      "\n",
      "\n",
      "Try to select samples [0,1] and features [0,2] simultaneously using\n",
      "dataset slicing.  Now apply the same slicing to the samples array itself\n",
      "(`ds.samples`) -- make sure that the result doesn't surprise you and find\n",
      "a pure NumPy way to achieve similar selection."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# you can use this cell for this exercise"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n",
      "\n",
      "\n",
      "One last interesting thing to look at, in the context of dataset slicing,\n",
      "are the attributes. What happens to them when a subset of samples and/or\n",
      "features is chosen? Our original dataset had both samples and feature attributes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print ds.sa.some_attr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print ds.fa.responsible"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Now let's look at what they became in the subset-dataset we previously\n",
      "created:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print subds.sa.some_attr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print subds.fa.responsible"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "We see that both attributes are still there and, moreover, also the\n",
      "corresponding subsets have been selected.  It makes it convenient to select\n",
      "subsets of the dataset matching specific values of sample or feature attributes,\n",
      "or both:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subds = ds[ds.sa.some_attr == 1., ds.fa.responsible == 'me']\n",
      "print subds.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "To simplify such selections based on the values of attributes, it is possible\n",
      "to specify the desired selection as a dictionary for either samples of features\n",
      "dimensions, where each key corresponds to an attribute name, and each value\n",
      "specifies a list of desired attribute values.  Specifying multiple keys for\n",
      "either dimension can be used to obtain the intersection of matching elements:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "subds = ds[{'some_attr': [1., 0.], 'literal': ['two']}, {'responsible': ['me', 'you']}]\n",
      "print subds.sa.some_attr, subds.sa.literal, subds.fa.responsible"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n",
      "**Exercise**\n",
      "\n",
      "\n",
      "\n",
      "Check the documentation of the \n",
      "[select()](http://pymvpa.org/generated/mvpa2.datasets.base.Dataset.select.html#mvpa2-datasets-base-dataset-select) method\n",
      "that can also be used to implement such a selection, but provides an\n",
      "additional argument `strict`.  Modify the example above to select\n",
      "non-existing elements via `[]`, and compare to the result to the output\n",
      "of `select()` with `strict=False`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# you can use this cell for this exercise"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Load fMRI data"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Enough theoretical foreplay -- let's look at a concrete example of loading an\n",
      "fMRI dataset. PyMVPA has several helper functions to load data from specialized\n",
      "formats, and the one for fMRI data is [fmri_dataset()](http://pymvpa.org/generated/mvpa2.datasets.mri.fmri_dataset.html#mvpa2-datasets-mri-fmri-dataset). The\n",
      "example dataset we are going to look at is a single subject from Haxby et al.\n",
      "(2001).  For more convenience and less typing, we have a short cut for the\n",
      "path of the directory with the fMRI data: `tutorial_data_path``.\n",
      "\n",
      "In the simplest case, we now let [fmri_dataset](http://pymvpa.org/generated/mvpa2.datasets.mri.fmri_dataset.html#mvpa2-datasets-mri-fmri-dataset) do its job,\n",
      "by just pointing it to the fMRI data file. The data is stored as a NIfTI file\n",
      "that has all volumes of one experiment concatenated into a single file."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bold_fname = os.path.join(tutorial_data_path, 'haxby2001', 'sub001',\n",
      "                          'BOLD', 'task001_run001', 'bold.nii.gz')\n",
      "ds = fmri_dataset(bold_fname)\n",
      "len(ds)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.nfeatures"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "We can notice two things. First -- *it worked!* Second, we obtained a\n",
      "two-dimensional dataset with 121 samples (these are volumes in the NIfTI file),\n",
      "and over 160k features (these are voxels in the volume). The voxels are\n",
      "represented as a one-dimensional vector, and it seems that they have lost their\n",
      "association with the 3D-voxel-space. However, this is not the case, as we will\n",
      "see later.  PyMVPA represents data in this simple format to make it compatible\n",
      "with a vast range of generic algorithms that expect data to be a simple matrix.\n",
      "\n",
      "We loaded all data from that NIfTI file, but usually we would be interested in\n",
      "a subset only, i.e. \"brain voxels\".  [fmri_dataset](http://pymvpa.org/generated/mvpa2.datasets.mri.fmri_dataset.html#mvpa2-datasets-mri-fmri-dataset) is\n",
      "capable of performing data masking. We just need to specify a mask image. Such\n",
      "a mask image is generated in pretty much any fMRI analysis pipeline -- may it\n",
      "be a full-brain mask computed during skull-stripping, or an activation map from\n",
      "a functional localizer. We are going to use the original GLM-based localizer\n",
      "mask of ventral temporal cortex from Haxby et al. (2001).  Let's reload the\n",
      "dataset:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mask_fname = os.path.join(tutorial_data_path, 'haxby2001', 'sub001',\n",
      "                          'masks', 'orig', 'vt.nii.gz')\n",
      "ds = fmri_dataset(bold_fname, mask=mask_fname)\n",
      "len(ds)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.nfeatures"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "As expected, we get the same number of samples, but now only 577 features\n",
      "-- voxels corresponding to non-zero elements in the mask image. Now, let's\n",
      "explore this dataset a little further."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n",
      "**Exercise**\n",
      "\n",
      "\n",
      "\n",
      "Explore the dataset attribute collections. What kind of information do they\n",
      "contain?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# you can use this cell for this exercise"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "- - -\n",
      "\n",
      "\n",
      "Besides samples, the dataset offers a number of attributes that enhance the\n",
      "data with information that is present in the NIfTI image file header.\n",
      "Each sample has information about its volume index in the time series and the\n",
      "actual acquisition time (relative to the beginning of the file). Moreover, the\n",
      "original voxel index (sometimes referred to as `ijk`) for each feature is\n",
      "available too.  Finally, the dataset also contains information about the\n",
      "dimensionality of the input volumes, voxel size, and any other NIfTI-specific\n",
      "information since it also includes a dump of the full NIfTI image header."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa.time_indices[:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.sa.time_coords[:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.fa.voxel_indices[:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.a.voxel_eldim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.a.voxel_dim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "'imghdr' in ds.a"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "In addition to all this information, the dataset also carries a key additional\n",
      "attribute: the *mapper*. A mapper is an important concept in PyMVPA, and hence\n",
      "has its own *tutorial chapter*."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print ds.a.mapper"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Having all these attributes being part of a dataset is often a useful thing to\n",
      "have, but in some cases (e.g. when it comes to efficiency, and/or very large\n",
      "datasets) one might want to have a leaner dataset with just the information\n",
      "that is really necessary. One way to achieve this, is to strip all unwanted\n",
      "attributes. The Dataset class' [AttrDataset.copy()](http://pymvpa.org/generated/mvpa2.base.dataset.AttrDataset.html#mvpa2.base.dataset.AttrDataset.copy)\n",
      "method can help with that."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "stripped = ds.copy(deep=False, sa=['time_coords'], fa=[], a=[])\n",
      "print stripped"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "We can see that all attributes besides `time_coords` have been filtered out.\n",
      "Setting the `deep` arguments to `False` causes the copy function to reuse\n",
      "the data from the source dataset to generate the new stripped one, without\n",
      "duplicating all data in memory -- meaning both datasets now share the sample\n",
      "data and any change done to `ds` will also affect `stripped`."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Intermediate storage"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Some data preprocessing can take a long time.  One would rather prevent having\n",
      "to do it over and over again, and instead just store the preprocessed data into\n",
      "a file for subsequent analyses. PyMVPA offers functionality to store a large\n",
      "variety of objects, including datasets, into [HDF5](http://en.wikipedia.org/wiki/Hierarchical_Data_Format) files. A variant of this\n",
      "format is also used by recent versions of Matlab to store data.\n",
      "\n",
      "For HDF5 support, PyMVPA depends on the [h5py](http://h5py.alfven.org) package. If it is available, any\n",
      "dataset can be saved to a file by simply calling\n",
      "[AttrDataset.save()](http://pymvpa.org/generated/mvpa2.base.dataset.AttrDataset.html#mvpa2.base.dataset.AttrDataset.save) with the desired filename."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import tempfile, shutil\n",
      "tempdir = tempfile.mkdtemp()\n",
      "ds.save(os.path.join(tempdir, 'mydataset.hdf5'))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "HDF5 is a flexible format that also supports, for example, data compression. To\n",
      "enable it, you can pass additional arguments to\n",
      "[AttrDataset.save()](http://pymvpa.org/generated/mvpa2.base.dataset.AttrDataset.html#mvpa2.base.dataset.AttrDataset.save) that are supported by h5py's\n",
      "`Group.create_dataset()`. Instead of using\n",
      "[AttrDataset.save()](http://pymvpa.org/generated/mvpa2.base.dataset.AttrDataset.html#mvpa2.base.dataset.AttrDataset.save) one can also use the\n",
      "[h5save()](http://pymvpa.org/generated/mvpa2.base.hdf5.h5save.html#mvpa2-base-hdf5-h5save) function in a similar way. Saving the same dataset\n",
      "with maximum gzip-compression looks like this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ds.save(os.path.join(tempdir, 'mydataset.gzipped.hdf5'), compression=9)\n",
      "h5save(os.path.join(tempdir, 'mydataset.gzipped.hdf5'), ds, compression=9)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Loading datasets from a file is easy too. [h5load()](http://pymvpa.org/generated/mvpa2.base.hdf5.h5load.html#mvpa2-base-hdf5-h5load) takes a\n",
      "filename as an argument and returns the stored dataset. Compressed data will be\n",
      "handled transparently."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "loaded = h5load(os.path.join(tempdir, 'mydataset.hdf5'))\n",
      "np.all(ds.samples == loaded.samples)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "shutil.rmtree(tempdir, ignore_errors=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "\n",
      "\n",
      "Note that this type of dataset storage is not appropriate from long-term\n",
      "archival of data, as it relies on a stable software environment. For long-term\n",
      "storage, use other formats."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}