mirror of
https://github.com/Richard-Sti/csiborgtools.git
synced 2024-12-23 06:08:01 +00:00
213 lines
193 KiB
Text
213 lines
193 KiB
Text
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# Lagrangian patch example"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 1,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"import csiborgtools\n",
|
||
|
|
"\n",
|
||
|
|
"import scienceplots\n",
|
||
|
|
"\n",
|
||
|
|
"%load_ext autoreload\n",
|
||
|
|
"%autoreload 2"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"FNAME overlap_09724_09844_13.25.npz\n",
|
||
|
|
"FNAME overlap_09844_09724_13.25.npz\n",
|
||
|
|
"FNAME overlap_09724_09844_13.25.npz\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
|
||
|
|
"bounds = {\"dist\": (0, 155)}\n",
|
||
|
|
"\n",
|
||
|
|
"nsim0 = 9844 - 24 * 5\n",
|
||
|
|
"nsnap_final = paths.get_snapshots(nsim0, \"csiborg1\")[-1]\n",
|
||
|
|
"snap0_initial = csiborgtools.read.CSiBORG1Snapshot(nsim0, nsnap=1)\n",
|
||
|
|
"cat0 = csiborgtools.read.CSiBORG1Catalogue(nsim0, paths, snapshot=snap0_initial, bounds=bounds)\n",
|
||
|
|
"\n",
|
||
|
|
"nsimx = 9844\n",
|
||
|
|
"nsnap_final = paths.get_snapshots(nsimx, \"csiborg1\")[-1]\n",
|
||
|
|
"snapx_initial = csiborgtools.read.CSiBORG1Snapshot(nsimx, nsnap=1)\n",
|
||
|
|
"catx = csiborgtools.read.CSiBORG1Catalogue(nsimx, paths, snapshot=snapx_initial, bounds=bounds)\n",
|
||
|
|
"\n",
|
||
|
|
"pair_overlapper = csiborgtools.summary.PairOverlap(cat0, catx, 13.25, )\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 3,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"15.431657\n",
|
||
|
|
"95.95069809055808\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"k = np.argsort(cat0[\"totmass\"])[::-1][1]\n",
|
||
|
|
"\n",
|
||
|
|
"print(np.log10(cat0[\"totmass\"][k]))\n",
|
||
|
|
"print(cat0[\"dist\"][k])\n",
|
||
|
|
"\n",
|
||
|
|
"hid0 = cat0[\"index\"][k]\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 4,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Overlaps are [7.2361535e-01 2.8234517e-02 2.3738865e-02 3.6706128e-03 3.3552472e-03\n",
|
||
|
|
" 1.9722811e-03 2.3486929e-03 1.5225195e-03 9.7722618e-04 6.2835851e-04\n",
|
||
|
|
" 1.6445392e-04 4.0640763e-04 3.2565400e-05]\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"print(\"Overlaps are \", pair_overlapper.overlap(True)[k])\n",
|
||
|
|
"hidx = catx[\"index\"][pair_overlapper[\"match_indxs\"][k][0]]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 5,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": []
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 7,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"overlapper = csiborgtools.match.ParticleOverlap(**{\"box_size\": 2048, \"bckg_halfsize\": 512})\n",
|
||
|
|
"\n",
|
||
|
|
"pos0_initial = snap0_initial.halo_coordinates(hid0)\n",
|
||
|
|
"pos0_initial = csiborgtools.match.pos2cell(pos0_initial / cat0.boxsize, 2048)\n",
|
||
|
|
"mass0_initial = snap0_initial.halo_masses(hid0)\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"posx_initial = snapx_initial.halo_coordinates(hidx)\n",
|
||
|
|
"posx_initial = csiborgtools.match.pos2cell(posx_initial / catx.boxsize, 2048)\n",
|
||
|
|
"massx_initial = snapx_initial.halo_masses(hidx)\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 9,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"smooth_kwargs = {\"sigma\": 1, \"mode\": \"constant\", \"cval\": 0}\n",
|
||
|
|
"# smooth_kwargs = None\n",
|
||
|
|
"\n",
|
||
|
|
"delta1_initial, delta2_initial, cellmins, __ = overlapper.make_deltas(\n",
|
||
|
|
" pos0_initial, posx_initial, mass0_initial, massx_initial, smooth_kwargs=smooth_kwargs, )\n",
|
||
|
|
"\n",
|
||
|
|
"x0 = cellmins[0] / 2048 * 677.7\n",
|
||
|
|
"xf = (cellmins[0] + delta1_initial.shape[0]) / 2048 * 677.7\n",
|
||
|
|
"\n",
|
||
|
|
"y0 = cellmins[1] / 2048 * 677.7\n",
|
||
|
|
"yf = (cellmins[1] + delta1_initial.shape[1]) / 2048 * 677.7"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 10,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAxgAAAGeCAYAAAAXC6eIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9y64lx3UujH4jInOtouS9/yKBc0AfwzxwqXdaRkluuWmyd5qk9hOIegMRfgKh+AaknkCb1fq7LDfdklRwZzdVBGwYUIsu4MASa82MGKcxLjEiMudaVWSxbowBrDXnzJmXyJyZY8Q3vnEhZmZMmTJlypQpU6ZMmTJlynOQ9LIHMGXKlClTpkyZMmXKlDdHJsCYMmXKlClTpkyZMmXKc5MJMKZMmTJlypQpU6ZMmfLcZAKMKVOmTJkyZcqUKVOmPDeZAGPKlClTpkyZMmXKlCnPTSbAmDJlypQpU6ZMmTJlynOTCTCmTJkyZcqUKVOmTJny3GQCjClTpkyZMmXKlClTpjw3mQBjyispn3/+OYgIH3300dl1PvnkExARPvnkk913jx8/xieffIKf/vSnICK8/fbb+Oijj/Do0aMXsv8pU6ZMmfLtZdqAKVNeb5kAY8orK7dv38b9+/fPfv/555/j9u3bu+WPHj3CT3/6Uzx8+BD37t3Df/3Xf+Ff/uVf8M477+Dx48cvbP+vsnz66af4yU9+grfffrsznvfv3wcRHf59/vnnAMSwfvTRR3j77bfxk5/85ND4mnzwwQcgou/9fKZMmfLmybQB35+cswEmn3zyCd5++228/fbb+OUvf3n2vG7S8dMG/HBlAowpr6y88847uHv37qEBePDgAd555x3cuXNn990HH3yAu3fv4ssvv8T777+P27dv4+7du/jss89w9+7dF7b/Z5EHDx7gJz/5ybfa9lnl008/xW9/+1t8+eWX+Oqrr/DgwQN8+umnAIAPP/wQf/zjH7u/P/zhD7h9+zbef/99AMAvfvELfPDBB/jqq6/w2Wef4f79+4dewPv37+Prr79+Iec0ZcqUN0+mDfh+5DobAMCZmK+++gr/9V//5ctGuUnHTxvww5YJMKa80vK//tf/wmeffbZb/sUXX+DDDz/cLf/888/x6NEj/OY3v3kl9v+qidH69+7dw507d3D79m3cu3cPv/71r32dO3fudH9x/UePHuHx48f4+OOPHXQYyBg9XJ988gn++Z//+QWf4ZQpU94kmTbg+crT2ID79+/jN7/5jbM3n332GR48ePDMOn7agB+2TIAx5ZWWDz/88FCxff755/jlL3+5Wz96lF6F/Y9itDMR4ac//SkePHiAjz76CB988AEePXrk4UhxPL/85S89HMnClGz5J5984t+//fbb19L9ADw+2NgIe//48WM8fPhwt755oD7++GMAAj5GY/yzn/0MAPD73//el3366ae4c+fOt/bmTZkyZQowbQDwcmzATczDTTp+2oApE2BMeaXFFNT//t//25fdv3/fveujPHz48HD5y9p/lAcPHuD+/fv46quvwMy4d+8e3nnnHXzxxRf44osvcOfOHTAzmNmNV6Sqv/zyS3zyySduBB49eoTPP/8cH330Eb766iv8/Oc/vzHJ8DqjcbTdL37xC9y7d69bNp6/AQsDGuYhO/IKTpkyZcqzyLQBL94G3Lt3z/MuHj9+jF/+8pfOWgM36/hpA6YAE2BMeQ3kl7/8Zaeofvvb3x56ll7V/Zs8fvwYX3/9tSvp999//1rvzqNHj3D//n188cUXuH37Nu7cuYN79+7ht7/9ra/z8ccfuzfts88+O2QYohgIiJ4qMyqj4bl//z7eeeedztN1JPfu3cOvfvUrP69f/OIX+Pjjj7+1EZ4yZcqUKNMGvFgb8Ktf/QoffPCBsyIAun3epOOnDZgCTIAx5TWQn//853j48KErwfv37x/GxgLA3bt3n7lM4Pe9f5P3338f77zzDogIH3zwwY1UthmAv/u7v3NFH71X545x3fgsb+KTTz7B48eP8ejRIzekozH47LPPzl4Hk48++gh37951luPBgwd48ODBjvWYMmXKlG8r0wa8WBvwy1/+Er/73e/wX//1X7sk75t0/LQBU0wmwJjyyospxPv37+Pzzz/H3bt3z3pGPvjgg8N42pe5/3icP/7xj/jss89w+/ZtfPTRR13ljiO5e/euK3n7+/LLL5/52FG++OILAGK0LIYX6AHG48eP8eDBA3zwwQdn9/PRRx+5R83kyy+/xOPHjz3G2Kqi3FRvfsqUKVPOybQBL84GPHz4EJ9//rkneRszcv/+fdy/f/9GHT9twBSTCTCmvBbyy1/+Er/97W/xxRdfXEtdGy37i1/84pXa/7iPL774Ap999llHdY9y9+5dPHz48JkM2YMHD/AP//AP165z+/ZtfPnll26oHj58uDOo5gEzOn0US0ocvVT37t3zGGJmxh/+8AcAADO7UZsyZcqUZ5VpA55OnpcNOCc36fhpA6aYTIAx5bWQDz/8EA8fPsSDBw/w85///Np1v/zyS6/MEZPhPvnkk7Peou97/4DQ7p9++qknzn355Zeu0GMJ2AcPHuDRo0e4c+cOPv744y5pz/Zh8vnnn7sBsvWs4tM5efjwoY/7/v37+PWvf70ruWjHO6qU8tFHH+Ef/uEf8POf/9zP5XVpLjVlypTXU6YNeDE24O7du7h79y5+8Ytf4NGjRx5CdefOnRtDZqdMiTIBxpTXRu7du4d79+7dWB7wzp07+Oqrr/DOO+/go48+8nKAjx8/vlZBft/7v3PnDr788kuPp338+PFOqf/d3/1dxwpY46af/vSnePvtt/HZZ5/tygv++te/xt/93d/h0aNH3hTvOnn06JGP+7PPPsMf/vCHXaLho0ePznawvX//ftfl1f4ePHhw7XGnTJky5bvItAEvxgZYV/Kf/vSn+OlPf4qvv/7amYgpU55WiJn5ZQ9iypQpzy7WTXYm002ZMmXKD0+mDZjyKstkMKZMmTJlypQpU6ZMmfLcZAKMKVOmTJkyZcqUKVOmPDeZAGPKlClTpkyZMmXKlCnPTWYOxpQpU6ZMmTJlypQpU56bTAZjypQpU6ZMmTJlypQpz00mwJgyZcqUKVOmTJkyZcpzkwkwpkyZMmXKlClTpkyZ8txkedkD+D7lr+h/4i381VOv/wR/wSXeeqZjfJtt5rHmsZ7nNvNY81jPc5un2e7/+f95G//n//yfZ97vi5ZpA+axXtY281jzWK/CsV6mDXijAcZb+Cv8Pf3jU6//b/yvz7T+t91mHmse63luM481j/U8t3ma7d76yZNn3ufLkGkD5rFe1jbzWPNYr8KxXqYNmCFS31Hexd++sO1e5LG+rbzq5zWvxXc/1reR1+G8XvXr8aZei9dd5r393Y/1OlzDbyPzWnz3Y73q1/Dbyg/BBrzRDMaLkHfpvRe23Ys81reVV/285rX47sf6NvI6nNerfj3e1Gvxusu8t7/7sV6Ja0gU3g++1xK+f4bK/q/ttfge5HU4r1f9erxu12ICjCDTa9dkXote5vVoMq9Fk3kt3iyZv2eTH9S1MPCgwIISdV+/W//fAXRUefkBtxD7Qd0bN8i8FudlAowg02vXZF6LXub1aDKvRZN5Ld4smb9nkx/EtQjAghIJiEgEoh5g/HX+CVAVUDCBKwMEgOu3P/ZrDFB+EPfGU8q8FudlAowpU6ZMmTJlyg9LFERQzg1YLDolSqkHGcxAKfK2FBCxgossy2oPFowBGZf3EsDJaww2pkw5JxNgTJkyZcqUKVN+eEJJ2IucgJwFbAQmw1erHNgOEkBQK1iBARF36/t2GY35iKLghKsCFdvnlClvkLzRAOMJ/oJ/43/Fu/jbSWNNmTJlyneUP/G/40/4D7z3n3/9sofyVDJtwJSdDMwFrYswFzkDF6swF0RAConetQLMYGbQtsmyUgR4AAFoNJDBETDUPVvBpQgAKUWZjjpBxpRXXp7FBrzRAOMSb32r+r9TpkyZMmUv79J7eBfv4a2/eT36YEwbMOWsGHNBJOAiJ1BKAizSUEVKGQaqFbwswko4k8HYcxe
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 830x415 with 2 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"with plt.style.context(\"science\"):\n",
|
||
|
|
" fig, axs = plt.subplots(1, 2, figsize=(8.3, 8.3 * 0.5), sharey=True, sharex=True)\n",
|
||
|
|
" fig.subplots_adjust(wspace=0.05, hspace=0.05)\n",
|
||
|
|
"\n",
|
||
|
|
" kwargs = {\"origin\": \"lower\", \"cmap\": \"viridis\", \"vmin\": 0, \"extent\": [x0, xf, y0, yf]}\n",
|
||
|
|
"\n",
|
||
|
|
" axs[0].imshow(np.mean(delta1_initial, axis=-1), **kwargs)\n",
|
||
|
|
" axs[1].imshow(np.mean(delta2_initial, axis=-1), **kwargs)\n",
|
||
|
|
"\n",
|
||
|
|
" axs[0].set_ylabel(r\"$y ~ [\\mathrm{cMpc} / h]$\")\n",
|
||
|
|
" for i in range(2):\n",
|
||
|
|
" axs[i].set_xlabel(r\"$x ~ [\\mathrm{cMpc} / h]$\")\n",
|
||
|
|
"\n",
|
||
|
|
" axs[0].set_title(f\"MCMC step {nsim0}\")\n",
|
||
|
|
" axs[1].set_title(f\"MCMC step {nsimx}\")\n",
|
||
|
|
"\n",
|
||
|
|
" fig.tight_layout()\n",
|
||
|
|
" fig.savefig(\"../../plots/lagrangian_patch_example.pdf\", bbox_inches=\"tight\", dpi=450)\n",
|
||
|
|
" fig.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.4"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 2
|
||
|
|
}
|