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

More flow preparation & Olympics (#143)

Browse source 
* Add more comments

* Add flow paths

* Simplify paths

* Update default arguemnts

* Update paths

* Update param names

* Update some of scipts for reading files

* Add the Mike method option

* Update plotting

* Update fnames

* Simplify things

* Make more default options

* Add print

* Update

* Downsample CF4

* Update numpyro selection

* Add selection fitting nb

* Add coeffs

* Update script

* Add nb

* Add label

* Increase number of steps

* Update default params

* Add more labels

* Improve file name

* Update nb

* Fix little bug

* Remove import

* Update scales

* Update labels

* Add script

* Update script

* Add more

* Add more labels

* Add script

* Add submit

* Update spacing

* Update submit scrips

* Update script

* Update defaults

* Update defaults

* Update nb

* Update test

* Update imports

* Add script

* Add support for Indranil void

* Add a dipole

* Update nb

* Update submit

* Update Om0

* Add final

* Update default params

* Fix bug

* Add option to fix to LG frame

* Add Vext label

* Add Vext label

* Update script

* Rm fixed LG

* rm LG stuff

* Update script

* Update bulk flow plotting

* Update nb

* Add no field option

* Update defaults

* Update nb

* Update script

* Update nb

* Update nb

* Add names to plots

* Update nb

* Update plot

* Add more latex names

* Update default

* Update nb

* Update np

* Add plane slicing

* Add nb with slices

* Update nb

* Update script

* Upddate nb

* Update nb
This commit is contained in:
Richard Stiskalek 2024-09-11 08:45:42 +02:00 committed by GitHub
parent 3d1e1c0ae3
commit 2b938c112c
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
20 changed files with 3106 additions and 379 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
csiborgtools/field/__init__.py
View file

@ -20,7 +20,7 @@ try:
VelocityField, radial_velocity, power_spectrum, # noqa VelocityField, radial_velocity, power_spectrum, # noqa
overdensity_field) # noqa overdensity_field) # noqa
from .interp import (evaluate_cartesian_cic, evaluate_los, field2rsp, # noqa from .interp import (evaluate_cartesian_cic, evaluate_los, field2rsp, # noqa
fill_outside, make_sky, # noqa fill_outside, make_sky, xy_supergalactic_slice, # noqa
observer_peculiar_velocity, smoothen_field, # noqa observer_peculiar_velocity, smoothen_field, # noqa
field_at_distance) # noqa field_at_distance) # noqa
except ImportError: except ImportError:

68
csiborgtools/field/interp.py
View file

