Add better diagnostics & plotting (#67)

* Add caching functions

* Add limts

* Add new mass runs

* Update .gitignore

* Edit which CDFs are loaded

* Stop saving cross hindxs

* Change dist to half precision

* New nearest path

* Add neighbour counting

* Add neighbour kwargs

* Update work in progress

* Add new counting

* Add how dist is built

* Collect dist to 1 file

* Update reading routine

* Delete Quijote files

* Remove file

* Back to float32

* Fix bugs

* Rename utils

* Remove neighbuor kwargs

* Rename file

* Fix bug

* Rename plt utils

* Change where nghb kwargs from

* line length

* Remove old notebooks

* Move survey

* Add white space

* Update TODO

* Update CDF calculation

* Update temporarily plotting

* Merge branch 'add_diagnostics' of github.com:Richard-Sti/csiborgtools into add_diagnostics

* Start adding documentation to plotting

* Remove comments

* Better code documentation

* Some work on tidal tensor

* Better plotting

* Add comment

* Remove nb

* Remove comment

* Add documentation

* Update plotting

* Update submission

* Update KL vs KS plots

* Update the plotting routine

* Update plotting

* Update plotting routines
This commit is contained in:
Richard Stiskalek 2023-06-16 14:33:27 +01:00 committed by GitHub
parent 004d9629a2
commit ccbbbd24b4
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
20 changed files with 1075 additions and 32121 deletions

1
.gitignore vendored
View file

@ -22,3 +22,4 @@ scripts_test/
scripts_plots/python.sh scripts_plots/python.sh
scripts_plots/submit.sh scripts_plots/submit.sh
scripts_plots/*.out scripts_plots/*.out
scripts_plots/*.sh

View file

@ -19,3 +19,26 @@ paths_glamdring = {"srcdir": "/mnt/extraspace/hdesmond/",
"postdir": "/mnt/extraspace/rstiskalek/CSiBORG/", "postdir": "/mnt/extraspace/rstiskalek/CSiBORG/",
"quijote_dir": "/mnt/extraspace/rstiskalek/Quijote", "quijote_dir": "/mnt/extraspace/rstiskalek/Quijote",
} }
neighbour_kwargs = {"rmax_radial": 155 / 0.705,
"nbins_radial": 50,
"rmax_neighbour": 100.,
"nbins_neighbour": 150,
"paths_kind": paths_glamdring}
###############################################################################
# Surveys #
###############################################################################
class SDSS:
@staticmethod
def steps(cls):
return [(lambda x: cls[x], ("IN_DR7_LSS",)),
(lambda x: cls[x] < 17.6, ("ELPETRO_APPMAG_r", )),
(lambda x: cls[x] < 155, ("DIST", ))
]
def __call__(self):
return read.SDSS(h=1, sel_steps=self.steps)

View file

@ -14,9 +14,6 @@
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
""" """
Density field and cross-correlation calculations. Density field and cross-correlation calculations.
TODO:
- [ ] Project the velocity field along the line of sight.
""" """
from abc import ABC from abc import ABC
@ -370,9 +367,9 @@ class TidalTensorField(BaseField):
box : :py:class:`csiborgtools.read.CSiBORGBox` box : :py:class:`csiborgtools.read.CSiBORGBox`
The simulation box information and transformations. The simulation box information and transformations.
MAS : str MAS : str
Mass assignment scheme. Options are Options are: 'NGP' (nearest grid Mass assignment scheme used to calculate the density field. Options
point), 'CIC' (cloud-in-cell), 'TSC' (triangular-shape cloud), 'PCS' are: 'NGP' (nearest grid point), 'CIC' (cloud-in-cell), 'TSC'
(piecewise cubic spline). (triangular-shape cloud), 'PCS' (piecewise cubic spline).
""" """
def __init__(self, box, MAS): def __init__(self, box, MAS):
self.box = box self.box = box
@ -384,8 +381,6 @@ class TidalTensorField(BaseField):
Calculate eigenvalues of the tidal tensor field, sorted in increasing Calculate eigenvalues of the tidal tensor field, sorted in increasing
order. order.
TODO: evaluate this on a grid instead.
Parameters Parameters
---------- ----------
tidal_tensor : :py:class:`MAS_library.tidal_tensor` tidal_tensor : :py:class:`MAS_library.tidal_tensor`
@ -396,20 +391,14 @@ class TidalTensorField(BaseField):
------- -------
eigvals : 3-dimensional array of shape `(grid, grid, grid)` eigvals : 3-dimensional array of shape `(grid, grid, grid)`
""" """
n_samples = tidal_tensor.T00.size # TODO needs to be checked further
# We create a array and then calculate the eigenvalues. grid = tidal_tensor.T00.shape[0]
Teval = numpy.full((n_samples, 3, 3), numpy.nan, dtype=numpy.float32) eigvals = numpy.full((grid, grid, grid, 3), numpy.nan,
Teval[:, 0, 0] = tidal_tensor.T00 dtype=numpy.float32)
Teval[:, 0, 1] = tidal_tensor.T01 dummy = numpy.full((3, 3), numpy.nan, dtype=numpy.float32)
Teval[:, 0, 2] = tidal_tensor.T02
Teval[:, 1, 1] = tidal_tensor.T11
Teval[:, 1, 2] = tidal_tensor.T12
Teval[:, 2, 2] = tidal_tensor.T22
eigvals = numpy.full((n_samples, 3), numpy.nan, dtype=numpy.float32) # FILL IN THESER ARGUMENTS
for i in range(n_samples): tidal_tensor_to_eigenvalues(eigvals, dummy, ...)
eigvals[i, :] = numpy.linalg.eigvalsh(Teval[i, ...], 'U')
eigvals[i, :] = numpy.sort(eigvals[i, :])
return eigvals return eigvals
@ -430,3 +419,23 @@ class TidalTensorField(BaseField):
""" """
return MASL.tidal_tensor(overdensity_field, self.box._omega_m, return MASL.tidal_tensor(overdensity_field, self.box._omega_m,
self.box._aexp, self.MAS) self.box._aexp, self.MAS)
@jit(nopython=True)
def tidal_tensor_to_eigenvalues(eigvals, dummy, T00, T01, T02, T11, T12, T22):
"""
TODO: needs to be checked further.
"""
grid = T00.shape[0]
for i in range(grid):
for j in range(grid):
for k in range(grid):
dummy[0, 0] = T00[i, j, k]
dummy[0, 1] = T01[i, j, k]
dummy[0, 2] = T02[i, j, k]
dummy[1, 1] = T11[i, j, k]
dummy[1, 2] = T12[i, j, k]
dummy[2, 2] = T22[i, j, k]
eigvals[i, j, k, :] = numpy.linalg.eigvalsh(dummy, 'U')
eigvals[i, j, k, :] = numpy.sort(eigvals[i, j, k, :])
return eigvals

View file

@ -260,7 +260,6 @@ def fill_outside(field, fill_value, rmax, boxsize):
N = imax N = imax
# Squared radial distance from the center of the box in box units. # Squared radial distance from the center of the box in box units.
rmax_box2 = (N * rmax / boxsize)**2 rmax_box2 = (N * rmax / boxsize)**2
# print("Box ", rmax_box2)
for i in range(N): for i in range(N):
idist2 = (i - 0.5 * (N - 1))**2 idist2 = (i - 0.5 * (N - 1))**2
@ -268,7 +267,6 @@ def fill_outside(field, fill_value, rmax, boxsize):
jdist2 = (j - 0.5 * (N - 1))**2 jdist2 = (j - 0.5 * (N - 1))**2
for k in range(N): for k in range(N):
kdist2 = (k - 0.5 * (N - 1))**2 kdist2 = (k - 0.5 * (N - 1))**2
# print(idist2 + jdist2 + kdist2 > rmax_box2)
if idist2 + jdist2 + kdist2 > rmax_box2: if idist2 + jdist2 + kdist2 > rmax_box2:
field[i, j, k] = fill_value field[i, j, k] = fill_value
return field return field

View file

@ -58,7 +58,6 @@ class BaseStructure(ABC):
@info.setter @info.setter
def info(self, info): def info(self, info):
# TODO turn this into a structured array and add some checks
self._info = info self._info = info
@property @property

View file

@ -19,10 +19,11 @@ the final snapshot.
from math import floor from math import floor
import numpy import numpy
from numba import jit from scipy.integrate import cumulative_trapezoid, quad
from scipy.integrate import quad
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
from scipy.stats import gaussian_kde, kstest from scipy.stats import gaussian_kde, kstest
from numba import jit
from tqdm import tqdm from tqdm import tqdm
@ -205,54 +206,57 @@ class NearestNeighbourReader:
Archive with keys `ndist`, `rdist`, `mass`, `cross_hindxs`` Archive with keys `ndist`, `rdist`, `mass`, `cross_hindxs``
""" """
assert simname in ["csiborg", "quijote"] assert simname in ["csiborg", "quijote"]
fpath = self.paths.cross_nearest(simname, run, nsim, nobs) fpath = self.paths.cross_nearest(simname, run, "dist", nsim, nobs)
return numpy.load(fpath) return numpy.load(fpath)
def build_dist(self, simname, run, kind, verbose=True): def count_neighbour(self, out, ndist, rdist):
""" """
Build the a PDF or a CDF for the nearest neighbour distribution. Count the number of neighbours for each halo as a function of its
Counts the binned number of neighbour for each halo as a funtion of its radial distance.
radial distance from the centre of the high-resolution region.
Parameters Parameters
---------- ----------
simname : str out : 2-dimensional array of shape `(nbins_radial, nbins_neighbour)`
Simulation name. Must be either `csiborg` or `quijote`. Output array to write to. Results are added to this array.
run : str ndist : 2-dimensional array of shape `(nhalos, ncross_simulations)`
Run name. Distance of each halo to its nearest neighbour from a cross
kind : str simulation.
Distribution kind. Either `pdf` or `cdf`. rdist : 1-dimensional array of shape `(nhalos, )`
verbose : bool, optional Distance of each halo to the centre of the high-resolution region.
Verbosity flag.
Returns
-------
out : 2-dimensional array of shape `(nbins_radial, nbins_neighbour)`
"""
return count_neighbour(out, ndist, rdist, self.radial_bin_edges,
self.rmax_neighbour, self.nbins_neighbour)
def build_dist(self, counts, kind):
"""
Build the a PDF or a CDF for the nearest neighbour distribution from
binned counts as a function of radial distance from the centre of the
high-resolution region.
Parameters
----------
counts : 2-dimensional array of shape `(nbins_radial, nbins_neighbour)`
Binned counts of the number of neighbours as a function of
radial distance.
Returns Returns
------- -------
dist : 2-dimensional array of shape `(nbins_radial, nbins_neighbour)` dist : 2-dimensional array of shape `(nbins_radial, nbins_neighbour)`
""" """
assert simname in ["csiborg", "quijote"]
assert kind in ["pdf", "cdf"] assert kind in ["pdf", "cdf"]
rbin_edges = self.radial_bin_edges
# We first bin the distances as a function of each reference halo
# radial distance and then its nearest neighbour distance.
fpaths = self.paths.cross_nearest(simname, run)
if simname == "quijote":
fpaths = fpaths
out = numpy.zeros((self.nbins_radial, self.nbins_neighbour),
dtype=numpy.float32)
for fpath in tqdm(fpaths) if verbose else fpaths:
data = numpy.load(fpath)
out = count_neighbour(
out, data["ndist"], data["rdist"], rbin_edges,
self.rmax_neighbour, self.nbins_neighbour)
if kind == "pdf": if kind == "pdf":
neighbour_bin_edges = self.neighbour_bin_edges neighbour_bin_edges = self.neighbour_bin_edges
dx = neighbour_bin_edges[1] - neighbour_bin_edges[0] dx = neighbour_bin_edges[1] - neighbour_bin_edges[0]
out /= numpy.sum(dx * out, axis=1).reshape(-1, 1) counts /= numpy.sum(dx * counts, axis=1).reshape(-1, 1)
else: else:
out = numpy.cumsum(out, axis=1, out=out) x = self.bin_centres("neighbour")
out /= out[:, -1].reshape(-1, 1) counts = cumulative_trapezoid(counts, x, axis=1, initial=0)
return out counts /= counts[:, -1].reshape(-1, 1)
return counts
def kl_divergence(self, simname, run, nsim, pdf, nobs=None, verbose=True): def kl_divergence(self, simname, run, nsim, pdf, nobs=None, verbose=True):
r""" r"""

View file

@ -410,7 +410,7 @@ class Paths:
fname = f"halo_counts_{simname}_{str(nsim).zfill(5)}.npz" fname = f"halo_counts_{simname}_{str(nsim).zfill(5)}.npz"
return join(fdir, fname) return join(fdir, fname)
def cross_nearest(self, simname, run, nsim=None, nobs=None): def cross_nearest(self, simname, run, kind, nsim=None, nobs=None):
""" """
Path to the files containing distance from a halo in a reference Path to the files containing distance from a halo in a reference
simulation to the nearest halo from a cross simulation. simulation to the nearest halo from a cross simulation.
@ -421,6 +421,9 @@ class Paths:
Simulation name. Must be one of: `csiborg`, `quijote`. Simulation name. Must be one of: `csiborg`, `quijote`.
run : str run : str
Run name. Run name.
kind : str
Whether raw distances or counts in bins. Must be one of `dist`,
`bin_dist` or `tot_counts`.
nsim : int, optional nsim : int, optional
IC realisation index. IC realisation index.
nobs : int, optional nobs : int, optional
@ -431,6 +434,7 @@ class Paths:
path : str path : str
""" """
assert simname in ["csiborg", "quijote"] assert simname in ["csiborg", "quijote"]
assert kind in ["dist", "bin_dist", "tot_counts"]
fdir = join(self.postdir, "nearest_neighbour") fdir = join(self.postdir, "nearest_neighbour")
if not isdir(fdir): if not isdir(fdir):
makedirs(fdir) makedirs(fdir)
@ -440,9 +444,9 @@ class Paths:
nsim = str(nsim).zfill(5) nsim = str(nsim).zfill(5)
else: else:
nsim = self.quijote_fiducial_nsim(nsim, nobs) nsim = self.quijote_fiducial_nsim(nsim, nobs)
return join(fdir, f"{simname}_nn_{nsim}_{run}.npz") return join(fdir, f"{simname}_nn_{kind}_{nsim}_{run}.npz")
files = glob(join(fdir, f"{simname}_nn_*")) files = glob(join(fdir, f"{simname}_nn_{kind}_*"))
run = "_" + run run = "_" + run
return [f for f in files if run in f] return [f for f in files if run in f]

View file

@ -36,7 +36,7 @@ class PKReader:
Output precision. By default `numpy.float32`. Output precision. By default `numpy.float32`.
""" """
def __init__(self, ics, hw, fskel=None, dtype=numpy.float32): def __init__(self, ics, hw, fskel=None, dtype=numpy.float32):
self.ics= ics self.ics = ics
self.hw = hw self.hw = hw
if fskel is None: if fskel is None:
fskel = "/mnt/extraspace/rstiskalek/csiborg/crosspk/out_{}_{}_{}.p" fskel = "/mnt/extraspace/rstiskalek/csiborg/crosspk/out_{}_{}_{}.p"

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

File diff suppressed because one or more lines are too long

View file

@ -19,12 +19,14 @@ MPI parallelized over the reference simulations.
from argparse import ArgumentParser from argparse import ArgumentParser
from datetime import datetime from datetime import datetime
from distutils.util import strtobool from distutils.util import strtobool
from os import remove
import numpy import numpy
import yaml import yaml
from mpi4py import MPI from mpi4py import MPI
from taskmaster import work_delegation from taskmaster import work_delegation
from tqdm import trange
from utils import open_catalogues from utils import open_catalogues
try: try:
@ -36,7 +38,7 @@ except ModuleNotFoundError:
import csiborgtools import csiborgtools
def find_neighbour(args, nsim, cats, paths, comm): def find_neighbour(args, nsim, cats, paths, comm, save_kind):
""" """
Find the nearest neighbour of each halo in the given catalogue. Find the nearest neighbour of each halo in the given catalogue.
@ -53,23 +55,78 @@ def find_neighbour(args, nsim, cats, paths, comm):
Paths object. Paths object.
comm : mpi4py.MPI.Comm comm : mpi4py.MPI.Comm
MPI communicator. MPI communicator.
save_kind : str
Kind of data to save. Must be either `dist` or `bin_dist`.
Returns Returns
------- -------
None None
""" """
assert save_kind in ["dist", "bin_dist"]
ndist, cross_hindxs = csiborgtools.match.find_neighbour(nsim, cats) ndist, cross_hindxs = csiborgtools.match.find_neighbour(nsim, cats)
mass_key = "totpartmass" if args.simname == "csiborg" else "group_mass" mass_key = "totpartmass" if args.simname == "csiborg" else "group_mass"
cat0 = cats[nsim] cat0 = cats[nsim]
mass = cat0[mass_key]
rdist = cat0.radial_distance(in_initial=False) rdist = cat0.radial_distance(in_initial=False)
fout = paths.cross_nearest(args.simname, args.run, nsim) # Distance is saved optionally, whereas binned distance is always saved.
if save_kind == "dist":
out = {"ndist": ndist,
"cross_hindxs": cross_hindxs,
"mass": cat0[mass_key],
"ref_hindxs": cat0["index"],
"rdist": rdist}
fout = paths.cross_nearest(args.simname, args.run, "dist", nsim)
if args.verbose:
print(f"Rank {comm.Get_rank()} writing to `{fout}`.", flush=True)
numpy.savez(fout, **out)
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
reader = csiborgtools.read.NearestNeighbourReader(
paths=paths, **csiborgtools.neighbour_kwargs)
counts = numpy.zeros((reader.nbins_radial, reader.nbins_neighbour),
dtype=numpy.float32)
counts = reader.count_neighbour(counts, ndist, rdist)
out = {"counts": counts}
fout = paths.cross_nearest(args.simname, args.run, "bin_dist", nsim)
if args.verbose: if args.verbose:
print(f"Rank {comm.Get_rank()} writing to `{fout}`.", flush=True) print(f"Rank {comm.Get_rank()} writing to `{fout}`.", flush=True)
numpy.savez(fout, ndist=ndist, cross_hindxs=cross_hindxs, mass=mass, numpy.savez(fout, **out)
ref_hindxs=cat0["index"], rdist=rdist)
def collect_dist(args, paths):
"""
Collect the binned nearest neighbour distances into a single file.
Parameters
----------
args : argparse.Namespace
Command line arguments.
paths : csiborgtools.paths.Paths
Paths object.
Returns
-------
"""
fnames = paths.cross_nearest(args.simname, args.run, "bin_dist")
if args.verbose:
print("Collecting counts into a single file.", flush=True)
for i in trange(len(fnames)) if args.verbose else range(len(fnames)):
fname = fnames[i]
data = numpy.load(fname)
if i == 0:
out = data["counts"]
else:
out += data["counts"]
remove(fname)
fout = paths.cross_nearest(args.simname, args.run, "tot_counts",
nsim=0, nobs=0)
if args.verbose:
print(f"Writing the summed counts to `{fout}`.", flush=True)
numpy.savez(fout, tot_counts=out)
if __name__ == "__main__": if __name__ == "__main__":
@ -87,16 +144,23 @@ if __name__ == "__main__":
with open("./match_finsnap.yml", "r") as file: with open("./match_finsnap.yml", "r") as file:
config = yaml.safe_load(file) config = yaml.safe_load(file)
if args.simname == "csiborg":
save_kind = "dist"
else:
save_kind = "bin_dist"
comm = MPI.COMM_WORLD comm = MPI.COMM_WORLD
rank = comm.Get_rank()
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
cats = open_catalogues(args, config, paths, comm) cats = open_catalogues(args, config, paths, comm)
def do_work(nsim): def do_work(nsim):
return find_neighbour(args, nsim, cats, paths, comm) return find_neighbour(args, nsim, cats, paths, comm, save_kind)
work_delegation(do_work, list(cats.keys()), comm, work_delegation(do_work, list(cats.keys()), comm,
master_verbose=args.verbose) master_verbose=args.verbose)
comm.Barrier() comm.Barrier()
if comm.Get_rank() == 0: if rank == 0:
collect_dist(args, paths)
print(f"{datetime.now()}: all finished. Quitting.") print(f"{datetime.now()}: all finished. Quitting.")

View file

@ -18,20 +18,77 @@ nbins_marks: 10
name: name:
- totpartmass - totpartmass
- group_mass - group_mass
min: 1.e+12 min: 12.4
max: 1.e+13 max: 12.8
islog: true
"mass002": "mass002":
primary: primary:
name: name:
- totpartmass - totpartmass
- group_mass - group_mass
min: 1.e+13 min: 12.6
max: 1.e+14 max: 13.0
islog: true
"mass003": "mass003":
primary: primary:
name: name:
- totpartmass - totpartmass
- group_mass - group_mass
min: 1.e+14 min: 12.8
max: 13.2
islog: true
"mass004":
primary:
name:
- totpartmass
- group_mass
min: 13.0
max: 13.4
islog: true
"mass005":
primary:
name:
- totpartmass
- group_mass
min: 13.2
max: 13.6
islog: true
"mass006":
primary:
name:
- totpartmass
- group_mass
min: 13.4
max: 13.8
islog: true
"mass007":
primary:
name:
- totpartmass
- group_mass
min: 13.6
max: 14.0
islog: true
"mass008":
primary:
name:
- totpartmass
- group_mass
min: 13.8
max: 14.2
islog: true
"mass009":
primary:
name:
- totpartmass
- group_mass
min: 14.0
islog: true

View file

@ -106,8 +106,12 @@ def read_single_catalogue(args, config, nsim, run, rmax, paths, nobs=None):
pname = _name pname = _name
if pname is None: if pname is None:
raise KeyError(f"Invalid names `{sel['name']}`.") raise KeyError(f"Invalid names `{sel['name']}`.")
xmin = sel.get("min", None)
cat.apply_bounds({pname: (sel.get("min", None), sel.get("max", None))}) xmax = sel.get("max", None)
if sel.get("islog", False):
xmin = 10**xmin if xmin is not None else None
xmax = 10**xmax if xmax is not None else None
cat.apply_bounds({pname: (xmin, xmax)})
# Now the secondary selection bounds. If needed transfrom the secondary # Now the secondary selection bounds. If needed transfrom the secondary
# property before applying the bounds. # property before applying the bounds.

View file

@ -21,7 +21,7 @@ import numpy
import healpy import healpy
import scienceplots # noqa import scienceplots # noqa
import utils import plt_utils
from cache_to_disk import cache_to_disk, delete_disk_caches_for_function # noqa from cache_to_disk import cache_to_disk, delete_disk_caches_for_function # noqa
from tqdm import tqdm from tqdm import tqdm
@ -35,7 +35,16 @@ except ModuleNotFoundError:
def open_csiborg(nsim): def open_csiborg(nsim):
""" """
Open a CSiBORG halo catalogue. Open a CSiBORG halo catalogue. Applies mass and distance selection.
Parameters
----------
nsim : int
Simulation index.
Returns
-------
cat : csiborgtools.read.HaloCatalogue
""" """
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
bounds = {"totpartmass": (None, None), "dist": (0, 155/0.705)} bounds = {"totpartmass": (None, None), "dist": (0, 155/0.705)}
@ -43,6 +52,20 @@ def open_csiborg(nsim):
def open_quijote(nsim, nobs=None): def open_quijote(nsim, nobs=None):
"""
Open a Quijote halo catalogue. Applies mass and distance selection.
Parameters
----------
nsim : int
Simulation index.
nobs : int, optional
Fiducial observer index.
Returns
-------
cat : csiborgtools.read.QuijoteHaloCatalogue
"""
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
cat = csiborgtools.read.QuijoteHaloCatalogue(nsim, paths, nsnap=4) cat = csiborgtools.read.QuijoteHaloCatalogue(nsim, paths, nsnap=4)
if nobs is not None: if nobs is not None:
@ -51,10 +74,24 @@ def open_quijote(nsim, nobs=None):
def plot_mass_vs_ncells(nsim, pdf=False): def plot_mass_vs_ncells(nsim, pdf=False):
"""
Plot the halo mass vs. number of occupied cells in the initial snapshot.
Parameters
----------
nsim : int
Simulation index.
pdf : bool, optional
Whether to save the figure as a PDF file.
Returns
-------
None
"""
cat = open_csiborg(nsim) cat = open_csiborg(nsim)
mpart = 4.38304044e+09 mpart = 4.38304044e+09
with plt.style.context(utils.mplstyle): with plt.style.context(plt_utils.mplstyle):
plt.figure() plt.figure()
plt.scatter(cat["totpartmass"], cat["lagpatch_ncells"], s=0.25, plt.scatter(cat["totpartmass"], cat["lagpatch_ncells"], s=0.25,
rasterized=True) rasterized=True)
@ -66,9 +103,9 @@ def plot_mass_vs_ncells(nsim, pdf=False):
plt.ylabel(r"$N_{\rm cells}$") plt.ylabel(r"$N_{\rm cells}$")
for ext in ["png"] if pdf is False else ["png", "pdf"]: for ext in ["png"] if pdf is False else ["png", "pdf"]:
fout = join(utils.fout, f"init_mass_vs_ncells_{nsim}.{ext}") fout = join(plt_utils.fout, f"init_mass_vs_ncells_{nsim}.{ext}")
print(f"Saving to `{fout}`.") print(f"Saving to `{fout}`.")
plt.savefig(fout, dpi=utils.dpi, bbox_inches="tight") plt.savefig(fout, dpi=plt_utils.dpi, bbox_inches="tight")
plt.close() plt.close()
@ -77,13 +114,19 @@ def plot_mass_vs_ncells(nsim, pdf=False):
############################################################################### ###############################################################################
def process_counts(counts):
mean = numpy.mean(counts, axis=0)
std = numpy.std(counts, axis=0)
return mean, std
def plot_hmf(pdf=False): def plot_hmf(pdf=False):
"""
Plot the (ultimate paretn) halo mass function of CSiBORG and Quijote.
Parameters
----------
pdf : bool, optional
Whether to save the figure as a PDF file.
Returns
-------
None
"""
print("Plotting the HMF...", flush=True) print("Plotting the HMF...", flush=True)
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
@ -114,35 +157,39 @@ def plot_hmf(pdf=False):
x = 10**(0.5 * (bins[1:] + bins[:-1])) x = 10**(0.5 * (bins[1:] + bins[:-1]))
# Edit lower limits # Edit lower limits
csiborg_counts[:, x < 1e12] = numpy.nan csiborg_counts[:, x < 1e12] = numpy.nan
quijote_counts[:, x < 8e12] = numpy.nan quijote_counts[:, x < 10**(12.4)] = numpy.nan
# Edit upper limits # Edit upper limits
csiborg_counts[:, x > 4e15] = numpy.nan csiborg_counts[:, x > 4e15] = numpy.nan
quijote_counts[:, x > 4e15] = numpy.nan quijote_counts[:, x > 4e15] = numpy.nan
with plt.style.context(utils.mplstyle): with plt.style.context(plt_utils.mplstyle):
cols = plt.rcParams["axes.prop_cycle"].by_key()["color"] cols = plt.rcParams["axes.prop_cycle"].by_key()["color"]
fig, ax = plt.subplots(nrows=2, sharex=True, fig, ax = plt.subplots(nrows=2, sharex=True,
figsize=(3.5, 2.625 * 1.25), figsize=(3.5, 2.625 * 1.25),
gridspec_kw={"height_ratios": [1, 0.65]}) gridspec_kw={"height_ratios": [1, 0.65]})
fig.subplots_adjust(hspace=0, wspace=0) fig.subplots_adjust(hspace=0, wspace=0)
mean_csiborg, std_csiborg = process_counts(csiborg_counts) # Upper panel data
mean_csiborg = numpy.mean(csiborg_counts, axis=0)
std_csiborg = numpy.std(csiborg_counts, axis=0)
ax[0].plot(x, mean_csiborg, label="CSiBORG") ax[0].plot(x, mean_csiborg, label="CSiBORG")
ax[0].fill_between(x, mean_csiborg - std_csiborg, ax[0].fill_between(x, mean_csiborg - std_csiborg,
mean_csiborg + std_csiborg, alpha=0.5) mean_csiborg + std_csiborg, alpha=0.5)
mean_quijote, std_quijote = process_counts(quijote_counts) mean_quijote = numpy.mean(quijote_counts, axis=0)
std_quijote = numpy.std(quijote_counts, axis=0)
ax[0].plot(x, mean_quijote, label="Quijote") ax[0].plot(x, mean_quijote, label="Quijote")
ax[0].fill_between(x, mean_quijote - std_quijote, ax[0].fill_between(x, mean_quijote - std_quijote,
mean_quijote + std_quijote, alpha=0.5) mean_quijote + std_quijote, alpha=0.5)
# Lower panel data
log_y = numpy.log10(mean_csiborg / mean_quijote) log_y = numpy.log10(mean_csiborg / mean_quijote)
err = numpy.sqrt((std_csiborg / mean_csiborg / numpy.log(10))**2 err = numpy.sqrt((std_csiborg / mean_csiborg / numpy.log(10))**2
+ (std_quijote / mean_quijote / numpy.log(10))**2) + (std_quijote / mean_quijote / numpy.log(10))**2)
ax[1].plot(x, 10**log_y, c=cols[2]) ax[1].plot(x, 10**log_y, c=cols[2])
ax[1].fill_between(x, 10**(log_y - err), 10**(log_y + err), alpha=0.5, ax[1].fill_between(x, 10**(log_y - err), 10**(log_y + err), alpha=0.5,
color=cols[2]) color=cols[2])
# Labels and accesories
ax[1].axhline(1, color="k", ls=plt.rcParams["lines.linestyle"], ax[1].axhline(1, color="k", ls=plt.rcParams["lines.linestyle"],
lw=0.5 * plt.rcParams["lines.linewidth"], zorder=0) lw=0.5 * plt.rcParams["lines.linewidth"], zorder=0)
ax[0].set_ylabel(r"$\frac{\mathrm{d} n}{\mathrm{d}\log M_{\rm h}}~\mathrm{dex}^{-1}$") # noqa ax[0].set_ylabel(r"$\frac{\mathrm{d} n}{\mathrm{d}\log M_{\rm h}}~\mathrm{dex}^{-1}$") # noqa
@ -156,14 +203,33 @@ def plot_hmf(pdf=False):
fig.tight_layout(h_pad=0, w_pad=0) fig.tight_layout(h_pad=0, w_pad=0)
for ext in ["png"] if pdf is False else ["png", "pdf"]: for ext in ["png"] if pdf is False else ["png", "pdf"]:
fout = join(utils.fout, f"hmf_comparison.{ext}") fout = join(plt_utils.fout, f"hmf_comparison.{ext}")
print(f"Saving to `{fout}`.") print(f"Saving to `{fout}`.")
fig.savefig(fout, dpi=utils.dpi, bbox_inches="tight") fig.savefig(fout, dpi=plt_utils.dpi, bbox_inches="tight")
plt.close() plt.close()
@cache_to_disk(7)
def load_field(kind, nsim, grid, MAS, in_rsp=False): def load_field(kind, nsim, grid, MAS, in_rsp=False):
"""
Load a single field.
Parameters
----------
kind : str
Field kind.
nsim : int
Simulation index.
grid : int
Grid size.
MAS : str
Mass assignment scheme.
in_rsp : bool, optional
Whether to load the field in redshift space.
Returns
-------
field : n-dimensional array
"""
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) 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))
@ -175,6 +241,30 @@ 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, MAS="PCS",
highres_only=True, pdf=False): highres_only=True, pdf=False):
"""
Plot the mean projected field.
Parameters
----------
kind : str
Field kind.
nsim : int
Simulation index.
grid : int
Grid size.
in_rsp : bool
Whether to load the field in redshift space.
MAS : str, optional
Mass assignment scheme.
highres_only : bool, optional
Whether to only plot the high-resolution region.
pdf : bool, optional
Whether to save the figure as a PDF.
Returns
-------
None
"""
print(f"Plotting projected field `{kind}`. ", flush=True) print(f"Plotting projected field `{kind}`. ", flush=True)
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
nsnap = max(paths.get_snapshots(nsim)) nsnap = max(paths.get_snapshots(nsim))
@ -190,14 +280,12 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS",
if highres_only: if highres_only:
csiborgtools.field.fill_outside(field, numpy.nan, rmax=155.5, csiborgtools.field.fill_outside(field, numpy.nan, rmax=155.5,
boxsize=677.7) boxsize=677.7)
# start = field.shape[0] // 4
start = round(field.shape[0] * 0.27) start = round(field.shape[0] * 0.27)
end = round(field.shape[0] * 0.73) end = round(field.shape[0] * 0.73)
# end = field.shape[0] - start
field = field[start:end, start:end, start:end] field = field[start:end, start:end, start:end]
labels = [r"$y-z$", r"$x-z$", r"$x-y$"] labels = [r"$y-z$", r"$x-z$", r"$x-y$"]
with plt.style.context(utils.mplstyle): with plt.style.context(plt_utils.mplstyle):
fig, ax = plt.subplots(figsize=(3.5 * 2, 2.625), ncols=3, sharey=True, fig, ax = plt.subplots(figsize=(3.5 * 2, 2.625), ncols=3, sharey=True,
sharex=True) sharex=True)
fig.subplots_adjust(hspace=0, wspace=0) fig.subplots_adjust(hspace=0, wspace=0)
@ -216,9 +304,10 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS",
fig.tight_layout(h_pad=0, w_pad=0) fig.tight_layout(h_pad=0, w_pad=0)
for ext in ["png"] if pdf is False else ["png", "pdf"]: for ext in ["png"] if pdf is False else ["png", "pdf"]:
fout = join(utils.fout, f"field_{kind}_{nsim}_rsp{in_rsp}.{ext}") fout = join(plt_utils.fout,
f"field_{kind}_{nsim}_rsp{in_rsp}.{ext}")
print(f"Saving to `{fout}`.") print(f"Saving to `{fout}`.")
fig.savefig(fout, dpi=utils.dpi, bbox_inches="tight") fig.savefig(fout, dpi=plt_utils.dpi, bbox_inches="tight")
plt.close() plt.close()
############################################################################### ###############################################################################
@ -227,6 +316,20 @@ def plot_projected_field(kind, nsim, grid, in_rsp, MAS="PCS",
def get_sky_label(kind, volume_weight): def get_sky_label(kind, volume_weight):
"""
Get the sky label for a given field kind.
Parameters
----------
kind : str
Field kind.
volume_weight : bool
Whether to volume weight the field.
Returns
-------
label : str
"""
if volume_weight: if volume_weight:
if kind == "density": if kind == "density":
label = r"$\log \int_{0}^{R} r^2 \rho(r, \mathrm{RA}, \mathrm{dec}) \mathrm{d} r$" # noqa label = r"$\log \int_{0}^{R} r^2 \rho(r, \mathrm{RA}, \mathrm{dec}) \mathrm{d} r$" # noqa
@ -255,8 +358,38 @@ def get_sky_label(kind, volume_weight):
def plot_sky_distribution(kind, nsim, grid, nside, MAS="PCS", plot_groups=True, def plot_sky_distribution(kind, nsim, grid, nside, MAS="PCS", plot_groups=True,
dmin=0, dmax=220, plot_halos=None, dmin=0, dmax=220, plot_halos=None,
volume_weight=True, pdf=False): volume_weight=True, pdf=False):
""" r"""
NOTE: add distance for groups. Plot the sky distribution of a given field kind on the sky along with halos
and selected observations.
TODO
----
- Add distance for groups.
Parameters
----------
field : str
Field kind.
nsim : int
Simulation index.
grid : int
Grid size.
nside : int
Healpix nside of the sky projection.
MAS : str, optional
Mass assignment scheme.
plot_groups : bool, optional
Whether to plot the 2M++ groups.
dmin : float, optional
Minimum projection distance in :math:`\mathrm{Mpc}/h`.
dmax : float, optional
Maximum projection distance in :math:`\mathrm{Mpc}/h`.
plot_halos : list, optional
Minimum halo mass to plot in :math:`M_\odot`.
volume_weight : bool, optional
Whether to volume weight the field.
pdf : bool, optional
Whether to save the figure as a pdf.
""" """
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
nsnap = max(paths.get_snapshots(nsim)) nsnap = max(paths.get_snapshots(nsim))
@ -274,7 +407,7 @@ def plot_sky_distribution(kind, nsim, grid, nside, MAS="PCS", plot_groups=True,
out = csiborgtools.field.make_sky(field, angpos=angpos, dist=dist, box=box, out = csiborgtools.field.make_sky(field, angpos=angpos, dist=dist, box=box,
volume_weight=volume_weight) volume_weight=volume_weight)
with plt.style.context(utils.mplstyle): with plt.style.context(plt_utils.mplstyle):
label = get_sky_label(kind, volume_weight) label = get_sky_label(kind, volume_weight)
if kind in ["density", "overdensity"]: if kind in ["density", "overdensity"]:
out = numpy.log10(out) out = numpy.log10(out)
@ -299,9 +432,9 @@ def plot_sky_distribution(kind, nsim, grid, nside, MAS="PCS", plot_groups=True,
plt.legend(markerscale=10) plt.legend(markerscale=10)
for ext in ["png"] if pdf is False else ["png", "pdf"]: for ext in ["png"] if pdf is False else ["png", "pdf"]:
fout = join(utils.fout, f"sky_{kind}_{nsim}_from_{dmin}_to_{dmax}_vol{volume_weight}.{ext}") # noqa fout = join(plt_utils.fout, f"sky_{kind}_{nsim}_from_{dmin}_to_{dmax}_vol{volume_weight}.{ext}") # noqa
print(f"Saving to `{fout}`.") print(f"Saving to `{fout}`.")
plt.savefig(fout, dpi=utils.dpi, bbox_inches="tight") plt.savefig(fout, dpi=plt_utils.dpi, bbox_inches="tight")
plt.close() plt.close()
@ -321,12 +454,19 @@ if __name__ == "__main__":
print(f"Cleaning cache for function {func}.") print(f"Cleaning cache for function {func}.")
delete_disk_caches_for_function(func) delete_disk_caches_for_function(func)
# plot_mass_vs_occupancy(7444) if False:
# plot_mass_vs_normcells(7444 + 24 * 4, pdf=False) plot_mass_vs_ncells(7444, pdf=False)
# plot_mass_vs_ncells(7444, pdf=True)
# plot_hmf(pdf=True)
# plot_sky_distribution("radvel", 7444, 256, nside=64,
# plot_groups=False, dmin=50, dmax=100,
# plot_halos=5e13, volume_weight=False)
plot_projected_field("potential", 7444, 256, in_rsp=True) if False:
plot_hmf(pdf=False)
if False:
plot_sky_distribution("radvel", 7444, 256, nside=64,
plot_groups=False, dmin=50, dmax=100,
plot_halos=5e13, volume_weight=False)
if True:
plot_projected_field("overdensity", 7444, 1024, in_rsp=True,
highres_only=False)
plot_projected_field("overdensity", 7444, 1024, in_rsp=False,
highres_only=False)

View file

@ -19,7 +19,7 @@ import matplotlib.pyplot as plt
import numpy import numpy
import scienceplots # noqa import scienceplots # noqa
import utils import plt_utils
try: try:
import csiborgtools import csiborgtools
@ -40,7 +40,7 @@ def plot_knn(runname):
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
reader = csiborgtools.read.kNNCDFReader(paths) reader = csiborgtools.read.kNNCDFReader(paths)
with plt.style.context(utils.mplstyle): with plt.style.context(plt_utils.mplstyle):
plt.figure() plt.figure()
# Quijote kNN # Quijote kNN
@ -92,9 +92,9 @@ def plot_knn(runname):
plt.ylabel(r"$P(k | V = 4 \pi r^3 / 3)$") plt.ylabel(r"$P(k | V = 4 \pi r^3 / 3)$")
for ext in ["png"]: for ext in ["png"]:
fout = join(utils.fout, f"knn_{runname}.{ext}") fout = join(plt_utils.fout, f"knn_{runname}.{ext}")
print("Saving to `{fout}`.".format(fout=fout)) print("Saving to `{fout}`.".format(fout=fout))
plt.savefig(fout, dpi=utils.dpi, bbox_inches="tight") plt.savefig(fout, dpi=plt_utils.dpi, bbox_inches="tight")
plt.close() plt.close()

File diff suppressed because it is too large Load diff