From 35ccfb5c673b73eef8c8a9067a377b702965b691 Mon Sep 17 00:00:00 2001 From: Richard Stiskalek Date: Sat, 17 Jun 2023 19:52:26 +0100 Subject: [PATCH] New matches (#69) * Remove old file * Add velocity plotting * add smooth scale * Fix bug * Improve paths * Edit plotting * Add smoothed density * Update boundaries * Add basics * Further docs * Remove blank * Better catalog broadcasting * Update high res size * Update plotting routines * Update routine * Update plotting * Fix field saving name * Add better colormap for environemnt --- csiborgtools/field/density.py | 9 +- csiborgtools/read/paths.py | 9 +- scripts/field_fromparts.py | 76 ------------- scripts/field_prop.py | 198 +++++++++++++++++++++++++++++----- scripts/match_finsnap.py | 1 - scripts/match_finsnap.yml | 8 -- scripts/utils.py | 23 +++- scripts_plots/plot_data.py | 97 ++++++++++++----- scripts_plots/plot_match.py | 91 ++++++++++++---- 9 files changed, 343 insertions(+), 169 deletions(-) delete mode 100644 scripts/field_fromparts.py diff --git a/csiborgtools/field/density.py b/csiborgtools/field/density.py index 78efc1c..1039766 100644 --- a/csiborgtools/field/density.py +++ b/csiborgtools/field/density.py @@ -305,7 +305,7 @@ class VelocityField(BaseField): vel *= mass.reshape(-1, 1) / mpart for i in range(3): - MASL.MA(pos, rho_vel[i], self.boxsize, self.MAS, W=vel[i, :], + MASL.MA(pos, rho_vel[i], self.boxsize, self.MAS, W=vel[:, i], verbose=False) if end == nparts: break @@ -417,11 +417,8 @@ class TidalTensorField(BaseField): Returns ------- environment : 3-dimensional array of shape `(grid, grid, grid)` - The environment of each grid cell. Possible values are: - - 0: void - - 1: sheet - - 2: filament - - 3: knot + The environment of each grid cell. Possible values are 0 (void), + 1 (sheet), 2 (filament), 3 (knot). """ environment = numpy.full(eigvals.shape[:-1], numpy.nan, dtype=numpy.float32) diff --git a/csiborgtools/read/paths.py b/csiborgtools/read/paths.py index f032c5d..51d9f6b 100644 --- a/csiborgtools/read/paths.py +++ b/csiborgtools/read/paths.py @@ -356,8 +356,8 @@ class Paths: fname = f"parts_{str(nsim).zfill(5)}.h5" return join(fdir, fname) - def field(self, kind, MAS, grid, nsim, in_rsp): - """ + def field(self, kind, MAS, grid, nsim, in_rsp, smooth_scale=None): + r""" Path to the files containing the calculated density fields in CSiBORG. Parameters @@ -373,6 +373,8 @@ class Paths: IC realisation index. in_rsp : bool Whether the calculation is performed in redshift space. + smooth_scale : float + Smoothing scale in :math:`\mathrm{Mpc}/h` Returns ------- @@ -387,6 +389,9 @@ class Paths: if in_rsp: kind = kind + "_rsp" fname = f"{kind}_{MAS}_{str(nsim).zfill(5)}_grid{grid}.npy" + if smooth_scale is not None and smooth_scale > 0: + smooth_scale = float(smooth_scale) + fname = fname.replace(".npy", f"smooth{smooth_scale}.npy") return join(fdir, fname) def halo_counts(self, simname, nsim): diff --git a/scripts/field_fromparts.py b/scripts/field_fromparts.py deleted file mode 100644 index 02ea934..0000000 --- a/scripts/field_fromparts.py +++ /dev/null @@ -1,76 +0,0 @@ -# Copyright (C) 2022 Richard Stiskalek -# This program is free software; you can redistribute it and/or modify it -# under the terms of the GNU General Public License as published by the -# Free Software Foundation; either version 3 of the License, or (at your -# option) any later version. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General -# Public License for more details. -# -# You should have received a copy of the GNU General Public License along -# with this program; if not, write to the Free Software Foundation, Inc., -# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -""" -MPI script to calculate the density fields on CSiBORG simulations in the final -snapshot. -""" -from argparse import ArgumentParser -from datetime import datetime -from distutils.util import strtobool - -import numpy -from mpi4py import MPI - -try: - import csiborgtools -except ModuleNotFoundError: - import sys - sys.path.append("../") - import csiborgtools - -comm = MPI.COMM_WORLD -rank = comm.Get_rank() -nproc = comm.Get_size() -verbose = nproc == 1 - -parser = ArgumentParser() -parser.add_argument("--ics", type=int, nargs="+", default=None, - help="IC realisations. If `-1` processes all simulations.") -parser.add_argument("--kind", type=str, choices=["density", "velocity"], - help="Calculate the density or velocity field?") -parser.add_argument("--MAS", type=str, choices=["NGP", "CIC", "TSC", "PCS"], - help="Mass assignment scheme.") -parser.add_argument("--grid", type=int, help="Grid resolution.") -parser.add_argument("--in_rsp", type=lambda x: bool(strtobool(x)), - help="Calculate the density field in redshift space?") -args = parser.parse_args() -paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) -mpart = 1.1641532e-10 # Particle mass in CSiBORG simulations. - -if args.ics is None or args.ics[0] == -1: - ics = paths.get_ics("csiborg") -else: - ics = args.ics - - -for i in csiborgtools.fits.split_jobs(len(ics), nproc)[rank]: - nsim = ics[i] - print(f"{datetime.now()}: rank {rank} working on simulation {nsim}.", - flush=True) - - nsnap = max(paths.get_snapshots(nsim)) - box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) - parts = csiborgtools.read.read_h5(paths.particles(nsim))["particles"] - - if args.kind == "density": - gen = csiborgtools.field.DensityField(box, args.MAS) - field = gen(parts, args.grid, in_rsp=args.in_rsp, verbose=verbose) - else: - gen = csiborgtools.field.VelocityField(box, args.MAS) - field = gen(parts, args.grid, mpart, verbose=verbose) - - fout = paths.field(args.kind, args.MAS, args.grid, nsim, args.in_rsp) - print(f"{datetime.now()}: rank {rank} saving output to `{fout}`.") - numpy.save(fout, field) diff --git a/scripts/field_prop.py b/scripts/field_prop.py index 7cf73df..a49a078 100644 --- a/scripts/field_prop.py +++ b/scripts/field_prop.py @@ -40,7 +40,23 @@ from utils import get_nsims ############################################################################### -def density_field(nsim, parser_args): +def density_field(nsim, parser_args, to_save=True): + """ + Calculate the density field in the CSiBORG simulation. + + Parameters + ---------- + nsim : int + Simulation index. + parser_args : argparse.Namespace + Parsed arguments. + to_save : bool, optional + Whether to save the output to disk. + + Returns + ------- + field : 3-dimensional array + """ paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) nsnap = max(paths.get_snapshots(nsim)) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) @@ -50,10 +66,47 @@ def density_field(nsim, parser_args): field = gen(parts, parser_args.grid, in_rsp=parser_args.in_rsp, verbose=parser_args.verbose) - fout = paths.field("density", parser_args.MAS, parser_args.grid, - nsim, parser_args.in_rsp) - print(f"{datetime.now()}: saving output to `{fout}`.") - numpy.save(fout, field) + if to_save: + fout = paths.field("density", parser_args.MAS, parser_args.grid, + nsim, parser_args.in_rsp) + print(f"{datetime.now()}: saving output to `{fout}`.") + numpy.save(fout, field) + return field + + +def density_field_smoothed(nsim, parser_args, to_save=True): + """ + Calculate the smoothed density field in the CSiBORG simulation. The + unsmoothed density field must already be precomputed. + + Parameters + ---------- + nsim : int + Simulation index. + parser_args : argparse.Namespace + Parsed arguments. + to_save : bool, optional + Whether to save the output to disk. + + Returns + ------- + smoothed_density : 3-dimensional array + """ + paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) + nsnap = max(paths.get_snapshots(nsim)) + box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) + + # Load the real space overdensity field + rho = numpy.load(paths.field("density", parser_args.MAS, parser_args.grid, + nsim, in_rsp=False)) + rho = csiborgtools.field.smoothen_field(rho, parser_args.smooth_scale, + box.boxsize, threads=1) + if to_save: + fout = paths.field("density", parser_args.MAS, parser_args.grid, + nsim, parser_args.in_rsp, parser_args.smooth_scale) + print(f"{datetime.now()}: saving output to `{fout}`.") + numpy.save(fout, rho) + return rho ############################################################################### @@ -61,9 +114,28 @@ def density_field(nsim, parser_args): ############################################################################### -def velocity_field(nsim, parser_args): +def velocity_field(nsim, parser_args, to_save=True): + """ + Calculate the velocity field in the CSiBORG simulation. + + Parameters + ---------- + nsim : int + Simulation index. + parser_args : argparse.Namespace + Parsed arguments. + to_save : bool, optional + Whether to save the output to disk. + + Returns + ------- + velfield : 4-dimensional array + """ if parser_args.in_rsp: raise NotImplementedError("Velocity field in RSP is not implemented.") + if parser_args.smooth_scale > 0: + raise NotImplementedError( + "Smoothed velocity field is not implemented.") paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) mpart = 1.1641532e-10 # Particle mass in CSiBORG simulations. nsnap = max(paths.get_snapshots(nsim)) @@ -73,10 +145,12 @@ def velocity_field(nsim, parser_args): gen = csiborgtools.field.VelocityField(box, parser_args.MAS) field = gen(parts, parser_args.grid, mpart, verbose=parser_args.verbose) - fout = paths.field("velocity", parser_args.MAS, parser_args.grid, - nsim, in_rsp=False) - print(f"{datetime.now()}: saving output to `{fout}`.") - numpy.save(fout, field) + if to_save: + fout = paths.field("velocity", parser_args.MAS, parser_args.grid, + nsim, in_rsp=False) + print(f"{datetime.now()}: saving output to `{fout}`.") + numpy.save(fout, field) + return field ############################################################################### @@ -84,7 +158,23 @@ def velocity_field(nsim, parser_args): ############################################################################### -def potential_field(nsim, parser_args): +def potential_field(nsim, parser_args, to_save=True): + """ + Calculate the potential field in the CSiBORG simulation. + + Parameters + ---------- + nsim : int + Simulation index. + parser_args : argparse.Namespace + Parsed arguments. + to_save : bool, optional + Whether to save the output to disk. + + Returns + ------- + potential : 3-dimensional array + """ paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) nsnap = max(paths.get_snapshots(nsim)) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) @@ -93,6 +183,9 @@ def potential_field(nsim, parser_args): density_gen = csiborgtools.field.DensityField(box, parser_args.MAS) rho = numpy.load(paths.field("density", parser_args.MAS, parser_args.grid, nsim, in_rsp=False)) + if parser_args.smooth_scale > 0: + rho = csiborgtools.field.smoothen_field(rho, parser_args.smooth_scale, + box.boxsize, threads=1) rho = density_gen.overdensity_field(rho) # Calculate the real space potentiel field gen = csiborgtools.field.PotentialField(box, parser_args.MAS) @@ -102,10 +195,12 @@ def potential_field(nsim, parser_args): parts = csiborgtools.read.read_h5(paths.particles(nsim))["particles"] field = csiborgtools.field.field2rsp(field, parts=parts, box=box, verbose=parser_args.verbose) - fout = paths.field(parser_args.kind, parser_args.MAS, parser_args.grid, - nsim, parser_args.in_rsp) - print(f"{datetime.now()}: saving output to `{fout}`.") - numpy.save(fout, field) + if to_save: + fout = paths.field(parser_args.kind, parser_args.MAS, parser_args.grid, + nsim, parser_args.in_rsp, parser_args.smooth_scale) + print(f"{datetime.now()}: saving output to `{fout}`.") + numpy.save(fout, field) + return field ############################################################################### @@ -113,9 +208,28 @@ def potential_field(nsim, parser_args): ############################################################################### -def radvel_field(nsim, parser_args): +def radvel_field(nsim, parser_args, to_save=True): + """ + Calculate the radial velocity field in the CSiBORG simulation. + + Parameters + ---------- + nsim : int + Simulation index. + parser_args : argparse.Namespace + Parsed arguments. + to_save : bool, optional + Whether to save the output to disk. + + Returns + ------- + radvel : 3-dimensional array + """ if parser_args.in_rsp: raise NotImplementedError("Radial vel. field in RSP not implemented.") + if parser_args.smooth_scale > 0: + raise NotImplementedError( + "Smoothed radial vel. field not implemented.") paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) nsnap = max(paths.get_snapshots(nsim)) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) @@ -124,11 +238,12 @@ def radvel_field(nsim, parser_args): nsim, parser_args.in_rsp)) gen = csiborgtools.field.VelocityField(box, parser_args.MAS) field = gen.radial_velocity(vel) - - fout = paths.field("radvel", parser_args.MAS, parser_args.grid, - nsim, parser_args.in_rsp) - print(f"{datetime.now()}: saving output to `{fout}`.") - numpy.save(fout, field) + if to_save: + fout = paths.field("radvel", parser_args.MAS, parser_args.grid, + nsim, parser_args.in_rsp) + print(f"{datetime.now()}: saving output to `{fout}`.") + numpy.save(fout, field) + return field ############################################################################### @@ -136,7 +251,23 @@ def radvel_field(nsim, parser_args): ############################################################################### -def environment_field(nsim, parser_args): +def environment_field(nsim, parser_args, to_save=True): + """ + Calculate the environmental classification in the CSiBORG simulation. + + Parameters + ---------- + nsim : int + Simulation index. + parser_args : argparse.Namespace + Parsed arguments. + to_save : bool, optional + Whether to save the output to disk. + + Returns + ------- + env : 3-dimensional array + """ if parser_args.in_rsp: raise NotImplementedError("Env. field in RSP not implemented.") paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) @@ -150,6 +281,9 @@ def environment_field(nsim, parser_args): print(f"{datetime.now()}: loading density field.") rho = numpy.load(paths.field("density", parser_args.MAS, parser_args.grid, nsim, in_rsp=False)) + if parser_args.smooth_scale > 0: + rho = csiborgtools.field.smoothen_field(rho, parser_args.smooth_scale, + box.boxsize, threads=1) rho = density_gen.overdensity_field(rho) # Calculate the real space tidal tensor field, delete overdensity. if parser_args.verbose: @@ -157,12 +291,16 @@ def environment_field(nsim, parser_args): tensor_field = gen(rho) del rho collect() + + # TODO: Optionally drag the field to RSP. + # Calculate the eigenvalues of the tidal tensor field, delete tensor field. if parser_args.verbose: print(f"{datetime.now()}: calculating eigenvalues.") eigvals = gen.tensor_field_eigvals(tensor_field) del tensor_field collect() + # Classify the environment based on the eigenvalues. if parser_args.verbose: print(f"{datetime.now()}: classifying environment.") @@ -170,10 +308,12 @@ def environment_field(nsim, parser_args): del eigvals collect() - fout = paths.field("environment", parser_args.MAS, parser_args.grid, - nsim, parser_args.in_rsp) - print(f"{datetime.now()}: saving output to `{fout}`.") - numpy.save(fout, env) + if to_save: + fout = paths.field("environment", parser_args.MAS, parser_args.grid, + nsim, parser_args.in_rsp, parser_args.smooth_scale) + print(f"{datetime.now()}: saving output to `{fout}`.") + numpy.save(fout, env) + return env ############################################################################### @@ -194,6 +334,7 @@ if __name__ == "__main__": parser.add_argument("--grid", type=int, help="Grid resolution.") parser.add_argument("--in_rsp", type=lambda x: bool(strtobool(x)), help="Calculate in RSP?") + parser.add_argument("--smooth_scale", type=float, default=0) parser.add_argument("--verbose", type=lambda x: bool(strtobool(x)), help="Verbosity flag for reading in particles.") parser_args = parser.parse_args() @@ -203,7 +344,10 @@ if __name__ == "__main__": def main(nsim): if parser_args.kind == "density": - density_field(nsim, parser_args) + if parser_args.smooth_scale > 0: + density_field_smoothed(nsim, parser_args) + else: + density_field(nsim, parser_args) elif parser_args.kind == "velocity": velocity_field(nsim, parser_args) elif parser_args.kind == "radvel": diff --git a/scripts/match_finsnap.py b/scripts/match_finsnap.py index 94a21df..4df0c4b 100644 --- a/scripts/match_finsnap.py +++ b/scripts/match_finsnap.py @@ -119,7 +119,6 @@ def collect_dist(args, paths): out = data["counts"] else: out += data["counts"] - remove(fname) fout = paths.cross_nearest(args.simname, args.run, "tot_counts", diff --git a/scripts/match_finsnap.yml b/scripts/match_finsnap.yml index 03edd89..5d9b7a6 100644 --- a/scripts/match_finsnap.yml +++ b/scripts/match_finsnap.yml @@ -19,7 +19,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 12.4 - max: 12.8 islog: true "mass002": @@ -28,7 +27,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 12.6 - max: 13.0 islog: true "mass003": @@ -37,7 +35,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 12.8 - max: 13.2 islog: true "mass004": @@ -46,7 +43,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 13.0 - max: 13.4 islog: true "mass005": @@ -55,7 +51,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 13.2 - max: 13.6 islog: true "mass006": @@ -64,7 +59,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 13.4 - max: 13.8 islog: true "mass007": @@ -73,7 +67,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 13.6 - max: 14.0 islog: true "mass008": @@ -82,7 +75,6 @@ nbins_marks: 10 - totpartmass - group_mass min: 13.8 - max: 14.2 islog: true "mass009": diff --git a/scripts/utils.py b/scripts/utils.py index ee2f186..8c9d88f 100644 --- a/scripts/utils.py +++ b/scripts/utils.py @@ -165,6 +165,8 @@ def open_catalogues(args, config, paths, comm): if args.verbose and rank == 0: print(f"{datetime.now()}: opening catalogues.", flush=True) + # We first load all catalogues on the zeroth rank and broadcast their + # names. if rank == 0: cats = {} if args.simname == "csiborg": @@ -182,12 +184,27 @@ def open_catalogues(args, config, paths, comm): name = paths.quijote_fiducial_nsim(nsim, nobs) cat = ref_cat.pick_fiducial_observer(nobs, rmax=args.Rmax) cats.update({name: cat}) - + names = list(cats.keys()) if nproc > 1: for i in range(1, nproc): - comm.send(cats, dest=i, tag=nproc + i) + comm.send(names, dest=i, tag=nproc + i) else: - cats = comm.recv(source=0, tag=nproc + rank) + names = comm.recv(source=0, tag=nproc + rank) + + comm.Barrier() + # We then broadcast the catalogues to all ranks, one-by-one as MPI can + # only pass messages smaller than 2GB. + if nproc == 1: + return cats + + if rank > 0: + cats = {} + for name in names: + if rank == 0: + for i in range(1, nproc): + comm.send(cats[name], dest=i, tag=nproc + i) + else: + cats.update({name: comm.recv(source=0, tag=nproc + rank)}) return cats diff --git a/scripts_plots/plot_data.py b/scripts_plots/plot_data.py index 3d8c9b7..c06ad23 100644 --- a/scripts_plots/plot_data.py +++ b/scripts_plots/plot_data.py @@ -16,6 +16,7 @@ from os.path import join from argparse import ArgumentParser +import matplotlib as mpl import matplotlib.pyplot as plt import numpy import healpy @@ -209,8 +210,8 @@ def plot_hmf(pdf=False): plt.close() -def load_field(kind, nsim, grid, MAS, in_rsp=False): - """ +def load_field(kind, nsim, grid, MAS, in_rsp=False, smooth_scale=None): + r""" Load a single field. Parameters @@ -225,13 +226,16 @@ def load_field(kind, nsim, grid, MAS, in_rsp=False): Mass assignment scheme. in_rsp : bool, optional Whether to load the field in redshift space. + smooth_scale : float, optional + Smoothing scale in :math:`\mathrm{Mpc} / h`. Returns ------- field : n-dimensional array """ paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) - return numpy.load(paths.field(kind, MAS, grid, nsim, in_rsp=in_rsp)) + return numpy.load(paths.field(kind, MAS, grid, nsim, in_rsp=in_rsp, + smooth_scale=smooth_scale)) ############################################################################### @@ -239,9 +243,9 @@ def load_field(kind, nsim, grid, MAS, in_rsp=False): ############################################################################### -def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS", +def plot_projected_field(kind, nsim, grid, in_rsp, smooth_scale, MAS="PCS", highres_only=True, slice_find=None, pdf=False): - """ + r""" Plot the mean projected field, however can also plot a single slice. Parameters @@ -254,6 +258,8 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS", Grid size. in_rsp : bool Whether to load the field in redshift space. + smooth_scale : float + Smoothing scale in :math:`\mathrm{Mpc} / h`. MAS : str, optional Mass assignment scheme. highres_only : bool, optional @@ -273,11 +279,16 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS", box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) if kind == "overdensity": - field = load_field("density", nsim, grid, MAS=MAS, in_rsp=in_rsp) + field = load_field("density", nsim, grid, MAS=MAS, in_rsp=in_rsp, + smooth_scale=smooth_scale) density_gen = csiborgtools.field.DensityField(box, MAS) - field = density_gen.overdensity_field(field) + 2 + field = density_gen.overdensity_field(field) + 1 else: - field = load_field(kind, nsim, grid, MAS=MAS, in_rsp=in_rsp) + field = load_field(kind, nsim, grid, MAS=MAS, in_rsp=in_rsp, + smooth_scale=smooth_scale) + + if kind == "velocity": + field = field[0, ...] if highres_only: csiborgtools.field.fill_outside(field, numpy.nan, rmax=155.5, @@ -286,10 +297,11 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS", end = round(field.shape[0] * 0.73) field = field[start:end, start:end, start:end] - if kind != "environment": - cmap = "viridis" + if kind == "environment": + cmap = mpl.colors.ListedColormap( + ['red', 'lightcoral', 'limegreen', 'khaki']) else: - cmap = "brg" + cmap = "viridis" labels = [r"$y-z$", r"$x-z$", r"$x-y$"] with plt.style.context(plt_utils.mplstyle): @@ -309,12 +321,15 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS", else: ax[i].imshow(img, vmin=vmin, vmax=vmax, cmap=cmap) - if not highres_only: + frad = 155.5 / 677.7 + if not highres_only and 0.5 - frad < slice_find < 0.5 + frad: theta = numpy.linspace(0, 2 * numpy.pi, 100) - rad = 155.5 / 677.7 * grid + z = (slice_find - 0.5) * grid + R = 155.5 / 677.7 * grid + rad = R * numpy.sqrt(1 - z**2 / R**2) ax[i].plot(rad * numpy.cos(theta) + grid // 2, rad * numpy.sin(theta) + grid // 2, - lw=plt.rcParams["lines.linewidth"], zorder=1, + lw=0.75 * plt.rcParams["lines.linewidth"], zorder=1, c="red", ls="--") ax[i].set_title(labels[i]) @@ -343,17 +358,32 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS", (xticks * size / ncells - size / 2).astype(int)) ax[i].set_xlabel(r"$x_j ~ [\mathrm{Mpc} / h]$") - cbar_ax = fig.add_axes([1.0, 0.1, 0.025, 0.8]) + cbar_ax = fig.add_axes([0.982, 0.155, 0.025, 0.75], + transform=ax[2].transAxes) if slice_find is None: clabel = "Mean projected field" else: clabel = "Sliced field" - fig.colorbar(im, cax=cbar_ax, label=clabel) + + if kind == "environment": + bounds = [0, 1, 2, 3, 4] + norm = mpl.colors.BoundaryNorm(bounds, cmap.N) + cbar = fig.colorbar( + mpl.cm.ScalarMappable(cmap=cmap, norm=norm), cax=cbar_ax, + ticks=[0.5, 1.5, 2.5, 3.5]) + cbar.ax.set_yticklabels(["knot", "filament", "sheet", "void"], + rotation=90, va="center") + else: + fig.colorbar(im, cax=cbar_ax, label=clabel) fig.tight_layout(h_pad=0, w_pad=0) for ext in ["png"] if pdf is False else ["png", "pdf"]: - fout = join(plt_utils.fout, - f"field_{kind}_{nsim}_rsp{in_rsp}.{ext}") + fout = join( + plt_utils.fout, + f"field_{kind}_{nsim}_rsp{in_rsp}_hres{highres_only}.{ext}") + if smooth_scale is not None and smooth_scale > 0: + smooth_scale = float(smooth_scale) + fout = fout.replace(f".{ext}", f"_smooth{smooth_scale}.{ext}") print(f"Saving to `{fout}`.") fig.savefig(fout, dpi=plt_utils.dpi, bbox_inches="tight") plt.close() @@ -404,8 +434,8 @@ def get_sky_label(kind, volume_weight): return label -def plot_sky_distribution(kind, nsim, grid, nside, MAS="PCS", plot_groups=True, - dmin=0, dmax=220, plot_halos=None, +def plot_sky_distribution(kind, nsim, grid, nside, smooth_scale, MAS="PCS", + plot_groups=True, dmin=0, dmax=220, plot_halos=None, volume_weight=True, pdf=False): r""" Plot the sky distribution of a given field kind on the sky along with halos @@ -425,6 +455,8 @@ def plot_sky_distribution(kind, nsim, grid, nside, MAS="PCS", plot_groups=True, Grid size. nside : int Healpix nside of the sky projection. + smooth_scale : float + Smoothing scale in :math:`\mathrm{Mpc} / h`. MAS : str, optional Mass assignment scheme. plot_groups : bool, optional @@ -445,11 +477,13 @@ def plot_sky_distribution(kind, nsim, grid, nside, MAS="PCS", plot_groups=True, box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) if kind == "overdensity": - field = load_field("density", nsim, grid, MAS=MAS, in_rsp=False) + field = load_field("density", nsim, grid, MAS=MAS, in_rsp=False, + smooth_scale=smooth_scale) density_gen = csiborgtools.field.DensityField(box, MAS) - field = density_gen.overdensity_field(field) + 2 + field = density_gen.overdensity_field(field) + 1 else: - field = load_field(kind, nsim, grid, MAS=MAS, in_rsp=False) + field = load_field(kind, nsim, grid, MAS=MAS, in_rsp=False, + smooth_scale=smooth_scale) angpos = csiborgtools.field.nside2radec(nside) dist = numpy.linspace(dmin, dmax, 500) @@ -519,7 +553,22 @@ if __name__ == "__main__": if True: kind = "environment" grid = 256 + smooth_scale = 0 # plot_projected_field("overdensity", 7444, grid, in_rsp=True, # highres_only=False) plot_projected_field(kind, 7444, grid, in_rsp=False, - slice_find=0.5, highres_only=False) + smooth_scale=smooth_scale, slice_find=0.5, + highres_only=False) + + if False: + paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) + + d = csiborgtools.read.read_h5(paths.particles(7444))["particles"] + + plt.figure() + plt.hist(d[:100000, 4], bins="auto") + + plt.tight_layout() + plt.savefig("../plots/velocity_distribution.png", dpi=450, + bbox_inches="tight") + diff --git a/scripts_plots/plot_match.py b/scripts_plots/plot_match.py index a8e92c9..b2d85ef 100644 --- a/scripts_plots/plot_match.py +++ b/scripts_plots/plot_match.py @@ -648,7 +648,32 @@ def plot_significance(simname, runs, nsim, nobs, kind, kwargs, runs_to_mass): plt.close() -def plot_significance_vs_mass(simname, runs, nsim, nobs, kind, kwargs): +def make_binlims(run, runs_to_mass): + """ + Make bin limits for the 1NN distance runs, corresponding to the first half + of the log-mass bin. + + Parameters + ---------- + run : str + Run name. + runs_to_mass : dict + Dictionary mapping run names to total halo mass range. + + Returns + ------- + xmin, xmax : floats + """ + xmin, xmax = runs_to_mass[run] + xmax = xmin + (xmax - xmin) / 2 + xmin, xmax = 10**xmin, 10**xmax + if run == "mass009": + xmax = numpy.infty + return xmin, xmax + + +def plot_significance_vs_mass(simname, runs, nsim, nobs, kind, kwargs, + runs_to_mass): """ Plot significance of the 1NN distance as a function of the total mass. @@ -667,6 +692,8 @@ def plot_significance_vs_mass(simname, runs, nsim, nobs, kind, kwargs): (Kolmogorov-Smirnov p-value). kwargs : dict Nearest neighbour reader keyword arguments. + runs_to_mass : dict + Dictionary mapping run names to total halo mass range. Returns ------- @@ -686,8 +713,12 @@ def plot_significance_vs_mass(simname, runs, nsim, nobs, kind, kwargs): y = make_kl(simname, run, nsim, nobs, kwargs) else: y = numpy.log10(make_ks(simname, run, nsim, nobs, kwargs)) - xs.append(x) - ys.append(y) + + xmin, xmax = make_binlims(run, runs_to_mass) + mask = (x >= xmin) & (x < xmax) + xs.append(x[mask]) + ys.append(y[mask]) + xs = numpy.concatenate(xs) ys = numpy.concatenate(ys) @@ -715,7 +746,7 @@ def plot_significance_vs_mass(simname, runs, nsim, nobs, kind, kwargs): plt.close() -def plot_kl_vs_ks(simname, runs, nsim, nobs, kwargs): +def plot_kl_vs_ks(simname, runs, nsim, nobs, kwargs, runs_to_mass): """ Plot Kullback-Leibler divergence vs Kolmogorov-Smirnov statistic p-value. @@ -731,6 +762,8 @@ def plot_kl_vs_ks(simname, runs, nsim, nobs, kwargs): Fiducial Quijote observer index. For CSiBORG must be set to `None`. kwargs : dict Nearest neighbour reader keyword arguments. + runs_to_mass : dict + Dictionary mapping run names to total halo mass range. Returns ------- @@ -741,9 +774,16 @@ def plot_kl_vs_ks(simname, runs, nsim, nobs, kwargs): xs, ys, cs = [], [], [] for run in runs: - cs.append(reader.read_single(simname, run, nsim, nobs)["mass"]) - xs.append(make_kl(simname, run, nsim, nobs, kwargs)) - ys.append(make_ks(simname, run, nsim, nobs, kwargs)) + c = reader.read_single(simname, run, nsim, nobs)["mass"] + x = make_kl(simname, run, nsim, nobs, kwargs) + y = make_ks(simname, run, nsim, nobs, kwargs) + + cmin, cmax = make_binlims(run, runs_to_mass) + mask = (c >= cmin) & (c < cmax) + xs.append(x[mask]) + ys.append(y[mask]) + cs.append(c[mask]) + xs = numpy.concatenate(xs) ys = numpy.log10(numpy.concatenate(ys)) cs = numpy.log10(numpy.concatenate(cs)) @@ -768,7 +808,7 @@ def plot_kl_vs_ks(simname, runs, nsim, nobs, kwargs): plt.close() -def plot_kl_vs_overlap(runs, nsim, kwargs): +def plot_kl_vs_overlap(runs, nsim, kwargs, runs_to_mass): """ Plot KL divergence vs overlap for CSiBORG. @@ -780,6 +820,8 @@ def plot_kl_vs_overlap(runs, nsim, kwargs): Simulation index. kwargs : dict Nearest neighbour reader keyword arguments. + runs_to_mass : dict + Dictionary mapping run names to total halo mass range. Returns ------- @@ -802,12 +844,15 @@ def plot_kl_vs_overlap(runs, nsim, kwargs): prob_nomatch = prob_nomatch[mask] mass = mass[mask] - kl = make_kl("csiborg", run, nsim, nobs=None, kwargs=kwargs) - - xs.append(kl) - ys1.append(1 - numpy.mean(prob_nomatch, axis=1)) - ys2.append(numpy.std(prob_nomatch, axis=1)) - cs.append(numpy.log10(mass)) + x = make_kl("csiborg", run, nsim, nobs=None, kwargs=kwargs) + y1 = 1 - numpy.mean(prob_nomatch, axis=1) + y2 = numpy.std(prob_nomatch, axis=1) + cmin, cmax = make_binlims(run, runs_to_mass) + mask = (mass >= cmin) & (mass < cmax) + xs.append(x[mask]) + ys1.append(y1[mask]) + ys2.append(y2[mask]) + cs.append(numpy.log10(mass[mask])) xs = numpy.concatenate(xs) ys1 = numpy.concatenate(ys1) @@ -877,7 +922,7 @@ if __name__ == "__main__": delete_disk_caches_for_function(func) # Plot 1NN distance distributions. - if False: + if True: for i in range(1, 10): run = f"mass00{i}" for pulled_cdf in [True, False]: @@ -886,12 +931,12 @@ if __name__ == "__main__": plot_dist(run, "pdf", neighbour_kwargs, runs_to_mass) # Plot 1NN CDF differences. - if False: + if True: runs = [f"mass00{i}" for i in range(1, 10)] for pulled_cdf in [True, False]: plot_cdf_diff(runs, neighbour_kwargs, pulled_cdf=pulled_cdf, runs_to_mass=runs_to_mass) - if False: + if True: runs = [f"mass00{i}" for i in range(1, 9)] for kind in ["kl", "ks"]: plot_significance("csiborg", runs, 7444, nobs=None, kind=kind, @@ -902,12 +947,14 @@ if __name__ == "__main__": runs = [f"mass00{i}" for i in range(1, 10)] for kind in ["kl", "ks"]: plot_significance_vs_mass("csiborg", runs, 7444, nobs=None, - kind=kind, kwargs=neighbour_kwargs) + kind=kind, kwargs=neighbour_kwargs, + runs_to_mass=runs_to_mass) - if False: + if True: runs = [f"mass00{i}" for i in range(1, 10)] - plot_kl_vs_ks("csiborg", runs, 7444, None, kwargs=neighbour_kwargs) + plot_kl_vs_ks("csiborg", runs, 7444, None, kwargs=neighbour_kwargs, + runs_to_mass=runs_to_mass) - if False: + if True: runs = [f"mass00{i}" for i in range(1, 10)] - plot_kl_vs_overlap(runs, 7444, neighbour_kwargs) + plot_kl_vs_overlap(runs, 7444, neighbour_kwargs, runs_to_mass)