Merge pull request #760 from LSSTDESC/k_NL_integral

Implemented integral at C level to compute scale k_NL for a non-linear cut.

CHANGELOG.md

@@ -23,6 +23,8 @@
 - Improved implementation of halo model quantities (n(M), b(M), c(M)) (#668, #655, #656, #657, #636).
 - Add straightforward single massive neutrino option (#730).
 - Updated MG support via mu / Sigma (scale-independent) with higher accuracy implementation (#737)
+- Added function to estimate the scale at which non-linear structure formation becomes 
+  important for the computation of the matter power spectrum (#760).
 
 ## C library
 - Deprecated old halo model (#736)

benchmarks/data/codes/kNL.txt

@@ -0,0 +1,11 @@
+# [0] a, [1] k_NL Mpc^-1
+1.000000000000000056e-01 9.589032014324436748e-01
+1.291549665014883885e-01 7.428033944455602056e-01
+1.668100537200058742e-01 5.757310448836822081e-01
+2.154434690031883370e-01 4.467698285641908407e-01
+2.782559402207124277e-01 3.475697112659149601e-01
+3.593813663804627523e-01 2.718037161313307526e-01
+4.641588833612778631e-01 2.147341521854057300e-01
+5.994842503189409255e-01 1.727980028153486558e-01
+7.742636826811269968e-01 1.431192890249111394e-01
+1.000000000000000000e+00 1.230450029799866479e-01

benchmarks/data/codes/kNL_benchmark.ipynb

@@ -0,0 +1 @@
+{"cells":[{"metadata":{"trusted":true},"cell_type":"code","source":"import pyccl as ccl\nimport numpy as np\nimport scipy.integrate\nimport math\nfrom functools import partial","execution_count":13,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cosmoin = ccl.Cosmology(\n    Omega_c=0.25,\n    Omega_b=0.05,\n    h=0.7,\n    sigma8=0.8,\n    n_s=0.96,\n    Neff=0,\n    m_nu=0.0,\n    w0=-1.,\n    wa=0.,\n    T_CMB=2.7,\n    m_nu_type='normal',\n    Omega_g=0,\n    Omega_k=0,\n    transfer_function='bbks',\n    matter_power_spectrum='linear')","execution_count":23,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"k_lin=np.logspace(-4., 3., 10000)\na_arr = np.logspace(-1., 0., 10)\nknl_arr = np.zeros((len(a_arr),))\nfor i in range(len(a_arr)):\n    a = a_arr[i]\n    #I will need linear power at z\n    pk_lin_z=ccl.linear_matter_power(cosmoin, k_lin, a)\n    interp_pk_lin_z=scipy.interpolate.interp1d(k_lin,pk_lin_z)\n\n\n    #---kNL prediction----\n    [intknl,errintknl]=scipy.integrate.quad(interp_pk_lin_z,min(k_lin),max(k_lin),epsabs=0,epsrel=1e-6)\n    print(\"Check error is small:\",errintknl/intknl)\n    knlemin2=1./(6*math.pi**2)*intknl\n    knl=1./math.sqrt(knlemin2)\n    print('knl[1/Mpc]=',knl)\n    knl_arr[i] = knl","execution_count":44,"outputs":[{"name":"stdout","output_type":"stream","text":"Check error is small: 1.5286495983617028e-07\nknl[1/Mpc]= 0.12304500297998665\n"}]},{"metadata":{"trusted":true},"cell_type":"code","source":"header_knl=\"[0] a, [1] k_NL Mpc^-1\"\n\nnp.savetxt(\"kNL.txt\", np.transpose(np.vstack((a_arr ,knl_arr))), header=header_knl)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.6.10","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":2}

benchmarks/data/codes/kNL_benchmark.py

@@ -0,0 +1,44 @@
+import pyccl as ccl
+import numpy as np
+import scipy.integrate
+import math
+from functools import partial
+
+cosmoin=ccl.Cosmology(
+    Omega_c=0.25,
+    Omega_b=0.05,
+    h=0.7,
+    sigma8=0.8,
+    n_s=0.96,
+    Neff=0,
+    m_nu=0.0,
+    w0=-1.,
+    wa=0.,
+    T_CMB=2.7,
+    m_nu_type='normal',
+    Omega_g=0,
+    Omega_k=0,
+    transfer_function='bbks',
+    matter_power_spectrum='linear')
+
+k_lin=np.logspace(-4., 3., 10000)
+a_arr = np.logspace(-1., 0., 10)
+knl_arr = np.zeros((len(a_arr),))
+for i in range(len(a_arr)):
+    a = a_arr[i]
+    #I will need linear power at z
+    pk_lin_z=ccl.linear_matter_power(cosmoin, k_lin, a)
+    interp_pk_lin_z=scipy.interpolate.interp1d(k_lin,pk_lin_z)
+
+
+    #---kNL prediction----
+    [intknl,errintknl]=scipy.integrate.quad(interp_pk_lin_z,min(k_lin),max(k_lin),epsabs=0,epsrel=1e-6)
+    print("Check error is small:",errintknl/intknl)
+    knlemin2=1./(6*math.pi**2)*intknl
+    knl=1./math.sqrt(knlemin2)
+    print('knl[1/Mpc]=',knl)
+    knl_arr[i] = knl
+
+header_knl="[0] a, [1] k_NL Mpc^-1"
+
+np.savetxt("kNL.txt", np.transpose(np.vstack((a_arr ,knl_arr))), header=header_knl)

benchmarks/data/kNL.txt

@@ -0,0 +1,11 @@
+# [0] a, [1] k_NL Mpc^-1
+1.000000000000000056e-01 9.589032014324436748e-01
+1.291549665014883885e-01 7.428033944455602056e-01
+1.668100537200058742e-01 5.757310448836822081e-01
+2.154434690031883370e-01 4.467698285641908407e-01
+2.782559402207124277e-01 3.475697112659149601e-01
+3.593813663804627523e-01 2.718037161313307526e-01
+4.641588833612778631e-01 2.147341521854057300e-01
+5.994842503189409255e-01 1.727980028153486558e-01
+7.742636826811269968e-01 1.431192890249111394e-01
+1.000000000000000000e+00 1.230450029799866479e-01

benchmarks/test_kNL.py

@@ -0,0 +1,31 @@
+import numpy as np
+import pyccl as ccl
+
+KNL_TOLERANCE = 1.0e-5
+
+
+def test_kNL():
+    cosmo = ccl.Cosmology(
+        Omega_c=0.25,
+        Omega_b=0.05,
+        h=0.7,
+        sigma8=0.8,
+        n_s=0.96,
+        Neff=0,
+        m_nu=0.0,
+        w0=-1.,
+        wa=0.,
+        T_CMB=2.7,
+        m_nu_type='normal',
+        Omega_g=0,
+        Omega_k=0,
+        transfer_function='bbks',
+        matter_power_spectrum='linear')
+
+    data = np.loadtxt('./benchmarks/data/kNL.txt')
+    a = data[:, 0]
+    kNL = data[:, 1]
+    kNL_ccl = ccl.kNL(cosmo, a)
+    for i in range(len(a)):
+        err = np.abs(kNL_ccl[i]/kNL[i] - 1)
+        assert np.allclose(err, 0, rtol=0, atol=KNL_TOLERANCE)

doc/0000-ccl_note/authors.csv

@@ -69,4 +69,5 @@ Valogiannis, Georgios, Georgios Valogiannis, Contributor, "Department of Astrono
 Villarreal,Antonio,Antonio Villarreal,Contributor,"Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260","Contributed to initial benchmarking, halo mass function code, and general code and issues review.",asv13@pitt.edu
 Vrastil,Michal,Michal Vrastil,Contributor,"Institute of Physics CAS, Prague, 182 21, CZ","Wrote documentation and example code, reviewed code.",vrastil@fzu.cz
 Wang,Kuan,Kuan Wang,Contributor,"Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, USA","Wrote halo profiles code.",kuw8@pitt.edu
+Wright,Bill,Bill Wright,Contributor,"School of Physics and Astronomy, Queen Mary University of London, London, E1 4NS","Implemented computation of the scale for making a non-linear scale cut.",w.wright@qmul.ac.uk
 Zuntz,Joe,Joe Zuntz,Contributor,"Institute for Astronomy, Royal Observatory Edinburgh, Edinburgh EH9 3HJ, UK","Wrote initial infrastructure, C testing setup, and reviewed code.",joezuntz@googlemail.com

doc/0000-ccl_note/main.bib

@@ -1319,7 +1319,7 @@ @article{baldauf12
 	Volume = 86,
 	Year = 2012,
 	Bdsk-Url-1 = {http://dx.doi.org/10.1103/PhysRevD.86.083540}}
-    
+
 @ARTICLE{desjacques18rev,
    author = {{Desjacques}, V. and {Jeong}, D. and {Schmidt}, F.},
     title = "{Large-scale galaxy bias}",
@@ -1349,3 +1349,48 @@ @article{bernardeau02
 	Volume = 367,
 	Year = 2002,
 	Bdsk-Url-1 = {http://dx.doi.org/10.1016/S0370-1573(02)00135-7}}
+
+@ARTICLE{cosmicvisions2018,
+       author = {{Cosmic Visions 21 cm Collaboration} and {Ansari}, R{\'e}za and
+         {Arena}, Evan J. and {Bandura}, Kevin and {Bull}, Philip and
+         {Castorina}, Emanuele and {Chang}, Tzu-Ching and {Chen}, Shi-Fan and
+         {Connor}, Liam and {Foreman}, Simon and {Frisch}, Josef and
+         {Green}, Daniel and {Johnson}, Matthew C. and {Karagiannis}, Dionysios and
+         {Liu}, Adrian and {Masui}, Kiyoshi W. and {Meerburg}, P. Daniel and
+         {M{\"u}nchmeyer}, Moritz and {Newburgh}, Laura B. and
+         {Obuljen}, Andrej and {O'Connor}, Paul and {Padmanabhan}, Hamsa and
+         {Shaw}, J. Richard and {Sheehy}, Christopher and {Slosar}, An{\v{z}}e and
+         {Smith}, Kendrick and {Stankus}, Paul and {Stebbins}, Albert and
+         {Timbie}, Peter and {Villaescusa-Navarro}, Francisco and
+         {Wallisch}, Benjamin and {White}, Martin},
+        title = "{Inflation and Early Dark Energy with a Stage II Hydrogen Intensity Mapping Experiment}",
+      journal = {arXiv e-prints},
+     keywords = {Astrophysics - Cosmology and Nongalactic Astrophysics, Astrophysics - Instrumentation and Methods for Astrophysics, High Energy Physics - Experiment},
+         year = 2018,
+        month = oct,
+          eid = {arXiv:1810.09572},
+        pages = {arXiv:1810.09572},
+archivePrefix = {arXiv},
+       eprint = {1810.09572},
+ primaryClass = {astro-ph.CO},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2018arXiv181009572C},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+@ARTICLE{white2015,
+       author = {{White}, Martin},
+        title = "{Reconstruction within the Zeldovich approximation}",
+      journal = {\mnras},
+     keywords = {gravitation, galaxies: haloes, galaxies: statistics, cosmological parameters, large-scale structure of Universe, Astrophysics - Cosmology and Nongalactic Astrophysics},
+         year = 2015,
+        month = jul,
+       volume = {450},
+       number = {4},
+        pages = {3822-3828},
+          doi = {10.1093/mnras/stv842},
+archivePrefix = {arXiv},
+       eprint = {1504.03677},
+ primaryClass = {astro-ph.CO},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2015MNRAS.450.3822W},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}

doc/0000-ccl_note/main.tex

@@ -424,13 +424,13 @@ \subsubsection{CAMB}
 \newpage
 
 \subsubsection{ISiTGR}
-\ccl can call the {\tt ISiTGR} software patch \citep{isitgr1,isitgr2} that is 
-built on the top of the {\tt CAMB} package to include modified gravity models. 
-This is invoked by specifying the transfer function option boltzmann\_isitgr. 
-\ccl then then calculates the GR or modified gravity linear matter power spectrum 
-at a given redshift. Here also, on initializing the cosmology object, we construct a 
-bi-dimensional spline in $k$ and the scale-factor which is then evaluated by the 
-relevant routines to obtain the matter power spectrum at the desired wavenumber 
+\ccl can call the {\tt ISiTGR} software patch \citep{isitgr1,isitgr2} that is
+built on the top of the {\tt CAMB} package to include modified gravity models.
+This is invoked by specifying the transfer function option boltzmann\_isitgr.
+\ccl then then calculates the GR or modified gravity linear matter power spectrum
+at a given redshift. Here also, on initializing the cosmology object, we construct a
+bi-dimensional spline in $k$ and the scale-factor which is then evaluated by the
+relevant routines to obtain the matter power spectrum at the desired wavenumber
 and redshift.
 
 \subsubsection{Halofit}
@@ -547,25 +547,25 @@ \subsubsection{Modified gravity ($\mu, \Sigma$)}
 scale-independent deviations from GR which arise at late times (currently with
 functional form given in equation \ref{muSigform}). Under these conditions,
 and the specified cosmology, the modified gravity matter power spectrum is
-calculated using {\tt ISiTGR}  (described further above) according to the values of the modified 
-gravity parameters. If these are zero or not specified then the GR matter power 
+calculated using {\tt ISiTGR}  (described further above) according to the values of the modified
+gravity parameters. If these are zero or not specified then the GR matter power
 spectrum is calculated.
 
-This implementation replaces the intilal one where the modified gravity power 
-spectrum is computed by re-scaling the GR matter power spectrum by 
-multiplying it by the squared ratio of the MG  growth factor $D_{MG}$ to the GR 
-one, $D_{GR}$. 
+This implementation replaces the intilal one where the modified gravity power
+spectrum is computed by re-scaling the GR matter power spectrum by
+multiplying it by the squared ratio of the MG  growth factor $D_{MG}$ to the GR
+one, $D_{GR}$.
 
-The new implementation solved accuracy issues for spectra and  correlations 
-encoutered in the initial implementation and provides the power  spectrum and 
-the 2D angular spectra and correlations at the same level of accuracy 
-as that of the LCDM standard model.   
+The new implementation solved accuracy issues for spectra and  correlations
+encoutered in the initial implementation and provides the power  spectrum and
+the 2D angular spectra and correlations at the same level of accuracy
+as that of the LCDM standard model.
 
-The splines for the power spectrum are then set as usual, but we note that 
-the $\mu, \Sigma$ MG implementation is valid only in the quasistatic limit 
-and is also restricted to the linear regime for now since no-reliable methods 
-are available to model such non-linear small scales nor the very 
-large scales where the quasistatic limit cease to be valid.   
+The splines for the power spectrum are then set as usual, but we note that
+the $\mu, \Sigma$ MG implementation is valid only in the quasistatic limit
+and is also restricted to the linear regime for now since no-reliable methods
+are available to model such non-linear small scales nor the very
+large scales where the quasistatic limit cease to be valid.
 
 \subsubsection{Extrapolation for the nonlinear power spectrum}
 \label{sec:NLextrapol}
@@ -686,6 +686,14 @@ \subsubsection{Normalization of the power spectrum}
 used, only the option to specify $\sigma_8$ is supported. The computation of
 $\sigma_8$ is described in Section \ref{sec:hmf}.
 
+\subsubsection{Non-linear scale cut}
+\label{sec:kNL}
+
+A conservative estimate of the scale at which non-linear structure formation becomes important for the computation of the matter power spectrum at a particular redshift can be obtained from the inverse of the RMS displacement $\sigma_{\eta}^2$ in the Zel'dovich approximation \cite{cosmicvisions2018}. Using the expression for $\sigma_{\eta}^2$ from \cite{white2015} we have
+\begin{equation}
+k_{\rm NL}(z) \approx \sigma_{\eta}(z)^{-1}=\left[ \frac{1}{6\pi^2} \int^{\infty}_0 dk P_{\rm lin}(k,z) \right]^{-1/2}~.
+\end{equation}
+
 \subsection{Nonlinear Power Spectra}
 \label{sec:Nonlinear_Pk}
 
@@ -745,7 +753,7 @@ \subsubsection{Intrinsic alignments}
 &+ (C_{1\delta,1} C_{2,2} + C_{1\delta,2} C_{2,1})D_{0E|E2}~,\\
 \label{eq:BBtot}
 P_{BB} =& C_{1\delta,1} C_{1\delta,2} A_{0B|0B} + C_{2,1} C_{2,2} A_{B2|B2} + (C_{1\delta,1} C_{2,2} + C_{1\delta,2} C_{2,1}) D_{0B|B2} ~,
-\end{align} 
+\end{align}
 where, as above, $P_{d1d1}$ corresponds to the chosen nonlinear matter power spectrum,
 and the other terms are one-loop PT expressions.
 

include/ccl_core.h

@@ -164,6 +164,8 @@ typedef struct ccl_gsl_params {
   double INTEGRATION_DISTANCE_EPSREL;
   // sigma_R integral
   double INTEGRATION_SIGMAR_EPSREL;
+  // k_NL integral
+  double INTEGRATION_KNL_EPSREL;
 
   // Root finding
   double ROOT_EPSREL;

include/ccl_power.h

@@ -110,6 +110,17 @@ double ccl_sigmaV(ccl_cosmology *cosmo, double R, double a, int * status);
  */
 double ccl_sigma8(ccl_cosmology *cosmo, int * status);
 
+/**
+ * Scale for the non-linear cut.
+ * Returns k_NL for specified cosmology at specified scale factor.
+ * @param cosmo Cosmology parameters and configurations
+ * @param a scale factor
+ * @param status Status flag. 0 if there are no errors, nonzero otherwise.
+ * For specific cases see documentation for ccl_error.c
+ * @return kNL.
+ */
+double ccl_kNL(ccl_cosmology *cosmo, double a, int * status);
+
 CCL_END_DECLS
 
 #endif

pyccl/__init__.py

@@ -27,9 +27,9 @@
 # Generalized power spectra
 from .pk2d import Pk2D
 
-# Power spectrum calculations and sigma8
+# Power spectrum calculations, sigma8 and kNL
 from .power import linear_matter_power, nonlin_matter_power, sigmaR, \
-    sigmaV, sigma8, sigmaM
+    sigmaV, sigma8, sigmaM, kNL
 
 # BCM stuff
 from .bcm import bcm_model_fka

pyccl/ccl_power.i

@@ -7,6 +7,7 @@
 // Enable vectorised arguments for arrays
 %apply (double* IN_ARRAY1, int DIM1) {(double* k, int nk)};
 %apply (double* IN_ARRAY1, int DIM1) {(double* R, int nR)};
+%apply (double* IN_ARRAY1, int DIM1) {(double* a, int na)};
 %apply (int DIM1, double* ARGOUT_ARRAY1) {(int nout, double* output)};
 
 %include "../include/ccl_power.h"
@@ -60,5 +61,26 @@ void sigmaV_vec(ccl_cosmology * cosmo, double a, double* R, int nR,
 
 %}
 
+
+/* The python code here will be executed before all of the functions that
+   follow this directive. */
+%feature("pythonprepend") %{
+    if numpy.shape(a) != (nout,):
+        raise CCLError("Input shape for `a` must match `(nout,)`!")
+%}
+
+%inline %{
+
+void kNL_vec(ccl_cosmology * cosmo,
+        double* a,  int na,
+        int nout, double* output, int *status) {
+    assert(nout == na);
+    for(int i=0; i < na; i++){
+      output[i] = ccl_kNL(cosmo, a[i], status);
+    }
+}
+
+%}
+
 /* The directive gets carried between files, so we reset it at the end. */
 %feature("pythonprepend") %{ %}

pyccl/power.py

@@ -1,5 +1,5 @@
 from . import ccllib as lib
-from .pyutils import _vectorize_fn2
+from .pyutils import _vectorize_fn2, _vectorize_fn
 import numpy as np
 from .core import check
 
@@ -108,3 +108,21 @@ def sigma8(cosmo):
     """
     cosmo.compute_linear_power()
     return sigmaR(cosmo, 8.0 / cosmo['h'])
+
+
+def kNL(cosmo, a):
+    """Scale for the non-linear cut.
+
+    .. note:: k_NL is calculated based on Lagrangian perturbation theory as the
+              inverse of the variance of the displacement field, i.e.
+              k_NL = 1/sigma_eta = [1/(6 pi^2) * int P_L(k) dk]^{-1/2}.
+
+    Args:
+        cosmo (:class:`~pyccl.core.Cosmology`): Cosmological parameters.
+        a (float or array_like): Scale factor(s), normalized to 1 today.
+
+    Returns:
+        float or array-like: Scale of non-linear cut-off; Mpc^-1.
+    """
+    cosmo.compute_linear_power()
+    return _vectorize_fn(lib.kNL, lib.kNL_vec, cosmo, a)

pyccl/tests/test_power.py

@@ -129,6 +129,17 @@ def test_sigma8_consistent():
     assert np.allclose(ccl.sigmaR(COSMO, 8 / COSMO['h'], 1), COSMO['sigma8'])
 
 
+@pytest.mark.parametrize('A', [
+    1,
+    1.0,
+    [0.3, 0.5, 1],
+    np.array([0.3, 0.5, 1])])
+def test_kNL(A):
+    knl = ccl.kNL(COSMO, A)
+    assert np.all(np.isfinite(knl))
+    assert np.shape(knl) == np.shape(A)
+
+
 @pytest.mark.parametrize('tf,pk,m_nu', [
     # ('boltzmann_class', 'emu', 0.06), - this case is slow and not needed
     (None, 'emu', 0.06),

pyccl/tests/test_swig_interface.py

@@ -198,6 +198,14 @@ def test_swig_power():
             3,
             status)
 
+    assert_raises(
+        CCLError,
+        ccllib.kNL_vec,
+        COSMO,
+        [0.5, 1.0],
+        3,
+        status)
+
 
 def test_swig_haloprofile():
     status = 0