Aquila-Consortium/csiborgtools
0
0
Fork 0
mirror of https://github.com/Richard-Sti/csiborgtools.git synced 2024-12-22 17:28:02 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Quijote snapshots support (#77)

Browse source 
* Renaming

* Edit docs

* Delete old function

* Add a blank space

* Rename particle reader

* Add comments

* Rename

* Rename

* edit get_snapshots

* More renaming

* Remove old correction

* Add import

* Add basics of the Quijote reader

* Add a blank space

* Fix paths

* Rename function

* Fix HID and path

* Add more FoF reading

* Move definition

* Adding arguments

* Renaming

* Add kwargs for backward comp

* FoF Quijote return only hids

* Add sorting of quijote

* Add path to CSiBORG ICs snapshot

* Add support for Quijote

* initmatch paths for quijote

* Add kwargs

* Fix blank lines

* Rename kwarg

* Remove unused import

* Remove hardcoded numbers

* Update for Quijote

* Do not store velocities in QUijote ICs

* Box units mass Quijote

* Fix typo

* Ensure particles are not right at the edge

* Add structfit paths for QUuijote

* Basic CSiBORG units

* Add more quijote halo reading

* Add Quijote fitting

* Docs changes

* Docs changes
This commit is contained in:
Richard Stiskalek 2023-07-27 18:41:00 +02:00 committed by GitHub
parent e08c741fc8
commit fb4b4edf19
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
18 changed files with 800 additions and 300 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
csiborgtools/fits/halo.py
View file

@ -117,7 +117,7 @@ class BaseStructure(ABC):
assert kind in ["crit", "matter"] assert kind in ["crit", "matter"]
rho = delta_mult * self.box.box_rhoc rho = delta_mult * self.box.box_rhoc
if kind == "matter": if kind == "matter":
rho *= self.box.box_Om rho *= self.box.Om
pos = self.pos pos = self.pos
mass = self["M"] mass = self["M"]
@ -216,8 +216,8 @@ class BaseStructure(ABC):
References References
---------- ----------
[1] A Universal Angular Momentum Profile for Galactic Halos; 2001; [1] A Universal Angular Momentum Profile for Galactic Halos; 2001;
Bullock, J. S.; Dekel, A.; Kolatt, T. S.; Kravtsov, A. V.; Bullock, J. S.; Dekel, A.; Kolatt, T. S.; Kravtsov, A. V.;
Klypin, A. A.; Porciani, C.; Primack, J. R. Klypin, A. A.; Porciani, C.; Primack, J. R.
""" """
pos = self.pos pos = self.pos
mask = periodic_distance(pos, ref, boxsize=1) < rad mask = periodic_distance(pos, ref, boxsize=1) < rad

148
csiborgtools/match/match.py
View file

