{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functional connectivity and resting state\n",
    "\n",
    "Functional connectivity and resting-state data can be studied in many different way. `Nilearn` provides tools to construct \"connectomes\" that capture functional interactions between regions or to extract regions and networks, via resting-state networks or parcellations. For a much more detailed guide, go to [Nilearn's Connectivity section](http://nilearn.github.io/connectivity/index.html), here we want to show you just a few basics.\n",
    "\n",
    "Speaking of which, we will be covering the following sections:\n",
    "\n",
    "1. Extracting times series to build a functional connectome\n",
    "1. Single subject maps of seed-to-voxel correlation\n",
    "1. Single subject maps of seed-to-seed correlation\n",
    "1. Group analysis of resting-state fMRI with ICA (CanICA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup\n",
    "\n",
    "Before we start with anything, let's set up the important plotting functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import plotting\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from nilearn import image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, let's specify which subjects we want to use for this notebook. So, who do we have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nib-ls /data/adhd/*/*.nii.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each subject we also have a regressor file, containing important signal confounds such motion parameters, compcor components and mean signal in CSF, GM, WM, Overal. Let's take a look at one of those regressor files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.read_table('/data/adhd/0010042/0010042_regressors.csv').head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So let's create two lists, containing the path to the resting-state images and confounds of all subjects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cp -R /data/adhd/ /home/neuro/workshop/notebooks/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Which subjects to consider\n",
    "sub_idx = ['0010042', '0010064', '0010128', '0021019', '0027018',\n",
    "           '0027034', '0027037', '1517058', '1562298', '2497695',\n",
    "           '2950754', '3007585', '3520880', '3994098', '4134561',\n",
    "           '6115230', '8409791', '8697774', '9744150', '9750701']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to resting state files\n",
    "rest_files = ['/data/adhd/%s/%s_rest_tshift_RPI_voreg_mni.nii.gz' % (sub, sub) for sub in sub_idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to counfound files\n",
    "confound_files = ['/data/adhd/%s/%s_regressors.csv' % (sub, sub) for sub in sub_idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perfect, now we're good to go!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Extracting times series to build a functional connectome\n",
    "\n",
    "So let's start with something simple: Extracting activation time series from regions defined by a parcellation atlas.\n",
    "\n",
    "## Brain parcellation\n",
    "\n",
    "As a first step, let's define the regions we want to extract the signal from:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Location of HarvardOxford parcellation atlas\n",
    "atlas_file = '/usr/share/fsl/data/atlases/HarvardOxford/HarvardOxford-cort-maxprob-thr25-2mm.nii.gz'\n",
    "\n",
    "# Visualize parcellation atlas\n",
    "plotting.plot_roi(atlas_file, draw_cross=False, annotate=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load labels for each atlas region\n",
    "with open('data/HarvardOxford_labels.txt', 'r') as f:\n",
    "    labels = f.read().split('\\n')\n",
    "labels[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting signals on a parcellation\n",
    "\n",
    "To extract signal on the parcellation, the easiest option is to use `NiftiLabelsMasker`. As any \"maskers\" in nilearn, it is a processing object that is created by specifying all the important parameters, but not the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiLabelsMasker\n",
    "masker = NiftiLabelsMasker(labels_img=atlas_file, standardize=True, verbose=1,\n",
    "                           memory=\"nilearn_cache\", memory_level=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Nifti data can then be turned to time-series by calling the `NiftiLabelsMasker` `fit_transform` method, that takes either filenames or NiftiImage objects.\n",
    "\n",
    "Let's do this now for the first subject:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_series = masker.fit_transform(rest_files[0], confounds=confound_files[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute and display the correlation matrix\n",
    "\n",
    "Now we're read to compute the functional connectome with `ConnectivityMeasure`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.connectome import ConnectivityMeasure\n",
    "correlation_measure = ConnectivityMeasure(kind='correlation')\n",
    "correlation_matrix = correlation_measure.fit_transform([time_series])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we can visualize this correlation matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mask the main diagonal for visualization:\n",
    "np.fill_diagonal(correlation_matrix, 0)\n",
    "\n",
    "# Plot correlation matrix - note: matrix is ordered for block-like representation\n",
    "plotting.plot_matrix(correlation_matrix, figure=(10, 8), labels=labels,\n",
    "                     vmax=0.8, vmin=-0.8, reorder=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Same thing without confounds, to stress the importance of confounds\n",
    "\n",
    "Let's do the same thing as before, but this time without using the confounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the signal from the regions\n",
    "time_series_bad = masker.fit_transform(rest_files[0]) # Note that we haven't specify confounds here\n",
    "\n",
    "# Compute the correlation matrix\n",
    "correlation_matrix_bad = correlation_measure.fit_transform([time_series_bad])[0]\n",
    "\n",
    "# Mask the main diagonal for visualization\n",
    "np.fill_diagonal(correlation_matrix_bad, 0)\n",
    "\n",
    "# Plot the correlation matrix\n",
    "plotting.plot_matrix(correlation_matrix_bad, figure=(10, 8), labels=labels,\n",
    "                     vmax=0.8, vmin=-0.8, title='No confounds', reorder=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, without any confounds all regions are connected to each other! One reference that discusses the importance of confounds is [Varoquaux & Craddock 2013](http://www.sciencedirect.com/science/article/pii/S1053811913003340)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probabilistic atlas\n",
    "\n",
    "Above we used a parcellation atlas. Now, with nilearn, you can do the same thing also with a probabilistic atlas. Let's use for example the [MSDL atlas](https://team.inria.fr/parietal/18-2/spatial_patterns/spatial-patterns-in-resting-state/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to MSDL atlas\n",
    "msdl_atlas = 'data/msdl_atlas/MSDL_rois/msdl_rois.nii'\n",
    "\n",
    "# Extract only default mode network nodes\n",
    "dmn_nodes = image.index_img(msdl_atlas, [3, 4, 5, 6])\n",
    "\n",
    "# Plot MSDL probability atlas\n",
    "plotting.plot_prob_atlas(dmn_nodes, cut_coords=(0, -60, 29), draw_cross=False,\n",
    "                         annotate=False, title=\"DMN nodes in MSDL atlas\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only difference to before is that we now need to use the `NiftiMapsMasker` function to create the masker that extracts the time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiMapsMasker\n",
    "masker = NiftiMapsMasker(maps_img=msdl_atlas, standardize=True, verbose=1,\n",
    "                         memory=\"nilearn_cache\", memory_level=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the signal from the regions\n",
    "time_series = masker.fit_transform(rest_files[0], confounds=confound_files[0])\n",
    "\n",
    "# Compute the correlation matrix\n",
    "correlation_matrix= correlation_measure.fit_transform([time_series])[0]\n",
    "\n",
    "# Mask the main diagonal for visualization\n",
    "np.fill_diagonal(correlation_matrix, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we plot the new correlation matrix, we also need to load the labels of the MSDL atlas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CSV containing label and coordinate of MSDL atlas\n",
    "msdl_info = 'data/msdl_atlas/MSDL_rois/msdl_rois_labels.csv'\n",
    "\n",
    "# Load the name and coordinates of the labels\n",
    "content = pd.read_csv(msdl_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read out label names\n",
    "labels = content['name'].tolist()\n",
    "labels[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read out label coordinats\n",
    "coords = np.array(content[['x', 'y', 'z']])\n",
    "coords[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perfect! Now as before, we can plot the correlation matrix as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the correlation matrix\n",
    "plotting.plot_matrix(correlation_matrix, figure=(10, 8), labels=labels,\n",
    "                     vmax=0.8, vmin=-0.8, reorder=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display corresponding graph on glass brain\n",
    "\n",
    "A square matrix, such as a correlation matrix, can also be seen as a \"graph\": a set of \"nodes\", connected by \"edges\". When these nodes are brain regions, and the edges capture interactions between them, this graph is a \"functional connectome\".\n",
    "\n",
    "As the MSDL atlas comes with (x, y, z) MNI coordinates for the different regions, we can visualize the matrix as a graph of interaction in a brain. To avoid having too dense a graph, we represent only the 20% edges with the highest values. For another atlas this information can be computed for each region with the `nilearn.plotting.find_xyz_cut_coords` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_connectome(correlation_matrix, coords, edge_threshold=\"80%\",\n",
    "                         colorbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the correlation matrix gives a very \"full\" graph: every node is connected to every other one. This is because it also captures indirect connections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From version `0.5.0` on, `nilearn` also provides an interactive plot for connectoms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.view_connectome(correlation_matrix, coords, threshold=\"80%\", cmap='bwr',\n",
    "                         symmetric_cmap=False, linewidth=6.0, marker_size=3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Single subject maps of seed-to-voxel correlation\n",
    "\n",
    "Above we computed the correlation between different regions. But what if we want to compute a seed-to-voxel correlation map for a single subject? The procedure is very similar to the one from before.\n",
    "\n",
    "## Time series extraction\n",
    "\n",
    "First, we need to extract the time series from the seed region. For this example, let's specify a sphere of radius 8 (in mm) located in the Posterior Cingulate Cortex. This sphere is considered to be part of the Default Mode Network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sphere radius in mm\n",
    "sphere_radius = 8\n",
    "\n",
    "# Sphere center in MNI-coordinate\n",
    "sphere_coords = [(0, -52, 18)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, we will use We use the `NiftiSpheresMasker` function to extract the time series within a given sphere. Before signal extraction, we can also directly detrend, standardize, and bandpass filter the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiSpheresMasker\n",
    "seed_masker = NiftiSpheresMasker(sphere_coords, radius=sphere_radius, detrend=True,\n",
    "                                 standardize=True, low_pass=0.1, high_pass=0.01,\n",
    "                                 t_r=2.0, verbose=1, memory=\"nilearn_cache\", memory_level=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're read to extract the mean time series within the seed region, while also regressing out the confounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_time_series = seed_masker.fit_transform(rest_files[0], confounds=confound_files[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to do a similar procedure for each voxel in the brain as well. For this, we can use the `NiftiMasker`, which the same arguments as before, plus additionally smoothing the signal with a smoothing kernel of 6mm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiMasker\n",
    "brain_masker = NiftiMasker(smoothing_fwhm=6, detrend=True, standardize=True,\n",
    "                           low_pass=0.1, high_pass=0.01, t_r=2., verbose=1,\n",
    "                           memory=\"nilearn_cache\", memory_level=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can extract the time series for every voxel while regressing out the confounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brain_time_series = brain_masker.fit_transform(rest_files[0], confounds=confound_files[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing the seed-based correlation analysis\n",
    "\n",
    "Now that we have two arrays (mean signal in seed region, signal for each voxel), we can correlate the two to each other. This can be done with the dot product between the two matrices.\n",
    "\n",
    "**Note**, that the signals have been variance-standardized during extraction. To have them standardized to\n",
    "norm unit, we further have to divide the result by the length of the time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_based_correlations = np.dot(brain_time_series.T, seed_time_series)\n",
    "seed_based_correlations /= seed_time_series.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the seed-based correlation map\n",
    "\n",
    "Finally, we can tranform the correlation array back to a Nifti image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_based_correlation_img = brain_masker.inverse_transform(seed_based_correlations.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this we can of course plot again to better investigate the correlation outcome:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display = plotting.plot_stat_map(seed_based_correlation_img, threshold=0.333,\n",
    "                                 cut_coords=sphere_coords[0])\n",
    "display.add_markers(marker_coords=sphere_coords, marker_color='black',\n",
    "                    marker_size=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The map above depicts the temporal correlation of a **seed region** with the **rest of the brain**. For a similar example but on the cortical surface, see [this example](http://nilearn.github.io/auto_examples/01_plotting/plot_surf_stat_map.html#seed-based-connectivity-on-the-surface)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Single subject maps of seed-to-seed correlation\n",
    "\n",
    "The next question is of course, how can compute the correlation between different seed regions?  It's actually very easy, even simpler than above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time series extraction\n",
    "\n",
    "First, we need to extract the time series from the seed regions. So as before, we need to define a sphere radius and centers of the seed regions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sphere radius in mm\n",
    "sphere_radius = 8\n",
    "\n",
    "# Sphere center in MNI-coordinate\n",
    "sphere_center = [(  0, -52, 18),\n",
    "                 (-46, -68, 32),\n",
    "                 ( 46, -68, 32),\n",
    "                 (  1,  50, -5)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can extract the time series from those spheres:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create masker object to extract average signal within spheres\n",
    "from nilearn.input_data import NiftiSpheresMasker\n",
    "masker = NiftiSpheresMasker(sphere_center, radius=sphere_radius, detrend=True,\n",
    "                            standardize=True, low_pass=0.1, high_pass=0.01,\n",
    "                            t_r=2.0, verbose=1, memory=\"nilearn_cache\", memory_level=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract average signal in spheres with masker object\n",
    "time_series = masker.fit_transform(rest_files[0], confounds=confound_files[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display mean signal per sphere\n",
    "\n",
    "If we want, we can even plot the average signal per sphere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 4))\n",
    "plt.plot(time_series)\n",
    "plt.xlabel('Scan number')\n",
    "plt.ylabel('Normalized signal')\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute partial correlation matrix\n",
    "\n",
    "Now that we have the average signal per sphere we can compute the partial correlation matrix, using the `ConnectivityMeasure` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.connectome import ConnectivityMeasure\n",
    "connectivity_measure = ConnectivityMeasure(kind='partial correlation')\n",
    "partial_correlation_matrix = connectivity_measure.fit_transform([time_series])[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the partical correlation matrix\n",
    "fig, ax = plt.subplots()\n",
    "plt.imshow(partial_correlation_matrix, cmap='Spectral')\n",
    "for (j,i),label in np.ndenumerate(partial_correlation_matrix):\n",
    "    ax.text(i, j, round(label, 2), ha='center', va='center', color='w')\n",
    "    ax.text(i, j, round(label, 2), ha='center', va='center', color='w')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display connectome\n",
    "\n",
    "Now that we have the correlation matrix, we can also plot it again on the glass brain with `plot_connectome`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.plotting import plot_connectome\n",
    "plot_connectome(partial_correlation_matrix, sphere_center,\n",
    "                display_mode='ortho', colorbar=True,  node_size=150,\n",
    "                title=\"Default Mode Network Connectivity\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And again with `nilearn`'s interactive connectome viewer function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.view_connectome(partial_correlation_matrix, sphere_center, cmap='bwr',\n",
    "                         symmetric_cmap=True, linewidth=6.0, marker_size=10.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Group analysis of resting-state fMRI with ICA (CanICA)\n",
    "\n",
    "This section demonstrates the use of multi-subject Independent Component Analysis (ICA) of resting-state fMRI data to extract brain networks in a data-driven way. Here we use the `CanICA` approach, that implements a multivariate random effects model across subjects. Afterward this, we will also show a newer technique, based on dictionary learning.\n",
    "\n",
    "## Multi-subject ICA: CanICA\n",
    "\n",
    "`CanICA` is a ready-to-use object that can be applied to multi-subject Nifti data, for instance, presented as filenames, and will perform a multi-subject ICA decomposition following the `CanICA` model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.decomposition import CanICA\n",
    "\n",
    "# Number of components to extract\n",
    "n_components = 20\n",
    "\n",
    "# Creating the CanICA object\n",
    "canica = CanICA(n_components=n_components,\n",
    "                smoothing_fwhm=6.,\n",
    "                threshold=3.,\n",
    "                random_state=0,\n",
    "                n_jobs=-1,\n",
    "                verbose=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with every object in nilearn, we give its parameters at construction, and then fit it on the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "canica.fit(rest_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once `CanICA`has finished we can retrieve the independent components directly in brain space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "components_img = canica.components_img_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing CanICA components\n",
    "\n",
    "To visualize the components we can plot the outline of all components in one figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot all ICA components together\n",
    "plotting.plot_prob_atlas(components_img, draw_cross=False, linewidths=None,\n",
    "                         cut_coords=[0, 0, 0], title='All ICA components');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can of course also plot the ICA components separately using the `plot_stat_map` and `iter_img`. Let's plot the first few components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract first few components\n",
    "first_few_comp = components_img.slicer[..., :4]\n",
    "\n",
    "# Plot first few components\n",
    "fig = plt.figure(figsize=(16, 3))\n",
    "for i, cur_img in enumerate(image.iter_img(first_few_comp)):\n",
    "    ax = fig.add_subplot(1, 4, i + 1)\n",
    "    plotting.plot_stat_map(cur_img, display_mode=\"z\", title=\"IC %d\" % i, \n",
    "                  cut_coords=1, colorbar=True, axes=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Beyond ICA : Dictionary learning\n",
    "\n",
    "Recent work has shown that **dictionary learning** based techniques outperform ICA in term of stability and constitutes a better first step in a statistical analysis pipeline. Dictionary learning in neuroimaging seeks to extract a few representative temporal elements along with their ***sparse spatial loadings***, which constitute good extracted maps.\n",
    "\n",
    "So let's do the same thing again as above, but this time with `DictLearning`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import dictionary learning algorithm\n",
    "from nilearn.decomposition import DictLearning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize DictLearning object\n",
    "dict_learn = DictLearning(n_components=n_components,\n",
    "                          smoothing_fwhm=6.,\n",
    "                          random_state=0,\n",
    "                          n_jobs=-1,\n",
    "                          verbose=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're ready and can apply the dictionar learning object on the functional data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Fit to the data\n",
    "dict_learn.fit(rest_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we can now retrieve the independent components directly in brain space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "components_img = dict_learn.components_img_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing DictLearning components\n",
    "\n",
    "To visualize the components we can plot the outline of all components in one figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all ICA components together\n",
    "plotting.plot_prob_atlas(components_img, draw_cross=False, linewidths=None,\n",
    "                         cut_coords=[0, 0, 0], title='Dictionary Learning maps');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract first few components\n",
    "first_few_comp = components_img.slicer[..., :4]\n",
    "\n",
    "# Plot first few components\n",
    "fig = plt.figure(figsize=(16, 3))\n",
    "for i, cur_img in enumerate(image.iter_img(first_few_comp)):\n",
    "    ax = fig.add_subplot(1, 4, i + 1)\n",
    "    plotting.plot_stat_map(cur_img, display_mode=\"z\", title=\"IC %d\" % i, \n",
    "                  cut_coords=1, colorbar=True, axes=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare ICA to DictLearning\n",
    "\n",
    "Now let's compare the two approaches by looking at some components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_stat_map(canica.components_img_.slicer[..., 3], display_mode='ortho',\n",
    "                       cut_coords=[45, -35, 50], colorbar=True, draw_cross=False,\n",
    "                       title='CanICA component - Motor Cortex')\n",
    "plotting.plot_stat_map(dict_learn.components_img_.slicer[..., 19], display_mode='ortho',\n",
    "                       cut_coords=[45, -35, 50], colorbar=True, draw_cross=False,\n",
    "                       title='DictLearning component - Motor Cortex')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_stat_map(canica.components_img_.slicer[..., 16],  display_mode='ortho',\n",
    "                       cut_coords=[50, -15, 12], colorbar=True, draw_cross=False,\n",
    "                       title='CanICA component - Auditory Cortex')\n",
    "plotting.plot_stat_map(dict_learn.components_img_.slicer[..., 16],  display_mode='ortho',\n",
    "                       cut_coords=[50, -15, 12], colorbar=True, draw_cross=False,\n",
    "                       title='DictLearning component - Auditory Cortex')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plotting.plot_stat_map(canica.components_img_.slicer[..., 0], display_mode='ortho',\n",
    "                       cut_coords=[0, -75, 4], colorbar=True, draw_cross=False,\n",
    "                       title='CanICA component - Visual Cortex')\n",
    "plotting.plot_stat_map(dict_learn.components_img_.slicer[..., 3], display_mode='ortho',\n",
    "                       cut_coords=[0, -75, 4], colorbar=True, draw_cross=False,\n",
    "                       title='DictLearning component - Visual Cortex')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the CanICA components looks much more noise, while the DictLearning components look sparser and more blobby. This becomes even more striking when we look at the corresponding glass brain plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_glass_brain(canica.components_img_.slicer[..., 16], black_bg=True,\n",
    "                          plot_abs=False, symmetric_cbar=False,\n",
    "                          title='CanICA component - Auditory Cortex', colorbar=True)\n",
    "plotting.plot_glass_brain(dict_learn.components_img_.slicer[..., 16], black_bg=True,\n",
    "                          title='DictLearning component - Auditory Cortex', colorbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maps obtained with dictionary leaning are often easier to exploit as they are less noisy than ICA maps, with blobs usually better defined. Typically, smoothing can be lower than when doing ICA. While dictionary learning computation time is comparable to CanICA, obtained atlases have been shown to outperform ICA in a variety of classification tasks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract functional connectome based on dictionary learning components\n",
    "\n",
    "Similar to the very first section of this notebook, we can now take the components from the dictionary learning and compute a correlation matrix between the regions defined by the components.\n",
    "\n",
    "We will be using nilearn's `RegionExtractor` to extract brain connected regions from the dictionary maps. We will be using the automatic thresholding strategy `ratio_n_voxels`. We use this thresholding strategy to first get foreground information present in the maps and then followed by robust region extraction on foreground information using Random Walker algorithm selected as `extractor='local_regions'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.regions import RegionExtractor\n",
    "extractor = RegionExtractor(dict_learn.components_img_, threshold=0.5,\n",
    "                            thresholding_strategy='ratio_n_voxels',\n",
    "                            extractor='local_regions', verbose=1,\n",
    "                            standardize=True, min_region_size=1350)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we control foreground extraction using parameter `threshold=0.5`, which represents the expected proportion of voxels included in the regions (i.e. with a non-zero value in one of the maps). If you need to keep more proportion of voxels then threshold should be tweaked according to the maps data.\n",
    "\n",
    "The parameter `min_region_size=1350 mm^3` is to keep the minimum number of extracted regions. We control the small spurious regions size by thresholding in voxel units to adapt well to the resolution of the image. Please see the documentation of [`nilearn.regions.connected_regions`](http://nilearn.github.io/modules/generated/nilearn.regions.connected_regions.html#nilearn.regions.connected_regions) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extractor.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So how many regions did we extract?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Total number of regions extracted\n",
    "n_regions_extracted = extractor.regions_img_.shape[-1]\n",
    "n_regions_extracted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, to get the average functional connectome over all subjects we need to compute the correlation matrix for each subject individually and than average those matrices into one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from nilearn.connectome import ConnectivityMeasure\n",
    "\n",
    "# Initializing ConnectivityMeasure object with kind='correlation'\n",
    "connectome_measure = ConnectivityMeasure(kind='correlation')\n",
    "\n",
    "# Iterate over the subjects and compute correlation matrix for each\n",
    "correlations = []\n",
    "for filename, confound in zip(rest_files, confound_files):\n",
    "    \n",
    "    timeseries_each_subject = extractor.transform(filename, confounds=confound)\n",
    "    \n",
    "    correlation = connectome_measure.fit_transform([timeseries_each_subject])\n",
    "    \n",
    "    correlations.append(correlation)\n",
    "\n",
    "# Get array in good numpy structure\n",
    "correlations = np.squeeze(correlations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that this is all computed, we can take the average correlation matrix and plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computing the mean correlation matrix\n",
    "mean_correlations = np.mean(correlations, axis=0)\n",
    "\n",
    "# Plot the average correlation matrix\n",
    "title = 'Correlation between %d regions' % n_regions_extracted\n",
    "plotting.plot_matrix(mean_correlations, vmax=1, vmin=-1, colorbar=True,\n",
    "                     labels=['IC %0d' % i for i in range(n_regions_extracted)],\n",
    "                     title=title, reorder=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And as a last step, let's plot the average function connectome, based on the dictionary learning components also on the glass brain.\n",
    "\n",
    "For this to work, we first need to find the center of the regions. Luckily nilearn provides a nice function, called `find_xyz_cut_coords` that does exactly that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the center of the regions with find_xyz_cut_coords\n",
    "coords_connectome = [plotting.find_xyz_cut_coords(img)\n",
    "                     for img in image.iter_img(extractor.regions_img_)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the functional connectome on the glass brain\n",
    "plotting.plot_connectome(mean_correlations, coords_connectome,\n",
    "                         edge_threshold='95%', title=title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.view_connectome(mean_correlations, coords_connectome,\n",
    "                          cmap='bwr', symmetric_cmap=True, linewidth=6.0, marker_size=10.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}