csiborgtools/notebooks/flow/quijote_Cell.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quijote $C_{\\ell}$ of the LOS pecululiar velocity covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import healpy\n",
    "from h5py import File\n",
    "\n",
    "import csiborgtools\n",
    "import healpy as hp\n",
    "from scipy.special import spherical_jn\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.integrate import quad\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "k, Pk = np.load(\"/mnt/users/rstiskalek/BayesianBulkFlows/data/pk_fiducial.npy\").T\n",
    "loglog_Pk = interp1d(np.log(k), np.log(Pk))\n",
    "\n",
    "def Pk(k):\n",
    "    return np.exp(loglog_Pk(np.log(k)))\n",
    "\n",
    "\n",
    "def djn(n, x):\n",
    "    return (n / x) * spherical_jn(n, x) - spherical_jn(n + 1, x)\n",
    "\n",
    "\n",
    "def C_ell_theory_equal_radius(ell, r):\n",
    "    dDdtau = 35.5\n",
    "    \n",
    "    norm = dDdtau**2 \n",
    "    integrand = lambda k: norm * Pk(k) * djn(ell, k * r)**2\n",
    "    return quad(integrand, 0.0001, 10, limit=1000, limlst=5)[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Num. simulations:   50\n",
      "Num. observers:     8\n",
      "Num. radii:         5\n",
      "Radii:              [100.  112.5 125.  137.5 150. ]\n",
      "Dividing the dipole by a factor of 9:\n"
     ]
    }
   ],
   "source": [
    "num_simulations = None\n",
    "load_skymaps = False\n",
    "\n",
    "with File(\"/mnt/extraspace/rstiskalek/BBF/Quijote_Cell/C_ell_fiducial.h5\", 'r') as f:\n",
    "    num_observers = f.attrs[\"num_observers\"]\n",
    "    num_radii = f.attrs[\"num_radii\"]\n",
    "    npix_per_map = f.attrs[\"npix_per_map\"]\n",
    "    if num_simulations is None:\n",
    "        num_simulations = f.attrs[\"num_simulations\"]\n",
    "\n",
    "    ell = f[\"ell\"][:]\n",
    "    radii = f[\"radii\"][:]\n",
    "\n",
    "    Cell = []\n",
    "    if load_skymaps:\n",
    "        skymap = []\n",
    "\n",
    "    for n in range(num_simulations):\n",
    "        grp_simulation = f[f\"nsim_{n}\"]\n",
    "\n",
    "        for i in range(num_observers):\n",
    "            grp_observer = grp_simulation[f\"observer_{i}\"]\n",
    "\n",
    "            Cell_ = np.full((num_radii, len(ell)), np.nan)\n",
    "            if load_skymaps:\n",
    "                skymap_ = np.full((num_radii, npix_per_map), np.nan)\n",
    "            else:\n",
    "                skymap_ = None\n",
    "\n",
    "            for j in range(num_radii):\n",
    "                Cell_[j] = grp_observer[f\"C_ell_{j}\"][...]\n",
    "                if load_skymaps:\n",
    "                    skymap_[j] = grp_observer[f\"skymap_{j}\"][...]\n",
    "\n",
    "\n",
    "            Cell.append(Cell_)\n",
    "            if load_skymaps:\n",
    "                skymap.append(skymap_)\n",
    "\n",
    "\n",
    "Cell = np.stack(Cell)\n",
    "if load_skymaps:\n",
    "    skymap = np.stack(skymap)\n",
    "\n",
    "\n",
    "print(f\"Num. simulations:   {num_simulations}\")\n",
    "print(f\"Num. observers:     {num_observers}\")\n",
    "print(f\"Num. radii:         {num_radii}\")\n",
    "print(f\"Radii:              {radii}\")\n",
    "\n",
    "\n",
    "\n",
    "print(\"Dividing the dipole by a factor of 9:\")\n",
    "Cell[..., 1] /= 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# healpy.mollview(skymap[2, 4], title=None, cmap=\"coolwarm\", unit=r\"$V_{\\rm rad} ~ [{\\rm km/s}]$\", )\n",
    "# plt.savefig(\"../../plots/skymap.png\", dpi=450)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1113107/3790696556.py:17: IntegrationWarning: The occurrence of roundoff error is detected, which prevents \n",
      "  the requested tolerance from being achieved.  The error may be \n",
      "  underestimated.\n",
      "  return quad(integrand, 0.0001, 10, limit=1000, limlst=5)[0]\n"
     ]
    }
   ],
   "source": [
    "k = -1\n",
    "\n",
    "\n",
    "ytheory = [C_ell_theory_equal_radius(i, radii[k]) for i in ell]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAHKCAYAAAANcD9RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAACDeElEQVR4nO3dd3TV9f3H8ee9N7k3e+9F2HsThqA4kKFFce89quIeP7VWsbUVrbV1kLqtW3GBoy5EQGXIBtlhhRCy97xJ7v3+/giJRFaAO0LyepyTc5r7Xe9LYu6rn2kyDMNAREREpIMze7sAERERkbZAoUhEREQEhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIpN169NFHMZlMFBYWersUj2t67/t64403MJlM7Ny584judSTXHei5x6q19+zIP28RV1EoEmlDNm7cyM0330y3bt0ICAggODiY0aNH88ILL9DQ0ODt8li0aBGPPvoopaWlR3xtW39v0jorVqxg4sSJhISEEBwczPjx41m9erW3yxJxCYUikTbipZdeYtCgQcyaNYtzzjmHZ599lmnTpuHj48Mtt9zCM8884+0SWbRoEX/5y1+OOBS1hfd2xRVXUFNTQ6dOnTxyXXu0cuVKxowZw/bt25k2bRqPPPIIGRkZjB07ls2bN3u7PJFj5uPtAkQE3nvvPW6++WbOPPNMPvjgAwIDA5uP3Xvvvbz33ntERUV5scKj11bem8ViwWKxeOy69ujhhx/G39+fxYsXExkZCcDll19Ojx49+NOf/sQnn3zi5QpFjo1aikS8LDc3l6lTp9KnTx8+/vjjFqGhyaWXXsr48eOP6v6FhYVceOGFhISEEBkZyR133EFtbW3z8auvvprU1NT9rvv9WJZHH32U++67D4DOnTtjMpkOO9bG3e8N4OeffyYtLQ0/Pz+6du3KSy+9dMDzfj826OOPP8ZkMrFgwYL9zn3ppZcwmUysW7fuoGOKWvtcgOzsbK699lpiY2Ox2Wz07duX119//ajfy6Ec6uc9b948TCYTs2bN2u+69957D5PJxOLFiw96759++olx48Y1ByKA+Ph4xo4dy5dffkllZeUR1yvSlqilSMTLnn32WUpLS3n//fex2Wwuv/+FF15Iamoq06dPZ8mSJTz33HOUlJTw1ltvHdF9zj33XLZs2cL777/Pv//97+bWnejo6INe4+739uuvvzJ+/Hiio6N59NFHaWhoYNq0acTGxh722jPPPJOgoCA+/PBDxo4d2+LYzJkz6du3L/369WP58uXH9Ny8vDxGjhyJyWTi1ltvJTo6mq+//prrrruO8vJy7rzzzmN+L/s61M/75JNPJjk5mXfffZdzzjmnxXXvvvsuXbt2ZdSoUQe9t91ux9/ff7/XAwICqKurY926dYwcOfKI6hVpUwwR8aqEhAQjOTnZcDgcLr3vtGnTDMA466yzWrx+yy23GICxZs0awzAM46qrrjI6dep00Ov39dRTTxmAsWPHjlbV4K731mTKlCmGn5+fkZmZ2fzahg0bDIvFsl/t//3vf/er/ZJLLjFiYmKMhoaG5tdycnIMs9ls/PWvfz3odUfy3Ouuu86Ij483CgsLW7x+8cUXG6GhoUZ1dfUR3/NAWvvzfvDBBw2bzWaUlpY2n5Ofn2/4+PgY06ZNO+Qz+vfvb/To0aPFv5fdbjdSUlIMwPj4448PW6dIW6buMxEv2rZtG3v27OHMM8/EbHbPf45Tp05t8f1tt90GwFdffeWW5zVx93tzOBx8++23TJkyhZSUlObXe/fuzYQJE1p1j4suuoj8/Hzmz5/f/NrHH3+M0+nkoosuOubnGobBJ598wuTJkzEMg8LCwuavCRMmUFZWxsqVK13yXpoc7ud95ZVXYrfb+fjjj5vPmTlzJg0NDVx++eWHvPctt9zCli1buO6669iwYQPr1q3jyiuvJCcnB4CampojqlWkrVEoEvGirKwsALp06dLqa6655hrS0tJwOp2tOr979+4tvu/atStms/mI1+s5Ukf63l544QWGDBmCr68vjz766GHPLygooKamZr/3B9CzZ89WPXPixImEhoYyc+bM5tdmzpzJoEGD6NGjxzE/t6CggNLSUl5++WWio6NbfF1zzTUA5Ofnu+S9NDncz7tXr16kpaXx7rvvNp/z7rvvMnLkSLp163bIe99000386U9/4r333qNv377079+fbdu28X//938ABAUFHVGtIm2NQpGIFzWtz1NfX9+q8xctWkRpaSkmk6nV1/ze7xcCPNjCgA6H46ju3+RI31t8fDyPPvoo55133jE990jYbDamTJnCrFmzaGhoIDs7m4ULFx60lehINQXXyy+/nDlz5hzwa/To0S551sEc6Od75ZVXsmDBAnbv3s22bdtYsmTJYVuJmvz9738nLy+Pn376ibVr17Js2bLm93mwIClyvNBAaxEvavp/5uvWrWvV+Z999hlnn302M2bMaPXA5YyMDDp37tz8/datW3E6nc0zzsLDww+47lBmZuZ+rx3Jas1H+t6mTJkCtL5bLzo6Gn9/fzIyMvY7diRr5lx00UW8+eabzJ07l40bN2IYxiFD0ZE8Nzo6muDgYBwOB+PGjTvoPR0Oh0veCxz+5w1w8cUXc/fdd/P+++9TU1ODr6/vEQXB8PBwxowZ0/z9999/T1JSEr169TqiWkXaGrUUiXhRamoqaWlpfPTRR6xcuXK/4xUVFTz99NPN38+fPx+Hw8Hpp5/e6mekp6e3+P75558HYNKkSUBj90pZWRlr165tPicnJ+eA07abptS3ZvHGI31vR8pisTBhwgRmz57Nrl27ml/fuHEj3377bavvM27cOCIiIpg5cyYzZ85k+PDhLULFsTzXYrFw3nnn8cknnxwwHBYUFLj0vcDhf94AUVFRTJo0iXfeeYd3332XiRMnHvVaUTNnzmTZsmXceeedbhsXJ+IxXh7oLdLhrV692ggNDTX8/PyMW265xXjhhReM559/3rjmmmuM0NBQ48wzzzQMwzAaGhqMpKQk47LLLjMyMjIOe9+m2Uj9+/c3Jk+ebKSnpxuXX365ARiXXnpp83mFhYVGYGCg0aVLF+OZZ54xHn/8cSM5OdkYMmTIfrOeli5dagDGGWecYbz11lvG+++/b1RWVh7ze9vXH//4x8POgmqyZs0aw8/Pz0hJSTGeeOIJ429/+5sRGxtrDBgwoFWzz5pcf/31RlBQkGEymYynn376sNcdyXNzc3ONTp06GQEBAcYdd9xhvPTSS8b06dONCy64wAgPDz+qex5Ia3/eTT7++GMDMABj5syZh72/YRjGggULjNNOO8148sknjVdffdW4/vrrDYvFYkycONGor69v1T1E2jKFIpE2YOvWrcY111xjJCYmGr6+vkZsbKwxcuRI49FHHzW2bt1qGIZhZGdnGykpKcbVV1/dqns2fUhu2LDBOP/8843g4GAjPDzcuPXWW42ampoW53733XdGv379DKvVavTs2dN45513Djgl3zAM47HHHjMSExMNs9ncqun5rXlv+zqSUGQYjR/UQ4cONaxWq9GlSxfjxRdfPGDthwpFc+bMMQDDZDIZWVlZrbqutc81DMPIy8szpk6daiQnJxu+vr5GXFyccdpppxkvv/zyUd/z947k520YjVPpw8PDjdDQ0AMeP5CtW7ca48ePN6KiogybzWb06tXLmD59umG321t1vUhbZzIMw/Bgw5SIHKUtW7YwaNAgNm3a1GLadntz0003ERcX16oZaHL0GhoaSEhIYPLkybz22mveLkekTVAHsMhxYuXKlfTo0YPk5GTmzp3r7XJcrqGhgdraWhwOR4v/Le4xe/ZsCgoKuPLKK71dikiboVAkchyoq6vjq6++Yvz48Zx88smUl5d7uySX+9vf/oa/vz+vvvoqf//73/H39+ftt9/2dlntzi+//MIrr7zC3XffzeDBg/fb4kSkI1P3mYhIB3L11VfzzjvvMGjQIN544w369evn7ZJE2gyFIhERERHUfSYiIiICKBS
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "plt.title(r\"$C_{\\ell}$ but $C_{1}$ divided by 9\")\n",
    "ylow, ymed, yhigh = np.percentile(Cell[:, k, :], [16, 50, 84], axis=0)\n",
    "plt.plot(ell, ymed, label=\"Quijote\")\n",
    "plt.fill_between(ell, ylow, yhigh, alpha=0.3)\n",
    "plt.plot(ell, ytheory, label=\"Linear theory\", color=\"black\", linestyle=\"--\")\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(r\"$\\ell$\")\n",
    "plt.ylabel(r\"$C_{\\rm ell}$\")\n",
    "plt.xlim(ell.min(), ell.max())\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.savefig(\"../../plots/C_ell.png\", dpi=450)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "fig, axs = plt.subplots(len(radii), 1, figsize=(5, 4 * len(radii)), sharex=True)\n",
    "\n",
    "for k in range(len(radii)):\n",
    "\n",
    "    ylow, ymed, yhigh = np.percentile(C_ell[:, k], [16, 50, 84], axis=0)\n",
    "    axs[k].plot(ell, ymed, label=rf\"$R = {radii[k]} ~ [\\mathrm{{Mpc}} / h]$\")\n",
    "    axs[k].fill_between(ell, ylow, yhigh, alpha=0.1)\n",
    "\n",
    "    Cell_theory = [C_ell_theory_equal_radius(ell_, radii[k]) for ell_ in ell]\n",
    "    axs[k].plot(ell, Cell_theory, color=\"black\", label=\"Theory\")\n",
    "\n",
    "    axs[k].set_yscale(\"log\")\n",
    "    axs[k].legend()\n",
    "\n",
    "    axs[k].set_ylabel(r\"$C_{\\ell}$\")\n",
    "    axs[k].set_xlabel(r\"$\\ell$\")\n",
    "\n",
    "axs[-1].set_xlim(0, ell.max())\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../../plots/C_ell_fiducial.png\", dpi=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(skymap[4], unit=r\"$V_{\\rm rad} ~ [\\mathrm{km} / \\mathrm{s}]$\",  cmap=\"coolwarm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 4\n",
    "\n",
    "Cell_theory = [C_ell_theory_equal_radius(ell_, radii[k]) for ell_ in ell]\n",
    "\n",
    "plt.figure()\n",
    "plt.title(f\"R = {radii[k]}\")\n",
    "\n",
    "\n",
    "for n in range(50):\n",
    "    plt.plot(ell, C_ell[n, k], alpha=0.1, color=\"black\")\n",
    "\n",
    "\n",
    "plt.plot(ell, C_ell[0, k], color=\"red\", label=\"Simulated\")\n",
    "plt.plot(ell, Cell_theory, color=\"blue\", label=\"Theory\")\n",
    "\n",
    "plt.yscale(\"log\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_csiborg",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}