Fast analytical calculation of the random pair counts for realistic survey geometry

Galaxy clustering is a standard cosmological probe that is commonly analysed through two-point statistics. In observations, the es- timation of the two-point correlation function crucially relies on counting pairs in a random catalogue. The latter contains a large number of randomly distributed points, which accounts for the survey window function. Random pair counts can also be advanta- geously used for modelling the window function in the observed power spectrum. Since pair counting scales as O ( N 2 ), where N is the number of points, the computational time to measure random pair counts can be very expensive for large surveys. In this work, we present an alternative approach for estimating those counts that does not rely on the use of a random catalogue. We derived an analytical expression for the anisotropic random-random pair counts that accounts for the galaxy radial distance distribution, survey geometry, and possible galaxy weights. We show that a prerequisite is the estimation of the two-point correlation function of the angular selection function, which can be obtained e ﬃ ciently using pixelated angular maps. Considering the cases of the VIPERS and SDSS-BOSS redshift surveys, we ﬁnd that the analytical calculation is in excellent agreement with the pair counts obtained from random catalogues. The main advantage of this approach is that the primary calculation only takes a few minutes on a single CPU and it does not depend on the number of random points. Furthermore, it allows for an accuracy on the monopole equivalent to what we would otherwise obtain when using a random catalogue with about 1500 times more points than in the data at hand. We also describe and test an approximate expression for data-random pair counts that is less accurate than for random-random counts, but still provides subpercent accuracy on the monopole. The presented formalism should be very useful in accounting for the window function in next-generation surveys, which will necessitate accurate two-point window function estimates over huge observed cosmological volumes.


