tutorial_mappers.ipynb

{
  "metadata": {
    "name": "Getting data in shape"
  }, 
  "nbformat": 3, 
  "nbformat_minor": 0, 
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "heading", 
          "level": 1, 
          "metadata": {}, 
          "source": [
            "Getting data in shape"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "In the tutorial part ", 
            "*chap_tutorial_datasets* we have discovered a magic\ningredient of datasets: a ", 
            "[Mapper](http://pymvpa.org/generated/mvpa2.mappers.base.Mapper.html#mvpa2-mappers-base-mapper). Mappers are\nprobably the most powerful concept in PyMVPA, and there is little one would do\nwithout them.\n\n", 
            "In general, a mapper is an algorithm that transforms data.\nThis transformation can be as simple as selecting a subset of data, or as\ncomplex as a multi-stage preprocessing pipeline. Some transformations are\nreversible, others are not. Some are simple one-step computations, others\nare iterative algorithms that have to be trained on data before they can be\nused. In PyMVPA, all these transformations are ", 
            "[mappers](http://pymvpa.org/generated/mvpa2.mappers.html#mvpa2-mappers)."
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n*Note*", 
            "\n\n", 
            "If you are an ", 
            "[MDP](http://mdp-toolkit.sourceforge.net/)-user you probably have realized the similarity of MDP's\nnodes and PyMVPA's mappers.", 
            "- - -\n"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Let's create a dummy dataset (5 samples, 12 features). This time we will use a\nnew method to create the dataset, the `dataset_wizard`. Here it is, fully\nequivalent to a regular constructor call (i.e.\n", 
            "[Dataset](http://pymvpa.org/generated/mvpa2.datasets.base.Dataset.html#mvpa2-datasets-base-dataset)), but we will shortly see some nice convenience\naspects."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "from mvpa2.tutorial_suite import *\n", 
            "ds = dataset_wizard(np.ones((5, 12)))\n", 
            "ds.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Some datasets (such as the ones ", 
            "[fmri_dataset()](http://pymvpa.org/generated/mvpa2.datasets.mri.fmri_dataset.html#mvpa2-datasets-mri-fmri-dataset) with a\nmask) contain mappers as a ", 
            "[dataset attribute](http://pymvpa.org/glossary.html#term-dataset-attribute) `.a.mapper`.\nHowever, not every dataset actually has\na mapper. For example, the simple one we have just created doesn't have any:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "'mapper' in ds.a"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now let's look at a very similar dataset that only differs in a tiny but\na very important detail:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "ds = dataset_wizard(np.ones((5, 4, 3)))\n", 
            "ds.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "'mapper' in ds.a"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print ds.a.mapper"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We see that the resulting dataset looks identical to the one above, but this\ntime it got created from a 3D samples array (i.e. five samples, where each is a\n4x3 matrix). Somehow this 3D array got transformed into a 2D samples array in\nthe dataset. This magic behavior is unveiled by observing that the dataset's\nmapper is a ", 
            "[FlattenMapper](http://pymvpa.org/generated/mvpa2.mappers.flatten.FlattenMapper.html#mvpa2-mappers-flatten-flattenmapper).\n\n", 
            "The purpose of this mapper is precisely what we have just observed: reshaping\ndata arrays into 2D. It does it by preserving the first axis (in PyMVPA\ndatasets this is the axis that separates the samples) and concatenates all\nother axis into the second one.\n\n", 
            "Since mappers represent particular transformations they can also be seen as a\nprotocol of what has been done. If we look at the dataset, we know that it had\nbeen flattened on the way from its origin to a samples array in a dataset. This\nfeature can become really useful, if the processing become more complex. Let's\nlook at a possible next step -- selecting a subset of interesting features:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "myfavs = [1, 2, 8, 10]\n", 
            "subds = ds[:, myfavs]\n", 
            "subds.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "'mapper' in subds.a"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print subds.a.mapper"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now the situation has changed: ", 
            "*two* new mappers appeared in the dataset -- a\n", 
            "[ChainMapper](http://pymvpa.org/generated/mvpa2.mappers.base.ChainMapper.html#mvpa2-mappers-base-chainmapper) and a\n", 
            "[StaticFeatureSelection](http://pymvpa.org/generated/mvpa2.featsel.base.StaticFeatureSelection.html#mvpa2-featsel-base-staticfeatureselection).  The latter describes (and\nactually performs) the slicing operation we just made, while the former\nencapsulates the two mappers into a processing pipeline.  We can see that the\nmapper chain represents the processing history of the dataset like a breadcrumb\ntrack.\n\n", 
            "As it has been mentioned, mappers not only can transform a single dataset, but\ncan be fed with other data (as long as it is compatible with the mapper)."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "fwdtest = np.arange(12).reshape(4,3)\n", 
            "print fwdtest"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "fmapped = subds.a.mapper.forward1(fwdtest)\n", 
            "fmapped.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print fmapped"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Although `subds` has less features than our input data, forward mapping\napplies the same transformation that had been done to the dataset itself also\nto our test 4x3 array. The procedure yields a feature vector of the same shape\nas the one in `subds`. By looking at the forward-mapped data, we can verify\nthat the correct features have been chosen."
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "Load real data"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We have pretty much all the pieces to start a first analysis.  We know how to\nload fMRI data from time series images, we know how to add and access\nattributes in a dataset, we know how to slice datasets, and we know that we can\nmanipulate datasets with mappers.\n\n", 
            "Now our goal is to combine all these little pieces into the code that produces\na dataset like the one used in the seminal work by ", 
            "\n\n> *Haxby et al. (2001)* -- a study where participants passively watched gray scale images of\neight object categories in a block-design experiment. From the raw BOLD time\nseries, of which we have the full 12 recording runs of the first subject, they\ncomputed:", 
            "A ", 
            "*pattern of activation* for each stimulus category in each half of the\ndata (split by odd vs. even runs; i.e. 16 samples), including the\nassociated ", 
            "[sample attribute](http://pymvpa.org/glossary.html#term-sample-attribute)s that are necessary to perform a\ncross-validated classification analysis of the data.", 
            "\n\n", 
            "We have already seen how fMRI data can be loaded from NIfTI images, but this\ntime we need more than just the EPI images. For a classification analysis we\nalso need to associate each sample with a corresponding experimental condition,\ni.e. a class label, also sometimes called ", 
            "[target](http://pymvpa.org/glossary.html#term-target) value.  Moreover, for\na cross-validation procedure we also need to partition the full dataset into,\npresumably, independent ", 
            "[chunk](http://pymvpa.org/glossary.html#term-chunk)s. Independence is critical to achieve\nan unbiased estimate of the generalization performance of a classifier, i.e.\nits accuracy in predicting the correct class label for new data, unseen during\ntraining. So, where do we get this information from?\n\n", 
            "Both, target values and chunks are defined by the design of the experiment.\nIn the simplest case the target value for an fMRI volume sample is the\nexperiment condition that has been present/active while the volume has been\nacquired. However, there are more complicated scenarios which we will look\nat later on. Chunks of independent data correspond to what fMRI volumes are\nassumed to be independent. The properties of the MRI acquisition process\ncause subsequently acquired volumes to be ", 
            "*very* similar, hence they cannot\nbe considered independent. Ideally, the experiment is split into several\nacquisition sessions, where the sessions define the corresponding data\nchunks.\n\n", 
            "There are many ways to import this information into PyMVPA. The most simple\none is to create a two-column text file that has the target value in the\nfirst column, and the chunk identifier in the second, with one line per\nvolume in the NIfTI image."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "data_path = os.path.join(tutorial_data_path, 'haxby2001')\n", 
            "attr_fname = os.path.join(data_path, 'sub001',\n                          'BOLD', 'task001_run001', 'attributes.txt')\n", 
            "attr = SampleAttributes(attr_fname)\n", 
            "len(attr.targets)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print np.unique(attr.targets)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "len(attr.chunks)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print np.unique(attr.chunks)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "[SampleAttributes](http://pymvpa.org/generated/mvpa2.misc.io.base.SampleAttributes.html#mvpa2-misc-io-base-sampleattributes) allows us to load this type of\nfile, and access its content. We got 121 labels and chunk values, one for each\nvolume. Moreover, we see that there are nine different conditions and all\nsamples are associated with the same chunk. The attributes file for a different\nscan/run would increment the chunk value.\n\n", 
            "Now we can load the fMRI data, as we have done before -- only loading voxels\ncorresponding to a mask of ventral temporal cortex, and assign the samples\nattributes to the dataset. ", 
            "[fmri_dataset()](http://pymvpa.org/generated/mvpa2.datasets.mri.fmri_dataset.html#mvpa2-datasets-mri-fmri-dataset) allows us to\npass them directly:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "bold_fname = os.path.join(data_path,\n                          'sub001', 'BOLD', 'task001_run001', 'bold.nii.gz')\n", 
            "mask_fname = os.path.join(tutorial_data_path, 'haxby2001',\n                          'sub001', 'masks', 'orig', 'vt.nii.gz')\n", 
            "fds = fmri_dataset(samples=bold_fname,\n                   targets=attr.targets, chunks=attr.chunks,\n                   mask=mask_fname)\n", 
            "fds.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print fds.sa"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "We got the dataset that we already know from the last part, but this time\nis also has information about chunks and targets."
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "More structure, less duplication of work"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Although one could craft individual attribute files for each fMRI scan, doing\nso would be suboptimal. Typically, stimulation is not synchronous with\nfMRI volume sampling rate, hence timing information would be lost. Moreover,\ninformation on stimulation, or experiment design in general, is most likely\navailable already in different form or shape.\n\n", 
            "To ease working with a broad range of datasets, PyMVPA comes with dedicated\nsupport for datasets following the specifications used by the ", 
            "[openfmri.org](http://www.openfmri.org)\ndata-sharing platform. These are simple guidelines for file name conventions\nand design specification that can be easily adopted for your own data."
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n**Exercise**", 
            "\n\n", 
            "The tutorial data you are working with is following the openfmri.org\nscheme. Open the dataset folder and inspect the structure and content\nof the files with meta data. Notice, that it is possible to run a standard\nanalysis using, for example, FSL's FEAT directly on this data in unmodified\nform."
          ]
        }, 
        {
          "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", 
            "Accessing such a dataset is done via a handler that simply needs to know\nwhere the dataset is stored on disk. This handler offers convenient access\nto basic information, such as the number of subjects, task descriptions,\nand other properties."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "dhandle = OpenFMRIDataset(data_path)\n", 
            "dhandle.get_subj_ids()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "dhandle.get_task_descriptions()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "More importantly, it supports access to information on experiment design:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "model = 1\n", 
            "subj = 1\n", 
            "run = 1\n", 
            "events = dhandle.get_bold_run_model(model, subj, run)\n", 
            "for ev in events[:2]:\n    print ev"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "As you can see, the stimulus design information is available in a list of\nstandard Python dictionaries for each event. This includes onset and duration\nof the stimulation, as well as literal condition labels, and task descriptions.\n\n", 
            "With a utility function it is straightforward to convert such an event list\ninto a sample attribute array like the one we have loaded from a file before.\n`events2sample_attr()` uses the sample acquisition time information stored in\nthe dataset's `time_coords` sample attribute to match stimulation events to\ndata samples."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "targets = events2sample_attr(events, fds.sa.time_coords,\n                             noinfolabel='rest', onset_shift=0.0)\n", 
            "print np.unique([attr.targets[i] == t for i, t in enumerate(targets)])"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print np.unique(attr.targets)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            " print len(fds), len(targets)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Note, that the conversion of stimulation events to attribute arrays is a rather\ncrude way of labeling fMRI data that only works well with block-design-like\nexperiments. We will see other approaches later in this tutorial.\n\n", 
            "In addition to experiment design information, the dataset handler also offers\nconvenient access to the actual BOLD fMRI data:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "task = 1\n", 
            "fds = dhandle.get_bold_run_dataset(subj, task, run, mask=mask_fname)\n", 
            "print fds"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "The method `get_bold_run_dataset()` works the same way as `fmri_dataset()`,\nwhich we have seen before, and also supports the same arguments. However,\ninstead of giving a custom filename, BOLD data is identified by subject, task,\nand acquisition run IDs."
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 3, 
          "metadata": {}, 
          "source": [
            "Multi-session data"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Many fMRI experiments involve multiple runs. Loading such data is best done\nin a loop. The following code snippet loads all available runs for the object\nviewing task from our example subject in the dataset."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "task = 1   # object viewing task\n", 
            "model = 1  # image stimulus category model\n", 
            "subj = 1\n", 
            "run_datasets = []\n", 
            "for run_id in dhandle.get_task_bold_run_ids(task)[subj]:\n    # load design info for this run\n    run_events = dhandle.get_bold_run_model(model, subj, run_id)\n    # load BOLD data for this run (with masking); add 0-based chunk ID\n    run_ds = dhandle.get_bold_run_dataset(subj, task, run_id,\n                                          chunks=run_id -1,\n                                          mask=mask_fname)\n    # convert event info into a sample attribute and assign as 'targets'\n    run_ds.sa['targets'] = events2sample_attr(\n                run_events, run_ds.sa.time_coords, noinfolabel='rest')\n    # additional time series preprocessing can go here\n    run_datasets.append(run_ds)\n", 
            "fds = vstack(run_datasets, a=0)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now it is a good time to obtain a ", 
            "[summary()](http://pymvpa.org/generated/mvpa2.datasets.miscfx.summary.html#mvpa2-datasets-miscfx-summary) overview\nof the dataset: basic statistics, balance in number of samples among targets\nper chunk, etc.:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print fds.summary()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "In ", 
            "*chap_tutorial_openfmri* you can take a look at an example on how the\nkind of data preparation described above can be perform in an even more compact\nway.\n\n", 
            "The next step is to extract the ", 
            "*patterns of activation* from the dataset that\nwe are interested in. But wait! We know that fMRI data is typically\ncontaminated with a lot of noise, or actually ", 
            "*information* that we are not\ninterested in. For example, there are temporal drifts in the data (the signal\ntends to increase when the scanner is warming up). We also know that the signal\nis not fully homogeneous throughout the brain.\n\n", 
            "All these artifacts carry a lot of variance that is (hopefully) unrelated\nto the experiment design, and we should try to remove it to present the\nclassifier with the cleanest signal possible. There are countless ways to\npre-process the data to try to achieve this goal. Some keywords are:\nhigh/low/band-pass filtering, de-spiking, motion-correcting, intensity\nnormalization, and so on. In this tutorial, we keep it simple. The data we\nhave just loaded is already motion corrected. For every experiment that is\nlonger than a few minutes, as in this case, temporal trend removal, or\n", 
            "[detrending](http://pymvpa.org/glossary.html#term-detrending), is crucial."
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "Basic preprocessing"
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 3, 
          "metadata": {}, 
          "source": [
            "Detrending"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "PyMVPA provides functionality to remove polynomial trends from the data (other\nmethods are available too), meaning that polynomials are fitted to the time\nseries and only what is not explained by them remains in the dataset. In the\ncase of linear detrending, this means fitting a straight line to the time\nseries of each voxel via linear regression and taking the residuals as the new\nfeature values. Detrending can be seen as a type of data transformation, hence\nin PyMVPA it is implemented as a mapper."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "detrender = PolyDetrendMapper(polyord=1, chunks_attr='chunks')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "What we have just created is a mapper that will perform chunk-wise linear\n(1st-order polynomial) detrending. Chunk-wise detrending is desirable,\nsince our data stems from 12 different runs, and the assumption of a\ncontinuous linear trend across all runs is not appropriate. The mapper is\ngoing to use the `chunks` attribute to identify the chunks in the\ndataset.\n\n", 
            "We have seen that we could simply forward-map our dataset with this mapper.\nHowever, if we want to have the mapper present in the datasets processing\nhistory breadcrumb track, we can use its\n", 
            "[get_mapped()](http://pymvpa.org/generated/mvpa2.datasets.base.Dataset.get_mapped.html#mvpa2-datasets-base-dataset-get-mapped) method. This method will cause\nthe dataset to map a shallow copy of itself with the given mapper, and\nreturn it. Let's try:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "detrended_fds = fds.get_mapped(detrender)\n", 
            "print detrended_fds.a.mapper"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "`detrended_fds` is easily identifiable as a dataset that has been\nflattened, sliced, and linearly detrended.\n\n", 
            "Note that detrending doesn't always have to be an explicit step. For example,\nif you plan on modelling your data with haemodynamic response functions in a\ngeneral linear model (like it is shown in ", 
            "*chap_tutorial_openfmri* with\nNiPy), polynomial detrending can be done simultaneously as part of the\nmodeling."
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 3, 
          "metadata": {}, 
          "source": [
            "Normalization"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "While this will hopefully have solved the problem of temporal drifts in the\ndata, we still have inhomogeneous voxel intensities that can be a problem for\nsome machine learning algorithms. For this tutorial, we are again following a\nsimple approach to address this issue, and perform a feature-wise, chunk-wise\nZ-scoring of the data.  This has many advantages. First, it is going to scale\nall features into approximately the same range, and also remove their mean.\nThe latter is quite important, since some classifiers are impaired when working\nwith data having large offsets.  However, we are not going to perform a very\nsimple Z-scoring removing the global mean, but use the ", 
            "*rest* condition samples\nof the dataset to estimate mean and standard deviation.  Scaling dataset\nfeatures using these parameters yields a score corresponding to the per\ntime-point voxel intensity difference from the ", 
            "*rest* average.\n\n", 
            "This type of data ", 
            "[normalization](http://pymvpa.org/glossary.html#term-normalization) is, you guessed it, also\nimplemented as a mapper:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "zscorer = ZScoreMapper(param_est=('targets', ['rest']))"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "This mapper configuration implements a ", 
            "[chunk](http://pymvpa.org/glossary.html#term-chunk)-wise (the default)\nZ-scoring, while estimating mean and standard deviation from samples targets\nwith 'rest' in the respective chunk of data.\n\n", 
            "Remember, all mappers return new datasets that only have copies of what has\nbeen modified. However, both detrending and Z-scoring have or will modify the\nsamples themselves. That means that the memory consumption will triple!  We\nwill have the original data, the detrended data, and the Z-scored data, but\ntypically we are only interested in the final processing stage. To reduce the\nmemory footprint, both mappers have siblings that perform the same processing,\nbut without copying the data. For ", 
            "[PolyDetrendMapper](http://pymvpa.org/generated/mvpa2.mappers.detrend.PolyDetrendMapper.html#mvpa2-mappers-detrend-polydetrendmapper)\nthis is ", 
            "[poly_detrend()](http://pymvpa.org/generated/mvpa2.mappers.detrend.poly_detrend.html#mvpa2-mappers-detrend-poly-detrend), and for\n", 
            "[ZScoreMapper](http://pymvpa.org/generated/mvpa2.mappers.zscore.ZScoreMapper.html#mvpa2-mappers-zscore-zscoremapper) this is ", 
            "[zscore()](http://pymvpa.org/generated/mvpa2.mappers.zscore.zscore.html#mvpa2-mappers-zscore-zscore).\nThe following call will do the same as the mapper we have created above, but\nusing less memory:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "zscore(detrended_fds, param_est=('targets', ['rest']))\n", 
            "fds = detrended_fds\n", 
            "print fds.a.mapper"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n**Exercise**", 
            "\n\n", 
            "Look at the ", 
            "*example_smellit* example. Using the techniques from\nthis example, explore the dataset we have just created and look at the\neffect of detrending and Z-scoring."
          ]
        }, 
        {
          "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", 
            "The resulting dataset is now both detrended and normalized. The information\nis nicely presented in the mapper. From this point on we have no use for\nthe samples of the ", 
            "*rest* category anymore, hence we remove them from the\ndataset:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "fds = fds[fds.sa.targets != 'rest']\n", 
            "print fds.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "heading", 
          "level": 3, 
          "metadata": {}, 
          "source": [
            "Computing ", 
            "*Patterns Of Activation*"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "The last preprocessing step is to compute the actual ", 
            "*patterns of activation*.\nIn the original study, Haxby and colleagues performed a GLM-analysis of odd vs.\neven runs of the data respectively and used the corresponding contrast\nstatistics (stimulus category vs. rest) as classifier input. In this tutorial,\nwe will use a much simpler shortcut and just compute ", 
            "*mean* samples per\ncondition for both odd and even run independently.\n\n", 
            "To achieve this, we first add a new sample attribute to assign a\ncorresponding label to each sample in the dataset that indicates which of\nboth run-types it belongs to:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "rnames = {0: 'even', 1: 'odd'}\n", 
            "fds.sa['runtype'] = [rnames[c % 2] for c in fds.sa.chunks]"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "The rest is trivial. For cases like this -- applying a function (i.e. mean)\nto a set of groups of samples (all combinations of stimulus category and\nrun-type) -- PyMVPA has ", 
            "[FxMapper](http://pymvpa.org/generated/mvpa2.mappers.fx.FxMapper.html#mvpa2-mappers-fx-fxmapper). it comes with a number\nof convenience functions. The one we need here is\n", 
            "[mean_group_sample()](http://pymvpa.org/generated/mvpa2.mappers.fx.mean_group_sample.html#mvpa2-mappers-fx-mean-group-sample). It takes a list of sample attributes,\ndetermines all possible combinations of its unique values, selects dataset\nsamples corresponding to these combinations, and averages them. Finally,\nsince this is also a mapper, a new dataset with mean samples is returned:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "averager = mean_group_sample(['targets', 'runtype'])\n", 
            "type(averager)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "fds = fds.get_mapped(averager)\n", 
            "fds.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print fds.sa.targets"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print fds.sa.chunks"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Here we go! We now have a fully-preprocessed dataset: masked, detrended, normalized,\nwith one sample per stimulus condition that is an average for odd and even runs\nrespectively. Now we could do some serious classification, and this will be\nshown in ", 
            "*chap_tutorial_classifiers*, but there is still an\nimportant aspect of mappers we have to look at first."
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "There and back again -- a Mapper's tale"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Let's take a look back at the simple datasets from the start of the tutorial\npart."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print ds"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print ds.a.mapper"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "A very important feature of mappers is that they allow to reverse a\ntransformation, if that is possible. In case of the simple dataset we can\nask the mapper to undo the flattening and to put our samples back into the\noriginal 3D shape."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "orig_data = ds.a.mapper.reverse(ds.samples)\n", 
            "orig_data.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "In interactive scripting sessions this is would be a relatively bulky command\nto type, although it might be quite frequently used. To make ones fingers\nsuffer less there is a little shortcut that does exactly the same:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "orig_data = ds.O\n", 
            "orig_data.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "It is important to realize that reverse-mapping not only works with a single\nmapper, but also with a ", 
            "[ChainMapper](http://pymvpa.org/generated/mvpa2.mappers.base.ChainMapper.html#mvpa2-mappers-base-chainmapper). Going back to our\ndemo dataset from the beginning we can see how it works:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print subds"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print subds.a.mapper"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "subds.nfeatures"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "revtest = np.arange(subds.nfeatures) + 10\n", 
            "print revtest"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "rmapped = subds.a.mapper.reverse1(revtest)\n", 
            "rmapped.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print rmapped"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Reverse mapping of a single sample (one-dimensional feature vector) through the\nmapper chain created a 4x3 array that corresponds to the dimensions of a sample\nin our original data space. Moreover, we see that each feature value is\nprecisely placed into the position that corresponds to the features selected\nin the previous dataset slicing operation.\n\n", 
            "But now let's look at our fMRI dataset again. Here the mapper chain is a little\nmore complex:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "print fds.a.mapper"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Initial flattening followed by mask, detrending, Z-scoring and finally\naveraging. We would reverse mapping do in this case? Let's test:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "fds.nfeatures"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "revtest = np.arange(100, 100 + fds.nfeatures)\n", 
            "rmapped = fds.a.mapper.reverse1(revtest)\n", 
            "rmapped.shape"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "What happens is exactly what we expect: The initial one-dimensional vector\nis passed backwards through the mapper chain. Reverting a group-based\naveraging doesn't make much sense for a single vector, hence it is ignored.\nSame happens for Z-Scoring and temporal detrending. However, for all\nremaining mappers the transformations are reverse. First unmasked, and\nthen reshaped into the original dimensionality -- the brain volume.\n\n", 
            "We can check that this is really the case by only reverse-mapping through\nthe first two mappers in the chain and compare the result:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "rmapped_partial = fds.a.mapper[:2].reverse1(revtest)\n", 
            "(rmapped == rmapped_partial).all()"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "In case you are wondering: The ", 
            "[ChainMapper](http://pymvpa.org/generated/mvpa2.mappers.base.ChainMapper.html#mvpa2-mappers-base-chainmapper) behaves\nlike a regular Python list. We have just selected the first two mappers in\nthe list as another ", 
            "[ChainMapper](http://pymvpa.org/generated/mvpa2.mappers.base.ChainMapper.html#mvpa2-mappers-base-chainmapper) and used that one for\nreverse-mapping."
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 3, 
          "metadata": {}, 
          "source": [
            "Back To NIfTI"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "One last interesting aspect in the context of reverse mapping: Whenever it\nis necessary to export data from PyMVPA, such as results, dataset mappers\nalso play a critical role. For example we can easily export the `revtest`\nvector into a NIfTI brain volume image. This is possible because the mapper\ncan put it back into 3D space, and because the dataset also stores\ninformation about the original source NIfTI image."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "'imghdr' in fds.a"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "PyMVPA offers ", 
            "[map2nifti()](http://pymvpa.org/generated/mvpa2.datasets.mri.map2nifti.html#mvpa2-datasets-mri-map2nifti), a function to combine these\ntwo things and convert any vector into the corresponding NIfTI image:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "nimg = map2nifti(fds, revtest)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "This image can now be stored as a file (e.g. `nimg.to_filename('mytest.nii.gz')`).\nIn this format it is now compatible with the vast majority of neuroimaging\nsoftware."
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "- - -\n**Exercise**", 
            "\n\n", 
            "Save the NIfTI image to some file, and use an MRI viewer to overlay it\non top of the anatomical image in the demo dataset. Does it match our\noriginal mask image of ventral temporal cortex?"
          ]
        }, 
        {
          "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", 
            "There are many more mappers in PyMVPA than we could cover in the tutorial\npart. Some more will be used in other parts, but even more can be found the\n", 
            "[mappers](http://pymvpa.org/generated/mvpa2.mappers.html#mvpa2-mappers) module. Even though they all implement different\ntransformations, they can all be used in the same way, and can all be\ncombined into a chain."
          ]
        }
      ], 
      "metadata": {}
    }
  ]
}