# Overview: Correlation functions

We provide methods for computing correlation functions in real and fourier space. So far this mainly include the two point correlation function in real space and $N$-point correlation functions in Fourier space. What is missing here and would be great to have is to implement the survey versions of all of these methods (which is straight forward, but just needs to be done).

# Paircounting

Paircounting is not always the easiest thing to parallelize effectively unless you have shared memory. This stems from the fact that we in principle need to correlate each pair of particles which in general resides on different tasks which would lead to tons of communication. Our paircounting algorithm is parallelized with MPI and OpenMP, but with the limitation that using MPI all tasks needs to have all the same particles to prevent having to do this communication! The way we do it is to make a grid, bin the particles to the grid and then correlate particles inside these boxes with neighboring boxes. This means we only consider the pairs we need to consider and not pairs that are separated by distances larger than the maximum radius we want. With MPI each task is responsible for correlating each part of the grid. For our general pair counting methods the user have to supply the binning function - the function that gets called for each pair and adds up the result - and can choose to bin up the data as they see fit so the algorithm itself don't have to be modified if one wants to bin up a quantity other than just how many pair we have. The fiducial thing implemented simply counts the number of pairs. The interface can be improved, see Paircount.h for how to use it.

For standard paircount, i.e. when we just want the number of pairs, we have the following methods:

// Do paircount of a set of particles
// Computes how many pairs we have
template<class T>
AutoPairCountData AutoPairCount(
std::vector<T> & particles,
int nbins,
double rmax,
bool periodic_box,
bool verbose);

// Do cross paircount of two sets of particles
// Computes how many pairs we have
template<class T, class U>
CrossPairCountData CrossPairCount(
std::vector<T> & particles1,
std::vector<U> & particles2,
int nbins,
double rmax,
bool periodic_box,
bool verbose);


The results are provided as follows:

// The result struct for auto pair counting
struct AutoPairCountData{
std::vector<double> r;
std::vector<double> r_edge;
std::vector<double> paircount;
double sum_weights;
double sum_weights_squared;
double norm;
};

// The result struct for cross pair counting
struct CrossPairCountData{
std::vector<double> r;
std::vector<double> r_edge;
std::vector<double> paircount;
double sum_weights;
double sum_weights_squared;
double sum2_weights;
double sum2_weights_squared;
double norm;
};


If you want to compute some more things for the pairs we provide the following methods:

// The general algorithm
// This is the method that does the hard work
// binning_func is the binning function telling us what to do
// with every pair
template<class T>
void AutoPairCountGridMethod(
FML::PARTICLE::ParticlesInBoxes<T> & grid,
std::function<void(int, double*, T&, T&)> &binning_func,
double rmax,
bool periodic,
bool verbose);

// The general algorithm
// This is the method that does the hard work
// binning_func is the binning function telling us what to do
// with every pair
template<class T, class U>
void CrossPairCountGridMethod(
FML::PARTICLE::ParticlesInBoxes<T> & grid,
FML::PARTICLE::ParticlesInBoxes<U> & grid2,
std::function<void(int, double*, T&, U&)> & binning_func,
double rmax,
bool periodic,
bool verbose);



For survey data, with random particles, we provide some common estimators for computing the correlation functions from the paircounts. For auto paircounts the availiable methods are LZ (Landy & Szalay 1993), HAM (Hamilton 1993), HEW (Hewett 1982), DP (Davis & Peebles 1983) and PH (Peebles & Hauser 1974) and for cross paircounts LZ (Landy & Szalay 1993), HAM (Hamilton 1993), DP_NO_R1 and DP_NO_R2 (Davis & Peebles 1983). See Blake et al. 2006 for more info.

// Some estimators for correlation functions in surveys
std::vector AutoCorrelationEstimator(
double *DD,
double *DR,
double *RR,
int nbins,
std::string estimator);

// Some estimators for cross correlation functions in surveys
std::vector CrossCorrelationEstimator(
double *D1D2,
double *D1R2,
double *R1D2,
double *R1R2,
int nbins,
std::string estimator);


# Power-spectrum

The power-spectrum is just binning up the norm of the Fourier modes so this is easy to do. For parallelization we just need to merge the binnings on each task in the end. We also provide a method using interlaced grids to reduce the effect of aliasing. The usual brute-force way of dealing with this is just using a very large grid, but with interlacing we can get away with a much smaller grid and thereby saving computational time and memory. For computing power-spectra we provide the following methods:

// Bin up the power-spectrum of a given fourier grid
template<int N>
void bin_up_power_spectrum(
FFTWGrid<N> & fourier_grid,
PowerSpectrumBinning<N> & pofk);

// Assign particles to grid using density_assignment_method = NGP,CIC,TSC,PCS,...
// Fourier transform, decolvolve the window function for the assignement above,
// With interlaing we assign particles to an interlaced grid (displaced by dx/2 in all directions)
// Fourier transform both and add the together to cancel the leading aliasing contributions
// bin up power-spectrum and subtract shot-noise 1/NumPartTotal
// Set pofk.subtract_shotnoise = false if this is not wanted
template<int N, class T>
void compute_power_spectrum(
int Ngrid,
T *part,
size_t NumPart,
size_t NumPartTotal,
PowerSpectrumBinning<N> & pofk,
std::string density_assignment_method,
bool interlacing);

// Brute force. Add particles to the grid using direct summation
// This gives alias free P(k), but scales as O(Npart)*O(Nmesh^N)
template<int N, class T>
void compute_power_spectrum_direct_summation(
int Ngrid,
T *part,
size_t NumPart,
PowerSpectrumBinning<N> & pofk);


PowerSpectrumBinning<N> pofk;

// ... compute it

// Put in physical units (scale k and P(k))
pofk.scale(boxsize);

// Fetch the arrays with k and P(k)
auto k  = pofk.k;
auto pk = pofk.pofk;


Direct summation above means computing the overdensity field as $\delta(k) = -1 + \frac{1}{N_{\rm part}}\sum_{i=1}^{N_{\rm part}} e^{-ik\cdot x_i}$. This gives us an alias-free estimate for the power-spectrum, however its insanely slow as it scales as $N_{\rm grid}^{\rm NDIM} \cdot N_{\rm part}$. This is only really useful for testing.

# Power-spectrum multipoles

Compute power-spectrum multipoles for when we have a fixed line of sight direction:

// Compute power-spectrum multipoles P0,P1,...,Pn-1 for the case
// where we have a fixed line_of_sight_direction (typical coordinate axes like (0,0,1))
// Pell contains P0,P1,P2,...Pell_max where ell_max = n-1 is the size of Pell
template<int N>
void compute_power_spectrum_multipoles(
FFTWGrid<N> & fourier_grid,
std::vector<PowerSpectrumBinning<N>> & Pell,
std::vector<double> line_of_sight_direction);

// Simple method to estimate multipoles from simulation data
// Take particles in realspace and use their velocity to put them into
// redshift space. Fourier transform and compute multipoles from this like in the method above.
// We do this for all coordinate axes and return the mean P0,P1,... we get from this
// velocity_to_displacement is factor to convert your velocity to a coordinate shift in [0,1)
// e.g. c/(a H Box) with H ~ 100 h km/s/Mpc and Box boxsize in Mpc/h if velocities are peculiar
// Subtracts shot-noise for the monopole. Set Pell.subtract_shotnoise = false if this is not wanted
template<int N, class T>
void compute_power_spectrum_multipoles(
int Ngrid,
MPIParticles<T> & part,
double velocity_to_displacement,
std::vector<PowerSpectrumBinning<N>> & Pell,
std::string density_assignment_method,
bool interlacing);


# Hankel transforms (FFTLog)

Hankel transforms and particulary the FFTLog algorithm have many great uses and particulart to convert between a 3D radial correlation function in real and fourier space. We provide a slightly modified version of the C++ library by A. Slosar to do just that:

//==========================================================================
// Compute the correlation function xi(r) from a power spectrum P(k), sampled
// at logarithmically spaced points k[j]
//==========================================================================
std::pair<DVector,DVector> ComputeCorrelationFunction(const DVector & k, const DVector & pk);

//==========================================================================
// Compute the power spectrum P(k) from a correlation function xi(r), sampled
// at logarithmically spaced points r[i]
//==========================================================================
std::pair<DVector,DVector> ComputePowerSpectrum(const DVector & r, const DVector & xi);

//==========================================================================
// Compute the discrete Hankel transform of the function a(r). See the FFTLog
// documentation for a description of exactly what this function computes.
//==========================================================================
std::pair<DVector, CVector> DiscreteHankelTransform(
const DVector & r,
const CVector & a,
double mu,
double q = 0,
double kcrc = 1,
bool noring = true,
CVector* u =  nullptr);


The first method computes $$\xi(r) = \int_0^\infty \frac{k^3 P(k)}{2\pi^2} \frac{\sin(kr)}{kr} \frac{dk}{k}$$ and returns $(r,\xi)$. The second one computes $$P(k) = \int_0^\infty \xi(r) \frac{\sin(kr)}{kr} 4\pi r^2 dr$$ and returns $(k,P(k))$. The third one computes $$a(k) = \int a(r) (kr)^qJ_\mu(kr) kdr$$ and returns $(k,a)$ where $J_\mu$ is a Bessel function. Read this for how to properly use this algorithm. You need to be careful to extend the input data in a sensible way to prevent ringing. You should always check that you have a good setup by applying it to test data where you know the result.

