tidy

Author: alexander-mead

.gitignore

@@ -5,6 +5,7 @@ dist/
 
 # Plots
 *.pdf
+*.png
 
 # VScode
 .vscode/

README.md

@@ -63,6 +63,12 @@ I compared it against the `CAMB-HMcode` version for 100 random sets of cosmologi
 
 These comparisons can be reproduced using the `comparisons/CAMB.py` script.
 
+The power-spectrum suppression caused by baryonic feedback has also been implemented, and the level of agreement between the *suppression* predicted by this code and the same implementation in `CAMB` is as follows:
+- LCDM: Mean error: 0.02%; Standard deviation of error: <0.01%; Worst-case error; 0.03%
+- nu-k-w(a)-CDM: Mean error: 0.06%; Standard deviation of error: 0.07%; Worst-case error; 0.49% (larger errors strongly correlated with neutrino mass)
+
+The comparisons can be reproduced using the `comparisons/CAMB_feedback.py` script.
+
 While the quoted accuracy of `HMcode-2020` relative to simulations is RMS ~2.5%, note that the accuracy is anti-correlated with neutrino masses (cf. Fig. 2 of [Mead et al. 2021](https://arxiv.org/abs/2009.01858)). The larger discrepancies between the codes for massive neutrinos (2% for ~1eV) may seem worrisome, but here are some reasons why I am not that worried:
 - Here, neutrinos are treated as completely cold matter when calculating the linear growth factors, whereas in `CAMB-HMcode` the transition from behaving like radiation to behaving like matter is accounted for in the linear growth.
 - Here the *cold* matter power spectrum is taken directly from `CAMB` whereas in `CAMB-HMcode` the *cold* spectrum is calculated approximately from the total matter power spectrum using approximations for the scale-dependent growth rate from [Eisenstein & Hu (1999)](https://arxiv.org/abs/astro-ph/9710252).

comparisons/CAMB.py

@@ -1,5 +1,7 @@
 # Standard imports
 import sys
+
+# Third-part imports
 import numpy as np
 
 # Project imports

comparisons/CAMB_Feedback.py

@@ -1,16 +1,17 @@
 # Standard imports
 import sys
+
+# Third-part imports
 import numpy as np
 import matplotlib.pyplot as plt
 
 # Project imports
-sys.path.append('../')
-
 import hmcode
 import hmcode.camb_stuff as camb_stuff
 import hmcode.utility as util
 
-plotting=True
+# Make plots or not
+plotting = False
 
 # Vary these parameters
 vary_Omega_k = True
@@ -30,7 +31,6 @@
 m_nu_min, m_nu_max = 0., 1.
 log10T_AGN_min, log10T_AGN_max = 7.6, 8.3
 
-
 # Number of cosmological models to test
 try:
     ncos = int(sys.argv[1])
@@ -69,7 +69,7 @@
         if w0+wa < 0.: break
     m_nu = rng.uniform(m_nu_min, m_nu_max) if vary_m_nu else 0.
 
-    #Gravity only spectra
+    ### Gravity only spectra ###
 
     # Get stuff from CAMB
     _, results_gravity, _, _, _ = camb_stuff.run(zs, Omega_c, Omega_b, Omega_k, h, ns, sigma_8, m_nu, w0, wa)
@@ -83,8 +83,9 @@
     # Get the new pyHMcode spectrum
     Pk_HMcode_gravity = hmcode.power(k, zs, results_gravity, verbose=verbose)
 
+    ### ###
 
-    # Spectra with feedback
+    ### Spectra with feedback ###
 
     # Get stuff from CAMB
     _, results_feedback, _, _, _ = camb_stuff.run(zs, Omega_c, Omega_b, Omega_k, h, ns, sigma_8, m_nu, w0, wa, log10_T_AGN=log10T_AGN)
@@ -98,11 +99,13 @@
     # Get the new pyHMcode spectrum
     Pk_HMcode_feedback = hmcode.power(k, zs, results_gravity, verbose=verbose, T_AGN=np.power(10, log10T_AGN))
 
-    # Calculate suppressions
+    ### ###
 
+    # Calculate suppressions
     Rk_CAMB=Pk_CAMB_feedback/Pk_CAMB_gravity
     Rk_HMcode=Pk_HMcode_feedback/Pk_HMcode_gravity
 
+    # Plotting
     if plotting:
         fig, axs=plt.subplots(nrows=2, sharex=True, figsize=(10,7))
         plt.subplots_adjust(hspace=0.001, bottom=0.2)
@@ -119,14 +122,13 @@
         for i,z in enumerate(zs):
             axs[0].plot(k, Rk_CAMB[i], label=f"CAMB, z={z}", color=f"C{i}", ls='--')
             axs[0].plot(k, Rk_HMcode[i], label=f"HMcode-Python, z={z}", color=f"C{i}")
-        
             axs[1].plot(k, 1-(Rk_CAMB[i]/Rk_HMcode[i]), label=f"z={z}", color=f"C{i}")
 
         axs[0].legend()
         axs[1].legend()
-        plt.savefig(f"Feedback_Cosmology_{icos}.png")
+        plt.savefig(f"plots/Feedback_Cosmology_{icos}.png")
+        plt.close()
 
-    
     # Calculate maximum deviation between pyHMcode and the version in CAMB
     max_error = np.max(np.abs(-1.+Rk_HMcode/Rk_CAMB))
     max_errors.append(max_error)

hmcode/camb_stuff.py

@@ -1,24 +1,24 @@
-# Standard imports
-import numpy as np
-
 # Third-party imports
 import camb
 
 def run(zs, Omega_c, Omega_b, Omega_k, h, ns, sigma_8, 
-    m_nu=0., w=-1., wa=0., As=2e-9, norm_sigma8=True, kmax_CAMB=200., log10_T_AGN=None, verbose=False):
+    m_nu=0., w=-1., wa=0., As=2e-9, norm_sigma8=True, kmax_CAMB=200., 
+    log10_T_AGN=None, verbose=False):
 
     # Sets cosmological parameters in camb to calculate the linear power spectrum
     if log10_T_AGN:
-        pars = camb.CAMBparams(NonLinearModel=camb.nonlinear.Halofit(halofit_version="mead2020_feedback", HMCode_logT_AGN=log10_T_AGN), WantCls=False)
+        non_linear_model = camb.nonlinear.Halofit(halofit_version="mead2020_feedback", 
+                                                  HMCode_logT_AGN=log10_T_AGN)
     else:
-        pars = camb.CAMBparams(WantCls=False,NonLinearModel=camb.nonlinear.Halofit(halofit_version="mead2020"))
+        non_linear_model = camb.nonlinear.Halofit(halofit_version="mead2020")
+    pars = camb.CAMBparams(WantCls=False, NonLinearModel=non_linear_model)
     wb, wc = Omega_b*h**2, Omega_c*h**2
 
     # This function sets standard and helium set using BBN consistency
     pars.set_cosmology(ombh2=wb, omch2=wc, H0=100.*h, mnu=m_nu, omk=Omega_k)
     pars.set_dark_energy(w=w, wa=wa, dark_energy_model='ppf')
     pars.InitPower.set_params(As=As, ns=ns, r=0.)
-    pars.set_matter_power(redshifts=zs, kmax=kmax_CAMB) # Setup the linear matter power spectrum
+    pars.set_matter_power(redshifts=zs, kmax=kmax_CAMB) # Linear matter power spectrum
     Omega_m = pars.omegam # Extract the matter density
 
     # Scale 'As' to be correct for the desired 'sigma_8' value if necessary

hmcode/constants.py

@@ -1,3 +1,4 @@
+# Standard imports
 from math import pi
 
 # Physical constants

hmcode/cosmology.py

@@ -1,4 +1,4 @@
-# Standard imports
+# Third-party imports
 import numpy as np
 import scipy.integrate as integrate
 

hmcode/hmcode.py

@@ -146,18 +146,19 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, T_AGN=None,
             print('RMS in matter displacement field: {:.4} Mpc/h'.format(sigmaV))
             print('Effective index at collapse scale: {:.4}'.format(neff))
 
-        if tweaks: # HMcode parameters (Table 2 of https://arxiv.org/pdf/2009.01858.pdf)
+        # HMcode parameters (Table 2 of https://arxiv.org/pdf/2009.01858.pdf)
+        ks = 0.05618*sigma8**-1.013  # One-halo damping wavenumber; equation (17)
+        if tweaks: 
             kd = 0.05699*sigma8**-1.089  # Two-halo damping wavenumber; equation (16)
             f = 0.2696*sigma8**0.9403    # Two-halo fractional damping; equation (16)
             nd = 2.853                   # Two-halo damping power; equation (16)
             ks = 0.05618*sigma8**-1.013  # One-halo damping wavenumber; equation (17)
             eta = 0.1281*sigma8**-0.3644 # Halo bloating parameter; equation (19)
-            B =   5.196                  # Minimum halo concentration; equation (20)
+            B = 5.196                  # Minimum halo concentration; equation (20)
             alpha = 1.875*(1.603)**neff  # Transition smoothing; equation (23)
-        else: # Use vanilla-ish halo model if no HMcode tweaks used (still uses 1-halo term suppression)
-            B, eta = 4., 0.
-            ks = 0.05618*sigma8**-1.013  # One-halo damping wavenumber; equation (17)
-
+        else: # Use vanilla-ish halo model if no HMcode tweaks used (still uses one-halo term suppression)
+            eta = 0.
+            B = 4.
 
         if T_AGN and not tweaks:
             B = feedback_params['B0']*np.power(10, z*feedback_params['Bz'])
@@ -201,16 +202,15 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, T_AGN=None,
         _, _Pk_1h, _ = hmod.power_spectrum(k, Pk_lin, M, sigmaM, {'m': profile}, simple_twohalo=True)
 
         # HMcode tweaks
+        # Still uses one-halo term dampening even if tweaks=False
+        Pk_1h = (k/ks)**4/(1.+(k/ks)**4)*_Pk_1h['m-m'] # One-halo term; equation (17)
         if tweaks:
             Pk_wig = _get_Pk_wiggle(k, Pk_lin, CAMB_results)   # Isolate spectral wiggle; footnote 7
             Pk_dwl = Pk_lin-(1.-np.exp(-(k*sigmaV)**2))*Pk_wig # Constuct linear spectrum with smoothed wiggle; equation (15)
             Pk_2h = Pk_dwl*(1.-f*(k/kd)**nd/(1.+(k/kd)**nd))   # Two-halo term; equation (16)
-            Pk_1h = (k/ks)**4/(1.+(k/ks)**4)*_Pk_1h['m-m']     # One-halo term; equation (17)
             Pk_hm = (Pk_2h**alpha+Pk_1h**alpha)**(1./alpha)    # Total prediction via smoothed sum; equation (23)
-        else: # Still uses 1halo term dampening even if tweaks=False
-            Pk_1h = (k/ks)**4/(1.+(k/ks)**4)*_Pk_1h['m-m']     # One-halo term; equation (17)
+        else: 
             Pk_hm = Pk_lin+Pk_1h
-
         Pk_HMcode[iz, :] = Pk_hm
 
     if T_AGN and tweaks:

hmcode/linear_growth.py

@@ -1,7 +1,5 @@
-# Standard imports
-import numpy as np
-
 # Third-party imports
+import numpy as np
 import camb
 
 # Parameters

hmcode/utility.py

@@ -1,4 +1,4 @@
-# Standard imports
+# Third-party imports
 import numpy as np
 
 def derivative_from_samples(x:float, xs:np.array, fs:np.array) -> float: