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 <dam.phys@gmail.com>

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
 
 

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
 
 

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
 
 

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
 
 

benchmarks/test_correlation_MG3_SD.py

@@ -17,15 +17,15 @@ 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,
                           mu_0=0.1, sigma_0=0.1,
                           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
 
 

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
 
 

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
 

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',

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)

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:

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);

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

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)

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