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": 2,
   "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": 3,
   "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": 4,
   "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": 5,
   "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": 6,
   "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": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzMAAAGhCAYAAABVmurXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOx9d5gcxZ32W9U9s1GbpFWOq5xhFRASSULCGDBOyBhHkiSfsY3DZwt8tnHgDos7Y85ZAgM2tjGWsE2wCbsSIFBAYVFEAWmFct4gbZ7pqu+PCl3d0zObFft9nt3p6VBd3TNTb7+/VIRzzhEiRIgQIUKECBEiRIgQ5xno2e5AiBAhQoQIESJEiBAhQrQFoZgJESJEiBAhQoQIESLEeYlQzIQIESJEiBAhQoQIEeK8RChmQoQIESJEiBAhQoQIcV4iFDMhQoQIESJEiBAhQoQ4LxGKmRAhQoQIESJEiBAhQpyXCMVMiBAhQoQIESJEiBAhzkuEYiZEiBAhQoQIESJEiBDnJUIxE+LCwqZFwM8I8OLs5Pssny/2WT4/cVtDlVj/9ASxz6/yRVtV5Wem/RAhQoQI0WYsWr4IZA7B7N8lH6PnL5kPModg/pLEMbqqrgrzl8zHhJ9MAJlDkH9vPmb/bjbKj5efkfZDhAjReoRiJsSFh7Q8YOeS5Ns3LRL7+FFVDvxpAnC0DLhqAXBPJTB7KZBeADRWnbn2z2Usny8E2K/ygZJ5QpyZWPMw8PhgsT1IzJnbS+a1vv0QIUKEaAZ5mXlYsj75GL3orUXIy8xLWF9+vBwTHpyAsn1lWPDJBaj8v0os/eZSFGQVoKqu6oy1fy7j4VcexuDvDkb+vfmBYm3+kvnIvzcf+ffmY97T85Je16xHZoHMIUnP09z2ECFMhGImxIWH9AKge3Gw4NhbKrbnFiVuWzIL6FEMzC4BBswE0vPE+1kLxeuZar812FsqxMGZwIuzgepy4O49wFcq3XUKax4GdjwL3FIi9tlbKtYpbFoEbFootn9uPXBknVfwNNd+iBAhQrQABVkFKO5fHCg4St8rRUFWAYq6JY7Rs34+C8X9i1HyzRLMHDUTeZl5KB5QjIWfX4jiAe4Y3dnttwal75Vi8HfPDAc8/MrDeHbtsyj5Rgn2PLQHpdtK8fAr7hg/+3ezUX6iHHse2oPK/6vU6/xYsn4JKmorkp6nue0hQvgRipkQFyaG3wpsXJi4fudiYNgties3LRIP0rMeOzfaPxexc4nof3qeeD9rIbCvVHhPGqqAt+YLj1NekdjnqgXAmofc45fPF8fkFYm/6x4D1j7csvZDhAgRohW4ddKtWPhm4hi9eP1i3FKcOEYvWr4I5cfL8dgXWjZGd3b75xqq6qow/7n5WPDJBSgqLEJeZh4WfHIBHnrZHeOXrF+Cx77wmPZKLfz8QpRuK03wzsx/bj7uv+H+pOdqbnuIEH6EYibEhYlhtwQ/CG9aBIwLCG/aWwL0n+k+SJ/t9v1Q4Vc/IyLfZm+p8FosmSVE0s+I+DP7UzJPHPP4YNEvc/3y+e72X+WnDpsz0ZDEWlYt470HzHTXDZgpwueOlokQu8Yq73bljdpb2nz7IUKECNEK3DLhlsAH6UXLF2He1YljdMl7JZg5cmZgeNjZaN8PFb5F5hBM+MkElL5Xitm/m41ZP5+F8uPlIHMIyBzi6c+8p+ch/958DP7uYCxavsizfv6S+Xp7/r35KcPmAOicnpmj3DF85qiZqKqrQtneMr2uOY/Kw688jKJuRSjuH+yJam57iBBBCMVMiAsTeUUyFOxv7rqdS0T4V15ACNjRsuD1Z6t9E3tLRdt37wG+xYXHI70A+Mhi8ZdbJNZ/i7tiSRUVuHuPCOtaPl/0ARDrNy0Chs0W24d/qmVFCK5c4OaxNFSJ5XFzxTlTiZDqclfs+JFb5G5L1X6IECFCtAJFheKB+G/r3DF6yfolKCosQlFh4lhctq8scP3Zat9E6XulWFK2BHse2gP+GMeCTy5AQVYBFn9pMRZ/aTGKCovAH+Pgj3EtllRRgT0P7UHJN0ow/7n5WnSUHy/HorcWYfaE2djz0B58auKnmi1CkEqklJ8Qxy345AKdJ1NVV4V5T8/D3Kvm6j4p787CzwdENbRge4gQyRCKmRAXLsbP84aC7Xg22Gtyrrav0FglxIJ6qB8wM3WOTVW5ED8fWSyOySsSAmjHs+4+4+a6eTuzFgpRsakZApn8HWDALODX+eIPEMcCQI+J4vWoa6HT4qil3pZU7YcIESJEKzHv6nmeULBn1z6LeVd13Bjd2e0rVNVXoaK2QouCmaNmpsyxKT9ejiXrl2DxlxYjLzMPRYVFWPDJBXh2rcsBc6+cq/N2Fn5+IYoKiwLD5hQmDhRjvOmFUeJHCZ3vXP8dzBo1S3t7AHiEyZw/zsHcq+YmFXXNbQ8RIhlCMRPiwsWwTwHHytyH6p1LgvNZACEOWlseubPbV+g/U3hifkaAxbOaDwk7Jsnm8UFuGJnpmQnCgJnJvScKJfOAo2tFFbZ7fAn66Xmin8vnC69KVblbrSyoGIKJ9ILm2w8RIkSIVuJTEz+Fsn1l+qF7yfoluGVC8Bhd3L+41eWRO7t9hZkjZ6IgqwBkDsGsR2Y1GxJWtk+M9YPuH6SFxfzn5uv1yc6hPCxByMvMw8yRMzH/ufmoqqtC+fFyzHtajPGq2MG8p+dh7QdrUfl/lQkFAErfK0XptlIs+OSCwPab2x4iRCrYZ7sDIUJ0GtQD9s4lYrl7cfJQrwGz3BCnloY1dXb75nnu3i1Cw/aWiAf8KxcIT0YydC8GPr++dedJhaNl4vz3VHoT9H9GXBH3kcWib48PEp6ayfNFXpEpZvzXX10uyli3pP0QIUKEaAXUA/iS9UtE5bD+xUmt/rNGzdIhUi3Na+ns9s3z7P7v3Vi0fBFK3ivB7N/NxoJPLsB3rk/OAcX9i7H++x3IAQAWf2kxZv9uNgbdPwgTB0zE/Ovno3RbKYoKi1C2twyLli9C5f9VegoAkDkES9Yvwdo9a1FVV6U9NgpkDsEtE25BUbeilNsXf2lxh15LiAsLoWcmxIWN8fNEeNWOxWI5GcbNFQ/dJXPOrfb9bXxksXjIN0PG/OheLLwzrakCtrcU6DGp7X0DhAiZXSLKKs8uEQJFCby8IiFa9hnJ/spTZBYFCBEiRIgOxLyr5+HZtc9i8brFgYn5Ciq8ac4fWzdGd3b7/jYWf2kxFn5+oSdkzI/i/sUo21fWqrlrSreVYtLA1ByQl5mHkm+WoPL/KlHyzRKU7StLKeBMLLhlgc7r4Y9xrP+eEFr8MY7FX1rc7PYQIVIhFDMhLmwMu0U82O8rFWFhqXBLiVslzEyWXz7fO1/KmWwfEJ6JNQ+7ifF7S1wPkEqgb6gSbVeVi23j5nqT+lUbCpsWiT40VLnzu4ybm7wPPYqFMCmZI9pUYWS5Ra7X5GiZe107l4iyzNcZZUivWiCutUr297U5wKTvtLz9ECFChGglbplwC8r2laF0Wyk+NTH1GF3yjRKUbhNVwsxk+flL5nvmUzmT7QMifO3hVx7WifUl75VoAVHUrQjlx8tRVVeF0vdKUX68HEWFRZh71VxPUr9qQ2HRW4tQtlcIHrXf3KtScABEvozq95L1S/DQyw/pUtPFA4pR3L8Yc/44B+XHy3UYWlFhUdLQuxAhOgqhmAlx4ePKBeKvufCuvCJR3Su9QDzg/4wAf5ogEvBTPVB3dvu5RULAPD5IJMY3VLnz1SgR8PggYI0Ra6wm4vzTBJEzs3FhYtnkNQ+J46rKxSSWzfV/9lLR9z9NEH8NFeI4hepy97o2LhTbzEIF4+aKAglLZonjB8wUAqel7YcIESJEG7Dgkwuw4JMLmg3vKioswp6H9qAgqwCzF84WZZAfnICq+qqUD+Sd3X5RtyKUvFeic2Cq6qoSRMSg+wdhwSvueLrw8wtR3L8YEx6cgPx787HwzYWYOdIoqzxyJh56+SEMun8Qyo+XY/331jfb//IT5brfC99ciPXfW+8pRLD0W0tRkFWACQ9OwIQHJ6CitkJ7WEKE6EwQzjk/250IESLEGcTiWUJ
      "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",
    "    X1 = np.mean(delta1_initial, axis=-1)\n",
    "    X2 = np.mean(delta2_initial, axis=-1)\n",
    "\n",
    "    img1 = axs[0].imshow(X1, **kwargs)\n",
    "    img2 = axs[1].imshow(X2, **kwargs)\n",
    "\n",
    "    kwargs = {\"extent\": [x0, xf, y0, yf]}\n",
    "    # contour_level_X1 = np.percentile(X1, [1, 99])  # This will give the value at the 95th percentile\n",
    "    # # contour_level_X2 = np.percentile(X2, [1, 99])  # You can adjust this to enclose more/less area\n",
    "    contour_level_X1 = 1\n",
    "    contour_level_X2 = 1\n",
    "\n",
    "    axs[0].contour(X1, colors='darkorange', levels=contour_level_X1, **kwargs)\n",
    "    axs[0].contour(X2, colors='darkgreen', levels=contour_level_X2, **kwargs)\n",
    "\n",
    "    axs[1].contour(X1, colors='darkorange', levels=contour_level_X1, **kwargs)\n",
    "    axs[1].contour(X2, colors='darkgreen', levels=contour_level_X2, **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}\", color=\"darkorange\")\n",
    "    axs[1].set_title(f\"MCMC step {nsimx}\", color=\"darkgreen\")\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
}