changes

Author: alexander-mead

pyhmcode/cosmology.py

@@ -102,6 +102,8 @@ def sigmaR_vec(R:float, Pk:callable):
     return sigma_func(R, Pk)
 
 
+def _sigmaV_integrand(k:float, Pk:callable):
+    return Pk(k)
 # def _sigmaV_integrand(t, Pk, alpha=3.):
 #     k = (-1.+1./t)**alpha
 #     return Pk(k)*k*alpha/(t*(1.-t))
@@ -119,28 +121,41 @@ def sigmaV(Pk:callable, eps=eps_sigmaV) -> float:
         eps: Integration accuracy
     '''
     kmin, kmax = 0., np.inf
-    sigmaV_squared, _ = integrate.quad(Pk, kmin, kmax, epsrel=eps, epsabs=eps)
+    #sigmaV_squared, _ = integrate.quad(Pk, kmin, kmax, epsrel=eps, epsabs=eps)
+    sigmaV_squared, _ = integrate.quad(lambda k: _sigmaV_integrand(k, Pk),  kmin, kmax, epsrel=eps, epsabs=eps)
     #sigmaV_squared, _ = integrate.quad(_sigmaV_integrand, 0., 1., args=(Pk,), epsabs=eps, epsrel=eps)
     sigmaV = np.sqrt(sigmaV_squared/(2.*np.pi**2))
     sigmaV /= np.sqrt(3.) # Convert from 3D displacement to 1D displacement
     return sigmaV
 
 
-def _dsigmaR_integrand(k:float, R:float, Pk) -> float:
+def _dsigmaR_integrand(k:float, R:float, Pk:callable) -> float:
     return Pk(k)*(k**3)*_Tophat_k(k*R)*_dTophat_k(k*R)
 
 
-def dlnsigma2_dlnR(R:float, Pk) -> float:
+def dlnsigma2_dlnR(R:float, Pk:callable) -> float:
     '''
     Calculates d(ln sigma^2)/d(ln R) by integration
+    3+neff = -d(ln sigma^2) / dR
     '''
-    def dsigmaR_vec(R, Pk):
-        kmin, kmax = 0., np.inf # Evaluate the integral and convert to a nicer form
-        dsigma, _ = integrate.quad(lambda k: _dsigmaR_integrand(k, R, Pk), kmin, kmax)
-        dsigma = R*dsigma/(np.pi*_sigmaR_quad(R, Pk))**2
-        return dsigma
-    dsigma_func = np.vectorize(dsigmaR_vec, excluded=['Pk'])
-    return dsigma_func(R, Pk)
+    # def dsigmaR_vec(R, Pk):
+    #     kmin, kmax = 0., np.inf # Evaluate the integral and convert to a nicer form
+    #     dsigma, _ = integrate.quad(lambda k: _dsigmaR_integrand(k, R, Pk), kmin, kmax)
+    #     dsigma = R*dsigma/(np.pi*_sigmaR_quad(R, Pk))**2
+    #     return dsigma
+    # dsigma_func = np.vectorize(dsigmaR_vec, excluded=['Pk'])
+    # return dsigma_func(R, Pk)
+    kmin, kmax = 0., np.inf # Evaluate the integral and convert to a nicer form
+    dsigma, _ = integrate.quad(lambda k: _dsigmaR_integrand(k, R, Pk), kmin, kmax)
+    dsigma = R*dsigma/(np.pi*_sigmaR_quad(R, Pk))**2
+    return dsigma
+
+
+def neff(R:float, Pk:callable) -> float:
+    '''
+    Effective index of the power spectrum at scale 'R'
+    '''
+    return -3.-dlnsigma2_dlnR(R, Pk)
 
 ### ###
 

pyhmcode/pyhmcode.py

@@ -1,5 +1,6 @@
 # Standard imports
 import numpy as np
+from time import time
 
 # Third-party imports
 import camb
@@ -26,7 +27,8 @@ def hmcode(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
     '''
 
     # Checks
-    if not util.is_array_monotonic(-zs):
+    if verbose: t_start = time()
+    if not util.is_array_monotonic(-np.array(zs)):
         raise ValueError('Redshift must be monotonically decreasing')
 
     # Halo mass range
@@ -107,8 +109,9 @@ def hmcode(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
         # Parameters of the linear spectrum pertaining to non-linear growth
         Rnl = _get_nonlinear_radius(R[0], R[-1], dc, iz, CAMB_results, cold=True) # Non-linear Lagrangian radius
         sigma8 = _get_sigmaR(8., iz, CAMB_results, cold=True)                     # RMS in the linear cold matter field at 8 Mpc/h
-        sigmaV = cosmology.sigmaV(Pk=lambda k: Pk_lin_interp(z, k))               # RMS in the linear displacement field
+        sigmaV = cosmology.sigmaV(lambda k: Pk_lin_interp(z, k))                  # RMS in the linear displacement field
         neff = _get_effective_index(Rnl, R, sigmaM)                               # Effective index of spectrum at collapse scale
+
         if verbose:
             print('Non-linear Lagrangian radius: {:.4} Mpc/h'.format(Rnl))
             print('RMS in matter field at 8 Mpc/h: {:.4}'.format(sigma8))
@@ -159,6 +162,10 @@ def hmcode(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
         Pk_HMcode[iz, :] = Pk_hm
 
     # Finish
+    if verbose:
+        t_finish = time()
+        print('HMcode predictions complete for {:} redshifts'.format(len(zs)))
+        print('Total HMcode run time: {:.3f}s'.format(t_finish-t_start))
     return Pk_HMcode