From 71f8d34c48441248e4d1bf0619b1679fd773d824 Mon Sep 17 00:00:00 2001 From: Nick Koukoufilippas Date: Tue, 2 Aug 2022 14:44:31 +0300 Subject: [PATCH] Tunable CCL global settings at the Python-level (#918) * first commit * ParamStruct methods keys, values, items * enable flake in constants.py * got rid of F401 in a Pythonic manner * coverage * reduced size of ParamStruct object * coverage * public public --> _public private * support immutable values; add physical constants * relax capitalized-only requirement * flake * coverage * added copy method * moved in separate module * expose CCLParameters * removed requirement for gsl_params and spline_params to be const at the C-level; makes the entire PR shorter and simpler * fix relative import * docstrings, frozen keys, reload option, tests * don't pickle swig spline types * singleton implementation moved under base class * implement get global na/nk sampling without cosmo; init Pk2D from pkfunc without cosmo * removed unnecessary test * new docstrings * option to call new funcs with or without cosmo * fix docs * updated warning in tracers.py * some minor improvements * updated readthedocs with new usage * fixed F401 * replace all swig lookups with python-level * clear up namespace & allow T_CMB changing * removed every trace of direct parameter assignment from code base * Refactor: direct interface with SWIG; deprecate direct assignment; no python-level copy; no singleton. * removed unecessary import * save spline and gsl params in cosmo * small improvement * tests to improve coverage * E303 * remove _type attribute * removed direct parameter assignment from benchmarks * reload frozen parameters (physical constants) * correct repr * verbose warning & added a test * no need to pass if docstring is provided * deprecation message Co-authored-by: David Alonso --- benchmarks/test_cls.py | 6 +- benchmarks/test_correlation.py | 6 +- benchmarks/test_correlation_MG.py | 6 +- benchmarks/test_correlation_MG2.py | 6 +- benchmarks/test_correlation_MG3_SD.py | 6 +- benchmarks/test_tracers.py | 8 +- include/ccl_core.h | 6 +- pyccl/__init__.py | 16 ++- pyccl/background.py | 2 +- pyccl/boltzmann.py | 5 +- pyccl/ccl_pk2d.i | 10 ++ pyccl/core.py | 16 ++- pyccl/errors.py | 23 ++++ pyccl/halos/halo_model.py | 3 +- pyccl/halos/hmfunc.py | 3 +- pyccl/neutrinos.py | 5 +- pyccl/nl_pt/tracers.py | 4 +- pyccl/parameters.py | 121 ++++++++++++++++++ pyccl/pk2d.py | 33 ++--- pyccl/pyutils.py | 62 +++++++++ pyccl/tests/test_cclerror.py | 29 +++++ pyccl/tests/test_cosmology.py | 77 +++++++++++ pyccl/tests/test_gsl_splines.py | 13 -- pyccl/tests/test_pk2d.py | 10 +- pyccl/tests/test_pt_power.py | 5 +- pyccl/tests/test_tracers.py | 29 +++-- pyccl/tracers.py | 15 ++- ...ion_and_other_cosmological_conventions.rst | 19 ++- src/ccl_core.c | 41 +++++- 29 files changed, 479 insertions(+), 106 deletions(-) create mode 100644 pyccl/parameters.py diff --git a/benchmarks/test_cls.py b/benchmarks/test_cls.py index 11d0f60bd..27606f3d3 100644 --- a/benchmarks/test_cls.py +++ b/benchmarks/test_cls.py @@ -11,13 +11,13 @@ def set_up(request): nztyp = request.param dirdat = os.path.dirname(__file__) + '/data/' + ccl.gsl_params.INTEGRATION_LIMBER_EPSREL = 1E-4 + ccl.gsl_params.INTEGRATION_EPSREL = 1E-4 cosmo = ccl.Cosmology(Omega_c=0.30, Omega_b=0.00, Omega_g=0, Omega_k=0, h=0.7, sigma8=0.8, n_s=0.96, Neff=0, m_nu=0.0, w0=-1, wa=0, T_CMB=2.7, transfer_function='bbks', mass_function='tinker', matter_power_spectrum='linear') - cosmo.cosmo.gsl_params.INTEGRATION_LIMBER_EPSREL = 1E-4 - cosmo.cosmo.gsl_params.INTEGRATION_EPSREL = 1E-4 # ell-arrays nls = 541 @@ -118,6 +118,8 @@ def read_bm(fname): bms['cc'] = read_bm(pre + 'log_cl_cc.txt') print('init and i/o time:', time.time() - t0) + ccl.gsl_params.reload() + return cosmo, trc, lfacs, bms diff --git a/benchmarks/test_correlation.py b/benchmarks/test_correlation.py index be4d7cbcf..27a988f44 100644 --- a/benchmarks/test_correlation.py +++ b/benchmarks/test_correlation.py @@ -20,13 +20,13 @@ def set_up(request): t0 = time.time() nztyp = request.param dirdat = os.path.dirname(__file__) + '/data/' + ccl.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 + ccl.gsl_params.INTEGRATION_EPSREL = 2.5E-5 cosmo = ccl.Cosmology(Omega_c=0.30, Omega_b=0.00, Omega_g=0, Omega_k=0, h=0.7, sigma8=0.8, n_s=0.96, Neff=0, m_nu=0.0, w0=-1, wa=0, T_CMB=2.7, transfer_function='bbks', mass_function='tinker', matter_power_spectrum='linear') - cosmo.cosmo.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 - cosmo.cosmo.gsl_params.INTEGRATION_EPSREL = 2.5E-5 # Ell-dependent correction factors # Set up array of ells @@ -173,6 +173,8 @@ def read_bm(fname): fill_value=d[3][0], bounds_error=False)(theta) print('setup time:', time.time() - t0) + + ccl.gsl_params.reload() return cosmo, trc, bms, ers, fl diff --git a/benchmarks/test_correlation_MG.py b/benchmarks/test_correlation_MG.py index edb156bd3..34ba6d5b8 100644 --- a/benchmarks/test_correlation_MG.py +++ b/benchmarks/test_correlation_MG.py @@ -16,14 +16,14 @@ def set_up(request): dirdat = os.path.dirname(__file__) + '/data/' h0 = 0.70001831054687500 logA = 3.05 # log(10^10 A_s) + ccl.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 + ccl.gsl_params.INTEGRATION_EPSREL = 2.5E-5 cosmo = ccl.Cosmology(Omega_c=0.12/h0**2, Omega_b=0.0221/h0**2, Omega_k=0, h=h0, A_s=np.exp(logA)/10**10, n_s=0.96, Neff=3.046, m_nu=0.0, w0=-1, wa=0, T_CMB=2.7255, mu_0=0.1, sigma_0=0.1, transfer_function='boltzmann_class', matter_power_spectrum='linear') - cosmo.cosmo.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 - cosmo.cosmo.gsl_params.INTEGRATION_EPSREL = 2.5E-5 # Ell-dependent correction factors # Set up array of ells @@ -118,6 +118,8 @@ def set_up(request): ers['ll_12_m'] = interp1d(d[0], d[3], fill_value=d[3][0], bounds_error=False)(theta) + + ccl.gsl_params.reload() return cosmo, trc, bms, ers, fl diff --git a/benchmarks/test_correlation_MG2.py b/benchmarks/test_correlation_MG2.py index 43a054332..44f0bdf88 100644 --- a/benchmarks/test_correlation_MG2.py +++ b/benchmarks/test_correlation_MG2.py @@ -16,14 +16,14 @@ def set_up(request): dirdat = os.path.dirname(__file__) + '/data/' h0 = 0.67702026367187500 logA = 3.05 # log(10^10 A_s) + ccl.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 + ccl.gsl_params.INTEGRATION_EPSREL = 2.5E-5 cosmo = ccl.Cosmology(Omega_c=0.12/h0**2, Omega_b=0.0221/h0**2, Omega_k=0, h=h0, A_s=np.exp(logA)/10**10, n_s=0.96, Neff=3.046, m_nu=0.0, w0=-1, wa=0, T_CMB=2.7255, mu_0=0.1, sigma_0=0.1, transfer_function='boltzmann_isitgr', matter_power_spectrum='linear') - cosmo.cosmo.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 - cosmo.cosmo.gsl_params.INTEGRATION_EPSREL = 2.5E-5 # Ell-dependent correction factors # Set up array of ells @@ -118,6 +118,8 @@ def set_up(request): ers['ll_12_m'] = interp1d(d[0], d[3], fill_value=d[3][0], bounds_error=False)(theta) + + ccl.gsl_params.reload() return cosmo, trc, bms, ers, fl diff --git a/benchmarks/test_correlation_MG3_SD.py b/benchmarks/test_correlation_MG3_SD.py index 65c3138a2..178c7e21b 100644 --- a/benchmarks/test_correlation_MG3_SD.py +++ b/benchmarks/test_correlation_MG3_SD.py @@ -17,6 +17,8 @@ def set_up(request): h0 = 0.67702026367187500 logA = 3.05 # log(10^10 A_s) # scale dependent MG cosmology + ccl.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 + ccl.gsl_params.INTEGRATION_EPSREL = 2.5E-5 cosmo = ccl.Cosmology(Omega_c=0.12/h0**2, Omega_b=0.0221/h0**2, Omega_k=0, h=h0, A_s=np.exp(logA)/10**10, n_s=0.96, Neff=3.046, m_nu=0.0, w0=-1, wa=0, T_CMB=2.7255, @@ -24,8 +26,6 @@ def set_up(request): c1_mg=1.1, c2_mg=1.1, lambda_mg=1, transfer_function='boltzmann_isitgr', matter_power_spectrum='linear') - cosmo.cosmo.gsl_params.INTEGRATION_LIMBER_EPSREL = 2.5E-5 - cosmo.cosmo.gsl_params.INTEGRATION_EPSREL = 2.5E-5 # Ell-dependent correction factors # Set up array of ells @@ -120,6 +120,8 @@ def set_up(request): ers['ll_12_m'] = interp1d(d[0], d[3], fill_value=d[3][0], bounds_error=False)(theta) + + ccl.gsl_params.reload() return cosmo, trc, bms, ers, fl diff --git a/benchmarks/test_tracers.py b/benchmarks/test_tracers.py index 6f0fd03b7..1bb466ea4 100644 --- a/benchmarks/test_tracers.py +++ b/benchmarks/test_tracers.py @@ -20,15 +20,17 @@ def get_prediction(ells, chi_i, chi_f, alpha, beta, gamma, @pytest.fixture(scope='module') def set_up(): + ccl.physical_constants.T_CMB = 2.7 + ccl.gsl_params.INTEGRATION_LIMBER_EPSREL = 1E-4 + ccl.gsl_params.INTEGRATION_EPSREL = 1E-4 cosmo = ccl.Cosmology(Omega_c=0.30, Omega_b=0.00, Omega_g=0, Omega_k=0, h=0.7, sigma8=0.8, n_s=0.96, Neff=0, m_nu=0.0, w0=-1, wa=0, transfer_function='bbks', mass_function='tinker', matter_power_spectrum='linear') - cosmo.cosmo.params.T_CMB = 2.7 - cosmo.cosmo.gsl_params.INTEGRATION_LIMBER_EPSREL = 1E-4 - cosmo.cosmo.gsl_params.INTEGRATION_EPSREL = 1E-4 + ccl.physical_constants.reload() + ccl.gsl_params.reload() return cosmo diff --git a/include/ccl_core.h b/include/ccl_core.h index 59afd66e6..7a09d65b5 100644 --- a/include/ccl_core.h +++ b/include/ccl_core.h @@ -149,7 +149,7 @@ typedef struct ccl_spline_params { const gsl_interp_type* CORR_SPLINE_TYPE; } ccl_spline_params; -extern const ccl_spline_params default_spline_params; +extern ccl_spline_params ccl_user_spline_params; /** * Struct that contains parameters that control the accuracy of various GSL @@ -195,7 +195,7 @@ typedef struct ccl_gsl_params { bool LENSING_KERNEL_SPLINE_INTEGRATION; } ccl_gsl_params; -extern const ccl_gsl_params default_gsl_params; +extern ccl_gsl_params ccl_user_gsl_params; /** * Struct containing the parameters defining a cosmology @@ -362,6 +362,8 @@ int ccl_get_pk_spline_na(ccl_cosmology *cosmo); int ccl_get_pk_spline_nk(ccl_cosmology *cosmo); void ccl_get_pk_spline_a_array(ccl_cosmology *cosmo,int ndout,double* doutput,int *status); void ccl_get_pk_spline_lk_array(ccl_cosmology *cosmo,int ndout,double* doutput,int *status); +void ccl_get_pk_spline_a_array_from_params(ccl_spline_params *spline_params, int ndout, double *doutput, int *status); +void ccl_get_pk_spline_lk_array_from_params(ccl_spline_params *spline_params, int ndout, double *doutput, int *status); CCL_END_DECLS diff --git a/pyccl/__init__.py b/pyccl/__init__.py index d7f7d5c25..b57a89545 100644 --- a/pyccl/__init__.py +++ b/pyccl/__init__.py @@ -20,6 +20,15 @@ from .errors import ( CCLError, CCLWarning, + CCLDeprecationWarning, +) + +# Constants and accuracy parameters +from .parameters import ( + CCLParameters, + gsl_params, + spline_params, + physical_constants, ) # Core data structures @@ -118,12 +127,11 @@ sigma2_B_from_mask, ) -# Parameters -physical_constants = lib.cvar.constants # Miscellaneous from .pyutils import debug_mode, resample_array + # Deprecated & Renamed modules from .halomodel import ( halomodel_matter_power, @@ -148,8 +156,8 @@ __all__ = ( 'lib', - 'physical_constants', - 'CCLError', 'CCLWarning', + 'CCLParameters', 'spline_params', 'gsl_params', 'physical_constants', + 'CCLError', 'CCLWarning', 'CCLDeprecationWarning', 'Cosmology', 'CosmologyVanillaLCDM', 'CosmologyCalculator', 'growth_factor', 'growth_factor_unnorm', 'growth_rate', 'comoving_radial_distance', 'angular_diameter_distance', diff --git a/pyccl/background.py b/pyccl/background.py index d36bd56cb..e715061f2 100644 --- a/pyccl/background.py +++ b/pyccl/background.py @@ -16,6 +16,7 @@ from . import ccllib as lib from .pyutils import _vectorize_fn, _vectorize_fn3 from .pyutils import _vectorize_fn4, _vectorize_fn5 +from .parameters import physical_constants species_types = { 'critical': lib.species_crit_label, @@ -299,7 +300,6 @@ def sigma_critical(cosmo, a_lens, a_source): float or array_like: :math:`\\Sigma_{\\mathrm{crit}}` in units of :math:`\\M_{\\odot}/Mpc^2` """ - physical_constants = lib.cvar.constants Ds = angular_diameter_distance(cosmo, a_source, a2=None) Dl = angular_diameter_distance(cosmo, a_lens, a2=None) Dls = angular_diameter_distance(cosmo, a_lens, a_source) diff --git a/pyccl/boltzmann.py b/pyccl/boltzmann.py index a7854288a..7c32339b9 100644 --- a/pyccl/boltzmann.py +++ b/pyccl/boltzmann.py @@ -9,6 +9,7 @@ from .pyutils import check from .pk2d import Pk2D from .errors import CCLError +from .parameters import physical_constants def get_camb_pk_lin(cosmo, nonlin=False): @@ -101,7 +102,7 @@ def get_camb_pk_lin(cosmo, nonlin=False): # thus we set T_i_eff = T_i = g^(1/4) * T_nu and solve for the right # value of g for CAMB. We get g = (TNCDM / (11/4)^(-1/3))^4 g = np.power( - lib.cvar.constants.TNCDM / np.power(11.0/4.0, -1.0/3.0), + physical_constants.TNCDM / np.power(11.0/4.0, -1.0/3.0), 4.0) if cosmo['N_nu_mass'] > 0: @@ -323,7 +324,7 @@ def get_isitgr_pk_lin(cosmo): # thus we set T_i_eff = T_i = g^(1/4) * T_nu and solve for the right # value of g for CAMB. We get g = (TNCDM / (11/4)^(-1/3))^4 g = np.power( - lib.cvar.constants.TNCDM / np.power(11.0/4.0, -1.0/3.0), + physical_constants.TNCDM / np.power(11.0/4.0, -1.0/3.0), 4.0) if cosmo['N_nu_mass'] > 0: diff --git a/pyccl/ccl_pk2d.i b/pyccl/ccl_pk2d.i index 74b45aac2..4ada1750b 100644 --- a/pyccl/ccl_pk2d.i +++ b/pyccl/ccl_pk2d.i @@ -32,11 +32,21 @@ void get_pk_spline_a(ccl_cosmology *cosmo,int ndout,double* doutput,int *status) ccl_get_pk_spline_a_array(cosmo,ndout,doutput,status); } +void get_pk_spline_a_from_params(ccl_spline_params *spline_params, int ndout, double *doutput, int *status) +{ + ccl_get_pk_spline_a_array_from_params(spline_params, ndout, doutput, status); +} + void get_pk_spline_lk(ccl_cosmology *cosmo,int ndout,double* doutput,int *status) { ccl_get_pk_spline_lk_array(cosmo,ndout,doutput,status); } +void get_pk_spline_lk_from_params(ccl_spline_params *spline_params, int ndout, double *doutput, int *status) +{ + ccl_get_pk_spline_lk_array_from_params(spline_params, ndout, doutput, status); +} + double pk2d_eval_single(ccl_f2d_t *psp,double lk,double a,ccl_cosmology *cosmo,int *status) { return ccl_f2d_t_eval(psp,lk,a,cosmo,status); diff --git a/pyccl/core.py b/pyccl/core.py index 00310dfe3..a8c272dd2 100644 --- a/pyccl/core.py +++ b/pyccl/core.py @@ -14,6 +14,7 @@ from .pyutils import check from .pk2d import Pk2D from .bcm import bcm_correct_pk2d +from .parameters import CCLParameters, physical_constants # Configuration types transfer_function_types = { @@ -200,9 +201,9 @@ class Cosmology(object): # an attribute of this class. from . import (background, bcm, boltzmann, cls, correlations, covariances, neutrinos, - pk2d, power, tk3d, tracers, halos, nl_pt) + pk2d, power, pyutils, tk3d, tracers, halos, nl_pt) subs = [background, boltzmann, bcm, cls, correlations, covariances, - neutrinos, pk2d, power, tk3d, tracers, halos, nl_pt] + neutrinos, pk2d, power, pyutils, tk3d, tracers, halos, nl_pt] funcs = [getmembers(sub, isfunction) for sub in subs] funcs = [func for sub in funcs for func in sub] for name, func in funcs: @@ -211,7 +212,7 @@ class Cosmology(object): vars()[name] = func # clear unnecessary locals del (background, boltzmann, bcm, cls, correlations, covariances, - neutrinos, pk2d, power, tk3d, tracers, halos, nl_pt, + neutrinos, pk2d, power, pyutils, tk3d, tracers, halos, nl_pt, subs, funcs, func, name, pars) def __init__( @@ -263,6 +264,9 @@ def _build_cosmo(self): self._build_parameters(**self._params_init_kwargs) self._build_config(**self._config_init_kwargs) self.cosmo = lib.cosmology_create(self._params, self._config) + self._spline_params = CCLParameters.get_params_dict("spline_params") + self._gsl_params = CCLParameters.get_params_dict("gsl_params") + self._accuracy_params = {**self._spline_params, **self._gsl_params} if self.cosmo.status != 0: raise CCLError( @@ -584,10 +588,10 @@ def _build_parameters( # Create new instance of ccl_parameters object # Create an internal status variable; needed to check massive neutrino # integral. - T_CMB_old = lib.cvar.constants.T_CMB + T_CMB_old = physical_constants.T_CMB try: if T_CMB is not None: - lib.cvar.constants.T_CMB = T_CMB + physical_constants.T_CMB = T_CMB status = 0 if nz_mg == -1: # Create ccl_parameters without modified growth @@ -605,7 +609,7 @@ def _build_parameters( df_mg, mnu_final_list, status) check(status) finally: - lib.cvar.constants.T_CMB = T_CMB_old + physical_constants.T_CMB = T_CMB_old if Omega_g is not None: total = self._params.Omega_g + self._params.Omega_l diff --git a/pyccl/errors.py b/pyccl/errors.py index 2a467df69..dba5ef8c5 100644 --- a/pyccl/errors.py +++ b/pyccl/errors.py @@ -1,3 +1,6 @@ +import warnings + + class CCLError(RuntimeError): """A CCL-specific RuntimeError""" def __repr__(self): @@ -20,3 +23,23 @@ def __eq__(self, other): def __hash__(self): return hash(repr(self)) + + +class CCLDeprecationWarning(FutureWarning): + """A CCL-specific deprecation warning.""" + def __repr__(self): + return 'pyccl.CCLDeprecationWarning(%r)' % (str(self)) + + def __eq__(self, other): + return repr(self) == repr(other) + + def __hash__(self): + return hash(repr(self)) + + @classmethod + def enable(cls): + warnings.simplefilter("always") + + @classmethod + def disable(cls): + warnings.filterwarnings(action="ignore", category=cls) diff --git a/pyccl/halos/halo_model.py b/pyccl/halos/halo_model.py index 0f1d04e7a..bc539b16e 100644 --- a/pyccl/halos/halo_model.py +++ b/pyccl/halos/halo_model.py @@ -12,10 +12,9 @@ from ..pyutils import _spline_integrate from .. import background from ..errors import CCLWarning +from ..parameters import physical_constants import numpy as np -physical_constants = lib.cvar.constants - class HMCalculator(object): """ This class implements a set of methods that can be used to diff --git a/pyccl/halos/hmfunc.py b/pyccl/halos/hmfunc.py index 00fe18a29..b9441053e 100644 --- a/pyccl/halos/hmfunc.py +++ b/pyccl/halos/hmfunc.py @@ -1,6 +1,7 @@ from .. import ccllib as lib from ..core import check from ..background import omega_x +from ..parameters import physical_constants from .massdef import MassDef, MassDef200m import numpy as np import functools @@ -149,7 +150,7 @@ def get_mass_function(self, cosmo, M, a, mdef_other=None): len(logM), status) check(status, cosmo=cosmo) - rho = (lib.cvar.constants.RHO_CRITICAL * + rho = (physical_constants.RHO_CRITICAL * cosmo['Omega_m'] * cosmo['h']**2) f = self._get_fsigma(cosmo, sigM, a, 2.302585092994046 * logM) mf = f * rho * dlns_dlogM / M_use diff --git a/pyccl/neutrinos.py b/pyccl/neutrinos.py index df0743c54..b55687857 100644 --- a/pyccl/neutrinos.py +++ b/pyccl/neutrinos.py @@ -1,6 +1,7 @@ import numpy as np from . import ccllib as lib from .core import check +from .parameters import physical_constants neutrino_mass_splits = { 'normal': lib.nu_normal, @@ -28,7 +29,7 @@ def Omeganuh2(a, m_nu, T_CMB=None): scalar = True if np.ndim(a) == 0 else False if T_CMB is None: - T_CMB = lib.cvar.constants.T_CMB + T_CMB = physical_constants.T_CMB # Convert to array if it's not already an array if not isinstance(a, np.ndarray): @@ -65,7 +66,7 @@ def nu_masses(OmNuh2, mass_split, T_CMB=None): status = 0 if T_CMB is None: - T_CMB = lib.cvar.constants.T_CMB + T_CMB = physical_constants.T_CMB if mass_split not in neutrino_mass_splits.keys(): raise ValueError( diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index 065da97fa..419fb6867 100644 --- a/pyccl/nl_pt/tracers.py +++ b/pyccl/nl_pt/tracers.py @@ -2,7 +2,7 @@ from scipy.interpolate import interp1d from ..pyutils import _check_array_params from ..background import growth_factor -from .. import ccllib as lib +from ..parameters import physical_constants def translate_IA_norm(cosmo, z, a1=1.0, a1delta=None, a2=None, @@ -50,7 +50,7 @@ def check_input_array(a, name): check_input_array(a1delta, 'a1delta') Om_m = cosmo['Omega_m'] - rho_crit = lib.cvar.constants.RHO_CRITICAL + rho_crit = physical_constants.RHO_CRITICAL c1 = c1delta = c2 = None gz = growth_factor(cosmo, 1./(1+z)) diff --git a/pyccl/parameters.py b/pyccl/parameters.py new file mode 100644 index 000000000..32df8db4a --- /dev/null +++ b/pyccl/parameters.py @@ -0,0 +1,121 @@ +from . import ccllib as lib +from .errors import warnings, CCLDeprecationWarning + + +class CCLParameters: + """Base for classes holding global CCL parameters and their values. + + Subclasses contain a reference to the C-struct with the collection + of parameters and their values (via SWIG). All subclasses act as proxies + to the CCL parameters at the C-level. + + Subclasses automatically store a backup of the initial parameter state + to enable ad hoc reloading. + """ + + def __init_subclass__(cls, instance=None, freeze=False): + """Routine for subclass initialization. + + Parameters: + instance (``pyccl.ccllib.cvar``): + Reference to ``cvar`` where the parameters are implemented. + freeze (``bool``): + Disable parameter mutation. + """ + super().__init_subclass__() + cls._instance = instance + cls._frozen = freeze + + def _new_setattr(self, key, value): + # Make instances of the SWIG-level class immutable + # so that everything is handled through this interface. + # SWIG only assigns `this` via the low level `_ccllib`; + # we therefore disable all other direct assignments. + if key == "this": + return object.__setattr__(self, key, value) + name = self.__class__.__name__ + # TODO: Deprecation cycle for fully immutable Cosmology objects. + # raise AttributeError(f"Direct assignment in {name} not supported.") # noqa + warnings.warn( + f"Direct assignment of {name} is deprecated " + "and an error will be raised in the next CCL release." + f"Set via `pyccl.{name}.{key}` before instantiation.", + CCLDeprecationWarning) + object.__setattr__(self, key, value) + + cls._instance.__class__.__setattr__ = _new_setattr + + def __init__(self): + # Keep a copy of the default parameters. + object.__setattr__(self, "_bak", CCLParameters.get_params_dict(self)) + + def __getattribute__(self, name): + get = object.__getattribute__ + try: + return get(get(self, "_instance"), name) + except AttributeError: + return get(self, name) + + def __setattr__(self, key, value): + if self._frozen and key != "T_CMB": + # `T_CMB` mutates in Cosmology. + name = self.__class__.__name__ + raise AttributeError(f"Instances of {name} are frozen.") + if not hasattr(self._instance, key): + raise KeyError(f"Parameter {key} does not exist.") + if (key, value) == ("A_SPLINE_MAX", 1.0): + # Setting `A_SPLINE_MAX` to its default value; do nothing. + return + object.__setattr__(self._instance, key, value) + + __getitem__ = __getattribute__ + + __setitem__ = __setattr__ + + def __repr__(self): + out = self._bak.copy() + for par in out: + out[par] = getattr(self, par) + return repr(out) + + def reload(self): + """Reload the C-level default CCL parameters.""" + frozen = self._frozen + if frozen: + object.__setattr__(self, "_frozen", False) + for param, value in self._bak.items(): + setattr(self, param, value) + object.__setattr__(self, "_frozen", frozen) + + @classmethod + def get_params_dict(cls, name): + """Get a dictionary of the current parameters. + + Arguments: + name (str or :obj:`CCLParameters`): + Name or instance of the parameters to look up. + """ + pars = eval(name) if isinstance(name, str) else name + out = {} + for par in dir(pars): + if not par.startswith("_") and par not in ["this", "thisown"]: + out[par] = getattr(pars, par) + return out + + +class SplineParams(CCLParameters, instance=lib.cvar.user_spline_params): + """Instances of this class hold the spline parameters.""" + + +class GSLParams(CCLParameters, instance=lib.cvar.user_gsl_params): + """Instances of this class hold the gsl parameters.""" + + +class PhysicalConstants(CCLParameters, instance=lib.cvar.constants, + freeze=True): + """Instances of this class hold the physical constants.""" + + +spline_params = SplineParams() +gsl_params = GSLParams() +physical_constants = PhysicalConstants() diff --git a/pyccl/pk2d.py b/pyccl/pk2d.py index cbf5b92ff..603eb2737 100644 --- a/pyccl/pk2d.py +++ b/pyccl/pk2d.py @@ -1,11 +1,10 @@ import warnings - import numpy as np from . import ccllib as lib - -from .errors import CCLError, CCLWarning -from .pyutils import check, _get_spline1d_arrays, _get_spline2d_arrays +from .errors import CCLWarning, CCLError +from .pyutils import (check, get_pk_spline_a, get_pk_spline_lk, + _get_spline1d_arrays, _get_spline2d_arrays) class Pk2D(object): @@ -56,10 +55,8 @@ class Pk2D(object): logarithm of the power spectrum. Otherwise, the true value of the power spectrum is expected. Note that arrays will be interpolated in log space if `is_logp` is set to `True`. - cosmo (:class:`~pyccl.core.Cosmology`): Cosmology object. The cosmology - object is needed in order if `pkfunc` is not `None`. The object is - used to determine the sampling rate in scale factor and - wavenumber. + cosmo (:class:`~pyccl.core.Cosmology`, optional): Cosmology object. + Used to determine sampling rates in scale factor and wavenumber. empty (bool): if True, just create an empty object, to be filled out later """ @@ -70,7 +67,6 @@ def __init__(self, pkfunc=None, a_arr=None, lk_arr=None, pk_arr=None, self.has_psp = False return - status = 0 if pkfunc is None: # Initialize power spectrum from 2D array # Make sure input makes sense if (a_arr is None) or (lk_arr is None) or (pk_arr is None): @@ -78,7 +74,7 @@ def __init__(self, pkfunc=None, a_arr=None, lk_arr=None, pk_arr=None, "you must provide arrays") # Check that `a` is a monotonically increasing array. - if not np.array_equal(a_arr, np.sort(a_arr)): + if not np.all((a_arr[1:] - a_arr[:-1]) > 0): raise ValueError("Input scale factor array in `a_arr` is not " "monotonically increasing.") @@ -93,27 +89,18 @@ def __init__(self, pkfunc=None, a_arr=None, lk_arr=None, pk_arr=None, except Exception: raise ValueError("Can't use input function") - if cosmo is None: - raise ValueError("A cosmology is needed if initializing " - "power spectrum from a function") - # Set k and a sampling from CCL parameters - nk = lib.get_pk_spline_nk(cosmo.cosmo) - na = lib.get_pk_spline_na(cosmo.cosmo) - a_arr, status = lib.get_pk_spline_a(cosmo.cosmo, na, status) - check(status, cosmo=cosmo) - lk_arr, status = lib.get_pk_spline_lk(cosmo.cosmo, nk, status) - check(status, cosmo=cosmo) + a_arr = get_pk_spline_a(cosmo=cosmo) + lk_arr = get_pk_spline_lk(cosmo=cosmo) # Compute power spectrum on 2D grid - pkflat = np.zeros([na, nk]) - for ia, a in enumerate(a_arr): - pkflat[ia, :] = pkfunc(k=np.exp(lk_arr), a=a) + pkflat = np.array([pkfunc(k=np.exp(lk_arr), a=a) for a in a_arr]) pkflat = pkflat.flatten() self.extrap_order_lok = extrap_order_lok self.extrap_order_hik = extrap_order_hik + status = 0 self.psp, status = lib.set_pk2d_new_from_arrays(lk_arr, a_arr, pkflat, int(extrap_order_lok), int(extrap_order_hik), diff --git a/pyccl/pyutils.py b/pyccl/pyutils.py index f76b129ab..fd0c0bc2c 100644 --- a/pyccl/pyutils.py +++ b/pyccl/pyutils.py @@ -2,6 +2,7 @@ well as wrappers to automatically vectorize functions.""" from . import ccllib as lib from ._types import error_types +from .parameters import spline_params from .errors import CCLError, CCLWarning import functools import warnings @@ -351,6 +352,67 @@ def _vectorize_fn6(fn, fn_vec, cosmo, x1, x2, returns_status=True): return f +def get_pk_spline_nk(cosmo=None): + """Get the number of sampling points in the wavenumber dimension. + + Arguments: + cosmo (``~pyccl.ccllib.cosmology`` via SWIG, optional): + Input cosmology. + """ + if cosmo is not None: + return lib.get_pk_spline_nk(cosmo.cosmo) + ndecades = np.log10(spline_params.K_MAX / spline_params.K_MIN) + return int(np.ceil(ndecades*spline_params.N_K)) + + +def get_pk_spline_na(cosmo=None): + """Get the number of sampling points in the scale factor dimension. + + Arguments: + cosmo (``~pyccl.ccllib.cosmology`` via SWIG, optional): + Input cosmology. + """ + if cosmo is not None: + return lib.get_pk_spline_na(cosmo.cosmo) + return spline_params.A_SPLINE_NA_PK + spline_params.A_SPLINE_NLOG_PK - 1 + + +def get_pk_spline_lk(cosmo=None): + """Get a log(k)-array with sampling rate defined by ``ccl.spline_params`` + or by the spline parameters of the input ``cosmo``. + + Arguments: + cosmo (``~pyccl.ccllib.cosmology`` via SWIG, optional): + Input cosmology. + """ + nk = get_pk_spline_nk(cosmo=cosmo) + if cosmo is not None: + lk_arr, status = lib.get_pk_spline_lk(cosmo.cosmo, nk, 0) + check(status, cosmo) + return lk_arr + lk_arr, status = lib.get_pk_spline_lk_from_params(spline_params, nk, 0) + check(status) + return lk_arr + + +def get_pk_spline_a(cosmo=None): + """Get an a-array with sampling rate defined by ``ccl.spline_params`` + or by the spline parameters of the input ``cosmo``. + + Arguments: + cosmo (``~pyccl.ccllib.cosmology`` via SWIG, optional): + Input cosmology. + """ + na = get_pk_spline_na(cosmo=cosmo) + if cosmo is not None: + a_arr, status = lib.get_pk_spline_a(cosmo.cosmo, na, 0) + check(status, cosmo) + return a_arr + a_arr, status = lib.get_pk_spline_a_from_params(spline_params, na, 0) + check(status) + return a_arr + + def resample_array(x_in, y_in, x_out, extrap_lo='none', extrap_hi='none', fill_value_lo=0, fill_value_hi=0): diff --git a/pyccl/tests/test_cclerror.py b/pyccl/tests/test_cclerror.py index fca894bfd..d01b73781 100644 --- a/pyccl/tests/test_cclerror.py +++ b/pyccl/tests/test_cclerror.py @@ -26,6 +26,11 @@ def test_cclwarning_repr(): assert_(str(w2) == str(w)) assert_(w2 == w) + v = pyccl.CCLDeprecationWarning("blah") + v2 = eval(repr(v)) + assert str(v2) == str(v) + assert v2 == v + def test_cclwarning_not_equal(): """Check that a CCLWarning can be built from its repr""" @@ -34,3 +39,27 @@ def test_cclwarning_not_equal(): assert_(w is not w2) assert_(w != w2) assert_(hash(w) != hash(w2)) + + v = pyccl.CCLDeprecationWarning("blah") + v2 = pyccl.CCLDeprecationWarning("blahh") + assert v is not v2 + assert v != v2 + assert hash(v) != hash(v2) + + +def test_ccl_deprecationwarning_switch(): + import warnings + + # check that warnings are enabled by default + with warnings.catch_warnings(record=True) as w: + warnings.warn("test", pyccl.CCLDeprecationWarning) + assert w[0].category == pyccl.CCLDeprecationWarning + + # switch off CCL (future) deprecation warnings + pyccl.CCLDeprecationWarning.disable() + with warnings.catch_warnings(record=True) as w: + warnings.warn("test", pyccl.CCLDeprecationWarning) + assert len(w) == 0 + + # switch back on + pyccl.CCLDeprecationWarning.enable() diff --git a/pyccl/tests/test_cosmology.py b/pyccl/tests/test_cosmology.py index 3273deb67..fac3e7545 100644 --- a/pyccl/tests/test_cosmology.py +++ b/pyccl/tests/test_cosmology.py @@ -243,3 +243,80 @@ def test_cosmology_context(): with pytest.raises(AttributeError): cosmo.has_growth + + +def test_pyccl_default_params(): + """Check that the Python-layer for setting the gsl and spline parameters + works on par with the C-layer. + """ + HM_MMIN = ccl.gsl_params["HM_MMIN"] + + # we will test with this parameter + assert HM_MMIN == 1e7 + + # can be accessed as an attribute and as a dictionary item + assert ccl.gsl_params.HM_MMIN == ccl.gsl_params["HM_MMIN"] + + # can be assigned as an attribute + ccl.gsl_params.HM_MMIN = 1e5 + assert ccl.gsl_params["HM_MMIN"] == 1e5 # cross-check + + ccl.gsl_params["HM_MMIN"] = 1e6 + assert ccl.gsl_params.HM_MMIN == 1e6 + + # does not accept extra assignment + with pytest.raises(KeyError): + ccl.gsl_params.test = "hello_world" + with pytest.raises(KeyError): + ccl.gsl_params["test"] = "hallo_world" + + # complains when we try to set A_SPLINE_MAX != 1.0 + ccl.spline_params.A_SPLINE_MAX = 1.0 + with pytest.raises(RuntimeError): + ccl.spline_params.A_SPLINE_MAX = 0.9 + + # complains when we try to change the spline type + ccl.spline_params.A_SPLINE_TYPE = None + with pytest.raises(TypeError): + ccl.spline_params.A_SPLINE_TYPE = "something_else" + + # complains when we try to change the physical constants + with pytest.raises(AttributeError): + ccl.physical_constants.CLIGHT = 1 + + # verify that this has changed + assert ccl.gsl_params.HM_MMIN != HM_MMIN + + # but now we reload it, so it should be the default again + ccl.gsl_params.reload() + assert ccl.gsl_params.HM_MMIN == HM_MMIN + + +def test_cosmology_default_params(): + """Check that the default params within Cosmology work as intended.""" + cosmo1 = ccl.CosmologyVanillaLCDM() + v1 = cosmo1.cosmo.gsl_params.HM_MMIN + + ccl.gsl_params.HM_MMIN = v1*10 + cosmo2 = ccl.CosmologyVanillaLCDM() + v2 = cosmo2.cosmo.gsl_params.HM_MMIN + assert v2 == v1*10 + assert v2 != v1 + + ccl.gsl_params.reload() + cosmo3 = ccl.CosmologyVanillaLCDM() + v3 = cosmo3.cosmo.gsl_params.HM_MMIN + assert v3 == v1 + + # warns when we try to mutate instantiated `cvar` objects + with pytest.warns(ccl.CCLDeprecationWarning): + cosmo1.cosmo.spline_params.A_SPLINE_MIN = 0.5 + + +def test_ccl_physical_constants_smoke(): + assert ccl.physical_constants.CLIGHT == ccl.ccllib.cvar.constants.CLIGHT + + +def test_ccl_global_parameters_repr(): + ccl.spline_params.reload() + assert eval(repr(ccl.spline_params)) == ccl.spline_params._bak diff --git a/pyccl/tests/test_gsl_splines.py b/pyccl/tests/test_gsl_splines.py index f428cb257..c8c2ecfa4 100644 --- a/pyccl/tests/test_gsl_splines.py +++ b/pyccl/tests/test_gsl_splines.py @@ -1,20 +1,7 @@ -import pytest import pyccl as ccl - import numpy as np -def test_spline_params(): - cosmo = ccl.Cosmology( - Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.8, n_s=0.96, - transfer_function='bbks', matter_power_spectrum='linear') - - assert cosmo.cosmo.spline_params.A_SPLINE_MAX == 1.0 - - with pytest.raises(RuntimeError): - cosmo.cosmo.spline_params.A_SPLINE_MAX = 0.9 - - def test_spline1d(): cosmo = ccl.CosmologyVanillaLCDM() cosmo.compute_distances() diff --git a/pyccl/tests/test_pk2d.py b/pyccl/tests/test_pk2d.py index be71072f6..efb703e41 100644 --- a/pyccl/tests/test_pk2d.py +++ b/pyccl/tests/test_pk2d.py @@ -45,9 +45,6 @@ def test_pk2d_init(): assert_raises(ValueError, ccl.Pk2D, pkfunc=pk1d) ccl.Pk2D(pkfunc=lpk2d, cosmo=cosmo) - # Input function but no cosmo - assert_raises(ValueError, ccl.Pk2D, pkfunc=lpk2d) - # Input arrays have incorrect sizes lkarr = -4.+6*np.arange(100)/99. aarr = 0.05+0.95*np.arange(100)/99. @@ -364,3 +361,10 @@ def test_pk2d_mul_pow(): pk2d_j = (pk2d_a + 0.5*pk2d_i)**1.5 _, _, zarr_j = pk2d_j.get_spline_arrays() assert np.allclose((zarr_a + 0.5*zarr_i)**1.5, zarr_j) + + +def test_pk2d_pkfunc_init_without_cosmo(): + cosmo = ccl.CosmologyVanillaLCDM(transfer_function="bbks") + arr1 = ccl.Pk2D(pkfunc=lpk2d, cosmo=cosmo).get_spline_arrays()[-1] + arr2 = ccl.Pk2D(pkfunc=lpk2d).get_spline_arrays()[-1] + assert np.allclose(arr1, arr2, rtol=0) diff --git a/pyccl/tests/test_pt_power.py b/pyccl/tests/test_pt_power.py index 472b82e2f..c487e8eea 100644 --- a/pyccl/tests/test_pt_power.py +++ b/pyccl/tests/test_pt_power.py @@ -1,7 +1,6 @@ import numpy as np import pyccl as ccl import pytest -from .. import ccllib as lib NZ = 128 ZZ = np.linspace(0., 1., NZ) @@ -27,8 +26,8 @@ gz_1 = ccl.growth_factor(COSMO, 1./(1+ZZ_1)) a_1_v = a_1*np.ones_like(ZZ) Om_m = COSMO['Omega_m'] -rho_crit = lib.cvar.constants.RHO_CRITICAL -rho_m = lib.cvar.constants.RHO_CRITICAL * COSMO['Omega_m'] +rho_crit = ccl.physical_constants.RHO_CRITICAL +rho_m = ccl.physical_constants.RHO_CRITICAL * COSMO['Omega_m'] Om_m_fid = 0.3 c_1_t = -1*a_1*5e-14*rho_crit*COSMO['Omega_m']/gz diff --git a/pyccl/tests/test_tracers.py b/pyccl/tests/test_tracers.py index efabbc11f..d959e9202 100644 --- a/pyccl/tests/test_tracers.py +++ b/pyccl/tests/test_tracers.py @@ -190,21 +190,24 @@ def test_tracer_nz_support(): _ = ccl.CMBLensingTracer(calculator_cosmo, z_source=2.0) +def new_simple_cosmo(): + return ccl.CosmologyVanillaLCDM(transfer_function='bbks', + matter_power_spectrum="linear") + + def test_tracer_nz_norm_spline_vs_gsl_intergation(): # Create a new Cosmology object so that we're not messing with the other # tests - cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, - sigma8=0.8, n_s=0.96, - transfer_function='bbks', - matter_power_spectrum='linear') - cosmo.cosmo.gsl_params.NZ_NORM_SPLINE_INTEGRATION = True + ccl.gsl_params.NZ_NORM_SPLINE_INTEGRATION = True + cosmo = new_simple_cosmo() tr_wl, _ = get_tracer("wl", cosmo) tr_nc, _ = get_tracer("nc", cosmo) w_wl_spline, _ = tr_wl.get_kernel(chi=None) w_nc_spline, _ = tr_nc.get_kernel(chi=None) - cosmo.cosmo.gsl_params.NZ_NORM_SPLINE_INTEGRATION = False + ccl.gsl_params.NZ_NORM_SPLINE_INTEGRATION = False + cosmo = new_simple_cosmo() tr_wl, _ = get_tracer("wl", cosmo) tr_nc, _ = get_tracer("nc", cosmo) @@ -216,6 +219,8 @@ def test_tracer_nz_norm_spline_vs_gsl_intergation(): for w_spline, w_gsl in zip(w_nc_spline, w_nc_gsl): assert np.allclose(w_spline, w_gsl, atol=0, rtol=1e-8) + ccl.gsl_params.reload() # reset to the default parameters + @pytest.mark.parametrize('z_min, z_max, n_z_samples', [(0.0, 1.0, 2000), (0.0, 1.0, 1000), @@ -226,10 +231,6 @@ def test_tracer_lensing_kernel_spline_vs_gsl_intergation(z_min, z_max, n_z_samples): # Create a new Cosmology object so that we're not messing with the other # tests - cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, - sigma8=0.8, n_s=0.96, - transfer_function='bbks', - matter_power_spectrum='linear') z = np.linspace(z_min, z_max, n_z_samples) n = dndz(z) @@ -237,11 +238,13 @@ def test_tracer_lensing_kernel_spline_vs_gsl_intergation(z_min, z_max, if z_min > 0: assert n[0] > 0 - cosmo.cosmo.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = True + ccl.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = True + cosmo = new_simple_cosmo() tr_wl = ccl.WeakLensingTracer(cosmo, dndz=(z, n)) w_wl_spline, _ = tr_wl.get_kernel(chi=None) - cosmo.cosmo.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = False + ccl.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION = False + cosmo = new_simple_cosmo() tr_wl = ccl.WeakLensingTracer(cosmo, dndz=(z, n)) w_wl_gsl, chi = tr_wl.get_kernel(chi=None) @@ -251,6 +254,8 @@ def test_tracer_lensing_kernel_spline_vs_gsl_intergation(z_min, z_max, else: assert np.allclose(w_wl_spline[0], w_wl_gsl[0], atol=5e-9, rtol=1e-5) + ccl.gsl_params.reload() # reset to the default parameters + @pytest.mark.parametrize('z_min, z_max, n_z_samples', [(0.0, 1.0, 2000), (0.0, 1.0, 1000), diff --git a/pyccl/tracers.py b/pyccl/tracers.py index e6b977fa8..7dcbe281b 100644 --- a/pyccl/tracers.py +++ b/pyccl/tracers.py @@ -4,11 +4,12 @@ from . import ccllib as lib from .core import check -from .background import comoving_radial_distance, growth_rate, \ - growth_factor, scale_factor_of_chi, h_over_h0 +from .background import (comoving_radial_distance, growth_rate, + growth_factor, scale_factor_of_chi, h_over_h0) from .errors import CCLWarning -from .pyutils import _check_array_params, NoneArr, _vectorize_fn6, \ - _get_spline1d_arrays +from .pyutils import (_check_array_params, NoneArr, _vectorize_fn6, + _get_spline1d_arrays) +from .parameters import physical_constants def _Sig_MG(cosmo, a, k=None): @@ -110,7 +111,7 @@ def get_lensing_kernel(cosmo, dndz, mag_bias=None, n_chi=None): f"The number of samples in the n(z) ({len(z_n)}) is smaller than " f"the number of samples in the lensing kernel ({n_chi}). Consider " f"disabling spline integration for the lensing kernel by setting " - f"cosmo.cosmo.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION " + f"pyccl.gsl_params.LENSING_KERNEL_SPLINE_INTEGRATION " f"= False", category=CCLWarning) @@ -689,7 +690,7 @@ def __init__(self, cosmo, dndz, has_shear=True, ia_bias=None, # Transfer # See Joachimi et al. (2011), arXiv: 1008.3491, Eq. 6. # and note that we use C_1= 5e-14 from arXiv:0705.0166 - rho_m = lib.cvar.constants.RHO_CRITICAL * cosmo['Omega_m'] + rho_m = physical_constants.RHO_CRITICAL * cosmo['Omega_m'] a = - tmp_a * 5e-14 * rho_m / D else: # use the raw input normalization. Normally, this will be 1 @@ -820,7 +821,7 @@ def __init__(self, cosmo, z_max=6., n_chi=1024): self.chi_max = comoving_radial_distance(cosmo, 1./(1+z_max)) chi = np.linspace(0, self.chi_max, n_chi) a_arr = scale_factor_of_chi(cosmo, chi) - H0 = cosmo['h'] / lib.cvar.constants.CLIGHT_HMPC + H0 = cosmo['h'] / physical_constants.CLIGHT_HMPC OM = cosmo['Omega_c']+cosmo['Omega_b'] Ez = h_over_h0(cosmo, a_arr) fz = growth_rate(cosmo, a_arr) diff --git a/readthedocs/source/notation_and_other_cosmological_conventions.rst b/readthedocs/source/notation_and_other_cosmological_conventions.rst index a5e807802..1443176ba 100644 --- a/readthedocs/source/notation_and_other_cosmological_conventions.rst +++ b/readthedocs/source/notation_and_other_cosmological_conventions.rst @@ -122,13 +122,20 @@ quantities (e.g., the transfer function). The various options are as follows. Controlling Splines and Numerical Accuracy ------------------------------------------ -The internal splines and integration accuracy are controlled by the -attributes of :obj:`Cosmology.cosmo.spline_params` and -:obj:`Cosmology.cosmo.gsl_params`. These should be set after instantiation, -but before the object is used. For example, you can set the generic relative +The internal splines and integration accuracy are controlled by the global +instances :obj:`pyccl.spline_params` and :obj:`pyccl.gsl_params`. +Upon instantiation, the :obj:`pyccl.Cosmology` object assumes the accuracy +parameters from these instances. For example, you can set the generic relative accuracy for integration by executing -``c = Cosmology(...); c.cosmo.gsl_params.INTEGRATION_EPSREL = 1e-5``. The -default values for these parameters are located in ``src/ccl_core.c``. +``pyccl.gsl_params.INTEGRATION_EPSREL = 1e-5`` or +``pyccl.gsl_params["INTEGRATION_EPSREL"] = 1e-5``. To reset the accuracy +parameters to their default valus listed in ``src/ccl_core.c``, you may run +``pyccl.gsl_params.reload()`` or ``pyccl.spline_params.reload()``. + +.. note:: + Previously, the indicated way of setting the accuracy parameters was + after instantiation of a :obj:`Cosmology` object. For mutation consinstency, + this functionality is now deprecated and will raise an error in the future. The internal splines are controlled by the following parameters. diff --git a/src/ccl_core.c b/src/ccl_core.c index a9d0cf4a4..48d79c41f 100644 --- a/src/ccl_core.c +++ b/src/ccl_core.c @@ -67,7 +67,7 @@ const ccl_configuration default_config = { */ #define GSL_EPSREL_DNDZ 1E-6 -const ccl_gsl_params default_gsl_params = { +ccl_gsl_params ccl_user_gsl_params = { GSL_N_ITERATION, // N_ITERATION GSL_INTEGRATION_GAUSS_KRONROD_POINTS,// INTEGRATION_GAUSS_KRONROD_POINTS GSL_EPSREL, // INTEGRATION_EPSREL @@ -100,7 +100,7 @@ const ccl_gsl_params default_gsl_params = { #undef GSL_EPSREL_DNDZ -const ccl_spline_params default_spline_params = { +ccl_spline_params ccl_user_spline_params = { // scale factor spline params 250, // A_SPLINE_NA @@ -259,8 +259,8 @@ ccl_cosmology * ccl_cosmology_create(ccl_parameters params, ccl_configuration co ccl_cosmology * cosmo = malloc(sizeof(ccl_cosmology)); cosmo->params = params; cosmo->config = config; - cosmo->gsl_params = default_gsl_params; - cosmo->spline_params = default_spline_params; + cosmo->gsl_params = ccl_user_gsl_params; + cosmo->spline_params = ccl_user_spline_params; cosmo->spline_params.A_SPLINE_TYPE = gsl_interp_akima; cosmo->spline_params.K_SPLINE_TYPE = gsl_interp_akima; cosmo->spline_params.M_SPLINE_TYPE = gsl_interp_akima; @@ -577,6 +577,24 @@ void ccl_get_pk_spline_a_array(ccl_cosmology *cosmo,int ndout,double* doutput,in free(d); } +void ccl_get_pk_spline_a_array_from_params(ccl_spline_params *spline_params, + int ndout, double *doutput, int *status) { + double *d = NULL; + if (*status == 0) { + d = ccl_linlog_spacing(spline_params->A_SPLINE_MINLOG_PK, + spline_params->A_SPLINE_MIN_PK, + spline_params->A_SPLINE_MAX, + spline_params->A_SPLINE_NLOG_PK, + spline_params->A_SPLINE_NA_PK); + if (d == NULL) + *status = CCL_ERROR_MEMORY; + } + if(*status==0) + memcpy(doutput, d, ndout*sizeof(double)); + + free(d); +} + int ccl_get_pk_spline_nk(ccl_cosmology *cosmo) { double ndecades = log10(cosmo->spline_params.K_MAX) - log10(cosmo->spline_params.K_MIN); return (int)ceil(ndecades*cosmo->spline_params.N_K); @@ -597,3 +615,18 @@ void ccl_get_pk_spline_lk_array(ccl_cosmology *cosmo,int ndout,double* doutput,i } free(d); } + +void ccl_get_pk_spline_lk_array_from_params(ccl_spline_params *spline_params, + int ndout, double *doutput, int *status) { + double *d = NULL; + if (*status == 0) { + d = ccl_log_spacing(spline_params->K_MIN, spline_params->K_MAX, ndout); + if (d == NULL) + *status = CCL_ERROR_MEMORY; + } + if (*status == 0) { + for(int ii=0; ii < ndout; ii++) + doutput[ii] = log(d[ii]); + } + free(d); +}