preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# default_exp preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing\n",
    "\n",
    "> Functions for preprocessing fmri data and preparing stimulus and fmri data for training voxel-wise encoding models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "from nbdev.showdoc import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "#export\n",
    "import os\n",
    "import warnings\n",
    "import numpy as np\n",
    "import joblib\n",
    "from nilearn.masking import unmask, apply_mask\n",
    "from nibabel import save, load, Nifti1Image\n",
    "from nilearn.signal import clean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing BOLD fMRI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def preprocess_bold_fmri(bold, mask=None, detrend=True, standardize='zscore', **kwargs):\n",
    "    '''Preprocesses BOLD data and returns ndarray of preprocessed data\n",
    "\n",
    "    Parameters\n",
    "\n",
    "        bold : path to bold nifti file or loaded bold nifti\n",
    "        mask : path to mask nifti file or loaded mask nifti, optional\n",
    "        detrend : bool, whether to linearly detrend the data, optional\n",
    "        standardize : {‘zscore’, ‘psc’, False}, default is ‘zscore’\n",
    "        kwargs : further arguments for nilearn's clean function\n",
    "\n",
    "    Returns\n",
    "        ndarray of the preprocessed bold data in (samples, voxels)\n",
    "    '''\n",
    "    if mask:\n",
    "        data = apply_mask(bold, mask)\n",
    "    else:\n",
    "        if not isinstance(bold, Nifti1Image):\n",
    "            data = load(bold).get_data()\n",
    "        else:\n",
    "            data = bold.get_data()\n",
    "        data = np.reshape(data, (-1, data.shape[-1])).T\n",
    "    return clean(data, detrend=detrend, standardize=standardize, **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`preprocess_bold_fmri` preprocessed a BOLD Nifti and returns a numpy ndarray of the optionally masked and preprocessed fMRI data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "def test_preprocess_bold_fmri():\n",
    "    test_nifti = Nifti1Image(np.ones((2, 2, 2, 2)), affine=np.eye(4))\n",
    "    assert np.allclose(preprocess_bold_fmri(test_nifti, standardize=True, detrend=False), 0.)\n",
    "    assert np.allclose(preprocess_bold_fmri(test_nifti, standardize=False, detrend=True), 0.)\n",
    "    assert np.allclose(preprocess_bold_fmri(test_nifti, standardize=False, detrend=False), 1.)\n",
    "    mask_array = np.ones((2, 2, 2))\n",
    "    mask_array[:,0,:] = 0.\n",
    "    mask = Nifti1Image(mask_array, affine=np.eye(4))\n",
    "    assert preprocess_bold_fmri(test_nifti, mask=mask).shape == (2,4)\n",
    "\n",
    "test_preprocess_bold_fmri()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def get_remove_idx(lagged_stimulus, remove_nan=True):\n",
    "    '''Returns indices of rows in lagged_stimulus to remove'''\n",
    "    if remove_nan is True:\n",
    "        return np.where(np.any(np.isnan(lagged_stimulus), axis=1))[0]\n",
    "    elif remove_nan is False:\n",
    "        # This will raise an error if it is supplied to np.delete\n",
    "        # which is what we want\n",
    "        return None\n",
    "    elif remove_nan <= 1. and remove_nan >= 0.:\n",
    "        return np.where(np.isnan(lagged_stimulus).mean(axis=1) > remove_nan)[0]\n",
    "    else:\n",
    "        raise ValueError('remove_nan needs to be either True, False, or a float between 0 and 1.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "def test_get_remove_idx():\n",
    "    lagged_data = np.ones((100, 10))\n",
    "    lagged_data[:5, :] = np.nan\n",
    "    lagged_data[5:10, :3] = np.nan\n",
    "    assert get_remove_idx(lagged_data, remove_nan=False) is None\n",
    "    assert np.all(get_remove_idx(lagged_data, remove_nan=True)  == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])) \n",
    "    assert np.all(get_remove_idx(lagged_data, remove_nan=0.5)  == np.array([0, 1, 2, 3, 4])) \n",
    "    assert np.all(get_remove_idx(lagged_data, remove_nan=0.1)  == np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])) \n",
    "\n",
    "test_get_remove_idx()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test = [1,2,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# export\n",
    "\n",
    "def make_lagged_stimulus(stimulus, n_lags, fill_value=np.nan):\n",
    "    '''Generates a lagged stimulus representation by adding nans'''\n",
    "    lagged_reps = [np.vstack([np.full((lag_i, stimulus.shape[1]), fill_value), stimulus[:-lag_i]])\n",
    "                                 for lag_i in range(1, n_lags)]\n",
    "    return np.hstack([stimulus]+lagged_reps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "\n",
    "def test_make_lagged_stimulus():\n",
    "    stimulus = np.arange(80)[:,None]\n",
    "    lagged_data = make_lagged_stimulus(stimulus, 3)\n",
    "    assert lagged_data.shape == (80, 3)\n",
    "    assert not np.any(np.isnan(lagged_data[:, 0]))\n",
    "    assert np.isnan(lagged_data[:, -1]).sum() == 2\n",
    "    assert np.isnan(lagged_data[0]).sum() == 2\n",
    "    lagged_data = make_lagged_stimulus(stimulus, 10)\n",
    "    assert np.isnan(lagged_data[0]).sum() == 9\n",
    "    assert np.isnan(lagged_data[:, -1]).sum() == 9\n",
    "    assert np.isnan(lagged_data[:, -2]).sum() == 8\n",
    "    lagged_data = make_lagged_stimulus(stimulus, 3, fill_value=1)\n",
    "    assert lagged_data.shape == (80, 3)\n",
    "    assert np.isnan(lagged_data[0]).sum() == 0\n",
    "\n",
    "test_make_lagged_stimulus()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aligning stimulus and fMRI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def generate_lagged_stimulus(stimulus, fmri_samples, TR, stim_TR,\n",
    "                             lag_time=None, start_time=0., offset_stim=0.,\n",
    "                             fill_value=np.nan):\n",
    "    '''Generates a lagged stimulus representation temporally aligned with the fMRI data\n",
    "\n",
    "    Parameters\n",
    "\n",
    "        stimuli : ndarray, stimulus representation of shape (samples, features)\n",
    "        fmri_samples : int, samples of corresponding fmri run\n",
    "        TR : int, float, repetition time of the fMRI data in seconds\n",
    "        stim_TR : int, float, repetition time of the stimulus in seconds\n",
    "        lag_time : int, float, or None, optional,\n",
    "               lag to introduce for stimuli in seconds,\n",
    "               if no lagging should be done set this to TR or None\n",
    "        start_time :  int, float, optional, default 0.\n",
    "                  starting time of the stimulus relative to fMRI recordings in seconds\n",
    "                  appends fill_value to stimulus representation to match fMRI and stimulus\n",
    "        offset_stim : int, float, optional, default 0.\n",
    "                  time to offset stimulus relative to fMRI in the lagged stimulus,\n",
    "                  i.e. when predicting fmri at time t use only stimulus features\n",
    "                  before t-offset_stim. This reduces the number of time points used\n",
    "                  in the model.\n",
    "        fill_value : int, float, or any valid numpy array element, optional, default np.nan\n",
    "                 appends fill_value to stimulus array to account for starting_time\n",
    "                 use np.nan here with remove_nans=True to remove fmri/stimulus samples where no stimulus was presented\n",
    "\n",
    "    Returns:\n",
    "        ndarray of the lagged stimulus of shape (samples, lagged features)\n",
    "    '''\n",
    "    # find out temporal alignment\n",
    "    stim_samples_per_TR = TR / stim_TR\n",
    "    if stim_samples_per_TR < 1:\n",
    "        raise ValueError('Stimulus TR is larger than fMRI TR')\n",
    "    # check if result is close to an integer\n",
    "    if not np.isclose(stim_samples_per_TR, np.round(stim_samples_per_TR)):\n",
    "        warnings.warn('Stimulus timing and fMRI timing do not align. '\n",
    "        'Stimulus samples per fMRI samples: {0} for stimulus TR {1} and fMRI TR {2}. '\n",
    "        'Proceeds by rounding stimulus samples '\n",
    "        'per TR.'.format(stim_samples_per_TR, stim_TR, TR), RuntimeWarning)\n",
    "    stim_samples_per_TR = int(np.round(stim_samples_per_TR))\n",
    "    if lag_time is None:\n",
    "        lag_time = TR\n",
    "    if np.isclose(lag_time, 0.):\n",
    "        warnings.warn('lag_time set to 0 - to disable lagging set lag_time either to None or TR ({}).'\n",
    "                      'Defaulting to lag_time=None.'.format( TR))\n",
    "        lag_time = TR\n",
    "        \n",
    "    if lag_time < TR:\n",
    "        warnings.warn('lag_time ({}) should not be smaller than TR ({}).'.format(lag_time, TR))\n",
    "    # check if lag time is multiple of TR\n",
    "    if not np.isclose(lag_time / TR, np.round(lag_time / TR)):\n",
    "        raise ValueError('lag_time should be a multiple of TR so '\n",
    "                'that stimulus/fMRI alignment does not change.')\n",
    "    if lag_time == TR:\n",
    "            warnings.warn('lag_time is None or equal to TR, no stimulus lagging will be done.', RuntimeWarning)\n",
    "    lag_TR = int(np.round(lag_time / TR))\n",
    "    offset_TR = int(np.round(offset_stim / TR))\n",
    "\n",
    "    n_features = stimulus.shape[1]\n",
    "    n_append = 0\n",
    "    # check if the stimulus start time is moved w.r.t. fmri\n",
    "    n_prepend = int(np.round(start_time / stim_TR))\n",
    "    stimulus = np.vstack([np.full((n_prepend, n_features), fill_value), stimulus])\n",
    "\n",
    "    # make reshapeable by appending filler\n",
    "    if stimulus.shape[0] % stim_samples_per_TR > 0:\n",
    "        # either remove part of the stimulus (if it is longer than fmri) or append filler\n",
    "        if stimulus.shape[0] / stim_samples_per_TR > fmri_samples:\n",
    "            stimulus = stimulus[:-(stimulus.shape[0] % stim_samples_per_TR)]\n",
    "        else:\n",
    "            n_append = stim_samples_per_TR - ((stimulus.shape[0]) % stim_samples_per_TR)\n",
    "            stimulus = np.vstack([np.full((n_append, n_features), fill_value), stimulus])\n",
    "\n",
    "    # now reshape and lag\n",
    "    stimulus = np.reshape(stimulus, (-1, stim_samples_per_TR * n_features))\n",
    "\n",
    "    # offset by appending filler values\n",
    "    if offset_stim > 0:\n",
    "        stimulus = np.vstack([np.full((offset_TR, stim_samples_per_TR * n_features), fill_value), stimulus])\n",
    "\n",
    "    # check if lagging should be done\n",
    "    if lag_time != TR:\n",
    "        stimulus = make_lagged_stimulus(stimulus, lag_TR, fill_value=fill_value)            \n",
    " \n",
    "    return stimulus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "def make_X_Y(stimuli, fmri, TR, stim_TR, lag_time=6.0, start_times=None, offset_stim=0., fill_value=np.nan, remove_nans=True):\n",
    "    '''Creates (lagged) features and fMRI matrices concatenated along runs\n",
    "\n",
    "    Parameters\n",
    "\n",
    "        stimuli : list, list of stimulus representations\n",
    "        fmri : list, list of fMRI ndarrays\n",
    "        TR : int, float, repetition time of the fMRI data in seconds\n",
    "        stim_TR : int, float, repetition time of the stimulus in seconds\n",
    "        lag_time : int, float, optional,\n",
    "                   lag to introduce for stimuli in seconds,\n",
    "                   if no lagging should be done set this to TR\n",
    "        start_times : list, list of int, float, optional,\n",
    "                      starting time of the stimuli relative to fMRI recordings in seconds\n",
    "                      appends fill_value to stimulus representation to match fMRI and stimulus\n",
    "        offset_stim : int, float, optional,\n",
    "                      time to offset stimulus relative to fMRI in the lagged stimulus,\n",
    "                      i.e. when predicting fmri at time t use only stimulus features\n",
    "                      before t-offset_stim. This reduces the number of time points used\n",
    "                      in the model.\n",
    "        fill_value : int, float, or any valid numpy array element, optional,\n",
    "                     appends fill_value to stimulus array to account for starting_time\n",
    "                     use np.nan here with remove_nans=True to remove fmri/stimulus samples where no stimulus was presented\n",
    "        remove_nans : bool, bool or float 0<=remove_nans<=1, optional\n",
    "                      True/False indicate whether to remove all or none\n",
    "                      stimulus/fmri samples that contain nans\n",
    "                      a proportion keeps all samples in the lagged stimulus that have\n",
    "                      lower number of nans than this proportion.\n",
    "                      Replace nans with zeros in this case.\n",
    "\n",
    "    Returns:\n",
    "    tuple of two ndarrays,\n",
    "    the first element are the (lagged) stimuli,\n",
    "    the second element is the aligned fMRI data\n",
    "    '''\n",
    "    if len(stimuli) != len(fmri):\n",
    "        raise ValueError('Stimulus and fMRI need to have the same number of runs. '\n",
    "        'Instead fMRI has {} and stimulus {} runs.'.format(len(fmri), len(stimuli)))\n",
    "    n_features = stimuli[0].shape[1]\n",
    "    if not np.all(np.array([stim.shape[1] for stim in stimuli]) == n_features):\n",
    "        raise ValueError('Stimulus has different number of features per run.')\n",
    "\n",
    "    lagged_stimuli = []\n",
    "    aligned_fmri = []\n",
    "    for i, (stimulus, fmri_run) in enumerate(zip(stimuli, fmri)):\n",
    "        stimulus = generate_lagged_stimulus(\n",
    "            stimulus, fmri_run.shape[0], TR, stim_TR, lag_time=lag_time,\n",
    "            start_time=start_times[i] if start_times else 0.,\n",
    "            offset_stim=offset_stim, fill_value=fill_value)\n",
    "        # remove nans in stim/fmri here\n",
    "        if remove_nans:\n",
    "            remove_idx = get_remove_idx(stimulus, remove_nans)\n",
    "            stimulus = np.delete(stimulus, remove_idx, axis=0)\n",
    "            fmri_run = np.delete(fmri_run, remove_idx, axis=0)\n",
    "\n",
    "        # remove fmri samples recorded after stimulus has ended\n",
    "        if fmri_run.shape[0] != stimulus.shape[0]:\n",
    "            # check if the difference is due to offsetting and warn if it is not\n",
    "            if np.round(offset_stim/TR) < abs(fmri_run.shape[0] - stimulus.shape[0]):\n",
    "                warnings.warn('fMRI data and stimulus samples differ.'\n",
    "                ' Removing additional fMRI/stimulus samples. This could mean that you recorded '\n",
    "                'after stimulus ended, stopped recording early, or that something went wrong in the '\n",
    "                'preprocessing. fMRI: {}s stimulus: {}s'.format(\n",
    "                    TR*fmri_run.shape[0], TR*stimulus.shape[0]), RuntimeWarning)\n",
    "            if fmri_run.shape[0] > stimulus.shape[0]:\n",
    "                fmri_run = fmri_run[:-(fmri_run.shape[0]-stimulus.shape[0])]\n",
    "            else:\n",
    "                stimulus = stimulus[:-(stimulus.shape[0]-fmri_run.shape[0])]\n",
    "        lagged_stimuli.append(stimulus)\n",
    "        aligned_fmri.append(fmri_run)\n",
    "    return np.vstack(lagged_stimuli), np.vstack(aligned_fmri)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "`make_X_Y` allows you to align the (preprocessed) fMRI and stimulus data by specifying fMRI `TR` and stimulus `stim_TR`, as well as the `lag_time` (how long a stimulus window should be in seconds to predict a single fMRI TR) and potential stimulus offsets.\n",
    "Since we potentially want to preprocess and concatenate multiple runs, both `fmri` and `stimuli` are supposed to be lists. To process only a single run, you can use a list of one element.\n",
    "\n",
    "Let's look at an example, where the stimulus is sample every 100 ms and fMRI every 2s, i.e. every fMRI sample corresponds to 20 stimulus samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stim_TR, TR = 0.1, 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create a simulated `stimulus` object of 80 samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(80, 1)\n"
     ]
    }
   ],
   "source": [
    "stimulus = np.tile(np.arange(80)[:, None], (1, 1))\n",
    "print(stimulus.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And an according `fmri` object of 4 samples and one voxel (since we TRs differ)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 1)\n"
     ]
    }
   ],
   "source": [
    "fmri = np.tile(np.arange(0, 4)[:, None], (1, 1))\n",
    "print(fmri.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first align `fMRI` and `stimulus` without any offset or lag:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mboos/anaconda3/envs/mne/lib/python3.7/site-packages/ipykernel_launcher.py:56: RuntimeWarning: lag_time is None or equal to TR, no stimulus lagging will be done.\n"
     ]
    }
   ],
   "source": [
    "X, y = make_X_Y([stimulus], [fmri], TR, stim_TR, lag_time=None, offset_stim=0, start_times=[0])\n",
    "assert X.shape == (4, 20)\n",
    "assert y.shape == (4, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We keep the original number of samples in fMRI, but represent stimulus (and hence X) by the number of samples per fmri TR: stimulus thus becomes a (4, 20) array."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lagging the stimulus\n",
    "\n",
    "We can now call `make_X_Y` with the stimulus and fMRI TRs and a specified `lag_time`.\n",
    "Here we want to use 4 seconds of the stimulus to predict fMRI, but do not want to shift `fmri` relative to `stimulus` (`offset_stim` is 0.).\n",
    "This means that our encoding model can approximate a hemodynamic response function (HRF) by estimating a finite impulse response (FIR) that is 4 seconds long."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = make_X_Y([stimulus], [fmri], TR, stim_TR, lag_time=4, offset_stim=0, start_times=[0])\n",
    "assert X.shape == (3, 40)\n",
    "assert y.shape == (3, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shifting the stimulus\n",
    "\n",
    "We could also shift `fmri` relative to `stimulus`, to account for the delayed onset of the hemodynamic response - this is different than estimating the hemodynamic response from the window given by `lag_time`.\n",
    "In practice this means we estimate an hemodynamic response function (HRF) by a FIR in the time period from -6s to -2s before each fMRI sample. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = make_X_Y([stimulus], [fmri], TR, stim_TR, lag_time=4, offset_stim=2, start_times=[0])\n",
    "assert X.shape == (2, 40)\n",
    "assert y.shape == (2, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Handling out-of-recording data\n",
    "\n",
    "Because of our shift we \"lose\" one sample, because by default `fill_value` fills values that lie outside the recording interval by NaNs and by default `remove_nans` specifies that all samples with NaNs are dropped.\n",
    "\n",
    "To check that behavior, we see what we get when we don't remove NaNs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X, y = make_X_Y([stimulus], [fmri], TR, stim_TR, lag_time=4, offset_stim=2, start_times=[0], remove_nans=False)\n",
    "assert X.shape == (4, 40)\n",
    "assert y.shape == (4, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We keep the original number of samples, but some are filled with NaNs now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan\n",
      "  nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan\n",
      "  nan nan nan nan]\n",
      " [ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16. 17.\n",
      "  18. 19. nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan nan\n",
      "  nan nan nan nan]\n",
      " [20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37.\n",
      "  38. 39.  0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15.\n",
      "  16. 17. 18. 19.]\n",
      " [40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57.\n",
      "  58. 59. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35.\n",
      "  36. 37. 38. 39.]]\n"
     ]
    }
   ],
   "source": [
    "assert np.isnan(X).sum() == 60\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the first samples completely consists of NaNs, because by lagging and offsetting we assume that the   fMRI sample at time point t can be predicted by the time period in the stimulus of t-6s to t-2s.\n",
    "However, we don't have any stimulus presented in that time!\n",
    "In the second sample we can see that the first half of the stimulus still consists of NaNs: that's because for t=2s, the time period in the stimulus from t-6s to t-2s has only data for t=0s but not t=4s.\n",
    "Keep in mind that the stimulus at t=0s corresponds to the first 2s of the stimulus (because we reshaped the stimulus TR to correspond to the 2s fmri TR)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "<h4 id=\"generate_lagged_stimulus\" class=\"doc_header\"><code>generate_lagged_stimulus</code><a href=\"__main__.py#L2\" class=\"source_link\" style=\"float:right\">[source]</a></h4>\n",
       "\n",
       "> <code>generate_lagged_stimulus</code>(**`stimulus`**, **`fmri_samples`**, **`TR`**, **`stim_TR`**, **`lag_time`**=*`None`*, **`start_time`**=*`0.0`*, **`offset_stim`**=*`0.0`*, **`fill_value`**=*`nan`*)\n",
       "\n",
       "Generates a lagged stimulus representation temporally aligned with the fMRI data\n",
       "\n",
       "Parameters\n",
       "\n",
       "    stimuli : ndarray, stimulus representation of shape (samples, features)\n",
       "    fmri_samples : int, samples of corresponding fmri run\n",
       "    TR : int, float, repetition time of the fMRI data in seconds\n",
       "    stim_TR : int, float, repetition time of the stimulus in seconds\n",
       "    lag_time : int, float, or None, optional,\n",
       "           lag to introduce for stimuli in seconds,\n",
       "           if no lagging should be done set this to TR or None\n",
       "    start_time :  int, float, optional, default 0.\n",
       "              starting time of the stimulus relative to fMRI recordings in seconds\n",
       "              appends fill_value to stimulus representation to match fMRI and stimulus\n",
       "    offset_stim : int, float, optional, default 0.\n",
       "              time to offset stimulus relative to fMRI in the lagged stimulus,\n",
       "              i.e. when predicting fmri at time t use only stimulus features\n",
       "              before t-offset_stim. This reduces the number of time points used\n",
       "              in the model.\n",
       "    fill_value : int, float, or any valid numpy array element, optional, default np.nan\n",
       "             appends fill_value to stimulus array to account for starting_time\n",
       "             use np.nan here with remove_nans=True to remove fmri/stimulus samples where no stimulus was presented\n",
       "\n",
       "Returns:\n",
       "    ndarray of the lagged stimulus of shape (samples, lagged features)"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_doc(generate_lagged_stimulus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`generate_lagged_stimulus` takes care of aligning fMRI and stimulus data, it is used internally by `make_X_Y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "<h4 id=\"get_remove_idx\" class=\"doc_header\"><code>get_remove_idx</code><a href=\"__main__.py#L2\" class=\"source_link\" style=\"float:right\">[source]</a></h4>\n",
       "\n",
       "> <code>get_remove_idx</code>(**`lagged_stimulus`**, **`remove_nan`**=*`True`*)\n",
       "\n",
       "Returns indices of rows in lagged_stimulus to remove"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show_doc(get_remove_idx)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (mne)",
   "language": "python",
   "name": "mne"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}