Calculator mode (#843)

* functions depend on psp * musigma just rescales * analytic standalone * linear analytic in python * non-linear almost decoupled * rescaling factorized * extricated linear * bcm factorized * nonlinear extricated, correlations generalized * flaked * test implemented, measure time on GHA * typo * actual timing * sigma_M generalized and sped up * mac is slow * started calculator creator * flaked * sigmaM splining wasn't good enough * bg and grth unit-tested * linear done * nonlinear implemented * precision * outdated test * bug kNL * tests removed * bcm simplified * core tests * pk2dtests * tests power * calculator types * growth normalization * outdated stuff * outdated * outdated * C cleanup * proper parsing * cleaned up C-level error messages * docstrings + power * classmethod to subclass * missing tests * memleak * easy PR comments * simplified api * get_power methods * simpler rescaling * test unnorm * missing rescaling * default name * nonlinear model * x to : * flake * catch warning
LSSTDESC · Jan 11, 2021 · 641ba87 · 641ba87
1 parent 593ed60
commit 641ba87
Show file tree

Hide file tree

Showing 56 changed files with 2,305 additions and 1,719 deletions.
diff --git a/benchmarks/ccl_test_f2d.c b/benchmarks/ccl_test_f2d.c
@@ -277,11 +277,13 @@ CTEST2(f2d,pk) {
   //Now verify that things scale with the CCL growth factor as intended
   //First initialize the cosmology object
   double gz;
+  double mnu=0.;
   ccl_configuration config = default_config;
   config.transfer_function_method = ccl_bbks;
   config.matter_power_spectrum_method = ccl_linear;
-  ccl_parameters params = ccl_parameters_create_flat_lcdm(data->Omega_c,data->Omega_b,data->h,
-							  data->A_s,data->n_s, &status);
+  ccl_parameters params = ccl_parameters_create(data->Omega_c,data->Omega_b, 0.0, 3.046, &mnu,
+                                                0, -1.0, 0.0, data->h, data->A_s,data->n_s,
+                                                -1, -1, -1, 0., 0., -1, NULL, NULL, &status);
   params.T_CMB=2.7;
   params.Omega_k=0;
   params.Omega_g=0;

diff --git a/benchmarks/data/dNdzs.npz b/benchmarks/data/dNdzs.npz
diff --git a/benchmarks/test_emu.py b/benchmarks/test_emu.py
@@ -36,7 +36,11 @@ def test_emu_nu(model):
 
     a = 1
     k = data[:, 0]
-    pk = ccl.nonlin_matter_power(cosmo, k, a)
+
+    # Catch warning about neutrino linear growth
+    pk = ccl.pyutils.assert_warns(ccl.CCLWarning,
+                                  ccl.nonlin_matter_power,
+                                  cosmo, k, a)
     err = np.abs(pk/data[:, 1]-1)
     assert np.allclose(err, 0, rtol=0, atol=EMU_TOLERANCE)
 

diff --git a/benchmarks/test_hod.py b/benchmarks/test_hod.py
@@ -33,17 +33,20 @@ def _nz_2mrs(z):
     dndz = _nz_2mrs(z_arr)
 
     # CCL prediction
-    cosmo = ccl.Cosmology(Omega_b=0.05,
-                          Omega_c=0.25,
-                          h=0.67,
-                          n_s=0.9645,
-                          A_s=2.0E-9,
-                          m_nu=0.00001,
-                          m_nu_type='equal')
     # Make sure we use the same P(k)
-    cosmo._set_linear_power_from_arrays(a_array=1/(1.+zs[::-1]),
-                                        k_array=ks,
-                                        pk_array=pks[::-1, :])
+    cosmo = ccl.CosmologyCalculator(
+        Omega_b=0.05,
+        Omega_c=0.25,
+        h=0.67,
+        n_s=0.9645,
+        A_s=2.0E-9,
+        m_nu=0.00001,
+        m_nu_type='equal',
+        pk_linear={'a': 1./(1.+zs[::-1]),
+                   'k': ks,
+                   'delta_matter:delta_matter': pks[::-1, :]})
+    cosmo.compute_growth()
+
     # Halo model setup
     mass_def = ccl.halos.MassDef(200, 'critical')
     cm = ccl.halos.ConcentrationDuffy08(mass_def)

diff --git a/benchmarks/test_timing.py b/benchmarks/test_timing.py
@@ -0,0 +1,39 @@
+import numpy as np
+import pyccl as ccl
+import time
+
+
+def test_timing():
+    ls = np.unique(np.geomspace(2, 2000, 128).astype(int)).astype(float)
+    nl = len(ls)
+    b_g = np.array([1.376695, 1.451179, 1.528404,
+                    1.607983, 1.689579, 1.772899,
+                    1.857700, 1.943754, 2.030887,
+                    2.118943])
+    dNdz_file = np.load('benchmarks/data/dNdzs.npz')
+    z_s = dNdz_file['z_sh']
+    dNdz_s = dNdz_file['dNdz_sh'].T
+    z_g = dNdz_file['z_cl']
+    dNdz_g = dNdz_file['dNdz_cl'].T
+    nt = len(dNdz_s) + len(dNdz_g)
+    nx = (nt*(nt+1))//2
+
+    cosmo = ccl.CosmologyVanillaLCDM(transfer_function='boltzmann_class')
+    cosmo.compute_nonlin_power()
+
+    start = time.time()
+    t_g = [ccl.NumberCountsTracer(cosmo, True,
+                                  (z_g, ng),
+                                  bias=(z_g, np.full(len(z_g), b)))
+           for ng, b in zip(dNdz_g, b_g)]
+    t_s = [ccl.WeakLensingTracer(cosmo, (z_s, ns)) for ns in dNdz_s]
+    t_all = t_g + t_s
+
+    cls = np.zeros([nx, nl])
+    ind1, ind2 = np.triu_indices(nt)
+    for ix, (i1, i2) in enumerate(zip(ind1, ind2)):
+        cls[ix, :] = ccl.angular_cl(cosmo, t_all[i1], t_all[i2], ls)
+    end = time.time()
+    t_seconds = end - start
+    print(end-start)
+    assert t_seconds < 3.
diff --git a/include/ccl_bcm.h b/include/ccl_bcm.h
@@ -17,6 +17,8 @@ CCL_BEGIN_DECLS
  */
 double ccl_bcm_model_fka(ccl_cosmology * cosmo, double k, double a, int *status);
 
+void ccl_bcm_correct(ccl_cosmology *cosmo, ccl_f2d_t *psp, int *status);
+
 CCL_END_DECLS
 
 #endif
diff --git a/include/ccl_cls.h b/include/ccl_cls.h
@@ -10,7 +10,7 @@ CCL_BEGIN_DECLS
  * @param cosmo Cosmological parameters
  * @param trc1 a ccl_cl_tracer_collection_t containing a bunch of individual contributions.
  * @param trc2 a ccl_cl_tracer_collection_t containing a bunch of individual contributions.
- * @param psp the p2d_t object representing the 3D power spectrum to integrate over. Pass null to use the non-linear matter power spectrum.
+ * @param psp the p2d_t object representing the 3D power spectrum to integrate over.
  * @param nl_out number of multipoles on which the power spectrum will be calculated.
  * @param l_out multipole values on which the power spectrum will be calculated.
  * @param cl_out will hold the calculated power spectrum values.
@@ -29,7 +29,7 @@ void ccl_angular_cls_limber(ccl_cosmology *cosmo,
  * @param cosmo Cosmological parameters
  * @param trc1 a ccl_cl_tracer_collection_t containing a bunch of individual contributions.
  * @param trc2 a ccl_cl_tracer_collection_t containing a bunch of individual contributions.
- * @param psp the p2d_t object representing the 3D power spectrum to integrate over. Pass null to use the non-linear matter power spectrum.
+ * @param psp the p2d_t object representing the 3D power spectrum to integrate over.
  * @param nl_out number of multipoles on which the power spectrum will be calculated.
  * @param l_out multipole values on which the power spectrum will be calculated.
  * @param cl_out will hold the calculated power spectrum values.

diff --git a/include/ccl_config.h b/include/ccl_config.h
@@ -33,8 +33,9 @@ typedef enum transfer_function_t
  */
 typedef enum matter_power_spectrum_t
 {
-    ccl_linear           = 0,
-    ccl_halofit          = 1,
+    ccl_pknl_none        = 0,
+    ccl_linear           = 1,
+    ccl_halofit          = 2,
     ccl_halo_model       = 3,
     ccl_emu              = 4,
     ccl_pknl_from_input  = 5
@@ -104,7 +105,6 @@ typedef struct ccl_configuration {
   mass_function_t          mass_function_method;
   halo_concentration_t     halo_concentration_method;
   emulator_neutrinos_t     emulator_neutrinos_method;
-  // TODO: Halo definition
 } ccl_configuration;
 
 /**

diff --git a/include/ccl_core.h b/include/ccl_core.h
@@ -107,6 +107,8 @@ typedef struct ccl_spline_params {
   double A_SPLINE_MIN;
   double A_SPLINE_MINLOG_PK;
   double A_SPLINE_MIN_PK;
+  double A_SPLINE_MINLOG_SM;
+  double A_SPLINE_MIN_SM;
   double A_SPLINE_MAX;
   double A_SPLINE_MINLOG;
   int A_SPLINE_NLOG;
@@ -118,6 +120,8 @@ typedef struct ccl_spline_params {
   double LOGM_SPLINE_MAX;
 
   //PS a and k spline
+  int A_SPLINE_NA_SM;
+  int A_SPLINE_NLOG_SM;
   int A_SPLINE_NA_PK;
   int A_SPLINE_NLOG_PK;
 
@@ -266,12 +270,7 @@ typedef struct ccl_data {
   gsl_spline * achi;
 
   // Function of Halo mass M
-  gsl_spline * logsigma;
-  gsl_spline * dlnsigma_dlogm;
-
-  // power spectrum splines
-  ccl_f2d_t * p_lin;
-  ccl_f2d_t * p_nl;
+  gsl_spline2d * logsigma;
 
   // real-space splines for RSD
   ccl_f1d_t* rsd_splines[3];
@@ -290,8 +289,6 @@ typedef struct ccl_cosmology {
 
   bool computed_distances;
   bool computed_growth;
-  bool computed_linear_power;
-  bool computed_nonlin_power;
   bool computed_sigma;
 
   int status;
@@ -338,13 +335,6 @@ ccl_parameters ccl_parameters_create(double Omega_c, double Omega_b, double Omeg
 				     double mu_0, double sigma_0, int nz_mgrowth, double *zarr_mgrowth,
 				     double *dfarr_mgrowth, int *status);
 
-/* ------- ROUTINE: ccl_parameters_create_flat_lcdm --------
-INPUT: some cosmological parameters needed to create a flat LCDM model
-TASK: call ccl_parameters_create to produce an LCDM model
-*/
-ccl_parameters ccl_parameters_create_flat_lcdm(double Omega_c, double Omega_b, double h,
-double norm_pk, double n_s, int *status);
-
 
 /**
  * Free a parameters struct

diff --git a/include/ccl_correlation.h b/include/ccl_correlation.h
@@ -44,35 +44,40 @@ void ccl_correlation(ccl_cosmology *cosmo,
 /**
  * Computes the 3dcorrelation function (wrapper)
  * @param cosmo :Cosmological parameters
+ * @param psp: power spectrum
  * @param a : scale factor
  * @param n_r : number of output values of distance r
  * @param r : values of the distance in Mpc
  * @param xi : the values of the correlation function at the distances above will be returned in this array, which should be pre-allocated
  * @param do_taper_pk : key for tapering (using cosine tapering by default)
  * @param taper_pk_limits: limits of tapering
  */
-void ccl_correlation_3d(ccl_cosmology *cosmo,double a,
-		     int n_r,double *r,double *xi,
-		     int do_taper_pk,double *taper_pk_limits,
-		     int *status);
+void ccl_correlation_3d(ccl_cosmology *cosmo,ccl_f2d_t *psp,
+                        double a,int n_r,double *r,double *xi,
+                        int do_taper_pk,double *taper_pk_limits,
+                        int *status);
 
-void ccl_correlation_multipole(ccl_cosmology *cosmo,double a,double beta,
-			   int l,int n_s,double *s,double *xi,
-			   int *status);
+void ccl_correlation_multipole(ccl_cosmology *cosmo,ccl_f2d_t *psp,
+                               double a,double beta,
+                               int l,int n_s,double *s,double *xi,
+                               int *status);
 
-void ccl_correlation_multipole_spline(ccl_cosmology *cosmo,double a,int *status);
+void ccl_correlation_multipole_spline(ccl_cosmology *cosmo,ccl_f2d_t *psp,
+                                      double a,int *status);
 
-void ccl_correlation_3dRsd(ccl_cosmology *cosmo,double a,
+void ccl_correlation_3dRsd(ccl_cosmology *cosmo,ccl_f2d_t *psp,double a,
 			   int n_s,double *s,double mu,double beta,double *xi,
 			   int use_spline, int *status);
 
-void ccl_correlation_3dRsd_avgmu(ccl_cosmology *cosmo, double a, int n_s, double *s,
+void ccl_correlation_3dRsd_avgmu(ccl_cosmology *cosmo, ccl_f2d_t *psp,
+                                 double a, int n_s, double *s,
                                  double beta, double *xi,
                                  int *status);
 
-void ccl_correlation_pi_sigma(ccl_cosmology *cosmo,double a,double beta,
-			   double pi,int n_sig,double *sig,double *xi,
-			   int use_spline,int *status);
+void ccl_correlation_pi_sigma(ccl_cosmology *cosmo,ccl_f2d_t *psp,
+                              double a,double beta,double pi,
+                              int n_sig,double *sig,double *xi,
+                              int use_spline,int *status);
 
 CCL_END_DECLS
 

diff --git a/include/ccl_halofit.h b/include/ccl_halofit.h
@@ -22,7 +22,8 @@ typedef struct halofit_struct {
  * @param cosmo Cosmological data
  * @return int, status of computations
  */
-halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo, int *status);
+halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo,
+                                       ccl_f2d_t *plin, int *status);
 
 /*
  * Free a halofit struct
@@ -39,7 +40,8 @@ void ccl_halofit_struct_free(halofit_struct *hf);
  * @param hf: halofit splines for evaluating the power spectrum
  * @return halofit_matter_power: halofit power spectrum, P(k), units of Mpc^{3}
  */
-double ccl_halofit_power(ccl_cosmology *cosmo, double k, double a, halofit_struct *hf, int *status);
+double ccl_halofit_power(ccl_cosmology *cosmo, ccl_f2d_t *plin,
+                         double k, double a, halofit_struct *hf, int *status);
 
 CCL_END_DECLS
 

diff --git a/include/ccl_massfunc.h b/include/ccl_massfunc.h
@@ -8,10 +8,11 @@ CCL_BEGIN_DECLS
  * Computes sigma(R), the power spectrum normalization, over log-spaced values of mass and radii
  * The result is attached to the cosmology object
  * @param cosmo Cosmological parameters
+ * @param psp linear matter power spectrum
  * @param status Status flag. 0 if there are no errors, nonzero otherwise.
  * For specific cases see documentation for ccl_error.
  */
-void ccl_cosmology_compute_sigma(ccl_cosmology *cosmo, int *status);
+void ccl_cosmology_compute_sigma(ccl_cosmology *cosmo, ccl_f2d_t *psp, int *status);
 
 /**
  * Calculate the standard deviation of density at smoothing mass M via interpolation.
@@ -30,11 +31,12 @@ double ccl_sigmaM(ccl_cosmology *cosmo, double log_halomass, double a, int *stat
  * via interpolation.
  * @param cosmo Cosmological parameters
  * @param log_halomass log10(Mass) to compute at, in units of Msun
+ * @param a scale factor, normalized to a=1 today
  * @param status Status flag. 0 if there are no errors, nonzero otherwise.
  * For specific cases see documentation for ccl_error.
  * @return sigmaM, the standard deviation of density at mass scale M
  */
-double ccl_dlnsigM_dlogM(ccl_cosmology *cosmo, double log_halomass, int *status);
+double ccl_dlnsigM_dlogM(ccl_cosmology *cosmo, double log_halomass, double a, int *status);
 
 /**
  * Fitting function for the spherical-model critical linear density for collapse

diff --git a/include/ccl_musigma.h b/include/ccl_musigma.h
@@ -10,7 +10,8 @@ CCL_BEGIN_DECLS
  ^ @param psp The linear power spectrum to spline.
  * @param status, integer indicating the status
  */
-void ccl_cosmology_spline_linpower_musigma(ccl_cosmology* cosmo, ccl_f2d_t *psp, int isitgr_flag, int* status);
+void ccl_rescale_musigma_s8(ccl_cosmology* cosmo, ccl_f2d_t *psp,
+                            int mg_rescale, int* status);
 
 CCL_END_DECLS
 

diff --git a/include/ccl_neutrinos.h b/include/ccl_neutrinos.h
@@ -31,15 +31,6 @@ typedef enum ccl_neutrino_mass_splits{
   ccl_nu_single=4
 } ccl_neutrino_mass_splits;
 
-/**
- * Spline for the phasespace integral required for getting the fractional energy density of massive neutrinos.
- * Returns a gsl spline for the phase space integral needed for massive neutrinos.
- * @param status Status flag. 0 if there are no errors, nonzero otherwise.
- * For specific cases see documentation for ccl_error.c
- * @return spl, the gsl spline for the phasespace integral required for massive neutrino calculations.
- */
-gsl_spline* calculate_nu_phasespace_spline(int *status);
-
 /**
  * Returns density of one neutrino species at a scale factor a.
  * Users are encouraged to access this quantity via the function ccl_omega_x.
@@ -51,8 +42,7 @@ gsl_spline* calculate_nu_phasespace_spline(int *status);
  * For specific cases see documentation for ccl_error.c
  * @return OmNuh2 Fractional energy density of neutrions with mass mnu, multiplied by h squared.
  */
-
-double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, int * status);
+double ccl_Omeganuh2(double a, int N_nu_mass, double* mnu, double T_CMB, int * status);
 
 /**
  * Returns mass of one neutrino species at a scale factor a.
@@ -64,7 +54,7 @@ double ccl_Omeganuh2 (double a, int N_nu_mass, double* mnu, double T_CMB, int *
  * For specific cases see documentation for ccl_error.c
  * @return Mnu Neutrino mass [eV].
  */
-double* ccl_nu_masses (double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB, int * status);
+double* ccl_nu_masses(double OmNuh2, ccl_neutrino_mass_splits mass_split, double T_CMB, int * status);
 
 CCL_END_DECLS
 #endif