Introduction
The spatial distribution of galaxies has a long history of providing cosmological parameter constraints (e.g., Strauss et al. 1992;Vogeley et al. 1992;Maddox et al. 1996;Peacock et al. 2001;Cole et al. 2005;Tegmark et al. 2006;Percival et al. 2010;Blake et al. 2012;de la Torre et al. 2013;Alam et al. 2017;eBOSS Collaboration 2020, and references therein). This arises from the fact that the statistical properties of galaxies, particularly spatial ones, can be predicted by cosmological models. When analysing galaxy clustering, we usually compress the information by using summary statistics, the most natural one being the two-point correlation function or its Fourier counterpart in the power spectrum. This is due to the nearly Gaussian nature of primordial matter perturbations, which are almost fully described by their two-point statistics. Although gravitational evolution leads to non-Gaussianity, and in turn, non-vanishing higher-order n-point statistics, two-point statistics continues to be very informative.
Despite the cosmological principle that implies that the correlation function is isotropic, meaning that it is only a function of the norm of the separation vector, because of the way the line-of-sight distance is measured in redshift surveys and the presence of peculiar velocities, the observed correlation function becomes anisotropic. These velocities are induced on large scales by the coherent convergence of matter towards overdensities as part of the general process of structure growth. This anisotropy makes observed galaxy n-point statistics sensitive to the strength of gravity acting on the large-scale structure (Kaiser 1987;Guzzo et al. 2008).
Formally, the two-point correlation function is the excess probability of finding a pair of objects at a given distance, with respect to the expectation in a random Poisson distribution of points. In practice, we rely on statistical estimators to measure the correlation function from galaxy survey data. The first estimator was proposed by Peebles & Hauser (1974), taking the form ξ PH (s) = DD(s)/RR(s) − 1, where DD and RR are the normalised number of distinct pairs separated by a vector s, in the data and random samples, respectively. The latter sample is constructed such that random points follow the same radial and angular selection functions as the data. Other estimators were later proposed (Hewett 1982;Davis & Peebles 1983;Hamilton 1993) to reduce the estimation variance, notably induced by discreteness and boundary effects. In particular, the Landy & Szalay (1993) minimum-variance estimator was designed such that for any survey geometry, its variance is nearly Poissonian. This estimator, defined as makes use of additional data-random pairs DR. To estimate the correlation function, we therefore need to compute the number of pairs as a function of the separation. To avoid introducing bias in the estimator and to minimise variance, the random catalogue must be much larger than the data catalogue (Landy & Szalay 1993;Keihänen et al. 2019). We usually consider that taking at least about 20−50 times more random points than objects in the data is enough to avoid introducing additional variance (e.g., Samushia et al. 2012;de la Torre et al. 2013;Sánchez et al. 2017;Bautista et al. 2021). A problem is that the computational time for direct pair counting scales as O(N 2 ), with N the number of elements in a given sample. Nonetheless, the complexity can be reduced, at best, to O(N) using appropriate algorithms and various efficient codes that implement them have been developed (e.g., Moore et al. 2001;Jarvis et al. 2004;Alonso 2012;Hearin et al. 2017;Marulli et al. 2016;Slepian & Eisenstein 2016;Sinha & Garrison 2020). For the random-random pairs calculation specifically, additional strategies can be used to speed up the computation beyond parallelisation, such as splitting the random sample and averaging the counts over subsamples (Keihänen et al. 2019). Still, for large surveys, the computational time for estimating the correlation function, especially random-random pairs, can become an issue. In particular, next-generation surveys such as Euclid (Laureijs et al. 2011) or DESI (DESI Collaboration 2016), will necessitate random samples as large as about 3 × 10 8 objects in several redshift bins, or even larger if we consider a single, wideredshift bin (e.g., Mueller et al. 2019;Castorina et al. 2019).
The role of random-random pair counts in the correlation function estimator is to account for the survey selection function, that is, the effective observed volume and its impact on the data-data pair counts. In Fourier space instead, common estimators for the power spectrum (e.g., Feldman et al. 1994;Yamamoto et al. 2006) provide a direct estimate of the windowconvolved power spectrum. To be able to compare theoretical predictions to observations, it is necessary to convolve the model power spectrum with the survey window function. This convolution is computationally expensive in the likelihood analysis, but it can be done efficiently by performing a multiplication in configuration space, as proposed by Wilson et al. (2017). The latter shows that the window-convolved power spectrum multipoles momentsP (k) can be written as: where H denotes the Hankel transform, A q p are coefficients, ξ (s) are model correlation function multipole moments, N is a normalisation factor, RR (s) are the multipole moments of the random-random pair counts, and ∆s is the bin size in s (Wilson et al. 2017;Beutler et al. 2017).
Random-random pairs counts are a purely geometrical quantity that depends on cosmology only through the radial selection function, which is defined in terms of the radial comoving distance. In the case of a simple geometry, such as a cubical volume with constant number density and periodic boundary conditions, the RR pair counts can be predicted from the appropriately normalised ratio between the spherical shell volume at s and the total volume. In the case of a realistic survey geometry, and taking advantage of the fact that radial and angular selection functions are usually assumed to be uncorrelated, Demina et al. (2018) developed a semi-analytical method to compute the RR and DR pair counts along the directions parallel and transverse to the line of sight, but still using a random sample to account for angular correlations. In fact, the process of spectroscopic observation can break this assumption, making radial and angular selection functions partially correlated. Nonetheless, in practice, as in the case of fibre or slit collision, for instance, we can assume independence in building RR but while applying an object or pairwise weighting scheme to account for those correlations (e.g., de la Torre et al. 2013;Bianchi & Percival 2017;Ross et al. 2017).
In this paper, we provide general expressions for the anisotropic RR and DR pair counts in the case of a realistic survey geometry, including the cases for the different definitions of the pair's line of sight. We apply this formalism to the VIMOS Public Extragalactic Redshift Survey (VIPERS, Guzzo et al. 2014;Garilli et al. 2014) and Sloan Digital Sky Survey Baryon Oscillation Spectroscopic Survey (SDSS-BOSS, Eisenstein et al. 2011;Dawson et al. 2013) and we perform an assessment of the accuracy of the method. This paper is organised as follows. Section 2 presents the formalism for random-random and data-random pair counts. This formalism is applied and its accuracy assessed in Sect. 3. We conclude in Sect. 4.

Formalism
In this section, we provide the analytical formalism for the random-random and data-random pair counts.