# Bispectrum

Computing the bispectrum $\left<\delta(k_1)\delta(k_2)\delta(k_3\right> \equiv \delta(k_1+k_2+k_3)B(k_1,k_2,k_3)$ requires us to bin up over closed triangled in Fourier space. Just as with paircounts in realspace doing this the straight forward way is not a parallelization friendly operation as fourier modes resides on different tasks. However there is a way to do the binning with just fourier transforms. With $N_b$ $k$-bins we can do the whole binning with $N_b$ fourier transforms and $N_b^3/8$ real-space integrals (and we are forced to store $N_b$ grids so it can get memory expensive for large $N_b$). What we compute is $$F_i = \int d^Dq_i \delta_{\rm shell~about~k_i}e^{iqx}$$ $$N_i = \int d^Dq_i 1_{\rm shell~about~k_i}e^{iqx}$$ where $X_{\rm shell~about~k_i}$ equals $X$ for $-\Delta k/2 < q - k_i< \Delta k/2$ and is zero otherwise (i.e. its $X$ in a shell about $k_i$). We can then estimate the bispectrum as $$B(k_1,k_2,k_3) = \frac{\int d^3 x F_{123}}{\int d^3 x N_{123}}$$ where $F_{123} = F_1F_2F_3$ and $N_{123} = N_1N_2N_2$. What this in effect is doing is binning up $\delta\delta\delta$ over closed triangles where $N_{123}$ counts the number of triangles in each bin (which we need to normalize the sum). One can estimate $N_{123}$ analytically and use this, but the fiducial option is to compute it numerically (which means another $N_b$ fourier transforms). If one were to compute many bispectra with the same settings then one would only need to compute $N_{123}$ once so its possible with just a few modifications to cut down the runtime by a factor of $2$. The same method easily generalizes to higher order correlators (e.g. for the tripspectrum we get $F_{1234}$ and $N_{1234}$ in the integrals above). The dimension $D$ has no big role in this algorithm (other than the fact that the complexity goes as $\propto D^3$) so in principle it works in any dimension $D\geq 2$.

It is straight forward to generalize this to compute bispectrum multipoles and do this for a survey setup (though I have not bothered doing this yet). See Scoccimarro 2015 for more info.

// Computes the bispectrum B(k1,k2,k3) from particles
// This deconvolves window function and does interlacing to reduce the alias contribution
template<int N, class T>
void compute_bispectrum(
int Ngrid,
T *part,
size_t NumPart,
size_t NumPartTotal,
BispectrumBinning<N> &bofk,
std::string density_assignment_method,
bool interlacing);

// Computes the bispectrum B(k1,k2,k3) from a fourier grid
template<int N>
void compute_bispectrum(
FFTWGrid<N> & density_k,
BispectrumBinning<N> & bofk);



# General polyspectra

The general version of the method described above, it returns the matrix $\left<\delta(k_1)\delta(k_2)\cdots\right> = P(k1,k2,k3,...)$. Below ORDER is how many fields to correlate. 2 is power-spectrum, 3 is bispectrum, 4 is trispectrum etc. The output has a memory footprint of $N_{\rm bins}^{\rm ORDER}$ elements. We have to do $N_b^{\rm ORDER}$ integrals so needless to say with a large number of bins this is quite expensive for large orders.

// Computes the polyspectrum P(k1,k2,k3,...) from particles
// This deconvolves window function and does interlacing to reduce the alias contribution if interlacing = true
// For the monospectrum and bispectrum we subtract shot-noise. Set polyofk.subtract_shotnoise = false if this is not wanted
template<int N, class T, int ORDER>
void compute_polyspectrum(
int Ngrid,
T *part,
size_t NumPart,
size_t NumPartTotal,
PolyspectrumBinning<N> & polyofk,
std::string density_assignment_method,
bool interlacing);

// Computes the polyspectrum P(k1,k2,k3,...) from a fourier grid
template<int N, int ORDER>
void compute_polyspectrum(
FFTWGrid<N> & density_k,
PolyspectrumBinning<N> & polyofk);

// We also have the wrapper methods (with same signature as the ones above)
// For computing P(k1,k2) (overkill, but useful for testing),
// B(k1,k2,k3) and T(k1,k2,k3,k4)
template<int N, class T, int ORDER>
void compute_monospectrum( ... )  // P(k1,k2) (diagonal is just the power-spectrum P(k))
void compute_bispectrum( ... )    // B(k1,k2,k3)
void compute_trispectrum( ... )   // T(k1,k2,k3,k4)