csiborgtools/notebooks/overlap/lagrangian_patch_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lagrangian patch example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import scienceplots\n",
    "\n",
    "import csiborgtools\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
    "bounds = {\"dist\": (0, 155)}\n",
    "\n",
    "nsim0 = 7444 + 24 * 61\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": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15.319891\n",
      "74.15419337073868\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": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overlaps are  [7.8512335e-01 9.1506150e-03 2.9659141e-03 3.0178351e-03 2.0852157e-04]\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": 52,
   "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": 55,
   "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": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_2722989/1176586521.py:10: RuntimeWarning: Mean of empty slice\n",
      "  axs[0].imshow(np.nanmean(delta1_initial, axis=-1), **kwargs)\n",
      "/tmp/ipykernel_2722989/1176586521.py:11: RuntimeWarning: Mean of empty slice\n",
      "  axs[1].imshow(np.nanmean(delta2_initial, axis=-1), **kwargs)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzMAAAGhCAYAAABVmurXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9zc4kR3IuDD5m7pGZVWRLRWIWBxBm8HVpeTaDYmsxGAwwwCkuP61INXQBIm/goAldgcC+A/JcgNAqrvQtWVoOZtEiZzZnyWpgAGFW0yx9UpP1Zoa7zcLM3M0jI9+3qlhk/TAMLL6ZkZERHpGZbv6YPfYYiYhgs80222yzzTbbbLPNNtvsNTN+2QPYbLPNNttss80222yzzTZ7HtvAzGabbbbZZpttttlmm232WtoGZjbbbLPNNttss80222yz19I2MLPZZpttttlmm2222WabvZa2gZnNNttss80222yzzTbb7LW0Dcxsttlmm2222WabbbbZZq+lbWBms80222yzzTbbbLPNNnstbQMzm2222WabbbbZZpttttlraRuY2eyNss8//xxEhA8//PDiPp988gmICJ988snZa48fP8Ynn3yC9957D0SEd955Bx9++CEePXr0kxx/s80222yz57fNB2y22c/PNjCz2Rtnd+7cwRdffHHx9c8//xx37tw52/7o0SO89957+Prrr/Hpp5/i22+/xb/8y7/g3XffxePHj3+y47/K9sknn+Cdd97BO++8g48//vhs3L/97W/xl3/5l3jnnXdWHXl8/eOPP37m42+22Wab3WSbD/jx7KY5/mnn8Pfffx9EdPE8N72+2WbRNjCz2Rtn7777Lu7du7fqbB4+fIh3330Xd+/ePXvt/fffx7179/Dll1/i/v37uHPnDu7du4fPPvsM9+7d+8mO/yz28OFD/OVf/uVzvfdZzaOHf/jDH/Dtt9+2bW6//e1v8bvf/Q5ffvkl/vCHP+Dhw4f47W9/217//PPP8dlnn+HLL7/EV199hX/9138dnOFNx99ss802exrbfMCPYzfN8U87h3/xxRf44x//ePE8N72+2WZnJptt9gbZZ599Jnfv3pVPP/1U7t+/f/b6Rx99JL/5zW/k3r178pvf/GZ4HwD59ttvX+rxn9W+/PJLuXv37gs95iVbG79v+/bbbwWAfPnll8PY7ty5057fuXNneP2rr76SOAVdd/zNNttss6exzQf8OPY0c/zTzuF3796VBw8eyKUl6E2vb7bZ0rbMzGZvpH3wwQd4+PDhWYr7888/X6U3xUjZq3D8pXnqnojw3nvv4eHDh/jwww/x/vvv49GjRyAiENEwno8//hjvvPMO/vIv/xKff/75sP2TTz5pr7/zzjvXUiaiXYqWOd/7/v37bdv9+/fx+PFjfP3113j06BEeP348vO6RyIcPH954/M0222yzZ7HNB7xYH3DTHO920xz+29/+Fnfv3r2Yibrp9c02W7MNzGz2RppPhv/0T//Utn3xxRe4e/fuavr/66+/Xt3+so4f7eHDh/jiiy/whz/8ASKCTz/9FO+++y4ePHiABw8e4O7duxARiEhzlDHd/+WXX+KTTz5pDufRo0f4/PPP8eGHH+IPf/gD/uZv/uapClA//fTTxoF+/PgxPv74Y3z00Ue4c+fOtQ7s0aNHF4999+7d9tp1x99ss802exbbfMCL9QE3zfHAzXO4ix989tlnq8e56fXNNrtkG5jZ7I21jz/+eJgUf/e7361GzF7V47s9fvwYf/zjH5tDuH///rVRq0ePHuGLL77AgwcPcOfOHdy9exeffvopfve737V9PvrooxYl/Oyzz3D37t0bHchvfvMbvP/++y2SB6C951e/+hUADBE6d3BPm2257vibbbbZZs9qmw94cT7gaeb4m+bwv/u7v8NHH310EdTd9Ppmm12yDcxs9sba3/zN3zSKE6BRsw8++GB133v37j2zNOaPfXy3+/fv49133wUR4f3337+RDuDO5pe//GVzKjEqd+kcN43v448/xu9//3t8++23Z8Wdd+7cwf379/HJJ5/g8ePHePToUXPqNzmmd99998bjb7bZZps9q20+4MX5gKeZ46+bwx8+fIiHDx/i008/XT3+Ta9vttl1ll/2ADbb7Mcyn3y/+OKLphpzaWH9/vvvt/T409Kafuzjx/N88803+Pzzz/Hll1/iww8/xKefforf/OY3F99z7949fPXVV890nuvs66+/xueff45vv/22jf+zzz4DETUH/uDBA3z44Yf45S9/iV/96lf45JNP8PDhw+GeLK//0aNHuHPnzlMdf7PNNtvsWWzzAS/OBwC4do6/aQ7//e9/j8ePH7eMjRsR4YMPPsDdu3evff3Bgwcv9Fo2e7Nsy8xs9kbbxx9/jN/97nd48ODBtel/T23/3d/93St1/OUxHjx4gM8++2ygCyzt3r17+Prrr5+pb8HDhw/xV3/1V889NkAd7pdffolvv/0WX375Jb7++uvm3O/evYs7d+4Mxf4eJYwFpZttttlmL9I2H/B09jQ+4Lo5/ib79NNPW12PiDSgJSJ48ODBja9vttm19jIk1Dbb7Mcyl82MBuBMHnIpmyki8s0338idO3fkgw8+kK+++qpt+81vfiOffvrpT3L8NXvw4IF8+umnTQL5gw8+kA8++EBEurzxt99+K19++aV88803IqLyoPfv32/P/RgiIvfv35c7d+7IV1991Y63HP+a3bt3Tz744AP55ptv5JtvvpGPPvpouBdfffVVu64HDx60c7j5vfvmm2/k22+/PbtHNx1/s8022+wm23zAj+cDbprjn2UOX0rzP+vrm20WbcvMbPbG26effopPP/30xtT+3bt38Yc//AHvvvsuPvzwwyaB+fjx42tpTj/28e/evYsvv/yy8Z8fP36M//E//gcAjcDdu3cPv/zlLweusTdhe++99/DOO+/gs88+O5PU/Id/+Af88pe/xKNHj/DVV1/dOH7vVP3ee+/hvffewx//+MeBxvDo0aN2XZ999hm++uqroUj1o48+wscff4z3338f7733Hu7fvz+M+abjb7bZZps9j20+4MX4gJvm+G0O3+xlGYmIvOxBbLbZZj+deRfqrdBys8022+znZ5sP2OxNsy0zs9lmm2222WabbbbZZpu9lraBmc0222yzzTbbbLPNNtvstbQNzGy22WabbbbZZpttttlmr6VtNTObbbbZZpttttlmm2222WtpW2Zms80222yzzTbbbLPNNnstbQMzm2222WabbbbZZpttttlraRuY2WyzzTbbbLPNNttss81eS8svewA/pv3X//pf8Zd/+Zf4t3/7N/zFX/zFDz7eizjOq3KMV2ks2/W8mOP8P/+3fz3bdoXvscetHzSOF3GMV2ksL/J6/u//6//tBx3jdf6+/VTH+Oabb/A//+f/fK5jbT7g9RjLdj0/7lj+P1//f3/QMX5qH/B/+V9/dfG1N/HzeRWO8SqN5Xl8wBsNZg6HQ3v8z//8zz/4eH/913/9g4/zqhzjVRrLdj0v5jjv84dn2/7f8v/A/5n+rz9oHC/iGK/SWF7k9bwK35U39ffzj//4j/jHf/zHYR5/Vtt8wOsxlu16ftyxvPv/uvuDjvFT+4B//ucHF197Ez+fV+EYr9JYnscHvNE0s7/4i7/AP//zP+O///f//kKO97d/+7cv5Dg/1F7UOF7EcV6VY7woe5Wu50Uc57/g//gCRvJi7EWM5VU5xos8zg+1V+n79iKP8bd/+7f453/+5x8U5dt8wI9/nFflGC/KXqXreZXG8iLsRcyZr9I9eVXu7at0PS/LB7zR0swvCq2+SHsVx/Sm2Jt8b9eyLj+Vvaio3M/NvqyXo4tub/J39kXZD7lHr+L9fRXH9KbYm3xvNx9w2Z5mrn1V7U3+zr4oe5p79EbTzF5Fe1WQ/Jto2719RiN6qt3+C/5Pui+9eolc4qe7hjWTKssN4fFPE+PZvrM/P9s+8x/Ptnv7gmzhG57JB8R5tG17Y2PmP9i27+yLsQ3M/MS2fXF/PNvu7VPaGoi5xkn9F/pffhBoeJZzrTrCZ3n
      "text/plain": [
       "<Figure size 830x498 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context(\"science\"):\n",
    "    # plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "    fig, axs = plt.subplots(1, 2, figsize=(8.3, 0.6* 8.3), sharey=True, sharex=True)\n",
    "    fig.subplots_adjust(wspace=0.05, hspace=0.05)\n",
    "\n",
    "    kwargs = {\"origin\": \"lower\", \"cmap\": \"viridis\", \"vmin\": 0,\n",
    "              \"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",
    "    for i in range(2):\n",
    "        axs[i].set_xlabel(r\"$x ~ [\\mathrm{cMpc} / h]$\")\n",
    "    axs[0].set_ylabel(r\"$y ~ [\\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
}