Random-random pairs
In a survey sample where sources are selected in redshift, the number of sources in a given radial distance interval [r min , r max ] is: with p(r) the number of sources as function of the radial distance r and where W(θ, ϕ) is the survey angular selection function in spherical coordinates. The latter encodes the probability of observing a source at any angular position on the sky and takes values from 0 to 1. Here, n(r) is the source number density given by with W the angular selection function averaged over the full sky. We note that radial weights, such as those of Feldman et al. (1994), can be included straightforwardly in the p(r), such that this becomes a weighted radial distribution in the equations. The total number of observed sources is therefore: In RR(s), we correlate points at two different positions r 1 and r 2 and it is convenient to write A40, page 2 of 10 with s = r 2 − r 1 and µ = r 1 · s/r 1 s. RR is obtained by integrating the angular and radial selection functions first over (r 1 , θ, ϕ) and then over the volume defined by the separation (s,θ,φ) as: where θ 1 , ϕ 1 (θ 2 , ϕ 2 ) are the angular positions at r 1 (r 2 ). Let us definer 1 = r 1 /r 1 andr 2 = r 2 /r 2 , we have then and the correlation function of the angular selection function is In our case, since we auto-correlate randoms points, RR will only depend on the angular separation (for cross-correlations with different angular selection functions, we would need to keep the angular dependence). This means that we can write ω We note that in the absence of an angular mask, that is, when we evenly probe the full sky, ω(φ) = 1. In compiling these results, we find that: We note that we have implicitly assumed an 'end-point' definition for the pair line of sight, that is, for every separation, the direction of the line of sight coincides with that of r 1 . With this definition, we can just use (µ * min , µ * max ) = (µ min , µ max ) for the integral limits in Eq. (12). In the case of the 'mid-point' definition for the pair line of sight, where µ = r · s/rs with r = 1 2 (r 1 + r 2 ), we can use the same equation, but we need to change the integral limits for each {r 1 , s, µ} as: We note that for the end-point definition it is important to compute pairs with µ < 0 since the correlation function in that case is not symmetric by pair exchange. For applications involving random-random multipole moments directly, the latter can be defined as: with L the Legendre polynomial of order . For the end-point definition, µ † = µ, while for the mid-point definition we have with r(r 1 , s, µ) = r 1 1 + µs/r + s 2 /(2r) 2 . Finally, we note that in order to cross-correlate tracers with different radial selection functions but the same angular selection function, we can use the same formalism but with different n A (r 1 ) and n B (r 2 ) in Eqs. (12)-(14).

Data-random pairs
The Landy & Szalay (1993) estimator includes data-random pairs to minimise variance. A similar formalism to that used for random-random pair counts can be employed to evaluate datarandom pair counts. However, contrarily to the RR case, we now have to cross-correlate a discrete set of sources with a continuous random distribution. The discrete limit of Eq. (3) is with δ D the Dirac delta function, r = (r, θ, ϕ), and r i the ith source position in the data vector. We can then use the same methodology as in Sect. 2.1. To make the computation of the data-random pair counts tractable we make the assumption: withr 1 = (θ, ϕ) andr i the angular position in the data vector. Under this assumption, the angular correlation at a particular point is given by that of the whole sample. For a large-enough N and a sufficiently homogeneous angular sampling of the data, the approximation should hold. We find that: with where we integrate over all the data and angular selection function positions. In practice, we can cross-correlate a map containing the galaxy number density per pixel with the angular selection function map. Introducing a generic data weight w i , for each source and assuming that the angular cross-correlation function does not depend on the pair orientation, we obtain DR(s min , s max , µ min , µ max ) = We can already see that the calculation involves a sum over all data sources, which is potentially more computationally expensive to evaluate than in the RR and direct pair-counting cases.
Nonetheless, to make the computation efficient, we can approximate this by taking the continuous limit on the sum and write it out similarly as we do for the RR case: DR(s min , s max , µ min , µ max ) = 8π 2 r max r min where the angular cross-correlation function, ω DR , can be written as ω DR = W 1 W 2 , with W 1 as a pixelated map containing (weighted) source counts and W 2 as the angular selection function as in Sect. 2.1. In this case, the normalisation on W 1 and W 2 does not matter since the final result is proportional to W 1 W 2 W 1 W 2 . In principle, Eq. (21) should be close to Eq. (20) when the p(r) is fine enough to faithfully reproduce data radial overdensities.

Application
In this section, we apply our formalism for RR and DR pair counts to the case of BOSS and VIPERS redshift surveys and test its accuracy.

