figure_kappa_maps_difference.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",
    "# Figure: comparison of mass maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os; \n",
    "import healpy as hp; import numpy as np; \n",
    "from IPython.core.display import display, HTML; display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "%pylab inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl = plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pathfig = './figures/'\n",
    "os.makedirs(pathfig, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_zoom(map_test1, cmap=cm.gray_r, title='', rot = (0,45,-20), dtheta=5, clim = [-0.025, 0.041]):\n",
    "    edge_arcmin = 2*dtheta*60\n",
    "    n_pix = 500\n",
    "    reso = edge_arcmin/float(n_pix)\n",
    "    hp.gnomview(map_test1, reso=reso, xsize=n_pix, notext=True, title=title, rot=(0,0,0), cmap=cmap, min=clim[0], max=clim[1], hold=True)\n",
    "    hp.graticule();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_ball(map_test1, cmap=cm.gray_r, clim=None, title='', rot = (0,45,-20), dtheta = 5):\n",
    "    cmap.set_over('w')\n",
    "    cmap.set_under('w')\n",
    "    cmap.set_bad('lightgrey')\n",
    "    dot_size=10\n",
    "    if clim==None:\n",
    "        hp.visufunc.orthview(map=map_test1, half_sky=True, title=title, rot=rot, cmap=cmap, cbar=True, hold=True, nest=False, notext=True);\n",
    "    else:\n",
    "        hp.visufunc.orthview(map=map_test1, half_sky=True, title=title, rot=rot, cmap=cmap, cbar=True, hold=True, nest=False, notext=True, min=clim[0], max=clim[1]);\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",
    "    \n",
    "    hp.graticule();\n",
    "    hp.graticule(dmer=360,dpar=360,alpha=1, rot=(0,0,180-rot[2]), local=True);\n",
    "    hp.graticule(dmer=360,dpar=360,alpha=1, rot=(0,0,-1*rot[2]), local=True);\n",
    "        \n",
    "    lw=2\n",
    "    \n",
    "    theta = np.linspace(90-dtheta, 90+dtheta, 100)*np.pi/180.\n",
    "    phi = np.linspace(-dtheta, -dtheta, 100)*np.pi/180.\n",
    "    hp.projplot(theta, phi, c='k', lw=lw);\n",
    "\n",
    "    theta = np.linspace(90-dtheta, 90+dtheta, 100)*np.pi/180.\n",
    "    phi = np.linspace(dtheta, dtheta, 100)*np.pi/180.\n",
    "    hp.projplot(theta, phi, c='k', lw=lw);\n",
    "\n",
    "    theta = np.linspace(90-dtheta, 90-dtheta, 100)*np.pi/180.\n",
    "    phi = np.linspace(-dtheta, dtheta, 100)*np.pi/180.\n",
    "    hp.projplot(theta, phi, c='k', lw=lw);\n",
    "\n",
    "    theta = np.linspace(90+dtheta, 90+dtheta, 100)*np.pi/180.\n",
    "    phi = np.linspace(-dtheta, dtheta, 100)*np.pi/180.\n",
    "    hp.projplot(theta, phi, c='k', lw=lw);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# map_kappa1 = np.load('./data/data_v5/kappa_omega_m_0p31_s_1.npy')\n",
    "# map_kappa1 -= np.mean(map_kappa1)\n",
    "# map_kappa2 = np.load('./data/data_v5/kappa_omega_m_0p26_s_1.npy')\n",
    "# map_kappa2 -= np.mean(map_kappa2)\n",
    "\n",
    "map_kappa1 = hp.read_map('data/same_psd/kappa_omega_m_0p31_s_2.fits')\n",
    "map_kappa2 = hp.read_map('data/same_psd/kappa_omega_m_0p26_s_2.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_hires, sigma_lores = 5, 60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_kappa1_lores = hp.ud_grade(map_kappa1, nside_out=512)\n",
    "sigma_rad = sigma_lores / 60. * np.pi / 180.\n",
    "map_kappa1_lores = hp.smoothing(map_kappa1_lores, fwhm=sigma_rad)\n",
    "\n",
    "map_kappa1_hires = hp.ud_grade(map_kappa1, nside_out=1024)\n",
    "sigma_rad = sigma_hires / 60. * np.pi / 180.\n",
    "map_kappa1_hires = hp.smoothing(map_kappa1_hires, fwhm=sigma_rad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_kappa2_lores = hp.ud_grade(map_kappa2, nside_out=512)\n",
    "sigma_rad = sigma_lores / 60. * np.pi / 180.\n",
    "map_kappa2_lores = hp.smoothing(map_kappa1_lores, fwhm=sigma_rad)\n",
    "\n",
    "map_kappa2_hires = hp.ud_grade(map_kappa2, nside_out=1024)\n",
    "sigma_rad = sigma_hires / 60. * np.pi / 180.\n",
    "map_kappa2_hires = hp.smoothing(map_kappa2_hires, fwhm=sigma_rad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_kappa1_hires.min(), map_kappa1_hires.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_kappa2_hires.min(), map_kappa2_hires.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clim = [-0.014, 0.029]\n",
    "clim_hires = [-0.03, 0.04]\n",
    "rot = (0,-20,-20) \n",
    "dtheta = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx, ny, fs = 2, 2, 3\n",
    "pl.figure(figsize=(nx*fs, ny*fs))\n",
    "make_ball(map_kappa1_lores, cmap=pl.cm.plasma, clim=clim, rot=rot, dtheta=dtheta, title=r'Model 1: $\\Omega_m=0.31$ $\\sigma_8=0.82$' + '\\nsmoothing 1 deg')\n",
    "filename_plot = os.path.join(pathfig,'figure_kappa_diff.part1.svg')\n",
    "pl.savefig(filename_plot, bbox_inches='tight', transparent=True); print('saved {}'.format(filename_plot))\n",
    "\n",
    "pl.figure(figsize=(nx*fs, ny*fs))\n",
    "make_zoom(map_kappa1_hires, cmap=cm.plasma, clim=clim_hires, rot = rot, dtheta=dtheta, title=r'zoom $10 \\times 10$ deg'+'\\nsmoothing 5 arcmin')\n",
    "filename_plot = os.path.join(pathfig,'figure_kappa_diff.part2.svg')\n",
    "pl.savefig(filename_plot, bbox_inches='tight', transparent=True); print('saved {}'.format(filename_plot))\n",
    "\n",
    "pl.figure(figsize=(nx*fs, ny*fs))\n",
    "make_ball(map_kappa2_lores, cmap=pl.cm.plasma, clim=clim, rot=rot, dtheta=dtheta, title=r'Model 2: $\\Omega_m=0.26$ $\\sigma_8=0.91$' + '\\nsmoothing 1 deg')\n",
    "filename_plot = os.path.join(pathfig,'figure_kappa_diff.part3.svg')\n",
    "pl.savefig(filename_plot, bbox_inches='tight', transparent=True); print('saved {}'.format(filename_plot))\n",
    "\n",
    "pl.figure(figsize=(nx*fs, ny*fs))\n",
    "make_zoom(map_kappa2_hires, cmap=cm.plasma, clim=clim_hires, rot = rot, dtheta=dtheta, title=r'zoom $10 \\times 10$ deg'+'\\nsmoothing 5 arcmin')\n",
    "filename_plot = os.path.join(pathfig,'figure_kappa_diff.part4.svg')\n",
    "pl.savefig(filename_plot, bbox_inches='tight', transparent=True); print('saved {}'.format(filename_plot))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import svgutils.compose as sc\n",
    "from IPython.display import SVG # /!\\ note the 'SVG' function also in svgutils.compose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lw=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svg_fig = sc.Figure(\"20cm\", \"10cm\",  \n",
    "                    sc.SVG(os.path.join(pathfig,\"figure_kappa_diff.part1.svg\")).scale(1).move(0,   0),\n",
    "                    sc.SVG(os.path.join(pathfig,\"figure_kappa_diff.part2.svg\")).scale(1).move(350, 0),\n",
    "                    sc.Line( ((191,116), (390, 40)), width=lw, color='black'),\n",
    "                    sc.Line( ((206,144), (390, 300)), width=lw, color='black')\n",
    "                   )\n",
    "svg_fig.save(os.path.join(pathfig,\"figure_kappa_diff_model1.svg\"))\n",
    "SVG(os.path.join(pathfig,'figure_kappa_diff_model1.svg'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svg_fig = sc.Figure(\"20cm\", \"10cm\",  \n",
    "                    sc.SVG(os.path.join(pathfig,\"figure_kappa_diff.part3.svg\")).scale(1).move(0,   0),\n",
    "                    sc.SVG(os.path.join(pathfig,\"figure_kappa_diff.part4.svg\")).scale(1).move(350, 0),\n",
    "                    sc.Line( ((191,116), (390, 40)), width=lw, color='black'),\n",
    "                    sc.Line( ((206,144), (390, 300)), width=lw, color='black')\n",
    "                   )\n",
    "svg_fig.save(os.path.join(pathfig,\"figure_kappa_diff_model2.svg\"))\n",
    "SVG(os.path.join(pathfig,'figure_kappa_diff_model2.svg'))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}