MiCM2022_Dim_Reduction_notebook_full.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f7041c7b",
   "metadata": {},
   "source": [
    "# Introduction to Dimensionality Reduction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30e36797",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import mne\n",
    "from mne.preprocessing import ICA\n",
    "from scipy.stats import zscore\n",
    "import pandas as pd\n",
    "from itertools import combinations\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from seaborn import lmplot\n",
    "import pickle\n",
    "import random\n",
    "import scanpy as sc\n",
    "import anndata\n",
    "from IPython.display import display, clear_output"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09d1ed79",
   "metadata": {},
   "source": [
    "#### A walk on the beach"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b9911dc",
   "metadata": {},
   "source": [
    "<center><img src=\"beach.jpeg\" width=\"600\"></center>\n",
    "<div style=\"text-align: center\"> source: https://www.publicdomainpictures.net/en/free-download.php?image=shadows-on-the-beach&id=177457 </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf3c5e1",
   "metadata": {},
   "source": [
    "#### Looking for repetition: a recipe book"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dcbed629",
   "metadata": {},
   "source": [
    "<center><img src=\"recipes2.png\" width=\"800\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dfd6c7b",
   "metadata": {},
   "source": [
    "## Principal Components Analysis (PCA)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f73b030e",
   "metadata": {},
   "source": [
    "With PCA, data from a high-dimensional space (e.g. 2D, a plane) can be projected onto a lower-dimensional space (e.g. 1D, a line).\n",
    "\n",
    "<center><img src=\"pca_proj.jpeg\" width=\"500\"></center>\n",
    "<div style=\"text-align: center\"> source: https://programmathically.com/principal-components-analysis-explained-for-dummies/ </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad360ea4",
   "metadata": {},
   "source": [
    "Here is an example of 3D data projected onto a 2D plane."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42c95f87",
   "metadata": {},
   "source": [
    "<center><img src=\"pca.png\" width=\"800\"></center>\n",
    "<div style=\"text-align: center\"> source: https://www.publicdomainpictures.net/en/free-download.php?image=shadows-on-the-beach&id=177457 </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a54091a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Import iris dataset\n",
    "from sklearn import datasets\n",
    "iris = datasets.load_iris()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36ceb21f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Store data in pandas DataFrame\n",
    "iris_df = pd.DataFrame(data= np.c_[iris['data'], iris['target']],\n",
    "                     columns= iris['feature_names'] + ['target'])\n",
    "\n",
    "## Rename columns\n",
    "iris_df['target'] = iris_df['target'].map({0:iris.target_names[0], 1:iris.target_names[1], 2:iris.target_names[2]})\n",
    "iris_df.rename(columns = {'target':'species'}, inplace=True)\n",
    "\n",
    "## Display data and list species' names\n",
    "display(iris_df)\n",
    "print(iris_df['species'].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab637e55",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Extract numerical values in arrays\n",
    "x = iris_df.iloc[:,:-1].values\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83e9f68a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Normalize data: zero mean & unit variance\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "scaler = StandardScaler()\n",
    "x = scaler.fit_transform(x)\n",
    "x_trunc = x[:,:-1]  ## first three features, for visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "960c01f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compare histograms of features before and after applying the standard scaler\n",
    "\n",
    "print('Before scaling')\n",
    "iris_df.hist(sharex=True); plt.show()\n",
    "\n",
    "print('After scaling')\n",
    "pd.DataFrame(x, columns=iris_df.columns[:-1].str.strip(' (cm)')).hist(sharex=True); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8210735",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Visualize truncated data containing first three features\n",
    "\n",
    "fig = plt.figure(1, figsize=(12, 8))\n",
    "ax = fig.add_subplot(111, projection=\"3d\", elev=-150, azim=110)\n",
    "\n",
    "ax.scatter(x_trunc[:,0], x_trunc[:,1], x_trunc[:,2]);\n",
    "\n",
    "ax.set_xlabel(iris_df.columns[0]); ax.set_ylabel(iris_df.columns[1]); ax.set_zlabel(iris_df.columns[2])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ad03174",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compare 2D projections of full vs truncated data\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=[6.,6.])\n",
    "ax.set_xlabel('PCA 1'); ax.set_ylabel('PCA 2')\n",
    "ax.set_xticks([]); ax.set_yticks([])\n",
    "\n",
    "## 2D PCA projection of truncated data\n",
    "pca = PCA(n_components=2)\n",
    "x_2d = pca.fit_transform(x_trunc)\n",
    "ax.scatter(x_2d[:,0], x_2d[:,1], label='truncated data\\n(3 features)')\n",
    "\n",
    "## 2D PCA projection of full data\n",
    "pca = PCA(n_components=2)\n",
    "x_2d = pca.fit_transform(x)\n",
    "ax.scatter(x_2d[:,0], x_2d[:,1], label='full data\\n(4 features)')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f12c5ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Show \"manual\" 2D PCA projection alongside PCA weights for different features\n",
    "\n",
    "fig, ax = plt.subplots(1,2,figsize=[12,6])\n",
    "ax[0].set_xlabel('PCA 1')  # proportion of variance explained by PCA component 1\n",
    "ax[0].set_ylabel('PCA 2')  # proportion of variance explained by PCA component 2\n",
    "ax[0].set_xticks([]); ax[0].set_yticks([])\n",
    "ax[0].set_title('PCA projection ($\\hat{X} = XU^T$)')\n",
    "\n",
    "## Matrix multiplication to obtain projection\n",
    "x_2d_manual = x @ pca.components_[:2,:].T\n",
    "ax[0].scatter(x_2d_manual[:,0], x_2d_manual[:,1]);\n",
    "\n",
    "xlim = ax[0].get_xlim()\n",
    "ylim = ax[0].get_ylim()\n",
    "\n",
    "## Barplot of components' weights\n",
    "bp = pd.DataFrame(pca.components_, columns=iris_df.columns[:-1].str.strip(' (cm)'), index=['PCA 1', 'PCA 2']).plot.bar(ax=ax[1], rot=0);\n",
    "\n",
    "bp.set_ylabel('weight');\n",
    "for p in ax[1].patches: ax[1].annotate(str(round(p.get_height(), 2)), (p.get_x()-0.005, p.get_height()*1.025))\n",
    "ax[1].set_title('PCA weights ($U$)')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e01d06fe",
   "metadata": {},
   "source": [
    "$$\n",
    "PCA.1 = 0.52 * sepal.length - 0.27 * sepal.width + 0.58 * petal.length + 0.56 * petal.width\\\\\n",
    "PCA.2 = 0.38 * sepal.length - 0.92 * sepal.width + 0.02 * petal.length + 0.07 * petal.width\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c49d3bcd",
   "metadata": {},
   "source": [
    "Performing a PCA projection amounts to applying the equations above to create PCA variables, which we can plot. In the code below, we will loop through each data point to compute PCA projections using the data features and PCA weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "579e0233",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create figure\n",
    "\n",
    "fig = plt.figure(figsize=[8,8])\n",
    "ax = fig.add_subplot(1, 1, 1) \n",
    "ax.set_xlim(xlim)\n",
    "ax.set_ylim(ylim)\n",
    "ax.set_xlabel('PCA 1')\n",
    "ax.set_ylabel('PCA 2')\n",
    "\n",
    "\n",
    "## Loop through data points\n",
    "for i in range(len(x)):\n",
    "    \n",
    "    xi = x[i,:].squeeze()\n",
    "    \n",
    "    c1 = pca.components_[0,:]\n",
    "    c2 = pca.components_[1,:]\n",
    "    \n",
    "    xi_pc1 = (xi * c1).sum()\n",
    "    xi_pc2 = (xi * c2).sum()\n",
    "    ax.scatter(xi_pc1, xi_pc2, c='blue')\n",
    "    \n",
    "    dpi = 'Data point ' + str(i)\n",
    "    pc1_1 = '({:.2f})*({:.2f})'.format(c1[0] , xi[0])\n",
    "    pc1_2 = '({:.2f})*({:.2f})'.format(c1[1] , xi[1])\n",
    "    pc1_3 = '({:.2f})*({:.2f})'.format(c1[2] , xi[2])\n",
    "    pc1_4 = '({:.2f})*({:.2f})'.format(c1[3] , xi[3])\n",
    "\n",
    "    pc2_1 = '({:.2f})*({:.2f})'.format(c2[0] , xi[0])\n",
    "    pc2_2 = '({:.2f})*({:.2f})'.format(c2[1] , xi[1])\n",
    "    pc2_3 = '({:.2f})*({:.2f})'.format(c2[2] , xi[2])\n",
    "    pc2_4 = '({:.2f})*({:.2f})'.format(c2[3] , xi[3])\n",
    "\n",
    "    pc_info = str(\n",
    "    '''\n",
    "    {:<16}: {:>15}    {:>15}    {:>15}    {:>15}\n",
    "    {:<16}: {:>15.2f}    {:>15.2f}    {:>15.2f}    {:>15.2f}\n",
    "    {:<16}: {:>16} + {:>16} + {:>16} + {:>16} = {:>5.2f}\n",
    "    {:<16}: {:>16} + {:>16} + {:>16} + {:>16} = {:>5.2f}\n",
    "    '''.format(\n",
    "    'features', 'sepal length', 'sepal width', 'petal length', 'petal width',\n",
    "    dpi, *xi,\n",
    "    'PCA 1' , pc1_1 , pc1_2 , pc1_3 , pc1_4, xi_pc1,\n",
    "    'PCA 2' , pc2_1 , pc2_2 , pc2_3 , pc2_4, xi_pc2))\n",
    "    \n",
    "    display(fig)\n",
    "    print(pc_info)\n",
    "  \n",
    "    clear_output(wait = True)\n",
    "    plt.pause(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26e5daaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "## \"Manual\" 2D projection - project dimensions individually for one data point\n",
    "\n",
    "xi = x[0,:]\n",
    "print('Feature values of xi:', xi, '\\n')\n",
    "print('Weights for PCA 1:', pca.components_[0,:])\n",
    "print('Weights for PCA 2:', pca.components_[1,:], '\\n')\n",
    "\n",
    "xi_1 = (xi * pca.components_[0,:]).sum()  # project onto 1st component\n",
    "xi_2 = (xi * pca.components_[1,:]).sum()  # project onto 2nd component\n",
    "\n",
    "print('Projection of xi onto PCA 1 and PCA 2 - manual computation:', np.float32([xi_1, xi_2]))\n",
    "print('Projection of xi onto PCA 1 and PCA 2 - automatic computation:', np.float32(x_2d[0,:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ca86f1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compare full \"manual\" 2D projection to automatic 2D projection\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=[6.,6.])\n",
    "ax.set_xlabel('PCA 1'); ax.set_ylabel('PCA 2')\n",
    "ax.set_xticks([]); ax.set_yticks([])\n",
    "\n",
    "## 2D PCA projection of full data\n",
    "pca = PCA(n_components=2)\n",
    "x_2d = pca.fit_transform(x)\n",
    "ax.scatter(x_2d[:,0], x_2d[:,1], label='automatic projection')\n",
    "\n",
    "## 2D PCA projection of full data - matrix multiplication\n",
    "x_2d_manual = x @ pca.components_.T\n",
    "ax.scatter(x_2d_manual[:,0], x_2d_manual[:,1], label='manual projection');\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d10cd0b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Color data points by species\n",
    "\n",
    "iris_df_2d = pd.DataFrame(data = np.c_[x_2d, iris_df['species']], columns=['PCA 1', 'PCA 2', 'species'])\n",
    "lmplot(x='PCA 1', y='PCA 2', data=iris_df_2d, hue='species', fit_reg=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e35edcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute 4 PCA components rather than just 2 components\n",
    "pca_4d = PCA(n_components=4)\n",
    "x_4d = pca_4d.fit_transform(x)\n",
    "\n",
    "## Plot all 6 possible 2D PCA projections\n",
    "p = np.arange(4)\n",
    "pcombs = list(combinations(p,2))\n",
    "\n",
    "plt_cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "colors = {'setosa':plt_cols[0], 'versicolor':plt_cols[1], 'virginica':plt_cols[2]}\n",
    "\n",
    "_, ax = plt.subplots(2,3, figsize=[15.,10.])\n",
    "for i, pcomb in enumerate(pcombs):\n",
    "    axi = np.unravel_index(i, ax.shape)\n",
    "    ax[axi].scatter(x_4d[:,pcomb[0]], x_4d[:,pcomb[1]], c=iris_df['species'].map(colors))\n",
    "    \n",
    "    ax[axi].set_xlabel('PCA ' + str(pcomb[0]+1)); ax[axi].set_ylabel('PCA ' + str(pcomb[1]+1))\n",
    "    ax[axi].set_xticks([]); ax[axi].set_yticks([])\n",
    "    \n",
    "\n",
    "## Plot singular values and related quantities\n",
    "fig, ax = plt.subplots(1,4,figsize=[15.,3.])\n",
    "\n",
    "sv_1 = pca_4d.singular_values_\n",
    "sv_2 = pca_4d.singular_values_**2\n",
    "sv_3 = pca_4d.singular_values_**2 / (pca_4d.singular_values_**2).sum()\n",
    "sv_4 = pca_4d.explained_variance_ratio_\n",
    "\n",
    "pd.DataFrame(sv_1, index=['PCA 1', 'PCA 2', 'PCA 3', 'PCA 4']).plot.bar(ax=ax[0], rot=0, legend=False, color='cyan', title='singular values');\n",
    "for p in ax[0].patches: ax[0].annotate(str(round(p.get_height(), 2)), (p.get_x()+0.01, p.get_height()*0.9))\n",
    "pd.DataFrame(sv_2, index=['PCA 1', 'PCA 2', 'PCA 3', 'PCA 4']).plot.bar(ax=ax[1], rot=0, legend=False, color='cyan', title='squared singular values');\n",
    "for p in ax[1].patches: ax[1].annotate(str(round(p.get_height(), 2)), (p.get_x()+0.01, p.get_height()*0.9))\n",
    "pd.DataFrame(sv_3, index=['PCA 1', 'PCA 2', 'PCA 3', 'PCA 4']).plot.bar(ax=ax[2], rot=0, legend=False, color='cyan', title='ratio of squared sing. values');\n",
    "for p in ax[2].patches: ax[2].annotate(str(round(p.get_height(), 2)), (p.get_x()+0.01, p.get_height()*0.9))\n",
    "pd.DataFrame(sv_4, index=['PCA 1', 'PCA 2', 'PCA 3', 'PCA 4']).plot.bar(ax=ax[3], rot=0, legend=False, color='cyan', title='variance explained');\n",
    "for p in ax[3].patches: ax[3].annotate(str(round(p.get_height(), 2)), (p.get_x()+0.01, p.get_height()*0.9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9b79418",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot cumulative sum of variances explained by the PCA components\n",
    "\n",
    "_, ax = plt.subplots(1,1)\n",
    "pd.DataFrame(sv_4.cumsum(), index=['PCA 1', 'PCA 2', 'PCA 3', 'PCA 4']).plot.bar(ax=ax, rot=0, legend=False, color='cyan', title='variance explained');\n",
    "for p in ax.patches: ax.annotate(str(round(p.get_height(), 2)), (p.get_x()+0.1, p.get_height()*0.9))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3face644",
   "metadata": {},
   "source": [
    "#### The biplot\n",
    "\n",
    "We've plotted the projections of individual samples onto the 2D PCA represenation, effectively achieving dimensionality reduction for the purposes of e.g. data visualization. We've also *separately* plotted the weights of the PCA components to visualize how different features are weighted when projecting data onto the 2D representation.\n",
    "\n",
    "Alternatively, one can visualize *both* the samples' projection onto the 2D representation as well as the features' weights simultaneously via a **biplot**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a6b659f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def biplot(score, coeff , y, variables, extend=1):\n",
    "    '''\n",
    "    \n",
    "    ***** Author: Serafeim Loukas, serafeim.loukas@epfl.ch *****\n",
    "    https://ai.plainenglish.io/how-to-implement-pca-with-python-and-scikit-learn-22f3de4e5983\n",
    "    \n",
    "    Inputs:\n",
    "       score: the projected data\n",
    "       coeff: the eigenvectors (PCs)\n",
    "       y: the class labels\n",
    "   '''\n",
    "    \n",
    "    xs = score[:,0] # projection on PC1\n",
    "    ys = score[:,1] # projection on PC2\n",
    "    n = coeff.shape[0] # number of variables\n",
    "    plt.figure(figsize=(8,6), dpi=100)\n",
    "    classes = np.unique(y)\n",
    "    for s,l in enumerate(classes):\n",
    "        plt.scatter(xs[y==l],ys[y==l], label=l) # color based on group\n",
    "    plt.legend()\n",
    "    for i in range(n):\n",
    "        #plot as arrows the variable scores (each variable has a score for PC1 and one for PC2)\n",
    "        plt.arrow(0, 0, coeff[i,0] * extend, coeff[i,1] * extend, color = 'k', alpha = 0.9,linestyle = '-',linewidth = 1.5, overhang=0.2)\n",
    "        plt.text(coeff[i,0]* extend, coeff[i,1] * extend, variables[i], color = 'k', ha = 'center', va = 'center',fontsize=10)\n",
    "    plt.xlabel(\"PCA {}\".format(1), size=14)\n",
    "    plt.ylabel(\"PCA {}\".format(2), size=14)\n",
    "    limx= int(xs.max()) + 1\n",
    "    limy= int(ys.max()) + 1\n",
    "    plt.xlim([-limx,limx])\n",
    "    plt.ylim([-limy,limy])\n",
    "    plt.grid()\n",
    "    plt.tick_params(axis='both', which='both', labelsize=14)\n",
    "    \n",
    "    \n",
    "def bisubplot(score, coeff , y, variables, pcomb, ax, extend=1):\n",
    "    '''\n",
    "    \n",
    "    ***** Author: Serafeim Loukas, serafeim.loukas@epfl.ch *****\n",
    "    https://ai.plainenglish.io/how-to-implement-pca-with-python-and-scikit-learn-22f3de4e5983\n",
    "    \n",
    "    Inputs:\n",
    "       score: the projected data\n",
    "       coeff: the eigenvectors (PCs)\n",
    "       y: the class labels\n",
    "   '''\n",
    "    \n",
    "    xs = score[:,0] # projection on PC1\n",
    "    ys = score[:,1] # projection on PC2\n",
    "    n = coeff.shape[0] # number of variables\n",
    "    classes = np.unique(y)\n",
    "    for s,l in enumerate(classes):\n",
    "        ax.scatter(xs[y==l],ys[y==l], label=l) # color based on group\n",
    "    ax.legend()\n",
    "    for i in range(n):\n",
    "        #plot as arrows the variable scores (each variable has a score for PC1 and one for PC2)\n",
    "        ax.arrow(0, 0, coeff[i,0] * extend, coeff[i,1] * extend, color = 'k', alpha = 0.9,linestyle = '-',linewidth = 1.5, overhang=0.2)\n",
    "        ax.text(coeff[i,0]* extend, coeff[i,1] * extend, variables[i], color = 'k', ha = 'center', va = 'center',fontsize=10)\n",
    "    ax.set_xlabel(\"PCA {}\".format(pcomb[0]+1), size=14)\n",
    "    ax.set_ylabel(\"PCA {}\".format(pcomb[1]+1), size=14)\n",
    "    limx= int(xs.max()) + 1\n",
    "    limy= int(ys.max()) + 1\n",
    "    ax.set_xlim([-limx,limx])\n",
    "    ax.set_ylim([-limy,limy])\n",
    "    ax.grid()\n",
    "    ax.tick_params(axis='both', which='both', labelsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4e196f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Barplot of components' weights\n",
    "\n",
    "fig, ax = plt.subplots(1,3,figsize=[14.,4.5])\n",
    "\n",
    "ax[0].scatter(x_2d[:,0], x_2d[:,1], c=iris_df['species'].map(colors))\n",
    "ax[0].set_xlabel('PCA 1'); ax[0].set_ylabel('PCA 2')\n",
    "\n",
    "## PCA weights as plotted previously\n",
    "bp = pd.DataFrame(pca.components_, columns=iris_df.columns[:-1].str.strip(' (cm)'), index=['PCA 1', 'PCA 2']).plot.bar(ax=ax[1], rot=0);\n",
    "\n",
    "bp.set_ylabel('weights');\n",
    "for p in ax[1].patches:\n",
    "    ax[1].annotate(str(round(p.get_height(), 2)), (p.get_x()-0.005, p.get_height()*1.025))\n",
    "    \n",
    "## PCA weights, plotted by feature rather than PCA component\n",
    "bp = pd.DataFrame(pca.components_.T, columns=['PCA 1','PCA 2'], index=iris_df.columns[:-1].str.strip(' (cm)')).plot.bar(ax=ax[2], rot=0, color=['pink','cyan']);\n",
    "\n",
    "bp.set_ylabel('weights');\n",
    "for p in ax[2].patches:\n",
    "    ax[2].annotate(str(round(p.get_height(), 2)), (p.get_x()-0.005, p.get_height()*1.025))\n",
    "    \n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a9af819",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create biplot\n",
    "biplot(x_2d, pca.components_.T, iris_df['species'], iris_df.columns[:-1].str.strip(' (cm)'), extend=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db5f4aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Histograms of original features, grouped by species\n",
    "fig, ax = plt.subplots(1,4,figsize=[18,3])\n",
    "for c,col in enumerate(iris_df.columns[:-1]):\n",
    "    h = iris_df.pivot(columns='species', values=col).plot.hist(alpha=0.8, title=col, ax=ax[c])\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3064c800",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot all 6 possible biplots arising from the 4 PCA components\n",
    "\n",
    "_, ax = plt.subplots(3,2, figsize=[15.,24.])\n",
    "for i, pcomb in enumerate(pcombs):\n",
    "    axi = np.unravel_index(i, ax.shape)\n",
    "    bisubplot(x_4d[:,pcomb], pca_4d.components_[pcomb,:].T, iris_df['species'], iris_df.columns[:-1].str.strip(' (cm)'), pcomb, ax[axi], extend=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1322543",
   "metadata": {},
   "source": [
    "#### Kernel PCA\n",
    "\n",
    "PCA provides *linear* components which can be unsuitable for data with non-linear structures. Non-linear functions called *kernels* can be used to represent data such that PCA can provide suitable linear components of that non-linear representation. This method is termed **Kernel PCA**.\n",
    "\n",
    "<img src=\"kpca.png\" width=\"900\">\n",
    "<div style=\"text-align: center\"> source: https://www.cs.mcgill.ca/~dprecup/courses/ML/Lectures/ml-lecture13.pdf </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80d8a29d",
   "metadata": {},
   "source": [
    "#### Probabilistic PCA\n",
    "PCA does not attempt to model noise that accompany its components. In **Probabilistic PCA**, a noise variance term $\\sigma^2$ accounts for random variations within its components, allowing projections against its components with varying spreads. Probabilistic PCA can thus model both the linear structure and noise components of its components. Due to its probabilistic nature, Probabilistic PCA can also be used as a generative model. Note that Probabilistic PCA can be viewed as a special case of **Factor Analysis**.\n",
    "\n",
    "<img src=\"ppca.png\" width=\"1000\">\n",
    "<div style=\"text-align: center\"> source: https://www.youtube.com/watch?v=lJ0cXPoEozg </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e264a72",
   "metadata": {},
   "source": [
    "### Independent Component Analysis (ICA)\n",
    "\n",
    "PCA provides *orthogonal* components / coordinates, which can limit the modelling capacities of PCA. On the other hand, such ICA components are not constrained to be orthogonal which provides ICA with more \"flexibility\".\n",
    "\n",
    "<center><img src=\"pca_vs_ica.png\" width=\"800\"></center>\n",
    "<div style=\"text-align: center\"> source: https://brainvision.com/applications/eeg/ </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c43870c",
   "metadata": {},
   "source": [
    "#### Electroencephalography (EEG) data analysis\n",
    " \n",
    " The EEG is a brain recording modality where electrodes are placed on the scalp/head to measure fluctuations in bioelectical potentials. Hence, every electrode measures a time-varying signal.\n",
    " \n",
    "  <center><img src=\"eeg_image.png\" width=\"400\"></center>\n",
    "  <div style=\"text-align: center\"> source: https://brainvision.com/applications/eeg/ </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba9d5c29",
   "metadata": {},
   "source": [
    "Import EEG data using [MNE Python](https://mne.tools/stable/index.html) analysis package. We will be using open-source data which has been reported in:\n",
    "\n",
    "[Babayan, A., Erbey, M., Kumral, D. et al. *A mind-brain-body dataset of MRI, EEG, cognition, emotion, and peripheral physiology in young and old adults.* Sci Data 6, 180308 (2019).](https://www.nature.com/articles/sdata2018308)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0a29b50",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load EEG data\n",
    "eeg = mne.io.read_raw_eeglab(\"sub-010321_EC_downsamp.set\")\n",
    "eeg.annotations.delete( np.arange( len(eeg.annotations.description) ) )  # remove annotations, not important"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "446de7d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot EEG electrode layout and EEG signals\n",
    "eeg.plot_sensors(show_names=True, sphere=90);\n",
    "eeg.plot(n_channels=12, duration=4);  # recommend plotting using \"notebook\" plotting backend, i.e. %matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be7ea1b6",
   "metadata": {},
   "source": [
    "From the time-series plots, it seems that many time-series are alike. Perhaps there's a way to find patterns and have a more succinct representation of our data. One way to find these patterns is to use **Independent Components Analysis** (ICA) on our EEG data.\n",
    "\n",
    "  <center><img src=\"eeg_ica.jpeg\" width=\"500\"></center>\n",
    "  <div style=\"text-align: center\"> source: https://sccn.ucsd.edu/~jung/Site/EEG_artifact_removal.html </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "898b56f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute ICA on EEG data\n",
    "num_comps = 15  # choose number of patterns, or \"components\"\n",
    "\n",
    "ica = ICA(n_components=num_comps, random_state=97)\n",
    "ica.fit(eeg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d15d2f26",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot ICA components\n",
    "ica.plot_sources(eeg, start=10, stop=16);\n",
    "ica.plot_components();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f7adf35",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute PCA on EEG data\n",
    "num_comps = 15  # choose number of patterns, or \"components\"\n",
    "\n",
    "eeg_ts = eeg.get_data().T\n",
    "pca = PCA(n_components=num_comps)\n",
    "pca_data_ts = pca.fit_transform(eeg_ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5818394",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot PCA components\n",
    "len_plot = 1000     # choose time-series length for plotting purposes\n",
    "\n",
    "num_plots = pca_data_ts.shape[1]\n",
    "plt.figure(figsize=[9.,9.])\n",
    "\n",
    "for i in range(num_plots):\n",
    "    ax = plt.subplot( num_plots , 1, i+1 )\n",
    "    plt.plot(zscore(pca_data_ts)[:len_plot,i])\n",
    "    \n",
    "    ax.set_ylabel(\"PCA {}\".format(i), fontsize=8, rotation=0)\n",
    "    ax.set_xticks([]); ax.set_yticks([])\n",
    "    ax.spines[\"top\"].set_visible(False)\n",
    "    ax.spines[\"right\"].set_visible(False)\n",
    "    ax.spines[\"left\"].set_visible(False)\n",
    "    ax.spines[\"bottom\"].set_visible(False)\n",
    "    \n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1750db33",
   "metadata": {},
   "source": [
    "### Nonlinear dimensionality reduction: t-SNE and UMAP\n",
    "\n",
    "Linear dimensionality reduction methods such as PCA may be incapable of properly displaying data in a low-dimensional space (e.g. 2D). As previously discussed, kernel PCA leverages nonlinear transformations (i.e. kernels) such that the resulting representation is ammenable to standard PCA. However, kernel PCA requires a pre-defined kernel function to obtain a suitable nonlinear representation.\n",
    "\n",
    "We will now explore two nonlinear dimensionality reduction methods which do not require such pre-defined kernels. These are:\n",
    "1. t-distributed Stochastic Neighbor Embedding (**t-SNE**)\n",
    "2. Uniform Manifold Approximation and Projection (**UMAP**)\n",
    "\n",
    "Essentially, t-SNE and UMAP aim to preserve the \"similarities\" between pairs of datapoints. In other words, points which are \"similar\" in high-dimensional space should remain \"similar\" in low-dimensional space, where similarities are quantified kernels computed over distances. The concept of a *graph* is thus relevant for t-SNE and UMAP.\n",
    "\n",
    "<center><img src=\"umap.png\" width=\"800\"></center>\n",
    ".\n",
    "<div style=\"text-align: center\"> source: https://umap-learn.readthedocs.io/en/latest/parametric_umap.html </div>\n",
    "\n",
    "Although UMAP and t-SNE are very similar methods for embedding high-dimensional data into low-dimensional space, they have their differences. These methods use different kernels to quantify distances, symmetrize distances differently, have different number of parameters, and initialize low-dimensional embeddings differently. To the last point, t-SNE initializes low-dimensional embeddings randomly, whereas UMAP does so using spectral embedding. UMAP should thus provide much more stable and reproducible results compared to t-SNE.\n",
    "\n",
    "Another (debated) difference between t-SNE and UMAP is the ability to preserve *global* structure in addition to *local* structure. One way of understanding this distinction is to view the graph structure in low-dimensional space as a set of attractive and repulsive forces between datapoints. The low-dimensional embedding is obtained once the forces equilibrate, similar to a set of springs left at its minimum energy equilibrium.\n",
    "\n",
    "<center><img src=\"spring_graph.png\" width=\"400\"></center>\n",
    "<div style=\"text-align: center\"> source: Dr Laleh Haghverdi - sc-sys-med Summer School (31 May, 2022) </div>\n",
    "\n",
    "#### Warning: speculative\n",
    "In t-SNE, the balance between attractive and repulsive forces are achieved via the minimization of the Kullback-Leibler divergence, whereas minimization in UMAP is achieved via cross-entropy. The latter directly accounts for repulsive forces between dis-similar points whereas the Kullback-Leibler divergence does not. This difference between cross-entropy (UMAP) and Kullback-Leibler divergence (t-SNE) perhaps accounts for the better preservation of global structure in UMAP, although other theoretically-grounded reasons are also at play.\n",
    "\n",
    "**Note**: t-SNE and UMAP are primarily used as *data visualization* methods and should generally not be used for other downstream purposes usually associated with dimensionality reduction such as with PCA.\n",
    "\n",
    "**Note**: Clustering is often performed separately using e.g. the Louvain or Leiden algorithm. Datapoints embedded within the low-dimensional t-SNE/UMAP space are then colour-coded by cluster label."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57fd668b",
   "metadata": {},
   "source": [
    "#### t-SNE and UMAP on scRNA-seq data\n",
    "t-SNE and UMAP are often used in the analysis of single-cell RNA sequencing (scRNA-seq) data which is very high dimensional, e.g. d = 20'000 genes. scRNA-seq data from the mouse prancreas will be used here for demonstration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d3f3b1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Import and format data\n",
    "\n",
    "random.seed(10)\n",
    "\n",
    "# mouse pancreas single-cell dataset\n",
    "# read in data and cell type labels\n",
    "with open('MP.pickle', 'rb') as f:\n",
    "    df = pickle.load(f)\n",
    "\n",
    "with open('MP_genes.pickle', 'rb') as f:\n",
    "    genes = pickle.load(f)\n",
    "\n",
    "df.set_index('Unnamed: 0', inplace=True)  # set first column (cell ID as the index column)\n",
    "sample_id = pickle.load(open('cell_IDs.pkl', 'rb'))\n",
    "df = df.loc[list(sample_id), :]\n",
    "\n",
    "X = df[genes].values  # extract the N x M cells-by-genes matrix\n",
    "\n",
    "sample_info = pd.read_csv('sample_info.csv')\n",
    "\n",
    "mp_anndata = anndata.AnnData(X=X)\n",
    "mp_anndata.obs['Celltype'] = sample_info['assigned_cluster'].values\n",
    "celltype_one_hot = pd.get_dummies(mp_anndata.obs['Celltype'], prefix='Celltype').values\n",
    "\n",
    "N = X.shape[0]  # number of single-cell samples\n",
    "K = len(sample_info['assigned_cluster'].unique())  # number of topics\n",
    "M = X.shape[1]  # number of genes\n",
    "\n",
    "display(mp_anndata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d5e0418",
   "metadata": {},
   "source": [
    "#### PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ff0d1c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute and plot 2D PCA projection\n",
    "_, ax = plt.subplots(figsize=(7, 6))\n",
    "\n",
    "sc.tl.pca(mp_anndata)\n",
    "sc.pl.pca(mp_anndata, color=[\"Celltype\"], ax=ax, title='two-dimensional PCA plot')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea5c1352",
   "metadata": {},
   "source": [
    "#### t-SNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b38c3bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute and plot 2D t-SNE projection\n",
    "_, ax = plt.subplots(figsize=(7, 6))\n",
    "\n",
    "sc.tl.tsne(mp_anndata, perplexity=30)\n",
    "sc.pl.tsne(mp_anndata, color=[\"Celltype\"], ax=ax, title='two-dimensional tSNE plot')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f9b2d2c",
   "metadata": {},
   "source": [
    "#### The perplexity parameter in t-SNE\n",
    "\n",
    "The *perplexity* parameter modulates the width of the Gaussian kernel used to quantify similarities between points based on their distance of separation. Larger perplexity values lead to wider Gaussian kernels, allowing greater similarity between points that are more distant from one another.\n",
    "\n",
    "<center><img src=\"tsne_perplexity.png\" width=\"800\"></center>\n",
    "<div style=\"text-align: center\"> source: Dr Laleh Haghverdi - sc-sys-med Summer School (31 May, 2022) </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4c858e0",
   "metadata": {},
   "source": [
    "#### UMAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fcfa0a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Compute neighbor affinities\n",
    "sc.pp.neighbors(mp_anndata, n_neighbors=30)\n",
    "\n",
    "## Compute and plot 2D UMAP projection\n",
    "_, ax = plt.subplots(figsize=(7, 6))\n",
    "\n",
    "sc.tl.umap(mp_anndata, min_dist=0.8)\n",
    "sc.pl.umap(mp_anndata, color=[\"Celltype\"], ax=ax, title='two-dimensional UMAP plot')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c80087e",
   "metadata": {},
   "source": [
    "#### Number of nearest neighbors and Minimum distance parameters in UMAP\n",
    "\n",
    "UMAP relies on the initial construction of a high-dimensional graph which is then embedded in low-dimensional space. This high-dimensional graph relies on the **number of nearest neighbors** parameter (n_neighbors) which controls how many adjacent datapoints are used to construct adaptive similarity kernels. To estimate the low-dimensional embedding, UMAP relies on the **minimum distance** parameter to \"re-calibrate\" the similarities' sensitivity to more distant points.\n",
    "\n",
    "In short, both parameters aim to balance local and global structures: *n_neighbors* does so during the initial contruction of the high-dimensional graph, whereas *min_distance* does so later during the search for a low-dimensional embedding. As such, altering *n_neighbors* has a more severe impact than *min_distance*.\n",
    "\n",
    "<center><img src=\"umap_n_neighbors_min_dist.png\" width=\"1000\"></center>\n",
    "<div style=\"text-align: center\"> source: https://pair-code.github.io/understanding-umap/ </div>\n",
    "\n",
    "**Note**: other parameters are at play for UMAP, but are considered less interpretable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cfea46b",
   "metadata": {},
   "source": [
    "#### Recommended reading and exercise:\n",
    "- Understanding UMAP (and t-SNE): https://pair-code.github.io/understanding-umap/\n",
    "\n",
    "- Fine-tuning UMAP visualizations: https://jlmelville.github.io/uwot/abparams.html"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0772f131",
   "metadata": {},
   "source": [
    "## Autoencoders\n",
    "\n",
    "#### A word about neural networks\n",
    "Neural networks combine linear and non-linear transformations to obtain powerful hidden representations of data. These hidden representations serve many purposes such as regression, classification, probability density estimation, image segmentation, etc.\n",
    "\n",
    "<center><img src=\"neuralnet.png\" width=\"400\"></center>\n",
    "<div style=\"text-align: center\"> source: https://en.wikipedia.org/wiki/Neural_network </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fce3d7a4",
   "metadata": {},
   "source": [
    "#### Autoencoders\n",
    "\n",
    "Within the context of dimensionality reduction, autoencoders are a subclass of neural networks which contains a hidden layer whose size is *smaller* than that of the input layer. Moreover, the output layer of an autoencoder aims at *reconstructing* the data provided at the input layer.\n",
    "\n",
    "<center><img src=\"autoencoder.png\" width=\"700\"></center>\n",
    "<div style=\"text-align: center\"> source: https://www.jeremyjordan.me/autoencoders/ </div>\n",
    "\n",
    "One can analyse the learned variables contained within the hidden layer. In this way, we achieve dimensionality reduction since these hidden variables are smaller in number than the original input variables."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34604ab2",
   "metadata": {},
   "source": [
    "#### Variational Autoencoders\n",
    "\n",
    "We've briefly discussed how PCA can be framed within a probabilistic setting, leading to Probabilistic PCA. Similarly, **variational autoencoders** (VAE) frame the autoencoder framework within a probabilistic setting. Rather than estimate hidden variables within the bottleneck layer, VAEs estimate *hidden probability distributions* from which hidden variables can be sampled from.\n",
    "\n",
    "<center><img src=\"vae.png\" width=\"700\"></center>\n",
    "<div style=\"text-align: center\"> source: https://www.jeremyjordan.me/autoencoders/ </div>\n",
    "\n",
    "These hidden probability distributions are often chosen to follow the Gaussian/Normal distribution. This design provides structure to the bottleneck layer while also accounting for variability which is inherent within the data. Defining probability distributions within the bottleneck layer also has interesting implications for data generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98170d3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import keras\n",
    "from keras import layers\n",
    "from keras.datasets import mnist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0c398b1",
   "metadata": {},
   "source": [
    "**Note**: exercises adapted from: https://www.theaidream.com/post/an-introduction-to-autoencoder-and-variational-autoencoder-vae"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c6a1644",
   "metadata": {},
   "source": [
    "#### MNIST data\n",
    "\n",
    "<center><img src=\"mnist.jpeg\" width=\"400\"></center>\n",
    "<div style=\"text-align: center\"> source: https://github.com/cazala/mnist </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2163369",
   "metadata": {},
   "source": [
    "#### Construct and train autoencoder model\n",
    "\n",
    "**Note**: There are three hidden layers in this model, where the middle bottleneck layer is called the *encoding* layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38d19162",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the size of our encoded representations\n",
    "encoding_dim = 2\n",
    "\n",
    "hidden_dim = 64\n",
    "\n",
    "# This is our input image\n",
    "input_img = keras.Input(shape=(784,))\n",
    "\n",
    "# \"encoded\" is the encoded representation of the input\n",
    "hidden_enc = layers.Dense(hidden_dim, activation='relu')(input_img)\n",
    "encoded = layers.Dense(encoding_dim, activation='relu')(hidden_enc)\n",
    "ae_encoder = keras.Model(input_img, encoded, name='encoder')\n",
    "\n",
    "# \"decoded\" is the lossy reconstruction of the input\n",
    "encoded_inputs = keras.Input(shape=(encoding_dim,), name='z_sampling')\n",
    "hidden_dec = layers.Dense(hidden_dim, activation='relu')(encoded_inputs)\n",
    "decoded = layers.Dense(784, activation='sigmoid')(hidden_dec)\n",
    "ae_decoder = keras.Model(encoded_inputs, decoded, name='decoder')\n",
    "\n",
    "output_img = ae_decoder(ae_encoder(input_img))\n",
    "\n",
    "autoencoder = keras.Model(input_img, output_img, name='ae')\n",
    "\n",
    "#Now let's train our autoencoder to reconstruct MNIST digits.\n",
    "#First, we'll configure our model to use a per-pixel binary crossentropy loss, and the Adam optimizer:\n",
    "autoencoder.compile(optimizer='adam', loss='binary_crossentropy')\n",
    "\n",
    "#Let's prepare our input data. We're using MNIST digits, and we're discarding the labels (since we're only interested in encoding/decoding the input images).\n",
    "(x_train, _), (x_test, y_test) = mnist.load_data()\n",
    "\n",
    "#We will normalize all values between 0 and 1 and we will flatten the 28x28 images into vectors of size 784.\n",
    "x_train = x_train.astype('float32') / 255.\n",
    "x_test = x_test.astype('float32') / 255.\n",
    "x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:])))\n",
    "x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))\n",
    "\n",
    "#Now let's train our autoencoder for 50 epochs:\n",
    "autoencoder.fit(x_train, x_train,\n",
    "                epochs=10,\n",
    "                batch_size=32,\n",
    "                shuffle=True,\n",
    "                validation_data=(x_test, x_test))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "270d9cc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Encode and decode some digits\n",
    "\n",
    "n = 10  # Number of digits to display\n",
    "\n",
    "encoded_imgs = ae_encoder.predict(x_test, verbose=0)\n",
    "decoded_imgs = ae_decoder.predict(encoded_imgs, verbose=0)\n",
    "\n",
    "plt.figure(figsize=(20, 4))\n",
    "for i in range(n):\n",
    "    # Display original\n",
    "    ax = plt.subplot(2, n, i + 1)\n",
    "    plt.imshow(x_test[i].reshape(28, 28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "\n",
    "    # Display reconstruction\n",
    "    ax = plt.subplot(2, n, i + 1 + n)\n",
    "    plt.imshow(decoded_imgs[i].reshape(28, 28))\n",
    "    plt.gray()\n",
    "    ax.get_xaxis().set_visible(False)\n",
    "    ax.get_yaxis().set_visible(False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82d1d56a",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot MNIST samples in bottleneck layer of autoencoder model\n",
    "\n",
    "x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))\n",
    "\n",
    "x_test_encoded = ae_encoder.predict(x_test, verbose=0)\n",
    "\n",
    "plt.figure(figsize=(7, 6))\n",
    "plt.scatter(x_test_encoded[:,0], x_test_encoded[:,1], c=y_test, cmap='tab10')\n",
    "plt.xlabel('autoencoder 1')\n",
    "plt.ylabel('autoencoder 2')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10614734",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Display a 2D manifold of the digits\n",
    "\n",
    "n = 15  # figure with 15x15 digits\n",
    "digit_size = 28\n",
    "figure = np.zeros((digit_size * n, digit_size * n))\n",
    "\n",
    "# We will sample n points\n",
    "grid_x = np.flip(np.linspace(0, 120, n))\n",
    "grid_y = np.linspace(0, 100, n)\n",
    "\n",
    "# Apply AE decoder along grid pattern\n",
    "for i, yi in enumerate(grid_x):\n",
    "    for j, xi in enumerate(grid_y):\n",
    "        z_sample = np.array([[xi, yi]])\n",
    "        x_decoded = ae_decoder.predict(z_sample, verbose=0)  # decoder\n",
    "        digit = x_decoded[0].reshape(digit_size, digit_size)\n",
    "        figure[i * digit_size: (i + 1) * digit_size,\n",
    "               j * digit_size: (j + 1) * digit_size] = digit\n",
    "        \n",
    "fig, ax = plt.subplots(1,1,figsize=(10, 10))\n",
    "ax.imshow(figure)\n",
    "\n",
    "ax.set_xticks( ticks=ax.get_xticks()[1:-1] , labels=np.linspace( grid_x[-1] , grid_x[0] , len(ax.get_xticks()[1:-1]) ))\n",
    "ax.set_yticks( ticks=ax.get_yticks()[1:-1] , labels=np.linspace( grid_y[-1] , grid_y[0] , len(ax.get_yticks()[1:-1]) ))\n",
    "\n",
    "ax.set_xlabel('autoencoder 1')\n",
    "ax.set_ylabel('autoencoder 2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cc01fe3",
   "metadata": {},
   "source": [
    "#### Construct and train variational autoencoder model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c267f3bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#First, here's our encoder network, mapping inputs to our latent distribution parameters:\n",
    "\n",
    "latent_dim = 2\n",
    "\n",
    "original_dim = 28 * 28\n",
    "intermediate_dim = 64\n",
    "\n",
    "inputs = keras.Input(shape=(original_dim,))\n",
    "h = layers.Dense(intermediate_dim, activation='relu')(inputs)\n",
    "z_mean = layers.Dense(latent_dim)(h)\n",
    "z_log_sigma = layers.Dense(latent_dim)(h)\n",
    "\n",
    "#We can use these parameters to sample new similar points from the latent space:\n",
    "from keras import backend as K\n",
    "\n",
    "def sampling(args):\n",
    "    z_mean, z_log_sigma = args\n",
    "    epsilon = K.random_normal(shape=(K.shape(z_mean)[0], latent_dim), mean=0., stddev=0.1)\n",
    "    return z_mean + K.exp(z_log_sigma) * epsilon\n",
    "\n",
    "z = layers.Lambda(sampling)([z_mean, z_log_sigma])\n",
    "\n",
    "#Finally, we can map these sampled latent points back to reconstructed inputs:\n",
    "# Create encoder\n",
    "vae_encoder = keras.Model(inputs, [z_mean, z_log_sigma, z], name='encoder')\n",
    "\n",
    "# Create decoder\n",
    "latent_inputs = keras.Input(shape=(latent_dim,), name='z_sampling')\n",
    "x = layers.Dense(intermediate_dim, activation='relu')(latent_inputs)\n",
    "outputs = layers.Dense(original_dim, activation='sigmoid')(x)\n",
    "\n",
    "vae_decoder = keras.Model(latent_inputs, outputs, name='decoder')\n",
    "\n",
    "# Instantiate VAE model\n",
    "outputs = vae_decoder(vae_encoder(inputs)[2])\n",
    "\n",
    "vae = keras.Model(inputs, outputs, name='vae_mlp')\n",
    "\n",
    "#We train the model using the end-to-end model, with a custom loss function: the sum of a reconstruction term, and the KL divergence regularization term.\n",
    "\n",
    "reconstruction_loss = keras.losses.binary_crossentropy(inputs, outputs)\n",
    "reconstruction_loss *= original_dim\n",
    "kl_loss = 1 + z_log_sigma - K.square(z_mean) - K.exp(z_log_sigma)\n",
    "kl_loss = K.sum(kl_loss, axis=-1)\n",
    "kl_loss *= -0.5\n",
    "vae_loss = K.mean(reconstruction_loss + kl_loss)\n",
    "vae.add_loss(vae_loss)\n",
    "\n",
    "vae.compile(optimizer='adam')\n",
    "\n",
    "(x_train, y_train), (x_test, y_test) = mnist.load_data()\n",
    "x_train = x_train.astype('float32') / 255.\n",
    "x_test = x_test.astype('float32') / 255.\n",
    "x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:])))\n",
    "x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:])))\n",
    "\n",
    "\n",
    "#We train our VAE on MNIST digits:\n",
    "vae.fit(x_train, x_train,\n",
    "        epochs=10,\n",
    "        batch_size=32,\n",
    "        validation_data=(x_test, x_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdf58822",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Encode and decode some digits\n",
    "\n",
    "x_test_encoded = vae_encoder.predict(x_test, verbose=0)\n",
    "x_test_encoded = x_test_encoded[0]\n",
    "\n",
    "plt.figure(figsize=(7, 6))\n",
    "plt.scatter(x_test_encoded[:,0], x_test_encoded[:,1], c=y_test, cmap='tab10')\n",
    "plt.xlabel('variational autoencoder 1')\n",
    "plt.ylabel('variational autoencoder 2')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0f5bd0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Display a 2D manifold of the digits\n",
    "\n",
    "n = 15  # figure with 15x15 digits\n",
    "digit_size = 28\n",
    "figure = np.zeros((digit_size * n, digit_size * n))\n",
    "\n",
    "# We will sample n points within [-15, 15]\n",
    "grid_x = np.flip(np.linspace(-3, 3, n))\n",
    "grid_y = np.linspace(-3, 3, n)\n",
    "\n",
    "# Apply VAE decoder along grid pattern\n",
    "for i, yi in enumerate(grid_x):\n",
    "    for j, xi in enumerate(grid_y):\n",
    "        z_sample = np.array([[xi, yi]])\n",
    "        x_decoded = vae_decoder.predict(z_sample, verbose=0)\n",
    "        digit = x_decoded[0].reshape(digit_size, digit_size)\n",
    "        figure[i * digit_size: (i + 1) * digit_size,\n",
    "               j * digit_size: (j + 1) * digit_size] = digit\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(10, 10))\n",
    "ax.imshow(figure)\n",
    "\n",
    "ax.set_xticks( ticks=ax.get_xticks()[1:-1] , labels=np.linspace( grid_x[-1] , grid_x[0] , len(ax.get_xticks()[1:-1]) ))\n",
    "ax.set_yticks( ticks=ax.get_yticks()[1:-1] , labels=np.linspace( grid_y[-1] , grid_y[0] , len(ax.get_yticks()[1:-1]) ))\n",
    "\n",
    "ax.set_xlabel('variational autoencoder 1')\n",
    "ax.set_ylabel('variational autoencoder 2')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86b1160e",
   "metadata": {},
   "source": [
    "#### What if the hidden layer of an autoencoder contains more than two hidden layers?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a331a60d",
   "metadata": {},
   "source": [
    "### Other topics in dimensionality reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dac234e",
   "metadata": {},
   "source": [
    "1. Dimensionality reduction via latent associations between datasets\n",
    "2. Tensor decompositions\n",
    "3. Dimensionality reduction of autocorrelated data\n",
    "4. Categorical and mixed data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a651d61",
   "metadata": {},
   "source": [
    "#### Latent associations between two datasets\n",
    "\n",
    "<center><img src=\"cca.png\" width=\"400\"></center>\n",
    "<div style=\"text-align: center\"> source: https://www.slideshare.net/blazf/canonical-correlation-analysis-319668 </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ded9f00d",
   "metadata": {},
   "source": [
    "#### Tensor decomposition\n",
    "\n",
    "<center><img src=\"tensor.png\" width=\"800\"></center>\n",
    "<div style=\"text-align: center\"> source: Hore et al. (2016) - Nature genetics </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "274e90b6",
   "metadata": {},
   "source": [
    "#### Autocorrelated data\n",
    "\n",
    "Many dimensionality reduction methods implicitly assume that the ordering of features is unimportant. However, this is not the case for autocorrelated signals such as EEG data where the ordering between features matters. Some methods such as Second-Order Blind Identification (SOBI) takes this ordering into account."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e2aa6e3",
   "metadata": {},
   "source": [
    "#### Categorical data\n",
    "<center><img src=\"mca.jpeg\" width=\"600\"></center>\n",
    "<div style=\"text-align: center\"> source: Bannister (2019) - PLoS ONE </div>\n",
    "\n",
    "\n",
    "For mixed data, where both numerical and categorical variables are present, one can use **Factor Analysis of Mixed Data** which combines PCA with Multiple Correspondance Analysis."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dimred_py",
   "language": "python",
   "name": "dimred_py"
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}