diff --git a/benchmarks/data/pt_bm_bacco_lbias_z0.txt b/benchmarks/data/pt_bm_bacco_lbias_z0.txt new file mode 100644 index 000000000..e6241e858 --- /dev/null +++ b/benchmarks/data/pt_bm_bacco_lbias_z0.txt @@ -0,0 +1,70 @@ +1.000000000000000048e-04 6.855035589030412957e+06 1.731322507686917903e+06 +1.125017868369810811e-04 6.305171508418726735e+06 1.582594949291167548e+06 +1.265665204151353231e-04 5.785967353791039437e+06 1.442792104142579250e+06 +1.423895970044196770e-04 5.296550671375462785e+06 1.311645761624566978e+06 +1.601908408999486275e-04 4.836049007399951108e+06 1.188887711121035740e+06 +1.802175583616277257e-04 4.403589908092526719e+06 1.074249742015381111e+06 +2.027479733508104082e-04 3.998300919680489227e+06 9.674636436908420874e+05 +2.280950927954279619e-04 3.619309588390721940e+06 8.682612055304173846e+05 +2.566110550823265716e-04 3.265743460452281870e+06 7.763742169183804654e+05 +2.886920221888471748e-04 2.936730082092395518e+06 6.915344672377756797e+05 +3.247836834182669485e-04 2.631396999538830481e+06 6.134737458721853327e+05 +3.653874472005138197e-04 2.348871759019009769e+06 5.419238422048285138e+05 +4.110674069786092737e-04 2.088281906761171995e+06 4.766165456194913713e+05 +4.624581779553800964e-04 1.848754988992715254e+06 4.172836454994115629e+05 +5.202737135735483647e-04 1.629418551941404119e+06 3.636569312280857121e+05 +5.853172242133589456e-04 1.429400141835035989e+06 3.154681921889922814e+05 +6.584923359046477350e-04 1.247827304901328404e+06 2.724492177655681735e+05 +7.408156440773042140e-04 1.083827587368055945e+06 2.343317973413098371e+05 +8.334308367548572648e-04 9.365285354628558271e+05 2.008477202996404667e+05 +9.376245833996174270e-04 8.050576954135099659e+05 1.717287760240163479e+05 +1.054844410147369503e-03 6.885426134477542946e+05 1.467067538979064557e+05 +1.186718809765804124e-03 5.861108357931714272e+05 1.255134433046632039e+05 +1.335079865717084041e-03 4.968899086777408957e+05 1.078806336279096286e+05 +1.501988704632487232e-03 4.200073783289733692e+05 9.354011425096351013e+04 +1.689764130801174212e-03 3.545907909746699152e+05 8.222367455732729286e+04 +1.901014840481703464e-03 2.997676928425814258e+05 7.366310393046194804e+04 +2.138675663578101793e-03 2.546656301604587352e+05 6.759019175383858965e+04 +2.406048336173026957e-03 2.184121491559549759e+05 6.373672741086473252e+04 +2.706847370356106099e-03 1.901347960568578856e+05 6.183450028503002977e+04 +3.045251658600457711e-03 1.689611170908559870e+05 6.161529975976930291e+04 +3.425962529608314513e-03 1.540186584857124253e+05 6.281091521855384053e+04 +3.854269062174795271e-03 1.444349664691298094e+05 6.515313604482416122e+04 +4.336121564451593177e-03 1.393375872688434902e+05 6.837375162203848595e+04 +4.878214239431705873e-03 1.378540671125804947e+05 7.220455133365021902e+04 +5.488078185096710929e-03 1.391119522280704114e+05 7.637732456311455462e+04 +6.174186021244368154e-03 1.422387888430423627e+05 8.062386069388614851e+04 +6.946069596539016129e-03 1.463621231852223282e+05 8.467594910941817216e+04 +7.814452411046684047e-03 1.506095014823389938e+05 8.826537919316513580e+04 +8.791398593953060961e-03 1.541084699621199979e+05 9.112394032858090941e+04 +9.890480506158424306e-03 1.562512351973984623e+05 9.303400141059812449e+04 +1.112696729619151899e-02 1.570776950108758174e+05 9.388496726643192233e+04 +1.251803702898198097e-02 1.560964727966866922e+05 9.339143319992923352e+04 +1.408301533451966799e-02 1.530083412248911336e+05 9.138866302858100971e+04 +1.584364389186067659e-02 1.478099369707513251e+05 8.778123262820804666e+04 +1.782438247843147125e-02 1.406713264526355197e+05 8.271700102806974610e+04 +2.005274878089318216e-02 1.318069295894066454e+05 7.639847539796734054e+04 +2.255970068843577217e-02 1.218050505921091826e+05 6.925515076792147011e+04 +2.538006637956496600e-02 1.115507349460478581e+05 6.193308795112917142e+04 +2.855302817742248062e-02 1.018606175297661539e+05 5.513927801614269265e+04 +3.212266689566698613e-02 9.343733965296312817e+04 4.933937925436145451e+04 +3.613857423731676316e-02 8.667749094105955737e+04 4.492700872193000396e+04 +4.065654175439026990e-02 8.118331846965101431e+04 4.146428609593502188e+04 +4.573933593981235196e-02 7.589140657157383976e+04 3.826747317999914230e+04 +5.145757021965836897e-02 6.968134210939257173e+04 3.432066468152480229e+04 +5.789068596000986056e-02 6.218479290794593544e+04 2.943168996150249950e+04 +6.512805611719643673e-02 5.468021221455463092e+04 2.454961468407437860e+04 +7.327022686403775187e-02 4.844426069067248318e+04 2.091755571477606281e+04 +8.243031444155220211e-02 4.383208375223904295e+04 1.867093857859003037e+04 +9.273557664208829932e-02 3.990856214840315806e+04 1.664717399621781078e+04 +1.043291807559273976e-01 3.563047891117057588e+04 1.425023996776120293e+04 +1.173721925428021368e-01 3.185891197906968591e+04 1.226188214779403825e+04 +1.320458138603942677e-01 2.915394661601902044e+04 1.099675414611043561e+04 +1.485539000363775841e-01 2.643317288311504672e+04 9.710702582257510585e+03 +1.671257919569474892e-01 2.400396424968933570e+04 8.614309544499872572e+03 +1.880195022170215446e-01 2.193334960257961939e+04 7.742739802899634014e+03 +2.115252995961465210e-01 2.007347703997483040e+04 6.967392010316811138e+03 +2.379697416579423519e-01 1.846920023339546606e+04 6.327384991037761210e+03 +2.677202114965329116e-01 1.708514059014760278e+04 5.777041934382074032e+03 +3.011900216573443756e-01 1.589338777559453592e+04 5.309448138988474966e+03 +3.388441561392025458e-01 1.486986077699188718e+04 4.904235139092390455e+03 diff --git a/benchmarks/data/pt_bm_bacco_lbias_z1.txt b/benchmarks/data/pt_bm_bacco_lbias_z1.txt new file mode 100644 index 000000000..6e81475cc --- /dev/null +++ b/benchmarks/data/pt_bm_bacco_lbias_z1.txt @@ -0,0 +1,70 @@ +1.000000000000000048e-04 -5.787538347034783801e+05 -5.760465749232481467e+05 +1.125017868369810811e-04 -5.399584957598655019e+05 -5.363730526740099303e+05 +1.265665204151353231e-04 -5.027504697382528102e+05 -4.984428364051764365e+05 +1.423895970044196770e-04 -4.670977722293614643e+05 -4.622170264950743876e+05 +1.601908408999486275e-04 -4.329684188237880589e+05 -4.276567233218124602e+05 +1.802175583616277257e-04 -4.003304251122052083e+05 -3.947230272636098671e+05 +2.027479733508104082e-04 -3.691518066854010685e+05 -3.633770386987907114e+05 +2.280950927954279619e-04 -3.394005791342977900e+05 -3.335798580057563377e+05 +2.566110550823265716e-04 -3.110447580495520961e+05 -3.052925855627732235e+05 +2.886920221888471748e-04 -2.840523590216873563e+05 -2.784763217479078448e+05 +3.247836834182669485e-04 -2.583913976416455116e+05 -2.530921669395852659e+05 +3.653874472005138197e-04 -2.340298895001568599e+05 -2.291012215160999040e+05 +4.110674069786092737e-04 -2.109358501878337993e+05 -2.064645858555978339e+05 +4.624581779553800964e-04 -1.890772952954104985e+05 -1.851433603363852890e+05 +5.202737135735483647e-04 -1.684222404137553240e+05 -1.650986453368044167e+05 +5.853172242133589456e-04 -1.489387011335245916e+05 -1.462915412350882252e+05 +6.584923359046477350e-04 -1.305946930453853856e+05 -1.286831484094547632e+05 +7.408156440773042140e-04 -1.133582317401100008e+05 -1.122345672382037301e+05 +8.334308367548572648e-04 -9.719733280848342110e+04 -9.690689809963956941e+04 +9.376245833996174270e-04 -8.208001184116789955e+04 -8.266124137198489916e+04 +1.054844410147369503e-03 -6.797428442893552710e+04 -6.945869743353693048e+04 +1.186718809765804124e-03 -5.484816616253430402e+04 -5.726036666257845354e+04 +1.335079865717084041e-03 -4.266967263262737106e+04 -4.602734943732642569e+04 +1.501988704632487232e-03 -3.140681943000727551e+04 -3.572074613609736844e+04 +1.689764130801174212e-03 -2.102762214539703564e+04 -2.630165713715392849e+04 +1.901014840481703464e-03 -1.150009636950433014e+04 -1.773118281874665990e+04 +2.138675663578101793e-03 -2.792257693070999721e+03 -9.970423559155035036e+03 +2.406048336173026957e-03 5.127878293163939816e+03 -2.980479736651973781e+03 +2.706847370356106099e-03 1.229229599847857935e+04 3.277548270499342379e+03 +3.045251658600457711e-03 1.873297983213460975e+04 8.842560084024857133e+03 +3.425962529608314513e-03 2.448191420340646437e+04 1.375345532565876601e+04 +3.854269062174795271e-03 2.957108352155875764e+04 1.804913361712914411e+04 +4.336121564451593177e-03 3.403247219586047140e+04 2.176849458016718927e+04 +4.878214239431705873e-03 3.789806463558069663e+04 2.495043783650359183e+04 +5.488078185096710929e-03 4.119984524998840061e+04 2.763386300786942229e+04 +6.174186021244368154e-03 4.396979844835220865e+04 2.985766971599531462e+04 +6.946069596539016129e-03 4.623990863994021493e+04 3.166075758261149167e+04 +7.814452411046684047e-03 4.804216023402116116e+04 3.308202622944871837e+04 +8.791398593953060961e-03 4.940853763986354897e+04 3.416037527823755954e+04 +9.890480506158424306e-03 5.034065830618274776e+04 3.490851339100255427e+04 +1.112696729619151899e-02 5.071714106207974692e+04 3.523463415399403311e+04 +1.251803702898198097e-02 5.043450745832359098e+04 3.505944959891346662e+04 +1.408301533451966799e-02 4.941743802003123710e+04 3.432452220551073697e+04 +1.584364389186067659e-02 4.762556696446789283e+04 3.299835495910294412e+04 +1.782438247843147125e-02 4.509617248215829750e+04 3.111436246673823553e+04 +2.005274878089318216e-02 4.190430272682777286e+04 2.872576903427655270e+04 +2.255970068843577217e-02 3.832294702413521009e+04 2.605077494000550723e+04 +2.538006637956496600e-02 3.459748322807742079e+04 2.326410160910756531e+04 +2.855302817742248062e-02 3.110758490381195224e+04 2.067548204131027524e+04 +3.212266689566698613e-02 2.810810286775376881e+04 1.846723907944000166e+04 +3.613857423731676316e-02 2.576981973068406296e+04 1.677732809927342532e+04 +4.065654175439026990e-02 2.393682114865565381e+04 1.547656622730936397e+04 +4.573933593981235196e-02 2.224548790820713839e+04 1.429676701664213215e+04 +5.145757021965836897e-02 2.012452510844651988e+04 1.278113062405563323e+04 +5.789068596000986056e-02 1.743829033173220159e+04 1.083803095979534010e+04 +6.512805611719643673e-02 1.462507993024422467e+04 8.802395076095741388e+03 +7.327022686403775187e-02 1.252046963015262190e+04 7.316186139403020206e+03 +8.243031444155220211e-02 1.128591522795189303e+04 6.518127387715445366e+03 +9.273557664208829932e-02 1.013164123800368179e+04 5.779984140146171740e+03 +1.043291807559273976e-01 8.633057393865612539e+03 4.757156171246198028e+03 +1.173721925428021368e-01 7.322863391620325274e+03 3.886577755426259500e+03 +1.320458138603942677e-01 6.576447003548580142e+03 3.434331748271198649e+03 +1.485539000363775841e-01 5.697144779268932325e+03 2.885211845994848773e+03 +1.671257919569474892e-01 4.946280733547519048e+03 2.429306315553752029e+03 +1.880195022170215446e-01 4.359654143806421416e+03 2.096614094600254248e+03 +2.115252995961465210e-01 3.807726728444853507e+03 1.783181421621180107e+03 +2.379697416579423519e-01 3.353523843515297358e+03 1.537292786463995981e+03 +2.677202114965329116e-01 2.956444373898936192e+03 1.327618530470946098e+03 +3.011900216573443756e-01 2.620988530749326856e+03 1.154202125870107466e+03 +3.388441561392025458e-01 2.334633412620704803e+03 1.010505217137959335e+03 diff --git a/benchmarks/test_ptpk_bacco_lbias.py b/benchmarks/test_ptpk_bacco_lbias.py new file mode 100644 index 000000000..4239e75a5 --- /dev/null +++ b/benchmarks/test_ptpk_bacco_lbias.py @@ -0,0 +1,57 @@ +import os +import numpy as np +import pyccl as ccl +import pyccl.nl_pt as pt +import pytest + +LPTPK_TOLERANCE = 1e-4 + +# Set cosmology +COSMO = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, + h=0.7, sigma8=0.8, n_s=0.96) + +# Redshifts +zs = np.array([0., 1.]) + +# Tracers +ptt = {} +ptt['g'] = pt.PTNumberCountsTracer(b1=1.3, b2=1.5, bs=1.7, bk2=0.1) +ptt['m'] = pt.PTMatterTracer() + +# Calculator +a_arr = 1./(1+np.array([0., 0.25, 0.5, 0.75, 1.]))[::-1] +ptc = pt.BaccoLbiasCalculator(log10k_min=-4, log10k_max=-0.47, + nk_per_decade=20, a_arr=a_arr) + +ptc.update_ingredients(COSMO) + +# Read data +dirdat = os.path.join(os.path.dirname(__file__), 'data') +data = [] +data.append(np.loadtxt(os.path.join(dirdat, 'pt_bm_bacco_lbias_z0.txt'), + unpack=True)) +data.append(np.loadtxt(os.path.join(dirdat, 'pt_bm_bacco_lbias_z1.txt'), + unpack=True)) +order = ['gg', 'gm'] + + +@pytest.mark.parametrize('comb', enumerate(order)) +def test_pt_pk(comb): + i_d, cc = comb + t1, t2 = cc + + ptt1 = ptt[t1] + ptt2 = ptt[t2] + + pk = ptc.get_biased_pk2d(ptt1, tracer2=ptt2) + + for iz, z in enumerate(zs): + a = 1./(1+z) + kin = data[iz][0] + ind = np.where((kin < 1.0) & (kin > 1E-3)) + k = kin[ind] + dpk = data[iz][i_d+1][ind] + tpk = pk(k, a, COSMO) + assert np.all( + np.fabs(tpk / dpk - 1) < LPTPK_TOLERANCE + ) diff --git a/pyccl/nl_pt/__init__.py b/pyccl/nl_pt/__init__.py index b7327d228..81e5e96ab 100644 --- a/pyccl/nl_pt/__init__.py +++ b/pyccl/nl_pt/__init__.py @@ -1,3 +1,4 @@ from .tracers import * from .ept import * from .lpt import * +from .bacco_lbias import * diff --git a/pyccl/nl_pt/bacco_lbias.py b/pyccl/nl_pt/bacco_lbias.py new file mode 100644 index 000000000..7ba04e761 --- /dev/null +++ b/pyccl/nl_pt/bacco_lbias.py @@ -0,0 +1,456 @@ +__all__ = ("BaccoLbiasCalculator",) + +import numpy as np + +from .. import (CCLAutoRepr, CCLError, Pk2D, + get_pk_spline_a, unlock_instance) + + +# All valid Pk pair labels and their aliases +# Note that bacco lbias has no b3nl, so terms +# with b3nl are automatically set to zero. +# TODO: make this common to all nl_pt models. +_PK_ALIAS = { + 'm:m': 'm:m', 'm:b1': 'm:b1', 'm:b2': 'm:b2', + 'm:b3nl': 'zero', 'm:bs': 'm:bs', 'm:bk2': 'm:bk2', + 'b1:b1': 'b1:b1', 'b1:b2': 'b1:b2', 'b1:b3nl': 'zero', + 'b1:bs': 'b1:bs', 'b1:bk2': 'b1:bk2', 'b2:b2': 'b2:b2', + 'b2:b3nl': 'zero', 'b2:bs': 'b2:bs', 'b2:bk2': 'b2:bk2', + 'b3nl:b3nl': 'zero', 'b3nl:bs': 'zero', + 'b3nl:bk2': 'zero', 'bs:bs': 'bs:bs', + 'bs:bk2': 'bs:bk2', 'bk2:bk2': 'bk2:bk2'} + + +class BaccoLbiasCalculator(CCLAutoRepr): + """ This class implements a set of methods that can be + used to compute the various components needed to estimate + hybrid Lagrangian bias expansion correlations using the + emulator baccoemu + (https://bitbucket.org/rangulo/baccoemu/src/master/baccoemu/). + + This is a hybrid model, featuring a second order + Lagrangian bias expansion coupled with advecting the + Lagrangian fields to Eulerian observables through + N-body simulations. It has been tested to reproduce the + galaxy-galaxy and galaxy-matter power spectra down + to scales of 0.7 h/Mpc. + + In the parametrisation used here, the bias function + in Lagrangian coordinates is expanded as (ignoring constant + terms): + + .. math:: + w_{\\rm g}(\\boldsymbol{q})=1 + b_1\\,\\delta+(b_2/2)\\,\\delta^2+ + (b_s/2)\\,s^2+(b_{k2}/2)\\nabla^2\\delta. + + This translates to Eulerian space + + .. math:: + \\delta(\\boldsymbol{x}) = \\int \\mathrm{d}^3\\boldsymbol{q} + w_{\\rm g}(\\boldsymbol{q}) \\delta^{\\rm D}(\\boldsymbol{x} - + \\boldsymbol{q} - \\boldsymbol{\\Psi}), + + where the displacement field :math:`\\Psi` is obtained from + simulations. + + .. note:: This calculator does not account for any form of + stochastic bias contribution to the power spectra. + If necessary, consider adding it in post-processing. + + Args: + cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object. + If present, internal PT power spectrum templates will + be initialized. If ``None``, you will need to initialize + them using the :meth:`update_ingredients` method. + log10k_min (:obj:`float`): decimal logarithm of the minimum + Fourier scale (in :math:`{\\rm Mpc}^{-1}`) for which you want to + calculate perturbation theory quantities. + log10k_max (:obj:`float`): decimal logarithm of the maximum + Fourier scale (in :math:`{\\rm Mpc}^{-1}`) for which you want to + calculate perturbation theory quantities. + nk_per_decade (:obj:`int` or :obj:`float`): number of k values per + decade. + a_arr (array): array of values of the scale factor at + which all power spectra will be evaluated. If ``None``, + the default sampling used internally by CCL will be + used. Note that this may be slower than a bespoke sampling + optimised for your particular application. + k_cutoff (:obj:`float`): exponential cutoff scale. All power + spectra will be multiplied by a cutoff factor of the + form :math:`\\exp(-(k/k_*)^n)`, where :math:`k_*` is + the cutoff scale. This may be useful when using the + resulting power spectra to compute correlation + functions if some of the PT contributions do not + fall sufficiently fast on small scales. If ``None`` + (default), no cutoff factor will be applied. + n_exp_cutoff (:obj:`float`): exponent of the cutoff factor (see + ``k_cutoff``). + """ + __repr_attrs__ = __eq_attrs__ = ('k_s', 'a_s', 'exp_cutoff') + + def __init__(self, *, cosmo=None, + log10k_min=-4, log10k_max=-0.47, nk_per_decade=20, + a_arr=None, k_cutoff=None, n_exp_cutoff=4): + # Load emulator + import warnings + with warnings.catch_warnings(): + warnings.filterwarnings('ignore', category=UserWarning) + import baccoemu + self.emu = baccoemu.Lbias_expansion() + + # k sampling + nk_total = int((log10k_max - log10k_min) * nk_per_decade) + self.k_s = np.logspace(log10k_min, log10k_max, nk_total) + + # a sampling + if a_arr is None: + a_arr = get_pk_spline_a().copy() + mask = a_arr >= self.emu.emulator['nonlinear']['bounds'][-1][0] + self.a_s = a_arr[mask] + self.z_s = 1/self.a_s-1 + + # Cutoff factor + if k_cutoff is not None: + self.exp_cutoff = np.exp(-(self.k_s/k_cutoff)**n_exp_cutoff) + self.exp_cutoff = self.exp_cutoff[None, :] + else: + self.exp_cutoff = 1 + + # Initialize all expensive arrays to ``None``. + self._cosmo = None + + # Fill them out if cosmo is present + if cosmo is not None: + self.update_ingredients(cosmo) + + # All valid Pk pair labels + self._pk_valid = list(_PK_ALIAS.keys()) + # List of Pk2Ds to fill out + self._pk2d_temp = {} + + def _check_init(self): + if self.initialised: + return + raise CCLError("PT templates have not been initialised " + "for this calculator. Please do so using " + "`update_ingredients`.") + + @property + def initialised(self): + return hasattr(self, "lpt_table") + + @unlock_instance + def update_ingredients(self, cosmo): + """ Update the internal PT arrays. + + Args: + cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object. + """ + if self.initialised and (cosmo == self._cosmo): + return + + # hubble + h = cosmo['h'] + + # get bacco parameters + emupars = dict( + omega_cold=cosmo['Omega_c'] + cosmo['Omega_b'], + omega_baryon=cosmo['Omega_b'], + ns=cosmo['n_s'], + hubble=cosmo['h'], + neutrino_mass=np.sum(cosmo['m_nu']), + w0=cosmo['w0'], + wa=cosmo['wa'], + expfactor=self.a_s + ) + + # if cosmo contains sigma8, we use it for baccoemu, otherwise we pass + # A_s to the emulator + if np.isnan(cosmo['A_s']): + sigma8tot = cosmo['sigma8'] + sigma8cold = self._sigma8tot_2_sigma8cold(emupars, sigma8tot) + emupars['sigma8_cold'] = sigma8cold + else: + emupars['A_s'] = cosmo['A_s'] + + # call bacco + k = self.k_s / h + lpt_table = self.emu.get_nonlinear_pnn(k=k, **emupars)[1] + lpt_table /= h ** 3 + + # save templates in a table + self.lpt_table = lpt_table + + # Reset template power spectra + self._pk2d_temp = {} + self._cosmo = cosmo + + def _sigma8tot_2_sigma8cold(self, emupars, sigma8tot): + """Use baccoemu to convert sigma8 total matter to sigma8 cdm+baryons + """ + if hasattr(emupars['omega_cold'], '__len__'): + _emupars = {} + for pname in emupars: + _emupars[pname] = emupars[pname][0] + else: + _emupars = emupars + A_s_fid = 2.1e-9 + sigma8tot_fid = self.emu.matter_powerspectrum_emulator.get_sigma8( + cold=False, A_s=A_s_fid, **_emupars) + A_s = (sigma8tot / sigma8tot_fid)**2 * A_s_fid + return self.emu.matter_powerspectrum_emulator.get_sigma8(cold=True, + A_s=A_s, + **_emupars) + + def _get_pgg(self, tr1, tr2): + """ Get the number counts auto-spectrum at the internal + set of wavenumbers and scale factors. + + Args: + tr1 (:class:`~pyccl.nl_pt.tracers.PTTracer`): first + tracer to correlate. + tr2 (:class:`~pyccl.nl_pt.tracers.PTTracer`): first + tracer to correlate. + + Returns: + array: 2D array of shape `(N_a, N_k)`, where `N_k` \ + is the size of this object's `k_s` attribute, and \ + `N_a` is the size of the object's `a_s` attribute. + """ + self._check_init() + # Get biases + b11 = tr1.b1(self.z_s) + b21 = tr1.b2(self.z_s) + bs1 = tr1.bs(self.z_s) + bk21 = tr1.bk2(self.z_s) + b12 = tr2.b1(self.z_s) + b22 = tr2.b2(self.z_s) + bs2 = tr2.bs(self.z_s) + bk22 = tr2.bk2(self.z_s) + + # Transform from Eulerian to Lagrangian biases + bL11 = b11 - 1 + bL12 = b12 - 1 + # Get Pk templates + Pdmdm = self.lpt_table[0, :, :] # 1 1 + Pdmd1 = self.lpt_table[1, :, :] # 1 d + Pdmd2 = self.lpt_table[2, :, :] * 0.5 # 1 d2 + Pdms2 = self.lpt_table[3, :, :] * 0.5 # 1 s2 + Pdmk2 = self.lpt_table[4, :, :] * 0.5 # 1 lap + Pd1d1 = self.lpt_table[5, :, :] # d d + Pd1d2 = self.lpt_table[6, :, :] * 0.5 # d d2 + Pd1s2 = self.lpt_table[7, :, :] * 0.5 # d s2 + Pd1k2 = self.lpt_table[8, :, :] * 0.5 # d k2 + Pd2d2 = self.lpt_table[9, :, :] * 0.25 # d2 d2 + Pd2s2 = self.lpt_table[10, :, :] * 0.25 # d2 s2 + Pd2k2 = self.lpt_table[11, :, :] * 0.25 # d2 k2 + Ps2s2 = self.lpt_table[12, :, :] * 0.25 # s2 s2 + Ps2k2 = self.lpt_table[13, :, :] * 0.25 # s2 k2 + Pk2k2 = self.lpt_table[14, :, :] * 0.25 # k2 k2 + + pgg = (Pdmdm + + (bL11+bL12)[:, None] * Pdmd1 + + (bL11*bL12)[:, None] * Pd1d1 + + (b21 + b22)[:, None] * Pdmd2 + + (bs1 + bs2)[:, None] * Pdms2 + + (bL11*b22 + bL12*b21)[:, None] * Pd1d2 + + (bL11*bs2 + bL12*bs1)[:, None] * Pd1s2 + + (b21*b22)[:, None] * Pd2d2 + + (b21*bs2 + b22*bs1)[:, None] * Pd2s2 + + (bs1*bs2)[:, None] * Ps2s2 + + (bk21 + bk22)[:, None] * Pdmk2 + + (bL12 * bk21 + bL11 * bk22)[:, None] * Pd1k2 + + (b22 * bk21 + b21 * bk22)[:, None] * Pd2k2 + + (bs2 * bk21 + bs1 * bk22)[:, None] * Ps2k2 + + (bk21 * bk22)[:, None] * Pk2k2) + + return pgg*self.exp_cutoff + + def _get_pgm(self, trg): + """ Get the number counts - matter cross-spectrum at the internal + set of wavenumbers and scale factors. + + Args: + trg (:class:`~pyccl.nl_pt.tracers.PTTracer`): number + counts tracer. + + Returns: + array: 2D array of shape `(N_a, N_k)`, where `N_k` \ + is the size of this object's `k_s` attribute, and \ + `N_a` is the size of the object's `a_s` attribute. + """ + self._check_init() + # Get biases + b1 = trg.b1(self.z_s) + b2 = trg.b2(self.z_s) + bs = trg.bs(self.z_s) + bk2 = trg.bk2(self.z_s) + + # Compute Lagrangian bias + bL1 = b1 - 1 + # Get Pk templates + Pdmdm = self.lpt_table[0, :, :] # 1 1 + Pdmd1 = self.lpt_table[1, :, :] # 1 d + Pdmd2 = self.lpt_table[2, :, :] * 0.5 # 1 d2 + Pdms2 = self.lpt_table[3, :, :] * 0.5 # 1 s2 + Pdmk2 = self.lpt_table[4, :, :] * 0.5 # 1 k2 + + pgm = (Pdmdm + + bL1[:, None] * Pdmd1 + + b2[:, None] * Pdmd2 + + bs[:, None] * Pdms2 + + bk2[:, None] * Pdmk2) + + return pgm*self.exp_cutoff + + def _get_pmm(self): + """ Get the one-loop matter power spectrum. + + Returns: + array: 2D array of shape `(N_a, N_k)`, where `N_k` \ + is the size of this object's `k_s` attribute, and \ + `N_a` is the size of the object's `a_s` attribute. + """ + self._check_init() + + pk = self.lpt_table[0, :, :] + + return pk*self.exp_cutoff + + def get_biased_pk2d(self, tracer1, *, tracer2=None, + extrap_order_lok=1, extrap_order_hik=2): + """Returns a :class:`~pyccl.pk2d.Pk2D` object containing + the PT power spectrum for two quantities defined by + two :class:`~pyccl.nl_pt.tracers.PTTracer` objects. + + Args: + tracer1 (:class:`~pyccl.nl_pt.tracers.PTTracer`): the first + tracer being correlated. + tracer2 (:class:`~pyccl.nl_pt.tracers.PTTracer`): the second + tracer being correlated. If ``None``, the auto-correlation + of the first tracer will be returned. + extrap_order_lok (:obj:`int`): extrapolation order to be used on + k-values below the minimum of the splines. See + :class:`~pyccl.pk2d.Pk2D`. + extrap_order_hik (:obj:`int`): extrapolation order to be used on + k-values above the maximum of the splines. See + :class:`~pyccl.pk2d.Pk2D`. + + Returns: + :class:`~pyccl.pk2d.Pk2D`: PT power spectrum. + """ + if tracer2 is None: + tracer2 = tracer1 + + t1 = tracer1.type + t2 = tracer2.type + + if t1 == 'IA' or t2 == 'IA': + raise ValueError("Intrinsic alignments not implemented in " + "BaccoLbiasCalculator.") + + if t1 == 'NC': + if t2 == 'NC': + pk = self._get_pgg(tracer1, tracer2) + else: # Must be matter + pk = self._get_pgm(tracer1) + else: # Must be matter + if t2 == 'NC': + pk = self._get_pgm(tracer2) + else: # Must be matter + pk = self._get_pmm() + + pk2d = Pk2D(a_arr=self.a_s, + lk_arr=np.log(self.k_s), + pk_arr=pk, + is_logp=False, + extrap_order_lok=extrap_order_lok, + extrap_order_hik=extrap_order_hik) + return pk2d + + def get_pk2d_template(self, kind, *, extrap_order_lok=1, + extrap_order_hik=2): + """Returns a :class:`~pyccl.pk2d.Pk2D` object containing + the power spectrum template for two of the PT operators. The + combination returned is determined by ``kind``, which must be + a string of the form ``'q1:q2'``, where ``q1`` and ``q2`` denote + the two operators whose power spectrum is sought. Valid + operator names are: ``'m'`` (matter overdensity), ``'b1'`` + (first-order overdensity), ``'b2'`` (:math:`\\delta^2` + term in galaxy bias expansion), ``'bs'`` (:math:`s^2` term + in galaxy bias expansion), ``'bk2'`` (non-local + :math:`\\nabla^2 \\delta` term in galaxy bias expansion) + + Args: + kind (:obj:`str`): string defining the pair of PT operators for + which we want the power spectrum. + extrap_order_lok (:obj:`int`): extrapolation order to be used on + k-values below the minimum of the splines. See + :class:`~pyccl.pk2d.Pk2D`. + extrap_order_hik (:obj:`int`): extrapolation order to be used on + k-values above the maximum of the splines. See + :class:`~pyccl.pk2d.Pk2D`. + + Returns: + :class:`~pyccl.pk2d.Pk2D`: PT power spectrum. + """ + if not (kind in _PK_ALIAS): + # Reverse order and check again + kind_reverse = ':'.join(kind.split(':')[::-1]) + if not (kind_reverse in _PK_ALIAS): + raise ValueError(f"Pk template {kind} not valid") + kind = kind_reverse + pk_name = _PK_ALIAS[kind] + + # If already built, return + if pk_name in self._pk2d_temp: + return self._pk2d_temp[pk_name] + + self._check_init() + + if pk_name == 'm:m': + pk = self._get_pmm() + elif pk_name == 'm:b1': + pk = self.lpt_table[1, :, :] + elif pk_name == 'm:b2': + pk = self.lpt_table[2, :, :] * 0.5 + elif pk_name == 'm:bs': + pk = self.lpt_table[3, :, :] * 0.5 + elif pk_name == 'm:bk2': + pk = self.lpt_table[4, :, :] * 0.5 + elif pk_name == 'b1:b1': + pk = self.lpt_table[5, :, :] + elif pk_name == 'b1:b2': + pk = self.lpt_table[6, :, :] * 0.5 + elif pk_name == 'b1:bs': + pk = self.lpt_table[7, :, :] * 0.5 + elif pk_name == 'b1:bk2': + pk = self.lpt_table[8, :, :] * 0.5 + elif pk_name == 'b2:b2': + pk = self.lpt_table[9, :, :] * 0.25 + elif pk_name == 'b2:bs': + pk = self.lpt_table[10, :, :] * 0.25 + elif pk_name == 'b2:bk2': + pk = self.lpt_table[11, :, :] * 0.25 + elif pk_name == 'bs:bs': + pk = self.lpt_table[12, :, :] * 0.25 + elif pk_name == 'bs:bk2': + pk = self.lpt_table[13, :, :] * 0.25 + elif pk_name == 'bk2:bk2': + pk = self.lpt_table[14, :, :] * 0.25 + elif pk_name == 'zero': + # If zero, store None and return + self._pk2d_temp[pk_name] = None + return None + + # Build interpolator + pk2d = Pk2D(a_arr=self.a_s, + lk_arr=np.log(self.k_s), + pk_arr=pk, + is_logp=False, + extrap_order_lok=extrap_order_lok, + extrap_order_hik=extrap_order_hik) + + # Store and return + self._pk2d_temp[pk_name] = pk2d + return pk2d diff --git a/pyccl/tests/test_bacco_lbias_power.py b/pyccl/tests/test_bacco_lbias_power.py new file mode 100644 index 000000000..21c727c9c --- /dev/null +++ b/pyccl/tests/test_bacco_lbias_power.py @@ -0,0 +1,245 @@ +import numpy as np +import pyccl as ccl +import pytest + + +COSMO = ccl.Cosmology( + Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.8, n_s=0.96, + transfer_function='bbks', matter_power_spectrum='linear') +TRS = {'TG': ccl.nl_pt.PTNumberCountsTracer(b1=2.0, b2=2.0, bs=2.0, bk2=2.0), + 'TM': ccl.nl_pt.PTMatterTracer()} +PTC = ccl.nl_pt.LagrangianPTCalculator(cosmo=COSMO) + + +def test_bacco_lbias_calculator_smoke(): + c = ccl.nl_pt.BaccoLbiasCalculator(log10k_min=-3, + log10k_max=1, + nk_per_decade=10) + assert len(c.k_s) == 40 + + +@pytest.mark.parametrize('tr1,tr2', + [['TG', 'TG'], + ['TG', 'TM'], + ['TM', 'TG'], + ['TM', 'TM']]) +def test_bacco_lbias_get_pk2d_smoke(tr1, tr2): + t2 = None if tr2 == tr1 else TRS[tr2] + ptc = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + + pk = ptc.get_biased_pk2d(TRS[tr1], tracer2=t2) + assert isinstance(pk, ccl.Pk2D) + + +def test_bacco_lbias_get_pk2d_nl(): + ptc = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + pk = ptc.get_biased_pk2d(TRS['TG']) + assert isinstance(pk, ccl.Pk2D) + + +@pytest.mark.parametrize('typ_nlin,typ_nloc', [('linear', 'nonlinear'), + ('nonlinear', 'linear')]) +def test_bacco_lbias_k2pk_types(typ_nlin, typ_nloc): + tg = ccl.nl_pt.PTNumberCountsTracer(1., 0., 0., bk2=1.) + tm = ccl.nl_pt.PTNumberCountsTracer(1., 0., 0.) + ptc1 = ccl.nl_pt.LagrangianPTCalculator( + b1_pk_kind=typ_nlin, bk2_pk_kind=typ_nloc, cosmo=COSMO) + ptc2 = ccl.nl_pt.LagrangianPTCalculator( + b1_pk_kind=typ_nloc, cosmo=COSMO) + pkmm = ptc1.get_biased_pk2d(tm, tracer2=tm) + pkmm2 = ptc2.get_biased_pk2d(tm, tracer2=tm) + pkgg = ptc1.get_biased_pk2d(tg, tracer2=tg) + ks = np.geomspace(1E-3, 1E1, 128) + p1 = pkgg(ks, 1., cosmo=COSMO) + p2 = pkmm(ks, 1., cosmo=COSMO)+ks**2*pkmm2(ks, 1., cosmo=COSMO) + assert np.allclose(p1, p2, atol=0, rtol=1E-4) + + +@pytest.mark.parametrize('kind', ['m:m', 'm:b1', 'm:b2', 'm:bs', 'm:bk2', + 'b1:b1', 'b1:b2', 'b1:bs', 'b1:bk2', 'b2:b2', + 'b2:bs', 'b2:bk2', 'bs:bs', 'bs:bk2', + 'bk2:bk2', 'b1:b3nl']) +def test_bacco_lbias_deconstruction(kind): + ptc = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + b_nc = ['b1', 'b2', 'bs', 'bk2', 'b3nl'] + pk1 = ptc.get_pk2d_template(kind) + + def get_tr(tn): + if tn == 'm': + return ccl.nl_pt.PTMatterTracer() + if tn in b_nc: + bdict = {b: 0.0 for b in b_nc} + # This is b1E = 1, so that internally it goes to b1L = 0 + if tn == 'b1': + bdict[tn] = 2.0 + else: + bdict['b1'] = 1.0 + bdict[tn] = 1.0 + return ccl.nl_pt.PTNumberCountsTracer( + b1=bdict['b1'], b2=bdict['b2'], + bs=bdict['bs'], bk2=bdict['bk2'], + b3nl=bdict['b3nl']) + + tn1, tn2 = kind.split(':') + t1 = get_tr(tn1) + t2 = get_tr(tn2) + + pkmm = ptc.get_pk2d_template('m:m') + pkx1 = ptc.get_pk2d_template(f'm:{tn1}') + pkx2 = ptc.get_pk2d_template(f'm:{tn2}') + + pk2 = ptc.get_biased_pk2d(t1, tracer2=t2) + if pk1 is None: + assert np.allclose(pk2(0.2, 1.0, cosmo=COSMO), + pkmm(0.2, 1.0, cosmo=COSMO) + + pkx1(0.2, 1.0, cosmo=COSMO), + atol=0, rtol=1e-6) + else: + if (t1.type == 'M') & (t2.type == 'M'): + v1 = pk1(0.2, 1.0, cosmo=COSMO) + elif (t1.type == 'M') | (t2.type == 'M'): + v1 = pkmm(0.2, 1.0, cosmo=COSMO) + pk1(0.2, 1.0, cosmo=COSMO) + else: + v1 = pkmm(0.2, 1.0, cosmo=COSMO) + pkx1(0.2, 1.0, cosmo=COSMO) \ + + pkx2(0.2, 1.0, cosmo=COSMO) + pk1(0.2, 1.0, cosmo=COSMO) + v2 = pk2(0.2, 1.0, cosmo=COSMO) + assert np.allclose(v1, v2, atol=0, rtol=1e-6) + # Check cached + pk3 = ptc._pk2d_temp[ccl.nl_pt.lpt._PK_ALIAS[kind]] + chached_pkmm = ptc._pk2d_temp[ccl.nl_pt.lpt._PK_ALIAS['m:m']] + chached_pkx1 = ptc._pk2d_temp[ccl.nl_pt.lpt._PK_ALIAS[f'm:{tn1}']] + chached_pkx2 = ptc._pk2d_temp[ccl.nl_pt.lpt._PK_ALIAS[f'm:{tn2}']] + if (t1.type == 'M') & (t2.type == 'M'): + v3 = pk3(0.2, 1.0, cosmo=COSMO) + elif (t1.type == 'M') | (t2.type == 'M'): + v3 = chached_pkmm(0.2, 1.0, cosmo=COSMO) \ + + pk3(0.2, 1.0, cosmo=COSMO) + else: + v3 = chached_pkmm(0.2, 1.0, cosmo=COSMO) \ + + chached_pkx1(0.2, 1.0, cosmo=COSMO) \ + + chached_pkx2(0.2, 1.0, cosmo=COSMO) \ + + pk3(0.2, 1.0, cosmo=COSMO) + assert np.allclose(v1, v3, atol=0, rtol=1e-6) + + +def test_bacco_lbias_pk_cutoff(): + # Tests the exponential cutoff + ks = np.geomspace(1E-2, 0.3, 30) + + t = ccl.nl_pt.PTNumberCountsTracer(b1=1.0) + ptc1 = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + pk2d1 = ptc1.get_biased_pk2d(t, tracer2=t) + pk1 = pk2d1(ks, 1.0, cosmo=COSMO) + ptc2 = ccl.nl_pt.BaccoLbiasCalculator(k_cutoff=0.1, + n_exp_cutoff=2., + cosmo=COSMO) + pk2d2 = ptc2.get_biased_pk2d(t, tracer2=t) + pk2 = pk2d2(ks, 1.0, cosmo=COSMO) + + expcut = np.exp(-(ks/0.1)**2) + assert np.allclose(pk1*expcut, pk2, atol=0, rtol=1E-2) + + +def test_bacco_lbias_matter_1loop(): + # Check P(k) for linear tracer with b1=1 is the same + # as matter P(k) + ptc = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + tg = ccl.nl_pt.PTNumberCountsTracer(b1=1.0) + tm = ccl.nl_pt.PTMatterTracer() + ks = np.geomspace(1E-3, 10, 64) + pk1 = ptc.get_biased_pk2d(tm)(ks, 1.0, cosmo=COSMO) + pk2 = ptc.get_biased_pk2d(tg)(ks, 1.0, cosmo=COSMO) + assert np.allclose(pk1, pk2, atol=0, rtol=1E-3) + + +def test_bacco_lbias_calculator_raises(): + # Uninitialized templates + with pytest.raises(ccl.CCLError): + ptc = ccl.nl_pt.BaccoLbiasCalculator() + ptc.get_biased_pk2d(TRS['TG']) + + # TODO: Discuss this test + # Wrong pair combination + with pytest.raises(ValueError): + ptc = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + ptc.get_pk2d_template('b1:b3') + + with pytest.raises(ValueError): + t1 = ccl.nl_pt.PTIntrinsicAlignmentTracer(c1=1) + ptc = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + ptc.get_biased_pk2d(t1, tracer2=t1) + + +def test_bacco_lbias_template_swap(): + # Test that swapping operator order gets you the same + # Pk + ks = np.array([0.01, 0.1, 1.0]) + ptc = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO) + pk1 = ptc.get_pk2d_template('b2:bs')(ks, 1.0, cosmo=COSMO) + pk2 = ptc.get_pk2d_template('bs:b2')(ks, 1.0, cosmo=COSMO) + assert np.all(pk1 == pk2) + + +def test_bacco_lbias_eq(): + ptc1 = ccl.nl_pt.BaccoLbiasCalculator() + # Should be the same + ptc2 = ccl.nl_pt.BaccoLbiasCalculator() + assert ptc1 == ptc2 + # Should still be the same + ptc2 = ccl.nl_pt.BaccoLbiasCalculator( + a_arr=ccl.pyutils.get_pk_spline_a()) + assert ptc1 == ptc2 + # Different a sampling + ptc2 = ccl.nl_pt.BaccoLbiasCalculator( + a_arr=np.linspace(0.5, 1., 30)) + assert ptc1 != ptc2 + # Should do nothing if cosmo is the same + ptc2 = ccl.nl_pt.BaccoLbiasCalculator( + cosmo=COSMO) + lpt_table_1 = ptc2.lpt_table + ptc2.update_ingredients(COSMO) + assert lpt_table_1 is ptc2.lpt_table + + +def test_bacco_lbias_sigma8_A_s(): + ks = np.logspace(-2, np.log10(0.3), 30) + + t = ccl.nl_pt.PTNumberCountsTracer(b1=1.0) + + COSMO_A_s = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, n_s=0.96, + A_s=2.1265e-9, m_nu=0.2) + COSMO_sigma8 = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, n_s=0.96, + sigma8=0.8059043572377348, m_nu=0.2) + + ptc1 = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO_A_s) + pk2d1 = ptc1.get_biased_pk2d(t, tracer2=t) + pk1 = pk2d1(ks, 1.0, cosmo=COSMO_A_s) + + ptc2 = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO_sigma8) + pk2d2 = ptc2.get_biased_pk2d(t, tracer2=t) + pk2 = pk2d2(ks, 1.0, cosmo=COSMO_sigma8) + + assert np.allclose(pk1, pk2, atol=0, rtol=1E-3) + + +def test_bacco_lbias_many_A_s(): + COSMO_sigma8 = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, n_s=0.96, + sigma8=0.81, m_nu=0.2) + pt = ccl.nl_pt.BaccoLbiasCalculator(cosmo=COSMO_sigma8) + emupars = dict( + omega_cold=COSMO_sigma8['Omega_c'] + COSMO_sigma8['Omega_b'], + omega_baryon=COSMO_sigma8['Omega_b'], + ns=COSMO_sigma8['n_s'], + hubble=COSMO_sigma8['h'], + neutrino_mass=np.sum(COSMO_sigma8['m_nu']), + w0=COSMO_sigma8['w0'], + wa=COSMO_sigma8['wa'], + expfactor=pt.a_s + ) + sigma8tot = COSMO_sigma8['sigma8'] + sigma8cold = pt._sigma8tot_2_sigma8cold(emupars, sigma8tot) + newemupars = {} + for key in emupars: + newemupars[key] = np.full(len(pt.a_s), emupars[key]) + sigma8cold_arr = pt._sigma8tot_2_sigma8cold(newemupars, sigma8tot) + assert np.allclose(np.mean(sigma8cold_arr), sigma8cold, atol=0, rtol=1E-4)