{
 "metadata": {
  "name": "",
  "signature": "sha256:3c8166417f7127847f677a068abc2d42a2488cdfb22325dd15b31bfccb030e00"
 },
 "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": [
      "## Range"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`range` is a built-in function in Python.\n",
      "\n",
      "\"built-in\" means you can use `range` without having to `import` it.\n",
      "\n",
      "In Python 2, `range` returns a list:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "range(10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "type(range(10))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In Python 3, range returns a special object that behaves rather like a list, but doesn't create the values until you need them.   This is more efficient in memory.\n",
      "\n",
      "Here is `range` in Python 3:\n",
      "\n",
      "```\n",
      "In [1]: range(10)\n",
      "Out[1]: range(0, 10)\n",
      "\n",
      "In [2]: type(range(10))\n",
      "Out[2]: range\n",
      "\n",
      "In [3]: list(range(10))\n",
      "Out[3]: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
      "```"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There are three ways of calling range - with 1, 2 or 3 arguments.\n",
      "\n",
      "With just one argument, the argument value is the length of the returned list, or, equivalently, the largest value in the list, plus 1 (as when you do slicing), or, equivalently, the value `stop` such that no list element is `>= stop`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "range(5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "range(2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you pass 2 arguments, the first value is the first value in the returned list, and the second value is the stop value:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "range(4, 10)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If you pass 3 arguments, the first value is the first returned value, the second is the stop value, and the third is the step size."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "range(4, 10, 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "range(4, 11, 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## np.arange"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`arange` in numpy is very like range with two important differences:\n",
      "\n",
      "* `arange` returns an array rather than a list (in Python 2 or Python 3);\n",
      "* `arange` arguments can be floating point values."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.arange(4, 11, 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.arange(4, 11, 0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Because `arange` returns arrays, you can also multiply by the step size and maybe add the start value to create equally spaced vectors:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.arange(10) * 0.5 + 4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Exponential floating point format"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have so far seen floating point numbers in this kind of format in text files : ``314.15926``.\n",
      "\n",
      "It is common for numpy and MATLAB to write floating point numbers in exponential format, like this ``3.1415926e02``.  You can read this as meaning `3.1415926 * 10**2` == `3.1415926 * 100` == `314.15926`.  Both `np.loadtxt` and the `float` function will happily convert these to floating point values in memory:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "float('2.2e02')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Legends on plots"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have already seen that we can plot two things onto the same plot:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(0, 2*np.pi, 0.1)\n",
      "y1 = np.sin(x)\n",
      "y2 = np.cos(x)\n",
      "plt.plot(x, y1)\n",
      "plt.plot(x, y2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It can be hard to remember which plot is which (is green the sin or the cos)?  `label` and `legend` are useful for that:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(x, y1, label='sin')\n",
      "plt.plot(x, y2, label='cos')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import gamma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Scipy"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "*Scipy* is a large library of scientific routines that builds on top of numpy. \n",
      "\n",
      "You can think of numpy as being a subset of MATLAB, and numpy + scipy as being as being roughly equivalent to MATLAB plus the MATLAB toolboxes."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Scipy has many sub-packages, for doing things like reading MATLAB `.mat` files (`scipy.io`) or working with sparse matrices (`scipy.sparse`). We are going to be using the functions and objects for working with statistical distributions in `scipy.stats`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.stats"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`scipy.stats` contains objects for working with many different distributions.  We are going to be working with `scipy.stats.gamma`, which implements the [gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution).\n",
      "\n",
      "In particular we are interested in the [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) (PDF) of the gamma distribution.\n",
      "\n",
      "Because this is a function, we need to pass it an array of values at which it will evaluate.\n",
      "\n",
      "We can also pass various parameters which change the shape, location and width of the gamma PDF.  The most important is the first parameter (after the input array) known as the [shape parameter](https://en.wikipedia.org/wiki/Shape_parameter) ($k$ in the [wikipedia page on gamma distributions](https://en.wikipedia.org/wiki/Gamma_distribution))."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we chose some x values at which to sample from the gamma PDF:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(0, 25, 0.1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next we plot the gamma PDF for shape values of 2, 4, 8, 12."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(x, gamma.pdf(x, 2), label='k=2')\n",
      "plt.plot(x, gamma.pdf(x, 4), label='k=4')\n",
      "plt.plot(x, gamma.pdf(x, 8), label='k=8')\n",
      "plt.plot(x, gamma.pdf(x, 12), label='k=12')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will see more from the statistics distributions next week."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Correlation coefficient"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By \"correlation coefficient\", we mean the [Pearson product moment correlation coefficient](https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient).\n",
      "\n",
      "Let's say I have two arrays, `x` and `y`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([ 20.57,  18.33,  20.8 ,  18.77,  18.46,  21.09,  19.96,  20.61,\n",
      "               18.73,  19.6 ,  18.15,  19.7 ,  20.36,  20.39,  19.67,  20.73])\n",
      "y = np.array([ 28.98,  30.13,  29.64,  29.87,  28.78,  33.48,  31.36,  31.04,\n",
      "               30.43,  31.69,  27.08,  29.25,  29.5 ,  30.04,  29.74,  30.06])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(x, y, 'o')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The numpy routine `np.corrcoef` gives me the correlation *matrix* between `x` and `y`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "corr_matrix = np.corrcoef(x, y)\n",
      "corr_matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`corr_matrix[0, 0]` is the correlation of `x` with itself, `corr_matrix[1, 0]` (and `corr_matrix[0, 1]`) is the correlation of `x` with `y`."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See the [correlation exercises](more_correlation.ipynb) for more on correlation."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Be careful - `np.correlate` is something else entirely."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## dot products"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If I have two vectors $\\mathbf{a}$ with elements $a_0, a_1, ... a_{n-1}$, and $\\mathbf{b}$ with elements $b_0, b_1, ... b_{n-1}$ then the [dot product](https://en.wikipedia.org/wiki/Dot_product) is defined as:\n",
      "\n",
      "$$\n",
      "\\mathbf{a}\\cdot \\mathbf{b} = \\sum_{i=0}^{n-1} a_ib_i = a_0b_0 + a_1b_1 + \\cdots + a_{n-1}b_{n-1}\n",
      "$$\n",
      "\n",
      "In code:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = np.arange(5)\n",
      "b = np.arange(10, 15)\n",
      "np.dot(a, b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The same thing as\n",
      "np.sum(a * b)  # Elementwise multiplication"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Matrix multiplication operates by taking dot products of the rows of the first array (matrix) with the columns of the second.\n",
      "\n",
      "Let's say I have a matrix $\\mathbf{X}$, and $X_{i,:}$ is row $i$ in $\\mathbf{X}$.  I have a matrix $\\mathbf{Y}$, and $Y_{:,j}$ is column $j$ in $\\mathbf{Y}$.  The output matrix $\\mathbf{Z} = \\mathbf{X} \\mathbf{Y}$ has entry $Z_{i,j} = X_{i,:} \\cdot Y_{:, j}$.   We will see this often over the next few weeks."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.array([[0, 1, 2], [3, 4, 5]])\n",
      "X"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y = np.array([[7, 8], [9, 10], [11, 12]])\n",
      "Y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X.dot(Y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[0, :].dot(Y[:, 0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[1, :].dot(Y[:, 0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Convolution and the tail values"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's say I have a neural vector with a couple of spikes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "times = np.arange(0, 60, 0.5)  # samples every 0.5 seconds\n",
      "neural_vector = np.zeros(times.shape)\n",
      "neural_vector[10] = 1  # At 5 seconds\n",
      "neural_vector[20] = 1  # At 10 seconds\n",
      "plt.plot(times, neural_vector)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Then I have my HRF function, sampled every half second, to match:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def hrf(t):\n",
      "    \"A hemodynamic response function\"\n",
      "    values = t ** 8.6 * np.exp(-t / 0.547)\n",
      "    # Scale max to 1\n",
      "    return values / np.max(values)\n",
      "\n",
      "hrf_times = np.arange(0, 20, 0.5)\n",
      "hrf_samples = hrf(hrf_times)\n",
      "plt.plot(hrf_times, hrf_samples)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The input neural vector is length 120, and the HRF vector is length 40:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(len(neural_vector))\n",
      "print(len(hrf_samples))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For reasons that may be familiar to you now, when we convolve the neural vector with the hrf signal, we get an output that is length 120 + 40 - 1 = 159:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hemodynamic_prediction = np.convolve(neural_vector, hrf_samples)\n",
      "len(hemodynamic_prediction)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is because of the HRF falling off the end of the input vector.  The value at index 120 in the new vector refers to time 60, and value 121 refers to time 60.5 seconds.  To retain only the values in the new hemodynamic vector that refer to times up to (not including) 60s, we can just drop the last `len(hrf_signal) - 1` values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hemodynamic_for_60s = hemodynamic_prediction[:len(neural_vector)]\n",
      "plt.plot(times, neural_vector, label='neural vector')\n",
      "plt.plot(times, hemodynamic_for_60s, label='hemodynamic prediction')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "##  Getting the time series from a voxel "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's say I have a 4D image:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nibabel as nib\n",
      "img = nib.load('ds114_sub009_t2r1.nii', mmap=False)\n",
      "data = img.get_data()\n",
      "data.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Each element slicing over the last dimension, is a 3D volume:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vol0 = data[..., 0]\n",
      "vol1 = data[..., 1]\n",
      "vol1.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The term *voxel* refers to a 3 dimensional pixel - a position in three dimensional space.\n",
      "\n",
      "These 3D volumes contain voxel values, for all of the 64, 64, 30 voxels, for a single time point.\n",
      "\n",
      "For example, let's say I was interested at the voxel at (i, j, k) coordinates (10, 23, 9):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vox_value = vol0[10, 23, 9]\n",
      "vox_value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This value is the value for this voxel (voxel at 10, 23, 9) for this time point (first time point).  This would be the value for the same voxel at the second time point (volume index 1):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vol1[10, 23, 9]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Of course I could get the value for this voxel, at every time point, by doing a loop:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vox_values = []\n",
      "for i in range(data.shape[-1]):\n",
      "    vox_values.append(data[10, 23, 9, i])\n",
      "len(vox_values)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vox_values[:15]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I can get the same effect by slicing over the last dimension:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sliced_vox_values = data[10, 23, 9, :]\n",
      "sliced_vox_values[:15]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.all(vox_values == sliced_vox_values)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}