@ -20,6 +20,8 @@ import numpy as np
import smoothing_library as SL import smoothing_library as SL
from numba import jit from numba import jit
from scipy.interpolate import RegularGridInterpolator from scipy.interpolate import RegularGridInterpolator
from astropy.coordinates import SkyCoord, Supergalactic, Galactic, ICRS
from astropy.coordinates import CartesianRepresentation
from tqdm import tqdm from tqdm import tqdm
from ..utils import periodic_wrap_grid, radec_to_cartesian from ..utils import periodic_wrap_grid, radec_to_cartesian
@ -351,6 +353,72 @@ def make_sky(field, angpos, rmax, dr, boxsize, interpolation_method="cic",
return finterp return finterp
###############################################################################
# Supergalactic plane slice #
###############################################################################
def xy_supergalactic_slice(field, boxsize, xmin, xmax, ngrid, field_frame,
z_value=0):
"""
Create a 2D slice of a scalar field in the x-y supergalactic plane.
Parameters
----------
field : 3-dimensional array of shape `(grid, grid, grid)`
Field to be interpolated.
boxsize : float
Box size in `Mpc / h`.
xmin, xmax : float
Minimum and maximum x and y values in supergalactic coordinates.
ngrid : int
Number of grid points along each axis.
field_frame : str
Frame of the field. Must be one of `galactic`, `supergalactic` or
`icrs`.
z_value : float, optional
Value of the z-coordinate in supergalactic coordinates.
Returns
-------
2-dimensional array of shape `(ngrid, ngrid)`
"""
# Coordinates of the 2D slice in supergalactic coordinates
xgrid = np.linspace(xmin, xmax, ngrid)
ygrid = np.copy(xgrid)
grid = np.stack(np.meshgrid(xgrid, ygrid))
grid = grid.reshape(2, -1).T
grid = np.hstack([grid, np.ones(ngrid**2).reshape(-1, 1) * z_value])
supergalactic_coord = SkyCoord(CartesianRepresentation(
grid[:, 0], grid[:, 1], grid[:, 2]), frame=Supergalactic)
# Create a Supergalactic SkyCoord object from Cartesian coordinates
if field_frame == "galactic":
original_frame = Galactic
elif field_frame == "supergalactic":
original_frame = Supergalactic
elif field_frame == "icrs":
original_frame = ICRS
else:
raise ValueError(f"Unknown field frame: {field_frame}")
# Convert to field frame
coords = supergalactic_coord.transform_to(original_frame).cartesian
pos = np.stack([coords.x, coords.y, coords.z]).value
pos = pos.T
# Convert to appropriate box units
pos /= boxsize
pos += 0.5
if np.any(pos <= 0) or np.any(pos >= 1):
raise ValueError("Some positions are outside the box.")
return evaluate_cartesian_cic(field, pos=pos).reshape(ngrid, ngrid)
############################################################################### ###############################################################################
# Average field at a radial distance # # Average field at a radial distance #
############################################################################### ###############################################################################

14
csiborgtools/field/utils.py
View file

@ -17,14 +17,14 @@ Utility functions used in the rest of the `field` module to avoid circular
imports. imports.
""" """
from numba import jit from numba import jit
import numpy import numpy as np
import healpy import healpy
def force_single_precision(x): def force_single_precision(x):
"""Attempt to convert an array `x` to float32.""" """Attempt to convert an array `x` to float32."""
if x.dtype != numpy.float32: if x.dtype != np.float32:
x = x.astype(numpy.float32) x = x.astype(np.float32)
return x return x
@ -46,10 +46,10 @@ def nside2radec(nside):
Generate RA [0, 360] deg and declination [-90, 90] deg for HEALPix pixel Generate RA [0, 360] deg and declination [-90, 90] deg for HEALPix pixel
centres at a given nside. centres at a given nside.
""" """
pixs = numpy.arange(healpy.nside2npix(nside)) pixs = np.arange(healpy.nside2npix(nside))
theta, phi = healpy.pix2ang(nside, pixs) theta, phi = healpy.pix2ang(nside, pixs)
ra = 180 / numpy.pi * phi ra = 180 / np.pi * phi
dec = 90 - 180 / numpy.pi * theta dec = 90 - 180 / np.pi * theta
return numpy.vstack([ra, dec]).T return np.vstack([ra, dec]).T

150
csiborgtools/flow/flow_model.py
View file

@ -22,6 +22,7 @@ References
[1] https://arxiv.org/abs/1912.09383. [1] https://arxiv.org/abs/1912.09383.
""" """
from abc import ABC, abstractmethod from abc import ABC, abstractmethod
from os.path import join
import numpy as np import numpy as np
from astropy import units as u from astropy import units as u
@ -100,7 +101,12 @@ class DataLoader:
d1, d2 = self._cat["RA"], self._cat["DEC"] d1, d2 = self._cat["RA"], self._cat["DEC"]
num_sims = len(self._los_density) num_sims = len(self._los_density)
radvel = np.empty((num_sims, nobject, len(self._field_rdist)), dtype) if "IndranilVoid" in simname:
self._los_radial_velocity = self._los_velocity
self._los_velocity = None
else:
radvel = np.empty(
(num_sims, nobject, len(self._field_rdist)), dtype)
for k in range(num_sims): for k in range(num_sims):
for i in range(nobject): for i in range(nobject):
radvel[k, i, :] = radial_velocity_los( radvel[k, i, :] = radial_velocity_los(
@ -182,6 +188,13 @@ class DataLoader:
if isinstance(ksims, int): if isinstance(ksims, int):
ksims = [ksims] ksims = [ksims]
# For no-field read in Carrick+2015 but then zero it.
if simname == "no_field":
simname = "Carrick2015"
to_wipe = True
else:
to_wipe = False
if not all(0 <= ksim < len(nsims) for ksim in ksims): if not all(0 <= ksim < len(nsims) for ksim in ksims):
raise ValueError(f"Invalid simulation index: `{ksims}`") raise ValueError(f"Invalid simulation index: `{ksims}`")
@ -189,6 +202,14 @@ class DataLoader:
fpath = paths.field_los(simname, "Pantheon+") fpath = paths.field_los(simname, "Pantheon+")
elif "CF4_TFR" in catalogue: elif "CF4_TFR" in catalogue:
fpath = paths.field_los(simname, "CF4_TFR") fpath = paths.field_los(simname, "CF4_TFR")
elif "IndranilVoid" in catalogue:
fdir = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_los" # noqa
if "exp" in catalogue:
fpath = join(fdir, "v_pec_EXP_IndranilVoid.dat")
elif "gauss" in catalogue:
fpath = join(fdir, "v_pec_GAUSS_IndranilVoid.dat")
else:
raise ValueError("Unknown `IndranilVoid` catalogue.")
else: else:
fpath = paths.field_los(simname, catalogue) fpath = paths.field_los(simname, catalogue)
@ -212,6 +233,10 @@ class DataLoader:
los_density = np.stack(los_density) los_density = np.stack(los_density)
los_velocity = np.stack(los_velocity) los_velocity = np.stack(los_velocity)
if to_wipe:
los_density = np.ones_like(los_density)
los_velocity = np.zeros_like(los_velocity)
return rdist, los_density, los_velocity return rdist, los_density, los_velocity
def _read_catalogue(self, catalogue, catalogue_fpath): def _read_catalogue(self, catalogue, catalogue_fpath):
@ -507,22 +532,18 @@ def e2_distmod_TFR(e2_mag, e2_eta, eta, b, c, e_mu_intrinsic):
def sample_TFR(e_mu_min, e_mu_max, a_mean, a_std, b_mean, b_std, def sample_TFR(e_mu_min, e_mu_max, a_mean, a_std, b_mean, b_std,
c_mean, c_std, alpha_min, alpha_max, sample_alpha, c_mean, c_std, alpha_min, alpha_max, sample_alpha,
sample_curvature, a_dipole_mean, a_dipole_std, sample_a_dipole, a_dipole_mean, a_dipole_std, sample_a_dipole, name):
name):
"""Sample Tully-Fisher calibration parameters.""" """Sample Tully-Fisher calibration parameters."""
e_mu = sample(f"e_mu_{name}", Uniform(e_mu_min, e_mu_max)) e_mu = sample(f"e_mu_{name}", Uniform(e_mu_min, e_mu_max))
a = sample(f"a_{name}", Normal(a_mean, a_std)) a = sample(f"a_{name}", Normal(a_mean, a_std))
if sample_a_dipole: if sample_a_dipole:
ax, ay, az = sample(f"a_dipole_{name}", Normal(0, 5).expand([3])) ax, ay, az = sample(f"a_dipole_{name}", Normal(a_dipole_mean, a_dipole_std).expand([3])) # noqa
else: else:
ax, ay, az = 0.0, 0.0, 0.0 ax, ay, az = 0.0, 0.0, 0.0
b = sample(f"b_{name}", Normal(b_mean, b_std)) b = sample(f"b_{name}", Normal(b_mean, b_std))
if sample_curvature:
c = sample(f"c_{name}", Normal(c_mean, c_std)) c = sample(f"c_{name}", Normal(c_mean, c_std))
else:
c = 0.0
alpha = sample_alpha_bias(name, alpha_min, alpha_max, sample_alpha) alpha = sample_alpha_bias(name, alpha_min, alpha_max, sample_alpha)
@ -571,8 +592,8 @@ def sample_calibration(Vext_min, Vext_max, Vmono_min, Vmono_max, beta_min,
beta_max, sigma_v_min, sigma_v_max, sample_Vmono, beta_max, sigma_v_min, sigma_v_max, sample_Vmono,
sample_beta): sample_beta):
"""Sample the flow calibration.""" """Sample the flow calibration."""
Vext = sample("Vext", Uniform(Vext_min, Vext_max).expand([3]))
sigma_v = sample("sigma_v", Uniform(sigma_v_min, sigma_v_max)) sigma_v = sample("sigma_v", Uniform(sigma_v_min, sigma_v_max))
Vext = sample("Vext", Uniform(Vext_min, Vext_max).expand([3]))
if sample_beta: if sample_beta:
beta = sample("beta", Uniform(beta_min, beta_max)) beta = sample("beta", Uniform(beta_min, beta_max))
@ -620,8 +641,8 @@ class PV_LogLikelihood(BaseFlowValidationModel):
Errors on the observed redshifts. Errors on the observed redshifts.
calibration_params: dict calibration_params: dict
Calibration parameters of each object. Calibration parameters of each object.
magmax_selection : float mag_selection : dict
Maximum magnitude selection if strict threshold. Magnitude selection parameters.
r_xrange : 1-dimensional array r_xrange : 1-dimensional array
Radial distances where the field was interpolated for each object. Radial distances where the field was interpolated for each object.
Omega_m : float Omega_m : float
@ -630,13 +651,11 @@ class PV_LogLikelihood(BaseFlowValidationModel):
Catalogue kind, either "TFR", "SN", or "simple". Catalogue kind, either "TFR", "SN", or "simple".
name : str name : str
Name of the catalogue. Name of the catalogue.
toy_selection : tuple of length 3, optional
Toy magnitude selection paramers `m1`, `m2` and `a`. Optional.
""" """
def __init__(self, los_density, los_velocity, RA, dec, z_obs, e_zobs, def __init__(self, los_density, los_velocity, RA, dec, z_obs, e_zobs,
calibration_params, maxmag_selection, r_xrange, Omega_m, calibration_params, mag_selection, r_xrange, Omega_m,
kind, name, toy_selection=None): kind, name):
if e_zobs is not None: if e_zobs is not None:
e2_cz_obs = jnp.asarray((SPEED_OF_LIGHT * e_zobs)**2) e2_cz_obs = jnp.asarray((SPEED_OF_LIGHT * e_zobs)**2)
else: else:
@ -657,8 +676,24 @@ class PV_LogLikelihood(BaseFlowValidationModel):
self.name = name self.name = name
self.Omega_m = Omega_m self.Omega_m = Omega_m
self.norm = - self.ndata * jnp.log(self.num_sims) self.norm = - self.ndata * jnp.log(self.num_sims)
self.maxmag_selection = maxmag_selection
self.toy_selection = toy_selection if mag_selection is not None:
self.mag_selection_kind = mag_selection["kind"]
if self.mag_selection_kind == "hard":
self.mag_selection_max = mag_selection["coeffs"]
fprint(f"catalogue {name} with selection mmax = {self.mag_selection_max}.") # noqa
elif self.mag_selection_kind == "soft":
self.m1, self.m2, self.a = mag_selection["coeffs"]
fprint(f"catalogue {name} with selection m1 = {self.m1}, m2 = {self.m2}, a = {self.a}.") # noqa
self.log_Fm = toy_log_magnitude_selection(
self.mag, self.m1, self.m2, self.a)
else:
self.mag_selection_kind = None
if mag_selection is not None and kind != "TFR":
raise ValueError("Magnitude selection is only implemented "
"for TFRs.")
if kind == "TFR": if kind == "TFR":
self.mag_min, self.mag_max = jnp.min(self.mag), jnp.max(self.mag) self.mag_min, self.mag_max = jnp.min(self.mag), jnp.max(self.mag)
@ -675,23 +710,13 @@ class PV_LogLikelihood(BaseFlowValidationModel):
else: else:
raise RuntimeError("Support most be added for other kinds.") raise RuntimeError("Support most be added for other kinds.")
if maxmag_selection is not None and self.maxmag_selection > self.mag_max: # noqa if self.mag_selection_kind == "hard" and self.mag_selection_max > self.mag_max: # noqa
raise ValueError("The maximum magnitude cannot be larger than the selection threshold.") # noqa raise ValueError("The maximum magnitude cannot be larger than "
"the selection threshold.")
if toy_selection is not None and self.maxmag_selection is not None:
raise ValueError("`toy_selection` and `maxmag_selection` cannot be used together.") # noqa
if toy_selection is not None:
self.m1, self.m2, self.a = toy_selection
self.log_Fm = toy_log_magnitude_selection(
self.mag, self.m1, self.m2, self.a)
if toy_selection is not None and self.kind != "TFR":
raise ValueError("Toy selection is only implemented for TFRs.")
def __call__(self, field_calibration_params, distmod_params, def __call__(self, field_calibration_params, distmod_params,
inference_method): inference_method):
if inference_method not in ["mike", "bayes"]: if inference_method not in ["mike", "bayes", "delta"]:
raise ValueError(f"Unknown method: `{inference_method}`.") raise ValueError(f"Unknown method: `{inference_method}`.")
ll0 = 0.0 ll0 = 0.0
@ -717,7 +742,7 @@ class PV_LogLikelihood(BaseFlowValidationModel):
"c", self.name, self.c_min, self.c_max) "c", self.name, self.c_min, self.c_max)
# NOTE: that the true variables are currently uncorrelated. # NOTE: that the true variables are currently uncorrelated.
with plate("true_SN", self.ndata): with plate(f"true_SN_{self.name}", self.ndata):
mag_true = sample( mag_true = sample(
f"mag_true_{self.name}", Normal(mag_mean, mag_std)) f"mag_true_{self.name}", Normal(mag_mean, mag_std))
x1_true = sample( x1_true = sample(
@ -726,7 +751,7 @@ class PV_LogLikelihood(BaseFlowValidationModel):
f"c_true_{self.name}", Normal(c_mean, c_std)) f"c_true_{self.name}", Normal(c_mean, c_std))
# Log-likelihood of the observed magnitudes. # Log-likelihood of the observed magnitudes.
if self.maxmag_selection is None: if self.mag_selection_kind is None:
ll0 += jnp.sum(normal_logpdf( ll0 += jnp.sum(normal_logpdf(
mag_true, self.mag, self.e_mag)) mag_true, self.mag, self.e_mag))
else: else:
@ -740,9 +765,12 @@ class PV_LogLikelihood(BaseFlowValidationModel):
mag_true = self.mag mag_true = self.mag
x1_true = self.x1 x1_true = self.x1
c_true = self.c c_true = self.c
if inference_method == "mike":
e2_mu = e2_distmod_SN( e2_mu = e2_distmod_SN(
self.e2_mag, self.e2_x1, self.e2_c, alpha_cal, beta_cal, self.e2_mag, self.e2_x1, self.e2_c, alpha_cal,
e_mu) beta_cal, e_mu)
else:
e2_mu = jnp.ones_like(mag_true) * e_mu**2
mu = distmod_SN( mu = distmod_SN(
mag_true, x1_true, c_true, mag_cal, alpha_cal, beta_cal) mag_true, x1_true, c_true, mag_cal, alpha_cal, beta_cal)
@ -761,22 +789,25 @@ class PV_LogLikelihood(BaseFlowValidationModel):
"mag", self.name, self.mag_min, self.mag_max) "mag", self.name, self.mag_min, self.mag_max)
eta_mean, eta_std = sample_gaussian_hyperprior( eta_mean, eta_std = sample_gaussian_hyperprior(
"eta", self.name, self.eta_min, self.eta_max) "eta", self.name, self.eta_min, self.eta_max)
corr_mag_eta = sample("corr_mag_eta", Uniform(-1, 1)) corr_mag_eta = sample(
f"corr_mag_eta_{self.name}", Uniform(-1, 1))
loc = jnp.array([mag_mean, eta_mean]) loc = jnp.array([mag_mean, eta_mean])
cov = jnp.array( cov = jnp.array(
[[mag_std**2, corr_mag_eta * mag_std * eta_std], [[mag_std**2, corr_mag_eta * mag_std * eta_std],
[corr_mag_eta * mag_std * eta_std, eta_std**2]]) [corr_mag_eta * mag_std * eta_std, eta_std**2]])
with plate("true_TFR", self.ndata): with plate(f"true_TFR_{self.name}", self.ndata):
x_true = sample("x_TFR", MultivariateNormal(loc, cov)) x_true = sample(
f"x_TFR_{self.name}", MultivariateNormal(loc, cov))
mag_true, eta_true = x_true[..., 0], x_true[..., 1] mag_true, eta_true = x_true[..., 0], x_true[..., 1]
# Log-likelihood of the observed magnitudes. # Log-likelihood of the observed magnitudes.
if self.maxmag_selection is not None: if self.mag_selection_kind == "hard":
ll0 += jnp.sum(upper_truncated_normal_logpdf( ll0 += jnp.sum(upper_truncated_normal_logpdf(
self.mag, mag_true, self.e_mag, self.maxmag_selection)) self.mag, mag_true, self.e_mag,
elif self.toy_selection is not None: self.mag_selection_max))
elif self.mag_selection_kind == "soft":
ll_mag = self.log_Fm ll_mag = self.log_Fm
ll_mag += normal_logpdf(self.mag, mag_true, self.e_mag) ll_mag += normal_logpdf(self.mag, mag_true, self.e_mag)
@ -805,8 +836,11 @@ class PV_LogLikelihood(BaseFlowValidationModel):
else: else:
eta_true = self.eta eta_true = self.eta
mag_true = self.mag mag_true = self.mag
if inference_method == "mike":
e2_mu = e2_distmod_TFR( e2_mu = e2_distmod_TFR(
self.e2_mag, self.e2_eta, eta_true, b, c, e_mu) self.e2_mag, self.e2_eta, eta_true, b, c, e_mu)
else:
e2_mu = jnp.ones_like(mag_true) * e_mu**2
mu = distmod_TFR(mag_true, eta_true, a, b, c) mu = distmod_TFR(mag_true, eta_true, a, b, c)
elif self.kind == "simple": elif self.kind == "simple":
@ -821,7 +855,10 @@ class PV_LogLikelihood(BaseFlowValidationModel):
raise NotImplementedError("Bayes for simple not implemented.") raise NotImplementedError("Bayes for simple not implemented.")
else: else:
mu_true = self.mu mu_true = self.mu
if inference_method == "mike":
e2_mu = e_mu**2 + self.e2_mu e2_mu = e_mu**2 + self.e2_mu
else:
e2_mu = jnp.ones_like(mag_true) * e_mu**2
mu = mu_true + dmu mu = mu_true + dmu
else: else:
@ -895,8 +932,7 @@ def PV_validation_model(models, distmod_hyperparams_per_model,
############################################################################### ###############################################################################
def get_model(loader, zcmb_min=0.0, zcmb_max=None, maxmag_selection=None, def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None):
toy_selection=None):
""" """
Get a model and extract the relevant data from the loader. Get a model and extract the relevant data from the loader.
@ -908,25 +944,20 @@ def get_model(loader, zcmb_min=0.0, zcmb_max=None, maxmag_selection=None,
Minimum observed redshift in the CMB frame to include. Minimum observed redshift in the CMB frame to include.
zcmb_max : float, optional zcmb_max : float, optional
Maximum observed redshift in the CMB frame to include. Maximum observed redshift in the CMB frame to include.
maxmag_selection : float, optional mag_selection : dict, optional
Maximum magnitude selection threshold. Magnitude selection parameters.
toy_selection : tuple of length 3, optional
Toy magnitude selection paramers `m1`, `m2` and `a` for TFRs of the
Boubel+24 model.
Returns Returns
------- -------
model : NumPyro model model : NumPyro model
""" """
zcmb_min = 0.0 if zcmb_min is None else zcmb_min
zcmb_max = np.infty if zcmb_max is None else zcmb_max zcmb_max = np.infty if zcmb_max is None else zcmb_max
los_overdensity = loader.los_density los_overdensity = loader.los_density
los_velocity = loader.los_radial_velocity los_velocity = loader.los_radial_velocity
kind = loader._catname kind = loader._catname
if maxmag_selection is not None and kind != "2MTF":
raise ValueError("Threshold magnitude selection implemented only for 2MTF.") # noqa
if kind in ["LOSS", "Foundation"]: if kind in ["LOSS", "Foundation"]:
keys = ["RA", "DEC", "z_CMB", "mB", "x1", "c", "e_mB", "e_x1", "e_c"] keys = ["RA", "DEC", "z_CMB", "mB", "x1", "c", "e_mB", "e_x1", "e_c"]
RA, dec, zCMB, mag, x1, c, e_mag, e_x1, e_c = ( RA, dec, zCMB, mag, x1, c, e_mag, e_x1, e_c = (
@ -941,7 +972,7 @@ def get_model(loader, zcmb_min=0.0, zcmb_max=None, maxmag_selection=None,
model = PV_LogLikelihood( model = PV_LogLikelihood(
los_overdensity[:, mask], los_velocity[:, mask], los_overdensity[:, mask], los_velocity[:, mask],
RA[mask], dec[mask], zCMB[mask], e_zCMB, calibration_params, RA[mask], dec[mask], zCMB[mask], e_zCMB, calibration_params,
maxmag_selection, loader.rdist, loader._Omega_m, "SN", name=kind) mag_selection, loader.rdist, loader._Omega_m, "SN", name=kind)
elif "Pantheon+" in kind: elif "Pantheon+" in kind:
keys = ["RA", "DEC", "zCMB", "mB", "x1", "c", "biasCor_m_b", "mBERR", 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",
@ -969,25 +1000,18 @@ def get_model(loader, zcmb_min=0.0, zcmb_max=None, maxmag_selection=None,
model = PV_LogLikelihood( model = PV_LogLikelihood(
los_overdensity[:, mask], los_velocity[:, mask], los_overdensity[:, mask], los_velocity[:, mask],
RA[mask], dec[mask], zCMB[mask], e_zCMB[mask], calibration_params, RA[mask], dec[mask], zCMB[mask], e_zCMB[mask], calibration_params,
maxmag_selection, loader.rdist, loader._Omega_m, "SN", name=kind) mag_selection, loader.rdist, loader._Omega_m, "SN", name=kind)
elif kind in ["SFI_gals", "2MTF", "SFI_gals_masked"]: elif kind in ["SFI_gals", "2MTF", "SFI_gals_masked"]:
keys = ["RA", "DEC", "z_CMB", "mag", "eta", "e_mag", "e_eta"] 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) RA, dec, zCMB, mag, eta, e_mag, e_eta = (loader.cat[k] for k in keys)
if kind == "SFI_gals" and toy_selection is not None:
if len(toy_selection) != 3:
raise ValueError("Toy selection must be a tuple with 3 elements.") # noqa
m1, m2, a = toy_selection
fprint(f"using toy selection with m1 = {m1}, m2 = {m2}, a = {a}.")
mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) mask = (zCMB < zcmb_max) & (zCMB > zcmb_min)
calibration_params = {"mag": mag[mask], "eta": eta[mask], calibration_params = {"mag": mag[mask], "eta": eta[mask],
"e_mag": e_mag[mask], "e_eta": e_eta[mask]} "e_mag": e_mag[mask], "e_eta": e_eta[mask]}
model = PV_LogLikelihood( model = PV_LogLikelihood(
los_overdensity[:, mask], los_velocity[:, mask], los_overdensity[:, mask], los_velocity[:, mask],
RA[mask], dec[mask], zCMB[mask], None, calibration_params, RA[mask], dec[mask], zCMB[mask], None, calibration_params,
maxmag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind, mag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind)
toy_selection=toy_selection)
elif "CF4_TFR_" in kind: elif "CF4_TFR_" in kind:
# The full name can be e.g. "CF4_TFR_not2MTForSFI_i" or "CF4_TFR_i". # The full name can be e.g. "CF4_TFR_not2MTForSFI_i" or "CF4_TFR_i".
band = kind.split("_")[-1] band = kind.split("_")[-1]
@ -1001,7 +1025,7 @@ def get_model(loader, zcmb_min=0.0, zcmb_max=None, maxmag_selection=None,
not_matched_to_2MTF_or_SFI = not_matched_to_2MTF_or_SFI.astype(bool) not_matched_to_2MTF_or_SFI = not_matched_to_2MTF_or_SFI.astype(bool)
# NOTE: fiducial uncertainty until we can get the actual values. # NOTE: fiducial uncertainty until we can get the actual values.
e_mag = 0.001 * np.ones_like(mag) e_mag = 0.05 * np.ones_like(mag)
z_obs /= SPEED_OF_LIGHT z_obs /= SPEED_OF_LIGHT
eta -= 2.5 eta -= 2.5
@ -1026,7 +1050,7 @@ def get_model(loader, zcmb_min=0.0, zcmb_max=None, maxmag_selection=None,
model = PV_LogLikelihood( model = PV_LogLikelihood(
los_overdensity[:, mask], los_velocity[:, mask], los_overdensity[:, mask], los_velocity[:, mask],
RA[mask], dec[mask], z_obs[mask], None, calibration_params, RA[mask], dec[mask], z_obs[mask], None, calibration_params,
maxmag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind) mag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind)
elif kind in ["CF4_GroupAll"]: elif kind in ["CF4_GroupAll"]:
# Note, this for some reason works terribly. # Note, this for some reason works terribly.
keys = ["RA", "DE", "Vcmb", "DMzp", "eDM"] keys = ["RA", "DE", "Vcmb", "DMzp", "eDM"]
@ -1042,7 +1066,7 @@ def get_model(loader, zcmb_min=0.0, zcmb_max=None, maxmag_selection=None,
model = PV_LogLikelihood( model = PV_LogLikelihood(
los_overdensity[:, mask], los_velocity[:, mask], los_overdensity[:, mask], los_velocity[:, mask],
RA[mask], dec[mask], zCMB[mask], None, calibration_params, RA[mask], dec[mask], zCMB[mask], None, calibration_params,
maxmag_selection, loader.rdist, loader._Omega_m, "simple", mag_selection, loader.rdist, loader._Omega_m, "simple",
name=kind) name=kind)
else: else:
raise ValueError(f"Catalogue `{kind}` not recognized.") raise ValueError(f"Catalogue `{kind}` not recognized.")

64
csiborgtools/flow/selection.py
View file

@ -13,56 +13,64 @@
# with this program; if not, write to the Free Software Foundation, Inc., # with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
"""Selection functions for peculiar velocities.""" """Selection functions for peculiar velocities."""
import numpy as np
from jax import numpy as jnp from jax import numpy as jnp
from scipy.integrate import quad from numpyro import factor, sample
from scipy.optimize import minimize from numpyro.distributions import Normal, Uniform
from quadax import simpson
class ToyMagnitudeSelection: class ToyMagnitudeSelection:
""" """
Toy magnitude selection according to Boubel et al 2024. Toy magnitude selection according to Boubel+2024 [1].
References
----------
[1] https://www.arxiv.org/abs/2408.03660
""" """
def __init__(self): def __init__(self):
pass self.mrange = jnp.linspace(0, 25, 1000)
def log_true_pdf(self, m, m1): def log_true_pdf(self, m, alpha, m1):
"""Unnormalized `true' PDF.""" """Unnormalized `true' PDF."""
return 0.6 * (m - m1) return alpha * (m - m1)
def log_selection_function(self, m, m1, m2, a): def log_selection_function(self, m, m1, m2, a):
return np.where(m <= m1, """Logarithm of the Boubel+2024 selection function."""
return jnp.where(m <= m1,
0, 0,
a * (m - m2)**2 - a * (m1 - m2)**2 - 0.6 * (m - m1)) a * (m - m2)**2 - a * (m1 - m2)**2 - 0.6 * (m - m1))
def log_observed_pdf(self, m, m1, m2, a): def log_observed_pdf(self, m, alpha, m1, m2, a):
# Calculate the normalization constant """
f = lambda m: 10**(self.log_true_pdf(m, m1) # noqa Logarithm of the unnormalized observed PDF, which is the product
+ self.log_selection_function(m, m1, m2, a)) of the true PDF and the selection function.
mmin, mmax = 0, 25 """
norm = quad(f, mmin, mmax)[0] y = 10**(self.log_true_pdf(self.mrange, alpha, m1)
+ self.log_selection_function(self.mrange, m1, m2, a)
)
norm = simpson(y, x=self.mrange)
return (self.log_true_pdf(m, m1) return (self.log_true_pdf(m, alpha, m1)
+ self.log_selection_function(m, m1, m2, a) + self.log_selection_function(m, m1, m2, a)
- np.log10(norm)) - jnp.log10(norm))
def fit(self, mag): def __call__(self, mag):
"""NumPyro model, uses an informative prior on `alpha`."""
alpha = sample("alpha", Normal(0.6, 0.1))
m1 = sample("m1", Uniform(0, 25))
m2 = sample("m2", Uniform(0, 25))
a = sample("a", Uniform(-10, 0))
def loss(x): ll = jnp.sum(self.log_observed_pdf(mag, alpha, m1, m2, a))
m1, m2, a = x factor("ll", ll)
if a >= 0:
return np.inf
return -np.sum(self.log_observed_pdf(mag, m1, m2, a))
x0 = [12.0, 12.5, -0.1]
return minimize(loss, x0, method="Nelder-Mead")
def toy_log_magnitude_selection(mag, m1, m2, a): def toy_log_magnitude_selection(mag, m1, m2, a):
"""JAX implementation of `ToyMagnitudeSelection` but natural logarithm.""" """
JAX implementation of `ToyMagnitudeSelection` but natural logarithm,
whereas the one in `ToyMagnitudeSelection` is base 10.
"""
return jnp.log(10) * jnp.where( return jnp.log(10) * jnp.where(
mag <= m1, mag <= m1,
0, 0,

3
csiborgtools/params.py
View file

@ -103,6 +103,9 @@ def simname2Omega_m(simname):
"CF4": 0.3, "CF4": 0.3,
"CF4gp": 0.3, "CF4gp": 0.3,
"Lilow2024": 0.3175, "Lilow2024": 0.3175,
"IndranilVoid_exp": 0.3,
"IndranilVoid_gauss": 0.3,
"no_field": 0.3,
} }
omega_m = d.get(simname, None) omega_m = d.get(simname, None)

61
csiborgtools/read/paths.py
View file

@ -15,11 +15,12 @@
""" """
CSiBORG paths manager. CSiBORG paths manager.
""" """
import datetime
from glob import glob from glob import glob
from os import makedirs, listdir from os import listdir, makedirs
from os.path import isdir, join from os.path import exists, getmtime, isdir, join
from warnings import warn
from re import search from re import search
from warnings import warn
import numpy import numpy
@ -117,15 +118,15 @@ class Paths:
files = glob(join(self.quijote_dir, "fiducial_processed", files = glob(join(self.quijote_dir, "fiducial_processed",
"chain_*")) "chain_*"))
files = [int(search(r'chain_(\d+)', f).group(1)) for f in files] files = [int(search(r'chain_(\d+)', f).group(1)) for f in files]
elif simname == "Carrick2015":
return [0]
elif simname == "CF4": elif simname == "CF4":
files = glob(join(self.CF4_dir, "CF4_new_128-z008_realization*_delta.fits")) # noqa files = glob(join(self.CF4_dir, "CF4_new_128-z008_realization*_delta.fits")) # noqa
files = [search(r'realization(\d+)_delta\.fits', file).group(1) files = [search(r'realization(\d+)_delta\.fits', file).group(1)
for file in files if search(r'realization(\d+)_delta\.fits', file)] # noqa for file in files if search(r'realization(\d+)_delta\.fits', file)] # noqa
files = [int(file) for file in files] files = [int(file) for file in files]
elif simname == "Lilow2024": # Downsample to only 20 realisations
return [0] files = files[::5]
elif simname in ["Carrick2015", "Lilow2024", "no_field"] or "IndranilVoid" in simname: # noqa
files = [0]
else: else:
raise ValueError(f"Unknown simulation name `{simname}`.") raise ValueError(f"Unknown simulation name `{simname}`.")
@ -635,6 +636,50 @@ class Paths:
try_create_directory(fdir) try_create_directory(fdir)
return join(fdir, f"los_{catalogue}_{simnname}.hdf5") return join(fdir, f"los_{catalogue}_{simnname}.hdf5")
def flow_validation(self, fdir, simname, catalogue, inference_method,
smooth=None, nsim=None, zcmb_min=None, zcmb_max=None,
mag_selection=None, sample_alpha=False,
sample_beta=False, sample_Vmono=False,
sample_mag_dipole=False, sample_curvature=False):
"""Flow validation file path."""
if isinstance(catalogue, list) and len(catalogue) == 1:
catalogue = catalogue[0]
if isinstance(catalogue, list):
catalogue = "_".join(catalogue)
if smooth == 0:
smooth = None
fname = f"samples_{simname}_{catalogue}_{inference_method}_"
keys = ["smooth", "nsim", "zcmb_min", "zcmb_max", "mag_selection",
"sample_alpha", "sample_beta", "sample_Vmono",
"sample_mag_dipole", "sample_curvature"]
values = [smooth, nsim, zcmb_min, zcmb_max, mag_selection,
sample_alpha, sample_beta, sample_Vmono, sample_mag_dipole,
sample_curvature]
for key, value in zip(keys, values):
if isinstance(value, bool):
if value:
fname += f"{key}_"
elif value is not None:
fname += f"{key}_{value}_"
fname = fname.strip("_")
fname = join(fdir, f"{fname}.hdf5")
# Print the last modified time of the file if it exists.
if exists(fname):
mtime = getmtime(fname)
mtime = datetime.datetime.fromtimestamp(mtime)
mtime = mtime.strftime("%d/%m/%Y %H:%M:%S")
print(f"File: {fname}")
print(f"Last modified: {mtime}")
return fname
def field_projected(self, simname, kind): def field_projected(self, simname, kind):
""" """
Path to the files containing the projected fields on the sky. Path to the files containing the projected fields on the sky.
@ -653,5 +698,3 @@ class Paths:
fdir = join(self.postdir, "field_projected") fdir = join(self.postdir, "field_projected")
try_create_directory(fdir) try_create_directory(fdir)
return join(fdir, f"{simname}_{kind}_volume_weighted.hdf5") return join(fdir, f"{simname}_{kind}_volume_weighted.hdf5")

588
notebooks/flow/PV_data.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

915
notebooks/flow/reconstruction_comparison.ipynb
View file

File diff suppressed because one or more lines are too long

248
notebooks/flow/reconstruction_comparison.py
View file

@ -13,8 +13,7 @@
# with this program; if not, write to the Free Software Foundation, Inc., # with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
"""Script to help with plots in `flow_calibration.ipynb`.""" """Script to help with plots in `flow_calibration.ipynb`."""
from copy import copy from copy import copy, deepcopy
from os.path import join
import numpy as np import numpy as np
from jax import numpy as jnp from jax import numpy as jnp
@ -41,25 +40,6 @@ def cartesian_to_radec(x, y, z):
return d, ra, dec return d, ra, dec
###############################################################################
# Get the filename of the samples #
###############################################################################
def get_fname(simname, catalogue, ksmooth=0, nsim=None, sample_beta=True):
"""Get the filename of the HDF5 file containing the posterior samples."""
FDIR = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/peculiar_velocity/" # noqa
fname = join(FDIR, f"samples_{simname}_{catalogue}_ksmooth{ksmooth}.hdf5")
if nsim is not None:
fname = fname.replace(".hdf5", f"_nsim{nsim}.hdf5")
if sample_beta:
fname = fname.replace(".hdf5", "_sample_beta.hdf5")
return fname
############################################################################### ###############################################################################
# Convert names to LaTeX # # Convert names to LaTeX #
############################################################################### ###############################################################################
@ -69,30 +49,96 @@ def names_to_latex(names, for_corner=False):
"""Convert the names of the parameters to LaTeX.""" """Convert the names of the parameters to LaTeX."""
ltx = {"alpha": "\\alpha", ltx = {"alpha": "\\alpha",
"beta": "\\beta", "beta": "\\beta",
"Vmag": "V_{\\rm ext}", "Vmag": "V_{\\rm ext} ~ [\\mathrm{km} / \\mathrm{s}]",
"sigma_v": "\\sigma_v", "Vx": "V_x ~ [\\mathrm{km} / \\mathrm{s}]",
"Vy": "V_y ~ [\\mathrm{km} / \\mathrm{s}]",
"Vz": "V_z ~ [\\mathrm{km} / \\mathrm{s}]",
"sigma_v": "\\sigma_v ~ [\\mathrm{km} / \\mathrm{s}]",
"alpha_cal": "\\mathcal{A}", "alpha_cal": "\\mathcal{A}",
"beta_cal": "\\mathcal{B}", "beta_cal": "\\mathcal{B}",
"mag_cal": "\\mathcal{M}", "mag_cal": "\\mathcal{M}",
"e_mu": "\\sigma_\\mu", "l": "\\ell ~ [\\mathrm{deg}]",
"aTF": "a_{\\rm TF}", "b": "b ~ [\\mathrm{deg}]",
"bTF": "b_{\\rm TF}",
} }
ltx_corner = {"alpha": r"$\alpha$", ltx_corner = {"alpha": r"$\alpha$",
"beta": r"$\beta$", "beta": r"$\beta$",
"Vmag": r"$V_{\rm ext}$", "Vmag": r"$V_{\rm ext}$",
"l": r"$\ell_{V_{\rm ext}}$", "l": r"$\ell$",
"b": r"$b_{V_{\rm ext}}$", "b": r"$b$",
"sigma_v": r"$\sigma_v$", "sigma_v": r"$\sigma_v$",
"alpha_cal": r"$\mathcal{A}$", "alpha_cal": r"$\mathcal{A}$",
"beta_cal": r"$\mathcal{B}$", "beta_cal": r"$\mathcal{B}$",
"mag_cal": r"$\mathcal{M}$", "mag_cal": r"$\mathcal{M}$",
"e_mu": r"$\sigma_\mu$",
"aTF": r"$a_{\rm TF}$",
"bTF": r"$b_{\rm TF}$",
} }
names = copy(names)
for i, name in enumerate(names):
if "SFI_gals" in name:
names[i] = names[i].replace("SFI_gals", "SFI")
if "CF4_GroupAll" in name:
names[i] = names[i].replace("CF4_GroupAll", "CF4Group")
if "CF4_TFR_i" in name:
names[i] = names[i].replace("CF4_TFR_i", "CF4,TFR")
for cat in ["2MTF", "SFI", "CF4,TFR"]:
ltx[f"a_{cat}"] = f"a_{{\\rm TF}}^{{\\rm {cat}}}"
ltx[f"b_{cat}"] = f"b_{{\\rm TF}}^{{\\rm {cat}}}"
ltx[f"c_{cat}"] = f"c_{{\\rm TF}}^{{\\rm {cat}}}"
ltx[f"corr_mag_eta_{cat}"] = f"\\rho_{{m,\\eta}}^{{\\rm {cat}}}"
ltx[f"eta_mean_{cat}"] = f"\\widehat{{\\eta}}^{{\\rm {cat}}}"
ltx[f"eta_std_{cat}"] = f"\\widehat{{\\sigma}}_\\eta^{{\\rm {cat}}}"
ltx[f"mag_mean_{cat}"] = f"\\widehat{{m}}^{{\\rm {cat}}}"
ltx[f"mag_std_{cat}"] = f"\\widehat{{\\sigma}}_m^{{\\rm {cat}}}"
ltx_corner[f"a_{cat}"] = rf"$a_{{\rm TF}}^{{\rm {cat}}}$"
ltx_corner[f"b_{cat}"] = rf"$b_{{\rm TF}}^{{\rm {cat}}}$"
ltx_corner[f"c_{cat}"] = rf"$c_{{\rm TF}}^{{\rm {cat}}}$"
ltx_corner[f"corr_mag_eta_{cat}"] = rf"$\rho_{{m,\eta}}^{{\rm {cat}}}$"
ltx_corner[f"eta_mean_{cat}"] = rf"$\widehat{{\eta}}^{{\rm {cat}}}$"
ltx_corner[f"eta_std_{cat}"] = rf"$\widehat{{\sigma}}_\eta^{{\rm {cat}}}$" # noqa
ltx_corner[f"mag_mean_{cat}"] = rf"$\widehat{{m}}^{{\rm {cat}}}$"
ltx_corner[f"mag_std_{cat}"] = rf"$\widehat{{\sigma}}_m^{{\rm {cat}}}$"
for cat in ["2MTF", "SFI", "Foundation", "LOSS", "CF4Group", "CF4,TFR"]:
ltx[f"alpha_{cat}"] = f"\\alpha^{{\\rm {cat}}}"
ltx[f"e_mu_{cat}"] = f"\\sigma_{{\\mu}}^{{\\rm {cat}}}"
ltx[f"a_dipole_mag_{cat}"] = f"\\epsilon_{{\\rm mag}}^{{\\rm {cat}}}"
ltx[f"a_dipole_l_{cat}"] = f"\\epsilon_{{\\ell}}^{{\\rm {cat}}} ~ [\\mathrm{{deg}}]" # noqa
ltx[f"a_dipole_b_{cat}"] = f"\\epsilon_{{b}}^{{\\rm {cat}}} ~ [\\mathrm{{deg}}]" # noqa
ltx["a_dipole_mag"] = "\\epsilon_{{\\rm mag}}"
ltx["a_dipole_l"] = "\\epsilon_{{\\ell}} ~ [\\mathrm{{deg}}]"
ltx["a_dipole_b"] = "\\epsilon_{{b}} ~ [\\mathrm{{deg}}]"
ltx_corner[f"alpha_{cat}"] = rf"$\alpha^{{\rm {cat}}}$"
ltx_corner[f"e_mu_{cat}"] = rf"$\sigma_{{\mu}}^{{\rm {cat}}}$"
ltx_corner[f"a_dipole_mag_{cat}"] = rf"$\epsilon_{{\rm mag}}^{{\rm {cat}}}$" # noqa
ltx_corner[f"a_dipole_l_{cat}"] = rf"$\epsilon_{{\ell}}^{{\rm {cat}}}$"
ltx_corner[f"a_dipole_b_{cat}"] = rf"$\epsilon_{{b}}^{{\rm {cat}}}$"
for cat in ["Foundation", "LOSS"]:
ltx[f"alpha_cal_{cat}"] = f"\\mathcal{{A}}^{{\\rm {cat}}}"
ltx[f"beta_cal_{cat}"] = f"\\mathcal{{B}}^{{\\rm {cat}}}"
ltx[f"mag_cal_{cat}"] = f"\\mathcal{{M}}^{{\\rm {cat}}}"
ltx_corner[f"alpha_cal_{cat}"] = rf"$\mathcal{{A}}^{{\rm {cat}}}$"
ltx_corner[f"beta_cal_{cat}"] = rf"$\mathcal{{B}}^{{\rm {cat}}}$"
ltx_corner[f"mag_cal_{cat}"] = rf"$\mathcal{{M}}^{{\rm {cat}}}$"
for cat in ["CF4Group"]:
ltx[f"dmu_{cat}"] = f"\\Delta\\mu^{{\\rm {cat}}}"
ltx[f"dmu_dipole_mag_{cat}"] = f"\\epsilon_\\mu_{{\\rm mag}}^{{\\rm {cat}}}" # noqa
ltx[f"dmu_dipole_l_{cat}"] = f"\\epsilon_\\mu_{{\\ell}}^{{\\rm {cat}}} ~ [\\mathrm{{deg}}]" # noqa
ltx[f"dmu_dipole_b_{cat}"] = f"\\epsilon_\\mu_{{b}}^{{\\rm {cat}}} ~ [\\mathrm{{deg}}]" # noqa
ltx_corner[f"dmu_{cat}"] = rf"$\Delta\mu_{{0}}^{{\rm {cat}}}$"
ltx_corner[f"dmu_dipole_mag_{cat}"] = rf"$\epsilon_{{\rm mag}}^{{\rm {cat}}}$" # noqa
ltx_corner[f"dmu_dipole_l_{cat}"] = rf"$\epsilon_{{\ell}}^{{\rm {cat}}}$" # noqa
ltx_corner[f"dmu_dipole_b_{cat}"] = rf"$\epsilon_{{b}}^{{\rm {cat}}}$" # noqa
labels = copy(names) labels = copy(names)
for i, label in enumerate(names): for i, label in enumerate(names):
if for_corner: if for_corner:
@ -113,21 +159,35 @@ def simname_to_pretty(simname):
} }
if isinstance(simname, list): if isinstance(simname, list):
return [ltx[s] if s in ltx else s for s in simname] names = [ltx[s] if s in ltx else s for s in simname]
return "".join([f"{n}, " for n in names]).rstrip(", ")
return ltx[simname] if simname in ltx else simname return ltx[simname] if simname in ltx else simname
def catalogue_to_pretty(catalogue):
ltx = {"SFI_gals": "SFI",
"CF4_TFR_not2MTForSFI_i": r"CF4 $i$-band",
"CF4_TFR_i": r"CF4 TFR $i$",
"CF4_TFR_w1": r"CF4 TFR W1",
}
if isinstance(catalogue, list):
names = [ltx[s] if s in ltx else s for s in catalogue]
return "".join([f"{n}, " for n in names]).rstrip(", ")
return ltx[catalogue] if catalogue in ltx else catalogue
############################################################################### ###############################################################################
# Read in goodness-of-fit # # Read in goodness-of-fit #
############################################################################### ###############################################################################
def get_gof(kind, simname, catalogue, ksmooth=0, nsim=None, sample_beta=True): def get_gof(kind, fname):
"""Read in the goodness-of-fit statistics `kind`.""" """Read in the goodness-of-fit statistics `kind`."""
if kind not in ["BIC", "AIC", "lnZ"]: if kind not in ["BIC", "AIC", "neg_lnZ_harmonic"]:
raise ValueError("`kind` must be one of 'BIC', 'AIC', 'lnZ'") raise ValueError("`kind` must be one of 'BIC', 'AIC', 'neg_lnZ_harmonic'.") # noqa
fname = get_fname(simname, catalogue, ksmooth, nsim, sample_beta)
with File(fname, 'r') as f: with File(fname, 'r') as f:
return f[f"gof/{kind}"][()] return f[f"gof/{kind}"][()]
@ -136,29 +196,48 @@ def get_gof(kind, simname, catalogue, ksmooth=0, nsim=None, sample_beta=True):
# Read in samples # # Read in samples #
############################################################################### ###############################################################################
def get_samples(simname, catalogue, ksmooth=0, nsim=None, sample_beta=True, def get_samples(fname, convert_Vext_to_galactic=True):
convert_Vext_to_galactic=True):
"""Read in the samples from the HDF5 file.""" """Read in the samples from the HDF5 file."""
fname = get_fname(simname, catalogue, ksmooth, nsim, sample_beta)
samples = {} samples = {}
with File(fname, 'r') as f: with File(fname, 'r') as f:
grp = f["samples"] grp = f["samples"]
for key in grp.keys(): for key in grp.keys():
samples[key] = grp[key][...] samples[key] = grp[key][...]
# Rename TF parameters
if "a" in samples:
samples["aTF"] = samples.pop("a")
if "b" in samples:
samples["bTF"] = samples.pop("b")
if convert_Vext_to_galactic: if convert_Vext_to_galactic:
Vext = samples.pop("Vext") Vext = samples.pop("Vext")
samples["Vmag"] = np.linalg.norm(Vext, axis=1) samples["Vmag"] = np.linalg.norm(Vext, axis=1)
Vext = csiborgtools.cartesian_to_radec(Vext) Vext = csiborgtools.cartesian_to_radec(Vext)
samples["l"], samples["b"] = csiborgtools.radec_to_galactic( samples["l"], samples["b"] = csiborgtools.radec_to_galactic(
Vext[:, 1], Vext[:, 2]) Vext[:, 1], Vext[:, 2])
else:
Vext = samples.pop("Vext")
samples["Vx"] = Vext[:, 0]
samples["Vy"] = Vext[:, 1]
samples["Vz"] = Vext[:, 2]
keys = list(samples.keys())
for key in keys:
if "dmu_dipole_" in key:
dmu = samples.pop(key)
dmu = csiborgtools.cartesian_to_radec(dmu)
dmu_mag = dmu[:, 0]
l, b = csiborgtools.radec_to_galactic(dmu[:, 1], dmu[:, 2])
samples[key.replace("dmu_dipole_", "dmu_dipole_mag_")] = dmu_mag
samples[key.replace("dmu_dipole_", "dmu_dipole_l_")] = l
samples[key.replace("dmu_dipole_", "dmu_dipole_b_")] = b
if "a_dipole" in key:
adipole = samples.pop(key)
adipole = csiborgtools.cartesian_to_radec(adipole)
adipole_mag = adipole[:, 0]
l, b = csiborgtools.radec_to_galactic(adipole[:, 1], adipole[:, 2])
samples[key.replace("a_dipole", "a_dipole_mag")] = adipole_mag
samples[key.replace("a_dipole", "a_dipole_l")] = l
samples[key.replace("a_dipole", "a_dipole_b")] = b
return samples return samples
@ -180,12 +259,20 @@ def get_bulkflow_simulation(simname, convert_to_galactic=True):
return r, B return r, B
def get_bulkflow(simname, catalogue, ksmooth=0, nsim=None, sample_beta=True, def get_bulkflow(fname, simname, convert_to_galactic=True, downsample=1,
convert_to_galactic=True, weight_simulations=True, Rmax=125):
downsample=1, Rmax=125): # Read in the samples
with File(fname, "r") as f:
Vext = f["samples/Vext"][...]
try:
beta = f["samples/beta"][...]
except KeyError:
beta = jnp.ones(len(Vext))
# Read in the bulk flow # Read in the bulk flow
f = np.load(f"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_{simname}.npz") # noqa f = np.load(f"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_{simname}.npz") # noqa
r = f["distances"] r = f["distances"]
# Shape of B_i is (nsims, nradial) # Shape of B_i is (nsims, nradial)
Bx, By, Bz = (f["cumulative_velocity"][..., i] for i in range(3)) Bx, By, Bz = (f["cumulative_velocity"][..., i] for i in range(3))
@ -197,38 +284,18 @@ def get_bulkflow(simname, catalogue, ksmooth=0, nsim=None, sample_beta=True,
By = By[:, mask] By = By[:, mask]
Bz = Bz[:, mask] Bz = Bz[:, mask]
# Read in the samples Vext = Vext[::downsample]
fname_samples = get_fname(simname, catalogue, ksmooth, nsim, sample_beta) beta = beta[::downsample]
with File(fname_samples, 'r') as f:
# Shape of Vext_i is (nsamples,)
Vext_x, Vext_y, Vext_z = (f["samples/Vext"][...][::downsample, i] for i in range(3)) # noqa
nsamples = len(Vext_x)
if weight_simulations: # Multiply the simulation velocities by beta.
simulation_weights = jnp.exp(f["simulation_weights"][...])[::downsample] # noqa
else:
nsims = len(Bx)
simulation_weights = jnp.ones((nsamples, nsims)) / nsims
if sample_beta:
beta = f["samples/beta"][...][::downsample]
else:
beta = jnp.ones(nsamples)
# Multiply the simulation velocities by beta
Bx = Bx[..., None] * beta Bx = Bx[..., None] * beta
By = By[..., None] * beta By = By[..., None] * beta
Bz = Bz[..., None] * beta Bz = Bz[..., None] * beta
# Shape of B_i is (nsims, nradial, nsamples) # Add V_ext, shape of B_i is `(nsims, nradial, nsamples)``
Bx = Bx + Vext_x Bx = Bx + Vext[:, 0]
By = By + Vext_y By = By + Vext[:, 1]
Bz = Bz + Vext_z Bz = Bz + Vext[:, 2]
simulation_weights = simulation_weights.T[:, None, :]
Bx = jnp.sum(Bx * simulation_weights, axis=0)
By = jnp.sum(By * simulation_weights, axis=0)
Bz = jnp.sum(Bz * simulation_weights, axis=0)
if convert_to_galactic: if convert_to_galactic:
Bmag, Bl, Bb = cartesian_to_radec(Bx, By, Bz) Bmag, Bl, Bb = cartesian_to_radec(Bx, By, Bz)
@ -237,6 +304,8 @@ def get_bulkflow(simname, catalogue, ksmooth=0, nsim=None, sample_beta=True,
else: else:
B = np.stack([Bx, By, Bz], axis=-1) B = np.stack([Bx, By, Bz], axis=-1)
# Stack over the simulations
B = np.hstack([B[i] for i in range(len(B))])
return r, B return r, B
############################################################################### ###############################################################################
@ -245,25 +314,30 @@ def get_bulkflow(simname, catalogue, ksmooth=0, nsim=None, sample_beta=True,
def samples_for_corner(samples): def samples_for_corner(samples):
samples = deepcopy(samples)
# Remove the true parameters of each galaxy.
keys = list(samples.keys())
for key in keys:
# Generally don't want to plot the true latent parameters..
if "x_TFR" in key or "_true_" in key:
samples.pop(key)
keys = list(samples.keys())
if any(x.ndim > 1 for x in samples.values()): if any(x.ndim > 1 for x in samples.values()):
raise ValueError("All samples must be 1D arrays.") raise ValueError("All samples must be 1D arrays.")
data = np.vstack([x for x in samples.values()]).T data = np.vstack([x for x in samples.values()]).T
labels = names_to_latex(list(samples.keys()), for_corner=True) labels = names_to_latex(list(samples.keys()), for_corner=True)
return data, labels return data, labels, keys
def samples_to_getdist(samples, simname, catalogue=None): def samples_to_getdist(samples, label):
data, __ = samples_for_corner(samples) data, __, keys = samples_for_corner(samples)
names = list(samples.keys())
if catalogue is None:
label = simname_to_pretty(simname)
else:
label = catalogue
return MCSamples( return MCSamples(
samples=data, names=names, samples=data, names=keys,
labels=names_to_latex(names, for_corner=False), labels=names_to_latex(keys, for_corner=False),
label=label) label=label)

506
notebooks/flow/reconstruction_slice.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

586
notebooks/flow/selection.ipynb Normal file
View file

@ -0,0 +1,586 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Selection fitting "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from tqdm import trange\n",
"from h5py import File\n",
"from jax.random import PRNGKey\n",
"from numpyro.infer import MCMC, NUTS, init_to_median\n",
"from astropy.cosmology import FlatLambdaCDM \n",
"from corner import corner\n",
"\n",
"import csiborgtools\n",
"\n",
"%matplotlib inline\n",
"%load_ext autoreload\n",
"%autoreload 2\n",
"\n",
"Om0 = 0.3\n",
"H0 = 100\n",
"cosmo = FlatLambdaCDM(H0=H0, Om0=Om0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit parameters of the toy selection model\n",
"\n",
"Choose either CF4 TFR or SFI."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# with File(\"/mnt/extraspace/rstiskalek/catalogs/PV_compilation.hdf5\", 'r') as f:\n",
"# grp = f[\"SFI_gals\"]\n",
"# # # print(grp.keys())\n",
"# mag = grp[\"mag\"][...]\n",
"\n",
"\n",
"# with File(\"/mnt/extraspace/rstiskalek/catalogs/PV/CF4/CF4_TF-distances.hdf5\", 'r') as f:\n",
" # mag = f[\"w1\"][...]\n",
"# mag = mag[mag > 3]\n",
"\n",
"model = csiborgtools.flow.ToyMagnitudeSelection()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nuts_kernel = NUTS(model, init_strategy=init_to_median(num_samples=5000))\n",
"mcmc = MCMC(nuts_kernel, num_warmup=15_000, num_samples=15_000)\n",
"mcmc.run(PRNGKey(42), extra_fields=(\"potential_energy\",), mag=mag)\n",
"samples = mcmc.get_samples()\n",
"\n",
"mcmc.print_summary()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"keys = [\"alpha\", \"a\", \"m1\", \"m2\"]\n",
"data = np.vstack([samples[key] for key in keys]).T\n",
"labels = [r\"$\\alpha$\", r\"$a$\", r\"$m_1$\", r\"$m_2$\"]\n",
"\n",
"fig = corner(data, labels=labels, show_titles=True, smooth=True)\n",
"# fig.savefig(\"../../plots/selection_corner_CF4.png\", dpi=450)\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for key in keys:\n",
" print(f\"{key}: {np.mean(samples[key]):.3f} +/- {np.std(samples[key]):.3f}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mrange = np.linspace(mag.min(), mag.max(), 1000)\n",
"nsamples = len(samples[\"m1\"])\n",
"\n",
"indx = np.random.choice(nsamples, 500)\n",
"\n",
"y = [model.log_observed_pdf(mrange, samples[\"alpha\"][i], samples[\"m1\"][i], samples[\"m2\"][i], samples[\"a\"][i]) for i in indx]\n",
"y = np.asarray(y)\n",
"y = 10**y"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.hist(mag, bins=\"auto\", density=True, histtype=\"step\", color=\"blue\",\n",
" label=\"Data\", zorder=1)\n",
"\n",
"for i in range(100):\n",
" plt.plot(mrange, y[i], color=\"black\", alpha=0.25, lw=0.25)\n",
"\n",
"plt.xlabel(r\"$m$\")\n",
"plt.ylabel(r\"$p(m)$\")\n",
"plt.tight_layout()\n",
"\n",
"plt.savefig(\"../../plots/CF4_selection.png\", dpi=450)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Hubble \n",
"\n",
"$p(m) \\propto 10^{0.6 m}$ ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from scipy.integrate import quad\n",
"from scipy.interpolate import interp1d\n",
"\n",
"zmin=0.00001\n",
"zmax=5\n",
"z_range = np.linspace(zmin, zmax, 100000)\n",
"r_range = cosmo.comoving_distance(z_range).value\n",
"distmod_range = cosmo.distmod(z_range).value\n",
"r2mu = interp1d(r_range, distmod_range, kind=\"cubic\")\n",
"\n",
"\n",
"def schechter_LF(M, M0=-20.83, alpha=-1):\n",
" return 10**(0.4 * (M0 - M) * (alpha + 1)) * np.exp(-10**(0.4 * (M0 - M)))\n",
"\n",
"\n",
"def sample_schechter_LF(M0=-20.83, alpha=-1, Mfaint=-16, Mbright=-30, npoints=1):\n",
" norm = quad(schechter_LF, Mbright, Mfaint, args=(M0, alpha))[0]\n",
"\n",
" samples = np.full(npoints, np.nan)\n",
" for i in trange(npoints):\n",
" while np.isnan(samples[i]):\n",
" M = np.random.uniform(Mbright, Mfaint)\n",
" if np.random.uniform(0, 1) < schechter_LF(M, M0, alpha) / norm:\n",
" samples[i] = M\n",
"\n",
" return samples\n",
"\n",
"\n",
"def sample_radial_distance(rmax, npoints):\n",
" return rmax * np.random.rand(npoints)**(1/3)\n",
"\n",
"\n",
"# z = np.linspace(0.001, 0.15, 100000)\n",
"# r = cosmo.comoving_distance(z).value\n",
"# mu = cosmo.distmod(z).value\n",
"# \n",
"# \n",
"# drdmu = np.gradient(r, mu)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"rmax = 300\n",
"npoints = 5000\n",
"\n",
"r_150 = sample_radial_distance(100, npoints)\n",
"r_300 = sample_radial_distance(300, npoints)\n",
"r_1000 = sample_radial_distance(5000, npoints)\n",
"\n",
"mu_150 = r2mu(r_150)\n",
"mu_300 = r2mu(r_300)\n",
"mu_1000 = r2mu(r_1000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def p_hubble(m, a, b):\n",
" norm = np.log10(- 5 / np.log(1000) * (10**(3 / 5 * a) - 10**(3 / 5 * b)))\n",
" return 10**(0.6 * m - norm)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"M_LF = sample_schechter_LF(npoints=npoints)\n",
"\n",
"M_LF2 = sample_schechter_LF(npoints=npoints, M0=-20.83, alpha=-1.5)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"M = -20.3\n",
"\n",
"# m = mu + M\n",
"# x = np.linspace(11, m.max(), 1000)\n",
"# plt.plot(x, p_hubble(x, m.min(), m.max()) * 5.5, color=\"black\")\n",
"\n",
"# plt.hist(m, bins=\"auto\", density=True, histtype=\"step\", color=\"blue\",)\n",
"\n",
"\n",
"cols = [\"red\", \"green\", \"blue\"]\n",
"rmax = [150, 300, 1000]\n",
"# for i, mu in enumerate([mu_150, mu_300, mu_1000]):\n",
"for i, mu in enumerate([mu_150, mu_300, mu_1000]):\n",
" plt.hist(mu + M_LF, bins=\"auto\", density=True,\n",
" histtype=\"step\", color=cols[i], label=rmax[i])\n",
"\n",
" plt.hist(mu + M_LF2, bins=\"auto\", density=True,\n",
" histtype=\"step\", color=cols[i], label=rmax[i], ls=\"--\")\n",
"\n",
"\n",
"plt.hist(mag, bins=\"auto\", density=True, histtype=\"step\", color=\"black\", label=\"Data\")\n",
"\n",
"plt.yscale(\"log\")\n",
"# plt.axvline(r2mu(rmax) + M, c=\"red\")\n",
"plt.legend()\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"M = sample_schechter_LF(npoints=10000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.hist(x, bins=\"auto\", density=True, histtype=\"step\", color=\"blue\",)\n",
"# plt.yscale(\"log\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"\n",
"\n",
"yeuclidean = 10**(0.6 * mu)\n",
"ycomoving = r**2 * drdmu\n",
"\n",
"\n",
"\n",
"k = np.argmin(np.abs(mu - 35)) \n",
"\n",
"yeuclidean /= yeuclidean[k]\n",
"ycomoving /= ycomoving[k]\n",
"\n",
"\n",
"\n",
"plt.figure()\n",
"plt.plot(z, yeuclidean, label=\"Euclidean\")\n",
"plt.plot(z, ycomoving, label=\"Comoving\")\n",
"\n",
"# plt.yscale('log')\n",
"plt.xlabel(r\"$z$\")\n",
"plt.ylabel(r\"$p(\\mu)$\")\n",
"\n",
"plt.legend()\n",
"plt.tight_layout()\n",
"plt.savefig(\"../../plots/pmu_comoving_vs_euclidean.png\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from scipy.interpolate import interp1d\n",
"from scipy.integrate import quad\n",
"from scipy.stats import norm\n",
"\n",
"z = np.linspace(0.001, 0.1, 100000)\n",
"r = cosmo.comoving_distance(z).value\n",
"mu = cosmo.distmod(z).value\n",
"\n",
"\n",
"drdmu = np.gradient(r, mu)\n",
"\n",
"\n",
"\n",
"mu2drdmu = interp1d(mu, drdmu, kind='cubic')\n",
"mu2r = interp1d(mu, r, kind='cubic')\n",
"\n",
"\n",
"\n",
"def schechter_LF(M):\n",
" M0 = -20.83\n",
" alpha = -1\n",
" return 10**(0.4 * (M0 - M) * (alpha + 1)) * np.exp(-10**(0.4 * (M0 - M)))\n",
" \n",
"\n",
"\n",
"\n",
"def phi(M):\n",
" # return 1\n",
" # return schechter_LF(M)# * norm.pdf(M, loc=-22, scale=1)\n",
" loc = -22\n",
" std = 0.1\n",
"\n",
" return norm.pdf(M, loc=loc, scale=std)\n",
"\n",
" # if -22 < M < -21:\n",
" # return 1\n",
" # else:\n",
" # return 0\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"xrange = np.linspace(-24, -18, 1000)\n",
"\n",
"plt.figure()\n",
"plt.plot(xrange, schechter_LF(xrange))\n",
"# plt.yscale(\"log\")\n",
"plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mu_min = mu.min()\n",
"mu_max = mu.max()\n",
"\n",
"\n",
"m = 12\n",
"\n",
"\n",
"m_range = np.linspace(10, 16, 100)\n",
"y = np.full_like(m_range, np.nan)\n",
"for i in trange(len(m_range)):\n",
" m = m_range[i]\n",
" # y[i] = quad(lambda x: mu2drdmu(x) * mu2r(x)**2 * phi(m - x), mu_min, mu_max)[0]\n",
" y[i] = quad(lambda x: 10**(0.6 * x) * phi(m - x), mu_min, mu_max)[0]\n",
"\n",
"\n",
"\n",
"y_hubble = 10**(0.6 * m_range)\n",
"ycomoving = r**2 * drdmu\n",
"\n",
"\n",
"k = np.argmin(np.abs(m_range - 12))\n",
"\n",
"y_hubble /= y_hubble[k]\n",
"y /= y[k]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"mu_max - 18"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.plot(m_range, y, label=\"Numerical\")\n",
"plt.plot(m_range, y_hubble, label=\"Hubble\")\n",
"# plt.plot(mu, ycomoving, label=\"Comoving\")\n",
"\n",
"plt.xlabel(r\"$m$\")\n",
"plt.ylabel(r\"$p(m)$\")\n",
"plt.legend()\n",
"\n",
"# plt.yscale(\"log\")\n",
"plt.tight_layout()\n",
"# plt.xlim(10, 14)\n",
"\n",
"plt.savefig(\"../../plots/pm.png\", dpi=450)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simple simulation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"npoints = 10000\n",
"rmax = 30000\n",
"\n",
"# pos = np.random.uniform(-boxsize, boxsize, (npoints, 3))\n",
"\n",
"\n",
"r = rmax * np.random.rand(npoints)**(1/3)\n",
"\n",
"mu = 5 * np.log10(r) + 25\n",
"\n",
"# M = np.ones(npoints) * -22\n",
"# M = np.random.normal(-22, 100, npoints)\n",
"M = np.random.uniform(-24, -18, npoints)\n",
"\n",
"\n",
"m = mu + M"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def f(m, a, b):\n",
" norm = np.log10(- 5 / np.log(1000) * (10**(3 / 5 * a) - 10**(3 / 5 * b)))\n",
" return 10**(0.6 * m - norm)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure()\n",
"plt.hist(m, bins=\"auto\", density=True, histtype=\"step\")\n",
"m_range = np.linspace(m.min(), m.max(), 100)\n",
"# plt.plot(m_range, f(m_range, m.min(), m.max()))\n",
"# plt.yscale(\"log\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "venv_csiborg",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

2
scripts/field_prop/clear.sh Executable file
View file

@ -0,0 +1,2 @@
cm="rm *.out"
$cm

11
scripts/field_prop/field_los.py
View file

@ -396,10 +396,17 @@ if __name__ == "__main__":
parser.add_argument("--grid", type=int, help="Grid resolution.") parser.add_argument("--grid", type=int, help="Grid resolution.")
args = parser.parse_args() args = parser.parse_args()
rmax = 300 rmax = 200
dr = 0.5 if args.catalogue == "CF4_GroupAll":
dr = 1
else:
dr = 0.75
# smooth_scales = [0, 2, 4, 6, 8]
smooth_scales = [0] smooth_scales = [0]
print(f"Running catalogue {args.catalogue} for simulation {args.simname}.")
comm = MPI.COMM_WORLD comm = MPI.COMM_WORLD
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
nsims = get_nsims(args, paths) nsims = get_nsims(args, paths)

15
scripts/field_prop/field_los.sh
View file

@ -1,17 +1,26 @@
nthreads=1 nthreads=1
memory=64 memory=64
on_login=1 on_login=${1}
queue="berg" queue="berg"
env="/mnt/users/rstiskalek/csiborgtools/venv_csiborg/bin/python" env="/mnt/users/rstiskalek/csiborgtools/venv_csiborg/bin/python"
file="field_los.py" file="field_los.py"
nsims="-1" nsims="-1"
# These are only for CB
MAS="SPH" MAS="SPH"
grid=1024 grid=1024
for simname in "Lilow2024"; do if [ "$on_login" != "1" ] && [ "$on_login" != "0" ]
for catalogue in "CF4_TFR"; do then
echo "'on_login' (1) must be either 0 or 1."
exit 1
fi
# for simname in "csiborg1" "csiborg2_main" "csiborg2X" "Lilow2024" "Carrick2015" "CF4"; do
for simname in "csiborg2_main"; do
for catalogue in "2MTF" "SFI_gals" "CF4_TFR"; do
pythoncm="$env $file --catalogue $catalogue --nsims $nsims --simname $simname --MAS $MAS --grid $grid" pythoncm="$env $file --catalogue $catalogue --nsims $nsims --simname $simname --MAS $MAS --grid $grid"
if [ $on_login -eq 1 ]; then if [ $on_login -eq 1 ]; then
echo $pythoncm echo $pythoncm

95
scripts/field_prop/field_los_indranil_void.py Normal file
View file

@ -0,0 +1,95 @@
# Copyright (C) 2024 Richard Stiskalek
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
"""
MPI script to interpolate the density and velocity fields along the line of
sight.
"""
from os.path import join
import csiborgtools
import numpy as np
from astropy.coordinates import SkyCoord, angular_separation
from h5py import File
from mpi4py import MPI
from scipy.interpolate import RegularGridInterpolator
from field_los import get_los
def interpolate_indranil_void(kind, RA, dec, rmax, dr, dump_folder, catalogue):
fdir = "/mnt/extraspace/rstiskalek/catalogs"
if kind == "exp":
fname = join(fdir, "v_pec_EXP_IndranilVoid.dat")
elif kind == "gauss":
fname = join(fdir, "v_pec_GAUSS_IndranilVoid.dat")
else:
raise ValueError("Invalid void kind.")
# These are only velocities.
data = np.loadtxt(fname)
fname_out = join(dump_folder, f"los_{catalogue}_IndranilVoid_{kind}.hdf5")
r_grid = np.arange(0, 251)
phi_grid = np.arange(0, 181)
r_eval = np.arange(0, rmax, dr).astype(float) / 0.674
model_axis = SkyCoord(l=117, b=4, frame='galactic', unit='deg').icrs
coords = SkyCoord(ra=RA, dec=dec, unit='deg').icrs
# Get angular separation in degrees
phi = angular_separation(coords.ra.rad, coords.dec.rad,
model_axis.ra.rad, model_axis.dec.rad)
phi *= 180 / np.pi
# Get the interpolator
f = RegularGridInterpolator((r_grid, phi_grid), data.T)
# Get the dummy x-values to evaluate for each LOS
x_dummy = np.ones((len(r_eval), 2))
x_dummy[:, 0] = r_eval
result = np.full((len(RA), len(r_eval)), np.nan)
for i in range(len(RA)):
x_dummy[:, 1] = phi[i]
result[i] = f(x_dummy)
# Write the output, homogenous density.
density = np.ones_like(result)
print(f"Writing to `{fname_out}`.")
with File(fname_out, 'w') as f_out:
f_out.create_dataset("rdist_0", data=r_eval * 0.674)
f_out.create_dataset("density_0", data=density)
f_out.create_dataset("velocity_0", data=result)
###############################################################################
# Command line interface #
###############################################################################
if __name__ == "__main__":
kind = "exp"
rmax = 165
dr = 1
comm = MPI.COMM_WORLD
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
out_folder = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_los"
for catalogue in ["LOSS", "Foundation", "2MTF", "SFI_gals", "CF4_TFR", "CF4_GroupAll"]: # noqa
print(f"Running kind `{kind}` for catalogue `{catalogue}`.")
RA, dec = get_los(catalogue, "", comm).T
interpolate_indranil_void(
kind, RA, dec, rmax, dr, out_folder, catalogue)

2
scripts/flow/clear.sh Executable file
View file

@ -0,0 +1,2 @@
cm="rm *.out"
$cm

124
scripts/flow/flow_validation.py
View file

@ -34,7 +34,7 @@ def parse_args():
help="Simulation name.") help="Simulation name.")
parser.add_argument("--catalogue", type=str, required=True, parser.add_argument("--catalogue", type=str, required=True,
help="PV catalogues.") help="PV catalogues.")
parser.add_argument("--ksmooth", type=int, default=1, parser.add_argument("--ksmooth", type=int, default=0,
help="Smoothing index.") help="Smoothing index.")
parser.add_argument("--ksim", type=none_or_int, default=None, parser.add_argument("--ksim", type=none_or_int, default=None,
help="IC iteration number. If 'None', all IC realizations are used.") # noqa help="IC iteration number. If 'None', all IC realizations are used.") # noqa
@ -61,6 +61,7 @@ import sys
from os.path import join # noqa from os.path import join # noqa
import csiborgtools # noqa import csiborgtools # noqa
from csiborgtools import fprint # noqa
import jax # noqa import jax # noqa
from h5py import File # noqa from h5py import File # noqa
from numpyro.infer import MCMC, NUTS, init_to_median # noqa from numpyro.infer import MCMC, NUTS, init_to_median # noqa
@ -72,7 +73,7 @@ def print_variables(names, variables):
print(flush=True) print(flush=True)
def get_models(get_model_kwargs, toy_selection, verbose=True): def get_models(get_model_kwargs, mag_selection, verbose=True):
"""Load the data and create the NumPyro models.""" """Load the data and create the NumPyro models."""
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
folder = "/mnt/extraspace/rstiskalek/catalogs/" folder = "/mnt/extraspace/rstiskalek/catalogs/"
@ -111,7 +112,7 @@ def get_models(get_model_kwargs, toy_selection, verbose=True):
cat, fpath, paths, cat, fpath, paths,
ksmooth=ARGS.ksmooth) ksmooth=ARGS.ksmooth)
models[i] = csiborgtools.flow.get_model( models[i] = csiborgtools.flow.get_model(
loader, toy_selection=toy_selection[i], **get_model_kwargs) loader, mag_selection=mag_selection[i], **get_model_kwargs)
print(f"\n{'Num. radial steps':<20} {len(loader.rdist)}\n", flush=True) print(f"\n{'Num. radial steps':<20} {len(loader.rdist)}\n", flush=True)
return models return models
@ -127,9 +128,15 @@ def get_harmonic_evidence(samples, log_posterior, nchains_harmonic, epoch_num):
data, log_posterior, return_flow_samples=False, epochs_num=epoch_num) data, log_posterior, return_flow_samples=False, epochs_num=epoch_num)
def run_model(model, nsteps, nburn, model_kwargs, out_folder, sample_beta, def run_model(model, nsteps, nburn, model_kwargs, out_folder,
calculate_harmonic, nchains_harmonic, epoch_num, kwargs_print): calculate_harmonic, nchains_harmonic, epoch_num, kwargs_print,
fname_kwargs):
"""Run the NumPyro model and save output to a file.""" """Run the NumPyro model and save output to a file."""
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
fname = paths.flow_validation(out_folder, ARGS.simname, ARGS.catalogue,
**fname_kwargs)
try: try:
ndata = sum(model.ndata for model in model_kwargs["models"]) ndata = sum(model.ndata for model in model_kwargs["models"])
except AttributeError as e: except AttributeError as e:
@ -159,13 +166,6 @@ def run_model(model, nsteps, nburn, model_kwargs, out_folder, sample_beta,
neg_ln_evidence = jax.numpy.nan neg_ln_evidence = jax.numpy.nan
neg_ln_evidence_err = (jax.numpy.nan, jax.numpy.nan) neg_ln_evidence_err = (jax.numpy.nan, jax.numpy.nan)
fname = f"samples_{ARGS.simname}_{'+'.join(ARGS.catalogue)}_ksmooth{ARGS.ksmooth}.hdf5" # noqa
if ARGS.ksim is not None:
fname = fname.replace(".hdf5", f"_nsim{ARGS.ksim}.hdf5")
if sample_beta:
fname = fname.replace(".hdf5", "_sample_beta.hdf5")
fname = join(out_folder, fname) fname = join(out_folder, fname)
print(f"Saving results to `{fname}`.") print(f"Saving results to `{fname}`.")
with File(fname, "w") as f: with File(fname, "w") as f:
@ -209,7 +209,7 @@ def run_model(model, nsteps, nburn, model_kwargs, out_folder, sample_beta,
def get_distmod_hyperparams(catalogue, sample_alpha, sample_mag_dipole): def get_distmod_hyperparams(catalogue, sample_alpha, sample_mag_dipole):
alpha_min = -1.0 alpha_min = -1.0
alpha_max = 3.0 alpha_max = 10.0
if catalogue in ["LOSS", "Foundation", "Pantheon+", "Pantheon+_groups", "Pantheon+_zSN"]: # noqa if catalogue in ["LOSS", "Foundation", "Pantheon+", "Pantheon+_groups", "Pantheon+_zSN"]: # noqa
return {"e_mu_min": 0.001, "e_mu_max": 1.0, return {"e_mu_min": 0.001, "e_mu_max": 1.0,
@ -224,7 +224,6 @@ def get_distmod_hyperparams(catalogue, sample_alpha, sample_mag_dipole):
"a_mean": -21., "a_std": 5.0, "a_mean": -21., "a_std": 5.0,
"b_mean": -5.95, "b_std": 4.0, "b_mean": -5.95, "b_std": 4.0,
"c_mean": 0., "c_std": 20.0, "c_mean": 0., "c_std": 20.0,
"sample_curvature": False,
"a_dipole_mean": 0., "a_dipole_std": 1.0, "a_dipole_mean": 0., "a_dipole_std": 1.0,
"sample_a_dipole": sample_mag_dipole, "sample_a_dipole": sample_mag_dipole,
"alpha_min": alpha_min, "alpha_max": alpha_max, "alpha_min": alpha_min, "alpha_max": alpha_max,
@ -242,14 +241,27 @@ def get_distmod_hyperparams(catalogue, sample_alpha, sample_mag_dipole):
raise ValueError(f"Unsupported catalogue: `{ARGS.catalogue}`.") raise ValueError(f"Unsupported catalogue: `{ARGS.catalogue}`.")
def get_toy_selection(toy_selection, catalogue): def get_toy_selection(catalogue):
if not toy_selection: """Toy magnitude selection coefficients."""
if catalogue == "SFI_gals":
kind = "soft"
# m1, m2, a
coeffs = [11.467, 12.906, -0.231]
elif "CF4_TFR" in catalogue and "_i" in catalogue:
kind = "soft"
coeffs = [13.043, 14.423, -0.129]
elif "CF4_TFR" in catalogue and "w1" in catalogue:
kind = "soft"
coeffs = [11.731, 14.189, -0.118]
elif catalogue == "2MTF":
kind = "hard"
coeffs = 11.25
else:
fprint(f"found no selection coefficients for {catalogue}.")
return None return None
if catalogue == "SFI_gals": return {"kind": kind,
return [1.221e+01, 1.297e+01, -2.708e-01] "coeffs": coeffs}
else:
raise ValueError(f"Unsupported catalogue: `{ARGS.catalogue}`.")
if __name__ == "__main__": if __name__ == "__main__":
@ -262,43 +274,63 @@ if __name__ == "__main__":
# Fixed user parameters # # Fixed user parameters #
########################################################################### ###########################################################################
nsteps = 1000 # `None` means default behaviour
nburn = 500 nsteps = 10_000
zcmb_min = 0 nburn = 2_000
zcmb_min = None
zcmb_max = 0.05 zcmb_max = 0.05
nchains_harmonic = 10 nchains_harmonic = 10
num_epochs = 50 num_epochs = 50
inference_method = "bayes" inference_method = "mike"
calculate_harmonic = True if inference_method == "mike" else False mag_selection = None
maxmag_selection = None sample_alpha = True
sample_alpha = False sample_beta = None
sample_beta = True
sample_Vmono = False sample_Vmono = False
sample_mag_dipole = False sample_mag_dipole = False
toy_selection = True calculate_harmonic = False if inference_method == "bayes" else True
if toy_selection and inference_method == "mike": fname_kwargs = {"inference_method": inference_method,
raise ValueError("Toy selection is not supported with `mike` inference.") # noqa "smooth": ARGS.ksmooth,
"nsim": ARGS.ksim,
if nsteps % nchains_harmonic != 0: "zcmb_min": zcmb_min,
raise ValueError( "zcmb_max": zcmb_max,
"The number of steps must be divisible by the number of chains.") "mag_selection": mag_selection,
"sample_alpha": sample_alpha,
"sample_beta": sample_beta,
"sample_Vmono": sample_Vmono,
"sample_mag_dipole": sample_mag_dipole,
}
main_params = {"nsteps": nsteps, "nburn": nburn, main_params = {"nsteps": nsteps, "nburn": nburn,
"zcmb_min": zcmb_min, "zcmb_min": zcmb_min,
"zcmb_max": zcmb_max, "zcmb_max": zcmb_max,
"maxmag_selection": maxmag_selection, "mag_selection": mag_selection,
"calculate_harmonic": calculate_harmonic, "calculate_harmonic": calculate_harmonic,
"nchains_harmonic": nchains_harmonic, "nchains_harmonic": nchains_harmonic,
"num_epochs": num_epochs, "num_epochs": num_epochs,
"inference_method": inference_method, "inference_method": inference_method,
"sample_mag_dipole": sample_mag_dipole, "sample_mag_dipole": sample_mag_dipole,
"toy_selection": toy_selection} }
print_variables(main_params.keys(), main_params.values()) print_variables(main_params.keys(), main_params.values())
calibration_hyperparams = {"Vext_min": -1000, "Vext_max": 1000, if sample_beta is None:
sample_beta = ARGS.simname == "Carrick2015"
if mag_selection and inference_method != "bayes":
raise ValueError("Magnitude selection is only supported with `bayes` inference.") # noqa
if inference_method != "bayes":
mag_selection = [None] * len(ARGS.catalogue)
elif mag_selection is None or mag_selection:
mag_selection = [get_toy_selection(cat) for cat in ARGS.catalogue]
if nsteps % nchains_harmonic != 0:
raise ValueError(
"The number of steps must be divisible by the number of chains.")
calibration_hyperparams = {"Vext_min": -3000, "Vext_max": 3000,
"Vmono_min": -1000, "Vmono_max": 1000, "Vmono_min": -1000, "Vmono_max": 1000,
"beta_min": -1.0, "beta_max": 3.0, "beta_min": -10.0, "beta_max": 10.0,
"sigma_v_min": 1.0, "sigma_v_max": 750., "sigma_v_min": 1.0, "sigma_v_max": 750.,
"sample_Vmono": sample_Vmono, "sample_Vmono": sample_Vmono,
"sample_beta": sample_beta, "sample_beta": sample_beta,
@ -315,15 +347,11 @@ if __name__ == "__main__":
kwargs_print = (main_params, calibration_hyperparams, kwargs_print = (main_params, calibration_hyperparams,
*distmod_hyperparams_per_catalogue) *distmod_hyperparams_per_catalogue)
########################################################################### ###########################################################################
get_model_kwargs = {"zcmb_min": zcmb_min, "zcmb_max": zcmb_max, get_model_kwargs = {"zcmb_min": zcmb_min, "zcmb_max": zcmb_max}
"maxmag_selection": maxmag_selection} models = get_models(get_model_kwargs, mag_selection)
toy_selection = [get_toy_selection(toy_selection, cat)
for cat in ARGS.catalogue]
models = get_models(get_model_kwargs, toy_selection)
model_kwargs = { model_kwargs = {
"models": models, "models": models,
"field_calibration_hyperparams": calibration_hyperparams, "field_calibration_hyperparams": calibration_hyperparams,
@ -334,5 +362,5 @@ if __name__ == "__main__":
model = csiborgtools.flow.PV_validation_model model = csiborgtools.flow.PV_validation_model
run_model(model, nsteps, nburn, model_kwargs, out_folder, run_model(model, nsteps, nburn, model_kwargs, out_folder,
calibration_hyperparams["sample_beta"], calculate_harmonic, calculate_harmonic, nchains_harmonic, num_epochs, kwargs_print,
nchains_harmonic, num_epochs, kwargs_print) fname_kwargs)

21
scripts/flow/flow_validation.sh
View file

@ -1,5 +1,5 @@
#!/bin/bash #!/bin/bash
memory=7 memory=14
on_login=${1} on_login=${1}
queue=${2} queue=${2}
ndevice=1 ndevice=1
@ -37,12 +37,19 @@ else
fi fi
# for simname in "Lilow2024" "CF4" "CF4gp" "csiborg1" "csiborg2_main" "csiborg2X"; do for simname in "Carrick2015" "csiborg2_main"; do
for simname in "Carrick2015"; do # for simname in "csiborg2_main" "csiborg2X" ; do
for catalogue in "SFI_gals"; do # for simname in "Carrick2015" "Lilow2024" "csiborg2_main" "csiborg2X" "CF4"; do
# for catalogue in "CF4_TFR_i"; do # for simname in "Carrick2015" "csiborg2X" "csiborg2_main"; do
# for ksim in 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20; do # for simname in "Carrick2015"; do
# for catalogue in "LOSS" "Foundation" "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do
for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i"; do
# for catalogue in "2MTF" "SFI" "CF4_TFR_not2MTForSFI_i"; do
# for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i"; do
# for catalogue in "CF4_TFR_w1"; do
# for catalogue in "CF4_GroupAll"; do
for ksim in "none"; do for ksim in "none"; do
for ksmooth in 1 2 3 4; do
pythoncm="$env $file --catalogue $catalogue --simname $simname --ksim $ksim --ksmooth $ksmooth --ndevice $ndevice --device $device" pythoncm="$env $file --catalogue $catalogue --simname $simname --ksim $ksim --ksmooth $ksmooth --ndevice $ndevice --device $device"
if [ "$on_login" == "1" ]; then if [ "$on_login" == "1" ]; then
@ -65,3 +72,5 @@ for simname in "Carrick2015"; do
done done
done done
done done
done

12
scripts/flow/test_harmonic.py
View file

@ -1,4 +1,4 @@
from argparse import ArgumentParser, ArgumentTypeError from argparse import ArgumentParser
def parse_args(): def parse_args():
@ -10,7 +10,7 @@ def parse_args():
ARGS = parse_args() ARGS = parse_args()
# This must be done before we import JAX etc. # This must be done before we import JAX etc.
from numpyro import set_host_device_count, set_platform # noqa from numpyro import set_platform # noqa
set_platform(ARGS.device) # noqa set_platform(ARGS.device) # noqa
@ -30,8 +30,8 @@ def get_harmonic_evidence(samples, log_posterior, nchains_harmonic, epoch_num):
data, log_posterior, return_flow_samples=False, epochs_num=epoch_num) data, log_posterior, return_flow_samples=False, epochs_num=epoch_num)
ndim = 250 ndim = 150
nsamples = 100_000 nsamples = 50_000
nchains_split = 10 nchains_split = 10
loc = jnp.zeros(ndim) loc = jnp.zeros(ndim)
cov = jnp.eye(ndim) cov = jnp.eye(ndim)
@ -42,10 +42,6 @@ X = gen.multivariate_normal(loc, cov, size=nsamples)
samples = {f"x_{i}": X[:, i] for i in range(ndim)} samples = {f"x_{i}": X[:, i] for i in range(ndim)}
logprob = multivariate_normal(loc, cov).logpdf(X) logprob = multivariate_normal(loc, cov).logpdf(X)
neg_lnZ_laplace, neg_lnZ_laplace_error = csiborgtools.laplace_evidence(
samples, logprob, nchains_split)
print(f"neg_lnZ_laplace: {neg_lnZ_laplace} +/- {neg_lnZ_laplace_error}")
neg_lnZ_harmonic, neg_lnZ_harmonic_error = get_harmonic_evidence( neg_lnZ_harmonic, neg_lnZ_harmonic_error = get_harmonic_evidence(
samples, logprob, nchains_split, epoch_num=30) samples, logprob, nchains_split, epoch_num=30)