@ -15,6 +15,7 @@
""" """
Support for matching halos between CSiBORG IC realisations. Support for matching halos between CSiBORG IC realisations.
""" """
from abc import ABC
from datetime import datetime from datetime import datetime
from functools import lru_cache from functools import lru_cache
from math import ceil from math import ceil
@ -26,20 +27,72 @@ from tqdm import tqdm, trange
from ..read import load_halo_particles from ..read import load_halo_particles
BCKG_HALFSIZE = 475
BOX_SIZE = 2048 class BaseMatcher(ABC):
"""
Base class for `RealisationsMatcher` and `ParticleOverlap`.
"""
_box_size = None
_bckg_halfsize = None
@property
def box_size(self):
"""
Number of cells in the box.
Returns
-------
box_size : int
"""
if self._box_size is None:
raise RuntimeError("`box_size` is not set.")
return self._box_size
@box_size.setter
def box_size(self, value):
assert isinstance(value, int)
assert value > 0
self._box_size = value
@property
def bckg_halfsize(self):
"""
Number of to each side of the centre of the box to calculate the
density field. This is because in CSiBORG we are only interested in the
high-resolution region.
Returns
-------
bckg_halfsize : int
"""
if self._bckg_halfsize is None:
raise RuntimeError("`bckg_halfsize` is not set.")
return self._bckg_halfsize
@bckg_halfsize.setter
def bckg_halfsize(self, value):
assert isinstance(value, int)
assert value > 0
self._bckg_halfsize = value
############################################################################### ###############################################################################
# Realisations matcher for calculating overlaps # # Realisations matcher for calculating overlaps #
############################################################################### ###############################################################################
class RealisationsMatcher: class RealisationsMatcher(BaseMatcher):
""" """
A tool to match haloes between IC realisations. A tool to match haloes between IC realisations.
Parameters Parameters
---------- ----------
box_size : int
Number of cells in the box.
bckg_halfsize : int
Number of to each side of the centre of the box to calculate the
density field. This is because in CSiBORG we are only interested in the
high-resolution region.
nmult : float or int, optional nmult : float or int, optional
Multiple of the sum of pair initial Lagrangian patch sizes Multiple of the sum of pair initial Lagrangian patch sizes
within which to return neighbours. By default 1. within which to return neighbours. By default 1.
@ -51,20 +104,22 @@ class RealisationsMatcher:
catalogue key. By default `totpartmass`, i.e. the total particle catalogue key. By default `totpartmass`, i.e. the total particle
mass associated with a halo. mass associated with a halo.
""" """
_nmult = None _nmult = None
_dlogmass = None _dlogmass = None
_mass_kind = None _mass_kind = None
_overlapper = None _overlapper = None
def __init__(self, nmult=1.0, dlogmass=2.0, mass_kind="totpartmass"): def __init__(self, box_size, bckg_halfsize, nmult=1.0, dlogmass=2.0,
mass_kind="totpartmass"):
assert nmult > 0 assert nmult > 0
assert dlogmass > 0 assert dlogmass > 0
assert isinstance(mass_kind, str) assert isinstance(mass_kind, str)
self.box_size = box_size
self.halfsize = bckg_halfsize
self._nmult = nmult self._nmult = nmult
self._dlogmass = dlogmass self._dlogmass = dlogmass
self._mass_kind = mass_kind self._mass_kind = mass_kind
self._overlapper = ParticleOverlap() self._overlapper = ParticleOverlap(box_size, bckg_halfsize)
@property @property
def nmult(self): def nmult(self):
@ -181,7 +236,7 @@ class RealisationsMatcher:
@lru_cache(maxsize=cache_size) @lru_cache(maxsize=cache_size)
def load_cached_halox(hid): def load_cached_halox(hid):
return load_processed_halo(hid, particlesx, halo_mapx, hid2mapx, return load_processed_halo(hid, particlesx, halo_mapx, hid2mapx,
nshift=0, ncells=BOX_SIZE) nshift=0, ncells=self.box_size)
if verbose: if verbose:
print(f"{datetime.now()}: calculating overlaps.", flush=True) print(f"{datetime.now()}: calculating overlaps.", flush=True)
@ -195,7 +250,8 @@ class RealisationsMatcher:
# Next, we find this halo's particles, total mass, minimum and # Next, we find this halo's particles, total mass, minimum and
# maximum cells and convert positions to cells. # maximum cells and convert positions to cells.
pos0, mass0, totmass0, mins0, maxs0 = load_processed_halo( pos0, mass0, totmass0, mins0, maxs0 = load_processed_halo(
k0, particles0, halo_map0, hid2map0, nshift=0, ncells=BOX_SIZE) k0, particles0, halo_map0, hid2map0, nshift=0,
ncells=self.box_size)
# We now loop over matches of this halo and calculate their # We now loop over matches of this halo and calculate their
# overlap, storing them in `_cross`. # overlap, storing them in `_cross`.
@ -267,7 +323,7 @@ class RealisationsMatcher:
@lru_cache(maxsize=cache_size) @lru_cache(maxsize=cache_size)
def load_cached_halox(hid): def load_cached_halox(hid):
return load_processed_halo(hid, particlesx, halo_mapx, hid2mapx, return load_processed_halo(hid, particlesx, halo_mapx, hid2mapx,
nshift=nshift, ncells=BOX_SIZE) nshift=nshift, ncells=self.box_size)
if verbose: if verbose:
print(f"{datetime.now()}: calculating smoothed overlaps.", print(f"{datetime.now()}: calculating smoothed overlaps.",
@ -277,7 +333,7 @@ class RealisationsMatcher:
for i, k0 in enumerate(tqdm(indxs) if verbose else indxs): for i, k0 in enumerate(tqdm(indxs) if verbose else indxs):
pos0, mass0, __, mins0, maxs0 = load_processed_halo( pos0, mass0, __, mins0, maxs0 = load_processed_halo(
k0, particles0, halo_map0, hid2map0, nshift=nshift, k0, particles0, halo_map0, hid2map0, nshift=nshift,
ncells=BOX_SIZE) ncells=self.box_size)
# Now loop over the matches and calculate the smoothed overlap. # Now loop over the matches and calculate the smoothed overlap.
_cross = numpy.full(match_indxs[i].size, numpy.nan, numpy.float32) _cross = numpy.full(match_indxs[i].size, numpy.nan, numpy.float32)
@ -327,13 +383,26 @@ def cosine_similarity(x, y):
return out return out
class ParticleOverlap: class ParticleOverlap(BaseMatcher):
r""" r"""
Class to calculate halo overlaps. The density field calculation is based on Halo overlaps calculator. The density field calculation is based on the
the nearest grid position particle assignment scheme, with optional nearest grid position particle assignment scheme, with optional Gaussian
Gaussian smoothing. smoothing.
Parameters
----------
box_size : int
Number of cells in the box.
bckg_halfsize : int
Number of to each side of the centre of the box to calculate the
density field. This is because in CSiBORG we are only interested in the
high-resolution region.
""" """
def __init__(self, box_size, bckg_halfsize):
self.box_size = box_size
self.bckg_halfsize = bckg_halfsize
def make_bckg_delta(self, particles, halo_map, hid2map, halo_cat, def make_bckg_delta(self, particles, halo_map, hid2map, halo_cat,
delta=None, verbose=False): delta=None, verbose=False):
""" """
@ -363,8 +432,8 @@ class ParticleOverlap:
------- -------
delta : 3-dimensional array delta : 3-dimensional array
""" """
cellmin = BOX_SIZE // 2 - BCKG_HALFSIZE cellmin = self.box_size // 2 - self.bckg_halfsize
cellmax = BOX_SIZE // 2 + BCKG_HALFSIZE cellmax = self.box_size // 2 + self.bckg_halfsize
ncells = cellmax - cellmin ncells = cellmax - cellmin
# We then pre-allocate the density field/check it is of the right shape # We then pre-allocate the density field/check it is of the right shape
if delta is None: if delta is None:
@ -379,7 +448,7 @@ class ParticleOverlap:
continue continue
pos, mass = pos[:, :3], pos[:, 3] pos, mass = pos[:, :3], pos[:, 3]
pos = pos2cell(pos, BOX_SIZE) pos = pos2cell(pos, self.box_size)
# We mask out particles outside the cubical high-resolution region # We mask out particles outside the cubical high-resolution region
mask = numpy.all((cellmin <= pos) & (pos < cellmax), axis=1) mask = numpy.all((cellmin <= pos) & (pos < cellmax), axis=1)
pos = pos[mask] pos = pos[mask]
@ -413,7 +482,7 @@ class ParticleOverlap:
delta : 3-dimensional array delta : 3-dimensional array
""" """
nshift = read_nshift(smooth_kwargs) nshift = read_nshift(smooth_kwargs)
cells = pos2cell(pos, BOX_SIZE) cells = pos2cell(pos, self.box_size)
# Check that minima and maxima are integers # Check that minima and maxima are integers
if not (mins is None and maxs is None): if not (mins is None and maxs is None):
assert mins.dtype.char in numpy.typecodes["AllInteger"] assert mins.dtype.char in numpy.typecodes["AllInteger"]
@ -421,12 +490,12 @@ class ParticleOverlap:
if subbox: if subbox:
if mins is None or maxs is None: if mins is None or maxs is None:
mins, maxs = get_halolims(cells, BOX_SIZE, nshift) mins, maxs = get_halolims(cells, self.box_size, nshift)
ncells = maxs - mins + 1 # To get the number of cells ncells = maxs - mins + 1 # To get the number of cells
else: else:
mins = [0, 0, 0] mins = [0, 0, 0]
ncells = (BOX_SIZE, ) * 3 ncells = (self.box_size, ) * 3
# Preallocate and fill the array # Preallocate and fill the array
delta = numpy.zeros(ncells, dtype=numpy.float32) delta = numpy.zeros(ncells, dtype=numpy.float32)
@ -472,8 +541,8 @@ class ParticleOverlap:
Calculated only if no smoothing is applied, otherwise `None`. Calculated only if no smoothing is applied, otherwise `None`.
""" """
nshift = read_nshift(smooth_kwargs) nshift = read_nshift(smooth_kwargs)
pos1 = pos2cell(pos1, BOX_SIZE) pos1 = pos2cell(pos1, self.box_size)
pos2 = pos2cell(pos2, BOX_SIZE) pos2 = pos2cell(pos2, self.box_size)
xc1, yc1, zc1 = [pos1[:, i] for i in range(3)] xc1, yc1, zc1 = [pos1[:, i] for i in range(3)]
xc2, yc2, zc2 = [pos2[:, i] for i in range(3)] xc2, yc2, zc2 = [pos2[:, i] for i in range(3)]
@ -490,7 +559,7 @@ class ParticleOverlap:
ymax = max(numpy.max(yc1), numpy.max(yc2)) + nshift ymax = max(numpy.max(yc1), numpy.max(yc2)) + nshift
zmax = max(numpy.max(zc1), numpy.max(zc2)) + nshift zmax = max(numpy.max(zc1), numpy.max(zc2)) + nshift
# Make sure shifting does not go beyond boundaries # Make sure shifting does not go beyond boundaries
xmax, ymax, zmax = [min(px, BOX_SIZE - 1) xmax, ymax, zmax = [min(px, self.box_size - 1)
for px in (xmax, ymax, zmax)] for px in (xmax, ymax, zmax)]
else: else:
xmin, ymin, zmin = [min(mins1[i], mins2[i]) for i in range(3)] xmin, ymin, zmin = [min(mins1[i], mins2[i]) for i in range(3)]
@ -573,12 +642,14 @@ class ParticleOverlap:
smooth_kwargs=smooth_kwargs) smooth_kwargs=smooth_kwargs)
if smooth_kwargs is not None: if smooth_kwargs is not None:
return calculate_overlap(delta1, delta2, cellmins, delta_bckg) return calculate_overlap(delta1, delta2, cellmins, delta_bckg,
self.box_size, self.bckg_halfsize)
# Calculate masses not given # Calculate masses not given
totmass1 = numpy.sum(mass1) if totmass1 is None else totmass1 totmass1 = numpy.sum(mass1) if totmass1 is None else totmass1
totmass2 = numpy.sum(mass2) if totmass2 is None else totmass2 totmass2 = numpy.sum(mass2) if totmass2 is None else totmass2
return calculate_overlap_indxs( return calculate_overlap_indxs(delta1, delta2, cellmins, delta_bckg,
delta1, delta2, cellmins, delta_bckg, nonzero, totmass1, totmass2) nonzero, totmass1, totmass2,
self.box_size, self.bckg_halfsize)
############################################################################### ###############################################################################
@ -725,7 +796,8 @@ def get_halolims(pos, ncells, nshift=None):
@jit(nopython=True) @jit(nopython=True)
def calculate_overlap(delta1, delta2, cellmins, delta_bckg): def calculate_overlap(delta1, delta2, cellmins, delta_bckg, box_size,
bckg_halfsize):
r""" r"""
Overlap between two halos whose density fields are evaluated on the Overlap between two halos whose density fields are evaluated on the
same grid. This is a JIT implementation, hence it is outside of the main same grid. This is a JIT implementation, hence it is outside of the main
@ -743,6 +815,12 @@ def calculate_overlap(delta1, delta2, cellmins, delta_bckg):
Summed background density field of the reference and cross simulations Summed background density field of the reference and cross simulations
calculated with particles assigned to halos at the final snapshot. calculated with particles assigned to halos at the final snapshot.
Assumed to only be sampled in cells :math:`[512, 1536)^3`. Assumed to only be sampled in cells :math:`[512, 1536)^3`.
box_size : int
Number of cells in the box.
bckg_halfsize : int
Number of to each side of the centre of the box to calculate the
density field. This is because in CSiBORG we are only interested in the
high-resolution region.
Returns Returns
------- -------
@ -751,8 +829,8 @@ def calculate_overlap(delta1, delta2, cellmins, delta_bckg):
totmass = 0.0 # Total mass of halo 1 and halo 2 totmass = 0.0 # Total mass of halo 1 and halo 2
intersect = 0.0 # Weighted intersecting mass intersect = 0.0 # Weighted intersecting mass
i0, j0, k0 = cellmins # Unpack things i0, j0, k0 = cellmins # Unpack things
bckg_size = 2 * BCKG_HALFSIZE bckg_size = 2 * bckg_halfsize
bckg_offset = BOX_SIZE // 2 - BCKG_HALFSIZE bckg_offset = box_size // 2 - bckg_halfsize
imax, jmax, kmax = delta1.shape imax, jmax, kmax = delta1.shape
for i in range(imax): for i in range(imax):
@ -777,7 +855,7 @@ def calculate_overlap(delta1, delta2, cellmins, delta_bckg):
@jit(nopython=True) @jit(nopython=True)
def calculate_overlap_indxs(delta1, delta2, cellmins, delta_bckg, nonzero, def calculate_overlap_indxs(delta1, delta2, cellmins, delta_bckg, nonzero,
mass1, mass2): mass1, mass2, box_size, bckg_halfsize):
r""" r"""
Overlap between two haloes whose density fields are evaluated on the Overlap between two haloes whose density fields are evaluated on the
same grid and `nonzero1` enumerates the non-zero cells of `delta1. This is same grid and `nonzero1` enumerates the non-zero cells of `delta1. This is
@ -801,6 +879,12 @@ def calculate_overlap_indxs(delta1, delta2, cellmins, delta_bckg, nonzero,
mass1, mass2 : floats, optional mass1, mass2 : floats, optional
Total masses of the two haloes, respectively. Optional. If not provided Total masses of the two haloes, respectively. Optional. If not provided
calculcated directly from the density field. calculcated directly from the density field.
box_size : int
Number of cells in the box.
bckg_halfsize : int
Number of to each side of the centre of the box to calculate the
density field. This is because in CSiBORG we are only interested in the
high-resolution region.
Returns Returns
------- -------
@ -808,8 +892,8 @@ def calculate_overlap_indxs(delta1, delta2, cellmins, delta_bckg, nonzero,
""" """
intersect = 0.0 # Weighted intersecting mass intersect = 0.0 # Weighted intersecting mass
i0, j0, k0 = cellmins # Unpack cell minimas i0, j0, k0 = cellmins # Unpack cell minimas
bckg_size = 2 * BCKG_HALFSIZE bckg_size = 2 * bckg_halfsize
bckg_offset = BOX_SIZE // 2 - BCKG_HALFSIZE bckg_offset = box_size // 2 - bckg_halfsize
for n in range(nonzero.shape[0]): for n in range(nonzero.shape[0]):
i, j, k = nonzero[n, :] i, j, k = nonzero[n, :]

4
csiborgtools/read/__init__.py
View file

@ -23,8 +23,8 @@ from .overlap_summary import (NPairsOverlap, PairOverlap, # noqa
binned_resample_mean, get_cross_sims) # noqa binned_resample_mean, get_cross_sims) # noqa
from .paths import Paths # noqa from .paths import Paths # noqa
from .pk_summary import PKReader # noqa from .pk_summary import PKReader # noqa
from .readsim import (MmainReader, ParticleReader, halfwidth_mask, # noqa from .readsim import (MmainReader, CSiBORGReader, QuijoteReader, halfwidth_mask, # noqa
load_halo_particles, read_initcm) # noqa load_halo_particles) # noqa
from .tpcf_summary import TPCFReader # noqa from .tpcf_summary import TPCFReader # noqa
from .utils import (M200_to_R200, cartesian_to_radec, # noqa from .utils import (M200_to_R200, cartesian_to_radec, # noqa
cols_to_structured, radec_to_cartesian, read_h5, cols_to_structured, radec_to_cartesian, read_h5,

129
csiborgtools/read/box_units.py
View file

@ -21,9 +21,9 @@ import numpy
from astropy import constants, units from astropy import constants, units
from astropy.cosmology import LambdaCDM from astropy.cosmology import LambdaCDM
from .readsim import ParticleReader from .readsim import CSiBORGReader, QuijoteReader
# Map of CSiBORG unit conversions
CSIBORG_CONV_NAME = { CSIBORG_CONV_NAME = {
"length": ["x", "y", "z", "peak_x", "peak_y", "peak_z", "Rs", "rmin", "length": ["x", "y", "z", "peak_x", "peak_y", "peak_z", "Rs", "rmin",
"rmax", "r200c", "r500c", "r200m", "r500m", "x0", "y0", "z0", "rmax", "r200c", "r500c", "r200m", "r500m", "x0", "y0", "z0",
@ -31,8 +31,14 @@ CSIBORG_CONV_NAME = {
"velocity": ["vx", "vy", "vz"], "velocity": ["vx", "vy", "vz"],
"mass": ["mass_cl", "totpartmass", "m200c", "m500c", "mass_mmain", "M", "mass": ["mass_cl", "totpartmass", "m200c", "m500c", "mass_mmain", "M",
"m200m", "m500m"], "m200m", "m500m"],
"density": ["rho0"]} "density": ["rho0"]
}
QUIJOTE_CONV_NAME = {
"length": ["x", "y", "z", "x0", "y0", "z0", "Rs", "r200c", "r500c",
"r200m", "r500m", "lagpatch_size"],
"mass": ["group_mass", "totpartmass", "m200c", "m500c", "m200m", "m500m"],
}
############################################################################### ###############################################################################
# Base box # # Base box #
@ -74,7 +80,7 @@ class BaseBox(ABC):
@property @property
def h(self): def h(self):
r""" r"""
The little 'h` parameter at the time of the snapshot. The little 'h' parameter at the time of the snapshot.
Returns Returns
------- -------
@ -157,12 +163,12 @@ class CSiBORGBox(BaseBox):
""" """
Read in the snapshot info file and set the units from it. Read in the snapshot info file and set the units from it.
""" """
partreader = ParticleReader(paths) partreader = CSiBORGReader(paths)
info = partreader.read_info(nsnap, nsim) info = partreader.read_info(nsnap, nsim)
pars = ["boxlen", "time", "aexp", "H0", "omega_m", "omega_l", pars = ["boxlen", "time", "aexp", "H0", "omega_m", "omega_l",
"omega_k", "omega_b", "unit_l", "unit_d", "unit_t"] "omega_k", "omega_b", "unit_l", "unit_d", "unit_t"]
for par in pars: for par in pars:
setattr(self, "_" + par, float(info[par])) setattr(self, "_" + par, info[par])
self._cosmo = LambdaCDM(H0=self._H0, Om0=self._omega_m, self._cosmo = LambdaCDM(H0=self._H0, Om0=self._omega_m,
Ode0=self._omega_l, Tcmb0=2.725 * units.K, Ode0=self._omega_l, Tcmb0=2.725 * units.K,
@ -226,7 +232,7 @@ class CSiBORGBox(BaseBox):
Returns Returns
------- -------
length : foat length : float
Length in :math:`\mathrm{ckpc}` Length in :math:`\mathrm{ckpc}`
""" """
return length * (self._unit_l / units.kpc.to(units.cm) / self._aexp) return length * (self._unit_l / units.kpc.to(units.cm) / self._aexp)
@ -243,7 +249,7 @@ class CSiBORGBox(BaseBox):
Returns Returns
------- -------
length : foat length : float
Length in box units. Length in box units.
""" """
return length / (self._unit_l / units.kpc.to(units.cm) / self._aexp) return length / (self._unit_l / units.kpc.to(units.cm) / self._aexp)
@ -260,7 +266,7 @@ class CSiBORGBox(BaseBox):
Returns Returns
------- -------
length : foat length : float
Length in box units. Length in box units.
""" """
return self.kpc2box(length * 1e3) return self.kpc2box(length * 1e3)
@ -277,7 +283,7 @@ class CSiBORGBox(BaseBox):
Returns Returns
------- -------
length : foat length : float
Length in :math:`\mathrm{ckpc}` Length in :math:`\mathrm{ckpc}`
""" """
return self.box2kpc(length) * 1e-3 return self.box2kpc(length) * 1e-3
@ -419,17 +425,110 @@ class QuijoteBox(BaseBox):
Empty keyword arguments. For backwards compatibility. Empty keyword arguments. For backwards compatibility.
""" """
def __init__(self, nsnap, **kwargs): def __init__(self, nsnap, nsim, paths):
zdict = {4: 0.0, 3: 0.5, 2: 1.0, 1: 2.0, 0: 3.0} zdict = {4: 0.0, 3: 0.5, 2: 1.0, 1: 2.0, 0: 3.0}
assert nsnap in zdict.keys(), f"`nsnap` must be in {zdict.keys()}." assert nsnap in zdict.keys(), f"`nsnap` must be in {zdict.keys()}."
self._aexp = 1 / (1 + zdict[nsnap]) self._aexp = 1 / (1 + zdict[nsnap])
self._cosmo = LambdaCDM(H0=67.11, Om0=0.3175, Ode0=0.6825, info = QuijoteReader(paths).read_info(nsnap, nsim)
Tcmb0=2.725 * units.K, Ob0=0.049) self._cosmo = LambdaCDM(H0=info["Hubble"], Om0=info["Omega_m"],
Ode0=info["Omega_l"], Tcmb0=2.725 * units.K)
self._info = info
@property @property
def boxsize(self): def boxsize(self):
return 1000. / (self._cosmo.H0.value / 100) return self._info["BoxSize"]
def box2mpc(self, length):
r"""
Convert length from box units to :math:`\mathrm{cMpc} / h`.
Parameters
----------
length : float
Length in box units.
Returns
-------
length : float
Length in :math:`\mathrm{cMpc} / h`
"""
return length * self.boxsize
def mpc2box(self, length):
r"""
Convert length from :math:`\mathrm{cMpc} / h` to box units.
Parameters
----------
length : float
Length in :math:`\mathrm{cMpc} / h`.
Returns
-------
length : float
Length in box units.
"""
return length / self.boxsize
def solarmass2box(self, mass):
r"""
Convert mass from :math:`M_\odot / h` to box units.
Parameters
----------
mass : float
Mass in :math:`M_\odot`.
Returns
-------
mass : float
Mass in box units.
"""
return mass / self._info["TotMass"]
def box2solarmass(self, mass):
r"""
Convert mass from box units to :math:`M_\odot / h`.
Parameters
----------
mass : float
Mass in box units.
Returns
-------
mass : float
Mass in :math:`M_\odot / h`.
"""
return mass * self._info["TotMass"]
def convert_from_box(self, data, names): def convert_from_box(self, data, names):
raise NotImplementedError("Conversion not implemented for Quijote.") names = [names] if isinstance(names, str) else names
transforms = {"length": self.box2mpc,
"mass": self.box2solarmass,
# "velocity": self.box2vel,
# "density": self.box2dens,
}
for name in names:
if name not in data.dtype.names:
continue
# Convert
found = False
for unittype, suppnames in QUIJOTE_CONV_NAME.items():
if name in suppnames:
data[name] = transforms[unittype](data[name])
found = True
continue
# If nothing found
if not found:
raise NotImplementedError(
f"Conversion of `{name}` is not defined.")
# # Center at the observer
# if name in ["x0", "y0", "z0"]:
# data[name] -= transforms["length"](0.5)
return data

86
csiborgtools/read/halo_cat.py
View file

@ -22,7 +22,6 @@ from copy import deepcopy
from functools import lru_cache from functools import lru_cache
from itertools import product from itertools import product
from math import floor from math import floor
from os.path import join
import numpy import numpy
from readfof import FoF_catalog from readfof import FoF_catalog
@ -30,14 +29,14 @@ from sklearn.neighbors import NearestNeighbors
from .box_units import CSiBORGBox, QuijoteBox from .box_units import CSiBORGBox, QuijoteBox
from .paths import Paths from .paths import Paths
from .readsim import ParticleReader from .readsim import CSiBORGReader
from .utils import (add_columns, cartesian_to_radec, cols_to_structured, from .utils import (add_columns, cartesian_to_radec, cols_to_structured,
flip_cols, radec_to_cartesian, real2redshift) flip_cols, radec_to_cartesian, real2redshift)
class BaseCatalogue(ABC): class BaseCatalogue(ABC):
""" """
Base (sub)halo catalogue. Base halo catalogue.
""" """
_data = None _data = None
_paths = None _paths = None
@ -125,11 +124,8 @@ class BaseCatalogue(ABC):
def position(self, in_initial=False, cartesian=True): def position(self, in_initial=False, cartesian=True):
r""" r"""
Position components. If Cartesian, then in :math:`\mathrm{cMpc}`. If Position components. If not Cartesian, then RA is in :math:`[0, 360)`
spherical, then radius is in :math:`\mathrm{cMpc}`, RA in degrees and DEC is in :math:`[-90, 90]` degrees.
:math:`[0, 360)` degrees and DEC in :math:`[-90, 90]` degrees. Note
that the position is defined as the minimum of the gravitationl
potential.
Parameters Parameters
---------- ----------
@ -399,33 +395,33 @@ class CSiBORGHaloCatalogue(BaseCatalogue):
names and the items are a len-2 tuple of (min, max) values. In case of names and the items are a len-2 tuple of (min, max) values. In case of
no minimum or maximum, use `None`. For radial distance from the origin no minimum or maximum, use `None`. For radial distance from the origin
use `dist`. use `dist`.
with_lagpatch : bool, optional
Whether to only load halos with a resolved Lagrangian patch.
load_fitted : bool, optional load_fitted : bool, optional
Whether to load fitted quantities. Whether to load fitted quantities.
load_initial : bool, optional load_initial : bool, optional
Whether to load initial positions. Whether to load initial positions.
with_lagpatch : bool, optional
Whether to only load halos with a resolved Lagrangian patch.
rawdata : bool, optional rawdata : bool, optional
Whether to return the raw data. In this case applies no cuts and Whether to return the raw data. In this case applies no cuts and
transformations. transformations.
""" """
def __init__(self, nsim, paths, bounds={"dist": (0, 155.5 / 0.705)}, def __init__(self, nsim, paths, bounds={"dist": (0, 155.5 / 0.705)},
with_lagpatch=True, load_fitted=True, load_initial=True, load_fitted=True, load_initial=True, with_lagpatch=True,
rawdata=False): rawdata=False):
self.nsim = nsim self.nsim = nsim
self.paths = paths self.paths = paths
reader = ParticleReader(paths) reader = CSiBORGReader(paths)
self._data = reader.read_fof_halos(self.nsim) self._data = reader.read_fof_halos(self.nsim)
if load_fitted: if load_fitted:
fits = numpy.load(paths.structfit(self.nsnap, nsim)) fits = numpy.load(paths.structfit(self.nsnap, nsim, "csiborg"))
cols = [col for col in fits.dtype.names if col != "index"] cols = [col for col in fits.dtype.names if col != "index"]
X = [fits[col] for col in cols] X = [fits[col] for col in cols]
self._data = add_columns(self._data, X, cols) self._data = add_columns(self._data, X, cols)
if load_initial: if load_initial:
fits = numpy.load(paths.initmatch(nsim, "fit")) fits = numpy.load(paths.initmatch(nsim, "csiborg", "fit"))
X, cols = [], [] X, cols = [], []
for col in fits.dtype.names: for col in fits.dtype.names:
if col == "index": if col == "index":
@ -465,7 +461,7 @@ class CSiBORGHaloCatalogue(BaseCatalogue):
@property @property
def nsnap(self): def nsnap(self):
return max(self.paths.get_snapshots(self.nsim)) return max(self.paths.get_snapshots(self.nsim, "csiborg"))
@property @property
def box(self): def box(self):
@ -497,28 +493,35 @@ class QuijoteHaloCatalogue(BaseCatalogue):
nsnap : int nsnap : int
Snapshot index. Snapshot index.
origin : len-3 tuple, optional origin : len-3 tuple, optional
Where to place the origin of the box. By default the centre of the box. Where to place the origin of the box. In units of :math:`cMpc / h`.
In units of :math:`cMpc`.
bounds : dict bounds : dict
Parameter bounds to apply to the catalogue. The keys are the parameter Parameter bounds to apply to the catalogue. The keys are the parameter
names and the items are a len-2 tuple of (min, max) values. In case of names and the items are a len-2 tuple of (min, max) values. In case of
no minimum or maximum, use `None`. For radial distance from the origin no minimum or maximum, use `None`. For radial distance from the origin
use `dist`. use `dist`.
load_initial : bool, optional
Whether to load initial positions.
with_lagpatch : bool, optional
Whether to only load halos with a resolved Lagrangian patch.
rawdata : bool, optional
Whether to return the raw data. In this case applies no cuts and
transformations.
**kwargs : dict **kwargs : dict
Keyword arguments for backward compatibility. Keyword arguments for backward compatibility.
""" """
_nsnap = None _nsnap = None
_origin = None _origin = None
def __init__(self, nsim, paths, nsnap, def __init__(self, nsim, paths, nsnap, origin=[0., 0., 0.],
origin=[500 / 0.6711, 500 / 0.6711, 500 / 0.6711], bounds=None, load_initial=True, with_lagpatch=True,
bounds=None, **kwargs): rawdata=False, **kwargs):
self.paths = paths self.paths = paths
self.nsnap = nsnap self.nsnap = nsnap
self.origin = origin self.origin = origin
self._boxwidth = 1000 / 0.6711 self._box = QuijoteBox(nsnap, nsim, paths)
self._boxwidth = self.box.boxsize
fpath = join(self.paths.quijote_dir, "halos", str(nsim)) fpath = self.paths.fof_cat(nsim, "quijote")
fof = FoF_catalog(fpath, self.nsnap, long_ids=False, swap=False, fof = FoF_catalog(fpath, self.nsnap, long_ids=False, swap=False,
SFR=False, read_IDs=False) SFR=False, read_IDs=False)
@ -529,18 +532,44 @@ class QuijoteHaloCatalogue(BaseCatalogue):
("index", numpy.int32)] ("index", numpy.int32)]
data = cols_to_structured(fof.GroupLen.size, cols) data = cols_to_structured(fof.GroupLen.size, cols)
pos = fof.GroupPos / 1e3 / self.box.h pos = self.box.mpc2box(fof.GroupPos / 1e3)
vel = fof.GroupVel * (1 + self.redshift) vel = fof.GroupVel * (1 + self.redshift)
for i, p in enumerate(["x", "y", "z"]): for i, p in enumerate(["x", "y", "z"]):
data[p] = pos[:, i] - self.origin[i] data[p] = pos[:, i] - self.origin[i]
data["v" + p] = vel[:, i] data["v" + p] = vel[:, i]
data["group_mass"] = fof.GroupMass * 1e10 / self.box.h data["group_mass"] = self.box.solarmass2box(fof.GroupMass * 1e10)
data["npart"] = fof.GroupLen data["npart"] = fof.GroupLen
data["index"] = numpy.arange(data.size, dtype=numpy.int32) # We want to start indexing from 1. Index 0 is reserved for
# particles unassigned to any FoF group.
data["index"] = 1 + numpy.arange(data.size, dtype=numpy.int32)
if load_initial:
fits = numpy.load(paths.initmatch(nsim, "quijote", "fit"))
X, cols = [], []
for col in fits.dtype.names:
if col == "index":
continue
if col in ['x', 'y', 'z']:
cols.append(col + "0")
else:
cols.append(col)
X.append(fits[col])
data = add_columns(data, X, cols)
self._data = data self._data = data
if bounds is not None: if not rawdata:
self.apply_bounds(bounds) if with_lagpatch:
mask = numpy.isfinite(self._data["lagpatch_size"])
self._data = self._data[mask]
names = ["x", "y", "z", "group_mass"]
self._data = self.box.convert_from_box(self._data, names)
if load_initial:
names = ["x0", "y0", "z0", "lagpatch_size"]
self._data = self.box.convert_from_box(self._data, names)
if bounds is not None:
self.apply_bounds(bounds)
@property @property
def nsnap(self): def nsnap(self):
@ -579,7 +608,7 @@ class QuijoteHaloCatalogue(BaseCatalogue):
------- -------
box : instance of :py:class:`csiborgtools.units.BaseBox` box : instance of :py:class:`csiborgtools.units.BaseBox`
""" """
return QuijoteBox(self.nsnap) return self._box
@property @property
def origin(self): def origin(self):
@ -629,6 +658,7 @@ class QuijoteHaloCatalogue(BaseCatalogue):
cat.apply_bounds({"dist": (0, rmax)}) cat.apply_bounds({"dist": (0, rmax)})
return cat return cat
############################################################################### ###############################################################################
# Utility functions for halo catalogues # # Utility functions for halo catalogues #
############################################################################### ###############################################################################

193
csiborgtools/read/paths.py
View file

@ -186,7 +186,7 @@ class Paths:
""" """
return join(self.borg_dir, "mcmc", f"mcmc_{nsim}.h5") return join(self.borg_dir, "mcmc", f"mcmc_{nsim}.h5")
def fof_membership(self, nsim, sorted=False): def fof_membership(self, nsim, simname, sorted=False):
""" """
Path to the file containing the FoF particle membership. Path to the file containing the FoF particle membership.
@ -194,10 +194,15 @@ class Paths:
---------- ----------
nsim : int nsim : int
IC realisation index. IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
sorted : bool, optional sorted : bool, optional
Whether to return path to the file that is sorted in the same Whether to return path to the file that is sorted in the same
order as the PHEW output. order as the PHEW output.
""" """
assert simname in ["csiborg", "quijote"]
if simname == "quijote":
raise RuntimeError("Quijote FoF membership is in the FoF cats..")
fdir = join(self.postdir, "FoF_membership", ) fdir = join(self.postdir, "FoF_membership", )
if not isdir(fdir): if not isdir(fdir):
mkdir(fdir) mkdir(fdir)
@ -207,20 +212,25 @@ class Paths:
fout = fout.replace(".npy", "_sorted.npy") fout = fout.replace(".npy", "_sorted.npy")
return fout return fout
def fof_cat(self, nsim): def fof_cat(self, nsim, simname):
""" r"""
Path to the FoF halo catalogue file. Path to the :math:`z = 0` FoF halo catalogue.
Parameters Parameters
---------- ----------
nsim : int nsim : int
IC realisation index. IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
""" """
fdir = join(self.postdir, "FoF_membership", ) assert simname in ["csiborg", "quijote"]
if not isdir(fdir): if simname == "csiborg":
mkdir(fdir) fdir = join(self.postdir, "FoF_membership", )
warn(f"Created directory `{fdir}`.", UserWarning, stacklevel=1) if not isdir(fdir):
return join(fdir, f"halo_catalog_{nsim}_FOF.txt") mkdir(fdir)
warn(f"Created directory `{fdir}`.", UserWarning)
return join(fdir, f"halo_catalog_{nsim}_FOF.txt")
return join(self.quijote_dir, "Halos_fiducial", str(nsim))
def mmain(self, nsnap, nsim): def mmain(self, nsnap, nsim):
""" """
@ -244,7 +254,7 @@ class Paths:
return join(fdir, return join(fdir,
f"mmain_{str(nsim).zfill(5)}_{str(nsnap).zfill(5)}.npz") f"mmain_{str(nsim).zfill(5)}_{str(nsnap).zfill(5)}.npz")
def initmatch(self, nsim, kind): def initmatch(self, nsim, simname, kind):
""" """
Path to the `initmatch` files where the halo match between the Path to the `initmatch` files where the halo match between the
initial and final snapshot of a CSiBORG realisaiton is stored. initial and final snapshot of a CSiBORG realisaiton is stored.
@ -253,19 +263,30 @@ class Paths:
---------- ----------
nsim : int nsim : int
IC realisation index. IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
kind : str kind : str
Type of match. Must be one of `["particles", "fit", "halomap"]`. Type of match. Must be one of `particles` or `fit`.
Returns Returns
------- -------
path : str path : str
""" """
assert kind in ["particles", "fit", "halomap"] assert simname in ["csiborg", "quijote"]
assert kind in ["particles", "fit"]
ftype = "npy" if kind == "fit" else "h5" ftype = "npy" if kind == "fit" else "h5"
fdir = join(self.postdir, "initmatch")
if simname == "csiborg":
fdir = join(self.postdir, "initmatch")
if not isdir(fdir):
mkdir(fdir)
warn(f"Created directory `{fdir}`.", UserWarning)
return join(fdir, f"{kind}_{str(nsim).zfill(5)}.{ftype}")
fdir = join(self.quijote_dir, "initmatch")
if not isdir(fdir): if not isdir(fdir):
mkdir(fdir) mkdir(fdir)
warn(f"Created directory `{fdir}`.", UserWarning, stacklevel=1) warn(f"Created directory `{fdir}`.", UserWarning)
return join(fdir, f"{kind}_{str(nsim).zfill(5)}.{ftype}") return join(fdir, f"{kind}_{str(nsim).zfill(5)}.{ftype}")
def get_ics(self, simname): def get_ics(self, simname):
@ -276,7 +297,7 @@ class Paths:
Parameters Parameters
---------- ----------
simname : str simname : str
Simulation name. Must be one of `["csiborg", "quijote"]`. Simulation name. Must be `csiborg` or `quijote`.
Returns Returns
------- -------
@ -295,75 +316,127 @@ class Paths:
except ValueError: except ValueError:
pass pass
return numpy.sort(ids) return numpy.sort(ids)
else:
# TODO here later read this from the catalogues instead.
return numpy.arange(100, dtype=int)
def ic_path(self, nsim, tonew=False): files = glob("/mnt/extraspace/rstiskalek/Quijote/Snapshots_fiducial/*")
files = [int(f.split("/")[-1]) for f in files]
return numpy.sort(files)
def snapshots(self, nsim, simname, tonew=False):
""" """
Path to a CSiBORG IC realisation folder. Path to an IC snapshots folder.
Parameters Parameters
---------- ----------
nsim : int nsim : int
IC realisation index. IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
tonew : bool, optional tonew : bool, optional
Whether to return the path to the '_new' IC realisation. Whether to return the path to the '_new' IC realisation of
CSiBORG. Ignored for Quijote.
Returns Returns
------- -------
path : str path : str
""" """
fname = "ramses_out_{}" assert simname in ["csiborg", "quijote"]
if tonew: if simname == "csiborg":
fname += "_new" fname = "ramses_out_{}"
return join(self.postdir, "output", fname.format(nsim)) if tonew:
fname += "_new"
return join(self.postdir, "output", fname.format(nsim))
return join(self.srcdir, fname.format(nsim))
return join(self.srcdir, fname.format(nsim)) return join(self.quijote_dir, "Snapshots_fiducial", str(nsim))
def get_snapshots(self, nsim): def get_snapshots(self, nsim, simname):
""" """
List of available snapshots of a CSiBORG IC realisation. List of available snapshots of simulation.
Parameters Parameters
---------- ----------
nsim : int nsim : int
IC realisation index. IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
Returns Returns
------- -------
snapshots : 1-dimensional array snapshots : 1-dimensional array
""" """
simpath = self.ic_path(nsim, tonew=False) simpath = self.snapshots(nsim, simname, tonew=False)
# Get all files in simpath that start with output_ if simname == "csiborg":
snaps = glob(join(simpath, "output_*")) # Get all files in simpath that start with output_
# Take just the last _00XXXX from each file and strip zeros snaps = glob(join(simpath, "output_*"))
snaps = [int(snap.split("_")[-1].lstrip("0")) for snap in snaps] # Take just the last _00XXXX from each file and strip zeros
snaps = [int(snap.split("_")[-1].lstrip("0")) for snap in snaps]
else:
snaps = glob(join(simpath, "snapdir_*"))
snaps = [int(snap.split("/")[-1].split("snapdir_")[-1])
for snap in snaps]
return numpy.sort(snaps) return numpy.sort(snaps)
def snapshot(self, nsnap, nsim): def snapshot(self, nsnap, nsim, simname):
""" """
Path to a CSiBORG IC realisation snapshot. Path to an IC realisation snapshot.
Parameters Parameters
---------- ----------
nsnap : int nsnap : int
Snapshot index. Snapshot index. For Quijote, `-1` indicates the IC snapshot.
nsim : int nsim : int
IC realisation index. IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
Returns Returns
------- -------
snappath : str snappath : str
""" """
tonew = nsnap == 1 simpath = self.snapshots(nsim, simname, tonew=nsnap == 1)
simpath = self.ic_path(nsim, tonew=tonew) if simname == "csiborg":
return join(simpath, f"output_{str(nsnap).zfill(5)}") return join(simpath, f"output_{str(nsnap).zfill(5)}")
else:
if nsnap == -1:
return join(simpath, "ICs", "ics")
nsnap = str(nsnap).zfill(3)
return join(simpath, f"snapdir_{nsnap}", f"snap_{nsnap}")
def structfit(self, nsnap, nsim): def particles(self, nsim, simname):
""" """
Path to the clump or halo catalogue from `fit_halos.py`. Only CSiBORG Path to the files containing all particles of a CSiBORG realisation at
is supported. :math:`z = 0`.
Parameters
----------
nsim : int
IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
Returns
-------
path : str
"""
assert simname in ["csiborg", "quijote"]
if simname == "csiborg":
fdir = join(self.postdir, "particles")
if not isdir(fdir):
makedirs(fdir)
warn(f"Created directory `{fdir}`.", UserWarning)
fname = f"parts_{str(nsim).zfill(5)}.h5"
return join(fdir, fname)
fdir = join(self.quijote_dir, "Particles_fiducial")
if not isdir(fdir):
makedirs(fdir)
warn(f"Created directory `{fdir}`.", UserWarning)
fname = f"parts_{str(nsim).zfill(5)}.h5"
return join(fdir, fname)
def structfit(self, nsnap, nsim, simname):
"""
Path to the halo catalogue from `fit_halos.py`.
Parameters Parameters
---------- ----------
@ -371,15 +444,26 @@ class Paths:
Snapshot index. Snapshot index.
nsim : int nsim : int
IC realisation index. IC realisation index.
simname : str
Simulation name. Must be one of `csiborg` or `quijote`.
Returns Returns
------- -------
path : str path : str
""" """
fdir = join(self.postdir, "structfit") assert simname in ["csiborg", "quijote"]
if simname == "csiborg":
fdir = join(self.postdir, "structfit")
if not isdir(fdir):
mkdir(fdir)
warn(f"Created directory `{fdir}`.", UserWarning, stacklevel=1)
fname = f"out_{str(nsim).zfill(5)}_{str(nsnap).zfill(5)}.npy"
return join(fdir, fname)
fdir = join(self.quijote_dir, "structfit")
if not isdir(fdir): if not isdir(fdir):
mkdir(fdir) mkdir(fdir)
warn(f"Created directory `{fdir}`.", UserWarning, stacklevel=1) warn(f"Created directory `{fdir}`.", UserWarning)
fname = f"out_{str(nsim).zfill(5)}_{str(nsnap).zfill(5)}.npy" fname = f"out_{str(nsim).zfill(5)}_{str(nsnap).zfill(5)}.npy"
return join(fdir, fname) return join(fdir, fname)
@ -409,27 +493,6 @@ class Paths:
fname = fname.replace("overlap", "overlap_smoothed") fname = fname.replace("overlap", "overlap_smoothed")
return join(fdir, fname) return join(fdir, fname)
def particles(self, nsim):
"""
Path to the files containing all particles of a CSiBORG realisation at
:math:`z = 0`.
Parameters
----------
nsim : int
IC realisation index.
Returns
-------
path : str
"""
fdir = join(self.postdir, "particles")
if not isdir(fdir):
makedirs(fdir)
warn(f"Created directory `{fdir}`.", UserWarning, stacklevel=1)
fname = f"parts_{str(nsim).zfill(5)}.h5"
return join(fdir, fname)
def field(self, kind, MAS, grid, nsim, in_rsp, smooth_scale=None): def field(self, kind, MAS, grid, nsim, in_rsp, smooth_scale=None):
r""" r"""
Path to the files containing the calculated density fields in CSiBORG. Path to the files containing the calculated density fields in CSiBORG.

337
csiborgtools/read/readsim.py
View file

@ -15,34 +15,28 @@
""" """
Functions to read in the particle and clump files. Functions to read in the particle and clump files.
""" """
from abc import ABC, abstractmethod
from datetime import datetime
from gc import collect
from os.path import isfile, join from os.path import isfile, join
from warnings import warn from warnings import warn
import numpy import numpy
import readfof
import readgadget
from scipy.io import FortranFile from scipy.io import FortranFile
from tqdm import tqdm, trange from tqdm import tqdm, trange
from .paths import Paths from .paths import Paths
from .utils import cols_to_structured from .utils import cols_to_structured
###############################################################################
# Fortran particle reader #
###############################################################################
class BaseReader(ABC):
class ParticleReader:
""" """
Object to read in particle files along with their corresponding haloes. Base class for all readers.
Parameters
----------
paths : py:class`csiborgtools.read.Paths`
""" """
_paths = None _paths = None
def __init__(self, paths):
self.paths = paths
@property @property
def paths(self): def paths(self):
""" """
@ -59,9 +53,10 @@ class ParticleReader:
assert isinstance(paths, Paths) assert isinstance(paths, Paths)
self._paths = paths self._paths = paths
@abstractmethod
def read_info(self, nsnap, nsim): def read_info(self, nsnap, nsim):
""" """
Read CSiBORG simulation snapshot info. Read simulation snapshot info.
Parameters Parameters
---------- ----------
@ -73,10 +68,63 @@ class ParticleReader:
Returns Returns
------- -------
info : dict info : dict
Dictionary of information paramaters. Note that both keys and Dictionary of information paramaters.
values are strings.
""" """
snappath = self.paths.snapshot(nsnap, nsim) pass
@abstractmethod
def read_particle(self, nsnap, nsim, pars_extract, return_structured=True,
verbose=True):
"""
Read particle files of a simulation at a given snapshot and return
values of `pars_extract`.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
pars_extract : list of str
Parameters to be extracted.
return_structured : bool, optional
Whether to return a structured array or a 2-dimensional array. If
the latter, then the order of the columns is the same as the order
of `pars_extract`. However, enforces single-precision floating
point format for all columns.
verbose : bool, optional
Verbosity flag while for reading in the files.
Returns
-------
out : structured array or 2-dimensional array
Particle information.
pids : 1-dimensional array
Particle IDs.
"""
pass
###############################################################################
# CSiBORG particle reader #
###############################################################################
class CSiBORGReader:
"""
Object to read in CSiBORG snapshots from the binary files and halo
catalogues.
Parameters
----------
paths : py:class`csiborgtools.read.Paths`
"""
def __init__(self, paths):
self.paths = paths
def read_info(self, nsnap, nsim):
snappath = self.paths.snapshot(nsnap, nsim, "csiborg")
filename = join(snappath, "info_{}.txt".format(str(nsnap).zfill(5))) filename = join(snappath, "info_{}.txt".format(str(nsnap).zfill(5)))
with open(filename, "r") as f: with open(filename, "r") as f:
info = f.read().split() info = f.read().split()
@ -87,7 +135,7 @@ class ParticleReader:
keys = info[eqs - 1] keys = info[eqs - 1]
vals = info[eqs + 1] vals = info[eqs + 1]
return {key: val for key, val in zip(keys, vals)} return {key: convert_str_to_num(val) for key, val in zip(keys, vals)}
def open_particle(self, nsnap, nsim, verbose=True): def open_particle(self, nsnap, nsim, verbose=True):
""" """
@ -109,7 +157,7 @@ class ParticleReader:
partfiles : list of `scipy.io.FortranFile` partfiles : list of `scipy.io.FortranFile`
Opened part files. Opened part files.
""" """
snappath = self.paths.snapshot(nsnap, nsim) snappath = self.paths.snapshot(nsnap, nsim, "csiborg")
ncpu = int(self.read_info(nsnap, nsim)["ncpu"]) ncpu = int(self.read_info(nsnap, nsim)["ncpu"])
nsnap = str(nsnap).zfill(5) nsnap = str(nsnap).zfill(5)
if verbose: if verbose:
@ -192,33 +240,6 @@ class ParticleReader:
def read_particle(self, nsnap, nsim, pars_extract, return_structured=True, def read_particle(self, nsnap, nsim, pars_extract, return_structured=True,
verbose=True): verbose=True):
"""
Read particle files of a simulation at a given snapshot and return
values of `pars_extract`.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
pars_extract : list of str
Parameters to be extracted.
return_structured : bool, optional
Whether to return a structured array or a 2-dimensional array. If
the latter, then the order of the columns is the same as the order
of `pars_extract`. However, enforces single-precision floating
point format for all columns.
verbose : bool, optional
Verbosity flag while for reading the CPU outputs.
Returns
-------
out : structured array or 2-dimensional array
Particle information.
pids : 1-dimensional array
Particle IDs.
"""
# Open the particle files # Open the particle files
nparts, partfiles = self.open_particle(nsnap, nsim, verbose=verbose) nparts, partfiles = self.open_particle(nsnap, nsim, verbose=verbose)
if verbose: if verbose:
@ -299,11 +320,11 @@ class ParticleReader:
""" """
nsnap = str(nsnap).zfill(5) nsnap = str(nsnap).zfill(5)
cpu = str(cpu + 1).zfill(5) cpu = str(cpu + 1).zfill(5)
fpath = join(self.paths.ic_path(nsim, tonew=False), f"output_{nsnap}", fpath = join(self.paths.snapshots(nsim, "csiborg", tonew=False),
f"unbinding_{nsnap}.out{cpu}") f"output_{nsnap}", f"unbinding_{nsnap}.out{cpu}")
return FortranFile(fpath) return FortranFile(fpath)
def read_clumpid(self, nsnap, nsim, verbose=True): def read_phew_clumpid(self, nsnap, nsim, verbose=True):
""" """
Read PHEW clump IDs of particles from unbinding files. This halo finder Read PHEW clump IDs of particles from unbinding files. This halo finder
was used when running the catalogue. was used when running the catalogue.
@ -337,7 +358,7 @@ class ParticleReader:
return clumpid return clumpid
def read_clumps(self, nsnap, nsim, cols=None): def read_phew_clups(self, nsnap, nsim, cols=None):
""" """
Read in a PHEW clump file `clump_xxXXX.dat`. Read in a PHEW clump file `clump_xxXXX.dat`.
@ -356,7 +377,7 @@ class ParticleReader:
out : structured array out : structured array
""" """
nsnap = str(nsnap).zfill(5) nsnap = str(nsnap).zfill(5)
fname = join(self.paths.ic_path(nsim, tonew=False), fname = join(self.paths.snapshots(nsim, "csiborg", tonew=False),
"output_{}".format(nsnap), "output_{}".format(nsnap),
"clump_{}.dat".format(nsnap)) "clump_{}.dat".format(nsnap))
if not isfile(fname): if not isfile(fname):
@ -387,7 +408,7 @@ class ParticleReader:
out[col] = data[:, clump_cols[col][0]] out[col] = data[:, clump_cols[col][0]]
return out return out
def read_fof_hids(self, nsim): def read_fof_hids(self, nsim, **kwargs):
""" """
Read in the FoF particle halo membership IDs that are sorted to match Read in the FoF particle halo membership IDs that are sorted to match
the PHEW output. the PHEW output.
@ -396,13 +417,16 @@ class ParticleReader:
---------- ----------
nsim : int nsim : int
IC realisation index. IC realisation index.
**kwargs : dict
Keyword arguments for backward compatibility.
Returns Returns
------- -------
hids : 1-dimensional array hids : 1-dimensional array
Halo IDs of particles. Halo IDs of particles.
""" """
return numpy.load(self.paths.fof_membership(nsim, sorted=True)) return numpy.load(self.paths.fof_membership(nsim, "csiborg",
sorted=True))
def read_fof_halos(self, nsim): def read_fof_halos(self, nsim):
""" """
@ -417,7 +441,7 @@ class ParticleReader:
------- -------
cat : structured array cat : structured array
""" """
fpath = self.paths.fof_cat(nsim) fpath = self.paths.fof_cat(nsim, "csiborg")
hid = numpy.genfromtxt(fpath, usecols=0, dtype=numpy.int32) hid = numpy.genfromtxt(fpath, usecols=0, dtype=numpy.int32)
pos = numpy.genfromtxt(fpath, usecols=(1, 2, 3), dtype=numpy.float32) pos = numpy.genfromtxt(fpath, usecols=(1, 2, 3), dtype=numpy.float32)
totmass = numpy.genfromtxt(fpath, usecols=4, dtype=numpy.float32) totmass = numpy.genfromtxt(fpath, usecols=4, dtype=numpy.float32)
@ -437,7 +461,7 @@ class ParticleReader:
############################################################################### ###############################################################################
# Summed substructure catalogue # # Summed substructure PHEW catalogue for CSiBORG #
############################################################################### ###############################################################################
@ -467,8 +491,8 @@ class MmainReader:
Parameters Parameters
---------- ----------
clumparr : structured array clumparr : structured array
Clump array. Read from `ParticleReader.read_clumps`. Must contain Clump array. Read from `CSiBORGReader.read_phew_clups`. Must
`index` and `parent` columns. contain `index` and `parent` columns.
verbose : bool, optional verbose : bool, optional
Verbosity flag. Verbosity flag.
@ -522,10 +546,10 @@ class MmainReader:
The ultimate parent halo index for every clump, i.e. referring to The ultimate parent halo index for every clump, i.e. referring to
its ultimate parent clump. its ultimate parent clump.
""" """
nsnap = max(self.paths.get_snapshots(nsim)) nsnap = max(self.paths.get_snapshots(nsim, "csiborg"))
partreader = ParticleReader(self.paths) partreader = CSiBORGReader(self.paths)
cols = ["index", "parent", "mass_cl", 'x', 'y', 'z'] cols = ["index", "parent", "mass_cl", 'x', 'y', 'z']
clumparr = partreader.read_clumps(nsnap, nsim, cols) clumparr = partreader.read_phew_clups(nsnap, nsim, cols)
ultimate_parent = self.find_parents(clumparr, verbose=verbose) ultimate_parent = self.find_parents(clumparr, verbose=verbose)
mask_main = clumparr["index"] == clumparr["parent"] mask_main = clumparr["index"] == clumparr["parent"]
@ -548,37 +572,162 @@ class MmainReader:
out["subfrac"] = 1 - clumparr["mass_cl"][mask_main] / out["M"] out["subfrac"] = 1 - clumparr["mass_cl"][mask_main] / out["M"]
return out, ultimate_parent return out, ultimate_parent
############################################################################### ###############################################################################
# Supplementary reading functions # # Quijote particle reader #
############################################################################### ###############################################################################
def read_initcm(nsim, srcdir, fname="clump_{}_cm.npy"): class QuijoteReader:
""" """
Read `clump_cm`, i.e. the center of mass of a clump at redshift z = 70. Object to read in Quijote snapshots from the binary files.
If the file does not exist returns `None`.
Parameters Parameters
---------- ----------
nsim : int paths : py:class`csiborgtools.read.Paths`
IC realisation index.
srcdir : str
Path to the folder containing the files.
fname : str, optional
File name convention. By default `clump_cm_{}.npy`, where the
substituted value is `nsim`.
Returns
-------
out : structured array
""" """
fpath = join(srcdir, fname.format(nsim))
try: def __init__(self, paths):
return numpy.load(fpath) self.paths = paths
except FileNotFoundError:
warn("File {} does not exist.".format(fpath), UserWarning, def read_info(self, nsnap, nsim):
stacklevel=1) snapshot = self.paths.snapshot(nsnap, nsim, "quijote")
return None header = readgadget.header(snapshot)
out = {"BoxSize": header.boxsize / 1e3, # Mpc/h
"Nall": header.nall[1], # Tot num of particles
"PartMass": header.massarr[1] * 1e10, # Part mass in Msun/h
"Omega_m": header.omega_m,
"Omega_l": header.omega_l,
"h": header.hubble,
"redshift": header.redshift,
}
out["TotMass"] = out["Nall"] * out["PartMass"]
out["Hubble"] = (100.0 * numpy.sqrt(
header.omega_m * (1.0 + header.redshift)**3 + header.omega_l))
return out
def read_particle(self, nsnap, nsim, pars_extract=None,
return_structured=True, verbose=True):
assert pars_extract in [None, "pids"]
snapshot = self.paths.snapshot(nsnap, nsim, "quijote")
info = self.read_info(nsnap, nsim)
ptype = [1] # DM in Gadget speech
if verbose:
print(f"{datetime.now()}: reading particle IDs.")
pids = readgadget.read_block(snapshot, "ID ", ptype)
if pars_extract == "pids":
return None, pids
if return_structured:
dtype = {"names": ['x', 'y', 'z', 'vx', 'vy', 'vz', 'M'],
"formats": [numpy.float32] * 7}
out = numpy.full(info["Nall"], numpy.nan, dtype=dtype)
else:
out = numpy.full((info["Nall"], 7), numpy.nan, dtype=numpy.float32)
if verbose:
print(f"{datetime.now()}: reading particle positions.")
pos = readgadget.read_block(snapshot, "POS ", ptype) / 1e3 # Mpc/h
pos /= info["BoxSize"] # Box units
for i, p in enumerate(['x', 'y', 'z']):
if return_structured:
out[p] = pos[:, i]
else:
out[:, i] = pos[:, i]
del pos
collect()
if verbose:
print(f"{datetime.now()}: reading particle velocities.")
# NOTE convert to box units.
vel = readgadget.read_block(snapshot, "VEL ", ptype) # km/s
vel *= (1 + info["redshift"])
for i, v in enumerate(['vx', 'vy', 'vz']):
if return_structured:
out[v] = vel[:, i]
else:
out[:, i + 3] = vel[:, i]
del vel
collect()
if verbose:
print(f"{datetime.now()}: reading particle masses.")
if return_structured:
out["M"] = info["PartMass"] / info["TotMass"]
else:
out[:, 6] = info["PartMass"] / info["TotMass"]
return out, pids
def read_fof_hids(self, nsnap, nsim, verbose=True, **kwargs):
"""
Read the FoF group membership of particles. Unassigned particles have
FoF group ID 0.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
verbose : bool, optional
Verbosity flag.
**kwargs : dict
Keyword arguments for backward compatibility.
Returns
-------
out : 1-dimensional array of shape `(nparticles, )`
Group membership of particles.
"""
redshift = {4: 0.0, 3: 0.5, 2: 1.0, 1: 2.0, 0: 3.0}.get(nsnap, None)
if redshift is None:
raise ValueError(f"Redshift of snapshot {nsnap} is not known.")
path = self.paths.fof_cat(nsim, "quijote")
cat = readfof.FoF_catalog(path, nsnap)
# Read the particle IDs of the snapshot
__, pids = self.read_particle(nsnap, nsim, pars_extract="pids",
verbose=verbose)
# Read the FoF particle membership. These are only assigned particles.
if verbose:
print(f"{datetime.now()}: reading the FoF particle membership.",
flush=True)
group_pids = cat.GroupIDs
group_len = cat.GroupLen
# Create a mapping from particle ID to FoF group ID.
if verbose:
print(f"{datetime.now()}: creating the particle to FoF ID to map.",
flush=True)
ks = numpy.insert(numpy.cumsum(group_len), 0, 0)
pid2hid = numpy.full((group_pids.size, 2), numpy.nan,
dtype=numpy.uint32)
for i, (k0, kf) in enumerate(zip(ks[:-1], ks[1:])):
pid2hid[k0:kf, 0] = i + 1
pid2hid[k0:kf, 1] = group_pids[k0:kf]
pid2hid = {pid: hid for hid, pid in pid2hid}
# Create the final array of hids matchign the snapshot array.
# Unassigned particles have hid 0.
if verbose:
print(f"{datetime.now()}: creating the final hid array.",
flush=True)
hids = numpy.full(pids.size, 0, dtype=numpy.uint32)
for i in trange(pids.size) if verbose else range(pids.size):
hids[i] = pid2hid.get(pids[i], 0)
return hids
###############################################################################
# Supplementary reading functions #
###############################################################################
def halfwidth_mask(pos, hw): def halfwidth_mask(pos, hw):
@ -627,3 +776,27 @@ def load_halo_particles(hid, particles, halo_map, hid2map):
return particles[k0:kf + 1, :] return particles[k0:kf + 1, :]
except KeyError: except KeyError:
return None return None
def convert_str_to_num(s):
"""
Convert a string representation of a number to its appropriate numeric type
(int or float).
Parameters
----------
s : str
The string representation of the number.
Returns
-------
num : int or float
"""
try:
return int(s)
except ValueError:
try:
return float(s)
except ValueError:
warn(f"Cannot convert string '{s}' to number", UserWarning)
return s

6
scripts/cluster_crosspk.py
View file

@ -51,7 +51,7 @@ MAS = "CIC" # mass asignment scheme
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
box = csiborgtools.read.CSiBORGBox(paths) box = csiborgtools.read.CSiBORGBox(paths)
reader = csiborgtools.read.ParticleReader(paths) reader = csiborgtools.read.CSiBORGReader(paths)
ics = paths.get_ics("csiborg") ics = paths.get_ics("csiborg")
nsims = len(ics) nsims = len(ics)
@ -66,8 +66,8 @@ jobs = csiborgtools.utils.split_jobs(nsims, nproc)[rank]
for n in jobs: for n in jobs:
print(f"Rank {rank} at {datetime.now()}: saving {n}th delta.", flush=True) print(f"Rank {rank} at {datetime.now()}: saving {n}th delta.", flush=True)
nsim = ics[n] nsim = ics[n]
particles = reader.read_particle(max(paths.get_snapshots(nsim)), nsim, particles = reader.read_particle(max(paths.get_snapshots(nsim, "csiborg")),
["x", "y", "z", "M"], verbose=False) nsim, ["x", "y", "z", "M"], verbose=False)
# Halfwidth -- particle selection # Halfwidth -- particle selection
if args.halfwidth < 0.5: if args.halfwidth < 0.5:
particles = csiborgtools.read.halfwidth_select( particles = csiborgtools.read.halfwidth_select(

22
scripts/field_prop.py
View file

@ -58,9 +58,10 @@ def density_field(nsim, parser_args, to_save=True):
field : 3-dimensional array field : 3-dimensional array
""" """
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, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
parts = csiborgtools.read.read_h5(paths.particles(nsim))["particles"] parts = csiborgtools.read.read_h5(paths.particles(nsim, "csiborg"))
parts = parts["particles"]
gen = csiborgtools.field.DensityField(box, parser_args.MAS) gen = csiborgtools.field.DensityField(box, parser_args.MAS)
if parser_args.kind == "density": if parser_args.kind == "density":
@ -114,9 +115,10 @@ def velocity_field(nsim, parser_args, to_save=True):
"Smoothed velocity field is not implemented.") "Smoothed velocity field is not implemented.")
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
mpart = 1.1641532e-10 # Particle mass in CSiBORG simulations. mpart = 1.1641532e-10 # Particle mass in CSiBORG simulations.
nsnap = max(paths.get_snapshots(nsim)) nsnap = max(paths.get_snapshots(nsim, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
parts = csiborgtools.read.read_h5(paths.particles(nsim))["particles"] parts = csiborgtools.read.read_h5(paths.particles(nsim, "csiborg"))
parts = parts["particles"]
gen = csiborgtools.field.VelocityField(box, parser_args.MAS) gen = csiborgtools.field.VelocityField(box, parser_args.MAS)
field = gen(parts, parser_args.grid, mpart, verbose=parser_args.verbose) field = gen(parts, parser_args.grid, mpart, verbose=parser_args.verbose)
@ -152,7 +154,7 @@ def potential_field(nsim, parser_args, to_save=True):
potential : 3-dimensional array potential : 3-dimensional array
""" """
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, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
# Load the real space overdensity field # Load the real space overdensity field
@ -168,7 +170,8 @@ def potential_field(nsim, parser_args, to_save=True):
field = gen(rho) field = gen(rho)
if parser_args.in_rsp: if parser_args.in_rsp:
parts = csiborgtools.read.read_h5(paths.particles(nsim))["particles"] parts = csiborgtools.read.read_h5(paths.particles(nsim, "csiborg"))
parts = parts["particles"]
field = csiborgtools.field.field2rsp(*field, parts=parts, box=box, field = csiborgtools.field.field2rsp(*field, parts=parts, box=box,
verbose=parser_args.verbose) verbose=parser_args.verbose)
if to_save: if to_save:
@ -207,7 +210,7 @@ def radvel_field(nsim, parser_args, to_save=True):
raise NotImplementedError( raise NotImplementedError(
"Smoothed radial vel. field not implemented.") "Smoothed radial vel. field not implemented.")
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, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
vel = numpy.load(paths.field("velocity", parser_args.MAS, parser_args.grid, vel = numpy.load(paths.field("velocity", parser_args.MAS, parser_args.grid,
@ -245,7 +248,7 @@ def environment_field(nsim, parser_args, to_save=True):
env : 3-dimensional array env : 3-dimensional array
""" """
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, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
density_gen = csiborgtools.field.DensityField(box, parser_args.MAS) density_gen = csiborgtools.field.DensityField(box, parser_args.MAS)
gen = csiborgtools.field.TidalTensorField(box, parser_args.MAS) gen = csiborgtools.field.TidalTensorField(box, parser_args.MAS)
@ -268,7 +271,8 @@ def environment_field(nsim, parser_args, to_save=True):
# Optionally drag the field to RSP. # Optionally drag the field to RSP.
if parser_args.in_rsp: if parser_args.in_rsp:
parts = csiborgtools.read.read_h5(paths.particles(nsim))["particles"] parts = csiborgtools.read.read_h5(paths.particles(nsim, "csiborg"))
parts = parts["particles"]
fields = (tensor_field.T00, tensor_field.T11, tensor_field.T22, fields = (tensor_field.T00, tensor_field.T11, tensor_field.T22,
tensor_field.T01, tensor_field.T02, tensor_field.T12) tensor_field.T01, tensor_field.T02, tensor_field.T12)

23
scripts/fit_halos.py
View file

@ -81,8 +81,8 @@ def _main(nsim, simname, verbose):
verbose : bool verbose : bool
Verbosity flag. Verbosity flag.
""" """
if simname == "quijote": # if simname == "quijote":
raise NotImplementedError("Quijote not implemented yet.") # raise NotImplementedError("Quijote not implemented yet.")
cols = [("index", numpy.int32), cols = [("index", numpy.int32),
("npart", numpy.int32), ("npart", numpy.int32),
@ -95,17 +95,22 @@ def _main(nsim, simname, verbose):
("m200c", numpy.float32), ("m200c", numpy.float32),
("lambda200c", numpy.float32),] ("lambda200c", numpy.float32),]
nsnap = max(paths.get_snapshots(nsim)) nsnap = max(paths.get_snapshots(nsim, simname))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) if simname == "csiborg":
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
cat = csiborgtools.read.CSiBORGHaloCatalogue(
nsim, paths, with_lagpatch=False, load_initial=False, rawdata=True,
load_fitted=False)
else:
box = csiborgtools.read.QuijoteBox(nsnap, nsim, paths)
cat = csiborgtools.read.QuijoteHaloCatalogue(
nsim, paths, nsnap, load_initial=False, rawdata=True)
# Particle archive # Particle archive
f = csiborgtools.read.read_h5(paths.particles(nsim)) f = csiborgtools.read.read_h5(paths.particles(nsim, simname))
particles = f["particles"] particles = f["particles"]
halo_map = f["halomap"] halo_map = f["halomap"]
hid2map = {hid: i for i, hid in enumerate(halo_map[:, 0])} hid2map = {hid: i for i, hid in enumerate(halo_map[:, 0])}
cat = csiborgtools.read.CSiBORGHaloCatalogue(
nsim, paths, with_lagpatch=False, load_initial=False, rawdata=True,
load_fitted=False)
out = csiborgtools.read.cols_to_structured(len(cat), cols) out = csiborgtools.read.cols_to_structured(len(cat), cols)
for i in trange(len(cat)) if verbose else range(len(cat)): for i in trange(len(cat)) if verbose else range(len(cat)):
@ -121,7 +126,7 @@ def _main(nsim, simname, verbose):
for key in _out.keys(): for key in _out.keys():
out[key][i] = _out[key] out[key][i] = _out[key]
fout = paths.structfit(nsnap, nsim) fout = paths.structfit(nsnap, nsim, simname)
if verbose: if verbose:
print(f"Saving to `{fout}`.", flush=True) print(f"Saving to `{fout}`.", flush=True)
numpy.save(fout, out) numpy.save(fout, out)

25
scripts/fit_init.py
View file

@ -50,9 +50,6 @@ def _main(nsim, simname, verbose):
verbose : bool verbose : bool
Verbosity flag. Verbosity flag.
""" """
if simname == "quijote":
raise NotImplementedError("Quijote not implemented yet.")
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
cols = [("index", numpy.int32), cols = [("index", numpy.int32),
("x", numpy.float32), ("x", numpy.float32),
@ -61,15 +58,26 @@ def _main(nsim, simname, verbose):
("lagpatch_size", numpy.float32), ("lagpatch_size", numpy.float32),
("lagpatch_ncells", numpy.int32),] ("lagpatch_ncells", numpy.int32),]
parts = csiborgtools.read.read_h5(paths.initmatch(nsim, "particles")) fname = paths.initmatch(nsim, simname, "particles")
parts = csiborgtools.read.read_h5(fname)
parts = parts['particles'] parts = parts['particles']
halo_map = csiborgtools.read.read_h5(paths.particles(nsim)) halo_map = csiborgtools.read.read_h5(paths.particles(nsim, simname))
halo_map = halo_map["halomap"] halo_map = halo_map["halomap"]
cat = csiborgtools.read.CSiBORGHaloCatalogue( if simname == "csiborg":
nsim, paths, rawdata=True, load_fitted=False, load_initial=False) cat = csiborgtools.read.CSiBORGHaloCatalogue(
nsim, paths, rawdata=True, load_fitted=False, load_initial=False)
else:
cat = csiborgtools.read.QuijoteHaloCatalogue(nsim, paths, nsnap=4)
hid2map = {hid: i for i, hid in enumerate(halo_map[:, 0])} hid2map = {hid: i for i, hid in enumerate(halo_map[:, 0])}
# Initialise the overlapper.
if simname == "csiborg":
kwargs = {"box_size": 2048, "bckg_halfsize": 475}
else:
kwargs = {"box_size": 512, "bckg_halfsize": 256}
overlapper = csiborgtools.match.ParticleOverlap(**kwargs)
out = csiborgtools.read.cols_to_structured(len(cat), cols) out = csiborgtools.read.cols_to_structured(len(cat), cols)
for i, hid in enumerate(tqdm(cat["index"]) if verbose else cat["index"]): for i, hid in enumerate(tqdm(cat["index"]) if verbose else cat["index"]):
out["index"][i] = hid out["index"][i] = hid
@ -88,12 +96,11 @@ def _main(nsim, simname, verbose):
out["lagpatch_size"][i] = numpy.percentile(distances, 99) out["lagpatch_size"][i] = numpy.percentile(distances, 99)
# Calculate the number of cells with > 0 density. # Calculate the number of cells with > 0 density.
overlapper = csiborgtools.match.ParticleOverlap()
delta = overlapper.make_delta(pos, mass, subbox=True) delta = overlapper.make_delta(pos, mass, subbox=True)
out["lagpatch_ncells"][i] = csiborgtools.fits.delta2ncells(delta) out["lagpatch_ncells"][i] = csiborgtools.fits.delta2ncells(delta)
# Now save it # Now save it
fout = paths.initmatch(nsim, "fit") fout = paths.initmatch(nsim, simname, "fit")
if verbose: if verbose:
print(f"{datetime.now()}: dumping fits to .. `{fout}`.", flush=True) print(f"{datetime.now()}: dumping fits to .. `{fout}`.", flush=True)
with open(fout, "wb") as f: with open(fout, "wb") as f:

15
scripts/match_singlematch.py
View file

@ -30,12 +30,15 @@ except ModuleNotFoundError:
def pair_match(nsim0, nsimx, sigma, smoothen, verbose): def pair_match(nsim0, nsimx, sigma, smoothen, verbose):
# TODO fix this.
simname = "csiborg"
overlapper_kwargs = {"box_size": 512, "bckg_halfsize": 475}
from csiborgtools.read import CSiBORGHaloCatalogue, read_h5 from csiborgtools.read import CSiBORGHaloCatalogue, read_h5
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
smooth_kwargs = {"sigma": sigma, "mode": "constant", "cval": 0.0} smooth_kwargs = {"sigma": sigma, "mode": "constant", "cval": 0.0}
overlapper = csiborgtools.match.ParticleOverlap() overlapper = csiborgtools.match.ParticleOverlap(**overlapper_kwargs)
matcher = csiborgtools.match.RealisationsMatcher() matcher = csiborgtools.match.RealisationsMatcher(**overlapper_kwargs)
# Load the raw catalogues (i.e. no selection) including the initial CM # Load the raw catalogues (i.e. no selection) including the initial CM
# positions and the particle archives. # positions and the particle archives.
@ -45,12 +48,12 @@ def pair_match(nsim0, nsimx, sigma, smoothen, verbose):
catx = CSiBORGHaloCatalogue(nsimx, paths, load_initial=True, bounds=bounds, catx = CSiBORGHaloCatalogue(nsimx, paths, load_initial=True, bounds=bounds,
with_lagpatch=True, load_clumps_cat=True) with_lagpatch=True, load_clumps_cat=True)
clumpmap0 = read_h5(paths.particles(nsim0))["clumpmap"] clumpmap0 = read_h5(paths.particles(nsim0, simname))["clumpmap"]
parts0 = read_h5(paths.initmatch(nsim0, "particles"))["particles"] parts0 = read_h5(paths.initmatch(nsim0, simname, "particles"))["particles"]
clid2map0 = {clid: i for i, clid in enumerate(clumpmap0[:, 0])} clid2map0 = {clid: i for i, clid in enumerate(clumpmap0[:, 0])}
clumpmapx = read_h5(paths.particles(nsimx))["clumpmap"] clumpmapx = read_h5(paths.particles(nsimx, simname))["clumpmap"]
partsx = read_h5(paths.initmatch(nsimx, "particles"))["particles"] partsx = read_h5(paths.initmatch(nsimx, simname, "particles"))["particles"]
clid2mapx = {clid: i for i, clid in enumerate(clumpmapx[:, 0])} clid2mapx = {clid: i for i, clid in enumerate(clumpmapx[:, 0])}
# We generate the background density fields. Loads halos's particles one by # We generate the background density fields. Loads halos's particles one by

12
scripts/mv_fofmembership.py
View file

@ -57,7 +57,7 @@ def copy_membership(nsim, verbose=True):
print(f"Loading from ... `{fpath}`.") print(f"Loading from ... `{fpath}`.")
data = numpy.genfromtxt(fpath, dtype=int) data = numpy.genfromtxt(fpath, dtype=int)
fout = paths.fof_membership(nsim) fout = paths.fof_membership(nsim, "csiborg")
if verbose: if verbose:
print(f"Saving to ... `{fout}`.") print(f"Saving to ... `{fout}`.")
numpy.save(fout, data) numpy.save(fout, data)
@ -77,7 +77,7 @@ def copy_catalogue(nsim, verbose=True):
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
source = join("/mnt/extraspace/jeg/greenwhale/Constrained_Sims", source = join("/mnt/extraspace/jeg/greenwhale/Constrained_Sims",
f"sim_{nsim}/halo_catalog_{nsim}_FOF.txt") f"sim_{nsim}/halo_catalog_{nsim}_FOF.txt")
dest = paths.fof_cat(nsim) dest = paths.fof_cat(nsim, "csiborg")
if verbose: if verbose:
print("Copying`{}` to `{}`.".format(source, dest)) print("Copying`{}` to `{}`.".format(source, dest))
copy(source, dest) copy(source, dest)
@ -96,14 +96,14 @@ def sort_fofid(nsim, verbose=True):
Verbosity flag. Verbosity flag.
""" """
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, "csiborg"))
fpath = paths.fof_membership(nsim) fpath = paths.fof_membership(nsim, "csiborg")
if verbose: if verbose:
print(f"{datetime.now()}: loading from ... `{fpath}`.") print(f"{datetime.now()}: loading from ... `{fpath}`.")
# Columns are halo ID, particle ID. # Columns are halo ID, particle ID.
fof = numpy.load(fpath) fof = numpy.load(fpath)
reader = csiborgtools.read.ParticleReader(paths) reader = csiborgtools.read.CSiBORGReader(paths)
pars_extract = ["x"] # Dummy variable pars_extract = ["x"] # Dummy variable
__, pids = reader.read_particle(nsnap, nsim, pars_extract, __, pids = reader.read_particle(nsnap, nsim, pars_extract,
return_structured=False, verbose=verbose) return_structured=False, verbose=verbose)
@ -123,7 +123,7 @@ def sort_fofid(nsim, verbose=True):
hid, pid = fof[i] hid, pid = fof[i]
fof_hids[pids_idx[pid]] = hid fof_hids[pids_idx[pid]] = hid
fout = paths.fof_membership(nsim, sorted=True) fout = paths.fof_membership(nsim, "csiborg", sorted=True)
if verbose: if verbose:
print(f"Saving the sorted data to ... `{fout}`") print(f"Saving the sorted data to ... `{fout}`")
numpy.save(fout, fof_hids) numpy.save(fout, fof_hids)

4
scripts/old/fit_profiles.py
View file

@ -58,10 +58,10 @@ for i, nsim in enumerate(nsims):
if rank == 0: if rank == 0:
now = datetime.now() now = datetime.now()
print(f"{now}: calculating {i}th simulation `{nsim}`.", flush=True) print(f"{now}: calculating {i}th simulation `{nsim}`.", flush=True)
nsnap = max(paths.get_snapshots(nsim)) nsnap = max(paths.get_snapshots(nsim, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
f = csiborgtools.read.read_h5(paths.particles(nsim)) f = csiborgtools.read.read_h5(paths.particles(nsim, "csiborg"))
particles = f["particles"] particles = f["particles"]
clump_map = f["clumpmap"] clump_map = f["clumpmap"]
clid2map = {clid: i for i, clid in enumerate(clump_map[:, 0])} clid2map = {clid: i for i, clid in enumerate(clump_map[:, 0])}

2
scripts/old/pre_mmain.py
View file

@ -38,7 +38,7 @@ mmain_reader = csiborgtools.read.MmainReader(paths)
def do_mmain(nsim): def do_mmain(nsim):
nsnap = max(paths.get_snapshots(nsim)) nsnap = max(paths.get_snapshots(nsim, "csiborg"))
# NOTE: currently works for highest snapshot anyway # NOTE: currently works for highest snapshot anyway
mmain, ultimate_parent = mmain_reader.make_mmain(nsim, verbose=False) mmain, ultimate_parent = mmain_reader.make_mmain(nsim, verbose=False)
numpy.savez(paths.mmain(nsnap, nsim), numpy.savez(paths.mmain(nsnap, nsim),

36
scripts/pre_dumppart.py
View file

@ -60,6 +60,11 @@ def minmax_halo(hid, halo_ids, start_loop=0):
return start, end return start, end
###############################################################################
# Sorting and dumping #
###############################################################################
def main(nsim, simname, verbose): def main(nsim, simname, verbose):
""" """
Read in the snapshot particles, sort them by their FoF halo ID and dump Read in the snapshot particles, sort them by their FoF halo ID and dump
@ -81,20 +86,21 @@ def main(nsim, simname, verbose):
None None
""" """
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
partreader = csiborgtools.read.ParticleReader(paths) if simname == "csiborg":
partreader = csiborgtools.read.CSiBORGReader(paths)
else:
partreader = csiborgtools.read.QuijoteReader(paths)
if simname == "quijote": nsnap = max(paths.get_snapshots(nsim, simname))
raise NotImplementedError("Not implemented for Quijote yet.") fname = paths.particles(nsim, simname)
# Keep "ID" as the last column!
pars_extract = ['x', 'y', 'z', 'vx', 'vy', 'vz', 'M', "ID"]
nsnap = max(paths.get_snapshots(nsim))
fname = paths.particles(nsim)
# We first read in the halo IDs of the particles and infer the sorting. # We first read in the halo IDs of the particles and infer the sorting.
# Right away we dump the halo IDs to a HDF5 file and clear up memory. # Right away we dump the halo IDs to a HDF5 file and clear up memory.
if verbose: if verbose:
print(f"{datetime.now()}: loading particles {nsim}.", flush=True) print(f"{datetime.now()}: loading PIDs of IC {nsim}.", flush=True)
part_hids = partreader.read_fof_hids(nsim) part_hids = partreader.read_fof_hids(
nsnap=nsnap, nsim=nsim, verbose=verbose)
if verbose:
print(f"{datetime.now()}: sorting PIDs of IC {nsim}.", flush=True)
sort_indxs = numpy.argsort(part_hids).astype(numpy.int32) sort_indxs = numpy.argsort(part_hids).astype(numpy.int32)
part_hids = part_hids[sort_indxs] part_hids = part_hids[sort_indxs]
with h5py.File(fname, "w") as f: with h5py.File(fname, "w") as f:
@ -106,6 +112,10 @@ def main(nsim, simname, verbose):
# Next we read in the particles and sort them by their halo ID. # Next we read in the particles and sort them by their halo ID.
# We cannot directly read this as an unstructured array because the float32 # We cannot directly read this as an unstructured array because the float32
# precision is insufficient to capture the halo IDs. # precision is insufficient to capture the halo IDs.
if simname == "csiborg":
pars_extract = ['x', 'y', 'z', 'vx', 'vy', 'vz', 'M', "ID"]
else:
pars_extract = None
parts, pids = partreader.read_particle( parts, pids = partreader.read_particle(
nsnap, nsim, pars_extract, return_structured=False, verbose=verbose) nsnap, nsim, pars_extract, return_structured=False, verbose=verbose)
# Now we in two steps save the particles and particle IDs. # Now we in two steps save the particles and particle IDs.
@ -129,11 +139,11 @@ def main(nsim, simname, verbose):
collect() collect()
if verbose: if verbose:
print(f"{datetime.now()}: creating halo map for {nsim}.", flush=True) print(f"{datetime.now()}: creating a halo map for {nsim}.", flush=True)
# Load clump IDs back to memory # Load clump IDs back to memory
with h5py.File(fname, "r") as f: with h5py.File(fname, "r") as f:
part_hids = f["halo_ids"][:] part_hids = f["halo_ids"][:]
# We loop over the unique clump IDs. # We loop over the unique halo IDs.
unique_halo_ids = numpy.unique(part_hids) unique_halo_ids = numpy.unique(part_hids)
halo_map = numpy.full((unique_halo_ids.size, 3), numpy.nan, halo_map = numpy.full((unique_halo_ids.size, 3), numpy.nan,
dtype=numpy.int32) dtype=numpy.int32)
@ -148,7 +158,7 @@ def main(nsim, simname, verbose):
start_loop = kf start_loop = kf
# We save the mapping to a HDF5 file # We save the mapping to a HDF5 file
with h5py.File(paths.particles(nsim), "r+") as f: with h5py.File(fname, "r+") as f:
f.create_dataset("halomap", data=halo_map) f.create_dataset("halomap", data=halo_map)
f.close() f.close()

43
scripts/pre_sortinit.py
View file

@ -50,34 +50,55 @@ def _main(nsim, simname, verbose):
verbose : bool verbose : bool
Verbosity flag. Verbosity flag.
""" """
if simname == "quijote":
raise NotImplementedError("Quijote not implemented yet.")
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
partreader = csiborgtools.read.ParticleReader(paths) if simname == "csiborg":
partreader = csiborgtools.read.CSiBORGReader(paths)
else:
partreader = csiborgtools.read.QuijoteReader(paths)
if verbose: if verbose:
print(f"{datetime.now()}: reading and processing simulation {nsim}.", print(f"{datetime.now()}: reading and processing simulation `{nsim}`.",
flush=True) flush=True)
# We first load the particle IDs in the final snapshot. # We first load the particle IDs in the final snapshot.
pidf = csiborgtools.read.read_h5(paths.particles(nsim)) pidf = csiborgtools.read.read_h5(paths.particles(nsim, simname))
pidf = pidf["particle_ids"] pidf = pidf["particle_ids"]
# Then we load the particles in the initil snapshot and make sure that # Then we load the particles in the initil snapshot and make sure that
# their particle IDs are sorted as in the final snapshot. Again, because of # their particle IDs are sorted as in the final snapshot. Again, because of
# precision this must be read as structured. # precision this must be read as structured.
# NOTE: ID has to be the last column. if simname == "csiborg":
pars_extract = ["x", "y", "z", "M", "ID"] pars_extract = ["x", "y", "z", "M", "ID"]
# CSiBORG's initial snapshot ID
nsnap = 1
else:
pars_extract = None
# Use this to point the reader to the ICs snapshot
nsnap = -1
part0, pid0 = partreader.read_particle( part0, pid0 = partreader.read_particle(
1, nsim, pars_extract, return_structured=False, verbose=verbose) nsnap, nsim, pars_extract, return_structured=False, verbose=verbose)
# Quijote's initial snapshot information also contains velocities but we
# don't need those.
if simname == "quijote":
part0 = part0[:, [0, 1, 2, 6]]
# In Quijote some particles are position precisely at the edge of the
# box. Move them to be just inside.
pos = part0[:, :3]
mask = pos >= 1
if numpy.any(mask):
spacing = numpy.spacing(pos[mask])
assert numpy.max(spacing) <= 1e-5
pos[mask] -= spacing
# First enforce them to already be sorted and then apply reverse # First enforce them to already be sorted and then apply reverse
# sorting from the final snapshot. # sorting from the final snapshot.
part0 = part0[numpy.argsort(pid0)] part0 = part0[numpy.argsort(pid0)]
del pid0 del pid0
collect() collect()
part0 = part0[numpy.argsort(numpy.argsort(pidf))] part0 = part0[numpy.argsort(numpy.argsort(pidf))]
fout = paths.initmatch(nsim, simname, "particles")
if verbose: if verbose:
print(f"{datetime.now()}: dumping particles for {nsim}.", flush=True) print(f"{datetime.now()}: dumping particles for `{nsim}` to `{fout}`",
with h5py.File(paths.initmatch(nsim, "particles"), "w") as f: flush=True)
with h5py.File(fout, "w") as f:
f.create_dataset("particles", data=part0) f.create_dataset("particles", data=part0)

9
scripts_plots/plot_data.py
View file

@ -352,7 +352,7 @@ def plot_projected_field(kind, nsim, grid, in_rsp, smooth_scale, MAS="PCS",
""" """
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, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
if kind == "overdensity": if kind == "overdensity":
@ -563,10 +563,10 @@ def plot_sky_distribution(field, nsim, grid, nside, smooth_scale=None,
Whether to save the figure as a pdf. 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, "csiborg"))
box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths) box = csiborgtools.read.CSiBORGBox(nsnap, nsim, paths)
if kind == "overdensity": if field== "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) smooth_scale=smooth_scale)
density_gen = csiborgtools.field.DensityField(box, MAS) density_gen = csiborgtools.field.DensityField(box, MAS)
@ -658,7 +658,8 @@ if __name__ == "__main__":
if False: if False:
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
d = csiborgtools.read.read_h5(paths.particles(7444))["particles"] d = csiborgtools.read.read_h5(paths.particles(7444, "csiborg"))
d = d["particles"]
plt.figure() plt.figure()
plt.hist(d[:100000, 4], bins="auto") plt.hist(d[:100000, 4], bins="auto")