SpeechModelTutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling the Representation of Speech in the Brain\n",
    "Originally given as a tutorial at EACL 2014 by Alex Huth.\n",
    "\n",
    "In this tutorial you will step through a voxel-wise modeling analysis. You will use computational models to extract semantic features from a natural speech stimulus. Then these features will be used to build linear models of fMRI data, and model weights and prediction performance will be visualized.\n",
    "\n",
    "If you so desire, you can step through this entire tutorial without modifying any code. But there are a few points where you will be able to make simply modifications and then see what effect those modifications have on the results. Additionally at the end you can re-run the model using phoneme features instead of semantic features.\n",
    "\n",
    "#### Acknowledgements\n",
    "This fMRI data used in this tutorial was collected by Alex Huth and Wendy de Heer at the University of California, Berkeley. All work was supervised by professors Jack Gallant and Frederic Theunissen of the UC Berkeley Psychology Department. Please do not redistribute the code or data used here. Visualization is done using [pycortex](https://pycortex.github.io/).\n",
    "\n",
    "#### Citation\n",
    "The analysis demonstrated in this tutorial forms the basis of this paper:\n",
    "[Huth, A. G. et al., \"Natural speech reveals the semantic maps that tile human cerebral cortex\" (2016) _Nature_.](https://www.nature.com/articles/nature17637)\n",
    "\n",
    "## The experiment\n",
    "In this experiment a subject underwent fMRI scanning while they listened to roughly 2 hours of natural narrative speech stimuli. These stimuli were 10-15 minute complete stories drawn from *The Moth Radio Hour*, a radio program where storytellers tell true, autobiographical stories in front of a live audience.\n",
    "\n",
    "## This notebook\n",
    "This notebook is written in python, and is presented using the iPython notebook. **To evaluate a block of code, click in it to select it, and then press shift-Enter.** You can go through the entire tutorial (until the optional section at the end) by just running each block in turn. But you'll learn more if you stop and mess around with each block as you go.\n",
    "\n",
    "* If you're totally new to python, you might check out this collection of [simple python programs](https://wiki.python.org/moin/SimplePrograms), which is a very quick intro to python's syntax and data structures. You can also find dozens of python beginner's guides on the internet.\n",
    "* If you're new to python but know MATLAB, check out [this reference](http://wiki.scipy.org/NumPy_for_Matlab_Users). It has python/numpy equivalents for most common MATLAB functions.\n",
    "* Plotting is accomplished here using matplotlib, a MATLAB-like plotting library for python. If you want to learn more about the basics of using matplotlib, check out [this notebook](http://nbviewer.ipython.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-4-Matplotlib.ipynb).\n",
    "* Finally, there are a number of points where the notebook imports libraries that I (Alex) have written. If you want to see the code for this, create a new cell (use the menu item `Insert > Insert cell below`), enter `%load libraryname.py`, and evaluate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell if your computer has a 'retina' or high DPI display. It will make the figures look much nicer.\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell imports libraries that you will need\n",
    "# Run this.\n",
    "from matplotlib.pyplot import figure, cm\n",
    "import numpy as np\n",
    "import logging\n",
    "logging.basicConfig(level=logging.DEBUG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The semantic model: English1000\n",
    "Here you will load a precomputed vector-space semantic model. This semantic model will transform any word (well, any word it knows about) into a 985-dimensional vector. This 985-dimensional semantic space has the nice property that words that are close together tend to have similar meanings. Although it would have been fun to have tried reconstructing this semantic model in this tutorial, it takes a really long time and it doesn't seem like the parameters matter that much. So today you're just going to work with the preconstructed semantic model.\n",
    "\n",
    "The semantic model was constructed using a decently large corpus of text (a couple billion words, comprising the stories used as stimuli here, 604 popular books, 2,405,569 wikipedia pages, and 36,333,459 user comments scraped from reddit.com) and a lexicon of roughly 10,000 words. We selected 985 \"basis words\" from the Wikipedia \"List of 1000 basic words\" (contrary to the title, this lost does not actually contain 1000 words, but this is where the title of the model comes from). These are common words that span many topics.\n",
    "\n",
    "We constructed a word co-occurrence matrix, $M$, with 985 rows and 10,470 columns. Iterating through the training corpus, we added 1 to $M_{ij}$ each time word $j$ appeared within 15 words of basis word $i$. The window size of 15 was selected to be large enough to suppress syntactic effects (word order) but no larger. Once the co-occurrence matrix was complete, we log transformed the counts, replacing $M_{ij}$ with $\\log(1 + M_{ij})$. Then each row of M was z-scored to correct for differences in basis word frequency, and finally each column of $M$ was z-scored to correct for word frequency. The resulting matrix is the one you're loading here.\n",
    "\n",
    "(As an aside, while I constructed this model in a totally ad hoc and unplanned way, it has properties that are very similar to Mikolov's [word2vec model](https://code.google.com/p/word2vec/) that's recently gained a lot of popularity.)\n",
    "\n",
    "Anyway, here you are going to load the model and then play with it a bit to see how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load semantic model\n",
    "# The SemanticModel class is something I wrote to make it easy to deal with vector-space semantic models.\n",
    "from SemanticModel import SemanticModel\n",
    "eng1000 = SemanticModel.load(\"data/english1000sm.hf5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can get the vector for a word by indexing the model with that word\n",
    "# For example, the vector for \"finger\":\n",
    "print(eng1000[\"finger\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing a word\n",
    "First let's plot the length 985 vector for one word to see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_word = \"finger\"\n",
    "\n",
    "f = figure(figsize=(15,5))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "ax.plot(eng1000[plot_word], 'k')\n",
    "ax.axis(\"tight\")\n",
    "ax.set_title(\"English1000 representation for %s\" % plot_word)\n",
    "ax.set_xlabel(\"Feature number\")\n",
    "ax.set_ylabel(\"Feature value\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing more than one word\n",
    "Next let's plot the vectors for three words: \"finger\", \"fingers\", and \"grief\". Here you will see that \"finger\" (in black) and \"fingers\" (in red) look very similar, but \"grief\" (in blue) looks very different. Neat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_words = [\"finger\", \"fingers\", \"language\"]\n",
    "colors = [\"k\", \"r\", \"b\"]\n",
    "\n",
    "f = figure(figsize=(15,5))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "wordlines = []\n",
    "\n",
    "for word, color in zip(plot_words, colors):\n",
    "    wordlines.append(ax.plot(eng1000[word], color)[0])\n",
    "\n",
    "ax.axis(\"tight\")\n",
    "ax.set_title(\"English1000 representations for some words\")\n",
    "ax.set_xlabel(\"Feature number\")\n",
    "ax.set_ylabel(\"Feature value\")\n",
    "ax.legend(wordlines, plot_words)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Semantic smoothness\n",
    "One nice test of a vector-space semantic model is whether it results in a \"semantically smooth\" representation of the words. That is, do nearby words in the space have intuitively similar meanings? Here you can test that using the method `find_words_like_word`. \n",
    "\n",
    "Give any word (that the model knows about), and it will print out the 10 closest words (that it knows about) and their cosine similarities (or correlations, same thing in this case). This includes the word you supplied.\n",
    "\n",
    "In this next example it prints the closest words to \"finger\". All of the 10 closest words are semantically related: 9 are nouns, and 1 is a verb (\"stick\"; of course this is also a noun, I'm just assuming that the sense of \"stick\" that's close to \"finger\" is probably the verb sense, but this brings up an important point: this model does nothing to disambiguate between different word senses!).\n",
    "\n",
    "You can put different words in here and see what the model comes up with. \n",
    "\n",
    "*(Be warned: the model knows some dirty words. It was trained using the internet, after all.)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test semantic model\n",
    "eng1000.find_words_like_word(\"finger\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is just another example, but this one an abstract noun, \"language\". Again the model does a pretty good job at finding related words."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eng1000.find_words_like_word(\"language\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eng1000.find_words_like_vec(eng1000[\"king\"] - eng1000[\"man\"] + eng1000[\"woman\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The stimuli: Moth stories\n",
    "Next we're going to load up the stimuli. We're not going to be dealing with the actual audio of the stories that were presented, but instead with aligned transcripts. These were generated using the UPenn forced aligner (P2FA), which figures out when each word was spoken given the transcript and the audio. The transcripts are stored in TextGrid format (native to Praat), which can be loaded directly into Python using some code from the natural language toolkit (NLTK).\n",
    "\n",
    "Here you will load the TextGrids for the stories, as well as 'TRfiles', which specify the time points relative to story onset when the fMRI data was collected (roughly every 2 seconds).\n",
    "\n",
    "Finally the TextGrids and TRfiles will be combined together into a representation I call a DataSequence. There is nothing interesting going on here scientifically, this is just something to make subsequent steps more manageable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These are lists of the stories\n",
    "# Rstories are the names of the training (or Regression) stories, which we will use to fit our models\n",
    "Rstories = ['alternateithicatom', 'avatar', 'howtodraw', 'legacy', \n",
    "            'life', 'myfirstdaywiththeyankees', 'naked', \n",
    "            'odetostepfather', 'souls', 'undertheinfluence']\n",
    "\n",
    "# Pstories are the test (or Prediction) stories (well, story), which we will use to test our models\n",
    "Pstories = ['wheretheressmoke']\n",
    "\n",
    "allstories = Rstories + Pstories\n",
    "\n",
    "# Load TextGrids\n",
    "from stimulus_utils import load_grids_for_stories\n",
    "grids = load_grids_for_stories(allstories)\n",
    "\n",
    "# Load TRfiles\n",
    "from stimulus_utils import load_generic_trfiles\n",
    "trfiles = load_generic_trfiles(allstories)\n",
    "\n",
    "# Make word and phoneme datasequences\n",
    "from dsutils import make_word_ds, make_phoneme_ds\n",
    "wordseqs = make_word_ds(grids, trfiles) # dictionary of {storyname : word DataSequence}\n",
    "phonseqs = make_phoneme_ds(grids, trfiles) # dictionary of {storyname : phoneme DataSequence}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before going on, let's play with the DataSequences a bit, both so you can see what the data structure looks like, and also so you can see what the stimuli look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "naked = wordseqs[\"naked\"]\n",
    "# The DataSequence stores a lot of information\n",
    "# naked.data is a list of all the words in the story\n",
    "print (\"There are %d words in the story called 'naked'\" % len(list(naked.data)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can print out the first 100 words like this\n",
    "print (list(naked.data)[:100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# or, if you want it to be more readable, like this\n",
    "print (\" \".join(list(naked.data)[:100]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the datasequence also stores when exactly each word was spoken (this time corresponds to the middle of each word)\n",
    "print (naked.data_times[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and it also stores the time of the middle of each fMRI acquisition (each acqusition takes 2.0045 seconds)\n",
    "# these times are relative to story start, so the fMRI scan started 10 seconds before the story\n",
    "print (naked.tr_times[:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# and it also makes it easy to, for example, find the words that were spoken during each fMRI acquisition\n",
    "# (the first few are empty because they came before the story started)\n",
    "print (naked.chunks()[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projecting the stimuli into the semantic space\n",
    "The next step in this analysis is that you need to project each word in the stimulus into the English1000 semantic feature space that you loaded above. I wrote a nice function to do this called `make_semantic_model` that simply takes the word DataSequence and the semantic model, and spits out a new DataSequence where each word is replaced by a 985-dimensional vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Project stimuli\n",
    "from dsutils import make_semantic_model\n",
    "semanticseqs = dict() # dictionary to hold projected stimuli {story name : projected DataSequence}\n",
    "for story in allstories:\n",
    "    semanticseqs[story] = make_semantic_model(wordseqs[story], eng1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a look at the projected stimuli\n",
    "naked_proj = semanticseqs[\"naked\"]\n",
    "\n",
    "print (naked_proj.data.shape) # prints the shape of 'data' as (rows, columns)\n",
    "print (naked_proj.data[:10]) # print the first 10 rows (this will be truncated)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Downsample the projected stimuli\n",
    "In order to build a model, you need to downsample the semantic representations of the stimuli to the same temporal scale as the fMRI responses that you will be modeling. The DataSequence provides a method that does this, called `chunksums`.\n",
    "\n",
    "For those of you who are interested, downsampling is accomplished here using a 3-lobe Lanczos filter (see [here](http://en.wikipedia.org/wiki/Lanczos_window) for details about the math). You can try changing the number of lobes, it shouldn't affect the results much."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downsample stimuli\n",
    "interptype = \"lanczos\" # filter type\n",
    "window = 3 # number of lobes in Lanczos filter\n",
    "\n",
    "downsampled_semanticseqs = dict() # dictionary to hold downsampled stimuli\n",
    "for story in allstories:\n",
    "    downsampled_semanticseqs[story] = semanticseqs[story].chunksums(interptype, window=window)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the downsampling\n",
    "Next you're going to visualize what the downsampling did. Here you're going to plot the value of one semantic feature (feature 2, which is actually the third feature: zero-based indexing) for each word, and also the downsampled vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the result\n",
    "s_words = wordseqs[\"naked\"]\n",
    "s_sem = semanticseqs[\"naked\"]\n",
    "s_semdown = downsampled_semanticseqs[\"naked\"]\n",
    "\n",
    "f = figure(figsize=(15,5))\n",
    "f.clf()\n",
    "schan = 2\n",
    "ax = f.add_subplot(1,1,1)\n",
    "wordstems = ax.stem(s_sem.data_times, \n",
    "                    s_sem.data[:,schan] / np.abs(s_sem.data[:,schan]).max(), \n",
    "                    linefmt=\"k-\", markerfmt=\"k.\", basefmt=\"k-\")\n",
    "interps = ax.plot(s_sem.tr_times, \n",
    "                  s_semdown[:,schan] / np.abs(s_semdown[:,schan]).max(), 'r.-')\n",
    "ax.set_xlim(-6, 60)\n",
    "ax.set_ylim(-1, 1)\n",
    "ax.set_xlabel(\"Time (seconds since story start)\")\n",
    "ax.set_ylabel(\"Semantic feature value\")\n",
    "ax.legend((wordstems, interps[0]), (\"Individual words\", \"Downsampled feature\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Concatenating features across stories\n",
    "Next you're going to combine together all the features from all the stories into one big matrix. Within this operation, you're also going to [z-score](http://en.wikipedia.org/wiki/Z-score) each feature within each story. This operation subtracts off the mean and then divides by the standard deviation. This might seem like a weird or incomprehensible thing to do, but I do it because the responses to each story are z-scored individually. Anyway not a big deal.\n",
    "\n",
    "The features for each story are also trimmed a bit (the variable `trim` determines how many time points are removed from the beginning and end of each story). The fMRI responses at the beginnings and ends of the stories are often noisier than at other times because of transients and problems with detrending (an fMRI preprocessing step that you don't need to worry about here aside from this point).\n",
    "\n",
    "The combined features are stored in big matrices called `Rstim` (with the training, or Regression stimuli) and `Pstim` (with the test, or Prediction stimuli)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine stimuli\n",
    "from npp import zscore\n",
    "trim = 5\n",
    "Rstim = np.vstack([zscore(downsampled_semanticseqs[story][5+trim:-trim]) for story in Rstories])\n",
    "Pstim = np.vstack([zscore(downsampled_semanticseqs[story][5+trim:-trim]) for story in Pstories])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "storylens = [len(downsampled_semanticseqs[story][5+trim:-trim]) for story in Rstories]\n",
    "print(storylens)\n",
    "\n",
    "print(np.cumsum(storylens))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the sizes of these matrices\n",
    "print (\"Rstim shape: \", Rstim.shape)\n",
    "print (\"Pstim shape: \", Pstim.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the combined stimuli\n",
    "Next you're going to plot some of the feature channels. This is just to see what the feature look like that are going to go into the regression model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot some of the combined stimuli\n",
    "f = figure(figsize=(20, 5))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "\n",
    "for ii in range(10):\n",
    "    # Plot each feature, offset by 5 vertically so they are easier to see\n",
    "    ax.plot(Rstim[:750,ii] - 5 * ii)\n",
    "\n",
    "ax.set_xlim(0, 750)\n",
    "ax.set_yticks([])\n",
    "ax.set_xticks(range(0, 750, 50))\n",
    "ax.set_xlabel(\"Time (fMRI volumes)\")\n",
    "ax.set_ylabel(\"Features 1-10\")\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Concatenate delayed stimuli for FIR model\n",
    "Next you are going to concatenate multiple delayed versions of the stimuli, in order to create a linear [finite impulse response (FIR) model](http://en.wikipedia.org/wiki/Fir_filter). This is a vitally important step, and is conceptually a bit difficult, so take a few minutes to make sure you understand what is going on here.\n",
    "\n",
    "#### Background: the hemodynamic response\n",
    "First you need to understand the problem that the FIR model is solving. fMRI measures the blood-oxygen level dependent (BOLD) signal, which is a complicated and nonlinear combination of blood oxygenation and blood volume. When neurons in an area of the brain become active, they start using up lots of energy. To compensate, nearby blood vessels dilate so that more oxygen and glucose become available to the neurons. The resulting changes in blood oxygenation (which increases) and volume (which also increases) create the magnetic signature that is recorded by fMRI. \n",
    "\n",
    "But this process is **slow**. It takes seconds after the neural activity begins for the blood vessels to dilate and for the BOLD response to become apparent. And then it takes more seconds for the response to go away. So although a neural response might only last milliseconds, the associated BOLD response will rise and fall over a span of maybe 10 seconds, orders of magnitude slower. The shape of this rise and fall is called the [hemodynamic response function (HRF)](http://en.wikipedia.org/wiki/Haemodynamic_response).\n",
    "\n",
    "Here is a pretty standard looking example of an HRF:\n",
    "\n",
    "<img src='http://www.brainmatters.nl/wp-content/uploads/bold.png' width=350px></img>\n",
    "\n",
    "#### FIR model\n",
    "To accurately model how the brain responds to these stimuli we must also model the HRF. There are many ways to do this. The most common is to assume that the HRF follows a canonical shape. But this approach turns out to not work very well: different parts of the brain have very different vasculature (blood vessels), so the HRF shape can vary a lot. \n",
    "\n",
    "Instead, what you are going to here is estimate a separate HRF for each semantic feature in each voxel that is being modeled. This estimate is going to take the form of a linear finite impulse response (FIR) model. The linear FIR form is particularly nice to use because it's very simple to estimate and powerful (if anything, it might be too powerful.. more on that later). To build a linear FIR model all you have to do is concatenate together multiple delayed copies of the stimulus. I usually use four delays: 1, 2, 3, and 4 time points. The resulting delayed features can be thought of as representing the stimulus 1, 2, 3, and 4 time points ago. So the regression weights for those features will represent how a particular voxel responds to a feature 1, 2, 3, or 4 time points in the past, and these regression weights are a 4-point estimate of the HRF for that feature in that voxel.\n",
    "\n",
    "The potential downside of the FIR model is that it may be too expressive. Each feature in each voxel is allowed to have any HRF, but this comes at the cost of multiplying the total number of regression weights that we must fit by the number of delays. In all likelihood the true HRFs vary, but they don't vary that much, so we probably don't need this many independent features. This cost becomes apparent if you increase the number of delays. This will slow down model fitting and likely decrease the stability of the regression weights, leading to decreased model performance. \n",
    "\n",
    "Feel free to play around with the number of delays and see how it affects the model results!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Delay stimuli\n",
    "from util import make_delayed\n",
    "ndelays = 4\n",
    "delays = range(1, ndelays+1)\n",
    "\n",
    "print (\"FIR model delays: \", delays)\n",
    "\n",
    "delRstim = make_delayed(Rstim, delays)\n",
    "delPstim = make_delayed(Pstim, delays)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the sizes of these matrices\n",
    "print (\"delRstim shape: \", delRstim.shape)\n",
    "print (\"delPstim shape: \", delPstim.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing FIR features\n",
    "Here you will visualize the first semantic feature at each of the delays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the same feature at different delays\n",
    "f = figure(figsize=(15, 4))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "for ii in range(ndelays):\n",
    "    ax.plot(delRstim[:500, ii * Rstim.shape[1]] - 5 * ii)\n",
    "ax.set_xlim(0, 500)\n",
    "ax.set_yticks([])\n",
    "ax.set_xticks(range(0, 500, 50))\n",
    "ax.set_xlabel(\"Time (fMRI volumes)\")\n",
    "ax.set_ylabel(\"Feature 1 across delays\")\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Response data\n",
    "Next you will load the fMRI data. This is totally the most exciting part! These responses have already been preprocessed (the 3D images were motion corrected and aligned to each other, detrended, and then z-scored within each stimulus) so you don't have to worry about that.\n",
    "\n",
    "You will load three different variables: `zRresp`, the responses to the regression dataset; `zPresp`, the responses to the prediction dataset; and `mask`, which is a 3D mask showing which voxels have been selected (we are not modeling every voxel in the scan, that would take forever, we are only modeling the voxels that overlap with the cerebral cortex)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load responses\n",
    "import tables\n",
    "resptf = tables.open_file(\"data/fmri-responses.hf5\")\n",
    "zRresp = resptf.root.zRresp.read()\n",
    "zPresp = resptf.root.zPresp.read()\n",
    "mask = resptf.root.mask.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print matrix shapes\n",
    "print (\"zRresp shape (num time points, num voxels): \", zRresp.shape)\n",
    "print (\"zPresp shape (num time points, num voxels): \", zPresp.shape)\n",
    "print (\"mask shape (Z, Y, X): \", mask.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize where the voxels are coming from (mask)\n",
    "Next you will visualize where the voxels are coming from in the brain. This will give you an idea of where the data come from.\n",
    "\n",
    "First you will plot a single slice through the mask in the Z dimension. This is called an \"axial\" slice. The top of the image is the front of the brain, the bottom is the back. The left side of the image is the right side of the brain, and the right side of the image is the left side of the brain (as if you are looking up at the brain from under the subject's chin; this left-right reversal is often referred to as \"radiological coordinates\", as opposed to \"neurological coordinates\" where you are looking down from the top).\n",
    "\n",
    "Then you will plot a mosaic of all the slices. This is done using the function `mosaic` from James Gao's pyCortex package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot one slice of the mask that was used to select the voxels\n",
    "f = figure()\n",
    "ax = f.add_subplot(1,1,1)\n",
    "ax.matshow(mask[16], interpolation=\"nearest\", cmap=cm.gray) # show the 17th slice of the mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot mask mosaic\n",
    "import cortex\n",
    "f = figure(figsize=(10,10))\n",
    "cortex.mosaic(mask, cmap=cm.gray, interpolation=\"nearest\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the responses of a few voxels over time\n",
    "Next you will visualize the responses of a few selected voxels over time. I selected these particular voxels because they are reasonably well explained by the semantic model, but have some differences in their responses across time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the response of a few voxels over time\n",
    "selvoxels = [20710, 27627, 24344, 34808, 22423, 25397]\n",
    "\n",
    "f = figure(figsize=(15, 5))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "for ii,vi in enumerate(selvoxels):\n",
    "    ax.plot(zRresp[:500, vi] - 5 * ii)\n",
    "ax.set_xlim(0, 500)\n",
    "ax.set_yticks([])\n",
    "ax.set_xticks(range(0, 500, 50))\n",
    "ax.set_xlabel(\"Time (fMRI volumes)\")\n",
    "ax.set_ylabel(\"Voxel responses\")\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression model\n",
    "Finally, the core of the analysis: you will fit a regression model that predicts the responses of each voxel as a weighted sum of the semantic features. This model will then be tested using a held out dataset (the Prediction dataset). And if the model proves to be reasonably predictive, then the weights of the regression model will tell us what semantic features each voxel responds to.\n",
    "\n",
    "This is a linear regression model, so if the response time course for voxel $j$ is $R_j$, the stimulus time course for semantic feature $i$ is $S_i$, and the regression weight for feature $i$ in voxel $j$ is $\\beta_{ij}$, then the model can be written as:\n",
    "\n",
    "$$\\hat{R}_j = \\beta_{0j} S_0 + \\beta_{1j} S_1 + \\cdots$$\n",
    "\n",
    "or:\n",
    "\n",
    "$$\\hat{R}_j = \\sum_i \\beta_{ij} S_i$$\n",
    "\n",
    "The trick, of course, is accurately estimating the $\\beta_j$ values. This is commonly done by minimizing the sum of the squared error (here across time, $t$):\n",
    "\n",
    "$$E_j(\\beta) = \\sum_t (R_{jt} - \\hat{R}_{jt})^2 = \\sum_t (R_{jt} - \\sum_i \\beta_{i} S_{it})^2$$\n",
    "\n",
    "$$\\beta_j = \\underset{\\beta}{\\operatorname{argmin}} E_j(\\beta)$$\n",
    "\n",
    "Computing $\\beta$ this way is called ordinary least squares (OLS), and this will not work in our case because the total number of features (3940) is smaller than the number of time points (3737). (It would be possible if the number of delays was smaller than 4, but it would give terrible results.. feel free to try it! OLS can be performed using the function `np.linalg.lstsq`.)\n",
    "\n",
    "In almost every case, linear regression can be improved by making some prior assumptions about the weights (or, equivalently, about the covariance structure of the stimuli). This is called **regularization**, or **regularized linear regression**. One way to do this is to penalize the error function by the sum of the squared weights. This is commonly known as **ridge regression**, and is a special case of [Tikhonov regularization](http://en.wikipedia.org/wiki/Ridge_regression). It finds the $\\beta$ that minimizes the following error function:\n",
    "\n",
    "$$E_j(\\beta) = \\sum_t (R_{jt} - \\sum_i \\beta_{i} S_{it})^2 + \\alpha \\sum_i \\beta_i^2$$\n",
    "\n",
    "(In practice we will use a different formulation that involves re-weighting the singular values of the matrix $S$ before computing its pseudoinverse. This method achieves the same results but is extremely efficient because it uses all the linear algebra machinery that computers are so good at to build many models in parallel.)\n",
    "\n",
    "### The hyperparameter: $\\alpha$\n",
    "You may have noticed in the equation above that we have introduced a new parameter, $\\alpha$, which controls the strength of the regularization. If $\\alpha$ is set to zero, then we get back to exactly the OLS formulation (above). As $\\alpha$ goes to infinity, the regularization forces all the weights to go to zero (in practice this also has the slightly weirder effect of making all the weights independent, as if each feature was regressed separately on the responses).\n",
    "\n",
    "So how do we choose $\\alpha$? We're going to do it here using cross-validation. First, we split the Regression dataset up into two parts. Then we estimate the weights for a given $\\alpha$ on the first part, and test how well we can predict responses on the second part. This is repeated for each possible $\\alpha$ that we want to test, and for a couple different splits of the Regression dataset. Then we find the $\\alpha^*$ that gave us the best predictions within the split Regression dataset. Finally we estimate the weights using the entire Regression dataset and the selected $\\alpha^*$.\n",
    "\n",
    "Because this is an annoying and laborious process, I've encapsulated it within the function `bootstrap_ridge`. You simply give this function your datasets, the possible $\\alpha$ values, and a few parameters for the cross-validation, and it does all the rest. The parameter `nboots` determines the number of cross-validation tests that will be run. \n",
    "\n",
    "To do cross-validation, `bootstrap_ridge` divides the Regression dataset into many small chunks, and then splits those chunks into the two groups that will be used to estimate weights and test $\\alpha$ values. This is better than just choosing individual time points because both the fMRI data and stimuli are autocorrelated (i.e. correlated across time). The parameter `chunklen` determines the length of the chunks, and the parameter `nchunks` determines the number of chunks in the $\\alpha$-testing dataset. By default I set `chunklen` to 40 time points (80-second chunks), and set `nchunks` to 20 (40 * 20 = 800 time points for testing $\\alpha$ values, 3737-800 = 2937 time points for estimating weights). These values should not matter too much.\n",
    "\n",
    "Running the regression will take a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run regression\n",
    "from ridge import bootstrap_ridge\n",
    "alphas = np.logspace(1, 3, 10) # Equally log-spaced alphas between 10 and 1000. The third number is the number of alphas to test.\n",
    "nboots = 1 # Number of cross-validation runs.\n",
    "chunklen = 40 # \n",
    "nchunks = 20\n",
    "\n",
    "wt, corr, alphas, bscorrs, valinds = bootstrap_ridge(delRstim, zRresp, delPstim, zPresp,\n",
    "                                                     alphas, nboots, chunklen, nchunks,\n",
    "                                                     singcutoff=1e-10, single_alpha=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = figure()\n",
    "ax = f.add_subplot(1,1,1)\n",
    "ax.semilogx( np.logspace(1, 3, 10), bscorrs.mean(2).mean(1), 'o-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variables returned by the regression\n",
    "Next let's have a look at the variables returned by the regression function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# wt is the regression weights\n",
    "print (\"wt has shape: \", wt.shape)\n",
    "\n",
    "# corr is the correlation between predicted and actual voxel responses in the Prediction dataset\n",
    "print (\"corr has shape: \", corr.shape)\n",
    "\n",
    "# alphas is the selected alpha value for each voxel, here it should be the same across voxels\n",
    "print (\"alphas has shape: \", alphas.shape)\n",
    "\n",
    "# bscorrs is the correlation between predicted and actual voxel responses for each round of cross-validation\n",
    "# within the Regression dataset\n",
    "print (\"bscorrs has shape (num alphas, num voxels, nboots): \", bscorrs.shape)\n",
    "\n",
    "# valinds is the indices of the time points in the Regression dataset that were used for each\n",
    "# round of cross-validation\n",
    "print (\"valinds has shape: \", np.array(valinds).shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing the regression models by predicting responses\n",
    "The `bootstrap_ridge` function already computed predictions and correlations for the Prediction dataset, but this is important so let's reproduce that step more explicitly.\n",
    "\n",
    "Remember that according to the linear model, the predicted responses for each voxel are a weighted sum of the semantic features. An easy way to compute that is by taking the dot product between the weights and semantic features: $$\\hat{R} = S \\beta$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predict responses in the Prediction dataset\n",
    "\n",
    "# First let's refresh ourselves on the shapes of these matrices\n",
    "print (\"zPresp has shape: \", zPresp.shape)\n",
    "print (\"wt has shape: \", wt.shape)\n",
    "print (\"delPstim has shape: \", delPstim.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Then let's predict responses by taking the dot product of the weights and stim\n",
    "pred = np.dot(delPstim, wt)\n",
    "\n",
    "print (\"pred has shape: \", pred.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualizing predicted and actual responses\n",
    "Next let's plot some predicted and actual responses side by side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = figure(figsize=(15,5))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "\n",
    "selvox = 20710 # a decent voxel\n",
    "\n",
    "realresp = ax.plot(zPresp[:,selvox], 'k')[0]\n",
    "predresp = ax.plot(pred[:,selvox], 'r')[0]\n",
    "\n",
    "ax.set_xlim(0, 291)\n",
    "ax.set_xlabel(\"Time (fMRI time points)\")\n",
    "\n",
    "ax.legend((realresp, predresp), (\"Actual response\", \"Predicted response\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualizing predicted and actual responses cont'd\n",
    "You might notice above that the predicted and actual responses look pretty different scale-wise, although the patterns of ups and downs are vaguely similar. But we don't really care about the scale -- for fMRI it's relatively arbitrary anyway, so let's rescale them both to have unit standard deviation and re-plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = figure(figsize=(15,5))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "\n",
    "selvox = 20710 # a good voxel\n",
    "\n",
    "realresp = ax.plot(zPresp[:,selvox], 'k')[0]\n",
    "predresp = ax.plot(zscore(pred[:,selvox]), 'r')[0]\n",
    "\n",
    "ax.set_xlim(0, 291)\n",
    "ax.set_xlabel(\"Time (fMRI time points)\")\n",
    "\n",
    "ax.legend((realresp, predresp), (\"Actual response\", \"Predicted response (scaled)\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you see that the actual and scaled predicted responses look very similar. We can quantify this similarity by computing the correlation between the two (correlation is scale-free, so it effectively automatically does the re-scaling that we did here). This voxel has high correlation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute correlation between single predicted and actual response\n",
    "# (np.corrcoef returns a correlation matrix; pull out the element [0,1] to get \n",
    "# correlation between the two vectors)\n",
    "voxcorr = np.corrcoef(zPresp[:,selvox], pred[:,selvox])[0,1]\n",
    "print (\"Correlation between predicted and actual responses for voxel %d: %f\" % (selvox, voxcorr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Computing correlations for all voxels\n",
    "Next let's compute this correlation for every voxel in the dataset. There are some very efficient ways to do this, but here I've written a for loop so that it's very explicit what's happening. (This should give exactly the same values as the variable `corr`, which was returned by `bootstrap_ridge`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxcorrs = np.zeros((zPresp.shape[1],)) # create zero-filled array to hold correlations\n",
    "for vi in range(zPresp.shape[1]):\n",
    "    voxcorrs[vi] = np.corrcoef(zPresp[:,vi], pred[:,vi])[0,1]\n",
    "print (voxcorrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing correlations across the brain\n",
    "Let's start with a supposition: the correlation should not be high everywhere, even if this is a good model of how the brain represents the semantic content of speech. There are parts of the brain that just don't respond to speech, so the correlation should be low in those areas. There are other parts of the brain that respond to speech, but maybe don't represent semantic information, so the correlation should be low in those areas as well.\n",
    "\n",
    "But let's begin by plotting a histogram of the correlations across the entire brain. This will show generally whether the model is working well or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histogram of correlations\n",
    "f = figure(figsize=(8,8))\n",
    "ax = f.add_subplot(1,1,1)\n",
    "ax.hist(voxcorrs, 100) # histogram correlations with 100 bins\n",
    "ax.set_xlabel(\"Correlation\")\n",
    "ax.set_ylabel(\"Num. voxels\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the semantic features didn't capture anything about brain activity, then we would expect the histogram to be symmetric and centered around zero. But here we see that it's highly skewed, with lots of positive values. This looks good! This model is working!\n",
    "\n",
    "Next, let's plot a mosaic of the correlations across the brain, as we plotted the mask earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot mosaic of correlations\n",
    "corrvolume = np.zeros(mask.shape)\n",
    "corrvolume[mask>0] = voxcorrs\n",
    "\n",
    "f = figure(figsize=(10,10))\n",
    "cortex.mosaic(corrvolume, vmin=0, vmax=0.5, cmap=cm.hot);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3D visualization of correlations\n",
    "In the mosaic we can see that there seem to be some concentrated areas of high correlation. But it's hard to say where in the brain those areas are based on the mosaic. So next you're going to create a fancy 3D visualization of the correlations using pyCortex.\n",
    "\n",
    "Once you've opened the viewer you'll be presented with a 3D view of the brain with colors showing the correlations. White outlines and labels show the locations of known brain areas (motor, somatosensory, visual, and some language areas). Drag around with your left mouse button to rotate the view, and the right mouse button to zoom in or out. \n",
    "\n",
    "By default you'll see a view of the cortex as it looks in reality: folded and convoluted. To better see parts of the brain that are hidden down in the folds, you can press \"i\" to see an inflated view (or drag the \"Mix\" slider at the bottom of the screen to the middle). This helps to see the data, but you will still need to rotate the brain to see all of it. To make the entire cortex visible at once, you can press \"f\" to see a flattened view. To create this view we cut the cortical surface at a few locations, and then flattened it out so that it can all be seen at once (but this introduces some distortions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot correlations on cortex\n",
    "import cortex\n",
    "corrvol = cortex.Volume(corr, \"S1\", \"fullhead\", mask=mask, vmin=0, vmax=0.5, cmap='hot')\n",
    "cortex.webshow(corrvol, port=8889, open_browser=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View 3D model\n",
    "# You will need to change where it says SERVERIP below to the IP you are connected to\n",
    "from IPython.display import HTML\n",
    "HTML(\"<a target='_blank' href='http://127.0.0.1/:8889'>Click here for viewer</a>\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simpler view of the correlations\n",
    "pyCortex also offers a simpler way to view the correlations. This method only shows the flat view, but can be embedded right here in the ipython notebook. This should look like the flat view in the 3D viewer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot correlation flatmap\n",
    "cortex.quickshow(corrvol, with_rois=False, with_labels=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What semantic features are the voxels responding to?\n",
    "Now that we have a working model, let's try to figure out what semantic features are making each voxel respond. One way to do this is to simulate how the voxel will respond to individual words, and then find the most preferred words for that voxel.\n",
    "\n",
    "But first we have an issue to contend with: we have separate weights for each delay. We could look at the weights for each delay, but instead here you will average the weights across delays to get a single set of weights for the voxel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Undelay voxel weights (average across delays)\n",
    "import operator\n",
    "from functools import reduce\n",
    "udwt = reduce(operator.add, np.split(wt/ndelays, ndelays))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "udwt.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next you will pick which voxel to visualize. Since many voxels are modeled poorly, we will pick from among the best modeled voxels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort voxels by correlation so that we can pick a good voxel\n",
    "# This will sort voxels in decreasing order of correlation\n",
    "corrsort = np.argsort(corr)[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define function that will print best words for a voxel\n",
    "import pprint\n",
    "\n",
    "def print_voxel_words(voxnum):\n",
    "    # find_words_like_vec returns 10 words most correlated with the given vector, and the correlations\n",
    "    voxwords = eng1000.find_words_like_vec(udwt[:,voxnum])\n",
    "    print (\"Best words for voxel %d (correlation %0.3f):\" % (voxnum, voxcorrs[voxnum]))\n",
    "    pprint.pprint(voxwords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print best words for some voxels\n",
    "print_voxel_words(corrsort[0]) # best voxel\n",
    "print_voxel_words(corrsort[14]) # 15th best voxel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### That's it!\n",
    "That's the semantic model! Since you made it this far, well done. If you're interested, you can go back and try changing some of the parameters and see how it affects the model. One easy parameter to change is the number of delays (or the delays themselves). Try using just one delay. Or try using 10 delays (that might be slow). You could also try pruning off some of the semantic features. How does the model work if you only use the first 100 semantic features?\n",
    "\n",
    "Alternatively, you can try using a different type of feature to model the fMRI responses: phonemes. Below are some blocks of code that will create stimulus vectors representing the number of times each different phoneme is spoken. The phoneme model will predict some voxels much better than the semantic model, and some voxels worse.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optional: phoneme model\n",
    "Another type of information that the brain extracts from speech is the phonemic content. The following blocks of code will extract phonemic features from the stimuli. Run this, and then go back and run the \"Regression Model\" block, above. You can compare the correlations of the semantic model and phoneme model to see which works best in each voxel.\n",
    "\n",
    "You can also visualize the phoneme stimuli and model as you go, building on the code blocks used above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create phoneme histogram DataSequences\n",
    "from dsutils import histogram_phonemes2, phonemes\n",
    "phonemehistseqs = dict() # dictionary to hold phoneme histograms {story name : DataSequence}\n",
    "for story in allstories:\n",
    "    phonemehistseqs[story] = histogram_phonemes2(phonseqs[story])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Phonemes were labeled using the ARPABET. The labeled phonemes are listed here.\n",
    "print (phonemes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downsample phoneme histograms\n",
    "interptype = \"lanczos\"\n",
    "window = 3\n",
    "\n",
    "downsampled_phonemehistseqs = dict()\n",
    "for story in allstories:\n",
    "    downsampled_phonemehistseqs[story] = phonemehistseqs[story].chunksums(interptype, window=window)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine phoneme stimuli\n",
    "trim = 5\n",
    "phRstim = np.vstack([np.nan_to_num(zscore(downsampled_phonemehistseqs[story][5+trim:-trim])) for story in Rstories])\n",
    "phPstim = np.vstack([np.nan_to_num(zscore(downsampled_phonemehistseqs[story][5+trim:-trim])) for story in Pstories])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Delay stimuli\n",
    "ndelays = 4\n",
    "delays = range(1, ndelays+1)\n",
    "\n",
    "delRstim = make_delayed(phRstim, delays)\n",
    "delPstim = make_delayed(phPstim, delays)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now go back to \"Regression Model\" and run that to fit a phoneme-based model!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}