Eisenstein & Hu power spectrum without wiggles (#898)

* implemented, benchmarked and tested

benchmarks/data/codes/ehpk_bm.py

@@ -8,10 +8,14 @@
 FS=16
 LKMAX=7
 
-def do_all(z_arr,k_arr,cpar,prefix) :
+def do_all(z_arr,k_arr,cpar,prefix,wig):
     pcs=csm.PcsPar()
     pcs.background_set(cpar['om'],cpar['ol'],cpar['ob'],cpar['w0'],cpar['wa'],cpar['hh'],TCMB)
-    pcs.set_linear_pk('EH',-3,LKMAX,0.01,cpar['ns'],cpar['s8'])
+    if wig:
+        typ = 'EH'
+    else:
+        typ = 'EH_smooth'
+    pcs.set_linear_pk(typ,-3,LKMAX,0.01,cpar['ns'],cpar['s8'])
 
     gf0=pcs.growth_factor(1)
     a_arr=1./(z_arr+1)
@@ -40,4 +44,5 @@ def do_all(z_arr,k_arr,cpar,prefix) :
 
 cpar_model1={'om': 0.3,'ol': 0.7,'ob':0.05,'hh': 0.7,'s8': 0.8,'ns': 0.96,'w0': -1.0, 'wa': 0.0}
 
-do_all(z_arr,k_arr,cpar_model1,"model1")
+do_all(z_arr,k_arr,cpar_model1,"model1",True)
+do_all(z_arr,k_arr,cpar_model1,"model1_nowig",False)

benchmarks/data/codes/model1_pk_eh.txt

@@ -21,7 +21,7 @@
 7.943282347242817953e-02 8.234247003818394660e+03 3.082130362809773487e+03 1.462537929152181505e+03 8.370869741684097107e+02 5.391416368610861127e+02 3.754642686884068326e+02
 1.000000000000000056e-01 5.480571979312075200e+03 2.051412508656723730e+03 9.734398773167208674e+02 5.571505703857434355e+02 3.588433218594377081e+02 2.499025046555387917e+02
 1.258925411794167559e-01 4.054255788034309262e+03 1.517533401305944153e+03 7.201026228306375288e+02 4.121524055006657932e+02 2.654545219254098356e+02 1.848654993983242321e+02
-1.584893192461114264e-01 2.694750256932136381e+03 1.008661943615306427e+03 4.786320423140367097e+02 2.739461589710427916e+02 1.764401849714377306e+02 1.228749191089040806e+02
+1.584893192461114264e-01 2.694750256932136836e+03 1.008661943615306541e+03 4.786320423140367666e+02 2.739461589710428484e+02 1.764401849714377590e+02 1.228749191089040949e+02
 1.995262314968880846e-01 1.883550438987185544e+03 7.050247576000597292e+02 3.345495899275985039e+02 1.914802333393953120e+02 1.233264518680408202e+02 8.588592105436052293e+01
 2.511886431509582351e-01 1.203641589366091921e+03 4.505306055071631022e+02 2.137865765669408518e+02 1.223612426926826231e+02 7.880906370479171130e+01 5.488351380578266259e+01
 3.162277660168379412e-01 7.651145793799045123e+02 2.863871918150734359e+02 1.358969547514811893e+02 7.778093708486780145e+01 5.009627796225638718e+01 3.488760853014236574e+01

benchmarks/data/model1_nowig_pk_eh.txt

@@ -0,0 +1,42 @@
+# [0] k (Mpc/h)^-1, [2] P(k,z=0.0) (Mpc/h)^3, [3] P(k,z=1.0) (Mpc/h)^3, [4] P(k,z=2.0) (Mpc/h)^3, [5] P(k,z=3.0) (Mpc/h)^3, [6] P(k,z=4.0) (Mpc/h)^3, [7] P(k,z=5.0) (Mpc/h)^3
+1.000000000000000021e-03 3.872331892535939460e+03 1.449438144778402147e+03 6.877899417485541562e+02 3.936581675769015192e+02 2.535429595496815693e+02 1.765701540742688849e+02
+1.258925411794167490e-03 4.792496131638794395e+03 1.793861397905921649e+03 8.512262705486640471e+02 4.872013292395770918e+02 3.137911952196702714e+02 2.185276995484048257e+02
+1.584893192461114082e-03 5.910682433207593022e+03 2.212405552549933873e+03 1.049834580105615942e+03 6.008752556232207098e+02 3.870050291821305564e+02 2.695146327532546593e+02
+1.995262314968878920e-03 7.255061260093071724e+03 2.715614989859850084e+03 1.288618408737105256e+03 7.375437334830405689e+02 4.750289372519654876e+02 3.308154672852602971e+02
+2.511886431509579437e-03 8.848484478846297861e+03 3.312043306990487508e+03 1.571636624432805547e+03 8.995298653220901315e+02 5.793591573647094037e+02 4.034721999850293059e+02
+3.162277660168379394e-03 1.070195359141481458e+04 4.005808435207850835e+03 1.900843275191089560e+03 1.087952055042062966e+03 7.007160186245508839e+02 4.879864761009121139e+02
+3.981071705534973415e-03 1.280579185658176357e+04 4.793288309506639962e+03 2.274519612345117821e+03 1.301826572859133648e+03 8.384659313302437340e+02 5.839170566753464300e+02
+5.011872336272724625e-03 1.511891414134865954e+04 5.659104506599821434e+03 2.685368239386540608e+03 1.536976736965109922e+03 9.899188248723943389e+02 6.893905464351411183e+02
+6.309573444801930275e-03 1.755708663227460056e+04 6.571727781146434609e+03 3.118427843275159830e+03 1.784839405158499403e+03 1.149559446182140391e+03 8.005660614296740505e+02
+7.943282347242813790e-03 1.998070882915169932e+04 7.478904789257882840e+03 3.548903075220979190e+03 2.031222902079402502e+03 1.308247379365496045e+03 9.110781137527912961e+02
+1.000000000000000021e-02 2.217964098431700768e+04 8.301978903752007682e+03 3.939469653934281723e+03 2.254764588807451673e+03 1.452223614342725114e+03 1.011344774827168521e+03
+1.258925411794167490e-02 2.384693804046715013e+04 8.926058662131945312e+03 4.235609080241578340e+03 2.424260675957304557e+03 1.561390762664801741e+03 1.087369998455200630e+03
+1.584893192461114125e-02 2.454790236558594916e+04 9.188434010927439886e+03 4.360111892944841202e+03 2.495520149427265096e+03 1.607286769118157508e+03 1.119332490907193687e+03
+1.995262314968879874e-02 2.381440928856386381e+04 8.913882946021058160e+03 4.229831437984250442e+03 2.420953829017739736e+03 1.559260925590682746e+03 1.085886756084762737e+03
+2.511886431509580825e-02 2.161125735039322535e+04 8.089229340248717563e+03 3.838515356287223312e+03 2.196983162519477446e+03 1.415008398110248891e+03 9.854276818195515943e+02
+3.162277660168379134e-02 1.863273109620549076e+04 6.974348258809088293e+03 3.309480021580228822e+03 1.894188562303696926e+03 1.219987831035708950e+03 8.496131767070349952e+02
+3.981071705534973415e-02 1.558134063389388393e+04 5.832193646697312033e+03 2.767502803054424021e+03 1.583986644882675364e+03 1.020196441757478169e+03 7.104762176282841892e+02
+5.011872336272724798e-02 1.269592331892837683e+04 4.752163825912453831e+03 2.255005149946303391e+03 1.290657425067242912e+03 8.312722312624035794e+02 5.789072834534885033e+02
+6.309573444801933051e-02 1.003285314162034592e+04 3.755359935044460599e+03 1.782000011710724721e+03 1.019931837689651047e+03 6.569063161029486082e+02 4.574769090518951771e+02
+7.943282347242817953e-02 7.657647497243417092e+03 2.866305546579246311e+03 1.360124357163718969e+03 7.784703288282803442e+02 5.013884820621140648e+02 3.491725492437656726e+02
+1.000000000000000056e-01 5.638006897142429807e+03 2.110341387057673273e+03 1.001402912502947061e+03 5.731552783977755325e+02 3.691514555914945390e+02 2.570812043304188705e+02
+1.258925411794167559e-01 4.007973007427742232e+03 1.500209465879108620e+03 7.118820384745520187e+02 4.074473349877110877e+02 2.624241326155997172e+02 1.827551023741333722e+02
+1.584893192461114264e-01 2.757356596682299823e+03 1.032095889737991683e+03 4.897519597087348870e+02 2.803106694698074648e+02 1.805393678688091654e+02 1.257296359468435583e+02
+1.995262314968880846e-01 1.841449505771730173e+03 6.892661139128096011e+02 3.270717705651537131e+02 1.872002860924097263e+02 1.205698712072184975e+02 8.396620743713356205e+01
+2.511886431509582351e-01 1.197855340214856597e+03 4.483647761134693610e+02 2.127588433877828606e+02 1.217730172251266936e+02 7.843020600994280755e+01 5.461967306782447196e+01
+3.162277660168379412e-01 7.615854316921858072e+02 2.850662096732802979e+02 1.352701199785535664e+02 7.742216674947029276e+01 4.986520516832649719e+01 3.472668685083784368e+01
+3.981071705534973137e-01 4.747917768167375243e+02 1.777175436516141076e+02 8.433084187564783463e+01 4.826695284113574047e+01 3.108724035148035014e+01 2.164950203450288768e+01
+5.011872336272724660e-01 2.910681209932793081e+02 1.089486255322880197e+02 5.169849370832840663e+01 2.958975272009248414e+01 1.905783772549292365e+01 1.327209144158199905e+01
+6.309573444801935826e-01 1.758843626604303267e+02 6.583462145934841203e+01 3.123996054725472149e+01 1.788026383890171545e+01 1.151612079878587735e+01 8.019955385040265838e+00
+7.943282347242821562e-01 1.049608296022576610e+02 3.928749765130995542e+01 1.864277259321896096e+01 1.067023411092923801e+01 6.872365311827218548e+00 4.785992102163725725e+00
+1.000000000000000000e+00 6.194885113802808974e+01 2.318784400627754039e+01 1.100313658489404567e+01 6.297670731554953605e+00 4.056133486004525501e+00 2.824736746157934686e+00
+1.258925411794167504e+00 3.620191123066213379e+01 1.355060271376682479e+01 6.430055870086385461e+00 3.680257383235608692e+00 2.370339105609538866e+00 1.650730676288923560e+00
+1.584893192461114042e+00 2.096506188486006295e+01 7.847354319532653832e+00 3.723740395376997636e+00 2.131291447574120035e+00 1.372698411434089660e+00 9.559625336665190654e-01
+1.995262314968880846e+00 1.203998631256167329e+01 4.506642485287637889e+00 2.138499930889724787e+00 1.223975392860779943e+00 7.883244120961442336e-01 5.489979416172608984e-01
+2.511886431509582351e+00 6.860775499928461585e+00 2.568031353801404570e+00 1.218586761775587890e+00 6.974609579990366237e-01 4.492128705214448181e-01 3.128368695427091173e-01
+3.162277660168379523e+00 3.881173185388950397e+00 1.452747496214680689e+00 6.893603010217548288e-01 3.945569663472847655e-01 2.541218478315767082e-01 1.769732983513079838e-01
+3.981071705534973137e+00 2.180747310048380072e+00 8.162674127699860449e-01 3.873366506206752935e-01 2.216930298966967861e-01 1.427855727153536369e-01 9.943747055217483366e-02
+5.011872336272724660e+00 1.217588219282982331e+00 4.557509166668581857e-01 2.162637278143884900e-01 1.237790459516281832e-01 7.972222660585616405e-02 5.551945066800999945e-02
+6.309573444801936049e+00 6.758351275758447407e-01 2.529693323492713808e-01 1.200394532098342792e-01 6.870485931705021254e-02 4.425065907793334502e-02 3.081665412897263354e-02
+7.943282347242821118e+00 3.730875792167216987e-01 1.396490238096840530e-01 6.626650077985186960e-02 3.792778533866170837e-02 2.442810472629340993e-02 1.701200547206999092e-02
+1.000000000000000000e+01 2.049207774317469144e-01 7.670313385062266232e-02 3.639730619281920238e-02 2.083208257476780109e-02 1.341729526939935448e-02 9.343954559754126532e-03

benchmarks/test_eh.py

@@ -1,11 +1,16 @@
 import numpy as np
 import pyccl as ccl
+import pytest
 
 EH_TOLERANCE = 1.0e-5
 
 
-def test_eh():
-    model = 1
+@pytest.mark.parametrize('transfer,fname',
+                         [('eisenstein_hu',
+                           'model1_pk_eh.txt'),
+                          ('eisenstein_hu_nowiggles',
+                           'model1_nowig_pk_eh.txt')])
+def test_eh(transfer, fname):
     cosmo = ccl.Cosmology(
         Omega_c=0.25,
         Omega_b=0.05,
@@ -19,10 +24,10 @@ def test_eh():
         m_nu_type='normal',
         Omega_g=0,
         Omega_k=0,
-        transfer_function='eisenstein_hu',
+        transfer_function=transfer,
         matter_power_spectrum='linear')
 
-    data = np.loadtxt('./benchmarks/data/model%d_pk_eh.txt' % model)
+    data = np.loadtxt('./benchmarks/data/' + fname)
 
     k = data[:, 0] * cosmo['h']
     for i in range(1):

include/ccl_config.h

@@ -22,6 +22,7 @@ typedef enum transfer_function_t
   ccl_boltzmann_camb    = 4,
   ccl_boltzmann_isitgr  = 5,
   ccl_pklin_from_input  = 6,
+  ccl_eisenstein_hu_nowiggles = 7,
 } transfer_function_t;
 
 /**

include/ccl_power.h

@@ -7,7 +7,7 @@ CCL_BEGIN_DECLS
 
 ccl_f2d_t *ccl_compute_linpower_bbks(ccl_cosmology *cosmo, int *status);
 
-ccl_f2d_t *ccl_compute_linpower_eh(ccl_cosmology *cosmo, int *status);
+ccl_f2d_t *ccl_compute_linpower_eh(ccl_cosmology *cosmo, int wiggled, int *status);
 
 ccl_f2d_t *ccl_compute_power_emu(ccl_cosmology * cosmo, int * status);
 

pyccl/core.py

@@ -19,6 +19,7 @@
 transfer_function_types = {
     None: lib.transfer_none,
     'eisenstein_hu': lib.eisenstein_hu,
+    'eisenstein_hu_nowiggles': lib.eisenstein_hu_nowiggles,
     'bbks': lib.bbks,
     'boltzmann_class': lib.boltzmann_class,
     'boltzmann_camb': lib.boltzmann_camb,
@@ -763,7 +764,7 @@ def compute_linear_power(self):
 
         if (self['N_nu_mass'] > 0 and
                 self._config_init_kwargs['transfer_function'] in
-                ['bbks', 'eisenstein_hu']):
+                ['bbks', 'eisenstein_hu', 'eisenstein_hu_nowiggles']):
             warnings.warn(
                 "The '%s' linear power spectrum model does not properly "
                 "account for massive neutrinos!" %
@@ -815,7 +816,7 @@ def compute_linear_power(self):
                                    "power spectrum using CAMB and specified"
                                    " sigma8 but the non-linear power spectrum "
                                    "cannot be consistenty rescaled.")
-        elif trf in ['bbks', 'eisenstein_hu']:
+        elif trf in ['bbks', 'eisenstein_hu', 'eisenstein_hu_nowiggles']:
             rescale_s8 = False
             rescale_mg = False
             pk = Pk2D.pk_from_model(self,

pyccl/pk2d.py

@@ -122,7 +122,8 @@ def pk_from_model(Pk2D, cosmo, model):
             cosmo (:class:`~pyccl.core.Cosmology`): A Cosmology object.
             model (:obj:`str`): model to use. Three models allowed:
                 `'bbks'` (Bardeen et al. ApJ 304 (1986) 15).
-                `'eisenstein_hu'` (Eisenstein & Hu astro-ph/9710252).
+                `'eisenstein_hu'` (Eisenstein & Hu astro-ph/9709112).
+                `'eisenstein_hu_nowiggles'` (Eisenstein & Hu astro-ph/9709112).
                 `'emu'` (arXiv:1508.02654).
         """
         pk2d = Pk2D(empty=True)
@@ -132,7 +133,10 @@ def pk_from_model(Pk2D, cosmo, model):
             ret = lib.compute_linpower_bbks(cosmo.cosmo, status)
         elif model == 'eisenstein_hu':
             cosmo.compute_growth()
-            ret = lib.compute_linpower_eh(cosmo.cosmo, status)
+            ret = lib.compute_linpower_eh(cosmo.cosmo, 1, status)
+        elif model == 'eisenstein_hu_nowiggles':
+            cosmo.compute_growth()
+            ret = lib.compute_linpower_eh(cosmo.cosmo, 0, status)
         elif model == 'emu':
             ret = lib.compute_power_emu(cosmo.cosmo, status)
         else:

pyccl/tests/test_pk2d.py

@@ -70,7 +70,8 @@ def test_pk2d_smoke():
     assert_(not np.isnan(psp.eval(1E-2, 0.5, cosmo)))
 
 
-@pytest.mark.parametrize('model', ['bbks', 'eisenstein_hu'])
+@pytest.mark.parametrize('model', ['bbks', 'eisenstein_hu',
+                                   'eisenstein_hu_nowiggles'])
 def test_pk2d_from_model(model):
     cosmo_fixed = ccl.Cosmology(
         Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.8, n_s=0.96)

src/ccl_power.c

@@ -153,11 +153,11 @@ ccl_f2d_t *ccl_compute_linpower_bbks(ccl_cosmology *cosmo, int *status)
   return psp;
 }
 
-ccl_f2d_t *ccl_compute_linpower_eh(ccl_cosmology *cosmo, int *status)
+ccl_f2d_t *ccl_compute_linpower_eh(ccl_cosmology *cosmo, int wiggled, int *status)
 {
   ccl_f2d_t *psp = NULL;
   eh_struct *eh = NULL;
-  eh = ccl_eh_struct_new(&(cosmo->params),1);
+  eh = ccl_eh_struct_new(&(cosmo->params),wiggled);
   if (eh != NULL) {
     psp=ccl_compute_linpower_analytic(cosmo, eh,
                                       eh_power,