diff --git a/csiborgtools/flow/__init__.py b/csiborgtools/flow/__init__.py index 4d8a2c5..b4cd399 100644 --- a/csiborgtools/flow/__init__.py +++ b/csiborgtools/flow/__init__.py @@ -17,6 +17,7 @@ from .io import (DataLoader, get_model, read_absolute_calibration, from .flow_model import (PV_LogLikelihood, PV_validation_model, dist2redshift, # noqa Observed2CosmologicalRedshift, predict_zobs, # noqa project_Vext, stack_pzosmo_over_realizations) # noqa +from .mocks import mock_Carrick2MTF # noqa from .selection import ToyMagnitudeSelection # noqa from .void_model import (load_void_data, interpolate_void, select_void_h, # noqa mock_void) # noqa diff --git a/csiborgtools/flow/flow_model.py b/csiborgtools/flow/flow_model.py index 7bb358d..7d15881 100644 --- a/csiborgtools/flow/flow_model.py +++ b/csiborgtools/flow/flow_model.py @@ -170,28 +170,6 @@ class BaseFlowValidationModel(ABC): self._setattr_as_jax(names, values) - def _set_abs_calibration_params(self, abs_calibration_params): - self.with_absolute_calibration = abs_calibration_params is not None - - if abs_calibration_params is None: - self.with_absolute_calibration = False - return - - self.calibration_distmod = jnp.asarray( - abs_calibration_params["calibration_distmod"][..., 0]) - self.calibration_edistmod = jnp.asarray( - abs_calibration_params["calibration_distmod"][..., 1]) - self.data_with_calibration = jnp.asarray( - abs_calibration_params["data_with_calibration"]) - self.data_wo_calibration = ~self.data_with_calibration - - # Calculate the log of the number of calibrators. Where there is no - # calibrator set the number of calibrators to 1 to avoid log(0) and - # this way only zeros are being added. - length_calibration = abs_calibration_params["length_calibration"] - length_calibration[length_calibration == 0] = 1 - self.log_length_calibration = jnp.log(length_calibration) - def _set_radial_spacing(self, r_xrange, Omega_m): cosmo = FlatLambdaCDM(H0=H0, Om0=Omega_m) @@ -356,7 +334,8 @@ 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, c_mean, c_std, alpha_min, alpha_max, sample_alpha, - a_dipole_mean, a_dipole_std, sample_a_dipole, name): + a_dipole_mean, a_dipole_std, sample_a_dipole, + sample_curvature, name): """Sample Tully-Fisher calibration parameters.""" e_mu = sample(f"e_mu_{name}", Uniform(e_mu_min, e_mu_max)) a = sample(f"a_{name}", Normal(a_mean, a_std)) @@ -367,7 +346,11 @@ def sample_TFR(e_mu_min, e_mu_max, a_mean, a_std, b_mean, b_std, ax, ay, az = 0.0, 0.0, 0.0 b = sample(f"b_{name}", Normal(b_mean, b_std)) - c = sample(f"c_{name}", Normal(c_mean, c_std)) + + if sample_curvature: + c = sample(f"c_{name}", Normal(c_mean, c_std)) + else: + c = 0. alpha = sample_alpha_bias(name, alpha_min, alpha_max, sample_alpha) @@ -495,8 +478,6 @@ class PV_LogLikelihood(BaseFlowValidationModel): Errors on the observed redshifts. calibration_params : dict Calibration parameters of each object. - abs_calibration_params : dict - Absolute calibration parameters. mag_selection : dict Magnitude selection parameters, optional. r_xrange : 1-dimensional array @@ -513,12 +494,17 @@ class PV_LogLikelihood(BaseFlowValidationModel): Whether to directly sample the distance without numerical marginalisation. in which case the tracers can be coupled by a covariance matrix. By default `False`. + with_homogeneous_malmquist : bool, optional + Whether to include the homogeneous Malmquist bias. By default `True`. + with_inhomogeneous_malmquist : bool, optional + Whether to include the inhomogeneous Malmquist bias. By default `True`. """ def __init__(self, los_density, los_velocity, RA, dec, z_obs, e_zobs, - calibration_params, abs_calibration_params, mag_selection, - r_xrange, Omega_m, kind, name, void_kwargs=None, - wo_num_dist_marginalisation=False): + calibration_params, mag_selection, r_xrange, Omega_m, kind, + name, void_kwargs=None, wo_num_dist_marginalisation=False, + with_homogeneous_malmquist=True, + with_inhomogeneous_malmquist=True): if e_zobs is not None: e2_cz_obs = jnp.asarray((SPEED_OF_LIGHT * e_zobs)**2) else: @@ -550,13 +536,14 @@ class PV_LogLikelihood(BaseFlowValidationModel): self._setattr_as_jax(names, values) self._set_calibration_params(calibration_params) - self._set_abs_calibration_params(abs_calibration_params) self._set_radial_spacing(r_xrange, Omega_m) self.kind = kind self.name = name self.Omega_m = Omega_m self.wo_num_dist_marginalisation = wo_num_dist_marginalisation + self.with_homogeneous_malmquist = with_homogeneous_malmquist + self.with_inhomogeneous_malmquist = with_inhomogeneous_malmquist self.norm = - self.ndata * jnp.log(self.num_sims) if mag_selection is not None: @@ -771,9 +758,14 @@ class PV_LogLikelihood(BaseFlowValidationModel): "marginalising the true distance.") # Calculate p(r) (Malmquist bias). Shape is (ndata, nxrange) - log_ptilde = log_ptilde_wo_bias( - self.mu_xrange[None, :], mu[:, None], e2_mu[:, None], - self.log_r2_xrange[None, :]) + if self.with_homogeneous_malmquist: + log_ptilde = log_ptilde_wo_bias( + self.mu_xrange[None, :], mu[:, None], e2_mu[:, None], + self.log_r2_xrange[None, :]) + else: + log_ptilde = log_ptilde_wo_bias( + self.mu_xrange[None, :], mu[:, None], e2_mu[:, None], + 0.) if self.is_void_data: rLG = field_calibration_params["rLG"] @@ -785,7 +777,9 @@ class PV_LogLikelihood(BaseFlowValidationModel): # Inhomogeneous Malmquist bias. Shape: (nsims, ndata, nxrange) alpha = distmod_params["alpha"] - log_ptilde = log_ptilde[None, ...] + alpha * log_los_density + log_ptilde = log_ptilde[None, ...] + if self.with_inhomogeneous_malmquist: + log_ptilde += alpha * log_los_density ptilde = jnp.exp(log_ptilde) # Normalization of p(r). Shape: (nsims, ndata) @@ -802,29 +796,51 @@ class PV_LogLikelihood(BaseFlowValidationModel): ptilde *= likelihood_zobs( self.z_obs[None, :, None], zobs, e2_cz[None, :, None]) - if self.with_absolute_calibration: - raise NotImplementedError( - "Absolute calibration not implemented for this model. " - "Use `PV_LogLikelihood_NoDistMarg` instead.") - # Integrate over the radial distance. Shape: (nsims, ndata) ll = jnp.log(simpson(ptilde, x=self.r_xrange, axis=-1)) ll -= jnp.log(pnorm) return ll0 + jnp.sum(logsumexp(ll, axis=0)) + self.norm else: - if field_calibration_params["sample_h"]: - raise NotImplementedError( - "Sampling of h is not yet implemented.") - e_mu = jnp.sqrt(e2_mu) - # True distance modulus, shape is `(n_data)`` + # True distance modulus, shape is `(n_data)`. If we have absolute + # calibration, then this distance modulus assumes a particular h. with plate("plate_mu", self.ndata): mu_true = sample("mu", Normal(mu, e_mu)) - # True distance and redshift, shape is `(n_data)`. - r_true = distmod2dist(mu_true, self.Omega_m) - z_true = distmod2redshift(mu_true, self.Omega_m) + # Likelihood of the true distance modulii given the calibration. + if field_calibration_params["sample_h"]: + raise RuntimeError( + "Sampling of 'h' has not yet been thoroughly tested.") + h = field_calibration_params["h"] + + # Now, the rest of the code except the calibration likelihood + # uses the distance modulus in units of h + mu_true_h = mu_true + 5 * jnp.log10(h) + + # Calculate the log-likelihood of the calibration, but the + # shape is `(n_calibrators, n_data)`. Where there is no data + # we set the likelihood to 0 (or the log-likelihood to -inf) + ll_calibration = jnp.where( + self.is_finite_calibrator, + normal_logpdf(self.mu_calibration, mu_true[None, :], + self.e_mu_calibration), + -jnp.inf) + + # Now average out over the calibrators, however only if the + # there is at least one calibrator. If there isn't, then we + # just assing a log-likelihood of 0. + ll_calibration = jnp.where( + self.any_calibrator, + logsumexp(ll_calibration, axis=0) - jnp.log(self.counts_calibrators), # noqa + 0.) + else: + mu_true_h = mu_true + + # True distance and redshift, shape is `(n_data)`. The distance + # here is in units of `Mpc / h``. + r_true = distmod2dist(mu_true_h, self.Omega_m) + z_true = distmod2redshift(mu_true_h, self.Omega_m) if self.is_void_data: raise NotImplementedError( @@ -840,26 +856,29 @@ class PV_LogLikelihood(BaseFlowValidationModel): alpha = distmod_params["alpha"] # Normalisation of p(mu), shape is `(n_sims, n_data, n_rad)` - pnorm = ( - + self.log_r2_xrange[None, None, :] - + alpha * log_los_density_grid - + normal_logpdf( - self.mu_xrange[None, :], mu[:, None], e_mu[:, None])[None, ...]) # noqa + pnorm = normal_logpdf( + self.mu_xrange[None, :], mu[:, None], e_mu[:, None])[None, ...] + if self.with_homogeneous_malmquist: + pnorm += self.log_r2_xrange[None, None, :] + if self.with_inhomogeneous_malmquist: + pnorm += alpha * log_los_density_grid + pnorm = jnp.exp(pnorm) # Now integrate over the radial steps. Shape is `(nsims, ndata)`. # No Jacobian here because I integrate over distance, not the # distance modulus. pnorm = simpson(pnorm, x=self.r_xrange, axis=-1) - # Jacobian |dr / dmu|_(mu_true), shape is `(n_data)`. - jac = jnp.abs(distmod2dist_gradient(mu_true, self.Omega_m)) + # Jacobian |dr / dmu|_(mu_true_h), shape is `(n_data)`. + jac = jnp.abs(distmod2dist_gradient(mu_true_h, self.Omega_m)) # Calculate unnormalized log p(mu). Shape is (nsims, ndata) - ll = ( - + jnp.log(jac)[None, :] - + (2 * jnp.log(r_true) - self.log_r2_xrange_mean)[None, :] - + alpha * log_density - ) + ll = 0.0 + if self.with_homogeneous_malmquist: + ll += (+ jnp.log(jac) + + (2 * jnp.log(r_true) - self.log_r2_xrange_mean)) + if self.with_inhomogeneous_malmquist: + ll += alpha * log_density # Subtract the normalization. Shape remains (nsims, ndata) ll -= jnp.log(pnorm) @@ -871,13 +890,13 @@ class PV_LogLikelihood(BaseFlowValidationModel): zobs *= 1 + vrad / SPEED_OF_LIGHT zobs -= 1. + # Add the log-likelihood of observed redshifts. Shape remains + # `(nsims, ndata)` ll += log_likelihood_zobs( self.z_obs[None, :], zobs, e2_cz[None, :]) - if self.with_absolute_calibration: - raise NotImplementedError( - "Absolute calibration not implemented for this model. " - "Use `PV_LogLikelihood_NoDistMarg` instead.") + if field_calibration_params["sample_h"]: + ll += ll_calibration[None, :] return ll0 + jnp.sum(logsumexp(ll, axis=0)) + self.norm diff --git a/csiborgtools/flow/io.py b/csiborgtools/flow/io.py index 694c657..d8c60ec 100644 --- a/csiborgtools/flow/io.py +++ b/csiborgtools/flow/io.py @@ -58,7 +58,8 @@ class DataLoader: def __init__(self, simname, ksim, catalogue, catalogue_fpath, paths, ksmooth=None, store_full_velocity=False, verbose=True): fprint("reading the catalogue,", verbose) - self._cat = self._read_catalogue(catalogue, catalogue_fpath) + self._cat, self._absmag_calibration = self._read_catalogue( + catalogue, catalogue_fpath) self._catname = catalogue fprint("reading the interpolated field.", verbose) @@ -132,6 +133,15 @@ class DataLoader: """The distance indicators catalogue (structured array).""" return self._cat[self._mask] + @property + def absmag_calibration(self): + """Returns the absolute magnitude calibration with masking applied.""" + if self._absmag_calibration is None: + return None + + return {key: val[:, self._mask] + for key, val in self._absmag_calibration.items()} + @property def catname(self): """Catalogue name.""" @@ -187,6 +197,8 @@ class DataLoader: fpath = paths.field_los(simname, "Pantheon+") elif "CF4_TFR" in catalogue: fpath = paths.field_los(simname, "CF4_TFR") + elif "Carrick2MTFmock" in catalogue: + fpath = paths.field_los(simname, "2MTF") else: fpath = paths.field_los(simname, catalogue) @@ -217,6 +229,8 @@ class DataLoader: return rdist, los_density, los_velocity def _read_catalogue(self, catalogue, catalogue_fpath): + absmag_calibration = None + if catalogue == "A2": with File(catalogue_fpath, 'r') as f: dtype = [(key, np.float32) for key in f.keys()] @@ -262,6 +276,20 @@ class DataLoader: arr = np.empty(len(mock_data["RA"]), dtype=dtype) for key in mock_data.keys(): arr[key] = mock_data[key] + elif "Carrick2MTFmock" in catalogue: + with File(catalogue_fpath, 'r') as f: + keys_skip = ["mu_calibration", "e_mu_calibration"] + dtype = [(key, np.float32) for key in f.keys() + if key not in keys_skip] + arr = np.empty(len(f["RA"]), dtype=dtype) + for key in f.keys(): + if key not in keys_skip: + arr[key] = f[key][:] + + absmag_calibration = { + "mu_calibration": f["mu_calibration"][...], + "e_mu_calibration": f["e_mu_calibration"][...]} + elif "UPGLADE" in catalogue: with File(catalogue_fpath, 'r') as f: dtype = [(key, np.float32) for key in f.keys()] @@ -286,7 +314,7 @@ class DataLoader: else: raise ValueError(f"Unknown catalogue: `{catalogue}`.") - return arr + return arr, absmag_calibration ############################################################################### @@ -322,7 +350,7 @@ def radial_velocity_los(los_velocity, ra, dec): def read_absolute_calibration(kind, data_length, calibration_fpath): """ Read the absolute calibration for the CF4 TFR sample from LEDA but - preprocessed by me. + preprocessed by me. Missing values are replaced with NaN. Parameters ---------- @@ -335,13 +363,13 @@ def read_absolute_calibration(kind, data_length, calibration_fpath): Returns ------- - data : 3-dimensional array of shape (data_length, max_calib, 2) + mu : 2-dimensional array of shape `(ncalib, ngalaxies)` Absolute calibration data. - with_calibration : 1-dimensional array of shape (data_length) - Whether the sample has a calibration. - length_calibration : 1-dimensional array of shape (data_length) - Number of calibration points per sample. + e_mu : 2-dimensional array of shape `(ncalib, ngalaxies)` + Uncertainties of the absolute calibration. """ + raise RuntimeError("The read-in functions are not guaranteed to work " + "properly.") data = {} with File(calibration_fpath, 'r') as f: for key in f[kind].keys(): @@ -355,14 +383,14 @@ def read_absolute_calibration(kind, data_length, calibration_fpath): max_calib = max(len(val) for val in data.values()) out = np.full((data_length, max_calib, 2), np.nan) - with_calibration = np.full(data_length, False) - length_calibration = np.full(data_length, 0) for i in data.keys(): out[int(i), :len(data[i]), :] = data[i] - with_calibration[int(i)] = True - length_calibration[int(i)] = len(data[i]) - return out, with_calibration, length_calibration + # Unpack from this the distsance modulus and its uncertainty. + mu = out[:, :, 0].T + e_mu = out[:, :, 1].T + + return mu, e_mu def mask_fields(density, velocity, mask, return_none): @@ -422,8 +450,9 @@ def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None, loader._field_rdist = rdist - if absolute_calibration is not None and "CF4_TFR_" not in kind: - raise ValueError("Absolute calibration supported only for the CF4 TFR sample.") # noqa + if absolute_calibration is not None and not ("CF4_TFR_" in kind or "Carrick2MTFmock" in kind): # noqa + raise ValueError("Absolute calibration supported only for either " + "the CF4 TFR sample or Carrick 2MTF mocks.") if kind in ["LOSS", "Foundation"]: keys = ["RA", "DEC", "z_CMB", "mB", "x1", "c", "e_mB", "e_x1", "e_c"] @@ -442,7 +471,7 @@ def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None, model = PV_LogLikelihood( los_overdensity, los_velocity, RA[mask], dec[mask], zCMB[mask], e_zCMB, calibration_params, - None, mag_selection, loader.rdist, loader._Omega_m, "SN", + mag_selection, loader.rdist, loader._Omega_m, "SN", name=kind, void_kwargs=void_kwargs, wo_num_dist_marginalisation=wo_num_dist_marginalisation) elif "Pantheon+" in kind: @@ -476,26 +505,62 @@ def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None, model = PV_LogLikelihood( los_overdensity, los_velocity, RA[mask], dec[mask], zCMB[mask], e_zCMB[mask], calibration_params, - None, mag_selection, loader.rdist, loader._Omega_m, "SN", + mag_selection, loader.rdist, loader._Omega_m, "SN", name=kind, void_kwargs=void_kwargs, wo_num_dist_marginalisation=wo_num_dist_marginalisation) - elif kind in ["SFI_gals", "2MTF", "SFI_gals_masked"] or "IndranilVoidTFRMock" in kind: # noqa + elif kind in ["SFI_gals", "2MTF", "SFI_gals_masked"] or "IndranilVoidTFRMock" in kind or "Carrick2MTFmock" in kind: # noqa 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) mask = (zCMB < zcmb_max) & (zCMB > zcmb_min) + if "Carrick2MTFmock" in kind: + # For the mock we only want to select objects with the '2M++' + # volume. + mask &= loader.cat["r"] < 130 + # The mocks are generated without Malmquist. + fprint("disabling homogeneous and inhomogeneous Malmquist bias for the mock.") # noqa + with_homogeneous_malmquist = False + with_inhomogeneous_malmquist = False + else: + with_homogeneous_malmquist = True + with_inhomogeneous_malmquist = True + calibration_params = {"mag": mag[mask], "eta": eta[mask], "e_mag": e_mag[mask], "e_eta": e_eta[mask]} + # Append the calibration data + if "Carrick2MTFmock" in kind: + absmag_calibration = loader.absmag_calibration + + # The shape of these is (`ncalibrators, nobjects`). + mu_calibration = absmag_calibration["mu_calibration"][:, mask] + e_mu_calibration = absmag_calibration["e_mu_calibration"][:, mask] + # Auxiliary parameters. + m = np.isfinite(mu_calibration) + + # NumPyro refuses to start if any inputs are not finite, so we + # replace with some ficutive mean and very large standard + # deviation. + mu_calibration[~m] = 0.0 + e_mu_calibration[~m] = 1000.0 + + calibration_params["mu_calibration"] = mu_calibration + calibration_params["e_mu_calibration"] = e_mu_calibration + calibration_params["is_finite_calibrator"] = m + calibration_params["counts_calibrators"] = np.sum(m, axis=0) + calibration_params["any_calibrator"] = np.any(m, axis=0) + los_overdensity, los_velocity = mask_fields( los_overdensity, los_velocity, mask, void_kwargs is not None) model = PV_LogLikelihood( los_overdensity, los_velocity, - RA[mask], dec[mask], zCMB[mask], None, calibration_params, None, + RA[mask], dec[mask], zCMB[mask], None, calibration_params, mag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind, void_kwargs=void_kwargs, - wo_num_dist_marginalisation=wo_num_dist_marginalisation) + wo_num_dist_marginalisation=wo_num_dist_marginalisation, + with_homogeneous_malmquist=with_homogeneous_malmquist, + with_inhomogeneous_malmquist=with_inhomogeneous_malmquist) elif "CF4_TFR_" in kind: # The full name can be e.g. "CF4_TFR_not2MTForSFI_i" or "CF4_TFR_i". band = kind.split("_")[-1] @@ -532,38 +597,39 @@ def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None, else: mask &= Qs == 5 - # Read the absolute calibration - if absolute_calibration is not None: - CF4_length = len(RA) - distmod, with_calibration, length_calibration = read_absolute_calibration( # noqa - "Cepheids", CF4_length, calibration_fpath) - - distmod = distmod[mask] - with_calibration = with_calibration[mask] - length_calibration = length_calibration[mask] - fprint(f"found {np.sum(with_calibration)} galaxies with absolute calibration.") # noqa - - distmod = distmod[with_calibration] - length_calibration = length_calibration[with_calibration] - - abs_calibration_params = { - "calibration_distmod": distmod, - "data_with_calibration": with_calibration, - "length_calibration": length_calibration} - else: - abs_calibration_params = None - calibration_params = {"mag": mag[mask], "eta": eta[mask], "e_mag": e_mag[mask], "e_eta": e_eta[mask]} + # Read the absolute calibration + mu_calibration, e_mu_calibration = read_absolute_calibration( + absolute_calibration, len(RA), calibration_fpath) + + # The shape of these is (`ncalibrators, nobjects`). + mu_calibration = mu_calibration[:, mask] + e_mu_calibration = e_mu_calibration[:, mask] + # Auxiliary parameters. + m = np.isfinite(mu_calibration) + + # NumPyro refuses to start if any inputs are not finite, so we + # replace with some ficutive mean and very large standard + # deviation. + mu_calibration[~m] = 0.0 + e_mu_calibration[~m] = 1000.0 + + calibration_params["mu_calibration"] = mu_calibration + calibration_params["e_mu_calibration"] = e_mu_calibration + calibration_params["is_finite_calibrator"] = m + calibration_params["counts_calibrators"] = np.sum(m, axis=0) + calibration_params["any_calibrator"] = np.any(m, axis=0) + los_overdensity, los_velocity = mask_fields( los_overdensity, los_velocity, mask, void_kwargs is not None) model = PV_LogLikelihood( los_overdensity, los_velocity, RA[mask], dec[mask], z_obs[mask], None, calibration_params, - abs_calibration_params, mag_selection, loader.rdist, - loader._Omega_m, "TFR", name=kind, void_kwargs=void_kwargs, + mag_selection, loader.rdist, loader._Omega_m, "TFR", name=kind, + void_kwargs=void_kwargs, wo_num_dist_marginalisation=wo_num_dist_marginalisation) elif kind in ["CF4_GroupAll"]: # Note, this for some reason works terribly. @@ -583,7 +649,7 @@ def get_model(loader, zcmb_min=None, zcmb_max=None, mag_selection=None, model = PV_LogLikelihood( los_overdensity, los_velocity, - RA[mask], dec[mask], zCMB[mask], None, calibration_params, None, + RA[mask], dec[mask], zCMB[mask], None, calibration_params, mag_selection, loader.rdist, loader._Omega_m, "simple", name=kind, void_kwargs=void_kwargs, wo_num_dist_marginalisation=wo_num_dist_marginalisation) diff --git a/csiborgtools/flow/mocks.py b/csiborgtools/flow/mocks.py new file mode 100644 index 0000000..f091083 --- /dev/null +++ b/csiborgtools/flow/mocks.py @@ -0,0 +1,150 @@ +# 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. +"""Mock data generators.""" +import numpy as np + +from ..field.interp import evaluate_cartesian_regular +from ..params import SPEED_OF_LIGHT +from ..utils import cartesian_to_radec, radec_to_cartesian, radec_to_galactic +from .cosmography import distmod2dist, distmod2redshift + +############################################################################### +# Mock Quijote observations # +############################################################################### + + +def mock_Carrick2MTF(velocity_field, boxsize, RA_2MTF, DEC_2MTF, + a_TF=-22.8, b_TF=-7.2, sigma_TF=0.25, sigma_v=100., + Vext=[150.0, 10.0, -100.0], h=1.0, beta=0.4, + mean_eta=0.069, std_eta=0.078, mean_e_eta=0.012, + mean_mag=10.31, std_mag=0.83, mean_e_mag=0.044, + sigma_calibration=0.05, calibration_max_percentile=10, + calibration_rand_fraction=0.5, nrepeat_calibration=1, + seed=42, Om0=0.3, verbose=True, **kwargs): + """ + Mock TFR catalogue build against the Carrick velocity field and the + 2MTF sky distribution to avoid recomputing the LOS velocities. + """ + nsamples = len(RA_2MTF) + + # Convert Vext from ICRS to Galactic coordinates. + Vext = np.asarray(Vext).reshape(1, 3) + Vext_mag, Vext_RA, Vext_DEC = cartesian_to_radec(Vext).reshape(-1, ) + Vext_l, Vext_b = radec_to_galactic(Vext_RA, Vext_DEC) + Vext_galactic = np.asanyarray([Vext_mag, Vext_l, Vext_b]).reshape(1, 3) + Vext = radec_to_cartesian(Vext_galactic).reshape(-1, ) + + truths = {"a": a_TF, "b": b_TF, "e_mu": sigma_TF, "sigma_v": sigma_v, + "Vext": Vext, + "mean_eta": mean_eta, "std_eta": std_eta, + "mean_mag": mean_mag, "std_mag": std_mag, + } + + gen = np.random.default_rng(seed) + + # The Carrick box is in the Galactic coordinates. + l, b = radec_to_galactic(RA_2MTF, DEC_2MTF) + gal_phi = np.deg2rad(l) + gal_theta = np.pi / 2 - np.deg2rad(b) + + # Sample the linewidth of each galaxy from a Gaussian distribution to mimic + # the MNR procedure. + eta_true = gen.normal(mean_eta, std_eta, nsamples) + eta_obs = gen.normal(eta_true, mean_e_eta) + + # Subtract the mean of the observed linewidths, so that they are + # centered around zero. For consistency subtract from both observed + # and true values. + eta_mean_sampled = np.mean(eta_obs) + eta_true -= eta_mean_sampled + eta_obs -= eta_mean_sampled + + # Sample the magnitude from some Gaussian distribution to replicate MNR. + mag_true = gen.normal(mean_mag, std_mag, nsamples) + mag_obs = gen.normal(mag_true, mean_e_mag) + + # Calculate the 'true' distance modulus and redshift from the TFR distance. + + # If h != 1, then these distance modulii are in physical units. + mu_TFR = mag_true - (a_TF + b_TF * eta_true) + mu_true = gen.normal(mu_TFR, sigma_TF) + # This is the distance modulus in units of little h. + mu_true_h = mu_true + 5 * np.log10(h) + + # Select a fraction of nearby galaxies. + mu_max = np.percentile(mu_true, calibration_max_percentile) + ks = np.where(mu_true < mu_max)[0] + nsel = int(calibration_rand_fraction * len(ks)) + if verbose: + print(f"Assigning calibration to {nsel}/{nsamples} galaxies.") + + mu_calibration = np.full((nrepeat_calibration, nsamples), np.nan) + e_mu_calibration = np.full((nrepeat_calibration, nsamples), np.nan) + + for n in range(nrepeat_calibration): + ks_n = gen.choice(ks, nsel, replace=False) + + mu_calibration[n, ks_n] = gen.normal(mu_true[ks_n], sigma_calibration) + e_mu_calibration[n, ks_n] = np.ones(len(ks_n)) * sigma_calibration + + # Convert the true distance modulus to true distance and cosmological + # redshift. The distance is in Mpc/h because the box is in Mpc / h. + r = distmod2dist(mu_true_h, Om0) + zcosmo = distmod2redshift(mu_true_h, Om0) + + # Calculate the Cartesian coordinates of each galaxy. This is initially + # centered at (0, 0, 0). + pos = r * np.asarray([ + np.sin(gal_theta) * np.cos(gal_phi), + np.sin(gal_theta) * np.sin(gal_phi), + np.cos(gal_theta)]) + pos = pos.T + pos_box = pos / boxsize + 0.5 + + vel = evaluate_cartesian_regular( + velocity_field[0], velocity_field[1], velocity_field[2], + pos=pos_box, smooth_scales=None, method="cubic") + vel = beta * np.vstack(vel).T + + for i in range(3): + vel[:, i] += Vext[i] + + # Compute the radial velocity. + Vr = np.sum(vel * pos, axis=1) / np.linalg.norm(pos, axis=1) + + # The true redshift of the source. + zCMB_true = (1 + zcosmo) * (1 + Vr / SPEED_OF_LIGHT) - 1 + zCMB_obs = gen.normal(zCMB_true, sigma_v / SPEED_OF_LIGHT) + + # These galaxies will be masked out when LOS is read it because they are + # too far away. + distance_mask = r < 125 + truths["distance_mask"] = distance_mask + + sample = {"RA": RA_2MTF, + "DEC": DEC_2MTF, + "z_CMB": zCMB_obs, + "eta": eta_obs, + "mag": mag_obs, + "e_eta": np.ones(nsamples) * mean_e_eta, + "e_mag": np.ones(nsamples) * mean_e_mag, + "mu_true": mu_true, + "mu_TFR": mu_TFR, + "mu_calibration": mu_calibration, + "e_mu_calibration": e_mu_calibration, + "r": r, + } + + return sample, truths diff --git a/notebooks/flow/mock_carrick.ipynb b/notebooks/flow/mock_carrick.ipynb new file mode 100644 index 0000000..25fa761 --- /dev/null +++ b/notebooks/flow/mock_carrick.ipynb @@ -0,0 +1,530 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Copyright (C) 2024 Richard Stiskalek\n", + "# This program is free software; you can redistribute it and/or modify it\n", + "# under the terms of the GNU General Public License as published by the\n", + "# Free Software Foundation; either version 3 of the License, or (at your\n", + "# option) any later version.\n", + "#\n", + "# This program is distributed in the hope that it will be useful, but\n", + "# WITHOUT ANY WARRANTY; without even the implied warranty of\n", + "# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General\n", + "# Public License for more details.\n", + "#\n", + "# You should have received a copy of the GNU General Public License along\n", + "# with this program; if not, write to the Free Software Foundation, Inc.,\n", + "# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from tqdm import trange\n", + "from joblib import dump\n", + "from h5py import File\n", + "\n", + "import csiborgtools\n", + "\n", + "%matplotlib inline\n", + "%load_ext autoreload\n", + "%autoreload 2\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "fpath = \"/mnt/extraspace/rstiskalek/catalogs/PV_compilation.hdf5\"\n", + "\n", + "with File(fpath, 'r') as f:\n", + " RA_2MTF = f[\"2MTF/RA\"][...]\n", + " DEC_2MTF = f[\"2MTF/DEC\"][...]" + ] + }, + { + "cell_type": "code", + "execution_count": 123, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/usr/local/shared/python/3.11.7/lib/python3.11/pty.py:89: RuntimeWarning: os.fork() was called. os.fork() is incompatible with multithreaded code, and JAX is multithreaded, so this will likely lead to a deadlock.\n", + " pid, fd = os.forkpty()\n", + " 0%| | 0/1 [00:00" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure()\n", + "plt.hist(mock[\"mu_calibration\"][0, m], bins=\"auto\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 126, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNfElEQVR4nO3df1iUZb4/8PeAzKDCDIHiwAqKP0IRsaTEOZWZQqIdM8XdMv2q5Wq64CZuu0ZnW7M9LbRe3zV32+jHMWu/RuzakVzriMdf0GaoiLGCFKssauUMrBIDQgzI3N8/2JllZAbmmRlgeHi/rmuunOe5n3vum9GLT/ePz60QQggQERERDXA+/d0AIiIiIk9gUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLAzp7wb0FbPZjKtXryIwMBAKhaK/m0NEREROEEKgsbER4eHh8PHpfixm0AQ1V69eRURERH83g4iIiFzw1VdfYfTo0d2WGTRBTWBgIICOH4pare7n1hAREZEzGhoaEBERYf093p1BE9RYppzUajWDGiIiogHGmaUjXChMREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZGDTJ94iIiMh17WaB09V1qG1sQWigP2ZEBcPXx7vOUmRQQ0RERN3KL9dj24EK6I0t1mvBw5V45I5wJMVovSbAUQghRH83oi80NDRAo9HAaDTymAQiIiIn5ZfrsWHPWXQXLIRp/LF1YQySY8M8/vlSfn9zTQ0RERHZ1W4W2HagotuABgD0xhas33MWvzxwHkVV19Fu7p/xEk4/ERERkV2nq+tsppx6suvEJew6calXR266w5EaIiIisqu20fmApjODsQUb9pxFfrnewy3qHoMaIiIisis00N+l5yyTT9sOVPTpVBSDGiIiIrJrRlQwwjT+cGVfk0DHWpvT1XWebpZDDGqIiIgGmXazQFHVdewv/abbhb2+PgpsXRgDAC4FNoDrU1iukBTUZGdnIy4uDmq1Gmq1GjqdDgcPHrTef/PNNzF79myo1WooFArU19f3WOfYsWOhUCi6vFJTU61lZs+e3eX++vXrpTSdiIiI0LFF+96Xj2HZWyfxdG4plr11Eve+fMzh+pfk2DBkr5gOzTA/lz7P1SksV0gKakaPHo2srCyUlJTgzJkzmDNnDhYtWoTz588DAJqbm5GcnIznnnvO6TqLi4uh1+utr8OHDwMAvv/979uUW7t2rU25X//611KaTkRENOhZcs7cuqNJ78TC3vrmNkmfpUBH/poZUcGuNNUlkrZ0L1y40Ob9Sy+9hOzsbJw8eRJTpkzBpk2bAAAFBQVO1zly5Eib91lZWRg/fjzuv/9+m+vDhg2DVquV0lwiIiL6p55yzggAGfvKkBSjtckObHlOCsvTWxfG9GmmYZfX1LS3tyM3NxdNTU3Q6XQeaUxrayv27NmDJ598EgqF7Q/hvffew4gRIxAbG4uMjAw0Nzd3W5fJZEJDQ4PNi4iIaLByJufMt81tePXYBcnP3Uqr8Uf2iul9nqdGcvK9srIy6HQ6tLS0ICAgAHl5eYiJifFIYz788EPU19dj9erVNtcff/xxjBkzBuHh4Th37hy2bNmCyspK7Nu3z2FdmZmZ2LZtm0faRURENNA5u2D37U+rkTZnonWExdnn0h4Yj4mjAvv1sEvJQU10dDRKS0thNBrxwQcfYNWqVSgsLPRIYLNr1y7Mnz8f4eHhNtfXrVtn/fPUqVMRFhaGuXPnoqqqCuPHj7dbV0ZGBjZv3mx939DQgIiICLfbSERENBA5u2DX2HITrx67iKcTJ0p67p4JI6EbH+Jy+zxBclCjVCoxYcIEAEB8fDyKi4uxc+dOvPHGG2415PLlyzhy5Ei3oy8WCQkJAICLFy86DGpUKhVUKpVbbSIiIhpI2s0Cp6vrUNvY0mXEZEZUMIKG+qH+u54X/O448jcAQNqcCdZcNd1NQfX1gmBH3D77yWw2w2Qyud2Q3bt3IzQ0FA899FCPZUtLSwEAYWF9O1dHRETkrfLL9dh2oMIm+Oh8BpOvjwJP3BNlDVh6suPI3/D+6ct44eEpeHhaGN74pNph2YenhfXLdNOtJAU1GRkZmD9/PiIjI9HY2IicnBwUFBTg0KFDAACDwQCDwYCLFy8C6Fh/ExgYiMjISAQHd0Rwc+fOxeLFi5GWlmat12w2Y/fu3Vi1ahWGDLFtUlVVFXJycrBgwQKEhITg3LlzSE9Px6xZsxAXF+dW54mIiOTAslX71p1NljOYsldMR1KMFmYh7cgCQ4MJ6/ecxXClb7fl3jv1FSZp1dBqhvbbehpAYlBTW1uLlStXQq/XQ6PRIC4uDocOHUJSUhIA4PXXX7dZnDtr1iwAHaMwlsW/VVVVuHbtmk29R44cwZUrV/Dkk092+UylUokjR47glVdeQVNTEyIiIpCSkoKf//znkjpKREQkR91t1Rbo2F6dsa8Mz/73OdR/d9Olz2hqbe/2/g3TTaT/6a8A0G8ndAOAQgiJYdsA1dDQAI1GA6PRCLVa3d/NISIi8oiiqutY9tbJ/m6GlWWMxlNbuqX8/ubZT0RERAOYwfhdfzfBRn+d0A0wqCEiIhqw8sv1+OXHX/R3M7rojxO6AQ/sfiIiIqK+52hxsKdZppNc+Zy+PKEb4EgNERHRgNPTOU6eYglo1s2Kcun5vjyhG2BQQ0RENOA4ex6Tws2d1ZYznDIWxOD1FdOhGercBE9/nNANMKghIiIacJyd1nF3f/N3rf/aAp4Uo8VQv56Dmv46oRtgUENERDTg9NW0Tv13N7F+z1nkl+txuroOhoaeg6ng4cp+OaEbYFBDRETkldrNAkVV17G/9BsUVV232R5tOY+pr8ZBtu4vx6cX/+FU2Z8/NLlfAhqAu5+IiIi8jjPnOG1dGIMNe85CAdd2JklR09iK3x+vcqqsVjO0l1vjGEdqiIiI+kh3oy8Wlq3aty4EtpzjlF+uBwAkx4Yhe8V0aDV9u8PIkf5aHNwZR2qIiIj6QE+jL4Bz5zhtO1CBpBgtfH0USI4NQ1KMFqer6/B6YRUK/+bcFJGn9efi4M44UkNERNTLnB196Wmrtr1Mvb4+CsyICkbZN0an2/PUrCiEeXCEx7L1u7/W0lhwpIaIiKgXSRl9cXar9tX677DrL3/H5bpmjAkehttHBaKuqbXH5xQAfv/4dCyIC8PPkifjdHUdjlQYsOvEJQk9spX2wASkJ93eryM0FgxqiIiIepGU0ZdL15qdqvOZD/5qk4PG2XBi9b+NwYK4jtEUXx8FdONDoBsfgrujgvHsvjLUN7c5WdO/3DNhhFcENACDGiIiol7l7OiLwfgd3j99xamytybVc3b3U9Awpd3rlrU5J/9+HUVV1wEIJIwNwU//+xxqGlrs1q9Ax7RTfy4MvhXX1BAREfUiZxPlXbvR6lRyO3e8f/qK3R1XQMfIzT0TRuCZedF4Zt4k3Bc9Ei88HAOg60iQtywMvhWDGiIiol7SbhYwC4GgoX4OyygABA3zw6vHL/Z6ewwNJptFxj1xtG3cWxYG34rTT0RERL3A3hbuW1kS57mylsVVzk6HWXTeNl7b2ILQwI4pJ28aobFgUENERORhli3cPa11GaVWoeWmuU+DGlfOjbIsKvZ2nH4iIiLyoO62cFsMU/riPxZMxval0yQFNO6MjXhDxt/exqCGiIjIg3rawg0Aza3teOl/vsDG9z+XVHdiTGi395NiQqHAwFnY62kMaoiIiDxIypqV+u+kjdKUf9OAtfdF4da4xEfRkSX4rZV3D6iFvZ7GNTVEREQe5MqaFWdYkvTNmTQKP503Cf+v6JI1o/D/0Y2FckjHOMVAWtjraQxqiIiIPGhGVDDCNP4wGO0nrXNXbWMLlEN8sOa+cQ7LDJSFvZ7G6SciIiIP8vVRYOtC+0nrPKG3RoLkgEENERGRhzlKWueu4SpfWe9echeDGiIiol6QHBuGT7fMQers8R6rs8nUjsMVBo/VJzcMaoiIiHqJr48C904c6dE6tx2ocHh+02DHoIaIiKgXzYgKRvBw+6dju0JvbJF0ftNgwqCGiIioG+1mgaKq69hf+g2Kqq5LHiXx9VHgPxfFerRNUs9vGiy4pZuIiMgBe4dShmn8sXVhDJJjw9BuFk7lg1kQF4anvo7CG59Ue6Rd3AFln6SRmuzsbMTFxUGtVkOtVkOn0+HgwYPW+2+++SZmz54NtVoNhUKB+vr6Hut84YUXoFAobF6TJk2yKdPS0oLU1FSEhIQgICAAKSkpqKmpkdJ0IiIiSSyHUt565IHB2IINe84iLecs4n95GMveOomnc0ux7K2TuPflY8gv19ut787I29xu02A4v8kdkoKa0aNHIysrCyUlJThz5gzmzJmDRYsW4fz58wCA5uZmJCcn47nnnpPUiClTpkCv11tfn376qc399PR0HDhwAHv37kVhYSGuXr2KJUuWSPoMIiIiZ3V3KKX45+ujc/ouxxxYAp5bAxtLfZ4g9/Ob3CFp+mnhwoU271966SVkZ2fj5MmTmDJlCjZt2gQAKCgokNaIIUOg1Wrt3jMajdi1axdycnIwZ84cAMDu3bsxefJknDx5EjNnzpT0WURERD1x5lBKewQ6RlO2HahAUozWGny4Wl9nwcP98KvFU2V/fpM7XF4o3N7ejtzcXDQ1NUGn07nViAsXLiA8PBzjxo3D8uXLceXKFeu9kpIStLW1ITEx0Xpt0qRJiIyMRFFRkcM6TSYTGhoabF5ERETOcGchruWMps47lI64mVtG7T8EJzMSGdD0QHJQU1ZWhoCAAKhUKqxfvx55eXmIiYlxuQEJCQl45513kJ+fj+zsbFRXV+O+++5DY2MjAMBgMECpVCIoKMjmuVGjRsFgcPyXJDMzExqNxvqKiIhwuY1ERDS4eGIhriUwyi/XY9eJS27V9f340dYDK8kxyT+h6OholJaW4tSpU9iwYQNWrVqFigrX5wnnz5+P73//+4iLi8O8efPwP//zP6ivr8ef/vQnl+sEgIyMDBiNRuvrq6++cqs+IiIaPCyHUrqzcuXStWaPraVJjLG/RINsSd7SrVQqMWHCBABAfHw8iouLsXPnTrzxxhseaVBQUBBuv/12XLx4EQCg1WrR2tqK+vp6m9Gampoah+twAEClUkGlUnmkTURENLhYDqVcv+esy3W8cuRvAODWWhoFAC13OznN7bEss9kMk8nkibYAAG7cuIGqqiqEhXXMG8bHx8PPzw9Hjx61lqmsrMSVK1fcXstDRETkSHJsGJ68Z6xbdez+zPm8NLeOClnec7eT8yQFNRkZGfjkk09w6dIllJWVISMjAwUFBVi+fDmAjvUvpaWl1lGWsrIylJaWoq7uX4ul5s6di1dffdX6/plnnkFhYSEuXbqEzz77DIsXL4avry+WLVsGANBoNFizZg02b96M48ePo6SkBE888QR0Oh13PhERUa9KcmPaRwCob27rsZyFZpifzXutxh/ZK6ZzcbAEkqafamtrsXLlSuj1emg0GsTFxeHQoUNISkoCALz++uvYtm2btfysWbMAdGzBXr16NQCgqqoK165ds5b5+uuvsWzZMly/fh0jR47Evffei5MnT2LkyH8dALZjxw74+PggJSUFJpMJ8+bNw2uvveZyp4mIiAD0mBHYsrbGYGyxm7PGGcOUvmhube+2jAKA/xAfvPfDBFy7Yeo2OzE5phBCDIqjPhsaGqDRaGA0GqFWq/u7OURE1IfsBS+HKwzdHoFgkV+ud2ttzXCVL5pM3Qc1Fu+vnQnd+BCXP0uOpPz+5tlPREQka/bObwoa5md3asiSEdiT0z5NpnbMj9XiYHnPuWp4UKV7uOmdiIhky9H5TY7WulimLrYdqEC7WXhsS/b4kcOdKseDKt3DoIaIiGSn3Sxw4sI1PPvfZZLXwnTOCOyJ4w0AQDduRLd5b3hQpWdw+omIiGTF3nSTKzw1FRSm8cfM8SHYujAGG/achQKwCbS4ddtzOFJDRESy4Wi6yRWhgf4emQ6yBCvJsWHIXjEdWo1tndy67TkcqSEiIlmwrH/x1Jbeb5taMS9W69aW7oVxWptgJTk2DEkx2m63kZPrOFJDRESy4Kn1Lxa//LhjgfDWhR2HNrsSdhw4Z0B+ud7mmq+PArrxIVh0x/egGx/CgMaDGNQQEZEseHo7tN7Ygh2H/wbNUCV+//idXaaNgm7JAGyPAv/aSUW9j9NPREQ0YLWbBU7+/TqKqq7j62+bPV7/q8cv4tXjFxGm8cfzD8XgtuFKm2mjV49dxI5/HlxpT+edVEyq1/sY1BAR0YCUX67Hs/vKJJ2v5CqDsQWpOR1J+Rbd8T3rdeN3rU49z6R6fYPTT0RENOBYji7oi4AG6JqUD+gYJfqw9KpTzzOpXt9gUENERF6h3SxQVHUd+0u/QVHVdYfrUNrNAi/82f0sv1J1nkoCOhYm1zX1PFITPNyPSfX6CKefiIio39lLmGfvcEmgI5gwNPTfdI5lKsnZKaXFd3yPO5z6CEdqiIioXzlKmGc5XPLWLdHurk+5bdgQl7ZnW1imkpydUkqM0brxaSQFgxoiIuo33SXMs7eOBXB/fcq2f48FID3vzK3nM82ICuZ5Tl6GQQ0REfWbnhLm3bqOBegIJrRq1wObfzSZ7B5X0B175zP5+igcJubjeU79g0ENERH1G2enkjqX8/VR4IWHY1z+zMt1zUiODcOnW+YgPfF2p55xdD4Tz3PyLlwoTERE/cbZqaRbyyXHhuH1FdOx+U9/RXNru6TPHBM8zPrn3OIr3ZYNGuqH3y+fjpnjHB9nwPOcvAeDGiIi6jeWdSmODoxUoGPUo/O6lHazwOnqOnzX2g5/P19JQY2PAvg/urEAnDsrqv67NvgoFD0GKJbznKh/MaghIqJ+4+ujwMPTwvDGJ9UOy3Rel2Jv67cUa++LgnJIx8oLV6a+yLsxqCEion6TX67Hm90ENOtmRVnXpVi2frtyNKSPoiOgyVjwr7U4rk59kffiQmEiIuoX3W3ntvjjma/RbhZOlXVkmJ8Pzm9LtgloAG7JliMGNURE1C+cWtPS3IYtH/wVJ/9+3eUpp+Y2M85e/rbLdW7Jlh8GNURE1C+cXavywdlv8KM9JW59VmpO18zEALdkyw3X1BARUb+QslbF2HLTrc+q/64NG/acdZhrhluy5YFBDRER9YsZUcEIGuqH+u/a+uwztx2oQFKMtkvAwi3Z8sDpJyIi6he+Pgqs/rcxffZ59o5cIHlhUENERP0iv1yP3OKvJD8XNNTP5n2Yxh9PzYrqct0R5p2RL04/ERERgH9l6nVmXYmUsvbKf9vUitQc13LO/P7x6fDxUXT57Fm3h2L5f53q8XnmnZEvBjVERGQ3U2+Yxh9bF8Z0WVgrpayj8j4KSA5oLEcmzBxv/xymmeNCJB+5QPLC6SciokHOkqn31jwwBmMLNuyx3QotpWx35c0SIxpn8sYw7wxJCmqys7MRFxcHtVoNtVoNnU6HgwcPWu+/+eabmD17NtRqNRQKBerr63usMzMzE3fffTcCAwMRGhqKRx55BJWVlTZlZs+eDYVCYfNav369lKYTEZEd3WXqtVzbdqCix6y+t5btqW6pnM0bw7wzg5uk6afRo0cjKysLEydOhBAC7777LhYtWoTPP/8cU6ZMQXNzM5KTk5GcnIyMjAyn6iwsLERqairuvvtu3Lx5E8899xwefPBBVFRUYPjw4dZya9euxYsvvmh9P2zYMHvVERGRBD1l9b11x5CzZXXjQ5zKGOyM5x+ajNX3RDk9wsK8M4OXpKBm4cKFNu9feuklZGdn4+TJk5gyZQo2bdoEACgoKHC6zvz8fJv377zzDkJDQ1FSUoJZs2ZZrw8bNgxarVZKc4mIqAe9cVK1payndhmNCFRJDkiYd2ZwcnlNTXt7O3Jzc9HU1ASdTuexBhmNRgBAcLDtQq733nsPI0aMQGxsLDIyMtDc3NxtPSaTCQ0NDTYvIiKyJeWkaqmnWntqlxF3K5GzJO9+Kisrg06nQ0tLCwICApCXl4eYmJieH3SC2WzGpk2bcM899yA2NtZ6/fHHH8eYMWMQHh6Oc+fOYcuWLaisrMS+ffsc1pWZmYlt27Z5pF1ERHJlOana2R1DUsr2VHdPuFuJpJI8UhMdHY3S0lKcOnUKGzZswKpVq1BRUeGRxqSmpqK8vBy5ubk219etW4d58+Zh6tSpWL58Of7whz8gLy8PVVVVDuvKyMiA0Wi0vr76SnqCJyIiuZOyY8iV3UWP3R3hMABSAFh731i77eJuJXKF5KBGqVRiwoQJiI+PR2ZmJqZNm4adO3e63ZC0tDR89NFHOH78OEaPHt1t2YSEBADAxYsXHZZRqVTWXVqWFxERdSVlx5CzZfPL9bj35WPYceSC3c+0lP+Ph6bg9RXTEcbdSuQBbiffM5vNMJlMLj8vhMDGjRuRl5eHgoICREVF9fhMaWkpACAsjH/ZiYg8QcqOoZ7KWnLTOJpySk+8HWlzJljLc7cSeYqkoCYjIwPz589HZGQkGhsbkZOTg4KCAhw6dAgAYDAYYDAYrCMoZWVlCAwMRGRkpHXh79y5c7F48WKkpaUB6JhyysnJwf79+xEYGAiDwQAA0Gg0GDp0KKqqqpCTk4MFCxYgJCQE586dQ3p6OmbNmoW4uDiP/SCIiAY7KTuGbi3bbhYoqroOg/E7/PLjLxwGNAoAucVXkDZngsufTeSIpKCmtrYWK1euhF6vh0ajQVxcHA4dOoSkpCQAwOuvv26zONeyJXv37t1YvXo1AKCqqgrXrl2zlsnOzgbQkWCvM8szSqUSR44cwSuvvIKmpiZEREQgJSUFP//5zyV3loiIPM/eMQiO3JrLhsiTFEIITyR79HoNDQ3QaDQwGo1cX0NE5CE9TTU5svOxO7Doju/1SptIXqT8/ubZT0RE5BJ3jkFg7hnqDTylm4iIXOLKMQjMPUO9iUENERH1qN0suuxOknoMAnPPUG9jUENERN2ytxA4TOOPx+6OkFSPVuOPrQtjmHuGeg2DGiIicsjRQmCDsQU7jlxA0DA/GJvbHK6rCR7uh+f/fQq0auaeod7HoIaIaJCwN4XUXZDR3UJgAdujEhT/vNb5PQD8avFUjsxQn2FQQ0Q0CDiaQupuOqinhcACQH1zG9ITb0du8RWbspxqov7AoIaISOa6m0LasOeswzOWnF0IPHbEMHy6ZQ6POaB+x6CGiEjGeppCAoBn95UhUOWHmeNDbAKREcNVTn1GaKA/jzkgr8Dke0REMuZMLpn65jYs33UK9758DPnlegAdozs/2fvXbp9ToGMKizlnyFtwpIaISMak5JKxTEetmxWFNz+p7jZTMHPOkDfiSA0RkYxJOY5A/PP11l+6D2gAYJRa5XAtDlF/YVBDRCRT7WYBs1kgaKifpOfMThzm9H9/cAcDGvI6nH4iIpIhe1u4PenaDVOv1EvkDgY1REQy42gLtyfxlG3yRpx+IiKSke62cFv0tKzXR+G4DHc8kTdjUENEJCPObOG2BDy3Bi6Kf77W3hfl8D7AHU/kvRjUEBHJiLNbuNfcMxZaje0Uklbjj+wV05GxIAbZK6Y7vM8FwuStuKaGiEhGnF3rkhijxXMPxTg82iA5NgxJMVoefUADCoMaIiIZmREVjDCNPwzGFrvrahToGHFxZk0Mjz6ggYZBDRHRANJuFj2Onjx2dyR2HPlbl2c7r4k5XGGQfGo3kbdjUENENEDYyz1jCUSSYrR49dgF7D5xCfXftdl9XvvPsgBcOrWbyNsxqCEiGgAc5Z4xGFuwfs9ZDFP6orm13eHz6YkTkTZnIgDg3pePOTy1WwFg24EKJMVouX6GBhzufiIi8nLd5Z6xXOsuoFEAyC3+CkDPW74FAL2xBaer61xuL1F/YVBDROTlnMk9053OgYqzW76lnO5N5C0Y1BAReTlPBRi1jS24dK3ZqbI8BoEGIgY1RERezlMBxojhKrx/+kqP5bRqFY9BoAGJQQ0RkZez5J5xZ9lumMYfUACGhp5HfZbNiOQiYRqQGNQQEXk5Xx+FdSu2q6HG1oUxuHbD5FTZsSOGu/gpRP2LQQ0RUR9oNwsUVV3H/tJvUFR1He3m7s7R7io5NszueUw9CRrmh9f/mXfG2WksrqehgYp5aoiIell3SfOkJLnrfB6ToaEFv/zoPOqa7CfaA4AAlS9OP5cI5ZCO/3/15BEKRN5I0khNdnY24uLioFaroVarodPpcPDgQev9N998E7Nnz4ZarYZCoUB9fb1T9f7+97/H2LFj4e/vj4SEBJw+fdrmfktLC1JTUxESEoKAgACkpKSgpqZGStOJiPqFJWnerVuyLdl788v1kuqznMekVft3G9AAwA1TO0ouf2vzrKNprM5HKHA9DQ1UkoKa0aNHIysrCyUlJThz5gzmzJmDRYsW4fz58wCA5uZmJCcn47nnnnO6zj/+8Y/YvHkztm7dirNnz2LatGmYN28eamtrrWXS09Nx4MAB7N27F4WFhbh69SqWLFkipelERH3OmaR52w5USJ6KApzf5n1rOUfTWFqNP49HoAFPIYSQ/q+pk+DgYGzfvh1r1qyxXisoKMADDzyAb7/9FkFBQd0+n5CQgLvvvhuvvvoqAMBsNiMiIgIbN27Es88+C6PRiJEjRyInJwdLly4FAHz55ZeYPHkyioqKMHPmTKfa2dDQAI1GA6PRCLVa7VpniYgkKKq6jmVvneyx3PtrZ0o+DdvZutMemIB7JozocvClMwdjEnkDKb+/XV4o3N7ejtzcXDQ1NUGn07lUR2trK0pKSpCYmPivBvn4IDExEUVFRQCAkpIStLW12ZSZNGkSIiMjrWXsMZlMaGhosHkREfUlZ0dTjlQYJNft7DbvV49fxLK3TuLel4/ZTHVZprEW3fE96MaHMKAhWZAc1JSVlSEgIAAqlQrr169HXl4eYmJiXPrwa9euob29HaNGjbK5PmrUKBgMHf/IDQYDlEpllxGfzmXsyczMhEajsb4iIiJcaiMRkas7l5zN3vvHM19JnoKSus3b1TU8RAOJ5KAmOjoapaWlOHXqFDZs2IBVq1ahoqKiN9rmloyMDBiNRuvrq6++6u8mEdEAlF+ux70vH8Oyt07i6dxSu6Me9oKe/HI9XjnyN6c+44apHa8euyC5bVK2ebu7hodoIJC8pVupVGLChAkAgPj4eBQXF2Pnzp144403JH/4iBEj4Ovr22UnU01NDbRaLQBAq9WitbUV9fX1NqM1ncvYo1KpoFKpJLeJiMjCsnPp1hDAMuqRvWI6AHTZrq1Vq9By02x3gbAju09cQtqciZKngTpv8z5x8R949XiVw7KdD7aUuoaHaCBwO/me2WyGyeRclspbKZVKxMfH4+jRozb1HT161LpOJz4+Hn5+fjZlKisrceXKFZfX8hAR9cSZnUsZ+8qw3t527QYT6pu73259q/rv2vDOiWqXkvNZ1sdMHBXoVHmewE1yJWmkJiMjA/Pnz0dkZCQaGxuRk5ODgoICHDp0CEDH+heDwYCLFy8C6Fh/ExgYiMjISAQHdyRzmjt3LhYvXoy0tDQAwObNm7Fq1SrcddddmDFjBl555RU0NTXhiSeeAABoNBqsWbMGmzdvRnBwMNRqNTZu3AidTuf0ziciIqlOV9d1CVY6EwC+lRi49OSXH39h/bMryfmYMZgGO0lBTW1tLVauXAm9Xg+NRoO4uDgcOnQISUlJAIDXX38d27Zts5afNWsWAGD37t1YvXo1AKCqqgrXrl2zlnn00Ufxj3/8A7/4xS9gMBhwxx13ID8/32bx8I4dO+Dj44OUlBSYTCbMmzcPr732msudJiLqSX+PZnSe4nI2sGHGYBrs3M5TM1AwTw0RSeFsHpjeZAlCPt0yx+m1NpZ1QABsAhvL00ywRwNNn+SpISKSs/gxt6G/U7d0XtjrLGYMpsGMB1oSEdlRcvlbeMvOZ6lTYZ13RDFjMA0mDGqIiOyQEkgoAEnbty3U/kPQ0HKzx3KuLOy17IgiGkw4/UREZIezgUR64kSnkt91lvrAeLy/dibO/Dyp26MOFOjYBcWFvUTOYVBDRGRHT2crWQKOtDkT8emWOXh/7UysmBnpVN0TRgZANz4EyiE+Do86sLzfujCG00ZETmJQQ0Sy4eoZTfZ0d7bSrQGHZaonKmS4U3XXNbVa/8yFvUSewzU1RCQL+eX6LscVuJLArjNLwNHlGAQH9QYHOHc0y63luLCXyDMY1BDRgOfMGU3uBDbOBhxatXNra+yV48JeIvcxqCGiAa2nM5oU6DhwMilG6/LIh7MBh2UdTnfHK3DhL1Hv4ZoaIhrQnDmjSWoCO1dZ1uEoYH8djgJc+EvUmxjUENGA5mw+mb46y4kLf4n6D6efiGhA88aTqbnwl6h/MKghogHNW0+m5sJfor7H6SciGtCk5JMhInljUENEXkdqEj131rF4MmEfEfUvTj8RkVdxNYmelHUs7WaB09V1OFJhQF7pN6hrapP0WUTknRRCiEHxvyUNDQ3QaDQwGo1Qq9X93RwissNREj2L1x6/Ewviwt3+jFuDps4sIRB3KhF5Bym/vzn9REReobskehZp73+O/zmnd/kzLEFTT3ltgI6EfZyKIhpYGNQQkVfoKYkeAJgF8KOcs8gvlx7YOBM0WfRlwj4i8hwGNUTkFaQkx3NlFMWZoMmdNhFR/2NQQ0ReQUpyPFdGUVwJUPoyYR8RuY9BDRF5BUsSPWdJDVKkBCgK8OBJooGIQQ0ReYXOSfScIXUUxRI0OZuCjwn7iAYeBjVE5DWSY8Pw2uN3ortYwtVRlO4yD3cWxoMniQYsJt8jIq+yIC4cr0KBH+Wc7XLP3WMPLJmHb81TEzJciUV3hCMpRsuDJ4kGMCbfIyKv5GpmYWdYMgrzBG0i7yfl9zeDGiLyWgw+iIgZhYlowGNAQ0RScU0NEXmd3px6IiL54kgNEXkVR+czGYwt2LDHtSMSiGhwYFBDRF6ju/OZeNAkEfVEUlCTnZ2NuLg4qNVqqNVq6HQ6HDx40Hq/paUFqampCAkJQUBAAFJSUlBTU9NtnQqFwu5r+/bt1jJjx47tcj8rK0tiV4nI2/V0PhMPmiSi7kgKakaPHo2srCyUlJTgzJkzmDNnDhYtWoTz588DANLT03HgwAHs3bsXhYWFuHr1KpYsWdJtnXq93ub19ttvQ6FQICUlxabciy++aFNu48aNErtKRP2l3SxQVHUd+0u/QVHVdYcjLc4efcCDJonIHkkLhRcuXGjz/qWXXkJ2djZOnjyJ0aNHY9euXcjJycGcOXMAALt378bkyZNx8uRJzJw5026dWq3W5v3+/fvxwAMPYNy4cTbXAwMDu5QlIu8nZdGvs0cfXGs0YX/pN9wVRUQ2XF5T097ejtzcXDQ1NUGn06GkpARtbW1ITEy0lpk0aRIiIyNRVFTkVJ01NTX4+OOPsWbNmi73srKyEBISgjvvvBPbt2/HzZs3u63LZDKhoaHB5kVEfUvqot/4Mbd1e0SCxS8//gJP55Zi2Vsnce/Lx7h4mIgAuBDUlJWVISAgACqVCuvXr0deXh5iYmJgMBigVCoRFBRkU37UqFEwGAxO1f3uu+8iMDCwy5TVj3/8Y+Tm5uL48eN46qmn8Ktf/Qo/+9nPuq0rMzMTGo3G+oqIiJDUTyJyjyuLfksufwupa4C5K4qILCTnqYmOjkZpaSmMRiM++OADrFq1CoWFhR5pzNtvv43ly5fD3992CHrz5s3WP8fFxUGpVOKpp55CZmYmVCqV3boyMjJsnmtoaGBgQ9SHpCz61Y0PAeDaWhmBjjOhth2oQFKM1vrZTNpHNPhIDmqUSiUmTJgAAIiPj0dxcTF27tyJRx99FK2traivr7cZrampqXFqLcxf/vIXVFZW4o9//GOPZRMSEnDz5k1cunQJ0dHRdsuoVCqHAQ8R9T5XFv06u6bmVpYA6dVjF5FbfIVJ+4gGKbfz1JjNZphMJsTHx8PPzw9Hjx613qusrMSVK1eg0+l6rGfXrl2Ij4/HtGnTeixbWloKHx8fhIaGutV2Iuo9zgYoncvNiApGmMYfro6r7DjyNybtIxrEJAU1GRkZ+OSTT3Dp0iWUlZUhIyMDBQUFWL58OTQaDdasWYPNmzfj+PHjKCkpwRNPPAGdTmez82nSpEnIy8uzqbehoQF79+7FD3/4wy6fWVRUhFdeeQV//etf8fe//x3vvfce0tPTsWLFCtx2220udpuIeltPAYoCHaMoM6KCrdd8fRTYujDGet8TmLSPaPCQFNTU1tZi5cqViI6Oxty5c1FcXIxDhw4hKSkJALBjxw78+7//O1JSUjBr1ixotVrs27fPpo7KykoYjUaba7m5uRBCYNmyZV0+U6VSITc3F/fffz+mTJmCl156Cenp6XjzzTel9pWI+lB3AYrl/daFMV3WuyTHhiF7xXRoNbYjPe4si2HSPqLBQSGEGBT/6yLl6HIi8hxXD6e89ZTub5tMSM35HABsdlQpbnnfnZ2P3YFFd3xPeieIqN9I+f3NU7qJqFclx4YhKUYreUeSr4/CuivKIttH0SVA0mr88djdEdhx5EKPbXF1ITIRDQwMaoio19kLUFzhKEACgNzir2AwttgdtVGgI/jpvH6HiOSHQQ0RDSiOAqStC2OwYc/ZLtNR3a3fISJ5cXtLNxGRN3C0wFir8Uf2iunMU0M0CHCkhohkw9X1O0QkDwxqiEhWPLV+h4gGHk4/ERERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLQ/q7AUTkvHazwOnqOtQ2tiA00B8zooLh66Po72YREXkFBjVEA0R+uR7bDlRAb2yxXgvT+GPrwhgkx4b1Y8uIiLwDp5+IBoD8cj027DlrE9AAgMHYgg17ziK/XN9PLSMi8h4Maoi8XLtZYNuBCgg79yzXth2oQLvZXgkiosGDQQ2RlztdXddlhKYzAUBvbMHp6rq+axQRkRdiUEPk5WobHQc0rpQjIpIrBjVEXi400N+j5YiI5IpBDZGXmxEVjDCNPxxt3FagYxfUjKjgvmwWEZHXYVBD5OV8fRTYujAGALoENpb3WxfGMF8NEQ16DGqIBoDk2DBkr5gOrcZ2ikmr8Uf2iunMU0NEBIlBTXZ2NuLi4qBWq6FWq6HT6XDw4EHr/ZaWFqSmpiIkJAQBAQFISUlBTU1Nt3WuXr0aCoXC5pWcnGxTpq6uDsuXL4darUZQUBDWrFmDGzduSGk60YCXHBuGT7fMwftrZ2LnY3fg/bUz8emWOQxoiIj+SVJG4dGjRyMrKwsTJ06EEALvvvsuFi1ahM8//xxTpkxBeno6Pv74Y+zduxcajQZpaWlYsmQJTpw40W29ycnJ2L17t/W9SqWyub98+XLo9XocPnwYbW1teOKJJ7Bu3Trk5ORIaT7RgOfro4BufEh/N4OIyCsphBBuZewKDg7G9u3bsXTpUowcORI5OTlYunQpAODLL7/E5MmTUVRUhJkzZ9p9fvXq1aivr8eHH35o9/4XX3yBmJgYFBcX46677gIA5OfnY8GCBfj6668RHh7uVDsbGhqg0WhgNBqhVquld5SIiIj6nJTf3y6vqWlvb0dubi6ampqg0+lQUlKCtrY2JCYmWstMmjQJkZGRKCoq6raugoIChIaGIjo6Ghs2bMD169et94qKihAUFGQNaAAgMTERPj4+OHXqlMM6TSYTGhoabF5EREQkX5KDmrKyMgQEBEClUmH9+vXIy8tDTEwMDAYDlEolgoKCbMqPGjUKBoPBYX3Jycn4wx/+gKNHj+Lll19GYWEh5s+fj/b2dgCAwWBAaGiozTNDhgxBcHBwt/VmZmZCo9FYXxEREVK7SkRERAOI5FO6o6OjUVpaCqPRiA8++ACrVq1CYWGhyw147LHHrH+eOnUq4uLiMH78eBQUFGDu3Lku15uRkYHNmzdb3zc0NDCwISIikjHJQY1SqcSECRMAAPHx8SguLsbOnTvx6KOPorW1FfX19TajNTU1NdBqtU7XP27cOIwYMQIXL17E3LlzodVqUVtba1Pm5s2bqKur67ZelUrVZcExERERyZfbeWrMZjNMJhPi4+Ph5+eHo0ePWu9VVlbiypUr0Ol0Ttf39ddf4/r16wgL69imqtPpUF9fj5KSEmuZY8eOwWw2IyEhwd3mExERkUxICmoyMjLwySef4NKlSygrK0NGRgYKCgqwfPlyaDQarFmzBps3b8bx48dRUlKCJ554Ajqdzmbn06RJk5CXlwcAuHHjBn7605/i5MmTuHTpEo4ePYpFixZhwoQJmDdvHgBg8uTJSE5Oxtq1a3H69GmcOHECaWlpeOyxx5ze+URERETyJ2n6qba2FitXroRer4dGo0FcXBwOHTqEpKQkAMCOHTvg4+ODlJQUmEwmzJs3D6+99ppNHZWVlTAajQAAX19fnDt3Du+++y7q6+sRHh6OBx98EL/85S9tpo7ee+89pKWlYe7cudb6f/vb37rbdyIiIpIRt/PUDBTMU0NERDTw9EmeGiIiIiJvwqCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREcnCkP5uANFA0W4WOF1dh9rGFoQG+mNGVDB8fRT93SwiIvonBjVETsgv12PbgQrojS3Wa2Eaf2xdGIPk2LB+bBkREVlw+omoB/nlemzYc9YmoAEAg7EFG/acRX65vp9aRkREnTGoIepGu1lg24EKCDv3LNe2HahAu9leCSIi6ksMaoi6cbq6rssITWcCgN7YgtPVdX3XKCIisotBDVE3ahsdBzSulCMiot7DoIaoG6GB/h4tR0REvYdBDVE3ZkQFI0zjD0cbtxXo2AU1Iyq4L5tFRER2MKgh6oavjwJbF8YAQJfAxvJ+68IY5qshIvICDGqIepAcG4bsFdOh1dhOMWk1/sheMZ15aoiIvAST7xE5ITk2DEkxWmYUJiLyYgxqiJzk66OAbnxIfzeDiIgckDT9lJ2djbi4OKjVaqjVauh0Ohw8eNB6v6WlBampqQgJCUFAQABSUlJQU1PjsL62tjZs2bIFU6dOxfDhwxEeHo6VK1fi6tWrNuXGjh0LhUJh88rKypLYVSIiIpIzSUHN6NGjkZWVhZKSEpw5cwZz5szBokWLcP78eQBAeno6Dhw4gL1796KwsBBXr17FkiVLHNbX3NyMs2fP4vnnn8fZs2exb98+VFZW4uGHH+5S9sUXX4Rer7e+Nm7cKLGrREREJGcKIYRb+d2Dg4Oxfft2LF26FCNHjkROTg6WLl0KAPjyyy8xefJkFBUVYebMmU7VV1xcjBkzZuDy5cuIjIwE0DFSs2nTJmzatMnldjY0NECj0cBoNEKtVrtcDxEREfUdKb+/Xd791N7ejtzcXDQ1NUGn06GkpARtbW1ITEy0lpk0aRIiIyNRVFTkdL1GoxEKhQJBQUE217OyshASEoI777wT27dvx82bN7utx2QyoaGhweZFRERE8iV5oXBZWRl0Oh1aWloQEBCAvLw8xMTEoLS0FEqlskswMmrUKBgMBqfqbmlpwZYtW7Bs2TKbaOzHP/4xpk+fjuDgYHz22WfIyMiAXq/Hb37zG4d1ZWZmYtu2bVK7R0RERAOU5KAmOjoapaWlMBqN+OCDD7Bq1SoUFha63ZC2tjb84Ac/gBAC2dnZNvc2b95s/XNcXByUSiWeeuopZGZmQqVS2a0vIyPD5rmGhgZERES43U4iIiLyTpKDGqVSiQkTJgAA4uPjUVxcjJ07d+LRRx9Fa2sr6uvrbUZrampqoNVqu63TEtBcvnwZx44d63HOLCEhATdv3sSlS5cQHR1tt4xKpXIY8BAREZH8uJ1R2Gw2w2QyIT4+Hn5+fjh69Kj1XmVlJa5cuQKdTufweUtAc+HCBRw5cgQhIT3nASktLYWPjw9CQ0PdbT4RERHJhKSRmoyMDMyfPx+RkZFobGxETk4OCgoKcOjQIWg0GqxZswabN29GcHAw1Go1Nm7cCJ1OZ7PzadKkScjMzMTixYvR1taGpUuX4uzZs/joo4/Q3t5uXX8THBwMpVKJoqIinDp1Cg888AACAwNRVFSE9PR0rFixArfddptnfxpEREQ0YEkKampra7Fy5Uro9XpoNBrExcXh0KFDSEpKAgDs2LEDPj4+SElJgclkwrx58/Daa6/Z1FFZWQmj0QgA+Oabb/DnP/8ZAHDHHXfYlDt+/Dhmz54NlUqF3NxcvPDCCzCZTIiKikJ6errNehkiIiIit/PUDBTMU0NERDTw9EmeGiIiIiJvwqCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBDREREckCgxoiIiKSBQY1REREJAsMaoiIiEgWGNQQERGRLEgKarKzsxEXFwe1Wg21Wg2dToeDBw9a77e0tCA1NRUhISEICAhASkoKampquq1TCIFf/OIXCAsLw9ChQ5GYmIgLFy7YlKmrq8Py5cuhVqsRFBSENWvW4MaNG1KaTkRERDInKagZPXo0srKyUFJSgjNnzmDOnDlYtGgRzp8/DwBIT0/HgQMHsHfvXhQWFuLq1atYsmRJt3X++te/xm9/+1u8/vrrOHXqFIYPH4558+ahpaXFWmb58uU4f/48Dh8+jI8++giffPIJ1q1b50J3B692s0BR1XXsL/0GRVXX0W4W/d0kIiIij1IIIdz67RYcHIzt27dj6dKlGDlyJHJycrB06VIAwJdffonJkyejqKgIM2fO7PKsEALh4eH4yU9+gmeeeQYAYDQaMWrUKLzzzjt47LHH8MUXXyAmJgbFxcW46667AAD5+flYsGABvv76a4SHhzvVzoaGBmg0GhiNRqjVane6PODkl+ux7UAF9MZ/BYphGn9sXRiD5NiwfmwZERFR96T8/nZ5TU17eztyc3PR1NQEnU6HkpIStLW1ITEx0Vpm0qRJiIyMRFFRkd06qqurYTAYbJ7RaDRISEiwPlNUVISgoCBrQAMAiYmJ8PHxwalTpxy2z2QyoaGhweY1GOWX67Fhz1mbgAYADMYWbNhzFvnl+n5qGRERkWdJDmrKysoQEBAAlUqF9evXIy8vDzExMTAYDFAqlQgKCrIpP2rUKBgMBrt1Wa6PGjXK4TMGgwGhoaE294cMGYLg4GCH9QJAZmYmNBqN9RURESG1qwNeu1lg24EK2BuKs1zbdqCCU1FERCQLkoOa6OholJaW4tSpU9iwYQNWrVqFioqK3mibWzIyMmA0Gq2vr776qr+b1OdOV9d1GaHpTADQG1twurqu7xpFRETUS4ZIfUCpVGLChAkAgPj4eBQXF2Pnzp149NFH0draivr6epvRmpqaGmi1Wrt1Wa7X1NQgLCzM5pk77rjDWqa2ttbmuZs3b6Kurs5hvQCgUqmgUqmkdk9WahsdBzSulCMiIvJmbuepMZvNMJlMiI+Ph5+fH44ePWq9V1lZiStXrkCn09l9NioqClqt1uaZhoYGnDp1yvqMTqdDfX09SkpKrGWOHTsGs9mMhIQEd5sva6GB/h4tR0RE5M0kjdRkZGRg/vz5iIyMRGNjI3JyclBQUIBDhw5Bo9FgzZo12Lx5M4KDg6FWq7Fx40bodDqbnU+TJk1CZmYmFi9eDIVCgU2bNuE///M/MXHiRERFReH5559HeHg4HnnkEQDA5MmTkZycjLVr1+L1119HW1sb0tLS8Nhjjzm982mwmhEVjDCNPwzGFrvrahQAtBp/zIgK7uumEREReZykoKa2thYrV66EXq+HRqNBXFwcDh06hKSkJADAjh074OPjg5SUFJhMJsybNw+vvfaaTR2VlZUwGo3W9z/72c/Q1NSEdevWob6+Hvfeey/y8/Ph7/+v0YP33nsPaWlpmDt3rrX+3/72t+70e1Dw9VFg68IYbNhzFgrAJrBR/PO/WxfGwNdHYedpIiKigcXtPDUDBfPUME8NERENPFJ+f0teKEy22s0Cp6vrUNvYgtDAjqkcbxv5SI4NQ1KM1uvbSURE5A4GNW4YSCMgvj4K6MaH9HcziIiIeg1P6XYRM/USERF5FwY1LmCmXiIiIu/DoMYFzNRLRETkfRjUuICZeomIiLwPgxoXMFMvERGR92FQ4wJLpl5HG6IV6NgFxUy9REREfYdBjQssmXoBdAlsmKmXiIiofzCocVFybBiyV0yHVmM7xaTV+CN7xXSvy1NDREQkd0y+5wZm6iUiIvIeDGrcxEy9RERE3oHTT0RERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFnhMgoe0mwXPgCIiIupHDGo8IL9cj20HKqA3tlivhWn8sXVhDE/rJiIi6iOcfnJTfrkeG/actQloAMBgbMGGPWeRX67vp5YRERENLgxq3NBuFth2oALCzj3LtW0HKtButleCiIiIPIlBjRtOV9d1GaHpTADQG1twurqu7xpFREQ0SDGocUNto+OAxpVyRERE5DoGNW4IDfT3aDkiIiJyHYMaN8yICkaYxh+ONm4r0LELakZUcF82i4iIaFCSFNRkZmbi7rvvRmBgIEJDQ/HII4+gsrLSpkxVVRUWL16MkSNHQq1W4wc/+AFqamq6rXfs2LFQKBRdXqmpqdYys2fP7nJ//fr1Uprvcb4+CmxdGAMAXQIby/utC2OYr4aIiKgPSApqCgsLkZqaipMnT+Lw4cNoa2vDgw8+iKamJgBAU1MTHnzwQSgUChw7dgwnTpxAa2srFi5cCLPZ7LDe4uJi6PV66+vw4cMAgO9///s25dauXWtT7te//rXU/npccmwYsldMh1ZjO8Wk1fgje8V05qkhIiLqI5KS7+Xn59u8f+eddxAaGoqSkhLMmjULJ06cwKVLl/D5559DrVYDAN59913cdtttOHbsGBITE+3WO3LkSJv3WVlZGD9+PO6//36b68OGDYNWq5XS5D6RHBuGpBgtMwoTERH1I7fW1BiNRgBAcHDHmhGTyQSFQgGVSmUt4+/vDx8fH3z66adO1dna2oo9e/bgySefhEJhGxS89957GDFiBGJjY5GRkYHm5mZ3mu9Rvj4K6MaHYNEd34NufAgDGiIioj7m8jEJZrMZmzZtwj333IPY2FgAwMyZMzF8+HBs2bIFv/rVryCEwLPPPov29nbo9c5l1v3www9RX1+P1atX21x//PHHMWbMGISHh+PcuXPYsmULKisrsW/fPrv1mEwmmEwm6/uGhgbXOkpEREQDgssjNampqSgvL0dubq712siRI7F3714cOHAAAQEB0Gg0qK+vx/Tp0+Hj49xH7dq1C/Pnz0d4eLjN9XXr1mHevHmYOnUqli9fjj/84Q/Iy8tDVVWV3XoyMzOh0Wisr4iICFe7SkRERAOASyM1aWlp+Oijj/DJJ59g9OjRNvcefPBBVFVV4dq1axgyZAiCgoKg1Woxbty4Huu9fPkyjhw54nD0pbOEhAQAwMWLFzF+/Pgu9zMyMrB582br+4aGBgY2REREMiYpqBFCYOPGjcjLy0NBQQGioqIclh0xYgQA4NixY6itrcXDDz/cY/27d+9GaGgoHnrooR7LlpaWAgDCwuzvLlKpVDZre4iIiEjeJAU1qampyMnJwf79+xEYGAiDwQAA0Gg0GDp0KICOwGTy5MkYOXIkioqK8PTTTyM9PR3R0dHWeubOnYvFixcjLS3Nes1sNmP37t1YtWoVhgyxbVZVVRVycnKwYMEChISE4Ny5c0hPT8esWbMQFxfncueJiIhIPiQFNdnZ2QA6EuF1tnv3buvC3srKSmRkZKCurg5jx47Ff/zHfyA9Pd2mvGV6qrMjR47gypUrePLJJ7t8rlKpxJEjR/DKK6+gqakJERERSElJwc9//nMpzSciIiIZUwghRH83oi80NDRAo9HAaDRac+gQERGRd5Py+5tnPxEREZEsMKghIiIiWXA5+d5AY5llYxI+IiKigcPye9uZ1TKDJqhpbGwEAOaqISIiGoAaGxuh0Wi6LTNoFgqbzWZcvXoVgYGBXc6UGugsiQW/+uqrQbEIejD1dzD1FWB/5Www9RUYXP3t7b4KIdDY2Ijw8PAeTycYNCM1Pj4+XbIfy41arZb9P57OBlN/B1NfAfZXzgZTX4HB1d/e7GtPIzQWXChMREREssCghoiIiGSBQY0MqFQqbN26ddCcdTWY+juY+gqwv3I2mPoKDK7+elNfB81CYSIiIpI3jtQQERGRLDCoISIiIllgUENERESywKCGiIiIZIFBTT/Lzs5GXFycNWmRTqfDwYMHrfdbWlqQmpqKkJAQBAQEICUlBTU1Nd3WKYTAL37xC4SFhWHo0KFITEzEhQsXbMrU1dVh+fLlUKvVCAoKwpo1a3Djxo1e6aOFp/va1taGLVu2YOrUqRg+fDjCw8OxcuVKXL161abc2LFjoVAobF5ZWVm91k+L3vhuV69e3aUvycnJNmXk8N0C6NJPy2v79u3WMt763b755puYPXs21Go1FAoF6uvrnar397//PcaOHQt/f38kJCTg9OnTNvdd+Tm6qzf6mpmZibvvvhuBgYEIDQ3FI488gsrKSpsys2fP7vLdrl+/3tPd66I3+vvCCy906cukSZNsyvTHdwv0Tn/t/btUKBRITU21lum171dQv/rzn/8sPv74Y/G3v/1NVFZWiueee074+fmJ8vJyIYQQ69evFxEREeLo0aPizJkzYubMmeLf/u3fuq0zKytLaDQa8eGHH4q//vWv4uGHHxZRUVHiu+++s5ZJTk4W06ZNEydPnhR/+ctfxIQJE8SyZcsGVF/r6+tFYmKi+OMf/yi+/PJLUVRUJGbMmCHi4+Ntyo0ZM0a8+OKLQq/XW183btzo1b4K0Tvf7apVq0RycrJNX+rq6mzKyOG7FULY9FGv14u3335bKBQKUVVVZS3jrd/tjh07RGZmpsjMzBQAxLfffttjnbm5uUKpVIq3335bnD9/Xqxdu1YEBQWJmpoaaxlXfo7u6o2+zps3T+zevVuUl5eL0tJSsWDBAhEZGWnz3d1///1i7dq1Nt+t0WjsrW5a9UZ/t27dKqZMmWLTl3/84x82ZfrjuxWid/pbW1tr09fDhw8LAOL48ePWMr31/TKo8UK33Xab+K//+i9RX18v/Pz8xN69e633vvjiCwFAFBUV2X3WbDYLrVYrtm/fbr1WX18vVCqVeP/994UQQlRUVAgAori42Frm4MGDQqFQiG+++aaXemWfO3215/Tp0wKAuHz5svXamDFjxI4dOzzZbJe5299Vq1aJRYsWObwv5+920aJFYs6cOTbXvPG77ez48eNO/yKYMWOGSE1Ntb5vb28X4eHhIjMzUwghPPZz9AR3+3qr2tpaAUAUFhZar91///3i6aefdrOlnuFuf7du3SqmTZvm8L43fbdCeP77ffrpp8X48eOF2Wy2Xuut75fTT16kvb0dubm5aGpqgk6nQ0lJCdra2pCYmGgtM2nSJERGRqKoqMhuHdXV1TAYDDbPaDQaJCQkWJ8pKipCUFAQ7rrrLmuZxMRE+Pj44NSpU73UO1ue6Ks9RqMRCoUCQUFBNtezsrIQEhKCO++8E9u3b8fNmzc91RWneLK/BQUFCA0NRXR0NDZs2IDr169b78n1u62pqcHHH3+MNWvWdLnnbd+tK1pbW1FSUmLzM/Lx8UFiYqL1Z+SpfyPu8ERf7TEajQCA4OBgm+vvvfceRowYgdjYWGRkZKC5udljn+kMT/b3woULCA8Px7hx47B8+XJcuXLFes8bvlugd77f1tZW7NmzB08++WSXw6R74/sdNAdaerOysjLodDq0tLQgICAAeXl5iImJQWlpKZRKZZdf0KNGjYLBYLBbl+X6qFGjHD5jMBgQGhpqc3/IkCEIDg52WK+neLKvt2ppacGWLVuwbNkym0PVfvzjH2P69OkIDg7GZ599hoyMDOj1evzmN7/xZNfs8nR/k5OTsWTJEkRFRaGqqgrPPfcc5s+fj6KiIvj6+sr2u3333XcRGBiIJUuW2Fz3xu/WFdeuXUN7e7vdf7dffvklgI5/t+7+HF3lyb7eymw2Y9OmTbjnnnsQGxtrvf74449jzJgxCA8Px7lz57BlyxZUVlZi3759Hvnc7ni6vwkJCXjnnXcQHR0NvV6Pbdu24b777kN5eTkCAwP79bsFevf7/fDDD1FfX4/Vq1fbXO+t75dBjReIjo5GaWkpjEYjPvjgA6xatQqFhYX93axe0Vt9bWtrww9+8AMIIZCdnW1zb/PmzdY/x8XFQalU4qmnnkJmZmavp/X2dH8fe+wx65+nTp2KuLg4jB8/HgUFBZg7d64nmuyy3vx7/Pbbb2P58uXw9/e3ue6N362nfhl4k97sa2pqKsrLy/Hpp5/aXF+3bp31z1OnTkVYWBjmzp2LqqoqjB8/3u3P7Y6n+zt//nzrn+Pi4pCQkIAxY8bgT3/6k93Rx77Wm9/vrl27MH/+fISHh9tc763vl9NPXkCpVGLChAmIj49HZmYmpk2bhp07d0Kr1aK1tbXLavOamhpotVq7dVmu37pqvvMzWq0WtbW1Nvdv3ryJuro6h/V6iif7amEJaC5fvozDhw/bjNLYk5CQgJs3b+LSpUtu9qZnvdHfzsaNG4cRI0bg4sWLAOT33QLAX/7yF1RWVuKHP/xhj2W94bt1xYgRI+Dr69vjv1tP/J1xhSf72llaWho++ugjHD9+HKNHj+62bEJCAgBY/673pt7qr0VQUBBuv/12m3+3/fXdAr3X38uXL+PIkSNO/9sF3P9+GdR4IbPZDJPJhPj4ePj5+eHo0aPWe5WVlbhy5YrD+c6oqChotVqbZxoaGnDq1CnrMzqdDvX19SgpKbGWOXbsGMxms/UvVl9xp6/AvwKaCxcu4MiRIwgJCenxM0tLS+Hj49NlmqYvuNvfW3399de4fv06wsLCAMjru7XYtWsX4uPjMW3atB7LesN36wqlUon4+Hibn5HZbMbRo0etPyNP/Z3xBHf6CnSknUhLS0NeXh6OHTuGqKioHp8pLS0FAOvf9b7kbn9vdePGDVRVVVn74k3fLeC5/u7evRuhoaF46KGHeizrse/X40uPSZJnn31WFBYWiurqanHu3Dnx7LPPCoVCIf73f/9XCNGxzS8yMlIcO3ZMnDlzRuh0OqHT6WzqiI6OFvv27bO+z8rKEkFBQWL//v3i3LlzYtGiRXa3dN95553i1KlT4tNPPxUTJ07s9W2/nu5ra2urePjhh8Xo0aNFaWmpzdZAk8kkhBDis88+Ezt27BClpaWiqqpK7NmzR4wcOVKsXLmyV/vaG/1tbGwUzzzzjCgqKhLV1dXiyJEjYvr06WLixImipaXF+owcvlsLo9Eohg0bJrKzs7t8pjd/t3q9Xnz++efirbfeEgDEJ598Ij7//HNx/fp1ax1z5swRv/vd76zvc3NzhUqlEu+8846oqKgQ69atE0FBQcJgMFjLOPNzHAh93bBhg9BoNKKgoMDm321zc7MQQoiLFy+KF198UZw5c0ZUV1eL/fv3i3HjxolZs2b1al97q78/+clPREFBgaiurhYnTpwQiYmJYsSIEaK2ttZapj++297qrxAdu/ciIyPFli1bunxmb36/DGr62ZNPPinGjBkjlEqlGDlypJg7d671L5MQQnz33XfiRz/6kbjtttvEsGHDxOLFi4Ver7epA4DYvXu39b3ZbBbPP/+8GDVqlFCpVGLu3LmisrLS5pnr16+LZcuWiYCAAKFWq8UTTzwhGhsbB1Rfq6urBQC7L0s+hJKSEpGQkCA0Go3w9/cXkydPFr/61a9sgoCB0t/m5mbx4IMPipEjRwo/Pz8xZswYsXbtWptfekLI47u1eOONN8TQoUNFfX19l8/05u9269atdv9edu7fmDFjxNatW23q/d3vficiIyOFUqkUM2bMECdPnrS578zP0dN6o6+O/t1anrly5YqYNWuWCA4OFiqVSkyYMEH89Kc/7ZM8Nb3R30cffVSEhYUJpVIpvve974lHH31UXLx40eZz++O7FaL3/i4fOnRIAOjyu0eI3v1+FUII4d5YDxEREVH/45oaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSwwqCEiIiJZYFBDREREssCghoiIiGSBQQ0RERHJAoMaIiIikgUGNURERCQLDGqIiIhIFhjUEBERkSz8f5TsShC908s9AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure()\n", + "# plt.hist(mock[\"mu_TFR\"] - mock[\"mu_calibration\"][1])\n", + "plt.scatter(mock[\"mu_true\"], mock[\"mu_calibration\"][0])\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 127, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([125])" + ] + }, + "execution_count": 127, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.sum(np.isfinite(mock[\"mu_calibration\"]), axis=1)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhOUlEQVR4nO3dfXBU5f2/8feGkBADuyEREiKJZJQRrAgKilEHETI8DoKmtVTaWmWIVUAhM2Li8FD7VYMMRQaLoFYjVKhKp6BCjbVRodYYIIgoRR4qSDRssIPZhSAhkvv3h8P+XIgJkLOce5PrNXNmmnNOls/tQnL17JPHGGMEAABgkRi3BwAAADgVgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOrFuD3AuGhoaVFVVpU6dOsnj8bg9DgAAOAPGGB0+fFjp6emKiWn6GklUBkpVVZUyMjLcHgMAAJyDyspKde/evclzojJQOnXqJOn7BXq9XpenAQAAZyIYDCojIyP0e7wpURkoJx/W8Xq9BAoAAFHmTJ6ewZNkAQCAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgnVi3BwDQ+vQoWNei7983d7RDkwCIVlxBAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1uHTjAFYh09DBsAVFAAAYJ2zDpQNGzZozJgxSk9Pl8fj0Zo1a0LH6uvr9dBDD6lPnz5KTExUenq6fv3rX6uqqirsNg4dOqQJEybI6/UqKSlJEydO1JEjR1q8GAAA0DqcdaDU1taqb9++Wrx48WnHjh49qi1btmjWrFnasmWL/va3v2nnzp265ZZbws6bMGGCtm/frrfffltr167Vhg0blJeXd+6rAAAArYrHGGPO+Zs9Hq1evVrjxo370XM2bdqka6+9Vl988YUyMzO1Y8cOXX755dq0aZMGDBggSSopKdGoUaP05ZdfKj09vdk/NxgMyufzKRAIyOv1nuv4ACKkpc8haSmegwLY6Wx+f0f8OSiBQEAej0dJSUmSpLKyMiUlJYXiRJJycnIUExOj8vLyRm+jrq5OwWAwbAMAAK1XRAPl2LFjeuihh/SLX/wiVEp+v19du3YNOy82NlbJycny+/2N3k5RUZF8Pl9oy8jIiOTYAADAZRELlPr6et1+++0yxmjJkiUtuq3CwkIFAoHQVllZ6dCUAADARhF5H5STcfLFF1/onXfeCXucKS0tTQcPHgw7/7vvvtOhQ4eUlpbW6O3Fx8crPj4+EqMCAAALOX4F5WSc7N69W//85z+VkpISdjw7O1s1NTWqqKgI7XvnnXfU0NCggQMHOj0OAACIQmd9BeXIkSPas2dP6Ou9e/dq69atSk5OVrdu3fTTn/5UW7Zs0dq1a3XixInQ80qSk5MVFxen3r17a8SIEZo0aZKWLl2q+vp6TZkyRePHjz+jV/AAAIDW76wDZfPmzbr55ptDX+fn50uS7rzzTv3ud7/T66+/Lknq169f2Pe9++67Gjx4sCRpxYoVmjJlioYOHaqYmBjl5uZq0aJF57gEAADQ2px1oAwePFhNvXXKmbytSnJyslauXHm2fzQAAGgj+CweAABgHQIFAABYh0ABAADWIVAAAIB1IvJGbQCim9sf9gcAXEEBAADWIVAAAIB1CBQAAGAdAgUAAFiHJ8kCrRBPcgUQ7biCAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsE6s2wMAOF2PgnVujwAAruIKCgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALDOWQfKhg0bNGbMGKWnp8vj8WjNmjVhx40xmj17trp166aEhATl5ORo9+7dYeccOnRIEyZMkNfrVVJSkiZOnKgjR460aCEAAKD1OOtAqa2tVd++fbV48eJGj8+bN0+LFi3S0qVLVV5ersTERA0fPlzHjh0LnTNhwgRt375db7/9ttauXasNGzYoLy/v3FcBAABalbP+LJ6RI0dq5MiRjR4zxmjhwoWaOXOmxo4dK0lavny5UlNTtWbNGo0fP147duxQSUmJNm3apAEDBkiSnnrqKY0aNUrz589Xenp6C5YDAABaA0efg7J37175/X7l5OSE9vl8Pg0cOFBlZWWSpLKyMiUlJYXiRJJycnIUExOj8vLyRm+3rq5OwWAwbAMAAK2Xo4Hi9/slSampqWH7U1NTQ8f8fr+6du0adjw2NlbJycmhc05VVFQkn88X2jIyMpwcGwAAWCYqXsVTWFioQCAQ2iorK90eCQAARJCjgZKWliZJqq6uDttfXV0dOpaWlqaDBw+GHf/uu+906NCh0Dmnio+Pl9frDdsAAEDr5WigZGVlKS0tTaWlpaF9wWBQ5eXlys7OliRlZ2erpqZGFRUVoXPeeecdNTQ0aODAgU6OAwAAotRZv4rnyJEj2rNnT+jrvXv3auvWrUpOTlZmZqamTZumRx99VD179lRWVpZmzZql9PR0jRs3TpLUu3dvjRgxQpMmTdLSpUtVX1+vKVOmaPz48byCBwAASDqHQNm8ebNuvvnm0Nf5+fmSpDvvvFMvvviiZsyYodraWuXl5ammpkY33nijSkpK1KFDh9D3rFixQlOmTNHQoUMVExOj3NxcLVq0yIHlAACA1sBjjDFuD3G2gsGgfD6fAoEAz0dBq9SjYJ3bI0S1fXNHuz0CgEacze/vqHgVDwAAaFsIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWiXV7AKA16lGwzu0R0AItvf/2zR3t0CRA28UVFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADW4X1QALQ6vA8NEP24ggIAAKxDoAAAAOsQKAAAwDqOB8qJEyc0a9YsZWVlKSEhQZdccon+7//+T8aY0DnGGM2ePVvdunVTQkKCcnJytHv3bqdHAQAAUcrxQHniiSe0ZMkS/fGPf9SOHTv0xBNPaN68eXrqqadC58ybN0+LFi3S0qVLVV5ersTERA0fPlzHjh1zehwAABCFHH8VzwcffKCxY8dq9OjvP82zR48e+stf/qKNGzdK+v7qycKFCzVz5kyNHTtWkrR8+XKlpqZqzZo1Gj9+vNMjAQCAKOP4FZTrr79epaWl2rVrlyTp448/1vvvv6+RI0dKkvbu3Su/36+cnJzQ9/h8Pg0cOFBlZWVOjwMAAKKQ41dQCgoKFAwG1atXL7Vr104nTpzQY489pgkTJkiS/H6/JCk1NTXs+1JTU0PHTlVXV6e6urrQ18Fg0OmxAQCARRy/gvLqq69qxYoVWrlypbZs2aJly5Zp/vz5WrZs2TnfZlFRkXw+X2jLyMhwcGIAAGAbxwPlwQcfVEFBgcaPH68+ffroV7/6laZPn66ioiJJUlpamiSpuro67Puqq6tDx05VWFioQCAQ2iorK50eGwAAWMTxQDl69KhiYsJvtl27dmpoaJAkZWVlKS0tTaWlpaHjwWBQ5eXlys7ObvQ24+Pj5fV6wzYAANB6Of4clDFjxuixxx5TZmamfvKTn+ijjz7SggULdPfdd0uSPB6Ppk2bpkcffVQ9e/ZUVlaWZs2apfT0dI0bN87pcQAAQBRyPFCeeuopzZo1S/fdd58OHjyo9PR03XPPPZo9e3bonBkzZqi2tlZ5eXmqqanRjTfeqJKSEnXo0MHpcQAAQBTymB++xWuUCAaD8vl8CgQCPNwDK/Fpum3bvrmj3R4BsNLZ/P7ms3gAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1Yt0eALANn0QMAO7jCgoAALAOgQIAAKzDQzwA4DAnHibcN3e0A5MA0YsrKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsE6s2wMATutRsM7tEQAALcQVFAAAYJ2IBMpXX32lX/7yl0pJSVFCQoL69OmjzZs3h44bYzR79mx169ZNCQkJysnJ0e7duyMxCgAAiEKOB8o333yjG264Qe3bt9ebb76p//znP/rDH/6gzp07h86ZN2+eFi1apKVLl6q8vFyJiYkaPny4jh075vQ4AAAgCjn+HJQnnnhCGRkZKi4uDu3LysoK/W9jjBYuXKiZM2dq7NixkqTly5crNTVVa9as0fjx450eCQAARBnHr6C8/vrrGjBggH72s5+pa9euuuqqq/Tcc8+Fju/du1d+v185OTmhfT6fTwMHDlRZWVmjt1lXV6dgMBi2AQCA1svxQPn888+1ZMkS9ezZU2+99Zbuvfde3X///Vq2bJkkye/3S5JSU1PDvi81NTV07FRFRUXy+XyhLSMjw+mxAQCARRwPlIaGBl199dV6/PHHddVVVykvL0+TJk3S0qVLz/k2CwsLFQgEQltlZaWDEwMAANs4HijdunXT5ZdfHravd+/e2r9/vyQpLS1NklRdXR12TnV1dejYqeLj4+X1esM2AADQejkeKDfccIN27twZtm/Xrl26+OKLJX3/hNm0tDSVlpaGjgeDQZWXlys7O9vpcQAAQBRy/FU806dP1/XXX6/HH39ct99+uzZu3Khnn31Wzz77rCTJ4/Fo2rRpevTRR9WzZ09lZWVp1qxZSk9P17hx45weBwAARCHHA+Waa67R6tWrVVhYqN///vfKysrSwoULNWHChNA5M2bMUG1trfLy8lRTU6Mbb7xRJSUl6tChg9PjAACAKOQxxhi3hzhbwWBQPp9PgUCA56PgNHwWD1qDfXNHuz0C4Liz+f3NZ/EAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArBPxQJk7d648Ho+mTZsW2nfs2DFNnjxZKSkp6tixo3Jzc1VdXR3pUQAAQJSIaKBs2rRJzzzzjK688sqw/dOnT9cbb7yhVatWaf369aqqqtJtt90WyVEAAEAUiVigHDlyRBMmTNBzzz2nzp07h/YHAgE9//zzWrBggYYMGaL+/furuLhYH3zwgT788MNIjQMAAKJIxAJl8uTJGj16tHJycsL2V1RUqL6+Pmx/r169lJmZqbKyskZvq66uTsFgMGwDAACtV2wkbvTll1/Wli1btGnTptOO+f1+xcXFKSkpKWx/amqq/H5/o7dXVFSkRx55JBKjAoCVehSsa9H375s72qFJAHc4fgWlsrJSDzzwgFasWKEOHTo4cpuFhYUKBAKhrbKy0pHbBQAAdnI8UCoqKnTw4EFdffXVio2NVWxsrNavX69FixYpNjZWqampOn78uGpqasK+r7q6WmlpaY3eZnx8vLxeb9gGAABaL8cf4hk6dKg++eSTsH133XWXevXqpYceekgZGRlq3769SktLlZubK0nauXOn9u/fr+zsbKfHQRRq6aVtAED0czxQOnXqpCuuuCJsX2JiolJSUkL7J06cqPz8fCUnJ8vr9Wrq1KnKzs7Wdddd5/Q4AAAgCkXkSbLNefLJJxUTE6Pc3FzV1dVp+PDhevrpp90YBQAAWMhjjDFuD3G2gsGgfD6fAoEAz0dphXiIB2g5XsUDG53N728+iwcAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHVi3R4ArU+PgnVujwAAiHJcQQEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANbh04wBoBVq6aeK75s72qFJgHPDFRQAAGAdAgUAAFiHQAEAANZxPFCKiop0zTXXqFOnTuratavGjRunnTt3hp1z7NgxTZ48WSkpKerYsaNyc3NVXV3t9CgAACBKOR4o69ev1+TJk/Xhhx/q7bffVn19vYYNG6ba2trQOdOnT9cbb7yhVatWaf369aqqqtJtt93m9CgAACBKOf4qnpKSkrCvX3zxRXXt2lUVFRUaNGiQAoGAnn/+ea1cuVJDhgyRJBUXF6t379768MMPdd111zk9EgAAiDIRfw5KIBCQJCUnJ0uSKioqVF9fr5ycnNA5vXr1UmZmpsrKyhq9jbq6OgWDwbANAAC0XhENlIaGBk2bNk033HCDrrjiCkmS3+9XXFyckpKSws5NTU2V3+9v9HaKiork8/lCW0ZGRiTHBgAALotooEyePFmffvqpXn755RbdTmFhoQKBQGirrKx0aEIAAGCjiL2T7JQpU7R27Vpt2LBB3bt3D+1PS0vT8ePHVVNTE3YVpbq6WmlpaY3eVnx8vOLj4yM1KgAAsIzjgWKM0dSpU7V69Wq99957ysrKCjvev39/tW/fXqWlpcrNzZUk7dy5U/v371d2drbT4+ActPQtsgEAaCnHA2Xy5MlauXKlXnvtNXXq1Cn0vBKfz6eEhAT5fD5NnDhR+fn5Sk5Oltfr1dSpU5Wdnc0reAAAgKQIBMqSJUskSYMHDw7bX1xcrN/85jeSpCeffFIxMTHKzc1VXV2dhg8frqefftrpUQAAQJSKyEM8zenQoYMWL16sxYsXO/3HAwCAVoDP4gEAANYhUAAAgHUIFAAAYB0CBQAAWCdib9QGAIheLX0/pH1zRzs0CdoqrqAAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA7vgwIAcBzvo4KW4goKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDp9m3Aq19FNEAQBwG1dQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB3e6h4A0Oq09CM/9s0d7dAkOFdcQQEAANYhUAAAgHV4iKcRXBoEAMBdXEEBAADWIVAAAIB1eIjHQi19iAkAoh0/B1vGif9+bj9dwdUrKIsXL1aPHj3UoUMHDRw4UBs3bnRzHAAAYAnXAuWVV15Rfn6+5syZoy1btqhv374aPny4Dh486NZIAADAEq49xLNgwQJNmjRJd911lyRp6dKlWrdunV544QUVFBS4NRYAAC3GQ1Qt50qgHD9+XBUVFSosLAzti4mJUU5OjsrKyk47v66uTnV1daGvA4GAJCkYDEZkvoa6oy36/pbO1dI/HwDQMvwcj8zv2JO3aYxp9lxXAuV///ufTpw4odTU1LD9qamp+uyzz047v6ioSI888shp+zMyMiI2Y0v4Fro9AQCgJfg5Htn/BocPH5bP52vynKh4FU9hYaHy8/NDXzc0NOjQoUNKSUmRx+NxZaZgMKiMjAxVVlbK6/W6MoMb2uq6JdbeFtfeVtcttd21t9V1S+dn7cYYHT58WOnp6c2e60qgXHjhhWrXrp2qq6vD9ldXVystLe208+Pj4xUfHx+2LykpKZIjnjGv19vm/hJLbXfdEmtvi2tvq+uW2u7a2+q6pcivvbkrJye58iqeuLg49e/fX6WlpaF9DQ0NKi0tVXZ2thsjAQAAi7j2EE9+fr7uvPNODRgwQNdee60WLlyo2tra0Kt6AABA2+VaoPz85z/X119/rdmzZ8vv96tfv34qKSk57YmztoqPj9ecOXNOe+iptWur65ZYe1tce1tdt9R2195W1y3Zt3aPOZPX+gAAAJxHfFggAACwDoECAACsQ6AAAADrECgAAMA6BEoTlixZoiuvvDL0pjXZ2dl68803Q8ePHTumyZMnKyUlRR07dlRubu5pbz4XrZpb+7PPPqvBgwfL6/XK4/GopqbGvWEd1NS6Dx06pKlTp+qyyy5TQkKCMjMzdf/994c+GyraNXef33PPPbrkkkuUkJCgLl26aOzYsY1+NEW0aW7dJxljNHLkSHk8Hq1Zs+b8DxoBza198ODB8ng8Ydtvf/tbFyd2zpnc72VlZRoyZIgSExPl9Xo1aNAgffvtty5N7Iym1r1v377T7u+T26pVq877rARKE7p37665c+eqoqJCmzdv1pAhQzR27Fht375dkjR9+nS98cYbWrVqldavX6+qqirddtttLk/tjObWfvToUY0YMUIPP/ywy5M6q6l1V1VVqaqqSvPnz9enn36qF198USUlJZo4caLbYzuiufu8f//+Ki4u1o4dO/TWW2/JGKNhw4bpxIkTLk/eMs2t+6SFCxe69tEakXIma580aZIOHDgQ2ubNm+fixM5pbu1lZWUaMWKEhg0bpo0bN2rTpk2aMmWKYmKi+9dmU+vOyMgIu68PHDigRx55RB07dtTIkSPP/7AGZ6Vz587mT3/6k6mpqTHt27c3q1atCh3bsWOHkWTKyspcnDByTq79h959910jyXzzzTfuDHUeNLbuk1599VUTFxdn6uvrz/NU50dTa//444+NJLNnz57zPFXknbrujz76yFx00UXmwIEDRpJZvXq1e8NF2A/XftNNN5kHHnjA3YHOox+ufeDAgWbmzJkuT3R+NPXvvF+/fubuu+8+zxN9L7pT8Dw6ceKEXn75ZdXW1io7O1sVFRWqr69XTk5O6JxevXopMzNTZWVlLk7qvFPX3lacyboDgYC8Xq9iY6PiczfPWHNrr62tVXFxsbKysqz9VPFz0di6jx49qjvuuEOLFy9u9LPCWosfu89XrFihCy+8UFdccYUKCwt19OhRF6eMjFPXfvDgQZWXl6tr1666/vrrlZqaqptuuknvv/++26M6qrl/5xUVFdq6dat7V4ldyaIosm3bNpOYmGjatWtnfD6fWbdunTHGmBUrVpi4uLjTzr/mmmvMjBkzzveYEfFja/+h1ngF5UzWbYwxX3/9tcnMzDQPP/zweZ4wcppb++LFi01iYqKRZC677LJWc/WkqXXn5eWZiRMnhr5WK7uC0tTan3nmGVNSUmK2bdtmXnrpJXPRRReZW2+91cVpnfVjay8rKzOSTHJysnnhhRfMli1bzLRp00xcXJzZtWuXy1O33Jn+jLv33ntN7969z/N0/x+B0oy6ujqze/dus3nzZlNQUGAuvPBCs3379jYRKD+29h9qjYFyJusOBALm2muvNSNGjDDHjx93aVLnNbf2mpoas2vXLrN+/XozZswYc/XVV5tvv/3WxYmd8WPrfu2118yll15qDh8+HDq3tQXKmfx9P6m0tLRVPaz3Y2v/97//bSSZwsLCsPP79OljCgoKXJrWOWdynx89etT4fD4zf/58l6YkUM7a0KFDTV5eXugf6qm/mDMzM82CBQvcGS7CTq79h1pjoJzq1HUHg0GTnZ1thg4d2ip+OTelsfv8pLq6OnPBBReYlStXnuepIu/kuh944AHj8XhMu3btQpskExMTY2666Sa3x4yIpu7zI0eOGEmmpKTkPE91fpxc++eff24kmT//+c9hx2+//XZzxx13uDRd5DR2ny9fvty0b9/eHDx40KWpeA7KWWtoaFBdXZ369++v9u3bq7S0NHRs586d2r9/f6t9nsbJtbc1P1x3MBjUsGHDFBcXp9dff10dOnRwebrIauo+N9//H5xW+Xfi5LoLCgq0bds2bd26NbRJ0pNPPqni4mJ3h4yQpu7zk+vv1q3beZzo/Dm59h49eig9PV07d+4MO75r1y5dfPHFLk0XOY3d588//7xuueUWdenSxaWpXPw042hQWFiokSNHKjMzU4cPH9bKlSv13nvv6a233pLP59PEiROVn5+v5ORkeb1eTZ06VdnZ2bruuuvcHr3Fmlq7JPn9fvn9fu3Zs0eS9Mknn6hTp07KzMxUcnKym6O3SFPrPhknR48e1UsvvaRgMKhgMChJ6tKli9q1a+fy9C3T1No///xzvfLKKxo2bJi6dOmiL7/8UnPnzlVCQoJGjRrl9ugt0tS609LSGn1ibGZmprKyslyY1llNrf2///2vVq5cqVGjRiklJUXbtm3T9OnTNWjQIF155ZVuj95iTa3d4/HowQcf1Jw5c9S3b1/169dPy5Yt02effaa//vWvbo/eIs39bJekPXv2aMOGDfr73//u4qTiSbJNufvuu83FF19s4uLiTJcuXczQoUPNP/7xj9Dxb7/91tx3332mc+fO5oILLjC33nqrOXDggIsTO6e5tc+ZM8dIOm0rLi52b2gHNLXukw9nNbbt3bvX3cEd0NTav/rqKzNy5EjTtWtX0759e9O9e3dzxx13mM8++8zlqVuuub/rp1Ireg5KU2vfv3+/GTRokElOTjbx8fHm0ksvNQ8++KAJBAIuT+2MM7nfi4qKTPfu3c0FF1xgsrOzzb/+9S+XpnXOmay7sLDQZGRkmBMnTrg05fc8xhjjThoBAAA0juegAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArPP/AH77U2J8Ne9jAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "m = np.isfinite(mock[\"mu_calibration\"])\n", + "\n", + "plt.figure()\n", + "plt.hist(mock[\"mu_calibration\"][m], bins=\"auto\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "31.757872992956177" + ] + }, + "execution_count": 101, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.percentile(mock[\"mu_true\"], 10)" + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhHElEQVR4nO3df1BVdf7H8ReIIKL3EiQgCcmUk7ZZlhZRjbnKqOiYFrutmzvrlqNtqaXMZNCorfu1MMctRzOptkg32cqdzUo3Whdb3XYJFTPLNdRNkw0vtmPcq5hI8vn+0Xinq4Q/OHA+F5+PmTuznHO4vj9di+d+7g8ijDFGAAAAFol0ewAAAIDTESgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArBPl9gAXoqmpSTU1NerevbsiIiLcHgcAAJwDY4yOHDmi1NRURUa2vEcSloFSU1OjtLQ0t8cAAAAXoLq6Wr169WrxmrAMlO7du0v6boEej8flaQAAwLkIBAJKS0sL/hxvSVgGyqmndTweD4ECAECYOZeXZ/AiWQAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWCfK7QEA2Kd3/rpWff/+BaMdmgTAxYodFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1+KA2ANbhg+IAsIMCAACsQ6AAAADrECgAAMA65x0omzZt0pgxY5SamqqIiAitWbMmeK6xsVGPPvqo+vfvr7i4OKWmpuqXv/ylampqQu7j8OHDmjBhgjwej+Lj4zVp0iQdPXq01YsBAAAdw3kHSn19va677jotW7bsjHPHjh3Ttm3bNGfOHG3btk1//vOfVVVVpTvuuCPkugkTJmjnzp1av3691q5dq02bNmnKlCkXvgoAANChnPe7eHJycpSTk9PsOa/Xq/Xr14cce/bZZ3XTTTfpwIEDSk9P165du1RaWqotW7Zo0KBBkqSlS5dq1KhRWrRokVJTUy9gGQAAoCNp89eg+P1+RUREKD4+XpJUXl6u+Pj4YJxIUnZ2tiIjI1VRUdHsfTQ0NCgQCITcAABAx9Wmn4Ny/PhxPfroo/r5z38uj8cjSfL5fEpKSgodIipKCQkJ8vl8zd5PYWGh5s2b15ajAnBQaz/HBADabAelsbFRd999t4wxWr58eavuq6CgQH6/P3irrq52aEoAAGCjNtlBORUnX3zxhTZs2BDcPZGklJQUHTp0KOT6b7/9VocPH1ZKSkqz9xcTE6OYmJi2GBUAAFjI8R2UU3GyZ88e/e1vf1NiYmLI+aysLNXV1amysjJ4bMOGDWpqalJmZqbT4wAAgDB03jsoR48e1d69e4Nf79u3T9u3b1dCQoJ69uypn/zkJ9q2bZvWrl2rkydPBl9XkpCQoOjoaPXr108jR47U5MmTVVRUpMbGRk2bNk3jx4/nHTwAAEDSBQTK1q1b9eMf/zj4dV5eniRp4sSJ+s1vfqO3335bkjRgwICQ73v//fc1ZMgQSdKqVas0bdo0DRs2TJGRkcrNzdWSJUsucAkAAKCjOe9AGTJkiIwxP3i+pXOnJCQkqKSk5Hz/aAAAcJHgd/EAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA60S5PQAA5/XOX+f2CADQKuygAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDp81D1gIT6qHsDFjh0UAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdc47UDZt2qQxY8YoNTVVERERWrNmTch5Y4zmzp2rnj17KjY2VtnZ2dqzZ0/INYcPH9aECRPk8XgUHx+vSZMm6ejRo61aCAAA6DjOO1Dq6+t13XXXadmyZc2eX7hwoZYsWaKioiJVVFQoLi5OI0aM0PHjx4PXTJgwQTt37tT69eu1du1abdq0SVOmTLnwVQAAgA4l6ny/IScnRzk5Oc2eM8Zo8eLFmj17tsaOHStJWrlypZKTk7VmzRqNHz9eu3btUmlpqbZs2aJBgwZJkpYuXapRo0Zp0aJFSk1NbcVyAABAR+Doa1D27dsnn8+n7Ozs4DGv16vMzEyVl5dLksrLyxUfHx+ME0nKzs5WZGSkKioqmr3fhoYGBQKBkBsAAOi4HA0Un88nSUpOTg45npycHDzn8/mUlJQUcj4qKkoJCQnBa05XWFgor9cbvKWlpTk5NgAAsExYvIunoKBAfr8/eKuurnZ7JAAA0IYcDZSUlBRJUm1tbcjx2tra4LmUlBQdOnQo5Py3336rw4cPB685XUxMjDweT8gNAAB0XI4GSkZGhlJSUlRWVhY8FggEVFFRoaysLElSVlaW6urqVFlZGbxmw4YNampqUmZmppPjAACAMHXe7+I5evSo9u7dG/x637592r59uxISEpSenq4ZM2Zo/vz56tOnjzIyMjRnzhylpqZq3LhxkqR+/fpp5MiRmjx5soqKitTY2Khp06Zp/PjxvIMHAABIuoBA2bp1q3784x8Hv87Ly5MkTZw4Ua+88opmzZql+vp6TZkyRXV1dbrttttUWlqqLl26BL9n1apVmjZtmoYNG6bIyEjl5uZqyZIlDiwHAAB0BBHGGOP2EOcrEAjI6/XK7/fzehR0SL3z17k9Qljbv2C02yMAaMb5/PwOi3fxAACAiwuBAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOtEuT0A0BH1zl/n9ggXtdb+89+/YLRDkwC4UOygAAAA6xAoAADAOgQKAACwDoECAACsw4tkAeA0vMgWcB87KAAAwDoECgAAsA6BAgAArEOgAAAA6zgeKCdPntScOXOUkZGh2NhYXXHFFfq///s/GWOC1xhjNHfuXPXs2VOxsbHKzs7Wnj17nB4FAACEKccD5amnntLy5cv17LPPateuXXrqqae0cOFCLV26NHjNwoULtWTJEhUVFamiokJxcXEaMWKEjh8/7vQ4AAAgDDn+NuN//etfGjt2rEaP/u5tdr1799Yf//hHbd68WdJ3uyeLFy/W7NmzNXbsWEnSypUrlZycrDVr1mj8+PFOjwQAAMKM4zsot9xyi8rKyrR7925J0scff6wPPvhAOTk5kqR9+/bJ5/MpOzs7+D1er1eZmZkqLy9v9j4bGhoUCARCbgAAoONyfAclPz9fgUBAffv2VadOnXTy5Ek98cQTmjBhgiTJ5/NJkpKTk0O+Lzk5OXjudIWFhZo3b57TowIAAEs5voPyxhtvaNWqVSopKdG2bdu0YsUKLVq0SCtWrLjg+ywoKJDf7w/eqqurHZwYAADYxvEdlEceeUT5+fnB15L0799fX3zxhQoLCzVx4kSlpKRIkmpra9WzZ8/g99XW1mrAgAHN3mdMTIxiYmKcHhUAAFjK8R2UY8eOKTIy9G47deqkpqYmSVJGRoZSUlJUVlYWPB8IBFRRUaGsrCynxwEAAGHI8R2UMWPG6IknnlB6erp+9KMf6aOPPtLTTz+t++67T5IUERGhGTNmaP78+erTp48yMjI0Z84cpaamaty4cU6PAwAAwpDjgbJ06VLNmTNHDz74oA4dOqTU1FTdf//9mjt3bvCaWbNmqb6+XlOmTFFdXZ1uu+02lZaWqkuXLk6PAwAAwlCE+f5HvIaJQCAgr9crv98vj8fj9jjAGXrnr3N7BLho/4LRbo8AWOl8fn7zu3gAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGCdKLcHAGzTO3+d2yMAwEWPHRQAAGAdAgUAAFiHp3gAwGFOPE24f8FoByYBwhc7KAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsw+egoMPho+oBIPyxgwIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA67RJoHz55Zf6xS9+ocTERMXGxqp///7aunVr8LwxRnPnzlXPnj0VGxur7Oxs7dmzpy1GAQAAYcjxQPn666916623qnPnznr33Xf173//W7/73e90ySWXBK9ZuHChlixZoqKiIlVUVCguLk4jRozQ8ePHnR4HAACEIcc/qO2pp55SWlqaiouLg8cyMjKC/9sYo8WLF2v27NkaO3asJGnlypVKTk7WmjVrNH78eKdHAgAAYcbxHZS3335bgwYN0k9/+lMlJSXp+uuv14svvhg8v2/fPvl8PmVnZwePeb1eZWZmqry83OlxAABAGHI8UD7//HMtX75cffr00XvvvacHHnhADz30kFasWCFJ8vl8kqTk5OSQ70tOTg6eO11DQ4MCgUDIDQAAdFyOP8XT1NSkQYMG6cknn5QkXX/99fr0009VVFSkiRMnXtB9FhYWat68eU6OCQAALOb4DkrPnj119dVXhxzr16+fDhw4IElKSUmRJNXW1oZcU1tbGzx3uoKCAvn9/uCturra6bEBAIBFHA+UW2+9VVVVVSHHdu/ercsvv1zSdy+YTUlJUVlZWfB8IBBQRUWFsrKymr3PmJgYeTyekBsAAOi4HH+KZ+bMmbrlllv05JNP6u6779bmzZv1wgsv6IUXXpAkRUREaMaMGZo/f7769OmjjIwMzZkzR6mpqRo3bpzT4wAAgDDkeKDceOONevPNN1VQUKDf/va3ysjI0OLFizVhwoTgNbNmzVJ9fb2mTJmiuro63XbbbSotLVWXLl2cHgcAAIShCGOMcXuI8xUIBOT1euX3+3m6B2fonb/O7RGAVtu/YLTbIwCOO5+f3/wuHgAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHWi3B4AAHCm3vnrWvX9+xeMdmgSwB3soAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACs0+aBsmDBAkVERGjGjBnBY8ePH9fUqVOVmJiobt26KTc3V7W1tW09CgAACBNtGihbtmzR888/r2uvvTbk+MyZM/XOO+9o9erV2rhxo2pqanTXXXe15SgAACCMtNlH3R89elQTJkzQiy++qPnz5weP+/1+vfTSSyopKdHQoUMlScXFxerXr58+/PBD3XzzzW01EsJEaz/iGwAQ/tpsB2Xq1KkaPXq0srOzQ45XVlaqsbEx5Hjfvn2Vnp6u8vLythoHAACEkTbZQXnttde0bds2bdmy5YxzPp9P0dHRio+PDzmenJwsn8/X7P01NDSooaEh+HUgEHB0XgAAYBfHd1Cqq6v18MMPa9WqVerSpYsj91lYWCiv1xu8paWlOXK/AADATo4HSmVlpQ4dOqQbbrhBUVFRioqK0saNG7VkyRJFRUUpOTlZJ06cUF1dXcj31dbWKiUlpdn7LCgokN/vD96qq6udHhsAAFjE8ad4hg0bpk8++STk2L333qu+ffvq0UcfVVpamjp37qyysjLl5uZKkqqqqnTgwAFlZWU1e58xMTGKiYlxelQAAGApxwOle/fuuuaaa0KOxcXFKTExMXh80qRJysvLU0JCgjwej6ZPn66srCzewQMAACS14duMW/LMM88oMjJSubm5amho0IgRI/Tcc8+5MQoAALBQhDHGuD3E+QoEAvJ6vfL7/fJ4PG6PA4fxOShA6+1fMNrtEYAznM/Pb34XDwAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACsQ6AAAADrECgAAMA6BAoAALBOlNsDoOPpnb/O7REAAGGOHRQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh4+6xxn4qHog/LX23+P9C0Y7NAlwYRzfQSksLNSNN96o7t27KykpSePGjVNVVVXINcePH9fUqVOVmJiobt26KTc3V7W1tU6PAgAAwpTjgbJx40ZNnTpVH374odavX6/GxkYNHz5c9fX1wWtmzpypd955R6tXr9bGjRtVU1Oju+66y+lRAABAmHL8KZ7S0tKQr1955RUlJSWpsrJSgwcPlt/v10svvaSSkhINHTpUklRcXKx+/frpww8/1M033+z0SAAAIMy0+Ytk/X6/JCkhIUGSVFlZqcbGRmVnZwev6du3r9LT01VeXt7sfTQ0NCgQCITcAABAx9WmgdLU1KQZM2bo1ltv1TXXXCNJ8vl8io6OVnx8fMi1ycnJ8vl8zd5PYWGhvF5v8JaWltaWYwMAAJe1aaBMnTpVn376qV577bVW3U9BQYH8fn/wVl1d7dCEAADARm32NuNp06Zp7dq12rRpk3r16hU8npKSohMnTqiuri5kF6W2tlYpKSnN3ldMTIxiYmLaalQAAGAZx3dQjDGaNm2a3nzzTW3YsEEZGRkh5wcOHKjOnTurrKwseKyqqkoHDhxQVlaW0+MAAIAw5PgOytSpU1VSUqK33npL3bt3D76uxOv1KjY2Vl6vV5MmTVJeXp4SEhLk8Xg0ffp0ZWVl8Q4eAAAgqQ0CZfny5ZKkIUOGhBwvLi7Wr371K0nSM888o8jISOXm5qqhoUEjRozQc8895/QoAAAgTEUYY4zbQ5yvQCAgr9crv98vj8fj9jgdDh91D6C1+Kh8NOd8fn7zywIBAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdQgUAABgHQIFAABYh0ABAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdaLcHgDO652/zu0RAABoFXZQAACAdQgUAABgHZ7iAQA4rrVPNe9fMNqhSRCu2EEBAADWIVAAAIB1CBQAAGAdAgUAAFiHQAEAANYhUAAAgHUIFAAAYB0CBQAAWIdAAQAA1iFQAACAdfioewCAdfiofLCDAgAArEOgAAAA6/AUDwCgw+EpovDHDgoAALAOgQIAAKzDUzzNYGsQAAB3sYMCAACsQ6AAAADr8BSPhVr7FBMAILzxUgOXd1CWLVum3r17q0uXLsrMzNTmzZvdHAcAAFjCtR2U119/XXl5eSoqKlJmZqYWL16sESNGqKqqSklJSW6NBQBAq7m9E+7En+/2LoxrOyhPP/20Jk+erHvvvVdXX321ioqK1LVrV7388stujQQAACzhyg7KiRMnVFlZqYKCguCxyMhIZWdnq7y8/IzrGxoa1NDQEPza7/dLkgKBQJvM19RwrFXf39q5WvvnAwBah/+Ot83P2FP3aYw567WuBMr//vc/nTx5UsnJySHHk5OT9dlnn51xfWFhoebNm3fG8bS0tDabsTW8i92eAADQGvx3vG3/GRw5ckRer7fFa8LiXTwFBQXKy8sLft3U1KTDhw8rMTFRERER7TpLIBBQWlqaqqur5fF42vXPdhtrZ+2s/eLB2ll7W6zdGKMjR44oNTX1rNe6EiiXXnqpOnXqpNra2pDjtbW1SklJOeP6mJgYxcTEhByLj49vyxHPyuPxXHR/cU9h7az9YsPaWfvFpi3Xfradk1NceZFsdHS0Bg4cqLKysuCxpqYmlZWVKSsry42RAACARVx7iicvL08TJ07UoEGDdNNNN2nx4sWqr6/Xvffe69ZIAADAEq4Fys9+9jN99dVXmjt3rnw+nwYMGKDS0tIzXjhrm5iYGD3++ONnPOV0MWDtrP1iw9pZ+8XGprVHmHN5rw8AAEA74pcFAgAA6xAoAADAOgQKAACwDoECAACsQ6D8gOXLl+vaa68NflhNVlaW3n333eD548ePa+rUqUpMTFS3bt2Um5t7xgfPhaOzrfuFF17QkCFD5PF4FBERobq6OveGdVhLaz98+LCmT5+uq666SrGxsUpPT9dDDz0U/L1Q4e5sj/v999+vK664QrGxserRo4fGjh3b7K+lCEdnW/spxhjl5OQoIiJCa9asaf9B28DZ1j5kyBBFRESE3H7961+7OLFzzuVxLy8v19ChQxUXFyePx6PBgwfrm2++cWli57S09v3795/xmJ+6rV69ul3nJFB+QK9evbRgwQJVVlZq69atGjp0qMaOHaudO3dKkmbOnKl33nlHq1ev1saNG1VTU6O77rrL5alb72zrPnbsmEaOHKnHHnvM5Umd19Laa2pqVFNTo0WLFunTTz/VK6+8otLSUk2aNMntsR1xtsd94MCBKi4u1q5du/Tee+/JGKPhw4fr5MmTLk/eemdb+ymLFy9u91+t0dbOZe2TJ0/WwYMHg7eFCxe6OLFzzrb28vJyjRw5UsOHD9fmzZu1ZcsWTZs2TZGR4f9js6W1p6WlhTzeBw8e1Lx589StWzfl5OS076AG5+ySSy4xv//9701dXZ3p3LmzWb16dfDcrl27jCRTXl7u4oRt49S6v+/99983kszXX3/tzlDtpLm1n/LGG2+Y6Oho09jY2M5TtY+W1v7xxx8bSWbv3r3tPFX7OH3tH330kbnsssvMwYMHjSTz5ptvujdcG/v+2m+//Xbz8MMPuztQO/r+2jMzM83s2bNdnqj9tPTv+4ABA8x9993XzhMZE/4p2A5Onjyp1157TfX19crKylJlZaUaGxuVnZ0dvKZv375KT09XeXm5i5M66/R1X0zOZe1+v18ej0dRUWHxOzfP2dnWXl9fr+LiYmVkZFj7G8UvVHNrP3bsmO655x4tW7as2d8V1lH80OO+atUqXXrppbrmmmtUUFCgY8eOuThl2zh97YcOHVJFRYWSkpJ0yy23KDk5Wbfffrs++OADt0d13Nn+fa+srNT27dvd2S1u9yQKIzt27DBxcXGmU6dOxuv1mnXr1hljjFm1apWJjo4+4/obb7zRzJo1q73HdNwPrfv7OuoOyrms3RhjvvrqK5Oenm4ee+yxdp6w7Zxt7cuWLTNxcXFGkrnqqqs61O5JS2ufMmWKmTRpUvBrdbAdlJbW/vzzz5vS0lKzY8cO8+qrr5rLLrvM3HnnnS5O66wfWnt5ebmRZBISEszLL79stm3bZmbMmGGio6PN7t27XZ7aGef637oHHnjA9OvXr52n+w6B0oKGhgazZ88es3XrVpOfn28uvfRSs3Pnzg4fKD+07u/rqIFyLmv3+/3mpptuMiNHjjQnTpxwaVLnnW3tdXV1Zvfu3Wbjxo1mzJgx5oYbbjDffPONixM754fW/tZbb5krr7zSHDlyJHhtRwuUc/k7f0pZWVmHemrvh9b+z3/+00gyBQUFIdf379/f5OfnuzSts87lcT927Jjxer1m0aJFrsxIoJyHYcOGmSlTpgT/JT39h3N6erp5+umn3RmuDZ1a9/d11EA53elrDwQCJisrywwbNqzD/HD+Ic097qc0NDSYrl27mpKSknaeqn2cWvvDDz9sIiIiTKdOnYI3SSYyMtLcfvvtbo/ZJlp63I8ePWokmdLS0naeqn2cWvvnn39uJJk//OEPIefvvvtuc88997g0Xdtq7nFfuXKl6dy5szl06JArM/EalPPQ1NSkhoYGDRw4UJ07d1ZZWVnwXFVVlQ4cONAhX6txat0Xo++vPRAIaPjw4YqOjtbbb7+tLl26uDxd22rpcTff/Z+bDvv34tTa8/PztWPHDm3fvj14k6RnnnlGxcXF7g7ZRlp63E+tv2fPnu04Ufs5tfbevXsrNTVVVVVVIed3796tyy+/3KXp2lZzj/tLL72kO+64Qz169HBlpo716j4HFRQUKCcnR+np6Tpy5IhKSkr097//Xe+99568Xq8mTZqkvLw8JSQkyOPxaPr06crKytLNN9/s9uit0tK6Jcnn88nn82nv3r2SpE8++UTdu3dXenq6EhIS3By91Vpa+6k4OXbsmF599VUFAgEFAgFJUo8ePdSpUyeXp2+dltb++eef6/XXX9fw4cPVo0cP/fe//9WCBQsUGxurUaNGuT16q7W09pSUlGZfGJuenq6MjAwXpnVWS2v/z3/+o5KSEo0aNUqJiYnasWOHZs6cqcGDB+vaa691e/RWa2ntEREReuSRR/T444/ruuuu04ABA7RixQp99tln+tOf/uT26K12tv/OS9LevXu1adMm/eUvf3FvUFf2bcLAfffdZy6//HITHR1tevToYYYNG2b++te/Bs9/88035sEHHzSXXHKJ6dq1q7nzzjvNwYMHXZzYGWdb9+OPP24knXErLi52b2iHtLT2U09pNXfbt2+fu4M7oKW1f/nllyYnJ8ckJSWZzp07m169epl77rnHfPbZZy5P7Yyz/Z0/nTrQa1BaWvuBAwfM4MGDTUJCgomJiTFXXnmleeSRR4zf73d5amecy+NeWFhoevXqZbp27WqysrLMP/7xD5emdda5rL2goMCkpaWZkydPujSlMRHGGONOGgEAADSP16AAAADrECgAAMA6BAoAALAOgQIAAKxDoAAAAOsQKAAAwDoECgAAsA6BAgAArEOgAAAA6xAoAADAOgQKAACwDoECAACs8/8fu0/PYQGRKAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure()\n", + "plt.hist(mock[\"mu_true\"], bins=\"auto\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": {}, + "outputs": [], + "source": [ + "from jax import numpy as jnp\n", + "from jax.scipy.special import logsumexp\n", + "\n", + "def normal_logpdf(x, loc, scale):\n", + " \"\"\"Log of the normal probability density function.\"\"\"\n", + " return (-0.5 * ((x - loc) / scale)**2\n", + " - jnp.log(scale) - 0.5 * jnp.log(2 * jnp.pi))" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": {}, + "outputs": [], + "source": [ + "mu_true = np.copy(mock[\"mu_true\"])\n", + "mu_calibration = np.copy(mock[\"mu_calibration\"])\n", + "e_mu = np.copy(mock[\"e_mu_calibration\"])\n", + "\n", + "mu_calibration = np.stack([mu_calibration, mu_calibration])\n", + "e_mu_calibration = np.stack([e_mu, e_mu])\n", + "\n", + "mu_calibration[0, 100:] = np.nan\n", + "e_mu_calibration[0, 100:] = np.nan\n", + "\n", + "mu_calibration[1, 50:] = np.nan\n", + "e_mu_calibration[1, 50:] = np.nan" + ] + }, + { + "cell_type": "code", + "execution_count": 105, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "h = 0.7\n", + "\n", + "# Now, the rest of the code except the calibration likelihood\n", + "# uses the distance modulus in units of h\n", + "mu_true_h = mu_true + 5 * jnp.log10(h)\n", + "\n", + "# Calculate the log-likelihood of the calibration, but the\n", + "# shape is `(n_calibrators, n_data)`.\n", + "ll_calibration = normal_logpdf(\n", + " mu_calibration, mu_true[None, :],\n", + " e_mu_calibration)\n", + "\n", + "# Create a mask for valid (non-NaN) log-likelihoods\n", + "calibration_mask = ~jnp.isnan(ll_calibration)\n", + "\n", + "# Replace NaN values with zero (or another neutral value) for safety\n", + "ll_calibration_clean = jnp.where(calibration_mask, ll_calibration, 0.0)\n", + "\n", + "# Count the number of valid calibrators for each galaxy (non-NaN entries)\n", + "counts = jnp.sum(calibration_mask, axis=0)\n", + "\n", + "# Now apply logsumexp only to the valid log-likelihoods\n", + "ll_calibration_sum = jnp.where(\n", + " counts > 0,\n", + " logsumexp(ll_calibration_clean, axis=0) - jnp.log(counts),\n", + " 0.0 # Return zero likelihood if no valid calibrators\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 120, + "metadata": {}, + "outputs": [], + "source": [ + "from jax.lax import cond\n", + "def ll_calibration(mu_calibration, mu_true, e_mu):\n", + " # Use jnp.where to apply element-wise conditional logic\n", + " return jnp.where(\n", + " jnp.isfinite(mu_calibration), # Check for finite values\n", + " normal_logpdf(mu_calibration, mu_true, e_mu), # Use valid values\n", + " 0.0 # Return 0 for invalid (non-finite) values\n", + " )\n" + ] + }, + { + "cell_type": "code", + "execution_count": 121, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "33.15014277245038" + ] + }, + "execution_count": 121, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mu_calibration[0, 0]" + ] + }, + { + "cell_type": "code", + "execution_count": 127, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Array([[1.5757915, 1.3846258, 2.0763617, ..., 0. , 0. ,\n", + " 0. ],\n", + " [1.5757915, 1.3846258, 2.0763617, ..., 0. , 0. ,\n", + " 0. ]], dtype=float32)" + ] + }, + "execution_count": 127, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ll_calibration(mu_calibration, mu_true[None, :], e_mu)" + ] + }, + { + "cell_type": "code", + "execution_count": 125, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([33.15014277, 33.36291487, 31.68284461, ..., nan,\n", + " nan, nan]),\n", + " array([33.10009268, 33.30408597, 31.68137438, ..., 31.3564011 ,\n", + " 32.63194483, 33.72852162]))" + ] + }, + "execution_count": 125, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mu_calibration[0], mu_true" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 69, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[1.5757915 1.3846258 2.0763617 ... nan nan nan]\n", + " [1.5757915 1.3846258 2.0763617 ... nan nan nan]]\n", + "\n", + "[[1.5757915 1.3846258 2.0763617 ... -inf -inf -inf]\n", + " [1.5757915 1.3846258 2.0763617 ... -inf -inf -inf]]\n", + "\n", + "[1.5757914 1.3846259 2.0763617 ... 0. 0. 0. ]\n", + "\n" + ] + } + ], + "source": [ + "ll = normal_logpdf(mu_calibration, mu_true[None, :], e_mu)\n", + "\n", + "print(ll)\n", + "print()\n", + "\n", + "\n", + "mask = ~jnp.isnan(ll)\n", + "ll = jnp.where(mask, ll, -jnp.inf)\n", + "\n", + "print(ll)\n", + "print()\n", + "\n", + "counts = jnp.sum(mask, axis=0)\n", + "ll = jnp.where(counts > 0, logsumexp(ll, axis=0) - jnp.log(counts), 0.)\n", + "\n", + "print(ll)\n", + "print()" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Array(100, dtype=int32)" + ] + }, + "execution_count": 68, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "jnp.sum(~jnp.isnan(ll))" + ] + }, + { + "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 +} diff --git a/notebooks/flow/flow_mock_pv.ipynb b/notebooks/flow/mock_csiborg.ipynb similarity index 100% rename from notebooks/flow/flow_mock_pv.ipynb rename to notebooks/flow/mock_csiborg.ipynb diff --git a/scripts/flow/flow_validation.py b/scripts/flow/flow_validation.py index 83d8669..4e4a21f 100644 --- a/scripts/flow/flow_validation.py +++ b/scripts/flow/flow_validation.py @@ -116,6 +116,9 @@ def get_models(ksim, get_model_kwargs, mag_selection, void_kwargs, "Pantheon+_groups", "Pantheon+_groups_zSN", "Pantheon+_zSN"]: fpath = join(folder, "PV_compilation.hdf5") + elif "Carrick2MTFmock" in cat: + ki = cat.split("_")[-1] + fpath =f"/mnt/extraspace/rstiskalek/csiborg_postprocessing/flow_mock/Carrick2MTFmock_seed{ki}.hdf5" # noqa elif "CF4_TFR" in cat: fpath = join(folder, "PV/CF4/CF4_TF-distances.hdf5") elif cat in ["CF4_GroupAll"]: @@ -186,7 +189,7 @@ def run_model(model, nsteps, nburn, model_kwargs, out_folder, neg_ln_evidence_err = (jax.numpy.nan, jax.numpy.nan) fname = join(out_folder, fname) - print(f"Saving results to `{fname}`.") + print(f"Saving results: `{fname}`.") with File(fname, "w") as f: # Write samples grp = f.create_group("samples") @@ -203,15 +206,22 @@ def run_model(model, nsteps, nburn, model_kwargs, out_folder, grp.create_dataset("neg_lnZ_harmonic", data=neg_ln_evidence) grp.create_dataset("neg_lnZ_harmonic_err", data=neg_ln_evidence_err) - fname_summary = fname.replace(".hdf5", ".txt") - print(f"Saving summary to `{fname_summary}`.") - with open(fname_summary, 'w') as f: + fname_config = fname.replace(".hdf5", "_config.txt") + print(f"Saving configuration: `{fname_config}`.") + with open(fname_config, 'w') as f: original_stdout = sys.stdout sys.stdout = f print("User parameters:") for kwargs in kwargs_print: print_variables(kwargs.keys(), kwargs.values()) + sys.stdout = original_stdout + + fname_summary = fname.replace(".hdf5", "_summary.txt") + print(f"Saving summary: `{fname_summary}`.") + with open(fname_summary, 'w') as f: + original_stdout = sys.stdout + sys.stdout = f print("HMC summary:") print(f"{'BIC':<20} {BIC}") @@ -238,7 +248,7 @@ def get_distmod_hyperparams(catalogue, sample_alpha, sample_mag_dipole): "alpha_min": alpha_min, "alpha_max": alpha_max, "sample_alpha": sample_alpha } - elif catalogue in ["SFI_gals", "2MTF"] or "CF4_TFR" in catalogue or "IndranilVoidTFRMock" in catalogue: # noqa + elif catalogue in ["SFI_gals", "2MTF"] or "CF4_TFR" in catalogue or "IndranilVoidTFRMock" in catalogue or "Carrick2MTFmock" in catalogue: # noqa return {"e_mu_min": 0.001, "e_mu_max": 1.0, "a_mean": -21., "a_std": 5.0, "b_mean": -5.95, "b_std": 4.0, @@ -247,6 +257,7 @@ def get_distmod_hyperparams(catalogue, sample_alpha, sample_mag_dipole): "sample_a_dipole": sample_mag_dipole, "alpha_min": alpha_min, "alpha_max": alpha_max, "sample_alpha": sample_alpha, + "sample_curvature": False if "Carrick2MTFmock" in catalogue else True, # noqa } elif catalogue in ["CF4_GroupAll"]: return {"e_mu_min": 0.001, "e_mu_max": 1.0, @@ -294,8 +305,8 @@ if __name__ == "__main__": ########################################################################### # `None` means default behaviour - nsteps = 2_000 - nburn = 2_000 + nsteps = 1_500 + nburn = 1_500 zcmb_min = None zcmb_max = 0.05 nchains_harmonic = 10 @@ -310,9 +321,16 @@ if __name__ == "__main__": sample_mag_dipole = False wo_num_dist_marginalisation = False absolute_calibration = None - calculate_harmonic = (False if inference_method == "bayes" else True) and (not wo_num_dist_marginalisation) # noqa + calculate_harmonic = (False if (inference_method == "bayes") else True) and (not wo_num_dist_marginalisation) # noqa sample_h = True if absolute_calibration is not None else False + # These mocks are generated without a density field, so there is no + # inhomogeneous Malmquist and we also do not need evidences. + for catalogue in ARGS.catalogue: + if "Carrick2MTFmock" in catalogue: + sample_alpha = False + calculate_harmonic = False + fname_kwargs = {"inference_method": inference_method, "smooth": ARGS.ksmooth, "nsim": ARGS.ksim, @@ -375,7 +393,7 @@ if __name__ == "__main__": "Vmono_min": -1000, "Vmono_max": 1000, "beta_min": -10.0, "beta_max": 10.0, "sigma_v_min": 1.0, "sigma_v_max": 1000 if "IndranilVoid_" in ARGS.simname else 750., # noqa - "h_min": 0.01, "h_max": 5.0, + "h_min": 0.25, "h_max": 1., "no_Vext": False if no_Vext is None else no_Vext, # noqa "sample_Vmag_vax": sample_Vmag_vax, "sample_Vmono": sample_Vmono, diff --git a/scripts/flow/flow_validation.sh b/scripts/flow/flow_validation.sh index f3dbd1f..4ad6873 100755 --- a/scripts/flow/flow_validation.sh +++ b/scripts/flow/flow_validation.sh @@ -38,10 +38,10 @@ fi # for simname in "IndranilVoid_exp" "IndranilVoid_gauss" "IndranilVoid_mb"; do -for simname in "IndranilVoid_exp"; do +for simname in "Carrick2015"; do # for simname in "Carrick2015" "Lilow2024" "csiborg2_main" "csiborg2X" "manticore_2MPP_N128_DES_V1" "CF4" "CLONES"; do # for catalogue in "LOSS" "Foundation" "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do - for catalogue in "2MTF"; do + for catalogue in "CF4_TFR_w1"; do # for catalogue in "CF4_TFR_i" "CF4_TFR_w1"; do # for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do for ksim in "none"; do