csiborgtools/notebooks/flow/reconstruction_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "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",
    "from os.path import exists\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from corner import corner\n",
    "from getdist import plots\n",
    "from astropy.coordinates import angular_separation\n",
    "import scienceplots\n",
    "from os.path import exists\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "from reconstruction_comparison import *\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
    "fdir = \"/mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogue = \"CF4_TFR_i\"\n",
    "simname = \"Carrick2015\"\n",
    "zcmb_max=0.05\n",
    "sample_beta = None\n",
    "sample_alpha = True\n",
    "\n",
    "fname_bayes = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"bayes\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "fname_mike = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"mike\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "\n",
    "X = []\n",
    "labels = [\"Full posterior\", \"Delta posterior\"]\n",
    "for i, fname in enumerate([fname_bayes, fname_mike]):\n",
    "    samples = get_samples(fname)\n",
    "    if i == 1:\n",
    "        print(samples.keys())\n",
    "\n",
    "    X.append(samples_to_getdist(samples, labels[i]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = [f\"a_{catalogue}\", f\"b_{catalogue}\", f\"c_{catalogue}\", f\"e_mu_{catalogue}\",\n",
    "          \"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"alpha_{catalogue}\"]\n",
    "# params = [\"beta\", f\"a_{catalogue}\", f\"b_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "# params = [\"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"mag_cal_{catalogue}\", f\"alpha_cal_{catalogue}\", f\"beta_cal_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "\n",
    "\n",
    "g = plots.get_subplot_plotter()\n",
    "g.settings.figure_legend_frame = False\n",
    "g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "plt.gcf().suptitle(catalogue_to_pretty(catalogue), y=1.025)\n",
    "plt.gcf().tight_layout()\n",
    "# plt.gcf().savefig(f\"../../plots/method_comparison_{simname}_{catalogue}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# catalogue = [\"LOSS\", \"Foundation\"]\n",
    "catalogue = \"CF4_TFR_i\"\n",
    "simname = \"IndranilVoid_exp\"\n",
    "zcmb_max = 0.05\n",
    "sample_alpha = False\n",
    "\n",
    "fname = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"mike\",\n",
    "    sample_mag_dipole=True,\n",
    "    sample_beta=False,\n",
    "    sample_alpha=sample_alpha, zcmb_max=zcmb_max)\n",
    "\n",
    "\n",
    "samples = get_samples(fname, convert_Vext_to_galactic=True)\n",
    "\n",
    "samples, labels, keys = samples_for_corner(samples)\n",
    "fig = corner(samples, labels=labels, show_titles=True,\n",
    "             title_kwargs={\"fontsize\": 12}, smooth=1)\n",
    "# fig.savefig(\"../../plots/test.png\", dpi=250)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paper plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 0. LOS velocity example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fpath = \"/mnt/extraspace/rstiskalek/catalogs/PV/CF4/CF4_TF-distances.hdf5\"\n",
    "\n",
    "loader_carrick = csiborgtools.flow.DataLoader(\"Carrick2015\", [0], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_lilow = csiborgtools.flow.DataLoader(\"Lilow2024\", [0], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_cb2 = csiborgtools.flow.DataLoader(\"csiborg2_main\", [i for i in range(20)], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_cb2X = csiborgtools.flow.DataLoader(\"csiborg2X\", [i for i in range(20)], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_CF4 = csiborgtools.flow.DataLoader(\"CF4\", [i for i in range(20)], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n",
    "loader_CLONES = csiborgtools.flow.DataLoader(\"CLONES\", [0], \"CF4_TFR_i\", fpath, paths, ksmooth=0, )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angdist = angular_separation(\n",
    "    np.deg2rad(loader_carrick.cat[\"RA\"]), np.deg2rad(loader_carrick.cat[\"DEC\"]),\n",
    "    np.deg2rad(csiborgtools.clusters[\"Virgo\"].spherical_pos[1]),\n",
    "    np.deg2rad(csiborgtools.clusters[\"Virgo\"].spherical_pos[2]))\n",
    "k = np.argmin(angdist)\n",
    "print([loader_carrick.cat[\"RA\"][k], loader_carrick.cat[\"DEC\"][k]])\n",
    "print(csiborgtools.clusters[\"Virgo\"].spherical_pos[1:])\n",
    "print(csiborgtools.clusters[\"Virgo\"].spherical_pos[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loaders = [loader_carrick, loader_lilow, loader_CF4, loader_cb2, loader_cb2X, loader_CLONES]\n",
    "simnames = [\"Carrick2015\", \"Lilow2024\", \"CF4\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\"]\n",
    "\n",
    "\n",
    "with plt.style.context(\"science\"):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    plt.figure()\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "    for i, (simname, loader) in enumerate(zip(simnames, loaders)):\n",
    "        r = loader.rdist\n",
    "        vrad = loader.los_radial_velocity[:, k, :]\n",
    "\n",
    "        if simname == \"Carrick2015\":\n",
    "            vrad *= 0.43\n",
    "\n",
    "        if len(vrad) > 1:\n",
    "            ylow, yhigh = np.percentile(vrad, [16, 84], axis=0)\n",
    "            plt.fill_between(r, ylow, yhigh, alpha=0.66, color=cols[i],\n",
    "                             label=simname_to_pretty(simname))\n",
    "        else:\n",
    "            plt.plot(r, vrad[0], label=simname_to_pretty(simname), c=cols[i])\n",
    "\n",
    "    plt.xlabel(r\"$r ~ [\\mathrm{Mpc} / h]$\")\n",
    "    plt.ylabel(r\"$V_{\\rm rad} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "    plt.xlim(0, 90)\n",
    "    plt.ylim(-1000, 1000)\n",
    "    plt.legend(ncols=2, fontsize=\"small\")\n",
    "    plt.axvline(12.045, zorder=0, c=\"k\", ls=\"--\", alpha=0.75)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/LOS_example.pdf\", dpi=450, bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Evidence comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\", \"CF4\",]\n",
    "catalogues = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "y_BIC = np.full((len(catalogues), len(sims)), np.nan)\n",
    "y_lnZ = np.full_like(y_BIC, np.nan)\n",
    "\n",
    "for i, catalogue in enumerate(catalogues):\n",
    "    for j, simname in enumerate(sims):\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, simname, catalogue, inference_method=\"mike\",\n",
    "            sample_alpha=simname != \"IndranilVoid_exp\",\n",
    "            zcmb_max=zcmb_max)\n",
    "\n",
    "        # y_BIC[i, j] = get_gof(\"BIC\", fname)z\n",
    "        y_lnZ[i, j] = get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "    y_lnZ[i] -= y_lnZ[i].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    figwidth = 8.3\n",
    "    fig, axs = plt.subplots(2, 3, figsize=(figwidth, 0.5 * figwidth))\n",
    "    fig.subplots_adjust(hspace=0)\n",
    "\n",
    "    x = np.arange(len(sims))\n",
    "    y = y_lnZ\n",
    "    for n in range(len(catalogues)):\n",
    "        i, j = n // 3, n % 3\n",
    "        ax = axs[i, j]\n",
    "        ax.text(0.1, 0.875, catalogue_to_pretty(catalogues[n]),\n",
    "            transform=ax.transAxes, #fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='left',\n",
    "            bbox=dict(facecolor='white', alpha=0.5),\n",
    "            )\n",
    "        ax.scatter(x, y[n], c=\"k\", s=7.5)\n",
    "\n",
    "        y_min, y_max = ax.get_ylim()\n",
    "        y_offset = (y_max - y_min) * 0.075  # Adjust the fraction (0.05) as needed\n",
    "\n",
    "        for k, txt in enumerate(y[n]):\n",
    "                ax.text(x[k], y[n, k] + y_offset, f\"({y[n, k]:.1f})\",\n",
    "                        ha='center', fontsize=\"small\")\n",
    "\n",
    "        ax.set_ylim(y_min, y_max + 2 * y_offset)\n",
    "\n",
    "    for i in range(3):\n",
    "        axs[1, i].set_xticks(\n",
    "            np.arange(len(sims)),\n",
    "            [simname_to_pretty(sim) for sim in sims], rotation=35)\n",
    "        axs[0, i].set_xticks([], [])\n",
    "\n",
    "    for i in range(2):\n",
    "        for j in range(3):\n",
    "            axs[i, j].set_xlim(-0.75, len(sims) - 0.25)\n",
    "\n",
    "            axs[i, j].tick_params(axis='x', which='major', top=False)\n",
    "            axs[i, j].tick_params(axis='x', which='minor', top=False, length=0)\n",
    "            axs[i, j].tick_params(axis='y', which='minor', length=0)\n",
    "\n",
    "        axs[i, 0].set_ylabel(r\"$-\\Delta \\ln \\mathcal{Z}$\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(f\"../../plots/lnZ_comparison.pdf\", dpi=500, bbox_inches='tight')\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Dependence of the evidence on smoothing scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 15:27:56\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:32:04\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_smooth_2_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:35:10\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_smooth_3_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:35:08\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_mike_smooth_4_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:34:31\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 15:28:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:37:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_smooth_2_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:39:21\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_smooth_3_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:41:15\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_mike_smooth_4_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:43:31\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 10:48:49\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:44:13\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_smooth_2_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:46:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_smooth_3_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:46:57\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_mike_smooth_4_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:48:46\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 12:03:40\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:52:24\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_smooth_2_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:53:39\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_smooth_3_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:58:54\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_mike_smooth_4_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 11:59:21\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 12:05:48\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:07:50\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_smooth_2_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:06:30\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_smooth_3_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:06:41\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_mike_smooth_4_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:14:03\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 11:13:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_smooth_1_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:20:28\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_smooth_2_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:27:13\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_smooth_3_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:32:47\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_mike_smooth_4_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 13/09/2024 12:36:56\n"
     ]
    }
   ],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "ksmooth = [0, 1, 2, 3, 4]\n",
    "scales = [0, 2, 4, 6, 8]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\"]\n",
    "catalogues = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\"]\n",
    "\n",
    "y = np.full((len(sims), len(catalogues), len(ksmooth)), np.nan)\n",
    "for i, simname in enumerate(sims):\n",
    "    for j, catalogue in enumerate(catalogues):\n",
    "        for n, k in enumerate(ksmooth):\n",
    "            fname = paths.flow_validation(\n",
    "                        fdir, simname, catalogue, inference_method=\"mike\",\n",
    "                        sample_alpha=True, smooth=k,\n",
    "                        zcmb_max=zcmb_max)\n",
    "            if not exists(fname):\n",
    "                raise FileNotFoundError(fname)\n",
    "\n",
    "            y[i, j, n] = get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "        y[i, j, :] -= y[i, j, :].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Carrick2015 2MTF 322.1943359375\n",
      "Carrick2015 SFI_gals 414.5947265625\n",
      "Carrick2015 CF4_TFR_i 835.421875\n",
      "csiborg2_main 2MTF 760.97265625\n",
      "csiborg2_main SFI_gals 800.328125\n",
      "csiborg2_main CF4_TFR_i 1914.80859375\n"
     ]
    }
   ],
   "source": [
    "for i, simname in enumerate(sims):\n",
    "    for j, catalogue in enumerate(catalogues):\n",
    "        print(simname, catalogue, y[i, j, -1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAD1CAYAAAA73MxWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABG3ElEQVR4nO2dd1hb59n/PxJ7C4lhg7FBeNvxkMF2dmKDM+0ssOOspkkNHW863iYQt+nK+2sdSNq3420TcNImaRsHQ5I2SbOQHWfbwch7WwIPbLOE2FM6vz8OkpHZGBCI53NdXHAenXN06xi+vp/nuYdCkiQJgUAgEPSL0tUGCAQCwXhBCKZAIBAMECGYAoFAMECEYAoEAsEA8RzJmxcUFGA2mzEajVgsFnJycgDQ6/WYTCbUajUmk4mMjIwhjQsEAsFoMmIepslkwmQykZaWRlZWFmazmezsbACysrJIS0sjJSUFkIV1KOMCgUAwmoyYYFosFvLy8hzHiYmJFBYWotfrUalUjnGdTkdeXt6gxwUCgWC0GbEpuU6no7i42HFsNBrR6XQYDAbUarVjXK1WYzAYBj3eG1OmTEGn012W7WVlZURHR4/b64UNw3O9sGHs2OCqz2A0Gjl06JDjeETXMO1YLBb0ej3FxcVs2rSpx3Oqq6sHNd4bDQ0NTsfr169n/fr1g7rHmjVrePvttwd1zVi6XtgwPNcLG8aODaP1GbZs2cKWLVscx76+vk6vj4pgbtiwgcLCQlQqFRqNBovF0u2cwY73hq+v72U/2MEK7Fi7fjgYDhvEcxgbz3E4cPVzGK1ncKmDtWbNGqfXRzysKDs7m6ysLLRaLSaTCZ1Oh9lsdrxuNpvR6XSDHh9JXP2P6w5/IMNxD3d4DmPhOQ4Hrn4OY+EZwAgLZkFBAUlJSWi1WkAOD0pKSsJkMjnOMRgMrFu3btDjvXG56xzuwlj5BXMl4hnIiOcgMxzPQTFSxTdMJhPx8fFOY1lZWWRkZKDX6zEYDA6vs2u85WDGe2I41joEAoEAuuvJiAmmqxCCKRAIhotL9USkRgoEAvemqRzeuQVe0sjfm8qHfCshmAKBwL3RPwBnC6HVLH/f9vCQbzUqYUUCgUAw6lQfhEPPw1n9xTHJCpVFQ76lEEyBQOB+HN4MO9LAfxIEx0N9qSyWCg8ITxzybd1uSl5WVsaaNWucovXHCwUFBeTm5pKZmUl6ejogRxukp6ejUCicwqu6XqNQKMjNzQXkiILk5GSys7PJzc0lNDSUzMxMcnNzSU1NJTs7G4PBQGpqKqGhoeTm5pKbm0t2djZLlixxFEgRCMYVtUb4MgP2/k4+nnYbrNoKD56Gu7+AKcngq5G/r3y539tt2bKFNWvWUFZW5vyC5GasXr3a1SYMCaPRKGVlZTmOU1JSHMfFxcVSSkqKlJGR0e26nJwcSaVSOR3X1NQ4jgHJaDQ6ju33LCwslLRabbf79fQeAsGYxNouSca3JOntmyTpz0jSi6GStPvXw/oWl+qJ23mYI0VFbTN3Pfsx075TwF3PfkxFbfOw3r+36k520tPTu5W1MxgMJCQkOI2p1Wqn6k6XYk8i6I1LY2cFgjGHrUP+XlEEH9wFbRZY8Tf4Rhks+cmIvvWEXcNsau3g+Pm6AZ//+Ku72W2sxipJfHzwAvf94TOeeyih/wuBmZOD8ffp+1H3Vt3JjlarRavVOrKlQE4TvVQw7TVDe6O313Nzc0lLS2Pt2rV9Xi8QuATJBme3waEX5LCguz+HyOWw7gBo5o+aGRNWMI+fr+Pan38wpGutksSuk1UDvv6zp29mUay6/xM76VrdqSvp6enk5OQ4BPNyMZvN5ObmYrFYMBqNAH16pwLBqNPeJO90H8qB2hOgngfzviMLqEI5qmIJE1gwZ04O5rOnbx7w+V09TA+FgoR4zaA8zMHQtbpTV1JSUtiwYQNwMS9/MJWcLkWtVpOWlobFYum17J5AMOpIEtQZIWS6LIr7fgdRN8CKv8Kkq0GhcJlpE1Yw/X08B+X1vfaDa0nP3YmhxIwuTk1O2nIiQvyG3a5Lqztduua4du1acnNz+12LHAwqlcqxKy8QuIz2Bjj+Tzj4PNQchofOgn8EPGACDx9XWwdMYMEcLBEhfrz1xI0j+h49VXdKS0sDcHiS6enppKamkp+fP6zvPZwCLBAMCkmCz38AR1+GjkaYdjss3wR+YfLrY0QswQ3jMMcrJpOJ1NRUlixZgkKhQKFQYLFYMBgMbNq0iU2bNjnqidq/TCYTmzZtwmKxOOIw7ej1eqemc1132A0GA/n5+ZhMJnJzc3uM7xQIRpSOFjixBazt8hTb1g4LfgAPlMCt/4Zpt8jT8TGGqFYkEAhGj1qjvIFz9K/QUg13bIfokZ25XQ5uX61oPGf6CARuzWffh39OhyMvwqyH4L5jY1Yse8v0ER6mQCAYGRrKZHGMuxPCFkLJv6HVAtPXgufwb5iOBJfqidj0EQgEw4c9wPzg81D6Nnj4ysUvwhZC3B2utu6yEYIpEAiGj+JN8PVToJ4P1/4RZj4A3oOLQx7LCMEUCARDQ5KgfKfsTYYvhoU/ktcmo693eYD5SCEEUyAQDI62ejjxmiyU1fsgOA6irpdfC4qRv9wUIZgCgWBgtDeBlz+U74JPvysHmF/5DMSsGpMxkyOBEEyBQNA7HS1gLJALYCi94c6PYcoKePAUBE5xtXWjjhDMMUTX/utmsxmj0UhWVpYj20ev15OVlQVcrJ+5bt06kpKSyMnJYevWrWRlZbF27VpRdUhwebSYwbAJjv5NDjCfslKuEgSyNzkBxRJwv4rrOp1OWr16tfTaa6+52pRBUVNTI6WkpDiNda2I3l+F9N5eFwgGjLVdkioM8s+tdZL0aqwkff7fklRzzLV2uYDXXntNWr16taTT6ZzG3c7DjI6OHpHA9fLWGh4++Cxf1x1jafAsXp7/BJE+ocN2/57yuQdSQUhUSBdcNg1lctOwI5vlwPKHL4B3EDxgnDBrk5eyfv161q9fz5o1a5zG3U4wB0qTtYWjjWcGfP5jR/7MrtqjWLFRWG3g7r2/4k9zvjega2cHxODv4dvnOTqdDr1eT2pqKunp6SQlJZGRkdHr+aJCuuCysXXAR+vkDBwPXzlmcv63ZbGECSuWfTFhBfNo4xmW7ByY4F2KFRtf1h4e8PXFy/+MLnhG/+cVF5OZmUlqaioWi4WcnBxHeTcQFdIFw0BLNZx4HeZ9G5SeEBDtlgHmI8WEFczZATEUL//zgM/v6mF6oGRZyOxBeZgDQavVOupcGgwGVq5c6VQfU1RIFwyKpnLY9jBUfA0hMyAwBk69IwecT74awhbJYikYMBNWMP09fAfk9dl5c9EvePjgsxTVHScxeOawr2Hq9XpHozOQp+hpaWmOXfOuiArpggGx7WE4WwiSFSp2QfV+WPo0zP4m+IW72rpxiVikGCCRPqG8v+Q3VN1YwPtLfjOsYmknJyfH6dhkMvXa8ExUSBf0SeN5WSQl68UxL39YnCHE8jIQgjmGiI+PJzc3l4KCAnJzc1m3bh0qlarfCumXvn45jdEE4xxJgiN/g9fngqc/KDzkcYUHhCe61jY3QNTDFAjchboS2JEGZ/Uw6xugexI+/xFUFsliufJl8I90tZXjClEPUyBwR9rqIF8HXsFw+wcw9SZ5fPX7rrXLzXC7KbloUSGYUNQcBWurHBKU9Brce/CiWAqGjGhRIRC4E9Y2MDwDxf8PrsySa1EKhh0xJRcIxjvlRfDxo2A5AoszLxbFEIw4QjAFgvFE9UF4c7kcdJ5SJH8XjBpCMAWC8UDVfghbAJr5sGqr3FBMKf58Rxu32/QRCNyKVoscKrR1IZzRy2Px9wixdBEj/tQtFotTgQiDwYBKpUKtVrN7924SEhJQqVTo9XpMJhNqtRqTyeSo1NPbuLtizxOPj49HrVZjNptJSEhwPKu+Cglf+mzS09O7ZQ/Zsd8jOTkZlUpFZmYmaWlpxMfHU1hYSGJiIklJSf0WLr70daPRSHx8vFPRkN4oKCjAbDYP6NwJScnb8Ol35B461z0vVzoXuJaRKsCZn58vZWRkdCv
      "text/plain": [
       "<Figure size 350x262.5 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "    plt.figure()\n",
    "\n",
    "    ls = [\"-\", \"--\", \"-.\", \":\"]\n",
    "    for i, simname in enumerate(sims):\n",
    "        for j, catalogue in enumerate(catalogues):\n",
    "            plt.plot(scales, y[i, j], marker='o', ms=2.5, ls=ls[i],\n",
    "                     label=catalogue_to_pretty(catalogue) if i == 0 else None, c=cols[j],)\n",
    "\n",
    "    plt.xlabel(r\"$R_{\\rm smooth} ~ [\\mathrm{Mpc} / h]$\")\n",
    "    plt.ylabel(r\"$-\\Delta \\ln \\mathcal{Z}$\")\n",
    "    plt.xlim(0)\n",
    "    plt.ylim(0)\n",
    "    plt.legend()\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/smoothing_comparison.pdf\", dpi=450)\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. External flow consistency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:26:54\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:28:43\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:20:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:19:39\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:27:13\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:28:48\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:19:59\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:26:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:54:29\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:10:19\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:27:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:58:15\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:59:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:16:15\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:39:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:01:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:12:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 14:07:26\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:23:47\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:19:35\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:49:16\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:50:20\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:52:12\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:51:03\n"
     ]
    }
   ],
   "source": [
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CF4\", \"CLONES\"]\n",
    "# sims = [\"Carrick2015\", \"Lilow2024\", \"CF4\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "# cats = [[\"LOSS\", \"Foundation\"], \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "# cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_not2MTForSFI_i\"]\n",
    "\n",
    "X = {}\n",
    "\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"bayes\",\n",
    "            sample_alpha=True, zcmb_max=0.05)\n",
    "\n",
    "        if not exists(fname):\n",
    "            raise FileNotFoundError(fname)\n",
    "\n",
    "        with File(fname, 'r') as f:\n",
    "            X[f\"{sim}_{cat}\"] = np.linalg.norm(f[f\"samples/Vext\"][...], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAFsCAYAAADsalOEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAACyN0lEQVR4nOy9d3gb553g/0EHSJAAwd6bumRbouTeEot0qp3Elqw4u9l2sZTs3e3d7saSlb377d1lE0vaJFvSTCnJJhvHjiW5xI6rKMlVkiWSVu8EeydBECRRCeD3xwBgA0CABEiCms/z4IEwM3jnJfTOd77zrRKv1+tFRERERCRhkM73BEREREREokMU3CIiIiIJhii4RURERBIMUXCLiIiIJBii4BYRERFJMETBLSIiIpJgiIJbREREJMEQBbeIiIhIgiEKbhEREZEEQxTcIiIiIgmGKLhFREREEgxRcIuIiIgkGPLZfLm1b4T+YUes5hKUdK2KwozkuJ5DJARDLWDvi+851BmQUhTfc9zgDPYNYR2yx/UcSSlqdBkpcT2HyBgzFtytfSOsf+qP2JzuWM5nChqljLpdX4xYeJvNZp544gnq6+sBKCsrY/fu3VRUVMzo/GazmbS0NAYGBtDr9SGPk0gkzLbQYk1NDdXV1Rw4cGDCuGVlZYHP1dXVVFZWzuo8ETHUAs+vhFFrfM8jT4LHL4nCO04M9g3xzPYXGHWOxvU8cqWcb+7ZErHwjvV1OlMivb4jYe/evezevRuATZs2Bf598OBBnn76acxmM5WVlVRXV8922jMX3P3DDmxON08+tJrCjKRZTyQYrX1W/vm1C/QPOyIS3EajkfXr17N79+6A8Kuvr58g+KJFr9dTV1c36//U6di8eTNmsxmj0ThlX0NDQ1zPHRR7nyC0K74DKcXxOcdQM9R/XzhXhII73AW/fv16zGYzZWVlGI1G9Ho9u3fvnnCjC3VxLVasQ3ZGnaPc/fA6UjO0cTmHpW+Yj179BOuQPSLBHY/rNNg5IhkvVte30Wjk0KFDNDQ0YDabWb9+PeXl5Tz22GO88MIL1NXVAVBeXs7BgwfZtGnTrM43K1MJQGFGEktyUmc7TEzYtm0bW7duZevWrYFtsbiDx0oLqKmpwWAwBB3vwIED1NfXs3nz5picK2akFIN+2XzPAojsgt+9e3fgovAfX1dXFxDmwS6u8etlsZKaoSU9Rz/f0wDid52OZ8eOHezcuTOicaM5d6hr2Gw2s3PnTkC4GWzbto1Dhw6xdevWCU/Qsbo5LSrnZE1NDdu2bQu5f/PmzZSXl1NeXk5NTU1gu//z+vXrMRqNUz5LJJIJ5/CP4dfwJrNt27YJ488Wv4DZsWNHzMZMREJd8KG0pbKyMjZs2BDQzkNdXCJzy3TX6cGDB1m/fj3r16+fcNz463DPnj0T9o2/Zh9//HEOHjzI5s2bWb9+fcTj+m/0/uvbv24ioaKiYoowNxgMU/6u2tramJg6F43g9psYwt3Rtm3bRkNDAwcOHJig2RqNRnbv3s3hw4cDmtn4z+OP27x5M3V1dYFxJguNzZs3o9frY2aH1uv1HDhwgIaGBmpqajh48GBMxk1Eprvgx2M2m9mzZw8mkymggUdycYnEl+mu0/r6ep5++mkOHz4cMGFEorCMv2aff/55KisrOXDgQMBEEcm4Bw8epLKykoaGBhoaGmblF6uurp4wflVVFU888QQ7d+6Midl11qaShYJ/IYSzbfmFaUVFBWazGbPZHPgRq6urJ/ygkz+D8B/72GOPBbZPPs+2bdswGo2BRyOz2czTTz8d2B9s0U5nYx0YGJgw/xdeeGHW9rFEJJIbMwiPyH5HkF6v5/Dhw0GP819cosY9t0x3nb7wwgts2bIlcI3t3LmT0tLSiHwRwa7ZaMatqKhg8+bNGI1Gtm3bRmVl5Yyu4c2bN3PgwIEJx/jXWVVVFXq9ftbmuUUjuIGAxzbUf/KePXs4depU0H2TF9FsbFE1NTVUVlYGnGPjt4eycYuEJ5IbM0y1cW/cuJGdO3dOudkFu7hE5obprtOZMtv/y8rKSurq6gJmln379k1xYE93DW/evDlsdExVVRXV1dWzFtyLxlQCwh1379697N27N7DNaDRSX1/PwYMHOXToEAcOHGDfvn0zGn/Tpk3s378/YNc2Go2Bf+v1+kAo37Zt24LavqPl4MGD7NmzZ8LnSE0Fi5FoQ6nKysrYuXPnlO9Md3GJxJdw1+mWLVuorq4OXD87duzgscceCxzn3z7dk5Jer59wfU43rn8OZWVlbN++nZ07d4ZU8kKxefPmKQ7RgwcPTjBvHjp0KCZm1Flr3K198YvzjXbssrIyGhsbeeKJJwJ3Sb1ez759+6isrGTHjh2sX7+esrKyGd2dy8rK2LdvX8DhUVZWNuXxrKysjB07dvDEE09M8CZPh9+h6XeQ+LXEbdu2UV5eHrDJzUkM93iGmhfM2NXV1VMiQfwXZzAh7DeHVFVVBbYFu7huBCx9wwtm7HDXaUVFBbt37w5cY+Nv1lu3bmXjxo2B6zecrdhvUy4rK+Pw4cNhx/VTX18fWCv+aztS/AJ6fFCC384+3nxXWVkZkycNiXeGWSMLNQFnIZNQppIFmoAzOY57/AXvvygNBgMmkwmz2cy2bdvYvn07QOARePwF77+4FisLNQEnUVko1/CMBTeIKe+LHjHlfVEgprwvPmYluEVERERE5p5F5ZwUERERuREQBbeIiIhIgiEKbhEREZEEQxTcIiIiIgmGKLhFREREEoxZJeC02Hrocw3Gai5ByVDoKNJkTXtceXn5hHRnP5MLpcei4UE0TFf/uaqqCpPJFCiGs6DobIGBOIcDpmVArhgOGE96enoYHIzvdarT6cjKmv46BSHR5YknnsBoNGIwGBZUXfRQ1+v4ksvB0uAnN0A5ePBgoMhUrJonjGfGgrvF1sPKj/4LVk9847iTpCou3f3LiIR3MOaqEUIwpqv/XF9fj8lkmvN5RURnCzy8EuxxTsBRJ8Grl6YV3lVVVdTW1k4ouuXHX143lg0nYtkZZT7p6enhiSeewOGI73WqUqnYt2/ftMLbbDazceNGDhw4ECjitH///rjOLVLCXa8bN24M1HWvqqoK1CMK1gDFnyTW2NgYKB+8d+/emNZ9n7Hg7nMNYvU4+E7p4xSrZyZUp6PZ3sP3G5+nzzU4Y8ENsS/SPpmZFFcHAmUex1cfWzAM9AlC+4nvQF6cOuB0NMO+7wvnikDrNhgMU7qHjK9HMRsmF6+azxt+LBkcHMThcPDVr341Yo04Wnp6evj973/P4ODgtOeora0Fxip1xqJSXrREe73667r718e2bdsCLQSDNUAxGo1s2LAhsHaqqqpiXoVy1rVKitVZLEsuiMVc4kYo84i/FxzAhg0bqK6uZvPmzYGSjjU1NWzevDmg5W3bto3169dHvNCC3TD89Z9rampmXDNlTskrhuKF0QFn06ZNVFdXTxDc1dXVgQtpNgTrmDLfac2xJCsri/z8/PmeRkBgj7/OxhPsmgTBFOqvcX3gwIFA6VUgUL5Xr9dPMFEEM52GI9T1OrnLkr+AVbhxTCZT4BqPRyXEG9Y5GaqwelVVVcBWdeDAgQkdVGpqaqZUFIuUycXVt23btmDseolCVVXVFA07WP++UJ2O0tLS2LNnT+AR2H/xbdu2LWjHlEg6HwXrqjK5g9LBgwcD372RG2H4aWxsBIT/p7S0tMBvEq7ZwfhGCf7vTm5oUl9fz44dOwLbZ3Mzn3y9Ti4kN91T3r59+wLrKVgDj9myqOpxR0OowuqNjY0BgVpbW8vOnTupqalBr9cHXrMtrn7w4EEqKiooKyuLqj2SyFhp3a1btwbsjJPxa3L19fVs3Lgx8MTkbyRcV1fH5s2bOXjwINu3b6e6ujogGIJdYP7OR36bpb8R8XhB4xcyO3bsmCBkjEZjQJjo9XqqqqpuyEYY4/F3dYKxwl8NDQ3TNjvwV+Lcu3dv0IYmL7zwAkDAbGE0GqmpqWHDhg2zboYwXlCPb8ASjMnrZdu2bezZsydQ7CwW3LCCOxR+4ex/zPE7IMa3I4u2QcL
      "text/plain": [
       "<Figure size 350x331.25 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(3.5, 2.65 * 1.25))\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "    for k, cat in enumerate(cats):\n",
    "        i, j = k // 2, k % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            sns.kdeplot(X[f\"{sim}_{cat}\"], fill=True, bw_adjust=0.75, ax=ax,\n",
    "                        label=simname_to_pretty(sim) if i == 0 else None)\n",
    "\n",
    "        ax.text(0.725, 0.85, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='center',\n",
    "            bbox=dict(facecolor='white', alpha=0.5, edgecolor='none'))\n",
    "\n",
    "        ax.set_ylabel(None)\n",
    "        ax.set_yticklabels([])\n",
    "        ax.set_xlim(0)\n",
    "\n",
    "    handles, labels = axs[0, 0].get_legend_handles_labels()\n",
    "    fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 1.1),\n",
    "               ncol=3)\n",
    "\n",
    "    for i in range(2):\n",
    "        axs[-1, i].set_xlabel(r\"$|\\mathbf{V}_{\\rm ext}| ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "        axs[i, 0].set_ylabel(\"Normalised PDF\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(f\"../../plots/Vext_comparison.pdf\", dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. What $\\beta$ is preferred by the data? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CF4\", \"CLONES\"]\n",
    "cats = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "# cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_not2MTForSFI_i\"]\n",
    "\n",
    "X = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"bayes\",\n",
    "            sample_alpha=True, zcmb_max=0.05, sample_beta=True)\n",
    "\n",
    "        if not exists(fname):\n",
    "            raise FileNotFoundError(fname)\n",
    "\n",
    "        with File(fname, 'r') as f:\n",
    "            X[f\"{sim}_{cat}\"] = f[f\"samples/beta\"][...]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "\n",
    "    fig, axs = plt.subplots(3, 2, figsize=(3.5, 2.65 * 1.8))\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "\n",
    "    for k, cat in enumerate(cats):\n",
    "        i, j = k // 2, k % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            sns.kdeplot(X[f\"{sim}_{cat}\"], fill=True, bw_adjust=0.75, ax=ax,\n",
    "                        label=simname_to_pretty(sim) if i == 0 else None)\n",
    "\n",
    "        ax.text(0.1, 0.85, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='left',\n",
    "            bbox=dict(facecolor='white', alpha=0.5, edgecolor='k')\n",
    "            )\n",
    "\n",
    "        ax.axvline(1, c=\"k\", ls=\"--\", alpha=0.75)\n",
    "        ax.set_ylabel(None)\n",
    "        ax.set_yticklabels([])\n",
    "\n",
    "    handles, labels = axs[0, 0].get_legend_handles_labels()\n",
    "    fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 1.075),\n",
    "               ncol=3)\n",
    "\n",
    "    # for i in range(3):\n",
    "    for j in range(2):\n",
    "        axs[-1, j].set_xlabel(r\"$\\beta$\")\n",
    "\n",
    "    for i in range(3):\n",
    "        axs[i, 0].set_ylabel(\"Normalised PDF\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fig.savefig(f\"../../plots/beta_comparison.pdf\", dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Bulk flow in the simulation rest frame "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAD1CAYAAADUHqdoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmAUlEQVR4nO2dd3hb53m3bywCBBcI7iWRoPayRFEe8Y4oO3ZsxbEpKYpjN8tiRpu0qS1aTtOsJjJZt03ar2lIO4njJFYkMrKteEWE97ZEalqbIDW4SRAcIDbO98dLgKS499C5r+tcAA7OeA9w8MM7nuf3KiRJkpCRkZGRGRHK6S6AjIyMzGxCFk0ZGRmZUSCLpoyMjMwokEVTRkZGZhSop7sA42H58uVkZmYO+F5NTQ0pKSmjPqa838TtNxvKKO83u/cb67kqKiq4dOnSqPcDQJrF3H333WN6b6zHlPebueeS97sy9xvruRISEsa0nyRJ0pxtnm/dunVO7zdWprKc8mcyM/YbK7Ph+qb6MwFQSNLsjdPcuHEj+/btm+5iTAtX8rUPhvyZDIz8ufQnMTGR+vr6Me07Z2uac53p+Ied6cifycDIn0t/xtIPGkAWzVmK/EPoj/yZDIz8ufRHFk0ZGRmZKUIWTRkZGZlRIIumjIyMzCiY1aJZU1PDxo0b2bVr13QXRUZGZhawa9cuNm7cSE1NzZiPMatFMyUlhX379skd3d3YbDby8/MpLi6mtLSU4uJiKioqKC4unpDjZ2ZmYrFYBn3fbDaTmZlJRUXFhJzPZrP1eV1RUYHFYsFms2E2m/u9P6l0tEFT3eQuHW0jKorFYiEvLw+FQjHg91FaWopCoZiw7703w90DM52tW7eyb9++cQ0Ezeo0SpkeLBYLGzZsoLy8HIPBEFy/adMm1q1bNyHnKCkpwWQyDfp+Tk4OWVlZ4z5PaWkpBw4cwGw2U15eHly/c+dOSktLMRgM7Nixg5ycnHGfa0R0tEHRT6G1eXLPEx0LeT+AiKghNzOZTOTl5WG1WikqKqKgoKDP+1arFYPBwLZt28ZVHLPZjNFo7POdDncPXAnM6pqmTA+bNm2ioKCgj2AC/X5QAzFQjW2gdRMhiAHMZvOgNdLc3Fzy8vL6rd+wYQOSJNHa2sr27dsnrCzD4uwSgqkLBUPM5Cy6UHEOZ9eIi5WXl0dpaWmfdRUVFWRnZ4/4GEPV1ge6d8Z7D0xp62CSmNU1TXubY7qLMCOw2WxUVFQMeEObTCZyc3MBKC4uxmQyUVZWxoYNG8jJycFsNrNp0yaefPJJysrKWLt2LSaTqd+67OxsNm3aRFFRUbCGFxA+k8mExWLpJ2TR0dHk5OTw5JNP9hPzsRJonk+kgI8YnR7CIibv+M7R3c8mkwmTyYTZbA5+J1ardUDRHOy7z8vLIz8/H6PRGHwvNzcXs9mM1Wpl9+7dHDx4kG3btlFRUTHie6C0tBSr1YrRaAyuH+he27ZtW3BbAKPRGLxfZyqzWjSt9bZpOW+Xy8uZuvYpOdeipEj02qG/pkAf02DNpsD6goICKisrycnJITo6mtbWVnJycsjJycFqtVJQUIDFYiErK2vQdb3PmZ+fH2w+B34AAcxmM08++eSE/wBsNlsfUZ8oMZ6t5OXl9RGxwRjqu7fZbGzbto2srCw2bdpEbm4uOTk5mEwmtmzZEvyDGuk9YLFY2L17NyUlJYC4F/Lz8ykoKOh3X1VUVFBUVERZWRlAUNBn8vc6q0XTaXfj7HKh02un9Lxn6tq58V9fnZJzvfOTz7A63TjkNgFRtFgsAwqnzWbDYDBQWVmJzWYbsCM/Ozsbg8HQpwY30LoApaWlfX5AvfsezWYzu3fv7rPOZrOxc+fO4OuBhH64roTeorxhwwby8/MpKioacp+5Tm5uLg899BBAsMY5UBN4qO8+8B0YjcZRNZ8Huwd2797dpx890EoJfL+976v8/HwMBgNmszm4feBPeqYyq0UTJC6cqmNRVvqUnnVRUiTv/OQzU3au4QjcgIFm0uWYzWZyc3PJz88nLy+PrKwsjMahhXg8mEwmsrOzKS4uDgqdwWDoI4oDDTIMhdlspqCgIFgjAYJNuiudzZs3B5vfgzHUdz9crc5ms2G1WidlAMhms7Fu3bqg+E7Z4N44mNWiqQ5Rc/5E7ZSLpl6rHrb2N9WUlJQM2LQJCFdxcXGwaQs9fYMHDx4c0/lycnKCNZzA8QJkZWWRm5tLZmZmsJk3XgIjxgHKy8vZsmXLuI87mwnUCvPy8ti0aVOwOXw5w333Q9UuA39MAwnrYPfAli1b+rQqDh48yObNmwc8fqBPNdAXOlxX00xgVo+e6yN0nD859iDVuYTJZKK8vJydO3f2idMM3KyBviSz2UxpaSnbt29n586ddHZ2BvuVAjdsoHO/97qKigoOHjwY/GFmZWWxY8cOCgsLKS0tDcZQVlRUBGuUgb7H0cZtms3m4LmLi4v7dDsUFxdTXFxMZmbmjB8wmCwqKirYuXMnO3fuDDZlA4vFYmHnzp3YbLZgnOZg371arQ52pQS6T6xWa3BEfseOHZSUlGA2mzGZTCO6B0wmE1lZWWzZsiV4H5aVlVFUVDTgfZWVlUVeXl6/Y8xkZrWf5i3Xf5oNGZv5p//7MvoI3XQXR2YUjLZ5Pq001cEvdoiwIJ1+cs7h7BKj5/+4E+KSJuccMkHG4zE6q5vnod1CeeFULUvWzex/J5m+zIa+qyA6vQg8b20edVjQqIiOnTxRlpkwZrVoNjY14DTaKXv+TVk0ZSaPiCiRqTOKwPMxodMPmw0kMz527drFrl27xpV7PqtFMyUlhWs/vZZLZxumuygyc52IKFnQ5gBbt25l69atbNy4cczHmNUDQQDpy1Jormmls22SawEyMjIyzHbRlPzMW5oMwPkTtdNcGBkZmSuB2S2aPhcR0WHEJBs4f0IOPZKRkZl8Zr1ogmiiV8uiKSMjMwXMbtH0uwGYvyyZ1oZ22ls6p7lAMjIyc51ZPXqOT4jmvCWiX/Pi6TqWf2rhdJZoWghkWezZsyfoJHN5VkVmZiZlZWVBK7GRuuNMFAH7r4BxxEBGG4EyyQxPb0u2wOc6Eu/UqWKw79tsNmOxWPpYxvUmYC4TIJBpZrVasdlsU+ujOhjSLObuG5cEn//qn5+VXn36nWkszfRSVlYmmUymQd8vLy/v83rbtm1SWVnZZBdLkiRJqqyslAoKCoKvc3Nz+7yWJEkqKCiQcnJypqQ8s53W1lYpNze3z7qhvvupZqjvu/d3XFBQIJWUlEiSJEklJSXS9u3bpaysrD7H6n2c8vLy4Pbj5e677x7zvrO7ee5zgt8HQOrCBC6drZ/mAs1cpiJdcTA3dpvNxu7du4Ov161b18etyGKxzGj/xJnGQPZuAzndX85Eu6aP9vs2m819vuesrKzgdoO59Qfy4mHg654OZrdoSj5wCReWlIWJNFxowe30THOhZh4VFRVkZmb28Sy8nIDBR2lpadA8ITo6muLi4uCEaQEDiICX5UjJysrq461ZWVnZR8R7O4/LDE9WVlbQBT3wnfZutl7+XYL4jDMyMigtLSUvLy/4XZrN5uD3Hdi29zECxxlt+Qb6visqKvrY0hmNxmHNXLKzs8nIyAgat8wEk5bZ3acp+cDRDKFxpC5MRPJL1FoaSV829pnmRkKXz8kp+8VJPUeAJWFp6FXjMyO53HH7cioqKvo5bZeVlQX7RgMTpgX6SfPy8sZ88wZmkgz8qMxmM5s3b57R3pgel4fmWtuUnCs22YBGqxl2u/LycvLz89m0aVOwzzAwJcVIXdOHcl6fKDf13t93b7u4kVJQUIDVah3XPTfRzHLR9EPnBTAuJS4lGm1oCJfO1k+6aJ6yX2Tth9+e1HMEKL/2f8mKHPvg1mBu7r0ZzGn7ySefpKSkJCiegecBxuLG/tBDD1FWVhb8AQY6/meyaDbX2vjtD/4yJef66k/vIykjbtjtAt8HiD+99evXk5OTMyrX9MLCwkGd1wdyUzeZTOP6vmNiYkbVRRCYkrqkpASLxcKmTZsoLCyc9sGgWS2aNS0uNj74z2zNs7J161ZSFiR
      "text/plain": [
       "<Figure size 350x262.5 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\", \"CF4\"]\n",
    "\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "    plt.figure()\n",
    "    for i, sim in enumerate(sims):\n",
    "        r, B = get_bulkflow_simulation(sim, convert_to_galactic=True)\n",
    "        B = B[..., 0]\n",
    "\n",
    "        if sim == \"Carrick2015\":\n",
    "            B *= 0.43\n",
    "\n",
    "        if sim in [\"Carrick2015\", \"Lilow2024\", \"CLONES\"]:\n",
    "            plt.plot(r, B[0], label=simname_to_pretty(sim), color=cols[i])\n",
    "        else:\n",
    "            ylow, yhigh = np.percentile(B, [16, 84], axis=0)\n",
    "            plt.fill_between(r, ylow, yhigh, alpha=0.5,\n",
    "                             label=simname_to_pretty(sim), color=cols[i])\n",
    "\n",
    "    plt.xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "    plt.ylabel(r\"$|\\mathbf{B}_{\\rm sim}| ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    plt.xlim(5, 200)\n",
    "    plt.legend(ncols=2)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/bulkflow_simulations_restframe.pdf\", dpi=450)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Bulk flow in the CMB frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:26:54\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:28:43\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:20:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Carrick2015_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:19:39\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:27:13\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:28:48\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:19:59\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_Lilow2024_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:26:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:54:29\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:10:19\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:27:52\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2_main_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 12:58:15\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 16:59:53\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:16:15\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:39:34\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_csiborg2X_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:01:55\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:49:16\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:50:20\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:52:12\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CLONES_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:51:03\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_2MTF_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 30/08/2024 17:12:11\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_SFI_gals_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 14:07:26\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_i_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:23:47\n",
      "File:          /mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/samples_CF4_CF4_TFR_w1_bayes_zcmb_max_0.05_sample_alpha.hdf5\n",
      "Last modified: 12/09/2024 13:19:35\n"
     ]
    }
   ],
   "source": [
    "sims = [\"Carrick2015\", \"Lilow2024\", \"csiborg2_main\", \"csiborg2X\", \"CLONES\", \"CF4\"]\n",
    "# cats = [[\"LOSS\", \"Foundation\"], \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "\n",
    "data = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"bayes\",\n",
    "            sample_alpha=True, zcmb_max=0.05)\n",
    "        data[f\"{sim}_{cat}\"] = get_bulkflow(fname, sim)\n",
    "\n",
    "def get_ax_centre(ax):\n",
    "    # Get the bounding box of the specific axis in relative figure coordinates\n",
    "    bbox = ax.get_position()\n",
    "\n",
    "    # Extract the position and size of the axis\n",
    "    x0, y0, width, height = bbox.x0, bbox.y0, bbox.width, bbox.height\n",
    "\n",
    "    # Calculate the center of the axis\n",
    "    center_x = x0 + width / 2\n",
    "    center_y = y0 + height / 2\n",
    "    return center_x, center_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAw0AAAOiCAYAAAAyqzvbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9eXgc533ni35q633HvoPgKkqiJJCybMuxZQm0HcdWFpNSlJk7Z24SkjPnzsy9M09sWidPztxzcjMymcxM7syciSH5TObcTKJIhO3YiRPbhGTLtlaKIEWJorg1QBL70ht676p67x+FbgIkAQIkQIBkffCAYFdXv/VWVVfV+31/mySEENjY2NjY2NjY2NjY2MyDvNodsLGxsbGxsbGxsbFZ29iiwcbGxsbGxsbGxsZmQWzRYGNjY2NjY2NjY2OzILZosLGxsbGxsbGxsbFZEFs02NjY2NjY2NjY2NgsiC0abGxsbGxsbGxsbGwWxBYNNjY2NjY2NjY2NjYLYosGGxsbGxsbGxsbG5sFsUWDjY2NjY2NjY2Njc2C2KLBxsbGxsbGxsbGxmZBbNFgY2NjY2NjY2NjY7MgtmiwsbGxsbGxsbGxsVkQdbU7YGNzJ5HMFskVjRXdhtuhEPQ4VnQbtxvJUoasmV/RbXhkF0HNu6LbuOUUkqBnV3YbqgecwZXdhs3VTCchv8Ln1uUBv31ubWzuFmzRYGOzTCSzRQ587wOmpgsrup0qv5P9v3qfLRxmSJYy/GH0L5ksJld0O9WOIH/Q8Y/uHOFQSMK7fwj5yZXdjqsadvyBLRxuJdNJ6P5DiK/wuQ1Xw74/sIWDjc1dgi0abGyWiVzRYGq6gFtT8DiVFdlGtmBtI1c0CHquv/7zzz/PgQMHANi1a1fl/zt37uTdd98lHo9f9Znt27eTSCQ4f/48Bw8epLu7u/JeLBYjEolUXh86dIjOzk62b99ONBqd815nZyeHDh260V1dNFkzz2QxiVt24FGcK7MNo8BkMUnWzBPk+qKhr6+PPXv2VI7J7GO/0LFav349Bw4cYNeuXSuyH3PQs5ZgUNyWNWAlt6FnFy0aEokEe/bsoa+vD4COjg4OHDiw6O/Zzp07icViHD169JrtL8d3unyNdHR0VNp49tlnr3veenp6iMVi7N27d1HH4obJZy3B4HJb1oCV3EY+u2jRcCee2927d/Pwww/zta99rbJs3759RKNRDh8+XFnW19fH/v37K8t6e3vp7u6+JfdIG5vlwhYNNjbLjMep4HNpK9Z+rrQ496fyQ+v8+fMkEgm2b9/O+vXrKwOWSCRCT0/PnIdhNBolkUhUXn/ta1+rPAz7+vrYvn37NYVGIpHghRdeuDWD3XnwKE78KzX4BXJmcVHrJRIJnnjiCQ4dOkRXVxeJRIKXX355zvurfazmoHrA4V+59o3coleNRqNs376dAwcOVAZTfX19lQHc9Y5dX18fsVhswW0s13d6trgrX18dHR10dnbOu+1bfs5dHvCu4LnN2+d2586dHDp0aI5oePfdd4lGo3PW6+3tZefOnYAlNBKJxFXr2NisdexAaBubO5REIsGzzz4LQCgUYt++fXNmvnbt2jVnVg6gu7ubffv23dJ+3mm8++67AHR1dQHWsV/xmeU7hH379rF37945x6uzs5NQKLSoz+/Zs6fynb+VhEIhurq66O3tveXbvl24U8/tU089Nee9aDRKR0cHXV1d9PT0VJa/9NJLFSFy6NChiuXRxuZ2whYNNjZ3KJ2dnVfNjM02x+/cufMqy8KVlgebpVMWC7t377YHkUukt7f3hkVrb28vHR0dlZnrW0lvby8vv/zygtfO888/f1cL8jv13IZCITo6OirXek9PDzt37mTnzp2VSZpEIjHH5cnG5nbFFg02NncBiUSC7u5u9u/fP2f5rl27Kq4zvb29lQHvjbB//362b99e+b2bTe/9/f2AJRzC4fCcGUewj9W1KB+D6w2s5jt2+/btW/bZ24XO0549e9i+fTvhcJgDBw7Q398/b997enrYu3fvXSsi7+RzC9Z9tOxydfjwYZ566im6urqW7d5qY7NWsGMabGzuAnbv3s2hQ4euevDt27eP3bt3s3fvXrq7u2/K/H/LAnhvA0KhUGUQ0dPTw+7duzl//nzl+NvH6mrKx6bs3jEf1zp2PT09dHZ20tHRUQmyXQ4WOk9ln/jyDPpCbja7du2a479/t3Enn1uAp59+mt27dwNW4HQoFCIUChGJRCqxZeX3bWxuZ2xLg43NHc7u3bsrGUqupPwA7+vro6+vb8EgTpsbY9euXXR2di7rgOdOpaur66o4m8Vw5MiRSqDt7t27K0Gwt4Kuri46Ojo4ePDgguvd7fFCd/K57ezsJBaL0dPTM8eiUI5rsC0NNncKtqXBxmaZyRZWrrjbUtvevXs3zz777IJioGxtuN0HNFlj5epjLKXt3t5e+vr62Lt3L6FQiN7eXqLR6NodNKxkcbcltt3d3X1Vlq9y3M1C3+HZriu9vb3s379/3rScK0G53+Vzfi3KKTaff/75WxcYv5LF3ZbY9p18bsESCM8999yc/u7evZv9+/fftRYmmzsPWzTY2CwTbodCld9p1VFYZFrUG6HK78TtuH4diJ6ensosV5murq6r8oLv3buXQ4cO3bYZfjyyi2pHkMlictFpUW+EakcQj+y67no7duzg8OHDrFu3DrCsOa+88sqis8TcMlSPVXgtP7mktKhLxlW96DoQHR0d9Pf3s2fPnsrgKxQK8cILL6xc/5aBjo4OnnrqKfbv3z/vbHp51vmpp55a+Q65PFbhtfjkktKiLplw9aLrQNzJ5xYuuyhdaWno6+u76nP79u2rTCZs3759UTU+bGzWApIQQqx2J2xs7hSS2SK54soJBrDEiV0Nei7JUoasmV/RbXhk151TDbpMIbmylgawBINdDfrWM51cWUsDWILBrgZtY3PXYIsGGxsbGxsbGxsbG5sFsQOhbWxsbGxsbGxsbGwWxBYNNjY2NjY2NjY2NjYLYosGGxsbGxsbGxsbG5sFsUWDjY2NjY2NjY2Njc2C2KLBxsbGxsbGxsbGxmZBbNFgY2NjY2NjY2NjY7MgtmiwsbGxsbGxsbGxsVkQWzTY2NjY2NjY2NjY2CyILRpsbGxsbGxsbGxsbBbEFg02NjY2NjY2NjY2NgtiiwYbGxsbGxsbGxsbmwWxRYONjY2NjY2NjY2NzYLYosHGxsbGxsbGxsbGZkFs0WBjY2NjY2NjY2NjsyDqanfgVnHvvfeyfv36BdcZGhqiqanpprZzs22shT7cSW2shT6spTb6+voYHBy8qTZuF653za+Vc7IW2lgLfbiT2lgLfYC763p3Op2Ew+F538/n87hcrpvaht3G2urDndTGjXw+n8+Ty+Uqr3VdxzCMG+7DohB3CV/+8peXZZ2VbmMt9OFOamMt9GEttVFXV3fTbdwuXO94rZVzshbaWAt9uJPaWAt9EOLuut6vt69r5ZzcKW2shT7cSW0sRx8cDsdNt3E9bPekWTzzzDOr3sZa6MNysRb2ZS30YblYjn7c7MzlncRa+W6slTbWQh/upDbWQh/s6/0ya+V7sVbaWAt9uJPaWAt9cLvdy9CThZGEEGLFt7IGePLJJ/n+97+/2t24bbCP1/KyVo7nWunHreBu2tebxT5Wy8taOZ5rpR+3gvr6ekZHR1e7G7cNd9N341awFo7nrbgGbEuDzTVZC8r7TsI+njZrGfv7ubzYx9NmrWN/R5eXu+V42qLB5prcLRfArcI+njZrGfv7ubzYx9NmrWN/R5eXu+V42qLBxsbmjmRoaIgnn3ySF198cbW7YmNzS3nxxRd58sknGRoaWu2u2NjY3EHcNSlXbWxs7i6amppW3cfUxmY1eOo3n+Lx3/gCT/36rtXuio2NzQrz4osv8uKLL5LP51d8W3eNaPjo3ABf/vKX+a3f+q27xoxkY1OmfFOxZx5tbO4simaJsWKckUKM0UKcC7kxorkRknqGs5m76HrXszB9Efytq90TG5tbyjPPPMMzzzxDfX39im/rrhENOcXHF//FH9P1ifbV7oq
      "text/plain": [
       "<Figure size 830x954.5 with 24 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "    nrows = len(sims)\n",
    "    ncols = 3\n",
    "\n",
    "    figwidth = 8.3\n",
    "    fig, axs = plt.subplots(nrows, ncols, figsize=(figwidth, 1.15 * figwidth), sharex=True, )\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "    cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "    # fig.suptitle(f\"Calibrated against {catalogue}\")\n",
    "\n",
    "    for i, sim in enumerate(sims):\n",
    "        for j, catalogue in enumerate(cats):\n",
    "            r, B = data[f\"{sim}_{catalogue}\"]\n",
    "            c = cols[j]\n",
    "            for n in range(3):\n",
    "                ylow, ymed, yhigh = np.percentile(B[..., n], [16, 50, 84], axis=-1)\n",
    "                axs[i, n].fill_between(\n",
    "                    r, ylow, yhigh, alpha=0.5, color=c, edgecolor=c,\n",
    "                    label=catalogue_to_pretty(catalogue) if i == 1 else None)\n",
    "\n",
    "\n",
    "    # CMB-LG velocity\n",
    "    kwargs = {\"color\": \"mediumblue\", \"alpha\": 0.5, \"zorder\": 10}\n",
    "    for n in range(len(sims)):\n",
    "        axs[n, 0].fill_between([r.min(), 15.], [627 - 22, 627 - 22], [627 + 22, 627 + 22], label=\"CMB-LG\" if n == 0 else None, **kwargs)\n",
    "        axs[n, 1].fill_between([r.min(), 15.], [276 - 3, 276 - 3], [276 + 3, 276 + 3], **kwargs)\n",
    "        axs[n, 2].fill_between([r.min(), 15.], [30 - 3, 30 - 3], [30 + 3, 30 + 3], **kwargs)\n",
    "\n",
    "    # LCDM expectation\n",
    "    Rs,mean,std,mode,p05,p16,p84,p95 = np.load(\"/mnt/users/rstiskalek/csiborgtools/data/BulkFlowPlot.npy\")\n",
    "    m = Rs < 175\n",
    "    kwargs = {\"color\": \"black\", \"zorder\": 0, \"alpha\": 0.25}\n",
    "    for n in range(len(sims)):\n",
    "        axs[n, 0].fill_between(\n",
    "            Rs[m], p16[m], p84[m],\n",
    "            label=r\"$\\Lambda\\mathrm{CDM}$\" if n == 0 else None, **kwargs)\n",
    "\n",
    "    for n in range(3):\n",
    "        axs[-1, n].set_xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "\n",
    "    for n in range(len(sims)):\n",
    "        axs[n, 0].set_ylabel(r\"$|\\mathbf{B}| ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "        axs[n, 1].set_ylabel(r\"$\\ell ~ [\\mathrm{deg}]$\")\n",
    "        axs[n, 2].set_ylabel(r\"$b ~ [\\mathrm{deg}]$\")\n",
    "\n",
    "    for i, sim in enumerate(sims):\n",
    "        ax = axs[i, -1].twinx()\n",
    "        ax.set_ylabel(simname_to_pretty(sim), rotation=270, labelpad=7.5)\n",
    "        ax.set_yticklabels([])\n",
    "\n",
    "    # Watkins numbers\n",
    "    # for n in range(len(sims)):\n",
    "    #     rx = 150\n",
    "\n",
    "    axs[0, 0].set_xlim(r.min(), r.max())\n",
    "\n",
    "    axs[0, 0].legend()\n",
    "    handles, labels = axs[1, 0].get_legend_handles_labels()  # get the labels from the first axis\n",
    "    fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 0.975), ncol=len(cats) + 2)\n",
    "\n",
    "    fig.tight_layout(rect=[0, 0, 0.95, 0.95], h_pad=0.01)\n",
    "    fig.savefig(f\"../../plots/bulkflow_CMB.pdf\", dpi=450)\n",
    "    fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Full vs Delta comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogue = \"CF4_TFR_i\"\n",
    "simname = \"csiborg2X\"\n",
    "zcmb_max=0.05\n",
    "sample_beta = True\n",
    "sample_alpha = True\n",
    "\n",
    "fname_bayes = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"bayes\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "fname_mike = paths.flow_validation(\n",
    "    fdir, simname, catalogue, inference_method=\"mike\",\n",
    "    sample_alpha=sample_alpha, sample_beta=sample_beta,\n",
    "    zcmb_max=zcmb_max)\n",
    "\n",
    "\n",
    "X = []\n",
    "labels = [\"Full posterior\", \"Delta posterior\"]\n",
    "for i, fname in enumerate([fname_bayes, fname_mike]):\n",
    "    samples = get_samples(fname)\n",
    "    if i == 1:\n",
    "        print(samples.keys())\n",
    "\n",
    "    X.append(samples_to_getdist(samples, labels[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = [f\"a_{catalogue}\", f\"b_{catalogue}\", f\"c_{catalogue}\", f\"e_mu_{catalogue}\",\n",
    "          \"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"alpha_{catalogue}\"]\n",
    "# params = [\"beta\", f\"a_{catalogue}\", f\"b_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "# params = [\"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\", f\"mag_cal_{catalogue}\", f\"alpha_cal_{catalogue}\", f\"beta_cal_{catalogue}\", f\"e_mu_{catalogue}\"]\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 11})\n",
    "    g = plots.get_subplot_plotter()\n",
    "    g.settings.figure_legend_frame = False\n",
    "    g.settings.alpha_filled_add = 0.75\n",
    "    g.settings.fontsize = 12\n",
    "\n",
    "    g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "    # plt.gcf().suptitle(catalogue_to_pretty(catalogue), y=1.025)\n",
    "    plt.gcf().tight_layout()\n",
    "    plt.gcf().savefig(f\"../../plots/method_comparison_{simname}_{catalogue}.pdf\", dpi=300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Guilhem plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manticore vs linear comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"Carrick2015\", \"csiborg2X\"]\n",
    "catalogues = [\"LOSS\", \"Foundation\", \"2MTF\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "y_lnZ = np.full((len(catalogues), len(sims)), np.nan)\n",
    "\n",
    "for i, catalogue in enumerate(catalogues):\n",
    "    for j, simname in enumerate(sims):\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, simname, catalogue, inference_method=\"mike\",\n",
    "            sample_alpha=simname != \"IndranilVoid_exp\",\n",
    "            zcmb_max=zcmb_max)\n",
    "\n",
    "        y_lnZ[i, j] = - get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "    # y_lnZ[i] -= y_lnZ[i].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayes_factor = y_lnZ[:, 1] - y_lnZ[:, 0]\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "    plt.figure()\n",
    "\n",
    "    sns.barplot(x=np.arange(len(catalogues)), y=bayes_factor / np.log(10), color=\"#21456D\")\n",
    "    plt.xticks(\n",
    "        np.arange(len(catalogues)),\n",
    "        [catalogue_to_pretty(cat) for cat in catalogues],\n",
    "        rotation=35, fontsize=\"small\", minor=False)\n",
    "    plt.ylabel(r\"$\\log \\left(\\mathcal{Z}_{\\rm Manticore} / \\mathcal{Z}_{\\rm linear}\\right)$\")\n",
    "    plt.tick_params(axis='x', which='both', bottom=False, top=False)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\"../../plots/manticore_vs_carrick.png\", dpi=450)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## All possible things"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dipole magnitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "sim = \"IndranilVoid_gauss\"\n",
    "\n",
    "X = []\n",
    "for cat in cats:\n",
    "    fname = paths.flow_validation(\n",
    "        fdir, sim, cat, inference_method=\"mike\",\n",
    "        sample_mag_dipole=False,\n",
    "        sample_alpha=False, zcmb_max=0.05)\n",
    "    \n",
    "    if not exists(fname):\n",
    "        raise FileNotFoundError(fname)\n",
    "\n",
    "    samples = get_samples(fname, convert_Vext_to_galactic=False)\n",
    "\n",
    "    # keys = list(samples.keys())\n",
    "    # for key in keys:\n",
    "    #     if cat in key:\n",
    "    #         value = samples.pop(key)\n",
    "    #         samples[key.replace(f\"_{cat}\",'')] = value\n",
    "    \n",
    "    samples = samples_to_getdist(samples, catalogue_to_pretty(cat))\n",
    "    X.append(samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# params = [\"Vmag\", \"l\", \"b\", \"a_dipole_mag\", \"a_dipole_l\", \"a_dipole_b\"]\n",
    "params = [\"Vx\", \"Vy\", \"Vz\"]\n",
    "# params = [\"Vmag\", \"l\", \"b\"]\n",
    "\n",
    "with plt.style.context('science'):\n",
    "    g = plots.get_subplot_plotter()\n",
    "    g.settings.figure_legend_frame = False\n",
    "    g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "    g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "    # plt.gcf().suptitle(catalogue_to_pretty(cat), y=1.025)\n",
    "    plt.gcf().tight_layout()\n",
    "    plt.gcf().savefig(f\"../../plots/vext_{sim}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow | catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogues = [\"LOSS\", \"Foundation\", \"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "params = [\"Vmag\", \"beta\", \"sigma_v\"]\n",
    "\n",
    "for catalogue in catalogues:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim)\n",
    "          for sim in sims]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{catalogue}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "     plt.gcf().savefig(f\"../../plots/calibration_{catalogue}.png\", dpi=500, bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow | simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogues = [\"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "params = [\"Vmag\", \"l\", \"b\", \"beta\", \"sigma_v\"]\n",
    "\n",
    "for sim in sims:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim, catalogue)\n",
    "          for catalogue in catalogues]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{sim}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "     plt.gcf().savefig(f\"../../plots/calibration_{sim}.png\", dpi=500, bbox_inches='tight')\n",
    "     plt.gcf().show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stacking vs marginalising CB boxes\n",
    "\n",
    "#### $V_{\\rm ext}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"csiborg2X\"\n",
    "catalogue = \"2MTF\"\n",
    "key = \"Vext\"\n",
    "\n",
    "X = [get_samples(sim, catalogue, nsim=nsim, convert_Vext_to_galactic=False)[key] for nsim in range(20)]\n",
    "Xmarg = get_samples(sim, catalogue, convert_Vext_to_galactic=False)[key]\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharey=True)\n",
    "fig.suptitle(f\"{simname_to_pretty(sim)}, {catalogue}\")\n",
    "fig.subplots_adjust(wspace=0.0, hspace=0)\n",
    "\n",
    "for i in range(3):\n",
    "    for n in range(20):\n",
    "        axs[i].hist(X[n][:, i], bins=\"auto\", alpha=0.25, histtype='step',\n",
    "                    color='black', linewidth=0.5, density=1, zorder=0,\n",
    "                    label=\"Individual box\" if (n == 0 and i == 0) else None)\n",
    "\n",
    "axs[i].hist(np.hstack([X[n][:, i] for n in range(20)]), bins=\"auto\",\n",
    "            histtype='step', color='blue', density=1,\n",
    "            label=\"Stacked individual boxes\" if i == 0 else None)\n",
    "axs[i].hist(Xmarg[:, i], bins=\"auto\", histtype='step', color='red',\n",
    "            density=1, label=\"Marginalised boxes\" if i == 0 else None)\n",
    "    \n",
    "axs[0].legend(fontsize=\"small\", loc='upper left', frameon=False)\n",
    "\n",
    "axs[0].set_xlabel(r\"$V_{\\mathrm{ext}, x} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[1].set_xlabel(r\"$V_{\\mathrm{ext}, y} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[2].set_xlabel(r\"$V_{\\mathrm{ext}, z} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[0].set_ylabel(\"Normalized PDF\")\n",
    "fig.tight_layout()\n",
    "fig.savefig(f\"../../plots/consistency_{sim}_{catalogue}_{key}.png\", dpi=450)\n",
    "fig.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### $\\beta$ and others"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"csiborg2_main\"\n",
    "catalogue = \"Pantheon+\"\n",
    "key = \"alpha\"\n",
    "\n",
    "X = [get_samples(sim, catalogue, nsim=nsim, convert_Vext_to_galactic=False)[key] for nsim in range(20)]\n",
    "Xmarg = get_samples(sim, catalogue, convert_Vext_to_galactic=False)[key]\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.title(f\"{simname_to_pretty(sim)}, {catalogue}\")\n",
    "for n in range(20):\n",
    "    plt.hist(X[n], bins=\"auto\", alpha=0.25, histtype='step',\n",
    "             color='black', linewidth=0.5, density=1, zorder=0,\n",
    "             label=\"Individual box\" if n == 0 else None)\n",
    "\n",
    "plt.hist(np.hstack([X[n] for n in range(20)]), bins=\"auto\",\n",
    "         histtype='step', color='blue', density=1,\n",
    "         label=\"Stacked individual boxes\")\n",
    "plt.hist(Xmarg, bins=\"auto\", histtype='step', color='red',\n",
    "         density=1, label=\"Marginalised boxes\")\n",
    "\n",
    "plt.legend(fontsize=\"small\", frameon=False, loc='upper left', ncols=3)\n",
    "plt.xlabel(names_to_latex([key], True)[0])\n",
    "plt.ylabel(\"Normalized PDF\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(f\"../../plots/consistency_{sim}_{catalogue}_{key}.png\", dpi=450)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SN/TFR Calibration consistency"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# catalogues = [\"LOSS\", \"Foundation\", \"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "catalogues = [\"Pantheon+\"]\n",
    "sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "\n",
    "for catalogue in catalogues:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim)\n",
    "          for sim in sims]\n",
    "\n",
    "     if \"Pantheon+\" in catalogue or catalogue in [\"Foundation\", \"LOSS\"]:\n",
    "          params = [\"alpha_cal\", \"beta_cal\", \"mag_cal\", \"e_mu\"]\n",
    "     else:\n",
    "          params = [\"aTF\", \"bTF\", \"e_mu\"]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{catalogue}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "     # plt.gcf().savefig(f\"../../plots/calibration_{catalogue}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $V_{\\rm ext}$ comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogues = [\"LOSS\"]\n",
    "# sims = [\"Carrick2015\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "sims = [\"Carrick2015\"]\n",
    "params = [\"Vmag\", \"l\", \"b\"]\n",
    "\n",
    "for sim in sims:\n",
    "     X = [samples_to_getdist(get_samples(sim, catalogue), sim, catalogue)\n",
    "          for catalogue in catalogues]\n",
    "\n",
    "     g = plots.get_subplot_plotter()\n",
    "     g.settings.figure_legend_frame = False\n",
    "     g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "     g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "     plt.gcf().suptitle(f'{simname_to_pretty(sim)}', y=1.025)\n",
    "     plt.gcf().tight_layout()\n",
    "    #  plt.gcf().savefig(f\"../../plots/calibration_{sim}.png\", dpi=500, bbox_inches='tight')\n",
    "     plt.gcf().show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Bulk flow in the simulation rest frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sims = [\"Carrick2015\", \"csiborg1\", \"csiborg2_main\", \"csiborg2X\"]\n",
    "convert_to_galactic = False\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n",
    "cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "for i, sim in enumerate(sims):\n",
    "    r, B = get_bulkflow_simulation(sim, convert_to_galactic=convert_to_galactic)\n",
    "    if sim == \"Carrick2015\":\n",
    "        if convert_to_galactic:\n",
    "            B[..., 0] *= 0.43\n",
    "        else:\n",
    "            B *= 0.43\n",
    "\n",
    "    for n in range(3):\n",
    "        ylow, ymed, yhigh = np.percentile(B[..., n], [16, 50, 84], axis=0)\n",
    "        axs[n].fill_between(r, ylow, yhigh, color=cols[i], alpha=0.5, label=simname_to_pretty(sim) if n == 0 else None)\n",
    "\n",
    "axs[0].legend()\n",
    "if convert_to_galactic:\n",
    "    axs[0].set_ylabel(r\"$B ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    axs[1].set_ylabel(r\"$\\ell_B ~ [\\degree]$\")\n",
    "    axs[2].set_ylabel(r\"$b_B ~ [\\degree]$\")\n",
    "else:\n",
    "    axs[0].set_ylabel(r\"$B_{x} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    axs[1].set_ylabel(r\"$B_{y} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "    axs[2].set_ylabel(r\"$B_{z} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "for n in range(3):\n",
    "    axs[n].set_xlabel(r\"$R ~ [\\mathrm{Mpc}]$\")\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../../plots/bulkflow_simulations_restframe.png\", dpi=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bulk flow in the CMB rest frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"csiborg2_main\"\n",
    "catalogues = [\"Pantheon+\", \"2MTF\", \"SFI_gals\"]\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(15, 5), sharex=True)\n",
    "cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "# fig.suptitle(f\"Calibrated against {catalogue}\")\n",
    "\n",
    "for i, catalogue in enumerate(catalogues):\n",
    "    r, B = get_bulkflow(sim, catalogue, sample_beta=True, convert_to_galactic=True,\n",
    "                        weight_simulations=True, downsample=3)\n",
    "    c = cols[i]\n",
    "    for n in range(3):\n",
    "        ylow, ymed, yhigh = np.percentile(B[..., n], [16, 50, 84], axis=-1)\n",
    "        axs[n].plot(r, ymed, color=c)\n",
    "        axs[n].fill_between(r, ylow, yhigh, alpha=0.5, color=c, label=catalogue)\n",
    "\n",
    "\n",
    "# CMB-LG velocity\n",
    "axs[0].fill_between([r.min(), 10.], [627 - 22, 627 - 22], [627 + 22, 627 + 22], color='black', alpha=0.5, zorder=0.5, label=\"CMB-LG\", hatch=\"x\")\n",
    "axs[1].fill_between([r.min(), 10.], [276 - 3, 276 - 3], [276 + 3, 276 + 3], color='black', alpha=0.5, zorder=0.5, hatch=\"x\")\n",
    "axs[2].fill_between([r.min(), 10.], [30 - 3, 30 - 3], [30 + 3, 30 + 3], color='black', alpha=0.5, zorder=0.5, hatch=\"x\")\n",
    "\n",
    "# LCDM expectation\n",
    "Rs,mean,std,mode,p05,p16,p84,p95 = np.load(\"/mnt/users/rstiskalek/csiborgtools/data/BulkFlowPlot.npy\")\n",
    "m = Rs < 175\n",
    "axs[0].plot(Rs[m], mode[m], color=\"violet\", zorder=0)\n",
    "axs[0].fill_between(Rs[m], p16[m], p84[m], alpha=0.25, color=\"violet\",\n",
    "                    zorder=0, hatch='//', label=r\"$\\Lambda\\mathrm{CDM}$\")\n",
    "\n",
    "for n in range(3):\n",
    "    axs[n].set_xlabel(r\"$r ~ [\\mathrm{Mpc} / h]$\")\n",
    "\n",
    "axs[0].legend()\n",
    "axs[0].set_ylabel(r\"$B ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[1].set_ylabel(r\"$\\ell_B ~ [\\mathrm{deg}]$\")\n",
    "axs[2].set_ylabel(r\"$b_B ~ [\\mathrm{deg}]$\")\n",
    "\n",
    "axs[0].set_xlim(r.min(), r.max())\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(f\"../../plots/bulkflow_{sim}_{catalogue}.png\", dpi=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Smoothing scale dependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simname = \"Carrick2015\"\n",
    "catalogue = \"Pantheon+\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Goodness-of-fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scales = [0, 4, 8, 16, 32]\n",
    "\n",
    "y = np.asarray([get_gof(\"BIC\", simname, catalogue, ksmooth=i)\n",
    "                for i in range(len(scales))])\n",
    "ymin = y.min()\n",
    "\n",
    "y -= ymin\n",
    "y_CF4 = get_gof(\"BIC\", \"CF4\", catalogue) - ymin\n",
    "y_CF4gp = get_gof(\"BIC\", \"CF4gp\", catalogue) - ymin\n",
    "\n",
    "plt.figure()\n",
    "plt.axhline(y[0], color='blue', label=\"Carrick+2015, no smoothing\")\n",
    "plt.plot(scales[1:], y[1:], marker=\"o\", label=\"Carrick+2015, smoothed\")\n",
    "\n",
    "plt.axhline(y_CF4, color='red', label=\"CF4, no smoothing\")\n",
    "\n",
    "plt.xlabel(r\"$R_{\\rm smooth} ~ [\\mathrm{Mpc}]$\")\n",
    "plt.ylabel(r\"$\\Delta \\mathrm{BIC}$\")\n",
    "plt.legend(ncols=1)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"../../plots/test_smooth.png\", dpi=450)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = \"Carrick2015\"\n",
    "catalogue = \"Pantheon+\"\n",
    "\n",
    "\n",
    "X = [samples_to_getdist(get_samples(sim, catalogue, ksmooth=ksmooth), ksmooth)\n",
    "     for ksmooth in [0, 1, 2]]\n",
    "\n",
    "params = [\"Vmag\", \"l\", \"b\", \"sigma_v\", \"beta\"]\n",
    "# if \"Pantheon+\" in catalogue or catalogue in [\"Foundation\", \"LOSS\"]:\n",
    "#      params += [\"alpha_cal\", \"beta_cal\", \"mag_cal\", \"e_mu\"]\n",
    "# else:\n",
    "#      params += [\"aTF\", \"bTF\", \"e_mu\"]\n",
    "\n",
    "\n",
    "\n",
    "g = plots.get_subplot_plotter()\n",
    "g.settings.figure_legend_frame = False\n",
    "g.settings.alpha_filled_add = 0.75\n",
    "\n",
    "g.triangle_plot(X, params=params, filled=True, legend_loc='upper right')\n",
    "plt.gcf().suptitle(f'{catalogue}', y=1.025)\n",
    "plt.gcf().tight_layout()\n",
    "plt.gcf().savefig(f\"../../plots/calibration_{catalogue}.png\", dpi=500, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Void testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evidence comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "\n",
    "sims = [\"no_field\", \"IndranilVoid_exp\"]\n",
    "cats = [\"LOSS\", \"Foundation\", \"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "neglnZ = {}\n",
    "kfound = []\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        sample_alpha = sim not in [\"IndranilVoid_exp\", \"no_field\"]\n",
    "        fname = paths.flow_validation(\n",
    "            fdir, sim, cat, inference_method=\"mike\",\n",
    "            sample_alpha=sample_alpha, zcmb_max=zcmb_max)\n",
    "        \n",
    "\n",
    "        neglnZ[f\"{sim}_{cat}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simA = sims[0]\n",
    "simB = sims[1]\n",
    "\n",
    "print(f\"lnZ_({simA}) - lnZ_({simB})\\n\")\n",
    "for cat in cats:\n",
    "    lnZ_A = - neglnZ[f\"{simA}_{cat}\"]\n",
    "    lnZ_B = - neglnZ[f\"{simB}_{cat}\"]\n",
    "    print(f\"{cat:15s} {lnZ_A - lnZ_B:.1f}\")\n",
    "\n",
    "\n",
    "print(f\"\\n(Positive -> preference for {simA})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Goodness-of-fit comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "no_Vext = True\n",
    "\n",
    "sims = [\"IndranilVoid_exp\", \"IndranilVoid_gauss\", \"IndranilVoid_mb\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "neglnZ = {}\n",
    "kfound = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        kfound[f\"{sim}_{cat}\"] = []\n",
    "        for ksim in range(500):\n",
    "            sample_alpha = False\n",
    "            fname = paths.flow_validation(\n",
    "                fdir, sim, cat, inference_method=\"mike\", nsim=ksim,\n",
    "                sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "                sample_beta=True,\n",
    "                no_Vext=no_Vext, verbose_print=False)\n",
    "\n",
    "            if not exists(fname):\n",
    "                continue\n",
    "\n",
    "            kfound[f\"{sim}_{cat}\"].append(ksim)\n",
    "            neglnZ[f\"{sim}_{cat}_{ksim}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "\n",
    "neglnZ_no_field = {}\n",
    "neglnZ_dipole = {}\n",
    "sim = \"no_field\"\n",
    "for cat in cats:\n",
    "    sample_alpha = False\n",
    "    fname = paths.flow_validation(\n",
    "        fdir, sim, cat, inference_method=\"mike\",\n",
    "        sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "        no_Vext=True, verbose_print=False)\n",
    "\n",
    "    if not exists(fname):\n",
    "        continue\n",
    "\n",
    "    neglnZ_no_field[f\"{cat}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n",
    "\n",
    "    fname = paths.flow_validation(\n",
    "        fdir, sim, cat, inference_method=\"mike\",\n",
    "        sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "        no_Vext=None, verbose_print=False)\n",
    "\n",
    "    if not exists(fname):\n",
    "        continue\n",
    "\n",
    "    neglnZ_dipole[f\"{cat}\"] = get_gof(\"neg_lnZ_harmonic\", fname)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving to `../../plots/void_goodness_of_fit_observer_no_Vext.png`.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAIUCAYAAAAzC0wXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd1yT1/7A8U8Swh5hizjBvQW0Vu1SsLVVuqTWerurdN/b21a0t+t3O/G2t7f7gh23kyp0aHfB2lZb2ypxTyS4ANlhr5Dz++MxUQQUJZAA591XXkieJ+f5JqR58n3OOd+jEkIIJEmSJEmSJEmSegG1vQOQJEmSJEmSJEnqKjIBkiRJkiRJkiSp15AJkCRJkiRJkiRJvYZMgCRJkiRJkiRJ6jVkAiRJkiRJkiRJUq8hEyBJkiRJkiRJknoNmQBJkiRJkiRJktRryARIkiRJkiRJkqReQyZAkiRJkiRJkiT1Gk72DsCe+vXrR0RExBn3y83NJTQ0tEv3s8cxZWwyNhmb48fm6K9HXV0du3btOmN7vUF7zjGO/veUscnYZGwyNkc4JoBer+fo0aNn3K9dRC8WHBzcrv3mzp3b5fvZ45jt3U/Gdm77ydjObT8Zm/2P2d795s6d2+72eoP2nGMc/e/Z1cds734ytnPbT8Z2bvvJ2Ox/TCHa/729PeQQuHZYsGBBl+/X3rbayx6x2eM52ONv1V72iq27Pwf5/0LnHtNe/z/3BnV1dcTGxpKSktLmPj3h7yk/FzqvrbMhz5Gdu589jumo/y909fNMSUkhNjaWurq6drXXHiohhLBZa91Mnz59OHbsmL3D6JDY2FjWrFlj7zA6RD4H++vu8YN8Do4gNjYWoFs/B1vq7ueY7v5+BPkcHEV3fw7dPX7oGc/Blp+pvboHqD3jDR1dT7gqK5+D/XX3+EE+B0fQ3eO3te5+jukJf0/5HBxDd38O3T1+6BnPwZafqb26B6gnZMOSJEmORH6uniBfC0mSJNux5Wdqr+4BkiRJkiRJkiSpd5EJkCRJkiRJkiRJvUavXgdIknqSvdWHMTZWM0U30t6h2I0Qgp1HjPTRuRHo7WrvcKRerq6uDr1e3+w+X19fBg8eTF1dHbt3727xGMu6Qfv27aO6urrZtkGDBuHn50dRURFHjhxpts3Ly4uhQ4fS1NTEtm3bWrQ7duxYtFot2dnZlJeXN9sWGhpKcHAwZWVl5OTkNNvm5ubGyJHKZ8qWLVs4ddT8yJEjcXNz49ChQ5SUlDTbFhwcTGhoKJWVlWRlZTXbptVqGTt2LAA7duygsbGx2fahQ4fi5eVFbm4uBQUFzbb5+/szcOBAamtr2bNnT7NtKpWKiRMnArBnzx5qa2ubbR88eDC+vr4UFBSQm5vbbJuPjw/h4eE0NjayY8cOTjV+/Hg0Gg1ZWVlUVlY229a/f38CAwMpLS3l4MGDzbZ5eHgwfPhwgBbvB4BRo0bh6upKTk4OZWVlzbaFhIQQEhJCRUUFBw4caLbNxcWF0aNHA7B9+3ZMJlOz7cOGDcPT05OjR49SWFjYbFtAQAADBgygpqaGvXv3NtumVquZMGECALt3725ReSssLAydTsexY8fIy8trtk2n0xEWFkZDQwM7d+5s8VwnTJiAWq1m//79VFVVNds2YMAAAgICKC4u5vDhw822eXp6MmzYMMxmM1u3bm3R7pgxY3B2dsZgMGA0Gptt69u3L3369MFoNGIwGJptc3V1ZdSoUQBs3boVs9ncbPuIESNwd3fn8OHDFBcXN9sWFBREv379qKqqYv/+/c22OTk5MW7cOAB27dpFfX19s+1DhgzB29ub/Px88vPzm22TnxGKM31GdJjNCmp3Q46+XgUgoqOjm926q+zs7G4df3sZjUaRl5fX5bfv9q0X3qtmCz6+UFz6/QNi/YHN7Xqc0Wi090tmE01NZvFV5hFxyZPfCc8bPxKeN34kpj36jXh85RaxN7dnPMfuwtE/V7vSsGHDBNDstnDhQiGEEFlZWS22nXxKnjJlSottH3zwgRBCiNdee63FtlmzZgkhhCgvL2+13cLCQiGE8vc5dduLL74ohBBi1apVLbZNnDjRGpOzs3OL7Tt37hRCCHH77be32LZ06VIhhBDr1q1rsS00NNTabmhoaIvt69atE0IIsXTp0hbbbr/9diGEEDt37myxzdnZ2druxIkTW2xftWqVEEKIF198scU2y3u3sLCw1dewvLxcCCHErFmzWmx77bXXhBBCfPDBBy22TZkyxRpTa+1mZWUJIYRYuHBhi21PPPGEEEKI7777rsW28PBwa7sBAQEttq9Zs0bk5eWJRYsWtdh28803K+eOVtr18PCwniOGDh3aYvu7774r8vLyWv3bzJs3TwghxJEjR1p9rnV1dUIIIS666KIW21asWCGEEGLFihUttl100UVCCCHq6upabffIkSNCCCHmzZvXYtuzzz4rhBBi9erVLbaNGjXK+hp6eXm12J6ZmSmEEOLuu+9use2BBx4QQgjx22+/tdgWEBBgbTc8PLzF9u+++04IIcQTTzzRYpv8jGj7M8KW5xdZBMGBJ6iqVKoWmbQjMRgMhIWFAcpVrUWLFpGZmdnmvvHx8aSnp3dliG1rqIRf7obpr4Crr02aLC8v57XXXrP9VYozKGmo4PuSTHyc3Bnu0Z+tFQeoNtczymMAUd7DUKlUbT5Wq9Vy77334uPj04UR29aeo0ZuffM3dh0xcv6wQO6fPZLKukbW7shn7Y58ahuaeOfuqVw+sZ+9Q+0VHP1ztSvNmjWL559/vtl98uquQvYAnWCLHqBff/2VlStXNvv7eHt7o9Vqqa6ubtGL4+LigqenJyaTqcX7QaVS4efnB4DRaKSpqanZdi8vL5ydnamtraW2tpZrr70WLy8vQPYAWcgeoBNs+Rnx5JNP2uz80qsToMjISEJDQ1mwYIFDlgd09AQoMjKyzYTnVA6XAOWug9UzYPZqGBxrkybz8/NJSkrimmuuITAw0CZtnklhfRmrCn5B5+TJvOALcNU4YzI3sal8H7+W7+bKoCkMdW/9i39RURGfffYZ8fHxhISEdEm8trY3t5zLn1tLkI8rL94UxbThQc2219SbuOO/v/G1PpfnF0Zw16zhdoq050tJSSElJYXc3Nx2fy70dDIZlLqKPc4/PeEcInUvtvxM7dVzgEJDQ8/4QtbUm9ifX9Epxx8W4o27y+n/BHFxcdZ/x8TEsHjxYhISEggPD2fx4sUsX74cnU7Hddddx8yZM4mKiiIjI4P4+HiWLFlCfHw8BoOB0tJSVqxYQUREBAaDgbi4OOu+iYmJzJs3D8C6v06nIzU1tc19ExISrNvi4+MJCwuzJjhGo5FFixZhNBqt7TicioPKz/Ks0+52LgIDA7vkZNBoNvHx4V8ZGDqAm0KicdU4W7f1D+1HdZ5guymP6X0iUat6Xr2T/fkVzHleSX6+WjqDAK+Wc37cXZz48L4LeGzlVpZ8mMnBwiqeXxhx2l4x6dxYLiRZFkOVJKnrddX5R5K6u16dALXH/vwKLnj8u05pe/0/L2PCIL/T7tNa8pCYmEhcXBxGo5Hs7GySkpKsXb5JSUmA0oWq0+nQ6XTWpCQyMpLs7GxA6ZHJzMy09szMmzeP5ORk6/4ZGRmkpaVZE6ZT901MTCQjI8Ma38ndyicnPXFxcWRkZFiHyjmMyoPKT6PtE6CusrUym6qmOm4NvbRZ8mMR7R9B0tGv2VqZTYT3UDtE2HkMBZXMeX4tvp4ufJnQevJjoVareGbBRAYGevDg+5sZO9CXv1zgYO9HqUfKzc0lNjbWYUcZSJIkdQcnjzCwFZkAncGwEG/W//OyTmv7XK1YsQJfX98W44UtoqKiWLlyJcuWLQOUpASwJkqWhMTPz4/S0lKAZkkOnOh9am3fM9Hr9WzevLnV8bgOoeL4GNZO6AHqCmZhZqNxD6M8BuCn9Wp1nxAXP8Z6DmJd6TbGeg5Gq+4Z/7s3mszc9NoG3J2d+CphRrurvS2OHsaWnFIeen8z5w8LJDy49ddNkmylPaMMJEmSpNPrjBEGPeMbUSdyd3E6Yy+NPVgKDix
      "text/plain": [
       "<Figure size 830x539.5 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "    figwidth = 8.3 \n",
    "    fig, axs = plt.subplots(2, 2, figsize=(figwidth, 0.65 * figwidth))\n",
    "\n",
    "    for n, cat in enumerate(cats):\n",
    "        i, j = n // 2, n % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            x = kfound[f\"{sim}_{cat}\"]\n",
    "            y = [neglnZ[f\"{sim}_{cat}_{ksim}\"] / np.log(10) for ksim in x]\n",
    "            x = np.array(x) * 0.674\n",
    "            ax.plot(x, y, label=simname_to_pretty(sim))\n",
    "        \n",
    "        # if no_Vext is None:\n",
    "        #     y_no_field = neglnZ_no_field[cat] / np.log(10)\n",
    "        #     if cat != \"CF4_TFR_w1\":\n",
    "        #         ax.axhline(y_no_field, color=\"black\", ls=\"--\", label=\"No peculiar velocity\")\n",
    "        y_no_field = neglnZ_no_field[cat] / np.log(10)\n",
    "        ax.axhline(y_no_field, color=\"black\", ls=\"--\", label=\"No peculiar velocity\")\n",
    "\n",
    "        y_dipole = neglnZ_dipole[cat] / np.log(10)\n",
    "        ax.axhline(y_dipole, color=\"black\", ls=\":\", label=\"Constant dipole\")\n",
    "\n",
    "        ax.text(0.5, 0.9, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, #fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='center',\n",
    "            bbox=dict(facecolor='white', alpha=0.5),\n",
    "            )\n",
    "\n",
    "        if n == 0:\n",
    "            ax.legend(fontsize=\"small\", loc=\"upper left\")\n",
    "\n",
    "        ax.set_ylabel(r\"$-\\Delta \\log \\mathcal{Z}$\")\n",
    "        ax.set_xlabel(r\"$R_{\\rm offset} ~ [\\mathrm{Mpc} / h]$\")\n",
    "        ax.set_xlim(0)\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fname = f\"../../plots/void_goodness_of_fit_observer.png\"\n",
    "    if no_Vext:\n",
    "        fname = fname.replace(\".png\", \"_no_Vext.png\")\n",
    "    print(f\"Saving to `{fname}`.\")\n",
    "    fig.savefig(fname, dpi=450)\n",
    "    fig.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Single parameter radial dependence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcmb_max = 0.05\n",
    "key = \"beta\"\n",
    "# key_label = r\"$\\sigma_v ~ [\\mathrm{km} / \\mathrm{s}]$\"\n",
    "# key_label = r\"$|\\mathbf{V}_{\\rm ext}| ~ [\\mathrm{km} / \\mathrm{s}]$\"\n",
    "key_label = r\"$\\beta$\"\n",
    "no_Vext = True\n",
    "\n",
    "sims = [\"IndranilVoid_exp\", \"IndranilVoid_gauss\", \"IndranilVoid_mb\"]\n",
    "cats = [\"2MTF\", \"SFI_gals\", \"CF4_TFR_i\", \"CF4_TFR_w1\"]\n",
    "\n",
    "data_mean = {}\n",
    "data_std = {}\n",
    "kfound = {}\n",
    "for sim in sims:\n",
    "    for cat in cats:\n",
    "        kfound[f\"{sim}_{cat}\"] = []\n",
    "        for ksim in range(500):\n",
    "            sample_alpha = False\n",
    "            fname = paths.flow_validation(\n",
    "                fdir, sim, cat, inference_method=\"mike\", nsim=ksim,\n",
    "                sample_alpha=sample_alpha, zcmb_max=zcmb_max,\n",
    "                sample_beta=True,\n",
    "                no_Vext=no_Vext, verbose_print=False)\n",
    "\n",
    "            if not exists(fname):\n",
    "                continue\n",
    "\n",
    "            kfound[f\"{sim}_{cat}\"].append(ksim)\n",
    "            with File(fname, 'r') as f:\n",
    "                x = f[f\"samples/{key}\"][...]\n",
    "            if key == \"Vext\":\n",
    "                x = np.linalg.norm(x, axis=-1)\n",
    "\n",
    "            data_mean[f\"{sim}_{cat}_{ksim}\"] = x.mean()\n",
    "            data_std[f\"{sim}_{cat}_{ksim}\"] = x.std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving to `../../plots/void_beta_per_observer_no_Vext.png`.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAIUCAYAAAAzC0wXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9aXAUabrnif7cPTz2TaEdCQQSkECSmyAza8vaUlSfrfL0qQOVp6an73T3WCVzZ86dD2NjRdeXa3fmSxtYm821mR4bg7rTPX2m+2TnSar6VJ2l6hSqLWvJBVAmkOwotO9S7LuHu98PHhFSIAkESEgRen9pYRHu/rrHG0DGG39/nuf/SKZpmggEAoFAIBAIBALBNkDe7AkIBAKBQCAQCAQCwdNCCCCBQCAQCAQCgUCwbRACSCAQCAQCgUAgEGwbhAASCAQCgUAgEAgE2wYhgAQCgUAgEAgEAsG2QQgggUAgEAgEAoFAsG0QAkggEAgEAoFAIBBsG4QAEggEAoFAIBAIBNsGIYAEAoFAIBAIBALBtsG22RN4WnR2dtLb2/vQcRMTE3R0dDzVcZvxnmLc1nlPMW7rvKcY92TjJiYmyOVyXL9+/aHXqye28vqy3uO28tzEuK3znmLc1nnPeho3MDDA+Pj4Q8etCXOb0NrauqZxX//615/6uM14TzFu67ynGLd13lOMe7JxX//619d8vXpiK68v6z1uK89NjNs67ynGbZ33rKdxa/2uXQsiBe4+vvWtbz31cWu91lrZjM/wKOPWivi72Phxa0X8XWz8uLVSL59jO1IPf3f18Bk24n238ueoh8+wmePWivi72Phx68q6SaktTiAQML/+9a+bf/mXf7nZU3ls6uXOaj18jnr4DKZZH5+jHj6Dadb+5/jLv/xLs7W11ezt7d3sqTx11vOu5GZS6/8GTbM+PoNp1sfnqIfPYJr18Tnq4TOY5vp+126bGqCenh5+9KMfbfY0noh6uQNbD5+jHj4D1MfnqIfPALX/Ocrzf/vttzd5Jk8fl8vFG2+8wbe+9a2a/nus5bmXqYfPAPXxOerhM0B9fI5a/wxvv/02b7/9Ni6Xa92uKZmmaa7b1bYwb7zxRs0LIIFAINjqbMfv2u34mQUCgeBps57ftaIGSCAQCAQCgUAgEGwbhAASCAQCgUAgEAgE24aaF0D9/f0MDAzQ39/PyZMnN3s6AoFAIBAIBAKBYAtT8yYIx44do1zGdPr0ac6dO8dbb721ybMSbCr5ODgCT+3t4vE4mUzmqb3feuB2uwkEnt6f0UYjSRJ9fX1V+y5cuLBJs3kywuEwJ0+erNn5CwQ1j66Bom72LNbEZqw/9bZ+CLYnNS+ABgcHK68jkQhHjx7dxNkItgTTv4WuP3gqbxWPx/k3/+bfoGnaU3m/9UJVVf78z/+8rhaxrSwYwuEw3d3dgNXJ+tvf/jaXL1/e5FkJBIIVyc2Dp32zZ/FQNmv9qcf1Q7D9qHkBVP5Rcf78efr6+ujt7d3kGQk2ndhtaD4K7pYNf6tMJoOmaXzjG9+gubl5w99vPZibm+MHP/gBmUxGLGBPiRMnTlQET29vrxA/dcbExERd2GALSmTnakIAbcb6I9YPwWZQtsGemJhYt2vWvAAC645qJBKhsbFx1THlBaqMWKjqGC0Dc5eeWhQIoLm5mfb2rb9gbhSZfJE7U4kNufb+dj9ux8O/qk6cOFF5fezYMd566y1OnTpFT08Pb731FmfOnCEYDPLNb36T119/naNHj1ZqB7/zne9w8uRJwuEwkUiE733ve/T29hIOhzlx4kRl7OnTpzl+/DhAZXwwGOTdd99ddeypU6cqx06ePEl3d3clxS0Wi/Htb3+bWCxWuU6tUV6YyqznAlUrdHR0CBvseiI3t9kzeCS2+/ojqH/Kv9mX/o5/UupCAPX29tLb21v5gXH27NllY8QCtU3QNTA0WLgKncdqJo+71rkzleC1//dPNuTav/6ff48Xd4ceOm4l8XD69GlOnDhBLBZjcHCQs2fPEovFACrfEw0NDQSDQYLBYEWUHDlypJJeGw6HuXz5cqU25/jx45w7d64yvr+/n/Pnz1cE0/1jT58+TX9/f2V+4XC4Mr+loufEiRP09/dXotq1wv03k9ZzgRIINoXs7GbPQCAQbDA1LYDOnz/P2bNnK7n/L7/8Mu+8884mz0qwqehZ67mYg8in0PzS5s5nm7C/3c+v/+ff27BrPwnf+973aGhoIBqNrnj86NGjvPPOO3z3u98FLFECVIRSWZCEQiEikQhAlciBxejTSmMfxsDAAJcuXSIcDlfeUyAQbCLZ2ooACQSCR6emBVAwGKxKe7l48SJvvvnmJs5IsOkUl7jhzF0SAugp4XbY1hSl2QzKhgPf/va3K9GWpULj0qVL/Pf//X/Pu+++S19fXyUNLhgMEolECIWWf64jR44QDAY5ffp0ZV84HF5x7P3vt5Rz585x+fJlTp8+XWXoIhAINglDh3wMTBMkabNnIxAINoia7gPU19dHKBTi/PnznDt3ju7ubr7zne9s9rQEm4m2RAClJiAzvWlTOXfuHD09PfT09HDq1KnK/mPHjtHQ0LDiOUeOHKGnpweAM2fOVM7v6emhoaGhantgYKByzv3Hlt4Y2C4cOXKk6hGLxTh58iQnT56kt7e38hosx8iTJ09y5MgRTp8+zf/0P/1PlWucOHGCn/3sZw98r7feeotYLMaxY8c4efIk/f39Dxzf19fHkSNHOHfuXNX+o0ePcunSJU6dOlX5+xQIBJtIPgqmAaa+2TN5YgYGBqrWh6Xr0IPWjZ6eHs6fP79Z0xYIngo1HQECKgXJAgFQHQECmP8Edm1MataDCIfDXLhwgcHBwUpNSbkYH6gI96X/fu9PgfrOd75TEfTlhWylNK5YLMb3vve9bf3/QrkX2P0srQfs6+ujr6+PWCxGd3f3slrBlWoHu7u7Kym2wWCwyr1tpfFrHVset5oj3Fa29BYI6pp8KXXV1Knln0ixWIzXX3+9EtmOxWL81V/9VdXx7b5uCLY3NR0BEgiWUcxWb2upTZlGLBarqim5v7Hl8ePHV/wBXo5QCAQCgWATyJUEkFHbEaBLly4BVBpEB4NB0SReIFiCEECC+uL+CJBR3JRplJ0Jl7K0PuTYsWPLIj73R4QEG8P90RmB4Ekpt1lYagcuqFGqIkC1S1n4lN0lBYJa5u233+aNN95Y1zYL20YAiQVqm7BFBNBSYrEYZ8+ercq/BisKVE5J6O/vryxYj8OpU6eq6l+WWi0LBE+DjVigaoVymwXRW64OKEeATGNz57EODA0NAZYIamhoWFbXI9YNQa3wrW99ix/96Ed0dHSs2zVrN8H1ERF9gLYJ9wugLXAX78SJE7z77rvL+rucPHmSEydO8NZbb3H27NlKytzjsLRBp0CwGWxEozqB4KlTJxEgqO4zdv78eU6cOMHg4GBlLRLrhmA7s20EkGCbsMUiQCdOnOD06dPL0uFgsWfMwMAAAwMDlUaaAoFg44nFYpWeTyvR399fsTYPh8PCYXQ7YBqWBTbUhQBayvHjx+nt7WVgYKDmmi0LBBvBtkmBE2wT7jdBMDdPAJ04cYLvfve7K4qfMuUokDA/eHLOnTvHsWPHKukcwlZasBLnz5/n1KlTvP766w8cd/r0ad56663KHXJhC7wNKFtgQ80LoP7+fs6cOVOpMy0L+idJtRYI6gkhgAT1xRaJAJ0/f57z58/z+uuv09DQQENDw4q9ed566y26u7uFO88TcurUKS5cuMCFCxe4fPkyP/vZzx54d/9xKFuRC2qb48ePP/SGQ39/f9W/n97eXt55550Nnplg0ynX/0DNu8AdPXqUhYUF9uzZQ0NDA6dOndqQ70WBoFYRKXCC+mKLCKDjx4+v2pvm/h4vS7e7u7sZHBxcdk5vb++q11tp/HYiFotx7ty5qh5JwWBw3Rf61Xr2COqPgYGBKtfGUCgkIorbgfwSAVTjEaBgMMjp06c5ffr0iscftG5s9zVFsD0QAkhQPxhF0AvV+zYxBW47kdFz3EqPbci1D3h24lacqx6/dOnSqjntsViMb3/
      "text/plain": [
       "<Figure size 830x539.5 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "with plt.style.context('science'):\n",
    "    plt.rcParams.update({'font.size': 9})\n",
    "\n",
    "    figwidth = 8.3\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(figwidth, 0.65 * figwidth))\n",
    "\n",
    "    for n, cat in enumerate(cats):\n",
    "        i, j = n // 2, n % 2\n",
    "        ax = axs[i, j]\n",
    "\n",
    "        for sim in sims:\n",
    "            x = kfound[f\"{sim}_{cat}\"]\n",
    "            y = [data_mean[f\"{sim}_{cat}_{ksim}\"] for ksim in x]\n",
    "            yerr = [data_std[f\"{sim}_{cat}_{ksim}\"] for ksim in x]\n",
    "            x = np.array(x) * 0.674\n",
    "\n",
    "            ax.plot(x, y, label=simname_to_pretty(sim))\n",
    "            ax.fill_between(x, np.array(y) - np.array(yerr), np.array(y) + np.array(yerr), alpha=0.5)\n",
    "\n",
    "        ax.text(0.5, 0.9, catalogue_to_pretty(cat),\n",
    "            transform=ax.transAxes, #fontsize=\"small\",\n",
    "            verticalalignment='center', horizontalalignment='center',\n",
    "            bbox=dict(facecolor='white', alpha=0.5),\n",
    "            )\n",
    "\n",
    "        if n == 0:\n",
    "            ax.legend(fontsize=\"small\", loc='upper right')\n",
    "\n",
    "        ax.set_ylabel(key_label)\n",
    "        ax.set_xlabel(r\"$R_{\\rm offset} ~ [\\mathrm{Mpc} / h]$\")\n",
    "        ax.set_xlim(0),\n",
    "\n",
    "    fig.tight_layout()\n",
    "    fname = f\"../../plots/void_{key}_per_observer.png\"\n",
    "    if no_Vext:\n",
    "        fname = fname.replace(\".png\", \"_no_Vext.png\")\n",
    "    print(f\"Saving to `{fname}`.\")\n",
    "    fig.savefig(fname, dpi=450)\n",
    "    fig.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": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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
}