Numerical implementation
The method takes as input the radial distance distribution of sources and the correlation function of the angular selection function. To compute the latter, we first produce a Healpix (Górski et al. 2005) full-sky map from the survey angular mask, which is generally in form of a list of distinct spherical polygons with associated weights. We infer the angular correlation function from the maps using Polspice (Szapudi et al. 2001;Chon et al. 2004), which takes advantage of Healpix isolatitude pixelation scheme to quickly evaluate the angular correlation function ω(θ) with fast spherical harmonic transforms. The CPU time to compute the angular correlation function depends on the map resolution, but it is generally quite efficient. The map resolution is controlled by the nside 1 parameter. For instance, for BOSS it takes 3, 20 and 130 min on a single CPU for nside = 2048, 4096, and 8192, respectively. For VIPERS, it takes 6, 40, and 320 min for nside = 8192, 16384, and 32768, respectively. We note that we compute the angular power spectra until max = 3 × nside; it is possible to reduce the run time by reducing max . The main limitation is memory since this operation needs to load the full-sky map, which can be difficult for highresolution maps. Once the correlation function of the angular mask, average map value W , and average map squared value W 2 are computed, we can estimate the pair counts. This can be done by numerically evaluating the multi-dimensional integrals of Eqs. (12) and (20) (or (21)), respectively, for RR and DR. We tested different integration schemes, namely: gsl (Gough 2009) integration algorithm cquad; cuba library (Hahn 2005) set of optimised algorithms for multi-dimensional integration: vegas, suave, divonne and cuhre. In each case, we need to specify a maximum tolerance on the integral relative error ε. The gsl cquad algorithm only performs one-dimensional integrals, and we thus implement nested integrals with same ε to perform the full three-dimensional integral.
There are three potential sources of error associated to the analytical pair count calculation: the p(r) estimation, the ω(θ) estimation, and the precision on numerical integrals. In our methodology, the error level associated to each source is controlled respectively by the binning in p(r), angular map resolution nside, and ε. In the last case however, we note that different algorithms can yield slightly different results even if ε is very small.
Regarding the performances of the RR implementation, for the two considered surveys and considering 6000 bins in (s, µ), the full RR(s, µ) based on Eq. (12) takes about 5−20 min on a single CPU using the cuba library, even for ε ≈ 10 −6 (gsl takes significantly more time when ε < 10 −3 ). In principle, we gain an additional factor of two when using the pair line-of-sight midpoint definition, as in that case, there is a symmetry along the line of sight for auto-correlation and we only need to compute µ > 0 pairs. The run time depends on the number of bins in (s, µ) but also in principle, on the shape of the integrand, as for complex p(r) or ω(θ), the integrals will take more time to converge (although, in our case, we did not see any noticeable difference). For the RR multipole moments in Eq. (14), the calculation only takes about 5 s with cuba algorithms using 30 bins in s, independently of the adopted value of ε.
Regarding DR, the run times for the approximation in Eq. (21) are similar to those for RR by definition. However, the evaluation of Eq. (20) leads to large computational times of up to several weeks on a single CPU for large datasets. The run times scale linearly with the number of objects in the data sample. In that case, direct data-random pair counting might be more efficient.
The C code that follows this implementation is publicly available 2 . It can be used with any input p(r) and ω(θ) to predict RR or DR pair counts.

Survey selection functions
We consider two realistic redshift survey selection functions, those of SDSS-BOSS DR12 CMASS (Alam et al. 2015) and VIPERS PDR2 (Scodeggio et al. 2018) galaxy samples. We use the public galaxy catalogues 3 and associated angular masks. Those two samples have complementary properties and thus allow testing the method in different conditions. Indeed, while BOSS survey is wide and has a low galaxy number density, VIPERS is much narrower and denser. Each survey is composed of two separated fields on the sky but here we only consider the BOSS North Galactic Cap (NGC) and VIPERS W1 fields and we focus on 0.5 < z < 0.75 and 0.7 < z < 1.2 redshift intervals for BOSS and VIPERS, respectively.
The radial selection functions are shown in Fig. 1. The BOSS p(r) is estimated from the data by taking the histogram of galaxy comoving distances and cubic-spline interpolating between the bins. In the case of VIPERS, we use the fitting function for the redshift distribution given in de la Torre et al. (2013). We assumed two cosmologies to convert redshift to comoving distance: flat ΛCDM with Ω m = 0.31 and Ω m = 0.25 for BOSS and VIPERS, respectively. Nonetheless, the choice of fiducial cosmology has no impact on the accuracy of the analytical predictions.
The angular selection functions that we used for VIPERS W1 and BOSS NGC are presented in Appendix A. In total, the surface covered by the BOSS NGC (VIPERS W1) field is  We test different map resolutions by varying the Healpix resolution parameter nside from 2048 to 8192. For BOSS, the correlation function is very smooth. The relative difference between nside = 2048, 4096 cases and nside = 8192 is roughly constant, at 0.1% and 0.05% respectively. In the case of VIPERS, the angular mask has more small-scale features but similarly, the relative differences between angular correlation functions based on different map resolutions are nearly constant in scale. The bias is larger than in the BOSS case, with a relative difference with respect to nside = 65536 of 4%, 2%, and 0.5%, respectively for nside = 8192, 16834, 32768. This difference is due to the fact that we need a significantly higher map resolution to correctly account for the angular selection function as illustrated in Fig. 4. The latter figure shows a detail of the VIPERS angular mask pixelated at nside = 8192 and nside = 65536. Overall, the  convergence of the angular correlation function with increasing nside, allows us to assess that resolutions of nside = 8192 and nside = 65536 are sufficiently accurate for BOSS and VIPERS respectively. We note that while the map resolution impacts the estimation of the angular selection function two-point correlation function, it also changes the estimation of n(r) through W , which partly compensates the bias from ω(θ) in the final RR and DR estimations.

