examples/reprsimil/bayesian_rsa_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This demo shows how to use the Bayesian Representational Similarity Analysis method in brainiak with a simulated dataset.\n",
    "### The brainik.reprsimil.brsa module has two estimators named BRSA and GBRSA. Both of them can be used to estimate representational similarity from a single participant, but with some differences in the assumptions of the models and fitting procedure. The basic usages are similar. We now generally recommend using GBRSA over BRSA for most of the cases. This document shows how to use BRSA for most of the part. At the end of the document, the usage of GBRSA is shown as well. You are encouranged to go through the example and try both estimators for your data.\n",
    "### The group_brsa_example.ipynb in the same directory demonstrates how to use GBRSA to estimate shared representational structure from multiple participants."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that the model assumes that the covariance matrix U which all $\\beta_i$ follow describe a multi-variate Gaussian distribution that is zero-meaned. This assumption does not imply that there must be both positive and negative responses across voxels.\n",
    "However, it means that (Group) Bayesian RSA treats the task-evoked activity against baseline BOLD level as signal, while in other RSA tools the deviation of task-evoked activity in each voxel from the average task-evoked activity level across voxels may be considered as signal of interest.\n",
    "Due to this assumption in (G)BRSA, relatively high degree of similarity may be expected when the activity patterns of two task conditions share a strong sensory driven components. When two task conditions elicit exactly the same activity pattern but only differ in their global magnitudes, under the assumption in (G)BRSA, their similarity is 1; under the assumption that only deviation of pattern from average patterns is signal of interest (which is currently not supported by (G)BRSA), their similarity would be -1 because the deviations of the two patterns from their average pattern are exactly opposite."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load some package which we will use in this demo.\n",
    "If you see error related to loading any package, you can install that package. For example, if you use Anaconda, you can use \"conda install matplotlib\" to install matplotlib.\n",
    "##### Notice that due to current implementation, you need to import either prior_GP_var_inv_gamma or prior_GP_var_half_cauchy from brsa module, in order to use the smooth prior imposed onto SNR in BRSA (see below). They are forms of priors imposed on the variance of Gaussian Process prior on log(SNR). (If you think these sentences are confusing, just import them like below and forget about this). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scipy.stats\n",
    "import scipy.spatial.distance as spdist\n",
    "import numpy as np\n",
    "from brainiak.reprsimil.brsa import BRSA, prior_GP_var_inv_gamma, prior_GP_var_half_cauchy\n",
    "from brainiak.reprsimil.brsa import GBRSA\n",
    "import brainiak.utils.utils as utils\n",
    "import matplotlib.pyplot as plt\n",
    "import logging\n",
    "np.random.seed(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### You might want to keep a log of the output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logging.basicConfig(\n",
    "    level=logging.DEBUG,\n",
    "    filename='brsa_example.log',\n",
    "    format='%(relativeCreated)6d %(threadName)s %(message)s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# We want to simulate some data in which each voxel responds to different task conditions differently, but following a common covariance structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load an example design matrix.\n",
    "The user should prepare their design matrix with their favorate software, such as using 3ddeconvolve of AFNI, or using SPM or FSL.\n",
    "The design matrix reflects your belief of how fMRI signal should respond to a task (if a voxel does respond).\n",
    "The common assumption is that a neural event that you are interested in will elicit a slow hemodynamic response in some voxels. The response peaks around 4-6 seconds after the event onset and dies down more than 12 seconds after the event. Therefore, typically you convolve a time series A, composed of delta (stem) functions reflecting the time of each neural event belonging to the same category (e.g. all trials in which a participant sees a face), with a hemodynamic response function B, to form the hypothetic response of any voxel to such type of neural event.\n",
    "For each type of event, such a convoluted time course can be generated. These time courses, put together, are called design matrix, reflecting what we believe a temporal signal would look like, if it exists in any voxel.\n",
    "Our goal is to figure out how the (spatial) response pattern of a population of voxels (in an Region of Interest, ROI) are similar or disimilar to different types of tasks (e.g., watching face vs. house, watching different categories of animals, different conditions of a cognitive task). So we need the design matrix in order to estimate the similarity matrix we are interested.\n",
    "\n",
    "We can use the utility called ReadDesign in brainiak.utils to read a design matrix generated from AFNI. For design matrix saved as Matlab data file by SPM or or other toolbox, you can use scipy.io.loadmat('YOURFILENAME') and extract the design matrix from the dictionary returned. Basically, the Bayesian RSA in this toolkit just needs a numpy array which is in size of {time points} * {condition}\n",
    "You can also generate design matrix using the function gen_design which is in brainiak.utils. It takes in (name of) event timing files in AFNI or FSL format (denoting onsets, duration, and weight for each event belongning to the same condition) and outputs the design matrix as numpy array.\n",
    "\n",
    "In typical fMRI analysis, some nuisance regressors such as head motion, baseline time series and slow drift are also entered into regression. In using our method, you should not include such regressors into the design matrix, because the spatial spread of such nuisance regressors might be quite different from the spatial spread of task related signal. Including such nuisance regressors in design matrix might influence the pseudo-SNR map, which in turn influence the estimation of the shared covariance matrix. \n",
    "\n",
    "### We concatenate the design matrix by 2 times, mimicking 2 runs of identical timing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "design = utils.ReadDesign(fname=\"example_design.1D\")\n",
    "\n",
    "n_run = 3\n",
    "design.n_TR = design.n_TR * n_run\n",
    "design.design_task = np.tile(design.design_task[:,:-1],\n",
    "                             [n_run, 1])\n",
    "# The last \"condition\" in design matrix\n",
    "# codes for trials subjects made and error.\n",
    "# We ignore it here.\n",
    "\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(12, 3),\n",
    "                 dpi=150, facecolor='w', edgecolor='k')\n",
    "plt.plot(design.design_task)\n",
    "plt.ylim([-0.2, 0.4])\n",
    "plt.title('hypothetic fMRI response time courses '\n",
    "          'of all conditions\\n'\n",
    "         '(design matrix)')\n",
    "plt.xlabel('time')\n",
    "plt.show()\n",
    "\n",
    "n_C = np.size(design.design_task, axis=1)\n",
    "# The total number of conditions.\n",
    "ROI_edge = 15\n",
    "# We simulate \"ROI\" of a rectangular shape\n",
    "n_V = ROI_edge**2 * 2\n",
    "# The total number of simulated voxels\n",
    "n_T = design.n_TR\n",
    "# The total number of time points,\n",
    "# after concatenating all fMRI runs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## simulate data: noise + signal\n",
    "### First, we start with noise, which is Gaussian Process in space and AR(1) in time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_bot = 0.5\n",
    "noise_top = 5.0\n",
    "noise_level = np.random.rand(n_V) * \\\n",
    "    (noise_top - noise_bot) + noise_bot\n",
    "# The standard deviation of the noise is in the range of [noise_bot, noise_top]\n",
    "# In fact, we simulate autocorrelated noise with AR(1) model. So the noise_level reflects\n",
    "# the independent additive noise at each time point (the \"fresh\" noise)\n",
    "\n",
    "# AR(1) coefficient\n",
    "rho1_top = 0.8\n",
    "rho1_bot = -0.2\n",
    "rho1 = np.random.rand(n_V) \\\n",
    "    * (rho1_top - rho1_bot) + rho1_bot\n",
    "\n",
    "noise_smooth_width = 10.0\n",
    "coords = np.mgrid[0:ROI_edge, 0:ROI_edge*2, 0:1]\n",
    "coords_flat = np.reshape(coords,[3, n_V]).T\n",
    "dist2 = spdist.squareform(spdist.pdist(coords_flat, 'sqeuclidean'))\n",
    "\n",
    "# generating noise\n",
    "K_noise = noise_level[:, np.newaxis] \\\n",
    "    * (np.exp(-dist2 / noise_smooth_width**2 / 2.0) \\\n",
    "       + np.eye(n_V) * 0.1) * noise_level\n",
    "# We make spatially correlated noise by generating\n",
    "# noise at each time point from a Gaussian Process\n",
    "# defined over the coordinates.\n",
    "plt.pcolor(K_noise)\n",
    "plt.colorbar()\n",
    "plt.xlim([0, n_V])\n",
    "plt.ylim([0, n_V])\n",
    "plt.title('Spatial covariance matrix of noise')\n",
    "plt.show()\n",
    "L_noise = np.linalg.cholesky(K_noise)\n",
    "noise = np.zeros([n_T, n_V])\n",
    "noise[0, :] = np.dot(L_noise, np.random.randn(n_V))\\\n",
    "    / np.sqrt(1 - rho1**2)\n",
    "for i_t in range(1, n_T):\n",
    "    noise[i_t, :] = noise[i_t - 1, :] * rho1 \\\n",
    "        + np.dot(L_noise,np.random.randn(n_V))\n",
    "# For each voxel, the noise follows AR(1) process:\n",
    "# fresh noise plus a dampened version of noise at\n",
    "# the previous time point.\n",
    "# In this simulation, we also introduced spatial smoothness resembling a Gaussian Process.\n",
    "# Notice that we simulated in this way only to introduce spatial noise correlation.\n",
    "# This does not represent the assumption of the form of spatial noise correlation in the model.\n",
    "# Instead, the model is designed to capture structured noise correlation manifested\n",
    "# as a few spatial maps each modulated by a time course, which appears as spatial noise correlation. \n",
    "fig = plt.figure(num=None, figsize=(12, 2), dpi=150,\n",
    "                 facecolor='w', edgecolor='k')\n",
    "plt.plot(noise[:, 0])\n",
    "plt.title('noise in an example voxel')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Then, we simulate signals, assuming the magnitude of response to each condition follows a common covariance matrix. \n",
    "#### Our model allows to impose a Gaussian Process prior on the log(SNR) of each voxels. \n",
    "What this means is that SNR turn to be smooth and local, but betas (response amplitudes of each voxel to each condition) are not necessarily correlated in space. Intuitively, this is based on the assumption that voxels coding for related aspects of a task turn to be clustered (instead of isolated)\n",
    "\n",
    "Our Gaussian Process are defined on both the coordinate of a voxel and its mean intensity.\n",
    "This means that voxels close together AND have similar intensity should have similar SNR level. Therefore, voxels of white matter but adjacent to gray matters do not necessarily have high SNR level.\n",
    "\n",
    "If you have an ROI saved as a binary Nifti file, say, with name 'ROI.nii'\n",
    "Then you can use nibabel package to load the ROI and the following example code to retrive the coordinates of voxels.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Note: the following code won't work if you just installed Brainiak and try this demo because ROI.nii does not exist. It just serves as an example for you to retrieve coordinates of voxels in an ROI. You can use the ROI_coords for the argument coords in BRSA.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import nibabel\n",
    "# ROI = nibabel.load('ROI.nii')\n",
    "# I,J,K = ROI.shape \n",
    "# all_coords = np.zeros((I, J, K, 3)) \n",
    "# all_coords[...,0] = np.arange(I)[:, np.newaxis, np.newaxis] \n",
    "# all_coords[...,1] = np.arange(J)[np.newaxis, :, np.newaxis] \n",
    "# all_coords[...,2] = np.arange(K)[np.newaxis, np.newaxis, :] \n",
    "# ROI_coords = nibabel.affines.apply_affine(\n",
    "#     ROI.affine, all_coords[ROI.get_data().astype(bool)])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Let's keep in mind of the pattern of the ideal covariance / correlation below and see how well BRSA can recover their patterns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ideal covariance matrix\n",
    "ideal_cov = np.zeros([n_C, n_C])\n",
    "ideal_cov = np.eye(n_C) * 0.6\n",
    "ideal_cov[8:12, 8:12] = 0.6\n",
    "for cond in range(8, 12):\n",
    "    ideal_cov[cond,cond] = 1\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
    "plt.pcolor(ideal_cov)\n",
    "plt.colorbar()\n",
    "plt.xlim([0, 16])\n",
    "plt.ylim([0, 16])\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('ideal covariance matrix')\n",
    "plt.show()\n",
    "\n",
    "std_diag = np.diag(ideal_cov)**0.5\n",
    "ideal_corr = ideal_cov / std_diag / std_diag[:, None]\n",
    "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
    "plt.pcolor(ideal_corr)\n",
    "plt.colorbar()\n",
    "plt.xlim([0, 16])\n",
    "plt.ylim([0, 16])\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('ideal correlation matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In the following, pseudo-SNR is generated from a Gaussian Process defined on a \"rectangular\" ROI, just for simplicity of code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_full = np.linalg.cholesky(ideal_cov)        \n",
    "\n",
    "# generating signal\n",
    "snr_level = 1.0\n",
    "# Notice that accurately speaking this is not SNR.\n",
    "# The magnitude of signal depends not only on beta but also on x.\n",
    "# (noise_level*snr_level)**2 is the factor multiplied\n",
    "# with ideal_cov to form the covariance matrix from which\n",
    "# the response amplitudes (beta) of a voxel are drawn from.\n",
    "\n",
    "tau = 1.0\n",
    "# magnitude of Gaussian Process from which the log(SNR) is drawn\n",
    "smooth_width = 3.0\n",
    "# spatial length scale of the Gaussian Process, unit: voxel\n",
    "inten_kernel = 4.0\n",
    "# intensity length scale of the Gaussian Process\n",
    "# Slightly counter-intuitively, if this parameter is very large,\n",
    "# say, much larger than the range of intensities of the voxels,\n",
    "# then the smoothness has much small dependency on the intensity.\n",
    "\n",
    "\n",
    "inten = np.random.rand(n_V) * 20.0\n",
    "# For simplicity, we just assume that the intensity\n",
    "# of all voxels are uniform distributed between 0 and 20\n",
    "# parameters of Gaussian process to generate pseuso SNR\n",
    "# For curious user, you can also try the following commond\n",
    "# to see what an example snr map might look like if the intensity\n",
    "# grows linearly in one spatial direction\n",
    "\n",
    "# inten = coords_flat[:,0] * 2\n",
    "\n",
    "\n",
    "inten_tile = np.tile(inten, [n_V, 1])\n",
    "inten_diff2 = (inten_tile - inten_tile.T)**2\n",
    "\n",
    "K = np.exp(-dist2 / smooth_width**2 / 2.0 \n",
    "           - inten_diff2 / inten_kernel**2 / 2.0) * tau**2 \\\n",
    "    + np.eye(n_V) * tau**2 * 0.001\n",
    "# A tiny amount is added to the diagonal of\n",
    "# the GP covariance matrix to make sure it can be inverted\n",
    "L = np.linalg.cholesky(K)\n",
    "snr = np.abs(np.dot(L, np.random.randn(n_V))) * snr_level\n",
    "sqrt_v = noise_level * snr\n",
    "betas_simulated = np.dot(L_full, np.random.randn(n_C, n_V)) * sqrt_v\n",
    "signal = np.dot(design.design_task, betas_simulated)\n",
    "\n",
    "\n",
    "Y = signal + noise + inten\n",
    "# The data to be fed to the program.\n",
    "\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
    "plt.pcolor(np.reshape(snr, [ROI_edge, ROI_edge*2]))\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('pseudo-SNR in a rectangular \"ROI\"')\n",
    "plt.show()\n",
    "\n",
    "idx = np.argmin(np.abs(snr - np.median(snr)))\n",
    "# choose a voxel of medium level SNR.\n",
    "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n",
    "                 facecolor='w', edgecolor='k')\n",
    "noise_plot, = plt.plot(noise[:,idx],'g')\n",
    "signal_plot, = plt.plot(signal[:,idx],'b')\n",
    "plt.legend([noise_plot, signal_plot], ['noise', 'signal'])\n",
    "plt.title('simulated data in an example voxel'\n",
    "          ' with pseudo-SNR of {}'.format(snr[idx]))\n",
    "plt.xlabel('time')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n",
    "                 facecolor='w', edgecolor='k')\n",
    "data_plot, = plt.plot(Y[:,idx],'r')\n",
    "plt.legend([data_plot], ['observed data of the voxel'])\n",
    "plt.xlabel('time')\n",
    "plt.show()\n",
    "\n",
    "idx = np.argmin(np.abs(snr - np.max(snr)))\n",
    "# display the voxel of the highest level SNR.\n",
    "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n",
    "                 facecolor='w', edgecolor='k')\n",
    "noise_plot, = plt.plot(noise[:,idx],'g')\n",
    "signal_plot, = plt.plot(signal[:,idx],'b')\n",
    "plt.legend([noise_plot, signal_plot], ['noise', 'signal'])\n",
    "plt.title('simulated data in the voxel with the highest'\n",
    "          ' pseudo-SNR of {}'.format(snr[idx]))\n",
    "plt.xlabel('time')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(12, 4), dpi=150,\n",
    "                 facecolor='w', edgecolor='k')\n",
    "data_plot, = plt.plot(Y[:,idx],'r')\n",
    "plt.legend([data_plot], ['observed data of the voxel'])\n",
    "plt.xlabel('time')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The reason that the pseudo-SNRs in the example voxels are not too small, while the signal looks much smaller is because we happen to have low amplitudes in our design matrix. The true SNR depends on both the amplitudes in design matrix and the pseudo-SNR. Therefore, be aware that pseudo-SNR does not directly reflects how much signal the data have, but rather a map indicating the relative strength of signal in differerent voxels.\n",
    "#### When you have multiple runs, the noise won't be correlated between runs. Therefore, you should tell BRSA when is the onset of each scan. \n",
    "#### Note that the data (variable Y above) you feed to BRSA is the concatenation of data from all runs along the time dimension, as a 2-D matrix of time x space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scan_onsets = np.int32(np.linspace(0, design.n_TR,num=n_run + 1)[: -1])\n",
    "print('scan onsets: {}'.format(scan_onsets))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit Bayesian RSA to our simulated data\n",
    "\n",
    "The nuisance regressors in typical fMRI analysis (such as head motion signal) are replaced by principal components estimated from residuals after subtracting task-related response. n_nureg tells the model how many principal components to keep from the residual as nuisance regressors, in order to account for spatial correlation in noise. \n",
    "If you prefer not using this approach based on principal components of residuals, you can set auto_nuisance=False, and optionally provide your own nuisance regressors as nuisance argument to BRSA.fit(). In practice, we find that the result is much better with auto_nuisance=True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "brsa = BRSA(GP_space=True, GP_inten=True)\n",
    "# Initiate an instance, telling it\n",
    "# that we want to impose Gaussian Process prior\n",
    "# over both space and intensity.\n",
    "\n",
    "brsa.fit(X=Y, design=design.design_task,\n",
    "         coords=coords_flat, inten=inten, scan_onsets=scan_onsets)\n",
    "\n",
    "# The data to fit should be given to the argument X.\n",
    "# Design matrix goes to design. And so on.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can have a look at the estimated similarity in matrix brsa.C_. \n",
    "#### We can also compare the ideal covariance above with the one recovered, brsa.U_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
    "plt.pcolor(brsa.C_, vmin=-0.1, vmax=1)\n",
    "plt.xlim([0, n_C])\n",
    "plt.ylim([0, n_C])\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Estimated correlation structure\\n shared between voxels\\n'\n",
    "         'This constitutes the output of Bayesian RSA\\n')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
    "plt.pcolor(brsa.U_)\n",
    "plt.xlim([0, 16])\n",
    "plt.ylim([0, 16])\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Estimated covariance structure\\n shared between voxels\\n')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In contrast, we can have a look of the similarity matrix based on Pearson correlation between point estimates of betas of different conditions.\n",
    "#### This is what vanila RSA might give"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "regressor = np.insert(design.design_task,\n",
    "                      0, 1, axis=1)\n",
    "betas_point = np.linalg.lstsq(regressor, Y)[0]\n",
    "point_corr = np.corrcoef(betas_point[1:, :])\n",
    "point_cov = np.cov(betas_point[1:, :])\n",
    "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
    "plt.pcolor(point_corr, vmin=-0.1, vmax=1)\n",
    "plt.xlim([0, 16])\n",
    "plt.ylim([0, 16])\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Correlation structure estimated\\n'\n",
    "         'based on point estimates of betas\\n')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(4, 4), dpi=100)\n",
    "plt.pcolor(point_cov)\n",
    "plt.xlim([0, 16])\n",
    "plt.ylim([0, 16])\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Covariance structure of\\n'\n",
    "         'point estimates of betas\\n')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can make a comparison between the estimated SNR map and the true SNR map (normalized)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(num=None, figsize=(5, 5), dpi=100)\n",
    "plt.pcolor(np.reshape(brsa.nSNR_, [ROI_edge, ROI_edge*2]))\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "ax.set_title('estimated pseudo-SNR')\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure(num=None, figsize=(5, 5), dpi=100)\n",
    "plt.pcolor(np.reshape(snr / np.exp(np.mean(np.log(snr))),\n",
    "                      [ROI_edge, ROI_edge*2]))\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "ax.set_title('true normalized pseudo-SNR')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RMS_BRSA = np.mean((brsa.C_ - ideal_corr)**2)**0.5\n",
    "RMS_RSA = np.mean((point_corr - ideal_corr)**2)**0.5\n",
    "print('RMS error of Bayesian RSA: {}'.format(RMS_BRSA))\n",
    "print('RMS error of standard RSA: {}'.format(RMS_RSA))\n",
    "print('Recovered spatial smoothness length scale: '\n",
    "      '{}, vs. true value: {}'.format(brsa.lGPspace_, smooth_width))\n",
    "print('Recovered intensity smoothness length scale: '\n",
    "      '{}, vs. true value: {}'.format(brsa.lGPinten_, inten_kernel))\n",
    "print('Recovered standard deviation of GP prior: '\n",
    "      '{}, vs. true value: {}'.format(brsa.bGP_, tau))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Empirically, the smoothness turns to be over-estimated when signal is weak."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can also look at how other parameters are recovered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "plt.scatter(rho1, brsa.rho_)\n",
    "plt.xlabel('true AR(1) coefficients')\n",
    "plt.ylabel('recovered AR(1) coefficients')\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(np.log(snr) - np.mean(np.log(snr)),\n",
    "            np.log(brsa.nSNR_))\n",
    "plt.xlabel('true normalized log SNR')\n",
    "plt.ylabel('recovered log pseudo-SNR')\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Even though the variation reduced in estimated pseudo-SNR (due to overestimation of smoothness of the GP prior under low SNR situation), betas recovered by the model has higher correlation with true betas than doing simple regression, shown below. Obiously there is shrinkage of the estimated betas, as a result of variance-bias tradeoff. But we think such shrinkage does preserve the patterns of betas, and therefore the result is suitable to be further used for decoding purpose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(betas_simulated, brsa.beta_)\n",
    "plt.xlabel('true betas (response amplitudes)')\n",
    "plt.ylabel('recovered betas by Bayesian RSA')\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.scatter(betas_simulated, betas_point[1:, :])\n",
    "plt.xlabel('true betas (response amplitudes)')\n",
    "plt.ylabel('recovered betas by simple regression')\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The singular decomposition of noise, and the comparison between the first two principal component of noise and the patterns of the first two nuisance regressors, returned by the model. \n",
    "The principal components may not look exactly the same. The first principal components both capture the baseline image intensities (although they may sometimes appear counter-phase)\n",
    "\n",
    "Apparently one can imagine that the choice of the number of principal components used as nuisance regressors can influence the result. If you just choose 1 or 2, perhaps only the global drift would be captured. But including too many nuisance regressors would slow the fitting speed and might have risk of overfitting. The users might consider starting in the range of 5-20. We do not have automatic cross-validation built in. But you can use the score() function to do cross-validation and select the appropriate number. The idea here is similar to that in GLMdenoise (http://kendrickkay.net/GLMdenoise/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u, s, v = np.linalg.svd(noise + inten)\n",
    "plt.plot(s)\n",
    "plt.xlabel('principal component')\n",
    "plt.ylabel('singular value of unnormalized noise')\n",
    "plt.show()\n",
    "\n",
    "plt.pcolor(np.reshape(v[0,:], [ROI_edge, ROI_edge*2]))\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Weights of the first principal component in unnormalized noise')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.pcolor(np.reshape(brsa.beta0_[0,:], [ROI_edge, ROI_edge*2]))\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Weights of the DC component in noise')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.pcolor(np.reshape(inten, [ROI_edge, ROI_edge*2]))\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('The baseline intensity of the ROI')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.pcolor(np.reshape(v[1,:], [ROI_edge, ROI_edge*2]))\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Weights of the second principal component in unnormalized noise')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "plt.pcolor(np.reshape(brsa.beta0_[1,:], [ROI_edge, ROI_edge*2]))\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('Weights of the first recovered noise pattern\\n not related to DC component in noise')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \"Decoding\" from new data\n",
    "### Now we generate a new data set, assuming signal is the same but noise is regenerated. We want to use the transform() function of brsa to estimate the \"design matrix\" in this new dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_new = np.zeros([n_T, n_V])\n",
    "noise_new[0, :] = np.dot(L_noise, np.random.randn(n_V))\\\n",
    "    / np.sqrt(1 - rho1**2)\n",
    "for i_t in range(1, n_T):\n",
    "    noise_new[i_t, :] = noise_new[i_t - 1, :] * rho1 \\\n",
    "        + np.dot(L_noise,np.random.randn(n_V))\n",
    "Y_new = signal + noise_new + inten\n",
    "\n",
    "\n",
    "ts, ts0 = brsa.transform(Y_new,scan_onsets=scan_onsets)\n",
    "# ts, ts0 = brsa.transform(Y_new,scan_onsets=scan_onsets)\n",
    "recovered_plot, = plt.plot(ts[:200, 8], 'b')\n",
    "design_plot, = plt.plot(design.design_task[:200, 8], 'g')\n",
    "plt.legend([design_plot, recovered_plot],\n",
    "           ['design matrix for one condition', 'recovered time course for the condition'])\n",
    "plt.show()\n",
    "# We did not plot the whole time series for the purpose of seeing closely how much the two\n",
    "# time series overlap\n",
    "\n",
    "c = np.corrcoef(design.design_task.T, ts.T)\n",
    "\n",
    "\n",
    "# plt.pcolor(c[0:n_C, n_C:],vmin=-0.5,vmax=1)\n",
    "plt.pcolor(c[0:16, 16:],vmin=-0.5,vmax=1)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('correlation between true design matrix \\nand the recovered task-related activity')\n",
    "plt.colorbar()\n",
    "plt.xlabel('recovered task-related activity')\n",
    "plt.ylabel('true design matrix')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# plt.pcolor(c[n_C:, n_C:],vmin=-0.5,vmax=1)\n",
    "plt.pcolor(c[16:, 16:],vmin=-0.5,vmax=1)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('correlation within the recovered task-related activity')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model selection by cross-validataion:\n",
    "You can compare different models by cross-validating the parameters of one model learnt from some training data\n",
    "on some testing data. BRSA provides a score() function, which provides you a pair of cross-validated log likelihood\n",
    "for testing data. The first value is the cross-validated log likelihood of the model you have specified. The second\n",
    "value is a null model which assumes everything else the same except that there is no task-related activity.\n",
    "### Notice that comparing the score of your model of interest against its corresponding null model is not the single way to compare models. You might also want to compare against a model using the same set of design matrix, but a different rank (especially rank 1, which means all task conditions have the same response pattern, only differing in the magnitude).\n",
    "In general, in the context of BRSA, a model means the timing of each event and the way these events are grouped, together with other trivial parameters such as the rank of the covariance matrix and the number of nuisance regressors. All these parameters can influence model performance.\n",
    "In future, we will provide interface to test the performance of a model with predefined similarity matrix or covariance matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[score, score_null] = brsa.score(X=Y_new, design=design.design_task, scan_onsets=scan_onsets)\n",
    "print(\"Score of full model based on the correct esign matrix, assuming {} nuisance\"\n",
    "     \" components in the noise: {}\".format(brsa.n_nureg_, score))\n",
    "print(\"Score of a null model with the same assumption except that there is no task-related response: {}\".format(\n",
    "        score_null))\n",
    "plt.bar([0,1],[score, score_null], width=0.5)\n",
    "plt.ylim(np.min([score, score_null])-100, np.max([score, score_null])+100)\n",
    "plt.xticks([0,1],['Model','Null model'])\n",
    "plt.ylabel('cross-validated log likelihood')\n",
    "plt.title('cross validation on new data')\n",
    "plt.show()\n",
    "\n",
    "[score_noise, score_noise_null] = brsa.score(X=noise_new+inten, design=design.design_task, scan_onsets=scan_onsets)\n",
    "print(\"Score of full model for noise, based on the correct design matrix, assuming {} nuisance\"\n",
    "     \" components in the noise: {}\".format(brsa.n_nureg_, score_noise))\n",
    "print(\"Score of a null model for noise: {}\".format(\n",
    "        score_noise_null))\n",
    "plt.bar([0,1],[score_noise, score_noise_null], width=0.5)\n",
    "plt.ylim(np.min([score_noise, score_noise_null])-100, np.max([score_noise, score_noise_null])+100)\n",
    "plt.xticks([0,1],['Model','Null model'])\n",
    "plt.ylabel('cross-validated log likelihood')\n",
    "plt.title('cross validation on noise')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen above, the model with the correct design matrix explains new data with signals generated from the true model better than the null model, but explains pure noise worse than the null model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# We can also try the version which marginalize SNR and rho for each voxel. \n",
    "### This version is intended for analyzing data of a group of participants and estimating their shared similarity matrix. But it also allows analyzing single participant. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "gbrsa = GBRSA(nureg_method='PCA', auto_nuisance=True, logS_range=1,\n",
    "              anneal_speed=20, n_iter=50)\n",
    "# Initiate an instance, telling it\n",
    "# that we want to impose Gaussian Process prior\n",
    "# over both space and intensity.\n",
    "\n",
    "gbrsa.fit(X=Y, design=design.design_task,scan_onsets=scan_onsets)\n",
    "\n",
    "# The data to fit should be given to the argument X.\n",
    "# Design matrix goes to design. And so on.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolor(np.reshape(gbrsa.nSNR_, (ROI_edge, ROI_edge*2)))\n",
    "plt.colorbar()\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('SNR map estimated by marginalized BRSA')\n",
    "plt.show()\n",
    "\n",
    "plt.pcolor(np.reshape(snr, (ROI_edge, ROI_edge*2)))\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.colorbar()\n",
    "plt.title('true SNR map')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(snr, gbrsa.nSNR_)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.xlabel('simulated pseudo-SNR')\n",
    "plt.ylabel('estimated pseudo-SNR')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(np.log(snr), np.log(gbrsa.nSNR_))\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.xlabel('simulated log(pseudo-SNR)')\n",
    "plt.ylabel('estimated log(pseudo-SNR)')\n",
    "plt.show()\n",
    "\n",
    "plt.pcolor(gbrsa.U_)\n",
    "plt.colorbar()\n",
    "plt.title('covariance matrix estimated by marginalized BRSA')\n",
    "plt.show()\n",
    "plt.pcolor(ideal_cov)\n",
    "plt.colorbar()\n",
    "plt.title('true covariance matrix')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(betas_simulated, gbrsa.beta_)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.xlabel('simulated betas')\n",
    "plt.ylabel('betas estimated by marginalized BRSA')\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(rho1, gbrsa.rho_)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.xlabel('simulated AR(1) coefficients')\n",
    "plt.ylabel('AR(1) coefficients estimated by marginalized BRSA')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also do \"decoding\" and cross-validating using the marginalized version in GBRSA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# \"Decoding\"\n",
    "ts, ts0 = gbrsa.transform([Y_new],scan_onsets=[scan_onsets])\n",
    "\n",
    "recovered_plot, = plt.plot(ts[0][:200, 8], 'b')\n",
    "design_plot, = plt.plot(design.design_task[:200, 8], 'g')\n",
    "plt.legend([design_plot, recovered_plot],\n",
    "           ['design matrix for one condition', 'recovered time course for the condition'])\n",
    "plt.show()\n",
    "# We did not plot the whole time series for the purpose of seeing closely how much the two\n",
    "# time series overlap\n",
    "\n",
    "c = np.corrcoef(design.design_task.T, ts[0].T)\n",
    "\n",
    "\n",
    "plt.pcolor(c[0:n_C, n_C:],vmin=-0.5,vmax=1)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('correlation between true design matrix \\nand the recovered task-related activity')\n",
    "plt.colorbar()\n",
    "plt.xlabel('recovered task-related activity')\n",
    "plt.ylabel('true design matrix')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.pcolor(c[n_C:, n_C:],vmin=-0.5,vmax=1)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect(1)\n",
    "plt.title('correlation within the recovered task-related activity')\n",
    "plt.colorbar()\n",
    "plt.show()\n",
    "\n",
    "# cross-validataion\n",
    "[score, score_null] = gbrsa.score(X=[Y_new], design=[design.design_task], scan_onsets=[scan_onsets])\n",
    "print(\"Score of full model based on the correct design matrix, assuming {} nuisance\"\n",
    "     \" components in the noise: {}\".format(gbrsa.n_nureg_, score))\n",
    "print(\"Score of a null model with the same assumption except that there is no task-related response: {}\".format(\n",
    "        score_null))\n",
    "plt.bar([0,1],[score[0], score_null[0]], width=0.5)\n",
    "plt.ylim(np.min([score[0], score_null[0]])-100, np.max([score[0], score_null[0]])+100)\n",
    "plt.xticks([0,1],['Model','Null model'])\n",
    "plt.ylabel('cross-validated log likelihood')\n",
    "plt.title('cross validation on new data')\n",
    "plt.show()\n",
    "\n",
    "[score_noise, score_noise_null] = gbrsa.score(X=[noise_new+inten], design=[design.design_task],\n",
    "                                              scan_onsets=[scan_onsets])\n",
    "print(\"Score of full model for noise, based on the correct design matrix, assuming {} nuisance\"\n",
    "     \" components in the noise: {}\".format(gbrsa.n_nureg_, score_noise))\n",
    "print(\"Score of a null model for noise: {}\".format(\n",
    "        score_noise_null))\n",
    "plt.bar([0,1],[score_noise[0], score_noise_null[0]], width=0.5)\n",
    "plt.ylim(np.min([score_noise[0], score_noise_null[0]])-100, np.max([score_noise[0], score_noise_null[0]])+100)\n",
    "plt.xticks([0,1],['Model','Null model'])\n",
    "plt.ylabel('cross-validated log likelihood')\n",
    "plt.title('cross validation on noise')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}