mirror of
https://github.com/Richard-Sti/csiborgtools_public.git
synced 2025-06-08 09:51:12 +00:00
Calculate upglade redshifts (#128)
* Update redshift reading * Add helio to CMB redshift * Update imports * Update nb * Run for Quijote * Add script * Update * Update .gitignore * Update imports * Add Peery estimator * Add bulk flow scripts * Update typs * Add comment * Add blank space * Update submission script * Update description * Add barriers * Update nb * Update nb * Rename script * Move to old * Update imports * Add nb * Update script * Fix catalogue key * Update script * Update submit * Update comment * Update .gitignore * Update nb * Update for stationary obsrevers * Update submission * Add nb * Add better verbose control * Update nb * Update submit * Update nb * Add SN errors * Add draft of the script * Update verbosity flags * Add submission script * Debug script * Quickfix * Remove comment * Update nb * Update submission * Update nb * Processed UPGLADE
This commit is contained in:
Richard Stiskalek
GitHub
parent
c447d2e7b0
commit
779f2e76ac
24 changed files with 2066 additions and 300 deletions
|
|
@ -20,7 +20,8 @@ from .utils import (center_of_mass, delta2ncells, number_counts,
|
|||
binned_statistic, cosine_similarity, fprint, # noqa
|
||||
hms_to_degrees, dms_to_degrees, great_circle_distance, # noqa
|
||||
radec_to_cartesian, cartesian_to_radec, # noqa
|
||||
thin_samples_by_acl, numpyro_gof, radec_to_galactic) # noqa
|
||||
thin_samples_by_acl, numpyro_gof, radec_to_galactic, # noqa
|
||||
heliocentric_to_cmb, calculate_acl) # noqa
|
||||
from .params import (paths_glamdring, simname2boxsize, simname2Omega_m, # noqa
|
||||
snap2redshift) # noqa
|
||||
|
||||
|
|
|
|||
|
|
@ -16,8 +16,9 @@ from .density import (DensityField, PotentialField, TidalTensorField,
|
|||
VelocityField, radial_velocity, power_spectrum, # noqa
|
||||
overdensity_field) # noqa
|
||||
from .enclosed_mass import (particles_enclosed_mass, # noqa
|
||||
particles_enclosed_momentum, field_enclosed_mass) # noqa
|
||||
from .interp import (evaluate_cartesian_cic, evaluate_sky, evaluate_los, # noqa
|
||||
particles_enclosed_momentum, field_enclosed_mass, # noqa
|
||||
bulkflow_peery2018) # noqa
|
||||
from .interp import (evaluate_cartesian_cic, evaluate_sky, evaluate_los, # noqa
|
||||
field2rsp, fill_outside, make_sky, # noqa
|
||||
observer_peculiar_velocity, smoothen_field, # noqa
|
||||
field_at_distance) # noqa
|
||||
|
|
|
|||
|
|
@ -12,10 +12,13 @@
|
|||
# You should have received a copy of the GNU General Public License along
|
||||
# with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
|
||||
|
||||
import numpy
|
||||
"""
|
||||
Code to calculate the enclosed mass, momentum or bulk flow from various
|
||||
radial-velocity based estimators (that may be flawed, be careful).
|
||||
"""
|
||||
import numpy as np
|
||||
from numba import jit
|
||||
from tqdm import tqdm
|
||||
|
||||
|
||||
###############################################################################
|
||||
|
|
@ -55,7 +58,7 @@ def particles_enclosed_mass(rdist, mass, distances):
|
|||
enclosed_mass : 1-dimensional array
|
||||
Enclosed mass at each distance.
|
||||
"""
|
||||
enclosed_mass = numpy.full_like(distances, 0.)
|
||||
enclosed_mass = np.full_like(distances, 0.)
|
||||
start_index = 0
|
||||
for i, dist in enumerate(distances):
|
||||
if i > 0:
|
||||
|
|
@ -121,9 +124,10 @@ def field_enclosed_mass(field, distances, boxsize):
|
|||
enclosed_volume : 1-dimensional array
|
||||
Enclosed grid-like volume at each distance.
|
||||
"""
|
||||
enclosed_mass = numpy.zeros_like(distances)
|
||||
enclosed_volume = numpy.zeros_like(distances)
|
||||
for i, dist in enumerate(distances):
|
||||
enclosed_mass = np.zeros_like(distances)
|
||||
enclosed_volume = np.zeros_like(distances)
|
||||
|
||||
for i, dist in enumerate(tqdm(distances)):
|
||||
enclosed_mass[i], enclosed_volume[i] = _field_enclosed_mass(
|
||||
field, dist, boxsize)
|
||||
|
||||
|
|
@ -137,7 +141,7 @@ def field_enclosed_mass(field, distances, boxsize):
|
|||
|
||||
@jit(nopython=True, boundscheck=False)
|
||||
def _enclosed_momentum(rdist, mass, vel, rmax, start_index):
|
||||
bulk_momentum = numpy.zeros(3, dtype=rdist.dtype)
|
||||
bulk_momentum = np.zeros(3, dtype=rdist.dtype)
|
||||
|
||||
for i in range(start_index, len(rdist)):
|
||||
if rdist[i] <= rmax:
|
||||
|
|
@ -169,7 +173,7 @@ def particles_enclosed_momentum(rdist, mass, vel, distances):
|
|||
bulk_momentum : 2-dimensional array
|
||||
Enclosed momentum at each distance.
|
||||
"""
|
||||
bulk_momentum = numpy.zeros((len(distances), 3))
|
||||
bulk_momentum = np.zeros((len(distances), 3))
|
||||
start_index = 0
|
||||
for i, dist in enumerate(distances):
|
||||
if i > 0:
|
||||
|
|
@ -179,4 +183,68 @@ def particles_enclosed_momentum(rdist, mass, vel, distances):
|
|||
start_index)
|
||||
bulk_momentum[i] += v
|
||||
|
||||
return bulk_momentum
|
||||
return bulk_momentum
|
||||
|
||||
|
||||
###############################################################################
|
||||
# Bulk flow estimators #
|
||||
###############################################################################
|
||||
|
||||
|
||||
def bulkflow_peery2018(rdist, mass, pos, vel, distances, weights,
|
||||
verbose=True):
|
||||
"""
|
||||
Calculate the bulk flow from a set of particles using the estimator from
|
||||
Peery+2018. Supports either `1/r^2` or `constant` weights. Particles
|
||||
are assumed to be sorted by distance from the center of the box.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rdist : 1-dimensional array
|
||||
Sorted distance of particles from the center of the box.
|
||||
mass : 1-dimensional array
|
||||
Sorted mass of particles.
|
||||
pos : 2-dimensional array
|
||||
Sorted position of particles.
|
||||
vel : 2-dimensional array
|
||||
Sorted velocity of particles.
|
||||
distances : 1-dimensional array
|
||||
Distances at which to calculate the bulk flow.
|
||||
weights : str
|
||||
Weights to use in the estimator, either `1/r^2` or `constant`.
|
||||
verbose : bool
|
||||
Verbosity flag.
|
||||
|
||||
Returns
|
||||
-------
|
||||
bulk_flow : 2-dimensional array
|
||||
"""
|
||||
# Select only the particles within the maximum distance to speed up the
|
||||
# calculation.
|
||||
if verbose:
|
||||
print("Selecting particles within the maximum distance...")
|
||||
kmax = np.searchsorted(rdist, np.max(distances))
|
||||
rdist = rdist[:kmax]
|
||||
mass = mass[:kmax]
|
||||
pos = pos[:kmax]
|
||||
vel = vel[:kmax]
|
||||
|
||||
if verbose:
|
||||
print("Computing the cumulative quantities...")
|
||||
if weights == "1/r^2":
|
||||
cumulative_x = np.cumsum(mass[:, np.newaxis] * np.sum(vel * pos, axis=1)[:, np.newaxis] * pos / rdist[:, np.newaxis]**4, axis=0) # noqa
|
||||
norm = lambda R: R**2 # noqa
|
||||
elif weights == "constant":
|
||||
cumulative_x = np.cumsum(mass[:, np.newaxis] * np.sum(vel * pos, axis=1)[:, np.newaxis] * pos / rdist[:, np.newaxis]**2, axis=0) # noqa
|
||||
norm = lambda R: 3. # noqa
|
||||
else:
|
||||
raise ValueError("Invalid weights.")
|
||||
cumulative_x /= np.cumsum(mass)[:, np.newaxis]
|
||||
|
||||
B = np.zeros((len(distances), 3))
|
||||
for i in range(3):
|
||||
for j, R in enumerate(distances):
|
||||
k = np.searchsorted(rdist, R)
|
||||
B[j, i] = norm(R) * cumulative_x[k - 1, i]
|
||||
|
||||
return B
|
||||
|
|
|
|||
|
|
@ -48,7 +48,6 @@ H0 = 100 # km / s / Mpc
|
|||
|
||||
|
||||
def t():
|
||||
"""Shortcut to get the current time."""
|
||||
return datetime.now().strftime("%H:%M:%S")
|
||||
|
||||
|
||||
|
|
@ -79,21 +78,26 @@ class DataLoader:
|
|||
store_full_velocity : bool, optional
|
||||
Whether to store the full 3D velocity field. Otherwise stores only
|
||||
the radial velocity.
|
||||
verbose : bool, optional
|
||||
Verbose flag.
|
||||
"""
|
||||
def __init__(self, simname, ksim, catalogue, catalogue_fpath, paths,
|
||||
ksmooth=None, store_full_velocity=False):
|
||||
print(f"{t()}: reading the catalogue.")
|
||||
ksmooth=None, store_full_velocity=False, verbose=True):
|
||||
if verbose:
|
||||
print(f"{t()}: reading the catalogue.", flush=True)
|
||||
self._cat = self._read_catalogue(catalogue, catalogue_fpath)
|
||||
self._catname = catalogue
|
||||
|
||||
print(f"{t()}: reading the interpolated field.")
|
||||
if verbose:
|
||||
print(f"{t()}: reading the interpolated field.", flush=True)
|
||||
self._field_rdist, self._los_density, self._los_velocity = self._read_field( # noqa
|
||||
simname, ksim, catalogue, ksmooth, paths)
|
||||
|
||||
if len(self._field_rdist) % 2 == 0:
|
||||
warn(f"The number of radial steps is even. Skipping the first "
|
||||
f"step at {self._field_rdist[0]} because Simpson's rule "
|
||||
"requires an odd number of steps.")
|
||||
if verbose:
|
||||
warn(f"The number of radial steps is even. Skipping the first "
|
||||
f"step at {self._field_rdist[0]} because Simpson's rule "
|
||||
"requires an odd number of steps.")
|
||||
self._field_rdist = self._field_rdist[1:]
|
||||
self._los_density = self._los_density[..., 1:]
|
||||
self._los_velocity = self._los_velocity[..., 1:]
|
||||
|
|
@ -102,7 +106,8 @@ class DataLoader:
|
|||
raise ValueError("The number of objects in the catalogue does not "
|
||||
"match the number of objects in the field.")
|
||||
|
||||
print(f"{t()}: calculating the radial velocity.")
|
||||
if verbose:
|
||||
print(f"{t()}: calculating the radial velocity.", flush=True)
|
||||
nobject = len(self._los_density)
|
||||
dtype = self._los_density.dtype
|
||||
|
||||
|
|
@ -231,9 +236,7 @@ class DataLoader:
|
|||
return rdist, los_density, los_velocity
|
||||
|
||||
def _read_catalogue(self, catalogue, catalogue_fpath):
|
||||
"""
|
||||
Read in the distance indicator catalogue.
|
||||
"""
|
||||
"""Read in the distance indicator catalogue."""
|
||||
if catalogue == "A2":
|
||||
with File(catalogue_fpath, 'r') as f:
|
||||
dtype = [(key, np.float32) for key in f.keys()]
|
||||
|
|
@ -256,10 +259,18 @@ class DataLoader:
|
|||
arr[key] = grp[key][:]
|
||||
elif "CB2_" in catalogue:
|
||||
with File(catalogue_fpath, 'r') as f:
|
||||
|
||||
dtype = [(key, np.float32) for key in f.keys()]
|
||||
arr = np.empty(len(f["RA"]), dtype=dtype)
|
||||
for key in f.keys():
|
||||
arr[key] = f[key][:]
|
||||
elif "UPGLADE" in catalogue:
|
||||
with File(catalogue_fpath, 'r') as f:
|
||||
dtype = [(key, np.float32) for key in f.keys()]
|
||||
arr = np.empty(len(f["RA"]), dtype=dtype)
|
||||
for key in f.keys():
|
||||
if key == "mask":
|
||||
continue
|
||||
|
||||
arr[key] = f[key][:]
|
||||
else:
|
||||
raise ValueError(f"Unknown catalogue: `{catalogue}`.")
|
||||
|
|
@ -916,6 +927,8 @@ class SN_PV_validation_model(BaseFlowValidationModel):
|
|||
Right ascension and declination in degrees.
|
||||
z_obs : 1-dimensional array of shape (n_objects)
|
||||
Observed redshifts.
|
||||
e_zobs : 1-dimensional array of shape (n_objects)
|
||||
Errors on the observed redshifts.
|
||||
mB, x1, c : 1-dimensional arrays of shape (n_objects)
|
||||
SALT2 light curve parameters.
|
||||
e_mB, e_x1, e_c : 1-dimensional arrays of shape (n_objects)
|
||||
|
|
@ -927,7 +940,7 @@ class SN_PV_validation_model(BaseFlowValidationModel):
|
|||
"""
|
||||
|
||||
def __init__(self, los_density, los_velocity, RA, dec, z_obs,
|
||||
mB, x1, c, e_mB, e_x1, e_c, r_xrange, Omega_m):
|
||||
e_zobs, mB, x1, c, e_mB, e_x1, e_c, r_xrange, Omega_m):
|
||||
dt = jnp.float32
|
||||
# Convert everything to JAX arrays.
|
||||
self._los_density = jnp.asarray(los_density, dtype=dt)
|
||||
|
|
@ -936,6 +949,11 @@ class SN_PV_validation_model(BaseFlowValidationModel):
|
|||
self._RA = jnp.asarray(np.deg2rad(RA), dtype=dt)
|
||||
self._dec = jnp.asarray(np.deg2rad(dec), dtype=dt)
|
||||
self._z_obs = jnp.asarray(z_obs, dtype=dt)
|
||||
if e_zobs is not None:
|
||||
self._e2_cz_obs = jnp.asarray(
|
||||
(SPEED_OF_LIGHT * e_zobs)**2, dtype=dt)
|
||||
else:
|
||||
self._e2_cz_obs = jnp.zeros_like(self._z_obs)
|
||||
|
||||
self._mB = jnp.asarray(mB, dtype=dt)
|
||||
self._x1 = jnp.asarray(x1, dtype=dt)
|
||||
|
|
@ -1121,6 +1139,7 @@ class SN_PV_validation_model(BaseFlowValidationModel):
|
|||
alpha_cal = numpyro.sample("alpha_cal", self._alpha_cal)
|
||||
beta_cal = numpyro.sample("beta_cal", self._beta_cal)
|
||||
|
||||
cz_err = jnp.sqrt(sigma_v**2 + self._e2_cz_obs)
|
||||
Vext_rad = project_Vext(Vx, Vy, Vz, self._RA, self._dec)
|
||||
|
||||
mu = self.mu(mag_cal, alpha_cal, beta_cal)
|
||||
|
|
@ -1137,7 +1156,7 @@ class SN_PV_validation_model(BaseFlowValidationModel):
|
|||
# Calculate p(z_obs) and multiply it by p(r)
|
||||
zobs_pred = self._f_zobs(beta, Vext_rad[i], self._los_velocity[i])
|
||||
ptilde *= calculate_likelihood_zobs(
|
||||
self._z_obs[i], zobs_pred, sigma_v)
|
||||
self._z_obs[i], zobs_pred, cz_err[i])
|
||||
|
||||
return ll + jnp.log(self._f_simps(ptilde) / pnorm), None
|
||||
|
||||
|
|
@ -1385,17 +1404,19 @@ def get_model(loader, zcmb_max=None, verbose=True):
|
|||
if kind in ["LOSS", "Foundation"]:
|
||||
keys = ["RA", "DEC", "z_CMB", "mB", "x1", "c", "e_mB", "e_x1", "e_c"]
|
||||
RA, dec, zCMB, mB, x1, c, e_mB, e_x1, e_c = (loader.cat[k] for k in keys) # noqa
|
||||
e_zCMB = None
|
||||
|
||||
mask = (zCMB < zcmb_max)
|
||||
model = SN_PV_validation_model(
|
||||
los_overdensity[mask], los_velocity[mask], RA[mask], dec[mask],
|
||||
zCMB[mask], mB[mask], x1[mask], c[mask], e_mB[mask], e_x1[mask],
|
||||
e_c[mask], loader.rdist, loader._Omega_m)
|
||||
zCMB[mask], e_zCMB, mB[mask], x1[mask], c[mask], e_mB[mask],
|
||||
e_x1[mask], e_c[mask], loader.rdist, loader._Omega_m)
|
||||
elif "Pantheon+" in kind:
|
||||
keys = ["RA", "DEC", "zCMB", "mB", "x1", "c", "biasCor_m_b", "mBERR",
|
||||
"x1ERR", "cERR", "biasCorErr_m_b", "zCMB_SN", "zCMB_Group"]
|
||||
"x1ERR", "cERR", "biasCorErr_m_b", "zCMB_SN", "zCMB_Group",
|
||||
"zCMBERR"]
|
||||
|
||||
RA, dec, zCMB, mB, x1, c, bias_corr_mB, e_mB, e_x1, e_c, e_bias_corr_mB, zCMB_SN, zCMB_Group = (loader.cat[k] for k in keys) # noqa
|
||||
RA, dec, zCMB, mB, x1, c, bias_corr_mB, e_mB, e_x1, e_c, e_bias_corr_mB, zCMB_SN, zCMB_Group, e_zCMB = (loader.cat[k] for k in keys) # noqa
|
||||
mB -= bias_corr_mB
|
||||
e_mB = np.sqrt(e_mB**2 + e_bias_corr_mB**2)
|
||||
|
||||
|
|
@ -1413,8 +1434,8 @@ def get_model(loader, zcmb_max=None, verbose=True):
|
|||
|
||||
model = SN_PV_validation_model(
|
||||
los_overdensity[mask], los_velocity[mask], RA[mask], dec[mask],
|
||||
zCMB[mask], mB[mask], x1[mask], c[mask], e_mB[mask], e_x1[mask],
|
||||
e_c[mask], loader.rdist, loader._Omega_m)
|
||||
zCMB[mask], e_zCMB[mask], mB[mask], x1[mask], c[mask], e_mB[mask],
|
||||
e_x1[mask], e_c[mask], loader.rdist, loader._Omega_m)
|
||||
elif kind in ["SFI_gals", "2MTF", "SFI_gals_masked"]:
|
||||
keys = ["RA", "DEC", "z_CMB", "mag", "eta", "e_mag", "e_eta"]
|
||||
RA, dec, zCMB, mag, eta, e_mag, e_eta = (loader.cat[k] for k in keys)
|
||||
|
|
|
|||
|
|
@ -15,7 +15,7 @@
|
|||
from .catalogue import (CSiBORG1Catalogue, CSiBORG2Catalogue, # noqa
|
||||
CSiBORG2SUBFINDCatalogue, # noqa
|
||||
CSiBORG2MergerTreeReader, QuijoteCatalogue, # noqa
|
||||
MDPL2Catalogue) # noqa
|
||||
MDPL2Catalogue, fiducial_observers) # noqa
|
||||
from .snapshot import (CSiBORG1Snapshot, CSiBORG2Snapshot, QuijoteSnapshot, # noqa
|
||||
CSiBORG1Field, CSiBORG2Field, QuijoteField, BORG2Field, # noqa
|
||||
BORG1Field, TNG300_1Field) # noqa
|
||||
|
|
|
|||
|
|
@ -1372,7 +1372,7 @@ class QuijoteCatalogue(BaseCatalogue):
|
|||
|
||||
@property
|
||||
def totmass(self):
|
||||
return self._read_fof_catalogue("group_mass")
|
||||
return self._read_fof_catalogue("GroupMass")
|
||||
|
||||
@property
|
||||
def index(self):
|
||||
|
|
|
|||
|
|
@ -413,6 +413,47 @@ def real2redshift(pos, vel, observer_location, observer_velocity, boxsize,
|
|||
return pos
|
||||
|
||||
|
||||
def heliocentric_to_cmb(z_helio, RA, dec, e_z_helio=None):
|
||||
"""
|
||||
Convert heliocentric redshift to CMB redshift using the Planck 2018 CMB
|
||||
dipole.
|
||||
"""
|
||||
# CMB dipole Planck 2018 values
|
||||
vsun_mag = 369 # km/s
|
||||
RA_sun = 167.942
|
||||
dec_sun = -6.944
|
||||
SPEED_OF_LIGHT = 299792.458 # km / s
|
||||
|
||||
theta_sun = np.pi / 2 - np.deg2rad(dec_sun)
|
||||
phi_sun = np.deg2rad(RA_sun)
|
||||
|
||||
# Convert to theat/phi in radians
|
||||
theta = np.pi / 2 - np.deg2rad(dec)
|
||||
phi = np.deg2rad(RA)
|
||||
|
||||
# Unit vector in the direction of each galaxy
|
||||
n = np.asarray([np.sin(theta) * np.cos(phi),
|
||||
np.sin(theta) * np.sin(phi),
|
||||
np.cos(theta)]).T
|
||||
# CMB dipole unit vector
|
||||
vsun_normvect = np.asarray([np.sin(theta_sun) * np.cos(phi_sun),
|
||||
np.sin(theta_sun) * np.sin(phi_sun),
|
||||
np.cos(theta_sun)])
|
||||
|
||||
# Project the CMB dipole onto the line of sight and normalize
|
||||
vsun_projected = vsun_mag * np.dot(n, vsun_normvect) / SPEED_OF_LIGHT
|
||||
|
||||
zsun_tilde = np.sqrt((1 - vsun_projected) / (1 + vsun_projected))
|
||||
zcmb = (1 + z_helio) / zsun_tilde - 1
|
||||
|
||||
# Optional linear error propagation
|
||||
if e_z_helio is not None:
|
||||
e_zcmb = np.abs(e_z_helio / zsun_tilde)
|
||||
return zcmb, e_zcmb
|
||||
|
||||
return zcmb
|
||||
|
||||
|
||||
###############################################################################
|
||||
# Statistics #
|
||||
###############################################################################
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue