diff --git a/LRG/lrg_cutsky_propagator.py b/LRG/lrg_cutsky_propagator.py index e4689c5..c69cbe3 100644 --- a/LRG/lrg_cutsky_propagator.py +++ b/LRG/lrg_cutsky_propagator.py @@ -1,32 +1,55 @@ -# Calculate the propagator defined as P_cross/(bias*P_init). -# We firstly generate meshes for the IC and the tracer, then we use MeshFFTCorrelator to calculate the correlator. -# The propagator information can be obtained from the correlator. +################################################################################################################## +# 1. Calculate the propagator defined as P_cross/(bias*P_init), given the galaxy mock and IC. +# 2. We firstly generate meshes for the IC and tracer, then we use MeshFFTCorrelator to calculate the correlator. +# *** Note that for the MultiGrid recon, it is better to set nmesh = a*2^n, e.g. 3*2^9, to avoid any weird performance. *** +# 3. The propagator information can be obtained from the correlator. You can check show_propagator.ipynb as an example. +# (To run the code, you can check the script propagator_nzweight_perlmutter.sl.) # import os, sys -import fitsio +import fitsio, asdf import numpy as np +from scipy import interpolate from pathlib import Path from cosmoprimo.fiducial import DESI from pyrecon import utils from pyrecon.metrics import MeshFFTCorrelator, MeshFFTPropagator, CatalogMesh +from pypower.mesh import ArrayMesh from pypower import setup_logging, mpi, MeshFFTPower import argparse def fkp_weights(nz, P0=10000): return 1 / (1 + nz * P0) + +def interpolate_nz(z, nz, zmin, zmax): + zbins = np.linspace(zmin, zmax, 41) + nz_list = [] + for z0, z1 in zip(zbins[0:-1], zbins[1:]): + zmask = (z>z0)&(zzmin)&(Z_iczmin)&(Z_ic" + "
" ] }, "metadata": { @@ -314,50 +237,37 @@ } ], "source": [ - "fig, lax = plt.subplots(nrows=1, ncols=3, figsize=(16,4), dpi=150)\n", - "fig.subplots_adjust(wspace=0.3)\n", - "lax = lax.flatten()\n", - "for i, convention in enumerate(convention_list[0:2]):\n", - " if recon_algo == \"MultiGrid\":\n", - " fname_appendix = f\"{cap.upper()}_{zmin}z{zmax}_shift_{recon_algo}_randoms20X_reso7.8_smooth{smooth_radius}_pad1.5_{convention}_f{growth_rate:.3f}_b{bias:.2f}\"\n", - " else:\n", - " fname_appendix = f\"{cap.upper()}_{zmin}z{zmax}_shift_{recon_algo}_randoms20X_reso7.8_smooth{smooth_radius}_pad1.5_niter{niter}_{convention}_f{growth_rate:.3f}_b{bias:.2f}\"\n", - " filename = f\"correlator_LRG_c000_ph000_{fname_appendix}_nmesh{Nmesh}.npy\"\n", - " correlator, transfer, propagator = read_correlator(data_dir, filename, 2.0)\n", - " \n", - " for imu, mu in enumerate(propagator.muavg[mu_id_list]):\n", - " lax[0].plot(k, correlator(k=k, mu=mu), color=f'C{imu}', ls=line_type_list[i], label=r'{0}, $\\mu = {1:.2f}$'.format(convention, mu))\n", - " lax[1].plot(k, transfer(k=k, mu=mu), color=f'C{imu}', ls=line_type_list[i])\n", - " lax[2].plot(k, propagator(k=k, mu=mu), color=f'C{imu}', ls=line_type_list[i])\n", - " \n", - "# plot pre-recon result\n", - "for imu, mu in enumerate(np.abs(propagator.muavg[mu_id_list])):\n", - " #print(\"mu:\", mu, transfer_pre(k=k, mu=mu))\n", - " lax[0].plot(k, correlator_pre(k=k, mu=mu), color=f'C{imu}', ls=\":\", label=r'{0}, $\\mu = {1:.2f}$'.format(\"pre-recon\", mu))\n", - " lax[1].plot(k, transfer_pre(k=k, mu=mu), color=f'C{imu}', ls=\":\")\n", - " lax[2].plot(k, propagator_pre(k=k, mu=mu), color=f'C{imu}', ls=\":\")\n", - " \n", - "for ax in lax:\n", - " #ax.legend(fontsize=14)\n", - " ax.grid(True)\n", - " ax.set_xscale(\"log\")\n", - " ax.set_ylim([0., 1.5])\n", - " ax.set_xlabel(r'$k$ [$\\mathrm{Mpc}/h$]', fontsize=14)\n", - " \n", - "lax[0].legend(fontsize=12)\n", - "lax[0].set_ylabel(r'$r(k) = P_{\\mathrm{rec},\\mathrm{init}}/\\sqrt{P_{\\mathrm{rec}}P_{\\mathrm{init}}}$', fontsize=14)\n", - "lax[1].set_ylabel(r'$t(k) = \\sqrt{P_{\\mathrm{rec}}/P_{\\mathrm{init}}}$', fontsize=14)\n", - "lax[2].set_ylabel(r'$g(k) = P_{\\mathrm{rec},\\mathrm{init}}/P_{\\mathrm{init}}$', fontsize=14)\n", - "plt.suptitle(recon_algo, fontsize=16)\n", - "plt.tight_layout()\n", - "# ofile = f\"figs/propagator/{rec_method}_BGS_propagator_z0.2_nmesh512_Rs15.0_bias1.63.png\"\n", - "# plt.savefig(ofile)" + "fig, axes = plt.subplots(nrows=1, ncols=2, sharey=True, figsize=(12,6), dpi=150)\n", + "plt.subplots_adjust(\n", + " left=0.15,\n", + " bottom=None,\n", + " right=None,\n", + " top=0.9,\n", + " wspace=0.0,\n", + " hspace=None,\n", + ")\n", + "for i, mu in enumerate(mu_list):\n", + " for j, smooth_radius in enumerate(sm_radius_list):\n", + " C_MG = propagator_list[j](k=k, mu=mu)\n", + " C_app = cal_propagator(k, mu, beta, smooth_radius, convention, Sigma_xy_list[j], Sigma_z_list[j])\n", + " axes[i].plot(k, C_MG, ls='-', color=f'C{j}', label=r\"$\\Sigma_{\\mathrm{sm}}=%.1f$, $\\Sigma_{\\perp}=%.2f$\"%(smooth_radius, Sigma_xy_list[j]) if i==0 else r\"$\\Sigma_{\\mathrm{sm}}=%.1f$, $\\Sigma_{\\parallel}=%.2f$\"%(smooth_radius, Sigma_z_list[j]))\n", + " axes[i].plot(k, C_app, ls=':', color=f'C{j}')\n", + " axes[i].set_xscale(\"log\")\n", + " axes[i].legend(loc=\"lower left\", fontsize=16)\n", + " axes[i].text(0.2, 1.1, r\"$\\mu=%.3f$\"%(mu) + f\"\\n{recon_algo}\", fontsize=16)\n", + " axes[i].set_xlabel(r\"$k\\;\\mathrm{[hMpc^{-1}]}$\", fontsize=16)\n", + "\n", + "\n", + "axes[0].set_ylabel(\"propagator\", fontsize=14)\n", + "fig.suptitle(r\"${0:.2f}