RR counts
We compare our analytical prediction for RR with the average random-random counts RR , obtained from 100 random samples constructed using the same radial and angular selection functions. Within the considered redshift intervals, there are 435185 BOSS and 24316 VIPERS galaxies and we generate 3 × 10 7 and 3.9 × 10 6 points per random sample, respectively (i.e. multiplicative factors of about 70 and 160 with respect to the data) with the radial distributions shown in Fig. 1. We compute the pair counts from the random samples using the fast Corrfunc pair-counting code (Sinha & Garrison 2020). Our method predicts anisotropic RR(s, µ) counts, but to simplify the comparison, we consider the first three even multipoles, that is, RR (s) = (2 + 1)/2 i RR(s, µ i )L (µ i )∆µ with = 0, 2, 4, where linear bins µ i extend from −1 to 1. Those comparisons are presented in Fig. 5 for BOSS and in Fig. 6 for VIPERS. We see that for both surveys, the relative difference between the analytical computation and RR is well within the variance of the random samples. We compare the results obtained with different numerical integration algorithms (see figure insets and Sect. 3.1) and find that cuhre tends to depart from the others algorithms, which is understandable since it is intrinsically different from the others. If we ignore cuhre, we see that, at most, the relative difference between the analytical computation and RR remains within 3 × 10 −5 for BOSS and 1.7 × 10 −4 for VIPERS on the monopole, and about 10 −3 for the quadrupole and hexadecapole for both surveys.
The variance on the random sample counts depends on the number of points in the sample and, thus, we may ask what is the number of random points needed to achieve the same accuracy as in the analytical method. Keihänen et al. (2019) showed that the relative variance on RR in a given bin is: (22) with N r the number of random points, and G p , G t terms are (Landy & Szalay 1993): with n p and n t the number of pairs and triplets averaged over several realisations. While G p can easily be estimated from the random samples, we directly solve for G t from the estimated var(RR). We can then deduce which N r give standard deviations similar to 3 × 10 −5 and 1.7 × 10 −4 for the monopole. We found that we need an additional factor of at least 20(10) for BOSS (VIPERS) in the number of random points. Therefore, the analytical method allows the achievement of the same accuracy as by using a random sample with about 20 × 70 (10 × 160) more points than data in BOSS (VIPERS). Finally, we note that cuba integration algorithms have parameters that can be potentially further fine-tuned to achieve better accuracy.

