diff --git a/environment.yml b/environment.yml index fb34e729..485010e0 100644 --- a/environment.yml +++ b/environment.yml @@ -23,6 +23,7 @@ dependencies: - pip: - autoclasstoc - cobaya + - pygobject-stubs - portalocker - pybobyqa - pyccl >= 2.8.0 diff --git a/examples/des_y1_3x2pt/numcosmo_run.py b/examples/des_y1_3x2pt/numcosmo_run.py index b4c26b02..1dd3e609 100755 --- a/examples/des_y1_3x2pt/numcosmo_run.py +++ b/examples/des_y1_3x2pt/numcosmo_run.py @@ -2,96 +2,230 @@ """Example of running the DES Y1 3x2pt likelihood using NumCosmo.""" import math +import os +from typing import Tuple +import argparse +from pathlib import Path + import yaml from numcosmo_py import Nc, Ncm +from numcosmo_py.sampling.esmcmc import create_esmcmc from firecrown.connector.numcosmo.numcosmo import MappingNumCosmo, NumCosmoFactory from firecrown.connector.numcosmo.model import define_numcosmo_model from firecrown.likelihood.likelihood import NamedParameters +# Any NumCosmo model should be loaded before NumCosmo cfg_init +# is called. Otherwise, the NumCosmo model will not be registered +# in the model set. This is required for the NumCosmo MPI support +# to work properly. -Ncm.cfg_init() -with open(r"numcosmo_firecrown_model.yml", "r", encoding="utf-8") as modelfile: +module_path = Path(os.path.dirname(__file__)) +with open( + module_path / r"numcosmo_firecrown_model.yml", "r", encoding="utf-8" +) as modelfile: ncmodel = yaml.load(modelfile, Loader=yaml.Loader) NcFirecrown = define_numcosmo_model(ncmodel) -cosmo = Nc.HICosmoDEXcdm(massnu_length=1) -cosmo.omega_x2omega_k() -cosmo.param_set_by_name("H0", 68.2) -cosmo.param_set_by_name("Omegak", 0.0) -cosmo.param_set_by_name("Omegab", 0.022558514 / 0.682**2) -cosmo.param_set_by_name("Omegac", 0.118374058 / 0.682**2) -cosmo.param_set_by_name("massnu_0", 0.06) -cosmo.param_set_by_name("ENnu", 2.0328) -cosmo.param_set_by_name("Yp", 0.2454) -cosmo.param_set_by_name("w", -1.0) - -prim = Nc.HIPrimPowerLaw.new() -prim.param_set_by_name("ln10e10ASA", math.log(1.0e10 * 2.0e-09)) -prim.param_set_by_name("n_SA", 0.971) - -reion = Nc.HIReionCamb.new() -reion.set_z_from_tau(cosmo, 0.0561) - -cosmo.add_submodel(prim) -cosmo.add_submodel(reion) - -p_ml = Nc.PowspecMLTransfer.new(Nc.TransferFuncEH.new()) -p_mnl = Nc.PowspecMNLHaloFit.new(p_ml, 3.0, 1.0e-5) -dist = Nc.Distance.new(6.0) -dist.comoving_distance_spline.set_reltol(1.0e-5) - -map_cosmo = MappingNumCosmo( - require_nonlinear_pk=True, - p_ml=p_ml, - p_mnl=p_mnl, - dist=dist, - model_list=["NcFirecrown"], -) - -nc_factory = NumCosmoFactory("des_y1_3x2pt.py", NamedParameters(), map_cosmo) - -fc = NcFirecrown() -# fc.params_set_default_ftype() - -mset = Ncm.MSet() -mset.set(cosmo) -mset.set(fc) - -fc_data = nc_factory.get_data() - -dset = Ncm.Dataset() -dset.append_data(fc_data) - -lh = Ncm.Likelihood(dataset=dset) - -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_delta_z", -0.001, 0.016) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src1_delta_z", -0.019, 0.013) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src2_delta_z", +0.009, 0.011) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src3_delta_z", -0.018, 0.022) - -lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens0_delta_z", +0.001, 0.008) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens1_delta_z", +0.002, 0.007) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens2_delta_z", +0.001, 0.007) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens3_delta_z", +0.003, 0.010) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens4_delta_z", +0.000, 0.010) - -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) -lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) - -fit = Ncm.Fit.new( - Ncm.FitType.NLOPT, - "ln-neldermead", - lh, - mset, - Ncm.FitGradType.NUMDIFF_FORWARD, -) - -mset.pretty_log() -fit.run_restart(Ncm.FitRunMsgs.FULL, 1.0e-3, 0.0, None, None) + +def setup_numcosmo_cosmology() -> Tuple[Nc.HICosmo, MappingNumCosmo]: + """Setup NumCosmo cosmology. + + Creates a NumCosmo cosmology object. + """ + + cosmo = Nc.HICosmoDEXcdm(massnu_length=1) + cosmo.omega_x2omega_k() + cosmo.param_set_by_name("H0", 68.2) + cosmo.param_set_by_name("Omegak", 0.0) + cosmo.param_set_by_name("Omegab", 0.022558514 / 0.682**2) + cosmo.param_set_by_name("Omegac", 0.118374058 / 0.682**2) + cosmo.param_set_by_name("massnu_0", 0.06) + cosmo.param_set_by_name("ENnu", 2.0328) + cosmo.param_set_by_name("Yp", 0.2454) + cosmo.param_set_by_name("w", -1.0) + + prim = Nc.HIPrimPowerLaw.new() + prim.param_set_by_name("ln10e10ASA", math.log(1.0e10 * 2.0e-09)) + prim.param_set_by_name("n_SA", 0.971) + + reion = Nc.HIReionCamb.new() + reion.set_z_from_tau(cosmo, 0.0561) + + cosmo.add_submodel(prim) + cosmo.add_submodel(reion) + + p_ml = Nc.PowspecMLTransfer.new(Nc.TransferFuncEH.new()) + p_mnl = Nc.PowspecMNLHaloFit.new(p_ml, 3.0, 1.0e-5) + dist = Nc.Distance.new(6.0) + dist.comoving_distance_spline.set_reltol(1.0e-5) + + map_cosmo = MappingNumCosmo( + require_nonlinear_pk=True, + p_ml=p_ml, + p_mnl=p_mnl, + dist=dist, + ) + + return cosmo, map_cosmo + + +def setup_firecrown(map_cosmo: MappingNumCosmo) -> Tuple[Ncm.Model, NumCosmoFactory]: + """Setup Firecrown object.""" + + nc_factory = NumCosmoFactory( + str(module_path / "des_y1_3x2pt.py"), + NamedParameters(), + map_cosmo, + [ncmodel.name], + ) + + fc = NcFirecrown() + + return fc, nc_factory + + +def setup_likelihood( + cosmo: Nc.HICosmo, + fc: Ncm.Model, + nc_factory: NumCosmoFactory, +) -> Tuple[Ncm.Likelihood, Ncm.MSet]: + """Setup the likelihood and model set objects.""" + + mset = Ncm.MSet() + mset.set(cosmo) + mset.set(fc) + + fc_data = nc_factory.get_data() + + dset = Ncm.Dataset() + dset.append_data(fc_data) + + lh = Ncm.Likelihood(dataset=dset) + + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_delta_z", -0.001, 0.016) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src1_delta_z", -0.019, 0.013) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src2_delta_z", +0.009, 0.011) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src3_delta_z", -0.018, 0.022) + + lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens0_delta_z", +0.001, 0.008) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens1_delta_z", +0.002, 0.007) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens2_delta_z", +0.001, 0.007) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens3_delta_z", +0.003, 0.010) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:lens4_delta_z", +0.000, 0.010) + + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) + lh.priors_add_gauss_param_name(mset, "NcFirecrown:src0_mult_bias", +0.012, 0.023) + + return lh, mset + + +def setup_fit(likelihood: Ncm.Likelihood, mset: Ncm.MSet) -> Ncm.Fit: + """Setup the fit object.""" + + fit = Ncm.Fit.new( + Ncm.FitType.NLOPT, + "ln-neldermead", + likelihood, + mset, + Ncm.FitGradType.NUMDIFF_FORWARD, + ) + + return fit + + +def setup_fit_all() -> Ncm.Fit: + """Setup all objects necessary to instantiate the fit object.""" + + cosmo, map_cosmo = setup_numcosmo_cosmology() + fc, nc_factory = setup_firecrown(map_cosmo) + likelihood, mset = setup_likelihood(cosmo, fc, nc_factory) + fit = setup_fit(likelihood, mset) + + return fit + + +def run_test(): + """Run the fit.""" + + fit = setup_fit_all() + mset = fit.peek_mset() + mset.param_set_all_ftype(Ncm.ParamType.FIXED) + mset.pretty_log() + + fit.run_restart(Ncm.FitRunMsgs.FULL, 1.0e-3, 0.0, None, None) + + +def run_compute_best_fit(): + """Run the fit.""" + + fit = setup_fit_all() + mset = fit.peek_mset() + + # Sets the default ftype for all model parameters + for i in range(mset.nmodels()): + model = mset.peek_array_pos(i) + model.params_set_default_ftype() + mset.prepare_fparam_map() + mset.pretty_log() + + fit.run_restart(Ncm.FitRunMsgs.FULL, 1.0e-3, 0.0, None, None) + + +def run_apes_sampler(ssize: int): + """Run the fit.""" + + fit = setup_fit_all() + mset = fit.peek_mset() + + # Sets the default ftype for all model parameters + for i in range(mset.nmodels()): + model = mset.peek_array_pos(i) + model.params_set_default_ftype() + mset.prepare_fparam_map() + + nwalkers = mset.fparam_len * 100 + esmcmc = create_esmcmc( + fit.lh, mset, "des_y1_3x2pt_apes", nwalkers=nwalkers, nthreads=1 + ) + + esmcmc.start_run() + esmcmc.run(ssize // nwalkers) + esmcmc.end_run() + + +if __name__ == "__main__": + Ncm.cfg_init() + + parser = argparse.ArgumentParser(description="Run DES Y1 3x2pt likelihood.") + + parser.add_argument( + "--run-mode", + choices=["test_likelihood", "compute_best_fit", "run_apes_sampler"], + default="test_likelihood", + help="Run mode.", + ) + + parser.add_argument( + "--apes-sampler-ssize", + type=int, + default=10, + help="Number of samples to draw from the APES sampler.", + ) + + args = parser.parse_args() + + if args.run_mode == "test_likelihood": + run_test() + elif args.run_mode == "compute_best_fit": + run_compute_best_fit() + elif args.run_mode == "run_apes_sampler": + run_apes_sampler(args.apes_sampler_ssize) + else: + raise ValueError(f"Unknown run mode: {args.run_mode}") diff --git a/examples/des_y1_3x2pt/numcosmo_run_PT.py b/examples/des_y1_3x2pt/numcosmo_run_PT.py index 66cea480..4935b00c 100755 --- a/examples/des_y1_3x2pt/numcosmo_run_PT.py +++ b/examples/des_y1_3x2pt/numcosmo_run_PT.py @@ -51,10 +51,14 @@ p_ml=p_ml, p_mnl=p_mnl, dist=dist, - model_list=["NcFirecrownPT"], ) -nc_factory = NumCosmoFactory("des_y1_3x2pt_PT.py", NamedParameters(), map_cosmo) +nc_factory = NumCosmoFactory( + "des_y1_3x2pt_PT.py", + NamedParameters(), + map_cosmo, + model_list=["NcFirecrownPT"], +) fc = NcFirecrownPT() # fc.params_set_default_ftype() diff --git a/examples/srd_sn/numcosmo_run.py b/examples/srd_sn/numcosmo_run.py index 7644d0e4..7bb6e6f9 100755 --- a/examples/srd_sn/numcosmo_run.py +++ b/examples/srd_sn/numcosmo_run.py @@ -56,14 +56,15 @@ dist = Nc.Distance.new(6.0) dist.comoving_distance_spline.set_reltol(1.0e-5) -map_cosmo = MappingNumCosmo( - require_nonlinear_pk=False, dist=dist, model_list=["NcFirecrownSNIa"] -) +map_cosmo = MappingNumCosmo(require_nonlinear_pk=False, dist=dist) sacc_file = os.path.expandvars("${FIRECROWN_DIR}/examples/srd_sn/srd-y1-converted.sacc") nc_factory = NumCosmoFactory( - "sn_srd.py", NamedParameters({"sacc_file": sacc_file}), map_cosmo + "sn_srd.py", + NamedParameters({"sacc_file": sacc_file}), + map_cosmo, + model_list=["NcFirecrownSNIa"], ) fc = NcFirecrownSNIa() diff --git a/firecrown/connector/mapping.py b/firecrown/connector/mapping.py index 90dd4e5f..a77b25ae 100644 --- a/firecrown/connector/mapping.py +++ b/firecrown/connector/mapping.py @@ -143,7 +143,7 @@ def set_params( self.Omega_b = Omega_b self.h = h - if A_s is not None and sigma8 is not None: + if (A_s is not None) == (sigma8 is not None): raise ValueError("Exactly one of A_s and sigma8 must be supplied") if sigma8 is None: self.A_s = A_s @@ -210,9 +210,6 @@ class MappingCLASS(Mapping): missing classes. """ - def get_params_names(self): - return [] - class MappingCosmoSIS(Mapping): """Mapping support for CosmoSIS.""" diff --git a/firecrown/connector/numcosmo/numcosmo.py b/firecrown/connector/numcosmo/numcosmo.py index 03567b81..e60402c0 100644 --- a/firecrown/connector/numcosmo/numcosmo.py +++ b/firecrown/connector/numcosmo/numcosmo.py @@ -6,10 +6,12 @@ """ from typing import Dict, Union, List, Any, Optional +import pickle +import base64 import numpy as np import pyccl as ccl -from numcosmo_py import Nc, Ncm +from numcosmo_py import Nc, Ncm, GObject from firecrown.likelihood.likelihood import load_likelihood from firecrown.likelihood.likelihood import Likelihood @@ -20,25 +22,107 @@ from firecrown.modeling_tools import ModelingTools -class MappingNumCosmo(Mapping): +class MappingNumCosmo(GObject.Object): """Mapping support for NumCosmo, this is a subclass of Mapping that provides a mapping from a NumCosmo Cosmological model to a CCL cosmology. It alsos convert NumCosmo models to `ParamsMap`s.""" + __gtype_name__ = "FirecrownMappingNumCosmo" + def __init__( self, - *, require_nonlinear_pk: bool = False, p_ml: Optional[Nc.PowspecML] = None, p_mnl: Optional[Nc.PowspecMNL] = None, - dist: Nc.Distance, - model_list: List[str], + dist: Optional[Nc.Distance] = None, ): - super().__init__(require_nonlinear_pk=require_nonlinear_pk) - self.p_ml = p_ml - self.p_mnl = p_mnl - self.dist = dist - self.model_list = model_list + """Initialize a MappingNumCosmo object.""" + super().__init__( # type: ignore + require_nonlinear_pk=require_nonlinear_pk, p_ml=p_ml, p_mnl=p_mnl, dist=dist + ) + self.mapping: Mapping + self._mapping_name: str + self._p_ml: Optional[Nc.PowspecML] + self._p_mnl: Optional[Nc.PowspecMNL] + self._dist: Nc.Distance + + def _get_mapping_name(self) -> str: + """Return the mapping name.""" + return self._mapping_name + + def _set_mapping_name(self, value: str): + """Set the mapping name.""" + self._mapping_name = value + self.mapping = Mapping() + + mapping_name = GObject.Property( + type=str, + default="default", + flags=GObject.ParamFlags.READWRITE | GObject.ParamFlags.CONSTRUCT_ONLY, + getter=_get_mapping_name, + setter=_set_mapping_name, + ) + + def _get_require_nonlinear_pk(self) -> bool: + """Return whether nonlinear power spectra are required.""" + return self.mapping.require_nonlinear_pk + + def _set_require_nonlinear_pk(self, value: bool): + """Set whether nonlinear power spectra are required.""" + self.mapping.require_nonlinear_pk = value + + require_nonlinear_pk = GObject.Property( + type=bool, + default=False, + flags=GObject.ParamFlags.READWRITE, + getter=_get_require_nonlinear_pk, + setter=_set_require_nonlinear_pk, + ) + + def _get_p_ml(self) -> Optional[Nc.PowspecML]: + """Return the NumCosmo PowspecML object.""" + return self._p_ml + + def _set_p_ml(self, value: Optional[Nc.PowspecML]): + """Set the NumCosmo PowspecML object.""" + self._p_ml = value + + p_ml = GObject.Property( + type=Nc.PowspecML, + flags=GObject.ParamFlags.READWRITE, + getter=_get_p_ml, + setter=_set_p_ml, + ) + + def _get_p_mnl(self) -> Optional[Nc.PowspecMNL]: + """Return the NumCosmo PowspecMNL object.""" + return self._p_mnl + + def _set_p_mnl(self, value: Optional[Nc.PowspecMNL]): + """Set the NumCosmo PowspecMNL object.""" + self._p_mnl = value + + p_mnl = GObject.Property( + type=Nc.PowspecMNL, + flags=GObject.ParamFlags.READWRITE, + getter=_get_p_mnl, + setter=_set_p_mnl, + ) + + def _get_dist(self) -> Optional[Nc.Distance]: + """Return the NumCosmo Distance object.""" + return self._dist + + def _set_dist(self, value: Nc.Distance): + """Set the NumCosmo Distance object.""" + self._dist = value + + dist = GObject.Property( + type=Nc.Distance, + flags=GObject.ParamFlags.READWRITE, + getter=_get_dist, + setter=_set_dist, + ) def set_params_from_numcosmo( self, mset: Ncm.MSet @@ -49,11 +133,11 @@ def set_params_from_numcosmo( hi_cosmo = mset.peek(Nc.HICosmo.id()) assert isinstance(hi_cosmo, Nc.HICosmo) - if self.p_ml is not None: - self.p_ml.prepare_if_needed(hi_cosmo) - if self.p_mnl is not None: - self.p_mnl.prepare_if_needed(hi_cosmo) - self.dist.prepare_if_needed(hi_cosmo) + if self._p_ml is not None: + self._p_ml.prepare_if_needed(hi_cosmo) + if self._p_mnl is not None: + self._p_mnl.prepare_if_needed(hi_cosmo) + self._dist.prepare_if_needed(hi_cosmo) h = hi_cosmo.h() Omega_b = hi_cosmo.Omega_b0() @@ -94,7 +178,7 @@ def set_params_from_numcosmo( n_s = hiprim.props.n_SA # pylint: disable=duplicate-code - self.set_params( + self.mapping.set_params( Omega_c=Omega_c, Omega_b=Omega_b, h=h, @@ -116,16 +200,16 @@ def calculate_ccl_args(self, mset: Ncm.MSet): # pylint: disable-msg=too-many-lo hi_cosmo = mset.peek(Nc.HICosmo.id()) assert isinstance(hi_cosmo, Nc.HICosmo) - if self.p_ml: - p_m_spline = self.p_ml.get_spline_2d(hi_cosmo) + if self._p_ml: + p_m_spline = self._p_ml.get_spline_2d(hi_cosmo) z = np.array(p_m_spline.xv.dup_array()) k = np.array(p_m_spline.yv.dup_array()) - scale = self.redshift_to_scale_factor(z) + scale = self.mapping.redshift_to_scale_factor(z) p_k = np.transpose( np.array(p_m_spline.zm.dup_array()).reshape(len(k), len(z)) ) - p_k = self.redshift_to_scale_factor_p_k(p_k) + p_k = self.mapping.redshift_to_scale_factor_p_k(p_k) ccl_args["pk_linear"] = { "a": scale, @@ -133,51 +217,59 @@ def calculate_ccl_args(self, mset: Ncm.MSet): # pylint: disable-msg=too-many-lo "delta_matter:delta_matter": p_k, } - if self.p_mnl: - p_mnl_spline = self.p_mnl.get_spline_2d(hi_cosmo) + if self._p_mnl: + p_mnl_spline = self._p_mnl.get_spline_2d(hi_cosmo) z = np.array(p_mnl_spline.xv.dup_array()) k = np.array(p_mnl_spline.yv.dup_array()) - scale_mpnl = self.redshift_to_scale_factor(z) + scale_mpnl = self.mapping.redshift_to_scale_factor(z) p_mnl = np.transpose( np.array(p_mnl_spline.zm.dup_array()).reshape(len(k), len(z)) ) - p_mnl = self.redshift_to_scale_factor_p_k(p_mnl) + p_mnl = self.mapping.redshift_to_scale_factor_p_k(p_mnl) ccl_args["pk_nonlin"] = { "a": scale_mpnl, "k": k, "delta_matter:delta_matter": p_mnl, } - elif self.require_nonlinear_pk: + elif self.mapping.require_nonlinear_pk: ccl_args["nonlinear_model"] = "halofit" else: ccl_args["nonlinear_model"] = None - d_spline = self.dist.comoving_distance_spline.peek_spline() + d_spline = self._dist.comoving_distance_spline.peek_spline() z_dist = np.array(d_spline.get_xv().dup_array()) c_dist = np.array(d_spline.get_yv().dup_array()) chi = np.flip(c_dist) * hi_cosmo.RH_Mpc() - scale_distances = self.redshift_to_scale_factor(z_dist) + scale_distances = self.mapping.redshift_to_scale_factor(z_dist) h_over_h0 = np.array([hi_cosmo.E(z) for z in reversed(z_dist)]) + # Too many points in the redshift spline can result in scale factors + # that are too close together for CCL to handle. This checks for + # duplicate scale factors and removes them. + a_unique, a_unique_indices = np.unique(scale_distances, return_index=True) + scale_distances = a_unique + chi = chi[a_unique_indices] + h_over_h0 = h_over_h0[a_unique_indices] + ccl_args["background"] = build_ccl_background_dict( a=scale_distances, chi=chi, h_over_h0=h_over_h0 ) + ccl_args["background"] = { "a": scale_distances, "chi": chi, "h_over_h0": h_over_h0, } - return ccl_args - def create_params_map(self, mset: Ncm.MSet) -> ParamsMap: + def create_params_map(self, model_list: List[str], mset: Ncm.MSet) -> ParamsMap: """Create a ParamsMap from a NumCosmo MSet.""" params_map = ParamsMap() - for model_ns in self.model_list: + for model_ns in model_list: mid = mset.get_id_by_ns(model_ns) if mid < 0: raise RuntimeError( @@ -195,7 +287,7 @@ def create_params_map(self, mset: Ncm.MSet) -> ParamsMap: raise RuntimeError( f"The following keys `{shared_keys}` appear " f"in more than one model used by the " - f"module {self.model_list}." + f"module {model_list}." ) params_map = ParamsMap({**params_map, **model_dict}) @@ -207,21 +299,100 @@ class NumCosmoData(Ncm.Data): object virtual methods using the prefix :python:`do_`. This class implement a general likelihood.""" - def __init__( - self, - likelihood: Likelihood, - tools: ModelingTools, - mapping: MappingNumCosmo, - ): + __gtype_name__ = "FirecrownNumCosmoData" + + def __init__(self): super().__init__() - self.likelihood: Likelihood = likelihood - self.tools: ModelingTools = tools + self.likelihood: Likelihood + self.tools: ModelingTools + self.ccl_cosmo: Optional[ccl.Cosmology] = None + self._model_list: List[str] + self._nc_mapping: MappingNumCosmo + self._likelihood_str: Optional[str] = None self.dof: int = 100 self.len: int = 100 - self.mapping = mapping - self.ccl_cosmo: Optional[ccl.Cosmology] = None self.set_init(True) + def _get_model_list(self) -> List[str]: + """Return the list of models.""" + return self._model_list + + def _set_model_list(self, value: List[str]): + """Set the list of models.""" + self._model_list = value + + model_list = GObject.Property( + type=GObject.TYPE_STRV, # type: ignore + flags=GObject.ParamFlags.READWRITE | GObject.ParamFlags.CONSTRUCT, + getter=_get_model_list, + setter=_set_model_list, + ) + + def _get_nc_mapping(self) -> MappingNumCosmo: + """Return the MappingNumCosmo object.""" + return self._nc_mapping + + def _set_nc_mapping(self, value: MappingNumCosmo): + """Set the MappingNumCosmo object.""" + self._nc_mapping = value + + nc_mapping = GObject.Property( + type=MappingNumCosmo, + flags=GObject.ParamFlags.READWRITE | GObject.ParamFlags.CONSTRUCT, + getter=_get_nc_mapping, + setter=_set_nc_mapping, + ) + + def _get_likelihood_str(self) -> Optional[str]: + """Return the likelihood string.""" + return self._likelihood_str + + def _set_likelihood_str(self, value: Optional[str]): + """Set the likelihood string.""" + self._likelihood_str = value + if value is not None: + likelihood_source, build_parameters = pickle.loads( + base64.b64decode(value.encode("ascii")) + ) + likelihood, tools = load_likelihood(likelihood_source, build_parameters) + assert isinstance(likelihood, Likelihood) + self.likelihood = likelihood + self.tools = tools + + likelihood_str = GObject.Property( + type=str, + flags=GObject.ParamFlags.READWRITE | GObject.ParamFlags.CONSTRUCT, + getter=_get_likelihood_str, + setter=_set_likelihood_str, + ) + + @classmethod + def new_from_likelihood( + cls, + likelihood: Likelihood, + model_list: List[str], + tools: ModelingTools, + nc_mapping: MappingNumCosmo, + likelihood_str: Optional[str] = None, + ): + """Initialize a NumCosmoGaussCov object representing a Gaussian likelihood + with a constant covariance.""" + + nc_data: NumCosmoData = GObject.new( + cls, + model_list=model_list, + nc_mapping=nc_mapping, + likelihood_str=None, + ) + + nc_data.likelihood = likelihood + nc_data.tools = tools + # pylint: disable=protected-access + nc_data._likelihood_str = likelihood_str + # pylint: enable=protected-access + + return nc_data + def do_get_length(self): # pylint: disable-msg=arguments-differ """ Implements the virtual Ncm.Data method get_length. @@ -253,10 +424,12 @@ def do_prepare(self, mset: Ncm.MSet): # pylint: disable-msg=arguments-differ self.likelihood.reset() self.tools.reset() - self.mapping.set_params_from_numcosmo(mset) - ccl_args = self.mapping.calculate_ccl_args(mset) - self.ccl_cosmo = ccl.CosmologyCalculator(**self.mapping.asdict(), **ccl_args) - params_map = self.mapping.create_params_map(mset) + self._nc_mapping.set_params_from_numcosmo(mset) + ccl_args = self._nc_mapping.calculate_ccl_args(mset) + self.ccl_cosmo = ccl.CosmologyCalculator( + **self._nc_mapping.mapping.asdict(), **ccl_args + ) + params_map = self._nc_mapping.create_params_map(self.model_list, mset) self.likelihood.update(params_map) self.tools.prepare(self.ccl_cosmo) @@ -276,32 +449,78 @@ class NumCosmoGaussCov(Ncm.DataGaussCov): object virtual methods using the prefix :python:`do_`. This class implement a Gaussian likelihood.""" - def __init__( - self, - likelihood: ConstGaussian, - tools: ModelingTools, - mapping: MappingNumCosmo, - ): - """Initialize a NumCosmoGaussCov object representing a Gaussian likelihood - with a constant covariance.""" - cov = likelihood.get_cov() - nrows, ncols = cov.shape - assert nrows == ncols + __gtype_name__ = "FirecrownNumCosmoGaussCov" - super().__init__(n_points=nrows) + def __init__(self): + """Initialize a NumCosmoGaussCov object. This class is a subclass of + Ncm.DataGaussCov and implements NumCosmo likelihood object virtual + methods using the prefix :python:`do_`. This class implement a Gaussian + likelihood. - self.likelihood: ConstGaussian = likelihood - self.tools: ModelingTools = tools - self.mapping = mapping - self.ccl_cosmo: Optional[ccl.Cosmology] = None + Due to the way GObject works, the constructor must have a `**kwargs` + argument, and the properties must be set after construction. + + In python one should use the `new_from_likelihood` method to create a + NumCosmoGaussCov object from a ConstGaussian object. This constuctor + has the correct signature for type checking. + """ + super().__init__() + self.likelihood: ConstGaussian + self.tools: ModelingTools + self.ccl_cosmo: ccl.Cosmology + self.dof: int + self.len: int + self._model_list: List[str] + self._nc_mapping: MappingNumCosmo + self._likelihood_str: Optional[str] = None + + def _get_model_list(self) -> List[str]: + """Return the list of models.""" + return self._model_list + + def _set_model_list(self, value: List[str]): + """Set the list of models.""" + self._model_list = value + + model_list = GObject.Property( + type=GObject.TYPE_STRV, # type: ignore + flags=GObject.ParamFlags.READWRITE | GObject.ParamFlags.CONSTRUCT, + getter=_get_model_list, + setter=_set_model_list, + ) + + def _get_nc_mapping(self) -> MappingNumCosmo: + """Return the MappingNumCosmo object.""" + return self._nc_mapping + + def _set_nc_mapping(self, value: MappingNumCosmo): + """Set the MappingNumCosmo object.""" + self._nc_mapping = value + + nc_mapping = GObject.Property( + type=MappingNumCosmo, + flags=GObject.ParamFlags.READWRITE | GObject.ParamFlags.CONSTRUCT, + getter=_get_nc_mapping, + setter=_set_nc_mapping, + ) + + def _configure_object(self): + """Configure the object.""" + assert self.likelihood is not None + + cov = self.likelihood.get_cov() + nrows, ncols = cov.shape + assert nrows == ncols + self.set_size(nrows) + self.ccl_cosmo = None self.dof = nrows self.len = nrows self.peek_cov().set_from_array( cov.flatten().tolist() ) # pylint: disable-msg=no-member - data_vector = likelihood.get_data_vector() + data_vector = self.likelihood.get_data_vector() assert len(data_vector) == ncols self.peek_mean().set_array( data_vector.tolist() @@ -309,6 +528,62 @@ def __init__( self.set_init(True) + def _get_likelihood_str(self) -> Optional[str]: + """Return the likelihood string.""" + return self._likelihood_str + + def _set_likelihood_str(self, value: Optional[str]): + """Set the likelihood string.""" + self._likelihood_str = value + if value is not None: + likelihood_source, build_parameters = pickle.loads( + base64.b64decode(value.encode("ascii")) + ) + likelihood, tools = load_likelihood(likelihood_source, build_parameters) + assert isinstance(likelihood, ConstGaussian) + self.likelihood = likelihood + self.tools = tools + self._configure_object() + + likelihood_str = GObject.Property( + type=str, + flags=GObject.ParamFlags.READWRITE | GObject.ParamFlags.CONSTRUCT, + getter=_get_likelihood_str, + setter=_set_likelihood_str, + ) + + @classmethod + def new_from_likelihood( + cls, + likelihood: ConstGaussian, + model_list: List[str], + tools: ModelingTools, + nc_mapping: MappingNumCosmo, + likelihood_str: Optional[str] = None, + ): + """Initialize a NumCosmoGaussCov object representing a Gaussian likelihood + with a constant covariance.""" + + cov = likelihood.get_cov() + nrows, ncols = cov.shape + assert nrows == ncols + + nc_gauss_cov: NumCosmoGaussCov = GObject.new( + cls, + model_list=model_list, + nc_mapping=nc_mapping, + likelihood_str=None, + ) + + nc_gauss_cov.likelihood = likelihood + nc_gauss_cov.tools = tools + # pylint: disable=protected-access + nc_gauss_cov._likelihood_str = likelihood_str + nc_gauss_cov._configure_object() + # pylint: enable=protected-access + + return nc_gauss_cov + def do_get_length(self): # pylint: disable-msg=arguments-differ """ Implements the virtual `Ncm.Data` method `get_length`. @@ -340,10 +615,13 @@ def do_prepare(self, mset: Ncm.MSet): # pylint: disable-msg=arguments-differ self.likelihood.reset() self.tools.reset() - self.mapping.set_params_from_numcosmo(mset) - ccl_args = self.mapping.calculate_ccl_args(mset) - self.ccl_cosmo = ccl.CosmologyCalculator(**self.mapping.asdict(), **ccl_args) - params_map = self.mapping.create_params_map(mset) + self._nc_mapping.set_params_from_numcosmo(mset) + ccl_args = self._nc_mapping.calculate_ccl_args(mset) + + self.ccl_cosmo = ccl.CosmologyCalculator( + **self._nc_mapping.mapping.asdict(), **ccl_args + ) + params_map = self._nc_mapping.create_params_map(self._model_list, mset) self.likelihood.update(params_map) self.tools.prepare(self.ccl_cosmo) @@ -360,6 +638,15 @@ def do_mean_func(self, _, mean_vector): mean_vector.set_array(theory_vector) +# These commands creates GObject types for the defined classes, enabling their use +# within the NumCosmo framework. It is essential to call these functions before +# initializing NumCosmo with the Ncm.init_cfg() function, as failure to do so +# will cause issues with MPI jobs using these objects. +GObject.type_register(MappingNumCosmo) +GObject.type_register(NumCosmoData) +GObject.type_register(NumCosmoGaussCov) + + class NumCosmoFactory: """NumCosmo likelihood class. This class provide the necessary factory methods to create NumCosmo+firecrown likelihoods.""" @@ -369,14 +656,23 @@ def __init__( likelihood_source: str, build_parameters: NamedParameters, mapping: MappingNumCosmo, + model_list: List[str], ): likelihood, tools = load_likelihood(likelihood_source, build_parameters) + likelihood_str = base64.b64encode( + pickle.dumps((likelihood_source, build_parameters)) + ).decode("ascii") + self.mapping: MappingNumCosmo = mapping if isinstance(likelihood, ConstGaussian): - self.data: Ncm.Data = NumCosmoGaussCov(likelihood, tools, mapping) + self.data: Ncm.Data = NumCosmoGaussCov.new_from_likelihood( + likelihood, model_list, tools, mapping, likelihood_str + ) else: - self.data = NumCosmoData(likelihood, tools, mapping) + self.data = NumCosmoData.new_from_likelihood( + likelihood, model_list, tools, mapping, likelihood_str + ) def get_data(self) -> Ncm.Data: """This method return the appropriated Ncm.Data class to be used by NumCosmo.""" diff --git a/tests/conftest.py b/tests/conftest.py index 219fc87b..139c3eb0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -10,6 +10,7 @@ from firecrown.likelihood.gauss_family.statistic.statistic import TrivialStatistic from firecrown.parameters import ParamsMap +from firecrown.connector.mapping import MappingCosmoSIS, mapping_builder def pytest_addoption(parser): @@ -68,3 +69,25 @@ def make_sacc_data(): result.add_data_point("count", (), -3.0) result.add_covariance([4.0, 9.0, 16.0]) return result + + +@pytest.fixture(name="mapping_cosmosis") +def fixture_mapping_cosmosis() -> MappingCosmoSIS: + """Return a MappingCosmoSIS instance.""" + mapping_cosmosis = mapping_builder(input_style="CosmoSIS") + assert isinstance(mapping_cosmosis, MappingCosmoSIS) + mapping_cosmosis.set_params( + Omega_c=0.26, + Omega_b=0.04, + h=0.72, + A_s=2.1e-9, + n_s=0.96, + Omega_k=0.0, + Neff=3.046, + m_nu=0.0, + m_nu_type="normal", + w0=-1.0, + wa=0.0, + T_CMB=2.7255, + ) + return mapping_cosmosis diff --git a/tests/connector/cosmosis/test_cosmosis_module.py b/tests/connector/cosmosis/test_cosmosis_module.py index cb276f2d..c3fbf655 100644 --- a/tests/connector/cosmosis/test_cosmosis_module.py +++ b/tests/connector/cosmosis/test_cosmosis_module.py @@ -261,3 +261,132 @@ def test_module_exec_working( firecrown_mod_with_const_gaussian: FirecrownLikelihood, sample_with_M: DataBlock ): assert firecrown_mod_with_const_gaussian.execute(sample_with_M) == 0 + + +def test_mapping_cosmosis_background(mapping_cosmosis): + block = DataBlock() + + background_d_m = np.geomspace(0.1, 10.0, 100) + background_z = np.linspace(0.0, 2.0, 10) + background_h = np.geomspace(0.1, 10.0, 100) + + block.put_double_array_1d("distances", "d_m", background_d_m) + block.put_double_array_1d("distances", "z", background_z) + block.put_double_array_1d("distances", "h", background_h) + + ccl_args = mapping_cosmosis.calculate_ccl_args(block) + + assert "background" in ccl_args + assert np.allclose( + ccl_args["background"]["a"], + mapping_cosmosis.redshift_to_scale_factor(background_z), + ) + assert np.allclose(ccl_args["background"]["chi"], np.flip(background_d_m)) + assert np.allclose( + ccl_args["background"]["h_over_h0"], + mapping_cosmosis.transform_h_to_h_over_h0(background_h), + ) + + assert "pk_linear" not in ccl_args + assert "pk_nonlin" not in ccl_args + + +def test_mapping_cosmosis_pk_linear(mapping_cosmosis): + block = DataBlock() + + block.put_double_array_1d("distances", "d_m", np.geomspace(0.1, 10.0, 100)) + block.put_double_array_1d("distances", "z", np.linspace(0.0, 2.0, 10)) + block.put_double_array_1d("distances", "h", np.geomspace(0.1, 10.0, 100)) + + matter_power_lin_k_h = np.geomspace(0.1, 10.0, 100) + matter_power_lin_z = np.linspace(0.0, 2.0, 10) + matter_power_lin_p_k = np.geomspace(1.0e-3, 1.0e3, 100) + + block.put_double_array_1d("matter_power_lin", "k_h", matter_power_lin_k_h) + block.put_double_array_1d("matter_power_lin", "z", matter_power_lin_z) + block.put_double_array_1d("matter_power_lin", "p_k", matter_power_lin_p_k) + + ccl_args = mapping_cosmosis.calculate_ccl_args(block) + + assert "background" in ccl_args + assert "pk_linear" in ccl_args + assert "pk_nonlin" not in ccl_args + + assert np.allclose( + ccl_args["pk_linear"]["a"], + mapping_cosmosis.redshift_to_scale_factor(matter_power_lin_z), + ) + assert np.allclose( + ccl_args["pk_linear"]["k"], + mapping_cosmosis.transform_k_h_to_k(matter_power_lin_k_h), + ) + assert np.allclose( + ccl_args["pk_linear"]["delta_matter:delta_matter"], + mapping_cosmosis.redshift_to_scale_factor_p_k( + mapping_cosmosis.transform_p_k_h3_to_p_k(matter_power_lin_p_k) + ), + ) + + +def test_mapping_cosmosis_pk_nonlin(mapping_cosmosis): + block = DataBlock() + + block.put_double_array_1d("distances", "d_m", np.geomspace(0.1, 10.0, 100)) + block.put_double_array_1d("distances", "z", np.linspace(0.0, 2.0, 10)) + block.put_double_array_1d("distances", "h", np.geomspace(0.1, 10.0, 100)) + + matter_power_nl_k_h = np.geomspace(0.1, 10.0, 100) + matter_power_nl_z = np.linspace(0.0, 2.0, 10) + matter_power_nl_p_k = np.geomspace(1.0e-3, 1.0e3, 100) + + block.put_double_array_1d("matter_power_nl", "k_h", matter_power_nl_k_h) + block.put_double_array_1d("matter_power_nl", "z", matter_power_nl_z) + block.put_double_array_1d("matter_power_nl", "p_k", matter_power_nl_p_k) + + ccl_args = mapping_cosmosis.calculate_ccl_args(block) + + assert "background" in ccl_args + assert "pk_linear" not in ccl_args + assert "pk_nonlin" in ccl_args + assert "nonlinear_model" not in ccl_args + + assert np.allclose( + ccl_args["pk_nonlin"]["a"], + mapping_cosmosis.redshift_to_scale_factor(matter_power_nl_z), + ) + assert np.allclose( + ccl_args["pk_nonlin"]["k"], + mapping_cosmosis.transform_k_h_to_k(matter_power_nl_k_h), + ) + assert np.allclose( + ccl_args["pk_nonlin"]["delta_matter:delta_matter"], + mapping_cosmosis.redshift_to_scale_factor_p_k( + mapping_cosmosis.transform_p_k_h3_to_p_k(matter_power_nl_p_k) + ), + ) + + +def test_mapping_cosmosis_pk_nonlin_nonlinear_model(mapping_cosmosis): + block = DataBlock() + + block.put_double_array_1d("distances", "d_m", np.geomspace(0.1, 10.0, 100)) + block.put_double_array_1d("distances", "z", np.linspace(0.0, 2.0, 10)) + block.put_double_array_1d("distances", "h", np.geomspace(0.1, 10.0, 100)) + + mapping_cosmosis.require_nonlinear_pk = True + ccl_args = mapping_cosmosis.calculate_ccl_args(block) + + assert "background" in ccl_args + assert "pk_linear" not in ccl_args + assert "pk_nonlin" not in ccl_args + assert "nonlinear_model" in ccl_args + assert ccl_args["nonlinear_model"] + + mapping_cosmosis.require_nonlinear_pk = False + ccl_args = mapping_cosmosis.calculate_ccl_args(block) + + assert "background" in ccl_args + assert "pk_linear" not in ccl_args + assert "pk_nonlin" not in ccl_args + assert "nonlinear_model" in ccl_args + assert not ccl_args["nonlinear_model"] diff --git a/tests/connector/numcosmo/test_numcosmo_connector.py b/tests/connector/numcosmo/test_numcosmo_connector.py index f362f79c..d3f0d20c 100644 --- a/tests/connector/numcosmo/test_numcosmo_connector.py +++ b/tests/connector/numcosmo/test_numcosmo_connector.py @@ -15,10 +15,12 @@ def test_numcosmo_connector_plain(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) nc_factory = NumCosmoFactory( - "tests/likelihood/lkdir/lkscript.py", NamedParameters(), map_cosmo + "tests/likelihood/lkdir/lkscript.py", + NamedParameters(), + map_cosmo, + model_list=["non_existing_model"], ) assert nc_factory.get_data() is not None @@ -31,12 +33,12 @@ def test_numcosmo_connector_const_gauss(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) nc_factory = NumCosmoFactory( "tests/likelihood/gauss_family/lkscript_const_gaussian.py", NamedParameters(), map_cosmo, + model_list=["non_existing_model"], ) assert nc_factory.get_data() is not None diff --git a/tests/connector/numcosmo/test_numcosmo_mapping.py b/tests/connector/numcosmo/test_numcosmo_mapping.py index 30a35fde..8d494a21 100644 --- a/tests/connector/numcosmo/test_numcosmo_mapping.py +++ b/tests/connector/numcosmo/test_numcosmo_mapping.py @@ -4,12 +4,14 @@ import pyccl as ccl from numcosmo_py import Ncm, Nc +from firecrown.likelihood.likelihood import Likelihood, NamedParameters from firecrown.likelihood.gauss_family.gaussian import ConstGaussian from firecrown.modeling_tools import ModelingTools from firecrown.connector.numcosmo.numcosmo import ( NumCosmoData, NumCosmoGaussCov, MappingNumCosmo, + NumCosmoFactory, ) from firecrown.connector.numcosmo.model import ( @@ -30,7 +32,6 @@ def test_numcosmo_mapping_create_params_map_non_existing_model(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) mset = Ncm.MSet() @@ -40,7 +41,7 @@ def test_numcosmo_mapping_create_params_map_non_existing_model(): RuntimeError, match="Model name non_existing_model was not found in the model set.", ): - map_cosmo.create_params_map(mset) + map_cosmo.create_params_map(["non_existing_model"], mset) def test_numcosmo_mapping_create_params_map_absent_model(): @@ -52,7 +53,6 @@ def test_numcosmo_mapping_create_params_map_absent_model(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["MyModel"], ) mset = Ncm.MSet() @@ -71,7 +71,7 @@ def test_numcosmo_mapping_create_params_map_absent_model(): RuntimeError, match="Model MyModel was not found in the model set.", ): - map_cosmo.create_params_map(mset) + map_cosmo.create_params_map(["MyModel"], mset) def test_numcosmo_mapping_create_params_map_two_models_sharing_parameters(): @@ -83,7 +83,6 @@ def test_numcosmo_mapping_create_params_map_two_models_sharing_parameters(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["MyModel1", "MyModel2"], ) mset = Ncm.MSet() @@ -120,7 +119,7 @@ def test_numcosmo_mapping_create_params_map_two_models_sharing_parameters(): RuntimeError, match="The following keys .* appear in more than one model used by the module", ): - map_cosmo.create_params_map(mset) + map_cosmo.create_params_map(["MyModel1", "MyModel2"], mset) def test_numcosmo_mapping_unsupported(): @@ -131,7 +130,6 @@ def test_numcosmo_mapping_unsupported(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) mset = Ncm.MSet() @@ -149,7 +147,6 @@ def test_numcosmo_mapping_missing_hiprim(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) mset = Ncm.MSet() @@ -171,7 +168,6 @@ def test_numcosmo_mapping_invalid_hiprim(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) mset = Ncm.MSet() @@ -194,7 +190,6 @@ def test_numcosmo_mapping_no_p_mnl_require_nonlinear_pk(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=True, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) mset = Ncm.MSet() @@ -218,7 +213,6 @@ def test_numcosmo_mapping_no_p_mnl(): map_cosmo = MappingNumCosmo( require_nonlinear_pk=False, dist=Nc.Distance.new(6.0), - model_list=["non_existing_model"], ) mset = Ncm.MSet() @@ -245,7 +239,6 @@ def test_numcosmo_mapping(numcosmo_cosmo_fixture, request): p_ml=numcosmo_cosmo["p_ml"], p_mnl=numcosmo_cosmo["p_mnl"], dist=numcosmo_cosmo["dist"], - model_list=["non_existing_model"], ) mset = Ncm.MSet() @@ -253,7 +246,7 @@ def test_numcosmo_mapping(numcosmo_cosmo_fixture, request): map_cosmo.set_params_from_numcosmo(mset) ccl_args = map_cosmo.calculate_ccl_args(mset) - ccl_cosmo = ccl.CosmologyCalculator(**map_cosmo.asdict(), **ccl_args) + ccl_cosmo = ccl.CosmologyCalculator(**map_cosmo.mapping.asdict(), **ccl_args) assert ccl_cosmo["H0"] == cosmo.param_get_by_name("H0") assert ccl_cosmo["Omega_c"] == cosmo.param_get_by_name("Omegac") @@ -261,6 +254,52 @@ def test_numcosmo_mapping(numcosmo_cosmo_fixture, request): assert ccl_cosmo["Omega_k"] == cosmo.param_get_by_name("Omegak") +@pytest.mark.parametrize( + "numcosmo_cosmo_fixture", + ["numcosmo_cosmo_xcdm", "numcosmo_cosmo_xcdm_no_nu", "numcosmo_cosmo_cpl"], +) +def test_numcosmo_serialize_mapping(numcosmo_cosmo_fixture, request): + """Test the NumCosmo mapping connector consistence.""" + numcosmo_cosmo = request.getfixturevalue(numcosmo_cosmo_fixture) + + cosmo = numcosmo_cosmo["cosmo"] + map_cosmo = MappingNumCosmo( + require_nonlinear_pk=True, + p_ml=numcosmo_cosmo["p_ml"], + p_mnl=numcosmo_cosmo["p_mnl"], + dist=numcosmo_cosmo["dist"], + ) + + assert map_cosmo.mapping_name == "default" + + ser = Ncm.Serialize.new(Ncm.SerializeOpt.NONE) + + map_cosmo_dup = ser.dup_obj(map_cosmo) + assert isinstance(map_cosmo_dup, MappingNumCosmo) + assert map_cosmo_dup.mapping_name == "default" + + mset = Ncm.MSet() + mset.set(cosmo) + + map_cosmo.set_params_from_numcosmo(mset) + map_cosmo_dup.set_params_from_numcosmo(mset) + + if map_cosmo_dup.p_ml is None: + assert map_cosmo_dup.p_ml is None + else: + assert id(map_cosmo_dup.p_ml) != id(map_cosmo.p_ml) + assert isinstance(map_cosmo_dup.p_ml, Nc.PowspecML) + + if map_cosmo_dup.p_mnl is None: + assert map_cosmo_dup.p_mnl is None + else: + assert id(map_cosmo_dup.p_mnl) != id(map_cosmo.p_mnl) + assert isinstance(map_cosmo_dup.p_mnl, Nc.PowspecMNL) + + assert id(map_cosmo_dup.dist) != id(map_cosmo.dist) + assert isinstance(map_cosmo_dup.dist, Nc.Distance) + + @pytest.mark.parametrize( "numcosmo_cosmo_fixture", ["numcosmo_cosmo_xcdm", "numcosmo_cosmo_xcdm_no_nu", "numcosmo_cosmo_cpl"], @@ -282,16 +321,16 @@ def test_numcosmo_data( p_ml=numcosmo_cosmo["p_ml"], p_mnl=numcosmo_cosmo["p_mnl"], dist=numcosmo_cosmo["dist"], - model_list=["NcFirecrownTrivial"], ) likelihood = ConstGaussian(statistics=trivial_stats) likelihood.read(sacc_data_for_trivial_stat) - fc_data = NumCosmoData( + fc_data = NumCosmoData.new_from_likelihood( likelihood=likelihood, tools=ModelingTools(), - mapping=map_cosmo, + nc_mapping=map_cosmo, + model_list=["NcFirecrownTrivial"], ) assert fc_data.get_length() > 0 @@ -334,16 +373,16 @@ def test_numcosmo_gauss_cov( p_ml=numcosmo_cosmo["p_ml"], p_mnl=numcosmo_cosmo["p_mnl"], dist=numcosmo_cosmo["dist"], - model_list=["NcFirecrownTrivial"], ) likelihood = ConstGaussian(statistics=trivial_stats) likelihood.read(sacc_data_for_trivial_stat) - fc_data = NumCosmoGaussCov( + fc_data = NumCosmoGaussCov.new_from_likelihood( likelihood=likelihood, tools=ModelingTools(), - mapping=map_cosmo, + nc_mapping=map_cosmo, + model_list=["NcFirecrownTrivial"], ) assert fc_data.get_length() > 0 @@ -363,3 +402,65 @@ def test_numcosmo_gauss_cov( nc_likelihood = Ncm.Likelihood(dataset=dset) assert nc_likelihood.m2lnL_val(mset) == 2.0 + + +@pytest.mark.parametrize( + "numcosmo_cosmo_fixture", + ["numcosmo_cosmo_xcdm", "numcosmo_cosmo_xcdm_no_nu", "numcosmo_cosmo_cpl"], +) +@pytest.mark.parametrize( + "likelihood_file", + [ + "tests/likelihood/lkdir/lkscript.py", + "tests/likelihood/gauss_family/lkscript_const_gaussian.py", + ], +) +def test_numcosmo_serialize_likelihood( + numcosmo_cosmo_fixture, + likelihood_file, + request, +): + """Test the NumCosmo data connector for NcmData with serialization.""" + + numcosmo_cosmo = request.getfixturevalue(numcosmo_cosmo_fixture) + + map_cosmo = MappingNumCosmo( + require_nonlinear_pk=True, + p_ml=numcosmo_cosmo["p_ml"], + p_mnl=numcosmo_cosmo["p_mnl"], + dist=numcosmo_cosmo["dist"], + ) + + nc_factory = NumCosmoFactory( + likelihood_file, + NamedParameters(), + map_cosmo, + model_list=["NcFirecrownTrivial"], + ) + + fc_data = nc_factory.get_data() + + assert isinstance(fc_data, (NumCosmoData, NumCosmoGaussCov)) + assert fc_data.get_length() > 0 + assert fc_data.get_dof() > 0 + + ser = Ncm.Serialize.new(Ncm.SerializeOpt.NONE) + + fc_data_dup = ser.dup_obj(fc_data) + assert isinstance(fc_data_dup, (NumCosmoData, NumCosmoGaussCov)) + + if fc_data.nc_mapping is None: + assert fc_data_dup.nc_mapping is None + else: + assert id(fc_data_dup.nc_mapping) != id(fc_data.nc_mapping) + assert isinstance(fc_data_dup.nc_mapping, MappingNumCosmo) + + assert id(fc_data_dup.likelihood) != id(fc_data.likelihood) + assert isinstance(fc_data_dup.likelihood, Likelihood) + + assert id(fc_data_dup.tools) != id(fc_data.tools) + assert isinstance(fc_data_dup.tools, ModelingTools) + + assert id(fc_data_dup.model_list) != id(fc_data.model_list) + assert isinstance(fc_data_dup.model_list, list) + assert fc_data_dup.model_list == fc_data.model_list diff --git a/tests/connector/test_mapping.py b/tests/connector/test_mapping.py index a816feaf..f14fd53c 100644 --- a/tests/connector/test_mapping.py +++ b/tests/connector/test_mapping.py @@ -1,8 +1,11 @@ """ Unit testing for the mapping module. """ +from typing import Any, Dict import pytest +import numpy as np from firecrown.connector import mapping +from firecrown.connector.mapping import Mapping, mapping_builder from firecrown.likelihood.likelihood import NamedParameters @@ -87,3 +90,103 @@ def test_conversion_from_cosmosis_camb_using_delta_neff(): p = mapping.mapping_builder(input_style="CosmoSIS") p.set_params_from_cosmosis(named_params) assert p.Neff == pytest.approx(3.171) + + +def test_get_params_names(): + fc_map = Mapping() + + with pytest.deprecated_call(): + params_names = fc_map.get_params_names() + assert not params_names + + +def test_transform_k_h_to_k(): + fc_map = Mapping() + + with pytest.deprecated_call(): + fc_map.transform_k_h_to_k([]) + + +def test_transform_p_k_h3_to_p_k(): + fc_map = Mapping() + + with pytest.deprecated_call(): + fc_map.transform_p_k_h3_to_p_k([]) + + +def test_transform_h_to_h_over_h0(): + fc_map = Mapping() + + with pytest.deprecated_call(): + fc_map.transform_h_to_h_over_h0([]) + + +def test_sigma8_and_A_s(): + fc_map = Mapping() + + params_dict: Dict[str, Any] = { + "Omega_c": 0.26, + "Omega_b": 0.04, + "h": 0.72, + "n_s": 0.96, + "Omega_k": 0.0, + "Neff": 3.046, + "m_nu": 0.0, + "m_nu_type": "normal", + "w0": -1.0, + "wa": 0.0, + "T_CMB": 2.7255, + } + + with pytest.raises( + ValueError, match="Exactly one of A_s and sigma8 must be supplied" + ): + fc_map.set_params( + **params_dict, + sigma8=0.8, + A_s=2.1e-9, + ) + + with pytest.raises( + ValueError, match="Exactly one of A_s and sigma8 must be supplied" + ): + fc_map.set_params(**params_dict) + + +def test_mappping_builder(): + with pytest.raises( + ValueError, match="input_style must be .* not invalid_input_style" + ): + mapping_builder(input_style="invalid_input_style") + + +def test_mapping_cosmosis(): + mapping_cosmosis = mapping_builder(input_style="CosmoSIS") + assert isinstance(mapping_cosmosis, Mapping) + + assert mapping_cosmosis.get_params_names() == [ + "h0", + "omega_b", + "omega_c", + "sigma_8", + "n_s", + "omega_k", + "delta_neff", + "omega_nu", + "w", + "wa", + ] + + +def test_mapping_cosmosis_k_h_to_h(mapping_cosmosis): + k_h_array = np.geomspace(0.1, 10.0, 10) + k_array = mapping_cosmosis.transform_k_h_to_k(k_h_array) + + assert np.allclose(k_h_array * mapping_cosmosis.h, k_array) + + +def test_mapping_cosmosis_p_k_h3_to_p_k(mapping_cosmosis): + p_k_h3_array = np.geomspace(0.1, 10.0, 10) + p_k_array = mapping_cosmosis.transform_p_k_h3_to_p_k(p_k_h3_array) + + assert np.allclose(p_k_h3_array / mapping_cosmosis.h**3, p_k_array) diff --git a/tests/likelihood/test_likelihood.py b/tests/likelihood/test_likelihood.py index 6ea10821..bb021ec7 100644 --- a/tests/likelihood/test_likelihood.py +++ b/tests/likelihood/test_likelihood.py @@ -74,10 +74,11 @@ def test_load_likelihood_submodule_returns_wrong_type_tools(): def test_load_likelihood_submodule_old(): dir_path = os.path.dirname(os.path.realpath(__file__)) - load_likelihood( - os.path.join(dir_path, "lkdir/lkscript_old.py"), - NamedParameters(), - ) + with pytest.deprecated_call(): + load_likelihood( + os.path.join(dir_path, "lkdir/lkscript_old.py"), + NamedParameters(), + ) def test_load_likelihood_correct_tools():