encoding.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# default_exp encoding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#hide\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "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 numpy as np\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.linear_model import RidgeCV\n",
    "import warnings\n",
    "import copy\n",
    "from joblib import Parallel, delayed\n",
    "from sklearn.multioutput import MultiOutputRegressor, _fit_estimator\n",
    "from sklearn.utils import check_X_y, check_array\n",
    "from sklearn.utils.validation import check_is_fitted\n",
    "from sklearn.base import RegressorMixin\n",
    "\n",
    "def product_moment_corr(x,y):\n",
    "    '''Product-moment correlation for two ndarrays x, y'''\n",
    "    from sklearn.preprocessing import StandardScaler\n",
    "    x = StandardScaler().fit_transform(x)\n",
    "    y = StandardScaler().fit_transform(y)\n",
    "    n = x.shape[0]\n",
    "    r = (1/(n-1))*(x*y).sum(axis=0)\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training and validating voxel-wise encoding models\n",
    "> Functions for training independent Ridge regressions for a large number of voxels and validating their performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "\n",
    "def get_model_plus_scores(X, y, estimator=None, cv=None, scorer=None,\n",
    "                          voxel_selection=True, validate=True, **kwargs):\n",
    "    '''Returns multiple estimator trained in a cross-validation on n_splits of the data and scores on the left-out folds\n",
    "\n",
    "    Parameters\n",
    "\n",
    "        X : ndarray of shape (samples, features)\n",
    "        y : ndarray of shape (samples, targets)\n",
    "        estimator : None or estimator object that implements fit and predict\n",
    "                    if None, uses RidgeCV per default\n",
    "        cv : int, None, or a cross-validation object that implements a split method, default is None, optional.\n",
    "             int specifies the number of cross-validation splits of a KFold cross validation\n",
    "             None defaults to a scikit-learn KFold cross-validation with default settings\n",
    "             a scikit-learn-like cross-validation object needs to implement a split method for X and y\n",
    "        scorer : None or any sci-kit learn compatible scoring function, optional\n",
    "                 default uses product moment correlation\n",
    "        voxel_selection : bool, optional, default True\n",
    "                          Whether to only use voxels with variance larger than zero.\n",
    "                          This will set scores for these voxels to zero.\n",
    "        validate : bool, optional, default True\n",
    "                     Whether to validate the model via cross-validation\n",
    "                     or to just train the estimator\n",
    "                     if False, scores will be computed on the training set\n",
    "        kwargs : additional parameters that will be used to initialize RidgeCV if estimator is None \n",
    "    Returns\n",
    "        tuple of n_splits estimators trained on training folds or single estimator if validation is False\n",
    "        and scores for all concatenated out-of-fold predictions'''\n",
    "    from sklearn.utils.estimator_checks import check_regressor_multioutput\n",
    "    if scorer is None:\n",
    "        scorer = product_moment_corr\n",
    "    if cv is None:\n",
    "        cv = KFold()\n",
    "    if isinstance(cv, int):\n",
    "        cv = KFold(n_splits=cv)\n",
    "    models = []\n",
    "    score_list = []\n",
    "    if estimator is None:\n",
    "        estimator = RidgeCV(**kwargs)\n",
    "        \n",
    "    if voxel_selection:\n",
    "        voxel_var = np.var(y, axis=0)\n",
    "        y = y[:, voxel_var > 0.]\n",
    "    if validate:\n",
    "        for train, test in cv.split(X, y):\n",
    "            models.append(copy.deepcopy(estimator).fit(X[train], y[train]))\n",
    "            if voxel_selection:\n",
    "                scores = np.zeros_like(voxel_var)\n",
    "                scores[voxel_var > 0.] =  scorer(y[test], models[-1].predict(X[test]))\n",
    "            else:\n",
    "                scores = scorer(y[test], models[-1].predict(X[test]))\n",
    "            score_list.append(scores[:, None])\n",
    "        score_list = np.concatenate(score_list, axis=-1)\n",
    "    else:\n",
    "        models = estimator.fit(X, y)\n",
    "        score_list = scorer(y, estimator.predict(X))\n",
    "    return models, score_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`get_model_plus_scores` is a convenience function that trains multiple Ridge regressions in a cross-validation scheme and evaluates their performance on the respective test set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples\n",
    "\n",
    "First, we create some simulated `stimulus` and `fmri` data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stimulus = np.random.randn(1000, 5)\n",
    "fmri = np.random.randn(1000, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the default Ridge regression\n",
    "\n",
    "We can now use `get_model_plus_scores` to estimate multiple [RidgeCV](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RidgeCV.html) regressions, one for each voxel (that maps the stimulus representation to this voxel) and one for each split (trained on a different training set and evaluated on the held-out set).\n",
    "Since sklearn's `RidgeCV` estimator allows multi-output, we get one `RidgeCV` object per split."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[RidgeCV(alphas=array([ 0.1,  1. , 10. ])),\n",
       " RidgeCV(alphas=array([ 0.1,  1. , 10. ])),\n",
       " RidgeCV(alphas=array([ 0.1,  1. , 10. ]))]"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridges, scores = get_model_plus_scores(stimulus, fmri, cv=3)\n",
    "assert len(ridges) == 3\n",
    "ridges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each `RidgeCV` estimator maps from the feature space to each voxel.\n",
    "In our example, that means it has 10 (the number of voxels-9 independently trained regression models with 5 coeficients each (the number of features)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.01930266  0.0350985  -0.04548384 -0.01058159 -0.07382483]\n",
      " [ 0.00183012 -0.00830046 -0.03604675 -0.00016843  0.03161116]\n",
      " [-0.04032306  0.01782385  0.02112695  0.01673908 -0.00645515]\n",
      " [ 0.0273047  -0.02382577 -0.06169262  0.06232742 -0.03331368]\n",
      " [ 0.01294108 -0.04825337 -0.04646228 -0.04701512 -0.00017405]\n",
      " [ 0.02008884 -0.07065883  0.01958404 -0.04115758 -0.02967363]\n",
      " [ 0.00502653 -0.02164034 -0.00419562 -0.05675778  0.00716245]\n",
      " [ 0.0080379   0.03230623  0.01527909 -0.02469508 -0.01681562]\n",
      " [ 0.01363082  0.02686557 -0.05923971  0.01392573 -0.00945206]\n",
      " [ 0.01665226 -0.01499506 -0.0043113  -0.01658976  0.06103525]]\n"
     ]
    }
   ],
   "source": [
    "assert ridges[0].coef_.shape == (10, 5)\n",
    "print(ridges[0].coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also get a set of scores (by default the [product moment correlation](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient), but you can supply your own via the `scorer` argument) that specifies how well we predict left-out data (with the usual caveats of using a correlation coefficient for evaluating it). In our case it is of shape (10, 3) because we predict 10 voxels and use a 3-fold cross-validation, i.e. we split 3 times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.11195214,  0.10260332,  0.02153565],\n",
       "       [-0.06450754, -0.07106461, -0.0447989 ],\n",
       "       [ 0.00559572, -0.03221425,  0.02866726],\n",
       "       [-0.04101258, -0.02197306, -0.04277958],\n",
       "       [ 0.02352969,  0.02008923, -0.02062713],\n",
       "       [ 0.01027339,  0.03074076,  0.01248573],\n",
       "       [-0.05974497, -0.03980094, -0.11293944],\n",
       "       [-0.01607721, -0.02264425, -0.07340733],\n",
       "       [-0.06009815, -0.05553956,  0.02102434],\n",
       "       [ 0.02388894, -0.01513094,  0.0904367 ]])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "assert scores.shape == (10, 3)\n",
    "scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also change the parameters of the `RidgeCV` function.\n",
    "For example, we can use pre-specified hyperparameters, like the values of the regularization parameter $\\alpha$ we want to perform a gridsearch over or whether we want to normalize features. If we want to use other parameters for the default `RidgeCV`, we can just pass the parameters as additional keyword arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphas = [100]\n",
    "ridges, scores = get_model_plus_scores(stimulus, fmri, alphas=alphas,\n",
    "                                       normalize=True, alpha_per_target=True)\n",
    "assert ridges[0].normalize\n",
    "assert ridges[0].alphas.shape == (1,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using your own estimator\n",
    "\n",
    "\n",
    "Additionally, we can use any other estimator that implements `fit` and `predict`.\n",
    "For example, we can use [CCA](https://scikit-learn.org/stable/modules/generated/sklearn.cross_decomposition.CCA.html) as an encoding model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import cross_decomposition\n",
    "\n",
    "our_estimator = cross_decomposition.CCA(n_components=2)\n",
    "\n",
    "ccas, scores = get_model_plus_scores(stimulus, fmri, our_estimator,\n",
    "                                     cv=3)\n",
    "assert type(ccas[0]) == cross_decomposition._pls.CCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your favorite estimator does not work in the multioutput regime, i.e. it cannot predict multiple targets/voxels, then `get_model_plus_scores` will wrap it into sklearn's [MultiOutputRegressor](https://scikit-learn.org/stable/modules/generated/sklearn.multioutput.MultiOutputRegressor.html) by default. However, for many voxels this can increase training time by a lot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[MultiOutputRegressor(estimator=Lasso()),\n",
       " MultiOutputRegressor(estimator=Lasso()),\n",
       " MultiOutputRegressor(estimator=Lasso())]"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.multioutput import MultiOutputRegressor\n",
    "\n",
    "our_estimator = MultiOutputRegressor(Lasso())\n",
    "\n",
    "lassos, scores = get_model_plus_scores(stimulus, fmri, our_estimator,\n",
    "                                       cv=3)\n",
    "lassos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training without validation\n",
    "\n",
    "We can also train an estimator without any validation, if, for example we want to test on a different dataset. In that case, the scores will be computed with the trained estimator on the training set, i.e. they will contain no information about the generalization performance of the estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "our_estimator = RidgeCV()\n",
    "\n",
    "model, scores = get_model_plus_scores(stimulus, fmri, our_estimator,\n",
    "                                       validate=False)\n",
    "assert type(model) == RidgeCV\n",
    "assert scores.shape == (10,)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using your own cross-validation method\n",
    "\n",
    "Instead of the default `KFold` cross-validation without shuffling, we can also use any sckit-learn compatible cross-validation iterators (e.g. [these](https://scikit-learn.org/stable/modules/cross_validation.html#cross-validation-iterators)).\n",
    "For example, we could use a `TimeSerisSplit` to test our predictions on only the most recent part of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import TimeSeriesSplit\n",
    "\n",
    "ts_cv = TimeSeriesSplit(n_splits=5)\n",
    "\n",
    "model, scores = get_model_plus_scores(stimulus, fmri, cv=ts_cv)\n",
    "assert scores.shape == (10, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Distributed training\n",
    "\n",
    "Voxel-wise encoding models can take a long time and a lot of memory to train, especially if we use the full brain or high resolution fMRI data.\n",
    "\n",
    "The `BlockMultiOutput` class can help distribute the load across multiple cores by splitting the fMRI data into multiple \"blocks\" (the `n_blocks` parameter) and training an estimator for each block.\n",
    "Without parallelization, this class allows one to train voxel-wise encoding models, even if training a single, large estimator takes up too much memory, by training the estimator for blocks of your data independently.\n",
    "\n",
    "This works even if the original fMRI data do not fit into memory, by using a [memmapped](https://numpy.org/doc/stable/reference/generated/numpy.memmap.html) Numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#export\n",
    "\n",
    "class BlockMultiOutput(MultiOutputRegressor, RegressorMixin):\n",
    "    \"\"\"Multi target regression with block-wise fit\n",
    "    This strategy consists of splitting the targets in blocks and fitting one regressor per block.\n",
    "    The estimator used needs to natively support multioutput.\n",
    "    \n",
    "    Parameters\n",
    "\n",
    "        estimator : estimator object\n",
    "            An estimator object implementing `fit` and `predict` and supporting multioutput.\n",
    "        n_blocks : int, optional, default=10\n",
    "            The number of blocks for the target variable.\n",
    "            This is a split along *targets* (columns of the array), not observations (rows of the array).\n",
    "        n_jobs : int, optional, default=1\n",
    "            The number of jobs to run in parallel for `fit`. If -1,\n",
    "            then the number of jobs is set to the number of cores.\n",
    "            When individual estimators are fast to train or predict\n",
    "            using `n_jobs>1` can result in slower performance due\n",
    "            to the overhead of spawning processes.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, estimator, n_blocks=10, n_jobs=1):\n",
    "        self.estimator = estimator\n",
    "        self.n_blocks = n_blocks\n",
    "        self.n_jobs = n_jobs\n",
    "\n",
    "    def fit(self, X, y, sample_weight=None):\n",
    "        \"\"\" Fit the model to data.\n",
    "        Fit a separate model for each chunk of output.\n",
    "        \n",
    "        Parameters\n",
    "        \n",
    "            X : (sparse) array-like, shape (n_samples, n_features)\n",
    "                Data.\n",
    "            y : (sparse) array-like, shape (n_samples, n_outputs)\n",
    "                Multi-output targets. An indicator matrix turns on multilabel\n",
    "                estimation.\n",
    "            sample_weight : array-like, shape = (n_samples) or None\n",
    "                Sample weights. If None, then samples are equally weighted.\n",
    "                Only supported if the underlying regressor supports sample\n",
    "                weights.\n",
    "                \n",
    "        Returns\n",
    "        \n",
    "            self : object\n",
    "                Returns self\n",
    "        \"\"\"\n",
    "        if not hasattr(self.estimator, \"fit\"):\n",
    "            raise ValueError(\"The base estimator should implement a fit method\")\n",
    "\n",
    "        if y.ndim == 1:\n",
    "            raise ValueError(\"y must have at least two dimensions for \"\n",
    "                             \"multi-output regression but has only one.\")\n",
    "\n",
    "        if (sample_weight is not None and\n",
    "                not has_fit_parameter(self.estimator, 'sample_weight')):\n",
    "            raise ValueError(\"Underlying estimator does not support\"\n",
    "                             \" sample weights.\")\n",
    "        kfold = KFold(n_splits=self.n_blocks)\n",
    "        smpl_X, smpl_y = np.zeros((y.shape[1],1)), np.zeros((y.shape[1],1))\n",
    "        self.estimators_ = Parallel(n_jobs=self.n_jobs)(\n",
    "            delayed(_fit_estimator)(\n",
    "                self.estimator, X, y[:, block], sample_weight)\n",
    "            for _, block in kfold.split(smpl_X, smpl_y))\n",
    "        return self\n",
    "\n",
    "    def partial_predict(self, X):\n",
    "        \"\"\"Predict multi-output variable using a model\n",
    "         trained for each target variable block and yields predictions for each block as an iterator.\n",
    "         \n",
    "        Parameters\n",
    "        \n",
    "        X : (sparse) array-like, shape (n_samples, n_features)\n",
    "            Data.\n",
    "            \n",
    "        Returns\n",
    "        \n",
    "        y : (sparse) array-like, shape (n_samples, n_outputs)\n",
    "            Multi-output targets predicted across multiple predictors.\n",
    "            Note: Separate models are generated for each predictor.\n",
    "        \"\"\"\n",
    "        check_is_fitted(self, 'estimators_')\n",
    "        if not hasattr(self.estimator, \"predict\"):\n",
    "            raise ValueError(\"The base estimator should implement a predict method\")\n",
    "\n",
    "        X = check_array(X, accept_sparse=True)\n",
    "\n",
    "        for estimator in self.estimators_:\n",
    "            yield estimator.predict(X)\n",
    "\n",
    "    def predict(self, X):\n",
    "        \"\"\"Predict multi-output variable using a model\n",
    "         trained for each target variable block.\n",
    "         \n",
    "        Parameters\n",
    "        \n",
    "            X : (sparse) array-like, shape (n_samples, n_features)\n",
    "                Data.\n",
    "                \n",
    "        Returns\n",
    "        \n",
    "            y : (sparse) array-like, shape (n_samples, n_outputs)\n",
    "                Multi-output targets predicted across multiple predictors.\n",
    "                Note: Separate models are generated for each predictor.\n",
    "        \"\"\"\n",
    "        check_is_fitted(self, 'estimators_')\n",
    "        if not hasattr(self.estimator, \"predict\"):\n",
    "            raise ValueError(\"The base estimator should implement a predict method\")\n",
    "\n",
    "        X = check_array(X, accept_sparse=True)\n",
    "\n",
    "        y = Parallel(n_jobs=self.n_jobs)(\n",
    "            delayed(e.predict)(X)\n",
    "            for e in self.estimators_)\n",
    "\n",
    "        return np.hstack(y)\n",
    "\n",
    "    def score(self, X, y):\n",
    "        \"\"\"Returns the correlation of the prediction with the target for each output.\n",
    "        \n",
    "        Parameters\n",
    "        \n",
    "            X : array-like, shape (n_samples, n_features)\n",
    "                Test samples.\n",
    "            y : array-like, shape (n_samples) or (n_samples, n_outputs)\n",
    "                True values for X.\n",
    "                \n",
    "        Returns\n",
    "        \n",
    "            score : float\n",
    "                Correlation of self.predict(X) wrt. y.\n",
    "        \"\"\"\n",
    "        from sklearn.preprocessing import StandardScaler\n",
    "        from itertools import izip\n",
    "        kfold = KFold(n_splits=self.n_blocks)\n",
    "        smpl_X, smpl_y = np.zeros((y.shape[1],1)), np.zeros((y.shape[1],1))\n",
    "        scores = []\n",
    "        for prediction, (_, block) in izip(self.partial_predict(X), kfold.split(smpl_X, smpl_y)):\n",
    "            mx = StandardScaler().fit_transform(prediction).astype('float32')\n",
    "            my = StandardScaler().fit_transform(y[:, block]).astype('float32')\n",
    "            n = mx.shape[0]\n",
    "            r = (1/(n-1))*(mx*my).sum(axis=0)\n",
    "            scores.append(r)\n",
    "        return np.concatenate(scores)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "Let's generate a larger data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "our_estimator = BlockMultiOutput(RidgeCV(alphas=[10,100]))\n",
    "\n",
    "estimators, scores = get_model_plus_scores(stimulus, fmri, our_estimator,\n",
    "                                     cv=3)\n",
    "assert len(estimators) == 3\n",
    "assert len(estimators[0].estimator.alphas) == 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each `BlockMultiOutput` estimator contains `n_blocks` estimators that are trained on different blocks of the target."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100])),\n",
       " RidgeCV(alphas=array([ 10, 100]))]"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "assert len(estimators[0].estimators_) == estimators[0].n_blocks\n",
    "estimators[0].estimators_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If `fmri` is of shape $(n\\_samples, n\\_targets)$, each of the `n_blocks` estimators in `BlockMultiOutput.estimators_` will contain the coefficients for ${n\\_targets}/{n\\_blocks}$ targets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fmri shape: (1000, 10) \n",
      "n_blocks: 10 \n",
      "coefficients of the estimator for one block: (1, 5)\n"
     ]
    }
   ],
   "source": [
    "assert estimators[0].estimators_[0].coef_.shape == (1, 5)\n",
    "print('fmri shape: {} \\nn_blocks: {} \\n'\n",
    "      'coefficients of the estimator for one block: {}'.format(\n",
    "          fmri.shape, our_estimator.n_blocks, estimators[0].estimators_[0].coef_.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `MultiBlockOutput` instance normally to predict data, i.e. it produces predictions of the full fmri data by concatenating the predictions of every block-estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 10)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "assert estimators[0].predict(stimulus).shape == (1000, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parallelizing voxel-wise encoding models\n",
    "\n",
    "We can use this to parallelize encoding models as well, by specifying the `n_jobs` parameter.\n",
    "Keep in mind that this requires copying the full `stimulus` data to every worker and can thus increase memory demand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "our_estimator = BlockMultiOutput(RidgeCV(alphas=[10,100]), n_jobs=10)\n",
    "\n",
    "estimators, scores = get_model_plus_scores(stimulus, fmri, our_estimator,\n",
    "                                     cv=3)\n",
    "assert len(estimators) == 3\n",
    "assert estimators[0].n_jobs == 10"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (mne)",
   "language": "python",
   "name": "mne"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}