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

@@ -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:
+            The environment of each grid cell. Possible values are 0 (void),
-            - 0: void
+            1 (sheet), 2 (filament), 3 (knot).
-            - 1: sheet
-            - 2: filament
-            - 3: knot
         """
         environment = numpy.full(eigvals.shape[:-1], numpy.nan,
                                  dtype=numpy.float32)

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):

scripts/field_fromparts.py

@@ -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)

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,
+    if to_save:
-                       nsim, parser_args.in_rsp)
+        fout = paths.field("density", parser_args.MAS, parser_args.grid,
-    print(f"{datetime.now()}: saving output to `{fout}`.")
+                           nsim, parser_args.in_rsp)
-    numpy.save(fout, field)
+        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,
+    if to_save:
-                       nsim, in_rsp=False)
+        fout = paths.field("velocity", parser_args.MAS, parser_args.grid,
-    print(f"{datetime.now()}: saving output to `{fout}`.")
+                           nsim, in_rsp=False)
-    numpy.save(fout, field)
+        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,
+    if to_save:
-                       nsim, parser_args.in_rsp)
+        fout = paths.field(parser_args.kind, parser_args.MAS, parser_args.grid,
-    print(f"{datetime.now()}: saving output to `{fout}`.")
+                           nsim, parser_args.in_rsp, parser_args.smooth_scale)
-    numpy.save(fout, field)
+        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)
-
+    if to_save:
-    fout = paths.field("radvel", parser_args.MAS, parser_args.grid,
+        fout = paths.field("radvel", parser_args.MAS, parser_args.grid,
-                       nsim, parser_args.in_rsp)
+                           nsim, parser_args.in_rsp)
-    print(f"{datetime.now()}: saving output to `{fout}`.")
+        print(f"{datetime.now()}: saving output to `{fout}`.")
-    numpy.save(fout, field)
+        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,
+    if to_save:
-                       nsim, parser_args.in_rsp)
+        fout = paths.field("environment", parser_args.MAS, parser_args.grid,
-    print(f"{datetime.now()}: saving output to `{fout}`.")
+                           nsim, parser_args.in_rsp, parser_args.smooth_scale)
-    numpy.save(fout, env)
+        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":

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",

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":

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
 
 

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":
+    if kind == "environment":
-        cmap = "viridis"
+        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,
+            fout = join(
-                        f"field_{kind}_{nsim}_rsp{in_rsp}.{ext}")
+                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,
+def plot_sky_distribution(kind, nsim, grid, nside, smooth_scale, MAS="PCS",
-                          dmin=0, dmax=220, plot_halos=None,
+                          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")
+

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"])
+        c = reader.read_single(simname, run, nsim, nobs)["mass"]
-        xs.append(make_kl(simname, run, nsim, nobs, kwargs))
+        x = make_kl(simname, run, nsim, nobs, kwargs)
-        ys.append(make_ks(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)
+        x = make_kl("csiborg", run, nsim, nobs=None, kwargs=kwargs)
-
+        y1 = 1 - numpy.mean(prob_nomatch, axis=1)
-        xs.append(kl)
+        y2 = numpy.std(prob_nomatch, axis=1)
-        ys1.append(1 - numpy.mean(prob_nomatch, axis=1))
+        cmin, cmax = make_binlims(run, runs_to_mass)
-        ys2.append(numpy.std(prob_nomatch, axis=1))
+        mask = (mass >= cmin) & (mass < cmax)
-        cs.append(numpy.log10(mass))
+        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)