figures.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [DeepSphere]: a spherical convolutional neural network\n",
    "[DeepSphere]: https://github.com/SwissDataScienceCenter/DeepSphere\n",
    "\n",
    "[Nathanaël Perraudin](https://perraudin.info), [Michaël Defferrard](http://deff.ch), Tomasz Kacprzak, Raphael Sgier\n",
    "\n",
    "# Figures for the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy import sparse\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib\n",
    "import healpy as hp\n",
    "import pygsp\n",
    "import svgutils.compose as sc\n",
    "import IPython.display as ipd\n",
    "\n",
    "import hyperparameters\n",
    "from deepsphere import utils, plot, models\n",
    "import experimental.cnn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"\"\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (17, 5)  # (9, 4) for matplotlib notebook\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pathfig = './figures/'\n",
    "os.makedirs(pathfig, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Full sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[8,6])\n",
    "\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "G = utils.healpix_graph(nside=8, nest=True)\n",
    "G.plotting.update(vertex_size=10)\n",
    "G.plot(ax=ax,edges=False)\n",
    "\n",
    "# Get rid of the ticks\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.set_zticks([])\n",
    "# Get rid of the panes\n",
    "ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Get rid of the spines\n",
    "ax.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "ax.set_title('Healpix sampling, Nside=8')\n",
    "plt.savefig(os.path.join(pathfig, \"healpix_sampling_8.pdf\"), bbox_inches='tight')\n",
    "\n",
    "fig = plt.figure(figsize=[8,6])\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "G = utils.healpix_graph(nside=4, nest=True)\n",
    "G.plotting.update(vertex_size=20)\n",
    "G.plot(ax=ax)\n",
    "\n",
    "# Get rid of the ticks\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.set_zticks([])\n",
    "# Get rid of the panes\n",
    "ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Get rid of the spines\n",
    "ax.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "ax.set_title('Graph, full sphere, Nside=4')\n",
    "plt.savefig(os.path.join(pathfig, \"healpix_graph_4.pdf\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Half the sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nside = 4\n",
    "npoints = hp.nside2npix(nside)\n",
    "indexes = hp.reorder(np.array(list(range(npoints))),n2r=True)[:npoints//2]\n",
    "\n",
    "G = utils.healpix_graph(nside=nside, nest=True, indexes=indexes)\n",
    "G.plotting['elevation']=90\n",
    "G.plotting['azimuth']=0\n",
    "G.plotting.update(vertex_size=50)\n",
    "\n",
    "fig = plt.figure(figsize=[8,8])\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "# plt.cm.Blues_r\n",
    "\n",
    "# Highlight the node with a degree of 7 on the full sphere\n",
    "G2 = utils.healpix_graph(nside=nside, nest=True)\n",
    "snode = np.arange(0,G2.N)[G2.d==7]\n",
    "sindex = set(indexes)\n",
    "snode2 = [el for el in snode if el in sindex]\n",
    "hl_index = [np.argmin(np.abs(indexes-el)) for el in snode2]\n",
    "sig = np.zeros([G.N])\n",
    "sig[hl_index]=1\n",
    "\n",
    "G.plot_signal(1-sig, ax=ax,colorbar=False)\n",
    "# G.plot(ax=ax)\n",
    "# Get rid of the ticks\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.set_zticks([])\n",
    "# Get rid of the panes\n",
    "ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Get rid of the spines\n",
    "ax.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Remove the title\n",
    "# ax.set_title('Graph, half sphere, Nside=4')\n",
    "ax.set_title('')\n",
    "\n",
    "# Zoom in\n",
    "c = 0.6\n",
    "plt.axis([c*min(G.coords[:,0]), c*max(G.coords[:,0]), c*min(G.coords[:,1]), c*max(G.coords[:,1]) ])\n",
    "\n",
    "fig.savefig(os.path.join(pathfig, \"half_graph_{}.pdf\").format(nside), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Pooling\n",
    "\n",
    "A better figure is made in the `figure_pooling` notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 4\n",
    "\n",
    "index = np.arange(hp.nside2npix(order)) + 1\n",
    "mask = np.zeros_like(index, dtype=np.bool)\n",
    "mask[:order**2] = 1\n",
    "index *= mask\n",
    "index = index.astype(np.float)\n",
    "index[index==0] = hp.UNSEEN\n",
    "hp.mollview(index, title='', nest=True, cbar=False,cmap=None, xsize=1600)\n",
    "plt.savefig(os.path.join(pathfig, \"pooling-order4.pdf\"), bbox_inches='tight')\n",
    "\n",
    "order = 2\n",
    "index = np.arange(hp.nside2npix(order)) + 1\n",
    "mask = np.zeros_like(index, dtype=np.bool)\n",
    "mask[:order**2] = 1\n",
    "index *= mask\n",
    "index = index.astype(np.float)\n",
    "index[index==0] = hp.UNSEEN\n",
    "hp.mollview(index, title='', nest=True, cbar=False,cmap=None, xsize=1600)\n",
    "plt.savefig(os.path.join(pathfig, \"pooling-order2.pdf\"), bbox_inches='tight')\n",
    "\n",
    "order = 1\n",
    "index = np.arange(hp.nside2npix(order)) + 1\n",
    "mask = np.zeros_like(index, dtype=np.bool)\n",
    "mask[:order**2] = 1\n",
    "index *= mask\n",
    "index = index.astype(np.float)\n",
    "index[index==0] = hp.UNSEEN\n",
    "hp.mollview(index, title='', nest=True, cbar=False,cmap=None, xsize=1600)\n",
    "plt.savefig(os.path.join(pathfig, \"pooling-order1.pdf\"), bbox_inches='tight')\n",
    "\n",
    "index = np.array(list(range(12)))\n",
    "hp.mollview(index, title='', nest=True, cbar=False,cmap=None, xsize=1600)\n",
    "plt.savefig(os.path.join(pathfig, \"12parts.pdf\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "orders = [1,2,4]\n",
    "order_max = max(orders)\n",
    "npix = hp.nside2npix(order_max)\n",
    "index = np.zeros([npix])\n",
    "for order in orders:\n",
    "    index[:order**2] = index[:order**2]+1\n",
    "index.astype(np.float)\n",
    "index[index==0] = hp.UNSEEN\n",
    "hp.mollview(index, title='', nest=True, cbar=False,cmap=None, xsize=1600)\n",
    "plt.savefig(os.path.join(pathfig, \"part_sphere.pdf\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_ball(map_test1, cmap=plt.cm.gray_r, sub=None, vmin =-0.5, vmax=1.5):\n",
    "    cmap.set_under('w')\n",
    "    cmap.set_bad('lightgray')\n",
    "    dot_size=10\n",
    "    rot = (0,30,345)\n",
    "    hp.visufunc.orthview(map=map_test1, half_sky=True, title='', rot=rot, cmap=cmap, cbar=False, hold=True, nest=True, min=vmin, max=vmax, notext=True, sub=sub);\n",
    "    theta, phi = hp.pix2ang(hp.npix2nside(len(map_test1)), range(len(map_test1)), nest=True);\n",
    "    hp.projscatter(theta, phi, c='k', s=dot_size);\n",
    "    hp.graticule();\n",
    "    hp.graticule(dmer=360,dpar=360,alpha=1, rot=(0,0,15), local=True);\n",
    "    hp.graticule(dmer=360,dpar=360,alpha=1, rot=(0,0,195), local=True);\n",
    "\n",
    "orders = [1,2,4]\n",
    "order_max = max(orders)\n",
    "npix = hp.nside2npix(order_max)\n",
    "index = np.zeros([npix])\n",
    "for order in orders:\n",
    "    index[:order**2] = index[:order**2]+1\n",
    "index.astype(np.float)\n",
    "index[index==0] = hp.UNSEEN\n",
    "make_ball(index, cmap=plt.cm.RdBu_r, vmin=0, vmax=np.max(index))\n",
    "plt.savefig(os.path.join(pathfig, \"part_sphere2.pdf\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Fourier basis\n",
    "\n",
    "Let us display a few Fourier modes on the healpix map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_eigenvectors = 16\n",
    "\n",
    "G = utils.healpix_graph(nside=16, lap_type='normalized', nest=True, dtype=np.float64)\n",
    "G.compute_fourier_basis(n_eigenvectors=n_eigenvectors)\n",
    "\n",
    "fig = plt.figure(figsize=(8, 5))\n",
    "cm = plt.cm.RdBu_r\n",
    "cm.set_under('w')\n",
    "\n",
    "l, m = 0, 0\n",
    "lm = []\n",
    "for idx in range(n_eigenvectors):\n",
    "    lm.append([l,m])\n",
    "    m += 1\n",
    "    if m > l:\n",
    "        l += 1\n",
    "        m = -l\n",
    "ind = np.array([ 0,  1,  3,  2,  4,  5,  7, 6,  8, 10, 12, 9,  15, 14, 11, 13])\n",
    "\n",
    "for idx in range(n_eigenvectors):\n",
    "    l,m = lm[ind[idx]]\n",
    "    hp.mollview(G.U[:, idx], \n",
    "                title='Mode {}: $\\ell$={}, $|m|$={}'.format(idx, l, np.abs(m)),\n",
    "                nest=True, \n",
    "                sub=(np.sqrt(n_eigenvectors), np.sqrt(n_eigenvectors), idx+1),\n",
    "                max=np.max(np.abs(G.U[:, :n_eigenvectors])),\n",
    "                min=-np.max(np.abs(G.U[:, :n_eigenvectors])),\n",
    "                cbar=False,\n",
    "                cmap=cm)\n",
    "\n",
    "hp.graticule(verbose=False)\n",
    "plt.savefig(os.path.join(pathfig, \"eigenvectors.pdf\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 Convolution on graphs: show filters in spectral and spatial domains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# taus = [5, 10, 20, 50]\n",
    "taus = [5, 20, 50]\n",
    "matplotlib.rcParams.update({'font.size': 14})\n",
    "# fig, ax = plt.subplots(1,len(taus), figsize=(17, 4))\n",
    "fig, ax = plt.subplots(1,len(taus), figsize=(12, 4))\n",
    "for i,tau in enumerate(taus):\n",
    "    hf = pygsp.filters.Heat(G, tau=tau)\n",
    "    hf.plot(eigenvalues=False, sum=False, ax=ax[i])\n",
    "    ax[i].set_xlabel('Graph eigenvalues', fontsize=18)\n",
    "    if i is not 0:\n",
    "        ax[i].set_ylabel('')\n",
    "    else:\n",
    "        ax[i].set_ylabel('Spectral response', fontsize=18)\n",
    "    ax[i].set_title('$t={}$'.format(tau), fontsize=22)\n",
    "fig.tight_layout(rect=[0, 0.05, 1, 0.92])\n",
    "\n",
    "# plt.suptitle('Filter response in the graph spectral domain', fontsize=24)\n",
    "plt.savefig(os.path.join(pathfig, \"gaussian_filters_spectral.pdf\"), bbox_inches='tight')\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hf = pygsp.filters.Heat(G,tau=taus)\n",
    "\n",
    "def arcmin2rad(x):\n",
    "    return x / 60 / 360 * 2 * np.pi\n",
    "\n",
    "def gaussian_smoothing(sig, sigma, nest=True):\n",
    "    if nest:\n",
    "        sig = hp.reorder(sig, n2r=True)\n",
    "    smooth = hp.sphtfunc.smoothing(sig, sigma=arcmin2rad(sigma))\n",
    "    if nest:\n",
    "        smooth = hp.reorder(smooth, r2n=True)\n",
    "    return smooth\n",
    "_, center = plot.get_index_equator(hp.npix2nside(G.N), radius=20)\n",
    "ind0 = center\n",
    "sig = np.zeros(G.N)\n",
    "sig[ind0] = 1\n",
    "conv = hf.analyze(sig)\n",
    "\n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "rel_diff = []\n",
    "matplotlib.rcParams.update({'font.size': 18})\n",
    "cm = plt.cm.seismic\n",
    "# cm = plt.cm.jet\n",
    "cm.set_under('w')\n",
    "m = 0\n",
    "#[315, 465, 670, 1080]\n",
    "for i, (tau, sigma) in enumerate(zip(taus, [315, 670, 1080])):\n",
    "    with utils.HiddenPrints():\n",
    "        smooth = gaussian_smoothing(sig, sigma, nest=True)\n",
    "    m = max(m, max(smooth))\n",
    "    hp.mollview(conv[:, i], \n",
    "                title='$t={}$'.format(tau), \n",
    "                nest=True, \n",
    "                min=-m, max=m,\n",
    "                cbar=False,\n",
    "                rot=(180,0,180),\n",
    "                sub=(2, len(taus), i+1),\n",
    "                cmap=cm)\n",
    "    \n",
    "    hp.mollview(smooth, \n",
    "                title='$\\sigma={}$'.format(sigma), \n",
    "                nest=True, \n",
    "                min=-m, max=m,\n",
    "                cbar=False,\n",
    "                rot=(180,0,180),\n",
    "                sub=(2,len(taus),i+len(taus)+1),\n",
    "                cmap=cm)\n",
    "    \n",
    "    diff = (conv[:, i]-smooth)\n",
    "    rel_diff.append(np.linalg.norm(diff)/np.linalg.norm(smooth))\n",
    "#     hp.mollview(diff, \n",
    "#                 title='', \n",
    "#                 nest=True,\n",
    "#                 cbar=False,\n",
    "#                 sub=(3, len(taus), i+2*len(taus)+1))\n",
    "with utils.HiddenPrints():\n",
    "    hp.graticule();\n",
    "\n",
    "print(rel_diff)\n",
    "\n",
    "plt.savefig(os.path.join(pathfig, \"gaussian_filters_sphere.pdf\"), bbox_inches='tight')\n",
    "\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hf = pygsp.filters.Heat(G,tau=taus)\n",
    "order = 20\n",
    "matplotlib.rcParams.update({'font.size': 20})\n",
    "fig = plt.figure( figsize=(12, 5.5))\n",
    "plot.plot_filters_gnomonic(filters=hf,order=order, title='', graticule=True)\n",
    "\n",
    "# plt.suptitle('Gnomonic projection of a convoluted delta', fontsize=27)\n",
    "\n",
    "plt.savefig(os.path.join(pathfig, \"gaussian_filters_gnomonic.pdf\"), bbox_inches='tight')\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 14})\n",
    "fig = plt.figure( figsize=(12, 4))\n",
    "plot.plot_filters_section(hf, order=order, xlabel='', ylabel='', title='', marker='o')\n",
    "# plt.suptitle('Section of a convoluted delta', fontsize=22)\n",
    "plt.savefig(os.path.join(pathfig, \"gaussian_filters_section.pdf\"), bbox_inches='tight')\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot.plot_index_filters_section(hf,order=order)\n",
    "plt.savefig(os.path.join(pathfig, \"index_plotting_order{}_nside16.pdf\".format(order)), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Experimental results: classification of convergence maps in two model classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma=3\n",
    "deepsphere_result_fcn = np.load('results/deepsphere/deepsphere_results_list_sigma{}_FCN.npz'.format(sigma))['data'][-15:]\n",
    "deepsphere_result_cnn = np.load('results/deepsphere/deepsphere_results_list_sigma{}_CNN.npz'.format(sigma))['data'][-15:]\n",
    "hist_result = np.load('results/histogram/histogram_results_list_sigma{}.npz'.format(sigma))['data'][-15:]\n",
    "psd_result = np.load('results/psd/psd_results_list_sigma{}.npz'.format(sigma))['data']\n",
    "\n",
    "cnn2d_result_fcn_reg = np.load('results/2dcnn/reg3_deepsphere_results_list_sigma{}_FCN-2d.npz'.format(sigma))['data'][:15]\n",
    "cnn2d_result_cnn_reg = np.load('results/2dcnn/reg3_deepsphere_results_list_sigma{}_CNN-2d.npz'.format(sigma))['data'][:15]\n",
    "cnn2d_result_fcn_big_reg = np.load('results/2dcnn/reg3_deepsphere_results_list_sigma{}_FCN-2d-big.npz'.format(sigma))['data'][:15]\n",
    "cnn2d_result_cnn_big_reg = np.load('results/2dcnn/reg3_deepsphere_results_list_sigma{}_CNN-2d-big.npz'.format(sigma))['data'][:15]\n",
    "\n",
    "\n",
    "def get_xy(result, order):\n",
    "    x = []\n",
    "    y = []\n",
    "    for d in result:\n",
    "        if d[0]==order:\n",
    "            x.append(d[1])\n",
    "            y.append(d[2])\n",
    "    x = np.array(x)\n",
    "    y = np.array(y)\n",
    "    a = np.argsort(x)\n",
    "    x = x[a]\n",
    "    y = y[a]\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for order in[1, 2, 4]:\n",
    "    x_hist, y_hist = get_xy(hist_result, order)\n",
    "    x_deepsphere_fcn, y_deepsphere_fcn = get_xy(deepsphere_result_fcn, order)\n",
    "    x_deepsphere_cnn, y_deepsphere_cnn = get_xy(deepsphere_result_cnn, order)\n",
    "    x_psd, y_psd = get_xy(psd_result, order)\n",
    "    \n",
    "    x_cnn2d_fcn_reg, y_cnn2d_fcn_reg = get_xy(cnn2d_result_fcn_reg, order)\n",
    "    x_cnn2d_cnn_reg, y_cnn2d_cnn_reg = get_xy(cnn2d_result_cnn_reg, order)\n",
    "    x_cnn2d_fcn_big_reg, y_cnn2d_fcn_big_reg = get_xy(cnn2d_result_fcn_big_reg, order)\n",
    "    x_cnn2d_cnn_big_reg, y_cnn2d_cnn_big_reg = get_xy(cnn2d_result_cnn_big_reg, order)    \n",
    "    \n",
    "    acc_hist = (1-y_hist)*100\n",
    "    acc_deepsphere_fcn = (1-y_deepsphere_fcn)*100\n",
    "    acc_deepsphere_cnn = (1-y_deepsphere_cnn)*100\n",
    "    acc_psd = (1-y_psd)*100\n",
    "\n",
    "    acc_cnn2d_fcn_reg = (1-y_cnn2d_fcn_reg)*100\n",
    "    acc_cnn2d_cnn_reg = (1-y_cnn2d_cnn_reg)*100\n",
    "    acc_cnn2d_fcn_big_reg = (1-y_cnn2d_fcn_big_reg)*100\n",
    "    acc_cnn2d_cnn_big_reg = (1-y_cnn2d_cnn_big_reg)*100\n",
    "\n",
    "    plt.figure(figsize=[6, 4])\n",
    "    plt.plot(x_deepsphere_fcn, acc_deepsphere_fcn,'g.-', label='DeepSphere (FCN variant)')\n",
    "    plt.plot(x_deepsphere_cnn, acc_deepsphere_cnn,'g.--', label='DeepSphere (CNN variant)')\n",
    "    plt.plot(x_cnn2d_fcn_reg, acc_cnn2d_fcn_reg,'c.-', label='2D ConvNet (FCN variant)')\n",
    "    plt.plot(x_cnn2d_cnn_reg, acc_cnn2d_cnn_reg,'c.--', label='2D ConvNet (CNN variant)')\n",
    "#     plt.plot(x_cnn2d_fcn_big_reg, acc_cnn2d_fcn_big_reg,'m.-', label='2DCNN (FCN variant)')\n",
    "#     plt.plot(x_cnn2d_cnn_big_reg, acc_cnn2d_cnn_big_reg,'m.--', label='2DCNN (CNN variant)')\n",
    "    plt.plot(x_hist, acc_hist,'r.-', label='linear SVM on histogram')\n",
    "    plt.plot(x_psd, acc_psd,'b.-', label='linear SVM on PSD')\n",
    "    plt.legend(loc=3, prop={'size': 12})\n",
    "    plt.xlabel('Relative noise level')\n",
    "    plt.ylabel('Accuracy in %')\n",
    "    npix = (1024//order)**2\n",
    "    part = 12 * order**2\n",
    "    plt.title('Order $o={}$: {:,} pixels per samples (1/{} sphere)'.format(order, npix, part))\n",
    "    if order==1:\n",
    "        plt.ylim( (77, 101))\n",
    "    \n",
    "    plt.savefig(os.path.join(pathfig, \"result_order{}.pdf\".format(order)), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# deepsphere_result_params = np.load('results/deepsphere/deepsphere_results_list_sigma{}_params.npz'.format(sigma))['data']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def make_tab(order, results):\n",
    "#     print('-'*48)\n",
    "#     print('| {} | {} |'.format('Network'.ljust(30),'Accuracy % '))\n",
    "#     print('-'*48)\n",
    "#     for result in results:\n",
    "#         if int(result[0])==int(order):\n",
    "#             print('| {} |    {:0.2f}    |'.format(result[3].ljust(30), 100*(1-float(result[2]))))\n",
    "#     print('-'*48)\n",
    "\n",
    "# make_tab(4, deepsphere_result_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make_tab(2, deepsphere_result_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6 Experimental data: show two convergence maps, one per model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img1 = hp.read_map('data/same_psd/kappa_omega_m_0p31_s_2.fits')\n",
    "img2 = hp.read_map('data/same_psd/kappa_omega_m_0p26_s_2.fits')\n",
    "img1 = hp.reorder(img1, r2n=True)\n",
    "img2 = hp.reorder(img2, r2n=True)\n",
    "Nside = 1024\n",
    "img1 = hp.ud_grade(img1, nside_out=Nside, order_in='NESTED')\n",
    "img2 = hp.ud_grade(img2, nside_out=Nside, order_in='NESTED')\n",
    "cmin = min(np.min(img1), np.min(img2))\n",
    "cmax = max(np.max(img1), np.max(img2))\n",
    "cmax = -2*cmin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# _ = plt.hist(img1,bins=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# hp.mollview(img1, title='Map 1, omega_m=0.31, pk_norm=0.82, h=0.7', nest=True, min=cmin, max=cmax)\n",
    "# hp.mollview(img2, title='Map 2, omega_m=0.26, sigma_8=0.91, h=0.7', nest=True, min=cmin, max=cmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def arcmin2rad(x):\n",
    "    return x / 60 / 360 * 2 * np.pi\n",
    "\n",
    "def gaussian_smoothing(sig, sigma, nest=True):\n",
    "    if nest:\n",
    "        sig = hp.reorder(sig, n2r=True)\n",
    "    smooth = hp.sphtfunc.smoothing(sig, sigma=arcmin2rad(sigma))\n",
    "    if nest:\n",
    "        smooth = hp.reorder(smooth, r2n=True)\n",
    "    return smooth\n",
    "sigma=3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot.zoom_mollview(img1, cmin=cmin, cmax=cmax)\n",
    "plt.suptitle('Sample from class 1, $\\Omega_m=0.31$, , $\\sigma_8=0.82$',y=0.78, fontsize=18);\n",
    "# omega_m=0.31, pk_norm=0.82, h=0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot.zoom_mollview(gaussian_smoothing(img1,sigma), cmin=cmin, cmax=cmax)\n",
    "plt.suptitle('Smoothed map from class 1, $\\Omega_m=0.31$, $\\sigma_8=0.82$',y=0.78, fontsize=18)\n",
    "plt.savefig(os.path.join(pathfig, \"smooth_map_class_1.pdf\"), bbox_inches='tight')\n",
    "# omega_m=0.31, pk_norm=0.82, h=0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot.zoom_mollview(img2, cmin=cmin, cmax=cmax)\n",
    "_ = plt.suptitle('Sample from class 2, $\\Omega_m=0.26$, $\\sigma_8=0.91$',y=0.78, fontsize=18)\n",
    "# omega_m=0.26, sigma_8=0.91, h=0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot.zoom_mollview(gaussian_smoothing(img2, sigma), cmin=cmin, cmax=cmax)\n",
    "_ = plt.suptitle('Smoothed map from class 2, $\\Omega_m=0.26$, $\\sigma_8=0.91$',y=0.78, fontsize=18)\n",
    "plt.savefig(os.path.join(pathfig, \"smooth_map_class_2.pdf\"), bbox_inches='tight')\n",
    "# omega_m=0.26, sigma_8=0.91, h=0.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7 Power spectral densities (PSD) of the convergence maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compute = False\n",
    "if compute:\n",
    "    def psd(x):\n",
    "        '''Spherical Power Spectral Densities'''\n",
    "        hatx = hp.map2alm(hp.reorder(x, n2r=True))\n",
    "        return hp.alm2cl(hatx)\n",
    "    data_path = 'data/same_psd/'\n",
    "    ds1 = np.load(data_path+'smoothed_class1_sigma{}.npz'.format(sigma))['arr_0']\n",
    "    ds2 = np.load(data_path+'smoothed_class2_sigma{}.npz'.format(sigma))['arr_0']\n",
    "    psds_img1 = [psd(img) for img in ds1]\n",
    "    psds_img2 = [psd(img) for img in ds2]\n",
    "    np.savez('results/psd_data_sigma{}'.format(sigma), psd_class1=psds_img1, psd_class2=psds_img2)\n",
    "else:\n",
    "    psds_img1 = np.load('results/psd_data_sigma{}.npz'.format(sigma))['psd_class1']\n",
    "    psds_img2 = np.load('results/psd_data_sigma{}.npz'.format(sigma))['psd_class2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 14})\n",
    "\n",
    "l = np.array(range(len(psds_img1[0])))\n",
    "\n",
    "plot.plot_with_std(l,np.stack(psds_img1)*l*(l+1), label='class 1, $\\Omega_m=0.31$, $\\sigma_8=0.82$, $h=0.7$', color='r')\n",
    "plot.plot_with_std(l,np.stack(psds_img2)*l*(l+1), label='class 2, $\\Omega_m=0.26$, $\\sigma_8=0.91$, $h=0.7$', color='b')\n",
    "\n",
    "\n",
    "plt.legend(fontsize=16);\n",
    "plt.xlim([11, np.max(l)])\n",
    "plt.ylim([1e-6, 5e-4])\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel('$\\ell$: spherical harmonic index', fontsize=18)\n",
    "plt.ylabel('$C_\\ell \\cdot \\ell \\cdot (\\ell+1)$', fontsize=18)\n",
    "plt.title('Power Spectrum Density, 3-arcmin smoothing, noiseless, Nside=1024', fontsize=18);\n",
    "plt.savefig(os.path.join(pathfig, \"psd_sigma{}.pdf\".format(sigma)), bbox_inches='tight')\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8 Checking convergence of SVM (w.r.t. the number of augmented training samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 3\n",
    "order = 2\n",
    "sigma_noise = 1.5\n",
    "\n",
    "# path = 'results/psd/'\n",
    "# name = '40sim_1024sides_{0}arcmin_{2:.1f}noise_{1}order.npz'.format(sigma, order, sigma_noise)\n",
    "\n",
    "path = 'results/histogram/'\n",
    "name = '40sim_1024sides_{2}noise_{1}order_{0}sigma.npz'.format(sigma, order, sigma_noise)\n",
    "\n",
    "filepath = os.path.join(path,name)\n",
    "data = np.load(filepath)['arr_0']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 24})\n",
    "\n",
    "plt.plot(data[0], data[1], linewidth=4)\n",
    "plt.plot(data[0], data[2], linewidth=4)\n",
    "plt.plot(data[0][-1], data[3],'x', markersize=10)\n",
    "plt.legend(['Training','Validation', 'Testing'])\n",
    "plt.xlabel('Number of training samples')\n",
    "plt.ylabel('Error rate in %')\n",
    "# plt.title('Error for the histogram + SVM, order: {}, noise level: {}'.format(order, sigma_noise))\n",
    "plt.savefig(os.path.join(pathfig, \"hist_error_order{}_noise{}.pdf\".format(order, sigma_noise)), bbox_inches='tight')\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9 Plotting learned graph convolutional filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nside = 1024\n",
    "order = 2  # 1,2,4,8 correspond to 12,48,192,768 parts of the sphere.\n",
    "sigma_noise = 2\n",
    "sigma = 3\n",
    "ntype = 'FCN'\n",
    "EXP_NAME = '40sim_{}sides_{:0.1f}noise_{}order_{}sigma_{}'.format(Nside, sigma_noise, order, sigma, ntype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = hyperparameters.get_params(12*40*0.8*order*order, EXP_NAME, order, Nside, ntype)\n",
    "\n",
    "model = models.deepsphere(**params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folder = 'figures/filters/{}/'.format(EXP_NAME)\n",
    "os.makedirs(folder, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "layer = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    model.plot_chebyshev_coeffs(layer, ind_in=range(5), ind_out=range(10))\n",
    "    plt.savefig('{}/layer{}_coefficients.png'.format(folder, layer), dpi=100)\n",
    "except ValueError:\n",
    "    raise ValueError('If checkpoints/{} is empty, run the demo_part_sphere notebook first.'.format(EXP_NAME))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.plot_filters_spectral(layer, ind_in=range(5), ind_out=range(10));\n",
    "\n",
    "plt.savefig('{}/layer{}_spectral.png'.format(folder, layer), dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "model.plot_filters_section(layer, ind_in=range(6), ind_out=range(4), title='');\n",
    "plt.savefig(os.path.join(pathfig, \"section_filter_last.pdf\".format(order)), bbox_inches='tight')\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 12)\n",
    "model.plot_filters_gnomonic(layer, ind_in=range(6), ind_out=range(4), title='');\n",
    "plt.savefig(os.path.join(pathfig, \"gnonomic_filter_last.pdf\".format(order)), bbox_inches='tight', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (17, 5)  # (9, 4) for matplotlib notebook\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "\n",
    "model.plot_filters_section(1, ind_out=range(4), title='');\n",
    "fig.savefig('{}/layer{}_section.png'.format(folder, layer), dpi=100)\n",
    "plt.savefig(os.path.join(pathfig, \"section_filter_first.pdf\".format(order)), bbox_inches='tight')\n",
    "matplotlib.rcParams.update({'font.size': 10})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10 Border effect of the graph convolution (part of the sphere)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams['image.cmap'] = 'RdBu_r'\n",
    "nside = 16\n",
    "indexes = range(nside**2)\n",
    "\n",
    "G = utils.healpix_graph(nside=nside, indexes=indexes)\n",
    "G.estimate_lmax()\n",
    "\n",
    "tau = 30\n",
    "hf = pygsp.filters.Heat(G, tau=tau)\n",
    "\n",
    "index1 = 170\n",
    "index2 = 64+2*16+2*4+2\n",
    "sig1 = np.zeros([nside**2])\n",
    "sig2 = np.zeros([nside**2])\n",
    "sig1[index1] = 1\n",
    "sig2[index2] = 1\n",
    "\n",
    "sig1 = hf.filter(sig1)\n",
    "sig2 = hf.filter(sig2)\n",
    "\n",
    "m = max(np.max(sig1), np.max(sig2))\n",
    "limits = [-m, m]\n",
    "\n",
    "# sig = np.arange(nside**2)\n",
    "fig = plt.figure(figsize=[12,6])\n",
    "\n",
    "ax1 = fig.add_subplot(121, projection='3d')\n",
    "\n",
    "G.plot_signal(sig1, ax=ax1, colorbar=False,limits=limits)\n",
    "\n",
    "# Get rid of the ticks\n",
    "ax1.set_xticks([])\n",
    "ax1.set_yticks([])\n",
    "ax1.set_zticks([])\n",
    "# Get rid of the panes\n",
    "ax1.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax1.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax1.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Get rid of the spines\n",
    "ax1.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax1.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax1.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Zoom\n",
    "a = 0.35\n",
    "ax1.set_xlim(-a,a)\n",
    "ax1.set_ylim(-a,a)\n",
    "\n",
    "# Remove the title\n",
    "# ax.set_title('Graph, half sphere, Nside=4')\n",
    "ax1.set_title('', FontSize=16)\n",
    "\n",
    "ax1.view_init(elev=10, azim=45)\n",
    "\n",
    "ax2 = fig.add_subplot(122, projection='3d')\n",
    "\n",
    "G.plot_signal(sig2, ax=ax2, limits=limits, colorbar=False)\n",
    "\n",
    "# Get rid of the ticks\n",
    "ax2.set_xticks([])\n",
    "ax2.set_yticks([])\n",
    "ax2.set_zticks([])\n",
    "# Get rid of the panes\n",
    "ax2.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax2.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax2.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Get rid of the spines\n",
    "ax2.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax2.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "ax2.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))\n",
    "\n",
    "# Zoom\n",
    "a = 0.35\n",
    "ax2.set_xlim(-a,a)\n",
    "ax2.set_ylim(-a,a)\n",
    "\n",
    "# Remove the title\n",
    "# ax.set_title('Graph, half sphere, Nside=4')\n",
    "ax2.set_title('', FontSize=16)\n",
    "\n",
    "ax2.view_init(elev=10, azim=45)\n",
    "\n",
    "plt.tight_layout(pad=0)\n",
    "plt.savefig(os.path.join(pathfig, \"border_effects.pdf\"), bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for nside in [16, 32, 64, 128, 256, 512, 1024, 2048]:\n",
    "#     print('Time for nside: {}'.format(nside))\n",
    "#     %timeit G = utils.healpix_graph(nside=nside)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11 Natural projection of Healpix to the plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nside = 8\n",
    "indexes = range(nside**2)\n",
    "G = utils.healpix_graph(nside=nside, indexes=indexes)\n",
    "sig1 = np.arange(nside**2)\n",
    "fig = plt.figure(figsize=[12,6])\n",
    "\n",
    "ax1 = fig.add_subplot(121, projection='3d')\n",
    "ax1.axis('off')\n",
    "\n",
    "G.plot_signal(sig1, ax=ax1, colorbar=False, edges=True)\n",
    "\n",
    "# Zoom\n",
    "a = 0.3\n",
    "ax1.set_xlim(-a, a)\n",
    "ax1.set_ylim(-a, a)\n",
    "\n",
    "ax1.view_init(elev=10, azim=45)\n",
    "ax1.set_title('Pixels on HEALPix (1/12 of the sphere)')\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "ind = experimental.cnn.build_index(3).flatten().astype(np.int)\n",
    "ax2.imshow(np.flipud(sig1[ind].reshape([nside, nside])))\n",
    "ax2.set_title('Pixels on the plane')\n",
    "ax2.get_xaxis().set_visible(False)\n",
    "ax2.get_yaxis().set_visible(False)\n",
    "\n",
    "fig.savefig(os.path.join(pathfig, 'projection_healpix_plane.pdf'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12 Filtering speed\n",
    "\n",
    "Numbers measured in the `spherical_vs_graph` notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = np.load('results/filtering_speed.npz')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5.9, 3.8))\n",
    "npix = [hp.nside2npix(nside) for nside in results['nsides']]\n",
    "\n",
    "ax.loglog(npix, results['times_sphere'][:, ::-1], '--', marker='.')\n",
    "ax.loglog(npix, results['times_graph'][:, ::-1], '-', marker='.')\n",
    "ax.loglog(npix, results['times_part'][:, 1, 0], '-.', marker='.')\n",
    "ax.loglog(npix, results['times_part'][:, 1, 2], '-.', marker='.')\n",
    "\n",
    "#ax.loglog(npix, np.array(npix)/1e6, ':', color='#808080')\n",
    "#ax.loglog(npix, (np.array(npix)/1e6)**1.5, ':', color='#808080')\n",
    "\n",
    "labels = ['Sph. harm., $\\ell_{{max}}$ = {}$N_{{side}}$'.format(lm) for lm in results['lmax'][::-1]]\n",
    "labels += ['Graph, poly. order K={}'.format(order) for order in results['orders'][::-1]]\n",
    "labels += ['Partial graph 1/{}, K={}'.format(v, 15) for v in [12, 192]]\n",
    "\n",
    "#labels += [r'Asymptotic $\\mathcal{O}(N_{side})$']\n",
    "#labels += [r'Asymptotic $\\mathcal{O}(N_{side}^{3/2})$']\n",
    "ax.legend(labels, loc='upper left')\n",
    "for i, nside in enumerate(results['nsides']):\n",
    "    x = npix[i]\n",
    "    y = results['times_sphere'][i, -1] * 2\n",
    "    ax.text(x, y, '$N_{{side}}$\\n{}'.format(nside), horizontalalignment='center')\n",
    "\n",
    "ax.set_ylim(0.6 * results['times_part'][:,1,:].min(), 8e4)#6e2 * results['times_sphere'].max())\n",
    "#ax.set_xlim(0.5 * min(npix), 2 * max(npix))\n",
    "ax.set_xlabel('Number of pixels')\n",
    "ax.set_ylabel('Processing time [s]')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(pathfig, 'filtering_speed.pdf'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 13 Eigenvalues of graph Laplacian are grouped (like the degree of the spherical harmonics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_eigenvalues = 50\n",
    "nside = 16\n",
    "\n",
    "graph = utils.healpix_graph(nside=nside, lap_type='normalized', nest=True, dtype=np.float64)\n",
    "graph.compute_fourier_basis(n_eigenvectors=n_eigenvalues)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 2.5))\n",
    "ax.plot(graph.e, '.-')\n",
    "idx = 1\n",
    "xticks = [idx]\n",
    "for l in range(1, 7):\n",
    "    ax.text(idx + l - 2.3, graph.e[idx + l] + 0.005, '$\\ell$ = {}'.format(l))\n",
    "    idx += 2*l + 1\n",
    "    xticks.append(idx)\n",
    "ax.set_xlabel('Eigenvalue $\\lambda$')\n",
    "ax.set_ylabel('Value')\n",
    "ax.set_xticks(xticks)\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(pathfig, 'graph_eigenvalues.pdf'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 14 Correspondance of subspaces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 14.1 Graph eigenvectors vs sampled spherical harmonics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nside = 16\n",
    "lmax = 8\n",
    "\n",
    "def compute_orhtogonality(nside, lmax, lap_type='normalized'):\n",
    "    n_harmonics = np.cumsum(np.arange(1, 2*lmax+2, 2))\n",
    "\n",
    "    harmonics = utils.compute_spherical_harmonics(nside, lmax=lmax)\n",
    "\n",
    "    graph = utils.healpix_graph(nside, lap_type=lap_type, nest=True, dtype=np.float64)\n",
    "    n_eigenvectors = min(n_harmonics[-1], graph.N)\n",
    "    graph.compute_fourier_basis(n_eigenvectors)\n",
    "\n",
    "    C = harmonics.T @ graph.U\n",
    "    \n",
    "    return C, n_harmonics\n",
    "\n",
    "C, n_harmonics = compute_orhtogonality(nside, lmax)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(5, 4))\n",
    "im = ax.imshow(np.abs(C), cmap=plt.cm.gist_heat_r, aspect='equal')\n",
    "ax.set_xlabel('Graph Fourier modes')\n",
    "ax.set_ylabel('Spherical harmonics')\n",
    "ax.set_xticks(n_harmonics - 0.5)\n",
    "ax.set_yticks(n_harmonics - 0.5)\n",
    "ax.set_xticklabels(n_harmonics)\n",
    "ax.set_yticklabels(n_harmonics)\n",
    "for l in range(4, lmax+1):\n",
    "    ax.text(n_harmonics[l-1] + l - 3.9, n_harmonics[l-1] - 1, '$\\ell={}$'.format(l))\n",
    "ax.grid(True)\n",
    "fig.colorbar(im)\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(pathfig, 'subspace_harmonics_eigenvectors.pdf'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 14.2 SHT of the graph eigenvectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsides = [4, 8, 16]\n",
    "\n",
    "spectral_content = dict()\n",
    "\n",
    "for nside in nsides:\n",
    "\n",
    "    lmax = 3 * nside - 1\n",
    "\n",
    "    n_harmonics = np.cumsum(np.arange(1, 2*lmax+2, 2))[-1]\n",
    "    graph = utils.healpix_graph(nside, lap_type='normalized', nest=True, dtype=np.float64)\n",
    "    graph.compute_fourier_basis(n_eigenvectors=n_harmonics)\n",
    "\n",
    "    cl = np.empty((n_harmonics, lmax+1))\n",
    "    for i in range(n_harmonics):\n",
    "        eigenvector = hp.reorder(graph.U[:, i], n2r=True)\n",
    "        cl[i] = hp.sphtfunc.anafast(eigenvector, lmax=lmax)\n",
    "\n",
    "    spectral_content[nside] = np.empty((lmax+1, lmax+1))\n",
    "    start = 0\n",
    "    for ell in range(lmax+1):\n",
    "        end = start + (2 * ell + 1)\n",
    "        spectral_content[nside][ell] = np.sum(cl[start:end,:], axis=0)\n",
    "        start = end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fontsize = matplotlib.rcParams['font.size']\n",
    "matplotlib.rcParams['font.size'] = 16 \n",
    "\n",
    "#fig, axes = plt.subplots(1, len(nsides), figsize=(6, 3))\n",
    "matplotlib.rcParams.update({'font.size': 16})\n",
    "fig, axes = plt.subplots(1, len(nsides), figsize=(12, 6))\n",
    "\n",
    "for ax, (nside, sc) in zip(axes.flatten(), spectral_content.items()):\n",
    "    lmax = 3*nside-1\n",
    "    sc = sc / sc[0, 0]\n",
    "    im = ax.imshow(sc, cmap=plt.cm.gist_heat_r)\n",
    "    #ax.set_title(rf'$N_{{side}}={nside}$')\n",
    "    ax.text(0.35, 0.8, rf'$N_{{side}}={nside}$', fontsize=20, transform=ax.transAxes)\n",
    "    #ax.set_xlabel('SH degree $\\ell$')\n",
    "    ax.text(0.5, -0.1, 'SH degree $\\ell$', horizontalalignment='center', transform=ax.transAxes)\n",
    "    #ax.yaxis.set_visible(False)\n",
    "    #ax.set_yticklabels([])\n",
    "    ax.set_yticks([])\n",
    "    #ax.set_yticks([0, lmax])\n",
    "    ax.set_xticks([0, lmax])\n",
    "\n",
    "axes[0].set_yticks([0])\n",
    "#axes[0].set_ylabel('subspaces of graph eigenvectors')\n",
    "axes[0].set_ylabel('graph eigenvectors', fontsize=20)\n",
    "\n",
    "fig.subplots_adjust(wspace=0.1, right=0.9)\n",
    "ax_cbar = fig.add_axes([0.92, 0.25, 0.02, 0.5])\n",
    "fig.colorbar(im, ax=ax, cax=ax_cbar, ticks=[0, 1])\n",
    "\n",
    "#fig.tight_layout()  # Doesn't work with color bar.\n",
    "fig.savefig(os.path.join(pathfig, 'subspace_harmonics_eigenvectors_v2.pdf'));\n",
    "\n",
    "matplotlib.rcParams['font.size'] = fontsize;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 15 Convolution basis: Chebyshev polynomials vs monomials\n",
    "\n",
    "1. Orthogonality of the basis in the spectral domain.\n",
    "1. Orthogonality of the basis in the vertex domain.\n",
    "1. Expected shape of the filters given a distribution over the coefficients.\n",
    "\n",
    "Todo:\n",
    "* compute the expectation analytically"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams.update({'font.size': 10})\n",
    "\n",
    "# Order of Chebyshev polynomials. Degree of monomials\n",
    "degree = 7\n",
    "\n",
    "n_points = 1000\n",
    "\n",
    "graph = pygsp.graphs.Path(64)\n",
    "\n",
    "# Irregular graph. Otherwise the Chebyshev polynomials are exactly orthogonal.\n",
    "graph.W.data = 0.5 + 0.1 * np.random.uniform(size=graph.W.data.shape)\n",
    "graph = pygsp.graphs.Graph(pygsp.utils.symmetrize(graph.W))\n",
    "#plt.imshow(graph.W.toarray())\n",
    "\n",
    "graph.estimate_lmax()\n",
    "graph.set_coordinates('line1D')\n",
    "\n",
    "fig = plt.figure(figsize=(8, 5))\n",
    "\n",
    "# Chebyshev\n",
    "\n",
    "#x = np.linspace(0, 1.05*graph.lmax, 1000)\n",
    "x = np.linspace(0, graph.lmax, n_points)\n",
    "\n",
    "coefficients = np.identity(degree)\n",
    "f = pygsp.filters.Chebyshev(graph, coefficients)\n",
    "Y = f.evaluate(x)\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (0, 0), colspan=2)\n",
    "ax.plot(x / graph.lmax * 2 - 1, Y.T)\n",
    "ax.legend(['k={}'.format(k) for k in range(degree)])\n",
    "ax.set_xlabel('Eigenvalue $\\lambda$')\n",
    "ax.set_ylabel('Polynomial $T_k(\\lambda)$')\n",
    "ax.set_title('Chebyshev basis (spectral domain)')\n",
    "ax.set_xticks([-1, 0, 1])\n",
    "ax.set_yticks([-1, 0, 1])\n",
    "ax.grid()\n",
    "\n",
    "C = Y @ Y.T\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (0, 2))\n",
    "im = ax.imshow(np.abs(C), cmap=plt.cm.gist_heat_r)\n",
    "fig.colorbar(im, ax=ax)\n",
    "ax.set_title('Cross-correlation')\n",
    "\n",
    "# Monomials\n",
    "\n",
    "x = np.linspace(-1, 1, n_points)\n",
    "\n",
    "Y = np.empty((degree, len(x)))\n",
    "for k in range(degree):\n",
    "    Y[k] = x**k\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (1, 0), colspan=2)\n",
    "plt.plot(x, Y.T)\n",
    "ax.legend(['k={}'.format(k) for k in range(degree)])\n",
    "ax.set_xlabel('Eigenvalue $\\lambda$')\n",
    "ax.set_ylabel('Monomial $\\lambda^k$')\n",
    "ax.set_title('Monomial basis (spectral domain)')\n",
    "ax.set_xticks([-1, 0, 1])\n",
    "ax.set_yticks([-1, 0, 1])\n",
    "ax.grid()\n",
    "\n",
    "C = Y @ Y.T\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (1, 2))\n",
    "im = ax.imshow(np.abs(C), cmap=plt.cm.gist_heat_r)\n",
    "fig.colorbar(im, ax=ax)\n",
    "ax.set_title('Cross-correlation')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(pathfig, 'polynomial_bases_spectrum.pdf'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 5))\n",
    "\n",
    "# Chebyshev\n",
    "\n",
    "Y = f.localize(graph.N // 2)\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (0, 0), colspan=2, fig=fig)\n",
    "for k in range(degree):\n",
    "    graph.plot_signal(Y[k], ax=ax)\n",
    "ax.legend(['k={}'.format(k) for k in range(degree)])\n",
    "ax.set_ylim(1.1*Y.min(), 1.1*Y.max())\n",
    "ax.set_xlim(graph.N // 2 - degree, graph.N // 2 + degree)\n",
    "ax.set_xticks(np.arange(graph.N // 2 - degree + 1, graph.N // 2 + degree, 2))\n",
    "ax.set_xticklabels('$v_{{{}}}$'.format(i) for i in range(- degree + 1, degree, 2))\n",
    "ax.set_title('Chebyshev basis (localized on vertex $v_0$)')\n",
    "ax.set_ylabel('($T_k(L) \\delta_0)_j$')\n",
    "\n",
    "C = Y @ Y.T\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (0, 2), fig=fig)\n",
    "im = ax.imshow(np.abs(C), cmap=plt.cm.gist_heat_r)\n",
    "fig.colorbar(im, ax=ax)\n",
    "ax.set_title('Cross-correlation')\n",
    "\n",
    "# Monomials\n",
    "\n",
    "Y = np.empty((degree, graph.N))\n",
    "s = np.zeros(graph.N)\n",
    "s[graph.N // 2] = 1\n",
    "L = graph.L / graph.lmax * 2 - sparse.identity(graph.N)\n",
    "for k in range(degree):\n",
    "    Y[k] = L**k @ s\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (1, 0), colspan=2, fig=fig)\n",
    "for k in range(degree):\n",
    "    graph.plot_signal(Y[k], ax=ax)\n",
    "ax.legend(['k={}'.format(k) for k in range(degree)])\n",
    "ax.set_ylim(1.1*Y.min(), 1.1*Y.max())\n",
    "ax.set_xlim(graph.N // 2 - degree, graph.N // 2 + degree)\n",
    "ax.set_xticks(np.arange(graph.N // 2 - degree + 1, graph.N // 2 + degree, 2))\n",
    "ax.set_xticklabels('$v_{{{}}}$'.format(i) for i in range(- degree + 1, degree, 2))\n",
    "ax.set_title('Monomial basis (localized on vertex $v_0$)')\n",
    "ax.set_ylabel('($L^k \\delta_0)_j$')\n",
    "\n",
    "C = Y @ Y.T\n",
    "\n",
    "ax = plt.subplot2grid((2, 3), (1, 2), fig=fig)\n",
    "im = ax.imshow(np.abs(C), cmap=plt.cm.gist_heat_r)\n",
    "fig.colorbar(im, ax=ax)\n",
    "ax.set_title('Cross-correlation')\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(pathfig, 'polynomial_bases_vertex.pdf'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "degrees = [5, 20, 100]\n",
    "n_realizations = int(1e4)\n",
    "n_points = 100\n",
    "\n",
    "x = np.linspace(-1, 1, n_points)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, sharey=True, figsize=(8.5, 3))\n",
    "\n",
    "for degree in degrees:\n",
    "    coefficients = np.random.normal(0, 1, size=(degree, n_realizations))\n",
    "    #coefficients = np.random.uniform(-1, 1, size=(order, n_realizations))\n",
    "\n",
    "    # Monomials.\n",
    "    y = np.zeros((n_realizations, n_points))\n",
    "    for k, c in enumerate(coefficients):\n",
    "        y += np.outer(c, x**k)\n",
    "    plot.plot_with_std(x, y, ax=axes[0])\n",
    "\n",
    "    # Chebyshev polynomials.\n",
    "    graph = pygsp.graphs.Path(n_points)\n",
    "    graph.estimate_lmax()\n",
    "    filters = pygsp.filters.Chebyshev(graph, coefficients)\n",
    "    y = filters.evaluate((x + 1) / 2 * graph.lmax)\n",
    "    plot.plot_with_std(x, y, ax=axes[1])\n",
    "\n",
    "legend = ['degree $K={}$'.format(degree) for degree in degrees]\n",
    "axes[0].legend(legend, loc='upper center')\n",
    "axes[1].legend(legend, loc='upper center')\n",
    "axes[0].set_xlabel('Scaled eigenvalue $x$')\n",
    "axes[1].set_xlabel('Scaled eigenvalue $x$')\n",
    "axes[0].set_ylabel(r'Expected filter value $\\mathbf{E}_\\theta[ g_\\theta(x) ]$')\n",
    "axes[0].set_title('Expected sum of monomials')\n",
    "axes[1].set_title('Expected sum of Chebyshev polynomials')\n",
    "axes[0].text(0, -7, r'$g_\\theta(x) = \\sum_{k=0}^K \\theta_k x^k$', horizontalalignment='center')\n",
    "axes[1].text(0, -6, r'$g_\\theta(x) = \\sum_{k=0}^K T_k(x)$', horizontalalignment='center')\n",
    "axes[1].text(0, -9.5, r'$T_k(x) = 2xT_{k-1}(x) - T_{k-2}(x), T_1(x) = x, T_0(x) = 1$', horizontalalignment='center')\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(pathfig, 'expected_filters.pdf'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x = np.arange(-1,1,0.001)\n",
    "\n",
    "# order = 20\n",
    "# c = np.random.randn(order,100)\n",
    "# f = []\n",
    "# for coeffs in c.T:\n",
    "#     s = 0*x\n",
    "#     for o, coeff in enumerate(coeffs):\n",
    "#         s += coeff*(x**o)\n",
    "#     f.append(s)\n",
    "# f = np.array(f)\n",
    "# ax = plot.plot_with_std(x, f)\n",
    "# ax.set_title('Monomial - order {}'.format(order));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# x = np.arange(-1,1,0.001)\n",
    "\n",
    "# order = 20\n",
    "# c = np.random.randn(order,100)\n",
    "# f = []\n",
    "# p = []\n",
    "# p.append(x**0)\n",
    "# p.append(x**1)\n",
    "# for o in range(2, order):\n",
    "#     p.append(2*x*p[o-1]-p[o-2])\n",
    "# for coeffs in c.T:\n",
    "#     s = x**0\n",
    "#     for o, coeff in enumerate(coeffs):\n",
    "#         s += coeff*p[o]\n",
    "#     f.append(s)\n",
    "# f = np.array(f)\n",
    "# ax = plot.plot_with_std(x, f)\n",
    "# ax.set_title('Chebyshev - order {}'.format(order));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(-1,1,0.001)\n",
    "order =20\n",
    "p = []\n",
    "p.append(x**0)\n",
    "p.append(x**1)\n",
    "for o in range(2,order):\n",
    "    p.append(2*x*p[o-1]-p[o-2])\n",
    "\n",
    "for o in range(order):\n",
    "    plt.plot(x, p[o])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for o in range(5,12):\n",
    "    plt.plot(x, np.sum(np.array(p[0:o])**2/(o+0.5)*2,axis=0))\n",
    "plt.plot(x, x**0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "o = 10\n",
    "plt.plot(x, np.sum(np.array(p[0:o])**2,axis=0))\n",
    "plt.plot(x, (o+0.5)/2*(x**0))"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}