DR counts
In the DR case, we need to rely on approximations. Under the approximation in Eq. (17), we have two possibilities to calculate DR counts: either a discrete sum over all source distances as in Eq. (20) or by further approximating the discrete sum by an integral as in Eq. (21). In the last case, we can already anticipate that the results will depend on the input p(r 1 ), particularly its ability to reproduce line-of-sight structures in the data. In Figs. 7 and 8, we show different estimations of the data p(r 1 ) in VIPERS and BOSS, varying the bin size in r 1 . In the limit where p(r 1 ) resembles a sum of Dirac delta functions, Eq. (21) should be equivalent Fig. 5. Relative difference between the analytical and random catalogue-based mean RR pair-count multipole moments ( = 0, 2, 4, in the top, middle, and bottom panels, respectively) for BOSS. The grey shaded area shows the standard deviation among the random catalogues, while blue, orange, green, red, and purple curves present the relative differences obtained with vegas, suave, divonne, cuhre, and gsl algorithms, respectively, when using ε = 10 −5 .
to Eq. (20). For the random part we use in p(r 2 ) the distributions provided in Fig. 1. It is worth noting that we use cubic splines to model the data distributions in Figs. 7 and 8. While other ways of estimating p(r 1 ) could have been chosen, we only focus on the relative importance of the binning, and therefore, the method used for the estimation is irrelevant here (however, it would be necessary for an accurate, in-depth characterisation of the radial selection function).
Following the same methodology as in Sect. 3.3, we compute DR for both surveys using the same data catalogue and 100 random samples, which we later compare to the predictions based on Eqs. (20) and (21) using different input data p(r 1 ). In the case of VIPERS, we find that when using Eq. (20), the discrepancy between the analytical prediction and direct pair counting is on the order of 1% for the monopole, up to several percent  for the quadrupole and hexadecapole, as shown in Fig. 9. Moreover, we see that the prescription in Eq. (21) leads to a systematic bias of up to about 2% on the monopole when using a large binning in the input p(r 1 ), but converges towards Eq. (20) result when a small binning is adopted, as expected.
In the case of BOSS, we find similar trends but with an higher accuracy, as shown in Fig. 10. We find at most a difference of 5 × 10 −4 on the monopole between the analytical solution, either Eqs. (20) or (21) with a fine p(r 1 ), and direct pair counting. Regarding the quadrupole and hexadecapole, the relative difference is about 1%. Here, the approximation in Eq. (17) is more appropriate since the data sample is larger. This explains the improved accuracy that is reached. We emphasise that the variance in Figs. 9 and 10 only comes from the random samples since a single data catalogue is used. Therefore, increasing the number of random points would reduce this variance. Overall, because of the approximation in Eq. (17), our analytical DR predictions remain biased, exceeding the typical variance introduced by random sampling in direct pair counting.

Conclusion
In this paper, we present general analytical expressions for the random-random and data-random pair counts in the case of a realistic survey geometry. The main results are given in Eqs. (12) (or (14) for the multipole moments) for RR and in Eq. (21) for in green is obtained with GSL using ε = 10 −3 .
DR. These expressions can be solved numerically in an efficient way. This method, which does not rely on generating random mocks, only takes as input the comoving radial distance distribution in an assumed cosmology and the angular selection function two-point correlation function, which only needs to be estimated a single time for a given survey. Once those quantities are provided, the full computation takes about a few minutes to obtain anisotropic pair counts RR(s, µ) and a few seconds for its multipole moments, using a single CPU and standard libraries for three-dimensional integration.
We tested this method in the context of the BOSS and VIPERS survey geometries and found excellent agreements with expected RR pair counts. The predicted counts exhibit a high accuracy for the cases investigated in this work, equivalent to that we would obtain by performing pair counting in random samples of about 1400−1600 more random points than data in those surveys for the monopole. The main advantage is that the method is fast and does not rely on any spatial sampling, while usually we need to generate a random catalogue with at least 50 times the number of objects in the data. We believe that this can be of some use for future surveys with large data samples and very expensive RR pair count calculations.
The DR pair counts can also be calculated analytically based on certain approximations. We found that the results are slightly biased with respect to the expected counts. For VIPERS and BOSS, we found a bias with respect to direct pair counts of 1% and 0.05%, respectively, for the monopole, up to several percents on the quadrupole and hexadecapole. This bias should decrease with the increasing number of data points. When estimating DR for several data samples, we need to compute, for each sample, its angular two-point correlation function with respect to the survey angular selection function.
Overall, the method presented in this paper for efficiently evaluating the survey window two-point function should be very useful when dealing with massive galaxy surveys. The formulae provided are fast in terms of the speed of the evaluation. With further efficient parallelisation (e.g., Hahn 2015), we should be able to compute RR and DR in an extremely small amount of time. In that case, we could imagine RR and DR being evaluated in different cosmologies at each step of a cosmological likelihood analysis. This opens up new horizons for the way we analyse galaxy survey data in the future.