{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nTutorial 2: Customizing and aligning gradients\n=================================================\nIn this tutorial you\u2019ll learn about the methods available within the\nGradientMaps class. The flexible usage of this class allows for the\ncustomization of gradient computation with different kernels and dimensionality\nreductions, as well as aligning gradients from different datasets. This\ntutorial will only show you how to apply these techniques.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Customizing gradient computation\n+++++++++++++++++++++++++++++++++\nAs before, we\u2019ll start by loading the sample data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from brainspace.datasets import load_group_fc, load_parcellation, load_conte69\n\n# First load mean connectivity matrix and Schaefer parcellation\nconn_matrix = load_group_fc('schaefer', scale=400)\nlabeling = load_parcellation('schaefer', scale=400, join=True)\n\nmask = labeling != 0\n\n# and load the conte69 hemisphere surfaces\nsurf_lh, surf_rh = load_conte69()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The GradientMaps object allows for many different kernels and dimensionality\nreduction techniques. Let\u2019s have a look at three different kernels.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nfrom brainspace.gradient import GradientMaps\nfrom brainspace.plotting import plot_hemispheres\nfrom brainspace.utils.parcellation import map_to_labels\n\nkernels = ['pearson', 'spearman', 'normalized_angle']\n\ngradients_kernel = [None] * len(kernels)\nfor i, k in enumerate(kernels):\n    gm = GradientMaps(kernel=k, approach='dm', random_state=0)\n    gm.fit(conn_matrix)\n\n    gradients_kernel[i] = map_to_labels(gm.gradients_[:, i], labeling, mask=mask,\n                                        fill=np.nan)\n\n\nlabel_text = ['Pearson', 'Spearman', 'Normalized\\nAngle']\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_kernel, size=(1200, 600),\n                 cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.45)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It seems the gradients provided by these kernels are quite similar although\ntheir scaling is quite different. Do note that the gradients are in arbitrary\nunits, so the smaller/larger axes across kernels do not imply anything.\nSimilar to using different kernels, we can also use different dimensionality\nreduction techniques.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# PCA, Laplacian eigenmaps and diffusion mapping\nembeddings = ['pca', 'le', 'dm']\n\ngradients_embedding = [None] * len(embeddings)\nfor i, emb in enumerate(embeddings):\n    gm = GradientMaps(kernel='normalized_angle', approach=emb, random_state=0)\n    gm.fit(conn_matrix)\n\n    gradients_embedding[i] = map_to_labels(gm.gradients_[:, 0], labeling, mask=mask,\n                                           fill=np.nan)\n\n\n# sphinx_gallery_thumbnail_number = 2\nlabel_text = ['PCA', 'LE', 'DM']\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_embedding, size=(1200, 600),\n                 cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.45)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Gradient alignment\n+++++++++++++++++++\n\nA more principled way of increasing comparability across gradients are\nalignment techniques. BrainSpace provides two alignment techniques:\nProcrustes analysis, and joint alignment. For this example we will load\nfunctional connectivity data of a second subject group and align it with the\nfirst group.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "conn_matrix2 = load_group_fc('schaefer', scale=400, group='holdout')\ngp = GradientMaps(kernel='normalized_angle', alignment='procrustes')\ngj = GradientMaps(kernel='normalized_angle', alignment='joint')\n\ngp.fit([conn_matrix, conn_matrix2])\ngj.fit([conn_matrix, conn_matrix2])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Here, `gp` contains the Procrustes aligned data and `gj` contains the joint\naligned data. Let\u2019s plot them, but in separate figures to keep things\norganized.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# First gradient from original and holdout data, without alignment\ngradients_unaligned = [None] * 2\nfor i in range(2):\n    gradients_unaligned[i] = map_to_labels(gp.gradients_[i][:, 0], labeling,\n                                           mask=mask, fill=np.nan)\n\nlabel_text = ['Unaligned\\nGroup 1', 'Unaligned\\nGroup 2']\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_unaligned, size=(1200, 400),\n                 cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.5)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# With procrustes alignment\ngradients_procrustes = [None] * 2\nfor i in range(2):\n    gradients_procrustes[i] = map_to_labels(gp.aligned_[i][:, 0], labeling, mask=mask,\n                                            fill=np.nan)\n\nlabel_text = ['Procrustes\\nGroup 1', 'Procrustes\\nGroup 2']\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_procrustes, size=(1200, 400),\n                 cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.5)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# With joint alignment\ngradients_joint = [None] * 2\nfor i in range(2):\n    gradients_joint[i] = map_to_labels(gj.aligned_[i][:, 0], labeling, mask=mask,\n                                       fill=np.nan)\n\nlabel_text = ['Joint\\nGroup 1', 'Joint\\nGroup 2']\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_joint, size=(1200, 400),\n                 cmap='viridis_r', color_bar=True, label_text=label_text, zoom=1.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Although in this example, we don't see any big differences, if the input data\nwas less similar, alignments may also resolve changes in the order of the\ngradients. However, you should always inspect the output of an alignment;\nif the input data are sufficiently dissimilar then the alignment may produce\nodd results.\n\n\nIn some instances, you may want to align gradients to an out-of-sample\ngradient, for example when aligning individuals to a hold-out group gradient.\nWhen performing a Procrustes alignemnt, a 'reference' can be specified.\nThe first alignment iteration will then be to the reference. For purposes of\nthis example, we will use the gradient of the hold-out group as the\nreference.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gref = GradientMaps(kernel='normalized_angle', approach='le')\ngref.fit(conn_matrix2)\n\ngalign = GradientMaps(kernel='normalized_angle', approach='le', alignment='procrustes')\ngalign.fit(conn_matrix, reference=gref.gradients_)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The gradients in `galign.aligned_` are now aligned to the reference\ngradients.\n\nGradient fusion\n+++++++++++++++++++\nWe can also fuse data across multiple modalities and build mutli-modal\ngradients. In this case we only look at one set of output gradients,\nrather than one per modality.\n\nFirst, let's load the example data of microstructural profile covariance\n`(Paquola et al., 2019) <https://journals.plos.org/plosbiology/article?\nid=10.1371/journal.pbio.3000284>`_ and functional connectivity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from brainspace.datasets import load_group_mpc\n\n# First load mean connectivity matrix and parcellation\nfc = load_group_fc('vosdewael', scale=200)\nmpc = load_group_mpc('vosdewael', scale=200)\n\nlabeling = load_parcellation('vosdewael', scale=200, join=True)\nmask = labeling != 0\n\nseeds = [None] * 2\nseeds[0] = map_to_labels(fc[0], labeling, mask=mask, fill=np.nan)\nseeds[1] = map_to_labels(mpc[0], labeling, mask=mask, fill=np.nan)\n\n# visualise the features from a seed region (seed 0)\nplot_hemispheres(surf_lh, surf_rh, array_name=seeds, label_text=['FC', 'MPC'],\n                 size=(1200, 400), color_bar=True, cmap='viridis', zoom=1.45)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In order to fuse the matrices, we simply pass the matrices to the fusion\ncommand which will rescale and horizontally concatenate the matrices.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Negative numbers are not allowed in fusion.\nfc[fc < 0] = 0\n\n\ndef fusion(*args):\n    from scipy.stats import rankdata\n    from sklearn.preprocessing import minmax_scale\n\n    max_rk = [None] * len(args)\n    masks = [None] * len(args)\n    for j, a in enumerate(args):\n        m = masks[j] = a != 0\n        a[m] = rankdata(a[m])\n        max_rk[j] = a[m].max()\n\n    max_rk = min(max_rk)\n    for j, a in enumerate(args):\n        m = masks[j]\n        a[m] = minmax_scale(a[m], feature_range=(1, max_rk))\n\n    return np.hstack(args)\n\n\n# fuse the matrices\nfused_matrix = fusion(fc, mpc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We then use this output in the fit function. This will convert the long\nhorizontal array into a square affinity matrix, and then perform embedding.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gm = GradientMaps(n_components=2, kernel='normalized_angle')\ngm.fit(fused_matrix)\n\n\ngradients_fused = [None] * 2\nfor i in range(2):\n    gradients_fused[i] = map_to_labels(gm.gradients_[:, i], labeling, mask=mask,\n                                       fill=np.nan)\n\nplot_hemispheres(surf_lh, surf_rh, array_name=gradients_fused,\n                 label_text=['Gradient 1', 'Gradient 2'], size=(1200, 400),\n                 color_bar=True, cmap='viridis', zoom=1.45)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>The mpc matrix presented here matches the subject cohort of `(Paquola et\n  al., 2019) <https://journals.plos.org/plosbiology/article?id=10.1371/\n  journal.pbio.3000284>`_. Other matrices in this package match the subject\n  groups used by `(Vos de Wael et al., 2018) <https://www.pnas.org/content/\n  115/40/10154.short>`_. We make direct comparisons in our tutorial for\n  didactic purposes only.</p></div>\n\nThat concludes the second tutorial. In the third tutorial we will consider\nnull hypothesis testing of comparisons between gradients and other markers.\n\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}