{
 "metadata": {
  "name": "",
  "signature": "sha256:3f1c0208726fd26279967156a4dd185e6c335c42d198f2f115b7d61e38b4f4ed"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - compatibility with Python 3\n",
      "from __future__ import print_function  # print('me') instead of print 'me'\n",
      "from __future__ import division  # 1/2 == 0.5, not 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import common modules\n",
      "import numpy as np  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Using our module"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Later on we are going to use our new module `spm_funcs.py`.\n",
      "\n",
      "In order to do this, we need to make sure Python can find the module.  Python stores all the directories that it will search for modules in the `path` variable in the `sys` module.  We can make sure Python can find our new module by appending the directory containing the module to the `sys.path` list.\n",
      "\n",
      "Import the `sys` module, and show the `path` variable from the `sys` module."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Show path in sys module"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Find the directory where you put the `spm_funcs.py` module.  For example, for me, that would be `/Users/mb312/pna_code`.   Append this directory to the path list."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Append the directory containing `spm_funcs.py` to sys.path"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Check you can import spm_funcs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## The task time-course"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We've previously had a look at the file `ds114_sub009_t2r1.nii`.  This is a 4D FMRI image.\n",
      "\n",
      "Now we want to see whether we can detect any signal in that image relating to the task.\n",
      "\n",
      "The task was a block design with 10 seconds rest followed by 7 repeats of (30 seconds when the subject thought of verbs followed by 30 seconds rest).\n",
      "\n",
      "There's a file in this directory called `ds114_sub009_t2r1_cond.txt` with the task information in it.   It comes from subject 9 task 2 and run 1 from OpenFMRI dataset https://openfmri.org/dataset/ds000114 .\n",
      "\n",
      "The file has one line of text per \"on\" block, giving the onset time (in seconds), the duration (in seconds) and the amplitude (expected amount of activation for this block).\n",
      "\n",
      "Here are the first four lines:\n",
      "\n",
      "```\n",
      "10\t30.000000\t1\n",
      "70\t30.000000\t1\n",
      "130\t30.000000\t1\n",
      "190\t30.000000\t1\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Read the file into an array called \"task\".\n",
      "# \"task\" should have 3 columns (onset, duration, amplitude)\n",
      "# HINT: np.loadtxt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The TR for this FMRI run was 2.5 seconds.  We need to convert the onsets and durations to TRs - so for example the first onset was at 10 seconds, which was the start of TR 4 (10 / 2.5).\n",
      "\n",
      "Select out the first two columns of task, and divide by the TR to convert the onset and duration times to be in terms of TRs:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Select first two columns and divide by TR"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Our next step is to make an on-off vector that is 0 when the subject is doing nothing and 1 when the subject is doing the covert verb task.\n",
      "\n",
      "The vector will have one value (either 0 or 1) for each TR."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First use nibabel to load the image `ds114_sub009_t2r1.nii`.  Check the image shape to find the number of TRs."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load the image and check the image shape to get the number of TRs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next make a new vector called `time_course` with one entry per TR, with all elements in the vector being zero:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make new zero vector"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Loop over the rows in the onsets / durations array to give you an onset / duration pair.  For each of these pairs, set the matching positions in the `time_course` vector to 1.  For example, the first pair will be `4, 12`.   That means the task started at the beginning of scan index 4, and lasted for 12 scans.  That means there should be 12 consecutive 1 values in `time_course`, starting at index 4.  That means that index 4 + 12 = 16 should be zero again, because there are 12 values starting at (including 4) going up to (but *not* including) 16."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - try running this if you don't believe me\n",
      "len(range(4, 16))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So, for the first row, you will want to set `time_course[4]` through `time_course[15]` equal to 1."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Fill in values of 1 for positions of on blocks in time course"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot the time course:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot the time course"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Comparing task to rest"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make a boolean array `is_task_tr` which is True when `time_course` is 1 and False otherwise.\n",
      "\n",
      "Make another array `is_rest_tr` that is the opposite - True when `time_course` is 0 and False otherwise."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make two boolean arrays encoding task, rest"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now read the image data into an array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Read the image data into an array."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remember that the 4D array consists of one volume (3D array) per TR.\n",
      "\n",
      "We want to select the volumes where the time course is 1 (task volumes). Do this by slicing, using the boolean array you just made."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create a new 4D array only containing the task volumes"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now select the volumes where the time course is 0 (rest volumes):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create a new 4D array only containing the rest volumes"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We want to know whether there is a difference in signal in the task volumes compared to the rest volumes.  Take the mean over the task volumes and mean over the rest volumes.  You should end up with two 3D volumes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create the mean volume across all the task volumes\n",
      "# Then create the mean volume across all the rest volumes\n",
      "# Hint: remember the `axis` keyword."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now subtract the rest mean from the task mean to get a difference volume."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create a difference volume"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Show a slice over the third dimension of the difference volume, from somewhere around the center of the third axis:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Show a slice over the third dimension"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the difference between activation and rest.  It looks a little strange.  Maybe there are some artefacts here."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Fixing an artefact"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's use our `spm_globals` function to see if there are any volumes with odd values for this estimate of the mean."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Use the `get_spm_globals` function in `spm_funcs` to get the SPM global values for each volume in `ds114_sub009_t2r1.nii`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get SPM global values for each volume"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot the SPM global values for each volume\n",
      "# Try using 'o' marker symbol to make the outlier easier to see"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There is one volume that is a large outlier in the SPM globals.\n",
      "\n",
      "Slice out the bad volume (to a 3D array) and the following volume.\n",
      "\n",
      "Use `plt.subplots` to show a matching slice from the bad volume and the following volume side-by-side:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Select bad volume and following volume\n",
      "# Use subplots to show slice from good and next"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In fact, the bad volume is one of the rest volumes.\n",
      "\n",
      "Use your slicing skills to remove this volume from your selection of rest volumes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use slicing to remove outlier volume from rest volumes"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make a new mean for the rest volumes, and subtract this mean from the mean for the task volumes to make a new difference image.\n",
      "\n",
      "Give the new difference image a new name, so we can compare to the old difference image later."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make new mean for rest volumes, subtract from task mean"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Use subplot to show the same slice from the old difference volume and the new difference volume, side by side:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use subplot to show same slice from old and new difference volume"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## If you are racing ahead"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can actually do this calculation of the mean image in another way, by doing an elementwise array multiplication, followed by a sum.\n",
      "\n",
      "See if you can work out how to do that.  This section is difficult, so don't be surprised if you have to stare at it for a long time, or you can't get it easily.  The exercise is to practice slicing and thinking about 4D arrays.\n",
      "\n",
      "Imagine that we are looking at a single voxel time series, for 3D voxel coordinate $i, j, k$. \n",
      "\n",
      "Call the time series vector $y_0, y_1 ... y_{T-1}$ where \n",
      "\n",
      "* $y_0 = image\\_data[i, j, k, 0]$,\n",
      "* $y_1 = image\\_data[i, j, k, 1]$\n",
      "* $...$\n",
      "* $y_{T-1} = image\\_data[i, j, k, T-1]$\n",
      "\n",
      "$T = 173$ - the number of TRs.\n",
      "\n",
      "Call the indices of the rest volumes sequence $R \\in 0, 1, 2, 3, 16 ... $ of length $N$, and the indices of the task volumes $K \\in 4, 5, 6 ...$ of length $M$.\n",
      "\n",
      "The average of the points from $y$ within the rest volumes is:\n",
      "\n",
      "$$\n",
      "\\frac{1}{N} \\sum_r{y_r} : r \\in R\n",
      "$$\n",
      "\n",
      "The average of the points from $y$ within the task volumes is:\n",
      "\n",
      "$$\n",
      "\\frac{1}{M} \\sum_k{y_k} : k \\in K\n",
      "$$\n",
      "\n",
      "Subtracting rest from task to get the difference in means at this voxel:\n",
      "\n",
      "$$\n",
      "\\frac{1}{M} \\sum_k{y_k} - \\frac{1}{N} \\sum_r{y_r}\n",
      "$$\n",
      "\n",
      "\n",
      "Moving the division inside the sums:\n",
      "\n",
      "$$\n",
      "\\sum_k{\\frac{1}{M} y_k} - \\sum_r{\\frac{1}{N} y_r}\n",
      "$$\n",
      "\n",
      "We can join these two sums by iterating over the entire time series $y_0, y_1 ... y_{T-1}$, selecting the weight to apply according to whether the index $t \\in 0, 1, ... T-1$ is a rest scan index or a task index:\n",
      "\n",
      "$$\n",
      "\\sum_t{ S y_t }\n",
      "\\begin{cases}\n",
      "    S = -\\frac{1}{N} ,& \\text{if } t \\in R \\\\\n",
      "    S = \\frac{1}{M} ,& \\text{if } t \\in K \\\\\n",
      "\\end{cases}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "* Make a new 4D array for `S` in the equation above, the same shape as the image data array (hint: use `np.zeros`);\n",
      "* Set the values for the array `S` to $-1/N$ for all the values in the rest volumes and $1/M$ for all the values in the task volumes (hint: use your boolean slicing arrays);\n",
      "* Do the $Sy_t$ calculation for every voxel and time point (hint: elementwise multiplication);\n",
      "* Do the sum (hint: `axis` again);\n",
      "* Check you get the same answer as you did when you subtracted the means before."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For extra points, can you change the `S` array so that `S` removes the bad volume from the difference calculation?   Do you see other ways of doing the same calculation, for example using the numpy .dot method ?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}