csiborgtools/notebooks/flow/flow_mapping.ipynb

539 lines
46 KiB
Plaintext
Raw Permalink Normal View History

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using a calibrated flow model to predict $z_{\\rm cosmo}$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The autoreload extension is already loaded. To reload it, use:\n",
" %reload_ext autoreload\n"
]
}
],
"source": [
"# Copyright (C) 2024 Richard Stiskalek\n",
"# This program is free software; you can redistribute it and/or modify it\n",
"# under the terms of the GNU General Public License as published by the\n",
"# Free Software Foundation; either version 3 of the License, or (at your\n",
"# option) any later version.\n",
"#\n",
"# This program is distributed in the hope that it will be useful, but\n",
"# WITHOUT ANY WARRANTY; without even the implied warranty of\n",
"# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General\n",
"# Public License for more details.\n",
"#\n",
"# You should have received a copy of the GNU General Public License along\n",
"# with this program; if not, write to the Free Software Foundation, Inc.,\n",
"# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from h5py import File\n",
"from tqdm import tqdm\n",
"\n",
"import csiborgtools\n",
"\n",
"%load_ext autoreload\n",
"%autoreload 2\n",
"%matplotlib inline\n",
"\n",
"paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def load_calibration(catalogue, simname, nsim, ksmooth):\n",
" fname = f\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/flow_samples_{catalogue}_{simname}_smooth_{ksmooth}.hdf5\" # noqa\n",
" keys = [\"Vext_x\", \"Vext_y\", \"Vext_z\", \"alpha\", \"beta\", \"sigma_v\"]\n",
"\n",
" # SN_keys = ['mag_cal', 'alpha_cal', 'beta_cal']\n",
" # SN_keys = []\n",
" calibration_samples = {}\n",
" with File(fname, 'r') as f:\n",
" for key in keys:\n",
" calibration_samples[key] = f[f\"sim_{nsim}/{key}\"][:][::10]\n",
"\n",
" # for key in SN_keys:\n",
" # calibration_samples[key] = f[f\"sim_{nsim}/{key}\"][:]\n",
"\n",
" return calibration_samples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Test running a model"
]
},
{
"cell_type": "code",
"execution_count": 135,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
" 0%| | 0/19 [00:00<?, ?it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:19: reading the catalogue.\n",
"10:32:19: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/mnt/users/rstiskalek/csiborgtools/csiborgtools/flow/flow_model.py:113: UserWarning: The number of radial steps is even. Skipping the first step at 0.0 because Simpson's rule requires an odd number of steps.\n",
" warn(f\"The number of radial steps is even. Skipping the first \"\n",
" 5%|▌ | 1/19 [00:00<00:05, 3.38it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:19: calculating the radial velocity.\n",
"10:32:20: reading the catalogue.\n",
"10:32:20: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 11%|█ | 2/19 [00:00<00:04, 3.43it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:20: calculating the radial velocity.\n",
"10:32:20: reading the catalogue.\n",
"10:32:20: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 16%|█▌ | 3/19 [00:00<00:04, 3.47it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:20: calculating the radial velocity.\n",
"10:32:20: reading the catalogue.\n",
"10:32:20: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 21%|██ | 4/19 [00:01<00:04, 3.44it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:20: calculating the radial velocity.\n",
"10:32:20: reading the catalogue.\n",
"10:32:20: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 26%|██▋ | 5/19 [00:01<00:04, 3.41it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:21: calculating the radial velocity.\n",
"10:32:21: reading the catalogue.\n",
"10:32:21: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 32%|███▏ | 6/19 [00:01<00:03, 3.44it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:21: calculating the radial velocity.\n",
"10:32:21: reading the catalogue.\n",
"10:32:21: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 37%|███▋ | 7/19 [00:02<00:03, 3.33it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:21: calculating the radial velocity.\n",
"10:32:21: reading the catalogue.\n",
"10:32:21: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 42%|████▏ | 8/19 [00:02<00:03, 3.31it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:22: calculating the radial velocity.\n",
"10:32:22: reading the catalogue.\n",
"10:32:22: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 47%|████▋ | 9/19 [00:02<00:03, 3.14it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:22: calculating the radial velocity.\n",
"10:32:22: reading the catalogue.\n",
"10:32:22: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 53%|█████▎ | 10/19 [00:03<00:02, 3.18it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:22: calculating the radial velocity.\n",
"10:32:22: reading the catalogue.\n",
"10:32:22: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 58%|█████▊ | 11/19 [00:03<00:02, 3.19it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:22: calculating the radial velocity.\n",
"10:32:23: reading the catalogue.\n",
"10:32:23: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 63%|██████▎ | 12/19 [00:03<00:02, 3.25it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:23: calculating the radial velocity.\n",
"10:32:23: reading the catalogue.\n",
"10:32:23: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 68%|██████▊ | 13/19 [00:03<00:01, 3.23it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:23: calculating the radial velocity.\n",
"10:32:23: reading the catalogue.\n",
"10:32:23: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 74%|███████▎ | 14/19 [00:04<00:01, 3.25it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:23: calculating the radial velocity.\n",
"10:32:23: reading the catalogue.\n",
"10:32:23: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 79%|███████▉ | 15/19 [00:04<00:01, 3.21it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:24: calculating the radial velocity.\n",
"10:32:24: reading the catalogue.\n",
"10:32:24: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 84%|████████▍ | 16/19 [00:04<00:00, 3.26it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:24: calculating the radial velocity.\n",
"10:32:24: reading the catalogue.\n",
"10:32:24: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 89%|████████▉ | 17/19 [00:05<00:00, 3.31it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:24: calculating the radial velocity.\n",
"10:32:24: reading the catalogue.\n",
"10:32:24: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" 95%|█████████▍| 18/19 [00:05<00:00, 3.33it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:25: calculating the radial velocity.\n",
"10:32:25: reading the catalogue.\n",
"10:32:25: reading the interpolated field.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 19/19 [00:05<00:00, 3.30it/s]"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"10:32:25: calculating the radial velocity.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"# fpath_data = \"/mnt/extraspace/rstiskalek/catalogs/PV_compilation_Supranta2019.hdf5\"\n",
"fpath_data = \"/mnt/extraspace/rstiskalek/catalogs/PV_mock_CB2_17417_large.hdf5\"\n",
"\n",
"simname = \"csiborg2_main\"\n",
"catalogue = \"CB2_large\"\n",
"\n",
"nsims = paths.get_ics(simname)[:-1]\n",
"ksmooth = 1\n",
"\n",
"loaders = []\n",
"models = []\n",
"zcosmo_mean = None\n",
"zobs = None\n",
"\n",
"for i, nsim in enumerate(tqdm(nsims)):\n",
" loader = csiborgtools.flow.DataLoader(simname, i, catalogue, fpath_data, paths, ksmooth=ksmooth)\n",
" calibration_samples = load_calibration(catalogue, simname, nsim, ksmooth)\n",
" model = csiborgtools.flow.Observed2CosmologicalRedshift(calibration_samples, loader.rdist, loader._Omega_m)\n",
"\n",
" if i == 0:\n",
" zcosmo_mean = loader.cat[\"zcosmo\"]\n",
" zobs = loader.cat[\"zobs\"]\n",
" vrad = loader.cat[\"vrad\"]\n",
"\n",
" loaders.append(loader)\n",
" models.append(model)\n"
]
},
{
"cell_type": "code",
"execution_count": 143,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Stacking: 0%| | 0/19 [00:00<?, ?it/s]"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Stacking: 100%|██████████| 19/19 [00:06<00:00, 3.06it/s]\n"
]
}
],
"source": [
"n = 400\n",
"zcosmo, pzcosmo = csiborgtools.flow.stack_pzosmo_over_realizations(\n",
" n, models, loaders, \"zobs\")"
]
},
{
"cell_type": "code",
"execution_count": 144,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgGUlEQVR4nO3deVhU9f4H8PcMy7AvimxuaIIbCCiLoGkWimWplWZer5rXa2Wr0eotNVsklxRTr5plYWZa3bSuuSWKG6CIoogbpiYKg7iwK8vM+f3hb+aCgLLMzHeYeb+e5zzBmTNnPscT+ua7HZkkSRKIiIiIqMWTiy6AiIiIiHSDwY6IiIjIRDDYEREREZkIBjsiIiIiE8FgR0RERGQiGOyIiIiITASDHREREZGJYLAjIiIiMhGWogtoCdRqNXJycuDo6AiZTCa6HCIiIjIjkiShuLgY3t7ekMvv3SbHYNcAOTk5aN++vegyiIiIyIxlZ2ejXbt29zyGwa4BHB0dAdz5A3VychJcDZHpK60ohffn3gCAnDdzYG9tL7giIoFKSwHvOz8PyMkB7PnzYG6KiorQvn17bR65Fwa7BtB0vzo5OTHYERmARYUFYHPnaycnJwY7Mm8WFv/72smJwc6MNWQ4GCdPEBEREZkIBjsiIiIiE8FgR0RERGQiOMaOiIjIRKlUKlRWVooug+7DysoKFtXHUjYDgx0REZGJkSQJSqUSBQUFokuhBnJxcYGnp2ez18tlsCMiIjIxmlDn7u4OOzs7Lq5vxCRJQllZGa5evQoA8PLyatb5GOyIiIhMiEql0oa61q1biy6HGsDW1hYAcPXqVbi7uzerW5aTJ4iIiEyIZkydnZ2d4EqoMTT3q7ljIhnsiIiITBC7X1sWXd0vBjsiIiIiE8FgR0RERGQiGOyIiIioRbl+/Trc3d1x8eJF0aU0yLPPPovPP//cIJ/FYEdEREQtyqeffooRI0bAx8dHdCkN8sEHH+DTTz9FYWGh3j+LwY6IiIhajLKyMnz99deYPHmy6FIazN/fHw888ADWrl2r989isCMiIiKjoVQqIZPJsHjxYgQHB8PGxgY9e/bE/v37AQBbtmyBQqFA3759a7xvzpw5kMlktba4uDgBV1HbE088gfXr1+v9c4wy2C1btgw+Pj6wsbFBeHg4Dh06VO+xv/zyC0JCQuDi4gJ7e3sEBQXhu+++q3GMJEmYOXMmvLy8YGtri6ioKGRlZen7MohIB959912cPXtWdBlEZCDp6ekAgNWrVyMuLg7p6eno0KEDxo0bB7VajX379qFPnz613vfqq68iNzdXu02ZMgUdO3bEqFGjDHwFdQsLC8OhQ4dQXl6u188xuidPbNiwATExMVixYgXCw8MRFxeH6OhonDlzBu7u7rWOb9WqFd5//31069YN1tbW2Lx5MyZNmgR3d3dER0cDAObNm4cvvvgC8fHx6NSpE2bMmIHo6GicPHkSNjY2hr5EImqEZcuW4ctlX2LatGn44IMP4OTkJLokohZH89gqERr7SLNjx47BysoKv/76q3YM3SeffIKQkBBcuXIFf/31F7y9vWu9z9HREY6OjgCAGTNmYMeOHUhMTES7du10ch3N5e3tjYqKCiiVSnTs2FF/HyQZmbCwMOnll1/Wfq9SqSRvb28pNja2wecIDg6WPvjgA0mSJEmtVkuenp7S/Pnzta8XFBRICoVC+uGHHxp0vsLCQgmAVFhY2OAaiKjpSspLJHwICR9CsrSzlABIACRPT0/p22+/lVQqlegSiQynpESSgDtbScl9D79165Z08uRJ6datW9VOUaL9OTL0VtKAmqt79tlnpTFjxtTYd/bsWQmAdOnSJWnIkCHSSy+9VO/7Z8yYIXXs2FG6cOFCoz5X3zTXcPLkyTpfr+u+aTQmhxhVV2xFRQXS0tIQFRWl3SeXyxEVFYXk5OT7vl+SJCQkJODMmTMYMGAAAODChQtQKpU1zuns7Izw8PB6z1leXo6ioqIaGxEZjkql0n49Z84cbN68Gb6+vlAqlXjuuecQERFxzyEaRNRypaenIygoqMa+5ORkuLm5oW3btnBzc8PNmzfrfO+sWbOwZs0aJCYm1pgxu3r1avTq1QuBgYF46623AABz586Fv78/AgIC8P333wMASkpKMHToUAQEBCAgIADbt2/HxYsXERgYiHHjxsHX1xdTp07Fpk2bEB4eDn9//xpDu+o6p8aNGzcAAG3atGnuH9E9GVVX7LVr16BSqeDh4VFjv4eHB06fPl3v+woLC9G2bVuUl5fDwsIC//73vzF48GAAdwZhas5x9zk1r90tNjYWs2fPbs6lEFEz/Hfzf7VfT5wwEe6u7oiKisLixYvx8ccf49ChQwgPD8ekSZMwZ84ceHp6CqyWyPjZ2dmhpKRE2Gc31K1bt5CVlVXjlzu1Wo24uDhMnDgRcrkcwcHBdc4unTVrFuLj42uFuoyMDCxatAj79u2Di4sLbty4gdTUVPz44484fPgwysrKEBoaikGDBiE5ORmtW7fGtm3bIEkSiouLcePGDZw6dQo//vgjunTpAn9/fzg4OODgwYNYuXIlli5disWLF9d7Tk238YkTJ9CuXTu4ubk1/Q+zAYyqxa6pHB0dkZ6ejtTUVHz66aeIiYlBYmJik883ffp0FBYWarfs7GzdFUtE9yRJEhYuXKj93t7eHgCgUCjwzjvv4OzZs5gwYQIA4JtvvoGfnx8+//xzVFRUCKmXqCWQyWSwt7cXsjVmfF1GRgZkMhnWrl2L5ORknDp1CmPGjEFBQQE++OADAEB0dDQyMzNrtNp98sknWL58OdavXw8bGxsolUoolUqUl5dj9+7dGDNmDFxcXADcGZt/4MABPP3007CxsUGrVq3wyCOPIDU1FQEBAdi7dy/eeecdpKSkaMf0du3aFV27doWFhQW6d++u7QUMCAjQLpJc3zk19u3bhyFDhjTnNjaIUQU7Nzc3WFhYIC8vr8b+vLy8e/5GLpfL0aVLFwQFBeHNN9/EqFGjEBsbCwDa9zXmnAqFAk5OTjU2IjKMvXv34kjakXpf9/LyQnx8PJKTkxESEoLi4mK89dZb6NWrF7Zu3WrASolI19LT09GtWzf861//wtNPP42QkBCoVCrs2bNHG8wCAgLQu3dv/PjjjwDu/DI4f/585OfnIyIiAl5eXtrt+PHjjfp8Pz8/pKeno2fPnoiJicHSpUsB3MkFGnK5XPu9XC6v0bpYn9u3b2PTpk2YMmVKo+ppCqMKdtbW1ujTpw8SEhK0+9RqNRISEhAREdHg86jVau104k6dOsHT07PGOYuKinDw4MFGnZOIDGPevHkNOq5v3744ePAgVq9eDXd3d5w5cwaPPfYYnnjiCZw7d07PVRKRPhw7dgwBAQEYP348cnJyUFpail9++QXt27evcdzMmTOxePFiqNVqyGQyFBYWQpKkWltoaCgefvhhbNiwQfvUhxs3bqB///745ZdfUF5ejps3b2LXrl0ICwtDTk4O7O3tMXHiREybNk279EpD1HdO4E7vQlhYWK219/TBqMbYAUBMTAwmTpyIkJAQhIWFIS4uDqWlpZg0aRIAYMKECWjbtq22RS42NhYhISF44IEHUF5eji1btuC7777D8uXLAdxpfp42bRo++eQT+Pr6apc78fb2xsiRI0VdJhHVISMjA1u2bAGsG3a8XC7HpEmT8NRTT+Gjjz7CF198gc2bN2PHjh1444038P7772uXPyAi45eeno4nnnjivscNGzYMWVlZuHLlSq3Qdzd/f3+8/vrr6NevHywtLTFkyBDMmzcPo0ePRp8+fSCTyTB79mx4eXlh+/bteOutt2BhYQFbW1t8/fXXDa49JCSkznMCgJWVFZYsWdLgczVLU6ft6tOSJUukDh06SNbW1lJYWJiUkpKifW3gwIHSxIkTtd+///77UpcuXSQbGxvJ1dVVioiIkNavX1/jfGq1WpoxY4bk4eEhKRQK6ZFHHpHOnDnT4Hq43AmRYUyYMEECII0cPVK73ElJecOXSjh16pQUHR2tXWbBy8tL+u677yS1Wq3Hqon0TAfLnbQEarVacnR
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"\n",
"# for i in range(len(nsims)):\n",
" # mask = pzcosmo[i] > 1e-5\n",
" # plt.plot(zcosmo[mask], pzcosmo[i][mask], color=\"black\", alpha=0.1)\n",
"\n",
"# mu = np.nanmean(pzcosmo, axis=0)\n",
"mask = pzcosmo > 1e-5\n",
"plt.plot(zcosmo[mask], pzcosmo[mask], color=\"black\", label=r\"$p(z_{\\rm cosmo})$\")\n",
"\n",
"plt.ylim(0)\n",
"plt.axvline(zcosmo_mean[n], color=\"green\", label=r\"$z_{\\rm cosmo}$\")\n",
"plt.axvline(zobs[n], color=\"red\", label=r\"$z_{\\rm CMB}$\")\n",
"\n",
"plt.xlabel(r\"$z$\")\n",
"plt.ylabel(r\"$p(z)$\")\n",
"plt.legend()\n",
"plt.tight_layout()\n",
"# plt.savefig(\"../plots/zcosmo_posterior_mock_example_B.png\", dpi=450)\n",
"plt.show()"
]
},
{
"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.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}