Merge pull request #899 from patricialarsen/emuissue_897

fix for issue 897, cosmic emulator omega_M missing massive neutrinos

CHANGELOG.md

@@ -6,6 +6,7 @@
 
 ## C library
 - Deprecate C-level yaml reader and writer (#852)
+- Fixed bug in emulator implementation (#899)
 
 # v2.2.0 Changes
 

benchmarks/data/emu_input_cosmologies.txt

@@ -0,0 +1,3 @@
+3.2759e-01 5.9450e-02 6.167e-01 8.778e-01 9.611e-01 -7.000e-01 6.722e-01 0.000e00
+1.9827e-01 4.1627e-02 7.352e-01 7.113e-01 9.851e-01 -8.971e-01 3.247e-01 6.9606e-03
+2.7631e-01 5.2441e-02 6.477e-01 8.631e-01 8.985e-01 -8.632e-01 8.285e-01 2.5792e-03

benchmarks/data/emu_nu_cosmologies.txt

@@ -1,4 +1,4 @@
-3.5595e-01 6.5074e-02 5.9020e-01 7.3252e-01 9.5623e-01 -8.0194e-01 3.6280e-01 2.0317e-02
-1.8371e-01 4.1553e-02 7.3266e-01 8.3109e-01 9.9143e-01 -7.7313e-01 4.8961e-01 3.6764e-03
-1.9260e-01 3.5647e-02 8.0155e-01 8.0444e-01 9.3760e-01 -9.5613e-01 -3.5868e-02 1.3898e-03
-2.6751e-01 6.0584e-02 6.1588e-01 7.7310e-01 8.5881e-01 -9.0428e-01 8.0951e-01 2.4240e-02
+3.3563e-01 6.5074e-02 5.9020e-01 7.3252e-01 9.5623e-01 -8.0194e-01 3.6280e-01 2.0317e-02
+1.8003e-01 4.1553e-02 7.3266e-01 8.3109e-01 9.9143e-01 -7.7313e-01 4.8961e-01 3.6764e-03
+1.9121e-01 3.5647e-02 8.0155e-01 8.0444e-01 9.3760e-01 -9.5613e-01 -3.5868e-02 1.3898e-03
+2.4327e-01 6.0584e-02 6.1588e-01 7.7310e-01 8.5881e-01 -9.0428e-01 8.0951e-01 2.4240e-02

benchmarks/test_emu.py

@@ -30,7 +30,7 @@ def test_emu_nu(model):
         Neff=3.04,
         Omega_g=0,
         Omega_k=0,
-        transfer_function='boltzmann_class',
+        transfer_function='boltzmann_camb',
         matter_power_spectrum='emu',
     )
 
@@ -73,3 +73,47 @@ def test_emu(model):
     pk = ccl.nonlin_matter_power(cosmo, k, a)
     err = np.abs(pk/data[:, 1]-1)
     assert np.allclose(err, 0, rtol=0, atol=EMU_TOLERANCE)
+
+
+@pytest.mark.parametrize('model', list(range(3)))
+def test_emu_lin(model):
+    cosmos = np.loadtxt("./benchmarks/data/emu_input_cosmologies.txt")
+
+    mnu = ccl.nu_masses(
+        cosmos[model, 7] * cosmos[model, 2]**2, 'equal', T_CMB=2.725)
+
+    cosmo = ccl.Cosmology(
+        Omega_c=cosmos[model, 0],
+        Omega_b=cosmos[model, 1],
+        h=cosmos[model, 2],
+        sigma8=cosmos[model, 3],
+        n_s=cosmos[model, 4],
+        w0=cosmos[model, 5],
+        wa=cosmos[model, 6],
+        m_nu=mnu,
+        m_nu_type='list',
+        Neff=3.04,
+        Omega_g=0,
+        Omega_k=0,
+        transfer_function='boltzmann_camb',
+        matter_power_spectrum='emu',
+    )
+
+    a = 1
+    k = np.logspace(-3, -2, 50)
+
+    # Catch warning about neutrino linear growth
+    if (np.sum(mnu) > 0):
+        pk = ccl.pyutils.assert_warns(ccl.CCLWarning,
+                                      ccl.nonlin_matter_power,
+                                      cosmo, k, a)
+    else:
+        pk = ccl.nonlin_matter_power(cosmo, k, a)
+
+    # Catch warning about linear matter power
+    pk_lin = ccl.pyutils.assert_warns(ccl.CCLWarning,
+                                      ccl.linear_matter_power,
+                                      cosmo, k, a)
+
+    err = np.abs(pk/pk_lin-1)
+    assert np.allclose(err, 0, rtol=0, atol=EMU_TOLERANCE)

src/ccl_power.c

@@ -320,7 +320,12 @@ ccl_f2d_t *ccl_compute_power_emu(ccl_cosmology * cosmo, int * status)
     //For each redshift:
     for (int j = 0; j < na; j++){
       //Turn cosmology into emu_par:
-      emu_par[0] = (cosmo->params.Omega_c+cosmo->params.Omega_b)*cosmo->params.h*cosmo->params.h;
+      if ((cosmo->params.N_nu_mass>0) &&
+          (cosmo->config.emulator_neutrinos_method == ccl_emu_equalize)){
+          emu_par[0] = (cosmo->params.Omega_c+cosmo->params.Omega_b)*cosmo->params.h*cosmo->params.h + Omeganuh2_eq;
+      }else{
+          emu_par[0] = (cosmo->params.Omega_c+cosmo->params.Omega_b+cosmo->params.Omega_nu_mass)*cosmo->params.h*cosmo->params.h;
+      }
       emu_par[1] = cosmo->params.Omega_b*cosmo->params.h*cosmo->params.h;
       emu_par[2] = cosmo->params.sigma8;
       emu_par[3] = cosmo->params.h;