FML
Namespaces | Data Structures | Typedefs | Functions
FML::CORRELATIONFUNCTIONS Namespace Reference

This namespace deals with computing correlations functions (N-point functions) in real and fourier space. More...

Namespaces

namespace  PAIRCOUNTS
 This namespace deals with paircounting and contains general auto and cross paircount methods in 1D, 2D or 3D for simulation boxes and survey data (with randoms).
 

Data Structures

class  PolyspectrumBinning
 Class for storing the results from general polyspectrum. More...
 
class  PowerSpectrumBinning
 Class for holding the results after binning up a power-spectrum. More...
 

Typedefs

template<int N>
using FFTWGrid = FML::GRID::FFTWGrid< N >
 
template<int N>
using MonospectrumBinning = PolyspectrumBinning< N, 2 >
 
template<int N>
using BispectrumBinning = PolyspectrumBinning< N, 3 >
 
template<int N>
using TrispectrumBinning = PolyspectrumBinning< N, 4 >
 

Functions

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)
 Assign particles to grid using density_assignment_method = NGP,CIC,TSC,PCS,... Fourier transform, decolvolve the window function for the assignement above. More...
 
template<int N, class T >
void compute_power_spectrum_direct_summation (int Ngrid, const T *part, size_t NumPart, PowerSpectrumBinning< N > &pofk)
 Brute force (but aliasing free) computation of the power spectrum. More...
 
template<int N>
void bin_up_power_spectrum (const FFTWGrid< N > &fourier_grid, PowerSpectrumBinning< N > &pofk)
 Compute the power-spectrum of a fourier grid. More...
 
template<int N>
void bin_up_deconvolved_power_spectrum (const FFTWGrid< N > &fourier_grid, PowerSpectrumBinning< N > &pofk, std::string density_assignment_method)
 Compute the power-spectrum of a fourier grid. More...
 
template<int N>
void bin_up_cross_power_spectrum (const FFTWGrid< N > &fourier_grid_1, const FFTWGrid< N > &fourier_grid_2, PowerSpectrumBinning< N > &pofk)
 Compute the cross power-spectrum of two fourier grids. More...
 
template<int N>
void compute_power_spectrum_multipoles_fourier (const FFTWGrid< N > &fourier_grid, std::vector< PowerSpectrumBinning< N > > &Pell, std::vector< double > line_of_sight_direction)
 Compute power-spectrum multipoles from a Fourier grid where we have a fixed line_of_sight_direction (typical coordinate axes like \( (0,0,1) \)). More...
 
template<int N, class T >
void compute_power_spectrum_multipoles (int Ngrid, FML::PARTICLE::MPIParticles< T > &part, double velocity_to_displacement, std::vector< PowerSpectrumBinning< N > > &Pell, std::string density_assignment_method, bool interlacing)
 A simple power-spectrum estimator for multipoles in simulations. More...
 
template<int N, class T , int ORDER>
void compute_polyspectrum (int Ngrid, T *part, size_t NumPart, size_t NumPartTotal, PolyspectrumBinning< N, ORDER > &polyofk, std::string density_assignment_method, bool interlacing)
 Computes the polyspectrum \( P(k_1,k_2,\ldots,k_{\rm ORDER}) = \left<\delta(k_1)\cdots\delta(k_{\rm ORDER})\right> \) from particles. More...
 
template<int N, int ORDER>
void compute_polyspectrum (const FFTWGrid< N > &fourier_grid, PolyspectrumBinning< N, ORDER > &polyofk)
 Computes the polyspectrum \( P(k_1,k_2,\ldots,k_{\rm ORDER}) = \left<\delta(k_1)\cdots\delta(k_{\rm ORDER})\right> \) from a fourier grid. More...
 
template<int N>
void compute_monospectrum (const FFTWGrid< N > &fourier_grid, MonospectrumBinning< N > &pofk)
 Computes the monospectrum \( P(k_1,k_2) = \left<\delta(k_1)\delta(k_2)\right> \) from a fourier grid. More...
 
template<int N, class T >
void compute_monospectrum (int Ngrid, T *part, size_t NumPart, size_t NumPartTotal, MonospectrumBinning< N > &pofk, std::string density_assignment_method, bool interlacing)
 Computes the monospectrum \( P(k_1,k_2) = \left<\delta(k_1)\delta(k_2)\right> \) from particles. More...
 
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(k_1,k_2,k_3) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\right> \) from particles. More...
 
template<int N>
void compute_bispectrum (const FFTWGrid< N > &fourier_grid, BispectrumBinning< N > &bofk)
 Computes the bispectrum \( B(k_1,k_2,k_3) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\right> \) from a fourier grid. More...
 
template<int N>
void compute_trispectrum (const FFTWGrid< N > &fourier_grid, TrispectrumBinning< N > &tofk)
 Computes the trispectrum \( T(k_1,k_2,k_3,k_4) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\delta(k_4)\right> \) from a fourier grid. More...
 
template<int N, class T >
void compute_trispectrum (int Ngrid, T *part, size_t NumPart, size_t NumPartTotal, TrispectrumBinning< N > &tofk, std::string density_assignment_method, bool interlacing)
 Computes the truspectrum \( T(k_1,k_2,k_3,k_4) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\delta(k_4)\right> \) from particles. More...
 
template<int N>
void bin_up_cross_power_spectrum (FFTWGrid< N > &fourier_grid_1, FFTWGrid< N > &fourier_grid_2, PowerSpectrumBinning< N > &pofk)
 
template<int N, class T >
void compute_power_spectrum_direct_summation (int Ngrid, T *part, size_t NumPart, PowerSpectrumBinning< N > &pofk)
 
template<int N, int ORDER>
void compute_polyspectrum_bincount (int Nmesh, PolyspectrumBinning< N, ORDER > &polyofk)
 This method is used by compute_polyspectrum. More...
 
template<int N>
void compute_monospectrum (const FFTWGrid< N > &fourier_grid, PolyspectrumBinning< N, 2 > &pofk)
 
template<int N>
void compute_bispectrum (const FFTWGrid< N > &fourier_grid, PolyspectrumBinning< N, 3 > &bofk)
 
template<int N>
void compute_trispectrum (const FFTWGrid< N > &fourier_grid, PolyspectrumBinning< N, 4 > &tofk)
 

Detailed Description

This namespace deals with computing correlations functions (N-point functions) in real and fourier space.

Typedef Documentation

◆ BispectrumBinning

Definition at line 61 of file ComputePowerSpectrum.h.

◆ FFTWGrid

template<int N>
using FML::CORRELATIONFUNCTIONS::FFTWGrid = typedef FML::GRID::FFTWGrid<N>

Definition at line 39 of file ComputePowerSpectrum.h.

◆ MonospectrumBinning

Definition at line 58 of file ComputePowerSpectrum.h.

◆ TrispectrumBinning

Definition at line 64 of file ComputePowerSpectrum.h.

Function Documentation

◆ bin_up_cross_power_spectrum() [1/2]

template<int N>
void FML::CORRELATIONFUNCTIONS::bin_up_cross_power_spectrum ( const FFTWGrid< N > &  fourier_grid_1,
const FFTWGrid< N > &  fourier_grid_2,
PowerSpectrumBinning< N > &  pofk 
)

Compute the cross power-spectrum of two fourier grids.

The result has no scales. Get scales by calling pofk.scale(boxsize) which does k \( k \to k/B \) and \( P(k) \to B^N P(k) \) where \( B\) is the boxsize once spectrum has been computed. The method assumes the two grids are fourier transforms of real grids (i.e. \( f(-k) = f^*(k) \)) and we only bin up the real part of \( f_1(k)f_2^*(k) \). The imaginary part is also binned up, but not returned.

Template Parameters
NDimension of the grid
Parameters
[in]fourier_grid_1Grid in fourier space
[in]fourier_grid_2Grid in fourier space
[out]pofkBinned cross power-spectrum

◆ bin_up_cross_power_spectrum() [2/2]

template<int N>
void FML::CORRELATIONFUNCTIONS::bin_up_cross_power_spectrum ( FFTWGrid< N > &  fourier_grid_1,
FFTWGrid< N > &  fourier_grid_2,
PowerSpectrumBinning< N > &  pofk 
)

Definition at line 630 of file ComputePowerSpectrum.h.

632 {
633
634 assert_mpi(fourier_grid_1.get_nmesh() > 0, "[bin_up_cross_power_spectrum] grid must have Nmesh > 0\n");
635 assert_mpi(fourier_grid_1.get_nmesh() == fourier_grid_2.get_nmesh(),
636 "[bin_up_cross_power_spectrum] Grids must have the same gridsize\n");
637 assert_mpi(pofk.n > 0 && pofk.kmax > pofk.kmin && pofk.kmin >= 0.0,
638 "[bin_up_cross_power_spectrum] Binning has inconsistent parameters\n");
639
640 const auto Nmesh = fourier_grid_1.get_nmesh();
641 const auto Local_nx = fourier_grid_1.get_local_nx();
642 const auto Local_x_start = fourier_grid_1.get_local_x_start();
643
644 // Initialize binning just in case
645 pofk.reset();
646
647 // Take a copy and bin up the imaginary part also
648 PowerSpectrumBinning<N> pofk_imag = pofk;
649 pofk_imag.reset();
650
651 // Bin up P(k)
652#ifdef USE_OMP
653#pragma omp parallel for
654#endif
655 for (int islice = 0; islice < Local_nx; islice++) {
656 [[maybe_unused]] double kmag;
657 [[maybe_unused]] std::array<double, N> kvec;
658 for (auto && fourier_index : fourier_grid_1.get_fourier_range(islice, islice + 1)) {
659 if (Local_x_start == 0 and fourier_index == 0)
660 continue; // DC mode( k=0)
661
662 // Special treatment of k = 0 plane (Safer way: fetch coord)
663 // auto coord = fourier_grid.get_fourier_coord_from_index(fourier_index);
664 // int last_coord = coord[N-1];
665 auto last_coord = fourier_index % (Nmesh / 2 + 1);
666 double weight = last_coord > 0 && last_coord < Nmesh / 2 ? 2.0 : 1.0;
667
668 auto delta_1 = fourier_grid_1.get_fourier_from_index(fourier_index);
669 auto delta_2 = fourier_grid_2.get_fourier_from_index(fourier_index);
670 auto delta12_real = delta_1.real() * delta_2.real() + delta_1.imag() * delta_2.imag();
671 auto delta12_imag = -delta_1.real() * delta_2.imag() + delta_1.imag() * delta_2.real();
672
673 // Add norm to bin
674 fourier_grid_1.get_fourier_wavevector_and_norm_by_index(fourier_index, kvec, kmag);
675 pofk.add_to_bin(kmag, delta12_real, weight);
676 pofk_imag.add_to_bin(kmag, delta12_imag, weight);
677 }
678 }
679
680 // Normalize to get P(k) (this communicates over tasks)
681 pofk.normalize();
682 pofk_imag.normalize();
683
684 // NB: we currently don't return the imaginary part. For real fields this should be zero
685 // so we just check this and give a warning if its large
686 const bool show_warning = false;
687 if (show_warning) {
688 for (int i = 0; i < pofk.n; i++) {
689 if (std::abs(pofk.pofk[i]) < 1e3 * std::abs(pofk_imag.pofk[i]) and
690 std::abs(pofk_imag.pofk[i]) > 1e-15) {
691 std::cout
692 << "Warning: the imaginary part of the cross spectrum is > 0.1%% times the real part [ k: "
693 << pofk.kbin[i] << " Real(d1d2): " << pofk.pofk[i] << " Imag(d1d2): " << pofk_imag.pofk[i]
694 << "]\n";
695 }
696 }
697 }
698 }
const int Nmesh
Definition: Main.cpp:48
Class for holding the results after binning up a power-spectrum.
void reset()
Reset everything. Call before starting to bim.
void add_to_bin(double kvalue, double power, double weight=1.0)
Add a new point to a bin.
void normalize()
Normalize (i.e. find mean in each bin) Do summation over MPI tasks.
std::vector< double > pofk
The power-spectrum in the bin.
int get_nmesh() const
Number of grid-cells per dimension.
Definition: FFTWGrid.h:1054
ptrdiff_t get_local_x_start() const
The x-slice the first local slice has in the global grid.
Definition: FFTWGrid.h:1084
ptrdiff_t get_local_nx() const
Number of local x-slices in the real grid.
Definition: FFTWGrid.h:1079
ComplexType get_fourier_from_index(const IndexIntType index) const
Fetch value in fourier grid by (local) index.
Definition: FFTWGrid.h:1031
FourierRange get_fourier_range(int islice_begin=0, int islice_end=0) const
Range iterator for going through all active cells in the main fourier grid by index [ e....
Definition: FFTWGrid.h:419
void get_fourier_wavevector_and_norm_by_index(const IndexIntType ind, std::array< double, N > &kvec, double &kmag) const
From index in the grid get the k-vector and the magnitude of it.
Definition: FFTWGrid.h:1103

◆ bin_up_deconvolved_power_spectrum()

template<int N>
void FML::CORRELATIONFUNCTIONS::bin_up_deconvolved_power_spectrum ( const FFTWGrid< N > &  fourier_grid,
PowerSpectrumBinning< N > &  pofk,
std::string  density_assignment_method 
)

Compute the power-spectrum of a fourier grid.

The result has no scales. Get scales by calling pofk.scale(boxsize) which does \( k \to k/B \) and \( P(k) \to B^N P(k) \) where \( B\) is the boxsize once spectrum has been computed. The method assumes the two grids are fourier transforms of real grids (i.e. \( f(-k) = f^*(k) \)). This method divides the power-spectrum by the square of the window function corresponding to the density assignement method for every mode we bin up.

Template Parameters
NDimension of the grid
Parameters
[in]fourier_gridGrid in fourier space
[out]pofkBinned power-spectrum
[in]density_assignment_methodThe density assignment method we use to create the grid

Definition at line 534 of file ComputePowerSpectrum.h.

536 {
537
538 assert_mpi(fourier_grid.get_nmesh() > 0, "[bin_up_deconvolved_power_spectrum] grid must have Nmesh > 0\n");
539 assert_mpi(pofk.n > 0 && pofk.kmax > pofk.kmin && pofk.kmin >= 0.0,
540 "[bin_up_deconvolved_power_spectrum] Binning has inconsistent parameters\n");
541
542 const auto Nmesh = fourier_grid.get_nmesh();
543 const auto Local_nx = fourier_grid.get_local_nx();
544 const auto Local_x_start = fourier_grid.get_local_x_start();
545 const auto window_function = FML::INTERPOLATION::get_window_function<N>(density_assignment_method, Nmesh);
546
547 // Initialize binning just in case
548 pofk.reset();
549
550 // Bin up P(k)
551#ifdef USE_OMP
552#pragma omp parallel for
553#endif
554 for (int islice = 0; islice < Local_nx; islice++) {
555 [[maybe_unused]] double kmag;
556 [[maybe_unused]] std::array<double, N> kvec;
557 for (auto && fourier_index : fourier_grid.get_fourier_range(islice, islice + 1)) {
558 if (Local_x_start == 0 and fourier_index == 0)
559 continue; // DC mode( k=0)
560
561 // Special treatment of k = 0 plane (Safer way: fetch coord)
562 // auto coord = fourier_grid.get_fourier_coord_from_index(fourier_index);
563 // int last_coord = coord[N-1];
564 auto last_coord = fourier_index % (Nmesh / 2 + 1);
565 double weight = last_coord > 0 && last_coord < Nmesh / 2 ? 2.0 : 1.0;
566
567 fourier_grid.get_fourier_wavevector_and_norm_by_index(fourier_index, kvec, kmag);
568
569 auto delta = fourier_grid.get_fourier_from_index(fourier_index);
570 auto delta_norm = std::norm(delta);
571 auto window = window_function(kvec);
572 delta_norm /= (window * window);
573
574 // Add norm to bin
575 pofk.add_to_bin(kmag, delta_norm, weight);
576 }
577 }
578
579 // Normalize to get P(k) (this communicates over tasks)
580 pofk.normalize();
581 }

◆ bin_up_power_spectrum()

template<int N>
void FML::CORRELATIONFUNCTIONS::bin_up_power_spectrum ( const FFTWGrid< N > &  fourier_grid,
PowerSpectrumBinning< N > &  pofk 
)

Compute the power-spectrum of a fourier grid.

The result has no scales. Get scales by calling pofk.scale(boxsize) which does \( k \to k/B \) and \( P(k) \to B^N P(k) \) where \( B\) is the boxsize once spectrum has been computed. The method assumes the two grids are fourier transforms of real grids (i.e. \( f(-k) = f^*(k) \)).

Template Parameters
NDimension of the grid
Parameters
[in]fourier_gridGrid in fourier space
[out]pofkBinned power-spectrum

Definition at line 585 of file ComputePowerSpectrum.h.

585 {
586
587 assert_mpi(fourier_grid.get_nmesh() > 0, "[bin_up_power_spectrum] grid must have Nmesh > 0\n");
588 assert_mpi(pofk.n > 0 && pofk.kmax > pofk.kmin && pofk.kmin >= 0.0,
589 "[bin_up_power_spectrum] Binning has inconsistent parameters\n");
590
591 const auto Nmesh = fourier_grid.get_nmesh();
592 const auto Local_nx = fourier_grid.get_local_nx();
593 const auto Local_x_start = fourier_grid.get_local_x_start();
594
595 // Initialize binning just in case
596 pofk.reset();
597
598 // Bin up P(k)
599#ifdef USE_OMP
600#pragma omp parallel for
601#endif
602 for (int islice = 0; islice < Local_nx; islice++) {
603 [[maybe_unused]] double kmag;
604 [[maybe_unused]] std::array<double, N> kvec;
605 for (auto && fourier_index : fourier_grid.get_fourier_range(islice, islice + 1)) {
606 if (Local_x_start == 0 and fourier_index == 0)
607 continue; // DC mode( k=0)
608
609 // Special treatment of k = 0 plane (Safer way: fetch coord)
610 // auto coord = fourier_grid.get_fourier_coord_from_index(fourier_index);
611 // int last_coord = coord[N-1];
612 auto last_coord = fourier_index % (Nmesh / 2 + 1);
613 double weight = last_coord > 0 && last_coord < Nmesh / 2 ? 2.0 : 1.0;
614
615 auto delta = fourier_grid.get_fourier_from_index(fourier_index);
616 auto delta_norm = std::norm(delta);
617
618 // Add norm to bin
619 fourier_grid.get_fourier_wavevector_and_norm_by_index(fourier_index, kvec, kmag);
620 pofk.add_to_bin(kmag, delta_norm, weight);
621 }
622 }
623
624 // Normalize to get P(k) (this communicates over tasks)
625 pofk.normalize();
626 }

◆ compute_bispectrum() [1/3]

template<int N>
void FML::CORRELATIONFUNCTIONS::compute_bispectrum ( const FFTWGrid< N > &  fourier_grid,
BispectrumBinning< N > &  bofk 
)

Computes the bispectrum \( B(k_1,k_2,k_3) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\right> \) from a fourier grid.

This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{3} \) integrals. This is just an alias for compute_polyspectrum<N, 3>

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]fourier_gridA fourier grid.
[out]bofkThe binned bispectrum. We required it to be initialized with the number of bins, kmin and kmax.

◆ compute_bispectrum() [2/3]

template<int N>
void FML::CORRELATIONFUNCTIONS::compute_bispectrum ( const FFTWGrid< N > &  fourier_grid,
PolyspectrumBinning< N, 3 > &  bofk 
)

Definition at line 1270 of file ComputePowerSpectrum.h.

1270 {
1271 compute_polyspectrum<N, 3>(fourier_grid, bofk);
1272 }

◆ compute_bispectrum() [3/3]

template<int N, class T >
void FML::CORRELATIONFUNCTIONS::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(k_1,k_2,k_3) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\right> \) from particles.

Note that with interlacing we change the particle positions, but when returning they should be in the same state as when we started. This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{3} \) integrals. This is just an alias for compute_polyspectrum<N, 3>. Shot-noise is subtracted.

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]NgridSize of the grid to use.
[in]partPointer to the first particle.
[in]NumPartNumber of particles on the local task.
[in]NumPartTotalTotal number of particles on all tasks.
[out]bofkThe binned bispectrum. We required it to be initialized with the number of bins, kmin and kmax.
[in]density_assignment_methodThe density assignment method (NGP, CIC, TSC, PCS or PQS)
[in]interlacingUse interlaced grids when computing density field for alias reduction

Definition at line 1043 of file ComputePowerSpectrum.h.

1049 {
1050
1051 // This will subtract shot-noise if bofk.subtract_shotnoise = true
1052 compute_polyspectrum<N, T, 3>(
1053 Ngrid, part, NumPart, NumPartTotal, bofk, density_assignment_method, interlacing);
1054 }

◆ compute_monospectrum() [1/3]

template<int N>
void FML::CORRELATIONFUNCTIONS::compute_monospectrum ( const FFTWGrid< N > &  fourier_grid,
MonospectrumBinning< N > &  pofk 
)

Computes the monospectrum \( P(k_1,k_2) = \left<\delta(k_1)\delta(k_2)\right> \) from a fourier grid.

This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{2} \) integrals. This is just an alias for compute_polyspectrum<N, 2>

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]fourier_gridA fourier grid.
[out]pofkThe binned monospectrum. We required it to be initialized with the number of bins, kmin and kmax.

◆ compute_monospectrum() [2/3]

template<int N>
void FML::CORRELATIONFUNCTIONS::compute_monospectrum ( const FFTWGrid< N > &  fourier_grid,
PolyspectrumBinning< N, 2 > &  pofk 
)

Definition at line 1265 of file ComputePowerSpectrum.h.

1265 {
1266 compute_polyspectrum<N, 2>(fourier_grid, pofk);
1267 }

◆ compute_monospectrum() [3/3]

template<int N, class T >
void FML::CORRELATIONFUNCTIONS::compute_monospectrum ( int  Ngrid,
T *  part,
size_t  NumPart,
size_t  NumPartTotal,
MonospectrumBinning< N > &  pofk,
std::string  density_assignment_method,
bool  interlacing 
)

Computes the monospectrum \( P(k_1,k_2) = \left<\delta(k_1)\delta(k_2)\right> \) from particles.

Note that with interlacing we change the particle positions, but when returning they should be in the same state as when we started. This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{2} \) integrals. This is just an alias for compute_polyspectrum<N, 2>. Shot-noise is subtracted

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]NgridSize of the grid to use.
[in]partPointer to the first particle.
[in]NumPartNumber of particles on the local task.
[in]NumPartTotalTotal number of particles on all tasks.
[out]pofkThe binned monospectrum. We required it to be initialized with the number of bins, kmin and kmax.
[in]density_assignment_methodThe density assignment method (NGP, CIC, TSC, PCS or PQS)
[in]interlacingUse interlaced grids when computing density field for alias reduction

Definition at line 1029 of file ComputePowerSpectrum.h.

1035 {
1036
1037 // This will subtract shot-noise if pofk.subtract_shotnoise = true
1038 compute_polyspectrum<N, T, 2>(
1039 Ngrid, part, NumPart, NumPartTotal, pofk, density_assignment_method, interlacing);
1040 }

◆ compute_polyspectrum() [1/2]

template<int N, int ORDER>
void FML::CORRELATIONFUNCTIONS::compute_polyspectrum ( const FFTWGrid< N > &  fourier_grid,
PolyspectrumBinning< N, ORDER > &  polyofk 
)

Computes the polyspectrum \( P(k_1,k_2,\ldots,k_{\rm ORDER}) = \left<\delta(k_1)\cdots\delta(k_{\rm ORDER})\right> \) from a fourier grid.

This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{\rm ORDER} \) integrals. If one is to compute many spectra with the same Ngrid and binning then one can precompute N123 in polyofk and set it using polyofk.set_bincount(N123). This speeds up the polyspectrum estimation by a factor of 2 by avoiding half of the fourier transforms.

Template Parameters
NThe dimension of the particles.
ORDERThe order. 2 is the power-spectrum, 3 is the bispectrum, 4 is the trispectrum.
Parameters
[in]fourier_gridGrid in fourier space
[out]polyofkThe binned polyspectrum.

Definition at line 1280 of file ComputePowerSpectrum.h.

1280 {
1281
1282 const auto Nmesh = fourier_grid.get_nmesh();
1283 const auto Local_nx = fourier_grid.get_local_nx();
1284 const auto nbins = polyofk.n;
1285
1286 assert_mpi(nbins > 0, "[compute_polyspectrum] nbins has to be >=0\n");
1287 assert_mpi(Nmesh > 0, "[compute_polyspectrum] grid is not allocated\n");
1288 assert_mpi(polyofk.P123.size() == polyofk.N123.size() and
1289 polyofk.N123.size() == size_t(FML::power(nbins, ORDER)),
1290 "[compute_polyspectrum] Binning is not good\n");
1291 static_assert(ORDER > 1);
1292
1293 // Get where to store the resuts in
1294 // We don't add to thes below so we don't need to clear the arrays
1295 std::vector<double> & P123 = polyofk.P123;
1296 std::vector<double> & N123 = polyofk.N123;
1297 std::vector<double> & pofk_bin = polyofk.pofk;
1298 std::vector<double> & kmean = polyofk.kmean;
1299 const size_t nbins_tot = P123.size();
1300
1301 // Get ranges for which we will compute F_k on
1302 const std::vector<double> & kbin = polyofk.kbin;
1303 const std::vector<double> & klow = polyofk.klow;
1304 const std::vector<double> & khigh = polyofk.khigh;
1305
1306 // Compute the bincount N123 if it does not already exist
1307 if (not polyofk.bincount_is_set)
1308 compute_polyspectrum_bincount<N, ORDER>(Nmesh, polyofk);
1309
1310 // Allocate grids
1311 std::vector<FFTWGrid<N>> F_k(nbins);
1312 for (int i = 0; i < nbins; i++) {
1313 F_k[i] = fourier_grid;
1314 F_k[i].add_memory_label("FFTWGrid::compute_polyspectrum::F_" + std::to_string(i));
1315 }
1316
1317 for (int i = 0; i < nbins; i++) {
1318#ifdef DEBUG_POLYSPECTRUM
1319 if (FML::ThisTask == 0)
1320 std::cout << "Computing polyspectrum<" << ORDER << "> " << i + 1 << " / " << nbins
1321 << " kbin: " << klow[i] / (2.0 * M_PI) << " -> " << khigh[i] / (2.0 * M_PI) << "\n";
1322#endif
1323
1324 FFTWGrid<N> & grid = F_k[i];
1325
1326 // For each bin get klow, khigh
1327 const double kmag2_max = khigh[i] * khigh[i];
1328 const double kmag2_min = klow[i] * klow[i];
1329
1330 // Loop over all cells
1331 double kmean_bin = 0.0;
1332 double nk = 0;
1333 double kmag2;
1334 pofk_bin[i] = 0.0;
1335 std::array<double, N> kvec;
1336 for (auto && fourier_index : grid.get_fourier_range()) {
1337 grid.get_fourier_wavevector_and_norm2_by_index(fourier_index, kvec, kmag2);
1338
1339 // Set to zero outside the bin
1340 if (kmag2 >= kmag2_max or kmag2 < kmag2_min) {
1341 grid.set_fourier_from_index(fourier_index, 0.0);
1342 } else {
1343 // Compute mean k in the bin
1344 kmean_bin += std::sqrt(kmag2);
1345 pofk_bin[i] += std::norm(grid.get_fourier_from_index(fourier_index));
1346 nk += 1.0;
1347 }
1348 }
1349 FML::SumOverTasks(&kmean_bin);
1350 FML::SumOverTasks(&pofk_bin[i]);
1351 FML::SumOverTasks(&nk);
1352
1353 // The mean k in the bin
1354 kmean[i] = (nk == 0) ? kbin[i] : kmean_bin / double(nk);
1355
1356 // Power spectrum in the bin
1357 pofk_bin[i] = (nk == 0) ? 0.0 : pofk_bin[i] / double(nk);
1358
1359#ifdef DEBUG_POLYSPECTRUM
1360 if (FML::ThisTask == 0)
1361 std::cout << "kmean: " << kmean[i] / (2.0 * M_PI) << "\n";
1362#endif
1363
1364 // Transform to real space
1365 grid.fftw_c2r();
1366 }
1367
1368 // We now have F_k and N_k for all bins
1369 for (size_t i = 0; i < nbins_tot; i++) {
1370#ifdef DEBUG_POLYSPECTRUM
1371 if (FML::ThisTask == 0)
1372 if ((i * 10) / nbins_tot != ((i + 1) * 10) / nbins_tot)
1373 std::cout << "Integrating up " << 100 * (i + 1) / nbins_tot << " %\n";
1374 ;
1375#endif
1376
1377 // Current values of ik1,ik2,ik3,...
1378 const auto ik = polyofk.get_coord_from_index(i);
1379
1380 // Symmetry: only do ik1 <= ik2 <= ... and don't need to do configurations that don't satisfy the
1381 // triangle inequality
1382 if (not polyofk.compute_this_configuration(ik)) {
1383 P123[i] = 0.0;
1384 continue;
1385 }
1386
1387 // Compute the sum over triangles by evaluating the integral Int dx^N/(2pi)^N
1388 // delta_k1(x)delta_k2(x)...delta_kORDER(x)
1389 double F123_current = 0.0;
1390#ifdef USE_OMP
1391#pragma omp parallel for reduction(+ : F123_current)
1392#endif
1393 for (int islice = 0; islice < Local_nx; islice++) {
1394 for (auto && real_index : F_k[0].get_real_range(islice, islice + 1)) {
1395 double Fproduct = 1.0;
1396 for (int ii = 0; ii < ORDER; ii++)
1397 Fproduct *= F_k[ik[ii]].get_real_from_index(real_index);
1398 F123_current += Fproduct;
1399 }
1400 }
1401 FML::SumOverTasks(&F123_current);
1402
1403 // Normalize by the integration measure dx^N / (2pi)^N
1404 F123_current *= std::pow(1.0 / double(Nmesh) / (2.0 * M_PI), N);
1405
1406 // Set the result
1407 const double N123_current = N123[i];
1408 P123[i] = (N123_current > 0.0) ? F123_current / N123_current : 0.0;
1409 }
1410
1411 // Set stuff not computed above which follows from symmetry
1412 for (size_t i = 0; i < nbins_tot; i++) {
1413
1414 // Current values of ik1,ik2,ik3
1415 auto ik = polyofk.get_coord_from_index(i);
1416
1417 // If its valid its already computed
1418 if (polyofk.compute_this_configuration(ik))
1419 continue;
1420
1421 // Find a cell given by symmetry that we have computed
1422 // by sorting ik in increasing order
1423 std::sort(ik.begin(), ik.end(), std::less<int>());
1424
1425 // Compute cell index of configuration we have computed
1426 size_t index = polyofk.get_index_from_coord(ik);
1427
1428 // Set value
1429 P123[i] = P123[index];
1430 }
1431 }
std::vector< double > pofk
The power-spectrum.
std::vector< double > N123
The volume (number of generalized triangles) of each bin.
std::vector< double > khigh
khigh is higher bin edge
std::vector< double > P123
The polyspectrum.
std::array< int, ORDER > get_coord_from_index(size_t index)
The coord (k1,k2,...) corresponding to a given index in the P123 and N123 array These methods determi...
std::vector< double > kmean
kmean is the mean of wavenumbers added to the bin
std::vector< double > klow
klow is lower bin edge
size_t get_index_from_coord(const std::array< int, ORDER > &ik)
bool bincount_is_set
If we have precomputed the volume factor (N123) or not and have it availiable for algorithms to use.
std::vector< double > kbin
kbin is center of bin
bool compute_this_configuration(const std::array< int, ORDER > &ik)
Symmetry: we only need to compute ik1 <= ik2 <= ... This function tells algorithms which configuratio...
Class for holding grids and performing real-to-complex and complex-to-real FFTs using FFTW with MPI.
Definition: FFTWGrid.h:88
void add_memory_label(std::string label)
For memory logging: add a label to the grid.
Definition: FFTWGrid.h:318
void grid(bool flag)
constexpr long long int power(int base, int exponent)
Simple integer a^b power-function by squaring.
Definition: Global.h:187
int ThisTask
The MPI task number.
Definition: Global.cpp:33
void SumOverTasks(T *value)
Inplace sum over tasks of a single value.
Definition: Global.h:147

◆ compute_polyspectrum() [2/2]

template<int N, class T , int ORDER>
void FML::CORRELATIONFUNCTIONS::compute_polyspectrum ( int  Ngrid,
T *  part,
size_t  NumPart,
size_t  NumPartTotal,
PolyspectrumBinning< N, ORDER > &  polyofk,
std::string  density_assignment_method,
bool  interlacing 
)

Computes the polyspectrum \( P(k_1,k_2,\ldots,k_{\rm ORDER}) = \left<\delta(k_1)\cdots\delta(k_{\rm ORDER})\right> \) from particles.

Note that with interlacing we change the particle positions, but when returning they should be in the same state as when we started. This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{\rm ORDER} \) integrals. If one is to compute many spectra with the same Ngrid and binning then one can precompute N123 in polyofk and set it using polyofk.set_bincount(N123). This speeds up the polyspectrum estimation by a factor of 2. Shot-noise term is subtracted if ORDER=2 (power) or ORDER=3 (bi). Doing this for higher order requires all lower order correlators and we currently only bin up the polyspectrum and the power-spectrum.

Template Parameters
NThe dimension of the particles.
ORDERThe order. 2 is the power-spectrum, 3 is the bispectrum, 4 is the trispectrum.
Parameters
[in]NgridSize of the grid to use.
[in]partPointer to the first particle.
[in]NumPartNumber of particles on the local task.
[in]NumPartTotalTotal number of particles on all tasks.
[out]polyofkThe binned polyspectrum. We required it to be initialized with the number of bins, kmin and kmax.
[in]density_assignment_methodThe density assignment method (NGP, CIC, TSC, PCS or PQS)
[in]interlacingUse interlaced grids when computing density field for alias reduction

Definition at line 1070 of file ComputePowerSpectrum.h.

1076 {
1077
1078 // Set how many extra slices we need for the density assignment to go smoothly
1079 const auto nleftright = get_extra_slices_needed_for_density_assignment(density_assignment_method);
1080 const int nleft = nleftright.first;
1081 const int nright = nleftright.second + 1;
1082
1083 FFTWGrid<N> density_k(Ngrid, nleft, nright);
1084 density_k.add_memory_label("FFTWGrid::compute_polyspectrum::density_k");
1085
1086 if (interlacing) {
1087
1088 // Bin particles to grid (use interlaced to reduce alias) and deconvolve window
1090 part, NumPart, NumPartTotal, density_k, density_assignment_method);
1091 deconvolve_window_function_fourier<N>(density_k, density_assignment_method);
1092
1093 } else {
1094
1095 // Bin to grid, fourier transform and deconvolve window function
1096 particles_to_grid<N, T>(part, NumPart, NumPartTotal, density_k, density_assignment_method);
1097 density_k.fftw_r2c();
1098 deconvolve_window_function_fourier<N>(density_k, density_assignment_method);
1099 }
1100
1101 // Compute polyspectrum
1102 compute_polyspectrum<N, ORDER>(density_k, polyofk);
1103
1104 // Subtract shotnoise if ORDER = 2 or 3
1105 if constexpr (ORDER == 2) {
1106 // Subtract shot-noise (along the diagonal)
1107 const double n = double(NumPartTotal);
1108 if (polyofk.subtract_shotnoise)
1109 for (int i = 0; i < polyofk.n; i++) {
1110 polyofk.pofk[i] -= 1.0 / n;
1111 std::array<int, 2> ik{i, i};
1112 auto index = polyofk.get_index_from_coord(ik);
1113 polyofk.P123[index] -= 1.0 / n;
1114 }
1115 } else if constexpr (ORDER == 3) {
1116 // Subtract shot-noise for both P(k) and B(k1,k2,k3)
1117 const double n = double(NumPartTotal);
1118 if (polyofk.subtract_shotnoise) {
1119 for (int i = 0; i < polyofk.n; i++) {
1120 polyofk.pofk[i] -= 1.0 / n;
1121 }
1122 for (size_t i = 0; i < polyofk.P123.size(); i++) {
1123 auto ik = polyofk.get_coord_from_index(i);
1124 polyofk.P123[i] -=
1125 1.0 / n * (polyofk.pofk[ik[0]] + polyofk.pofk[ik[1]] + polyofk.pofk[ik[2]]) - 1.0 / (n * n);
1126 }
1127 }
1128 }
1129 }
bool subtract_shotnoise
Subtract shotnoise? Algorithms (compute_polyspectrum) using this class checks this value.
std::pair< int, int > get_extra_slices_needed_for_density_assignment(std::string density_assignment_method)
Compute how many extra slices we need in the FFTWGrid for a given density assignement / interpolation...
void particles_to_fourier_grid_interlacing(T *part, size_t NumPart, size_t NumPartTot, FFTWGrid< N > &density_grid_fourier, std::string density_assignment_method)
Internal method.

◆ compute_polyspectrum_bincount()

template<int N, int ORDER>
void FML::CORRELATIONFUNCTIONS::compute_polyspectrum_bincount ( int  Nmesh,
PolyspectrumBinning< N, ORDER > &  polyofk 
)

This method is used by compute_polyspectrum.

It computes the number of generalized triangles of the bins needed to normalize the polyspectra up to symmetry (i.e. we only compute it for \( k_1 \leq k_2 \leq k_3 \) and only for valid triangle configurations \( k_1 + k_2 \geq k_3 \) and then set rest using symmetry). If one is to compute many spectra with the same Nmesh and binning then one can precompute N123 and set it using polyofk.set_bincount(N123) (which sets polyofk.bincount_is_set = true and avoid a call to this function) This speeds up the polyspectrum estimation by a factor of 2.

Template Parameters
NThe dimension we work in
ORDERThe order (mono = 2, bi = 3, tri = 4)
Parameters
[in]NmeshThe size of the grid we us
[out]polyofkThe binning (we compute and store the volumes of each bin, N123, in this binning).

Definition at line 1147 of file ComputePowerSpectrum.h.

1147 {
1148
1149 const auto nbins = polyofk.n;
1150 const auto klow = polyofk.klow;
1151 const auto khigh = polyofk.khigh;
1152 const auto kbin = polyofk.kbin;
1153 auto & N123 = polyofk.N123;
1154 const size_t nbins_tot = N123.size();
1155
1156 // Allocate grids
1157 std::vector<FFTWGrid<N>> N_k(nbins);
1158 for (int i = 0; i < nbins; i++) {
1159 N_k[i] = FFTWGrid<N>(Nmesh);
1160 N_k[i].add_memory_label("FFTWGrid::compute_polyspectrum_bincount::N_" + std::to_string(i));
1161 N_k[i].set_grid_status_real(false);
1162 N_k[i].fill_fourier_grid(0.0);
1163 }
1164
1165 // Set the grids
1166#ifdef USE_OMP
1167#pragma omp parallel for
1168#endif
1169 for (int i = 0; i < nbins; i++) {
1170 FFTWGrid<N> & grid = N_k[i];
1171 const double kmag2_max = khigh[i] * khigh[i];
1172 const double kmag2_min = klow[i] * klow[i];
1173 double kmag2;
1174 std::array<double, N> kvec;
1175 for (auto && fourier_index : grid.get_fourier_range()) {
1176 grid.get_fourier_wavevector_and_norm2_by_index(fourier_index, kvec, kmag2);
1177 if (not(kmag2 > kmag2_max or kmag2 < kmag2_min)) {
1178 grid.set_fourier_from_index(fourier_index, 1.0);
1179 }
1180 }
1181 }
1182
1183 // Transform to real space
1184 for (int i = 0; i < nbins; i++) {
1185 N_k[i].fftw_c2r();
1186 }
1187
1188 // We now have N_k for all bins, integrate up
1189 for (size_t i = 0; i < nbins_tot; i++) {
1190
1191 // Current values of ik1,ik2,ik3,...
1192 std::array<int, ORDER> ik = polyofk.get_coord_from_index(i);
1193
1194 // Symmetry: only do ik1 <= ik2 <= ...
1195 if (not polyofk.compute_this_configuration(ik)) {
1196 N123[i] = 0.0;
1197 continue;
1198 }
1199
1200 // Compute number of triangles in current bin
1201 // Norm represents integration measure dx^N / (2pi)^N
1202 double N123_current = 0.0;
1203 const double norm = std::pow(1.0 / double(Nmesh) / (2.0 * M_PI), N);
1204 const auto Local_nx = N_k[0].get_local_nx();
1205#ifdef USE_OMP
1206#pragma omp parallel for reduction(+ : N123_current)
1207#endif
1208 for (int islice = 0; islice < Local_nx; islice++) {
1209 for (auto && real_index : N_k[0].get_real_range(islice, islice + 1)) {
1210 if constexpr (ORDER == 2) {
1211 double N1 = N_k[ik[0]].get_real_from_index(real_index);
1212 double N2 = N_k[ik[1]].get_real_from_index(real_index);
1213 N123_current += N1 * N2;
1214 } else if constexpr (ORDER == 3) {
1215 double N1 = N_k[ik[0]].get_real_from_index(real_index);
1216 double N2 = N_k[ik[1]].get_real_from_index(real_index);
1217 double N3 = N_k[ik[2]].get_real_from_index(real_index);
1218 N123_current += N1 * N2 * N3;
1219 } else if constexpr (ORDER == 4) {
1220 double N1 = N_k[ik[0]].get_real_from_index(real_index);
1221 double N2 = N_k[ik[1]].get_real_from_index(real_index);
1222 double N3 = N_k[ik[2]].get_real_from_index(real_index);
1223 double N4 = N_k[ik[3]].get_real_from_index(real_index);
1224 N123_current += N1 * N2 * N3 * N4;
1225 } else {
1226 double Nproduct = 1.0;
1227 for (int ii = 0; ii < ORDER; ii++)
1228 Nproduct *= N_k[ik[ii]].get_real_from_index(real_index);
1229 N123_current += Nproduct;
1230 }
1231 }
1232 }
1233 FML::SumOverTasks(&N123_current);
1234 N123[i] = N123_current * norm;
1235
1236 // We cannot have less than 1 generalized triangle so put to zero if small
1237 // due to numerical noise
1238 if (N123[i] < 1.0)
1239 N123[i] = 0.0;
1240 }
1241
1242 // Set stuff not computed
1243 for (size_t i = 0; i < nbins_tot; i++) {
1244
1245 // Current values of ik1,ik2,ik3,...
1246 std::array<int, ORDER> ik = polyofk.get_coord_from_index(i);
1247
1248 // If its already computed we don't need to set it
1249 if (polyofk.compute_this_configuration(ik))
1250 continue;
1251
1252 // Find a cell given by symmetry that we have computed
1253 // by sorting ik in increasing order
1254 std::sort(ik.begin(), ik.end(), std::less<int>());
1255
1256 // Compute cell index of cell we have computed
1257 size_t index = polyofk.get_index_from_coord(ik);
1258
1259 // Set value
1260 N123[i] = N123[index];
1261 }
1262 }

◆ compute_power_spectrum()

template<int N, class T >
void FML::CORRELATIONFUNCTIONS::compute_power_spectrum ( int  Ngrid,
T *  part,
size_t  NumPart,
size_t  NumPartTotal,
PowerSpectrumBinning< N > &  pofk,
std::string  density_assignment_method,
bool  interlacing 
)

Assign particles to grid using density_assignment_method = NGP,CIC,TSC,PCS,... Fourier transform, decolvolve the window function for the assignement above.

With interlacing bin to a grid and 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. Note that with interlacing we change the particle positions, but when returning they should be in the same state as when we started.

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]NgridSize of the grid to use.
[in]partPointer to the first particle.
[in]NumPartNumber of particles on the local task.
[in]NumPartTotalTotal number of particles on all tasks.
[out]pofkThe binned power-spectrum. We required it to be initialized with the number of bins, kmin and kmax.
[in]density_assignment_methodThe density assignment method (NGP, CIC, TSC, PCS or PQS)
[in]interlacingUse interlaced grids for alias reduction. Twice as expensive, but allows us to use a smaller grid (so its actually faster if used correctly).

Definition at line 891 of file ComputePowerSpectrum.h.

897 {
898
899 // Set how many extra slices we need for the density assignment to go smoothly
900 const auto nleftright = get_extra_slices_needed_for_density_assignment(density_assignment_method);
901 const int nleft = nleftright.first;
902 const int nright = nleftright.second + (interlacing ? 1 : 0);
903
904 // Bin particles to grid
905 FFTWGrid<N> density_k(Ngrid, nleft, nright);
906 density_k.add_memory_label("FFTWGrid::compute_power_spectrum::density_k");
907
908 if (interlacing) {
909
910 // Bin to grid using interlaced grids to produce a fourier space density field
912 part, NumPart, NumPartTotal, density_k, density_assignment_method);
913 deconvolve_window_function_fourier<N>(density_k, density_assignment_method);
914
915 } else {
916
917 // Bin to grid, fourier transform and deconvolves the window function
918 particles_to_grid<N, T>(part, NumPart, NumPartTotal, density_k, density_assignment_method);
919 density_k.fftw_r2c();
920 deconvolve_window_function_fourier<N>(density_k, density_assignment_method);
921 }
922
923 // Bin up power-spectrum
924 bin_up_power_spectrum<N>(density_k, pofk);
925
926 // Subtract shotnoise
927 if (pofk.subtract_shotnoise)
928 for (int i = 0; i < pofk.n; i++) {
929 pofk.pofk[i] -= 1.0 / double(NumPartTotal);
930 }
931 }

◆ compute_power_spectrum_direct_summation() [1/2]

template<int N, class T >
void FML::CORRELATIONFUNCTIONS::compute_power_spectrum_direct_summation ( int  Ngrid,
const T *  part,
size_t  NumPart,
PowerSpectrumBinning< N > &  pofk 
)

Brute force (but aliasing free) computation of the power spectrum.

Loop over all grid-cells and all particles and add up contribution and subtracts shot-noise term. Since we need to combine all particles with all cells this is not easiy parallelizable with MPI so we assume all CPUs have exactly the same particles when this is run on more than 1 MPI tasks (so best run just using OpenMP). This method scales as \( \mathcal{O}({\rm Npart}\cdot {\rm Nmesh}^N) \) so will be slow!

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]NgridSize of the grid to use.
[in]partPointer to the first particle.
[in]NumPartNumber of particles on the local task. With MPI assumes all tasks have the same particles.
[out]pofkThe binned power-spectrum. We required it to be initialized with the number of bins, kmin and kmax.

◆ compute_power_spectrum_direct_summation() [2/2]

template<int N, class T >
void FML::CORRELATIONFUNCTIONS::compute_power_spectrum_direct_summation ( int  Ngrid,
T *  part,
size_t  NumPart,
PowerSpectrumBinning< N > &  pofk 
)

Definition at line 704 of file ComputePowerSpectrum.h.

704 {
705
706 static_assert(
707 FML::PARTICLE::has_get_pos<T>(),
708 "[compute_power_spectrum_direct_summation] Particle class needs to have positions to use this method");
709
710 // Very simple check to see if all tasks do have the same particles
711 if (FML::NTasks > 1) {
712 long long int tmp1 = NumPart;
713 FML::MinOverTasks(&tmp1);
714 long long int tmp2 = NumPart;
715 FML::MaxOverTasks(&tmp2);
716 double x = FML::PARTICLE::GetPos(part[0])[0];
717 double y = x;
720 assert_mpi(tmp1 == tmp2 and std::abs(x - y) < 1e-10,
721 "[direct_summation_power_spectrum] All tasks must have the same particles for this method "
722 "to work\n");
723 }
724
725 assert_mpi(Ngrid > 0, "[direct_summation_power_spectrum] Ngrid > 0 required\n");
726 if (NTasks > 1 and ThisTask == 0)
727 std::cout << "[direct_summation_power_spectrum] Warning: this method assumes all tasks have the same "
728 "particles\n";
729
730 const std::complex<FML::GRID::FloatType> I(0, 1);
731 const FML::GRID::FloatType norm = 1.0 / FML::GRID::FloatType(NumPart);
732
733 FFTWGrid<N> density_k(Ngrid, 1, 1);
734 density_k.add_memory_label("FFTWGrid::compute_power_spectrum_direct_summation::density_k");
735 density_k.set_grid_status_real(false);
736
737 for (auto && fourier_index : density_k.get_fourier_range()) {
738 auto kvec = density_k.get_fourier_wavevector_from_index(fourier_index);
739 FML::GRID::FloatType real = 0.0;
740 FML::GRID::FloatType imag = 0.0;
741#ifdef USE_OMP
742#pragma omp parallel for reduction(+ : real, imag)
743#endif
744 for (size_t i = 0; i < NumPart; i++) {
745 auto * x = FML::PARTICLE::GetPos(part[i]);
746 FML::GRID::FloatType kx = 0.0;
747 for (int idim = 0; idim < N; idim++) {
748 kx += kvec[idim] * x[idim];
749 }
750 auto val = std::exp(-kx * I);
751 real += val.real();
752 imag += val.imag();
753 }
754
755 std::complex<FML::GRID::FloatType> sum = {real, imag};
756 if (ThisTask == 0 and fourier_index == 0)
757 sum -= 1.0;
758 density_k.set_fourier_from_index(fourier_index, sum * norm);
759 }
760
761 // Bin up the power-spectrum
762 bin_up_power_spectrum<N>(density_k, pofk);
763
764 // Subtract shot-noise
765 if (pofk.subtract_shotnoise)
766 for (int i = 0; i < pofk.n; i++)
767 pofk.pofk[i] -= 1.0 / double(NumPart);
768 }
float FloatType
Definition: FFTWGlobal.h:29
get_vel constexpr double * GetPos(...)
void MinOverTasks(T *value)
Inplace min over tasks of a single value.
Definition: Global.h:134
void MaxOverTasks(T *value)
Inplace max over tasks of a single value.
Definition: Global.h:121
return y
Definition: Global.h:265
const double exp
Definition: Global.h:259
int NTasks
Total number of MPI tasks.
Definition: Global.cpp:34
x
Definition: test.py:18

◆ compute_power_spectrum_multipoles()

template<int N, class T >
void FML::CORRELATIONFUNCTIONS::compute_power_spectrum_multipoles ( int  Ngrid,
FML::PARTICLE::MPIParticles< T > &  part,
double  velocity_to_displacement,
std::vector< PowerSpectrumBinning< N > > &  Pell,
std::string  density_assignment_method,
bool  interlacing 
)

A simple power-spectrum estimator for multipoles in simulations.

Displacing particles from realspace to redshift-space using their velocities along each of the coordinate axes. Result is the mean of this. Deconvolving the window-function and subtracting shot-noise (1/NumPartTotal) for monopole.

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]NgridSize of the grid to use.
[in]partParticles in the form of a MPIParticle container
[in]velocity_to_displacementFactor to convert a velocity to a displacement. This is \( \frac{c}{aH(a)B} \) for peculiar and \( \frac{c}{H(a)B} \) (where \( B\) is the boxsize) for comoving velocities At \( z = 0 \) velocity_to_displacement = \( \frac{1}{100 {\rm Boxsize}} \) when Boxsize is in Mpc/h
[out]PellVector of power-spectrum binnings. The size of Pell is the maximum ell to compute. All binnings has to have nbins, kmin and kmax set. At the end Pell[ ell ] is a binning of P_ell(k).
[in]density_assignment_methodThe density assignment method (NGP, CIC, TSC, PCS or PQS) to use.
[in]interlacingUse interlaced grids for alias reduction when computing the density field

Definition at line 777 of file ComputePowerSpectrum.h.

782 {
783
784 // Sanity check
785 static_assert(FML::PARTICLE::has_get_pos<T>(),
786 "[compute_power_spectrum_multipoles] Particle class needs to have positions to use "
787 "this method");
788 static_assert(
789 FML::PARTICLE::has_get_vel<T>(),
790 "[compute_power_spectrum_multipoles] Particle class needs to have velocity to use this method");
791
792 // Set how many extra slices we need for the density assignment to go smoothly
793 // One extra slice if we use interlacing due to displacing particles by one half cell to the right
794 const auto nleftright = get_extra_slices_needed_for_density_assignment(density_assignment_method);
795 const int nleft = nleftright.first;
796 const int nright = nleftright.second + (interlacing ? 1 : 0);
797
798 // Initialize binning just in case
799 for (size_t ell = 0; ell < Pell.size(); ell++)
800 Pell[ell].reset();
801
802 // Set a binning for each axes
803 std::vector<std::vector<PowerSpectrumBinning<N>>> Pell_all(N);
804 for (int idim = 0; idim < N; idim++) {
805 Pell_all[idim] = Pell;
806 }
807
808 // Allocate density grid
809 FFTWGrid<N> density_k(Ngrid, nleft, nright);
810 density_k.add_memory_label("FFTWGrid::compute_power_spectrum_multipoles::density_k");
811
812 // Loop over all the N axes we are going to put the particles
813 // into redshift space
814 for (int idim = 0; idim < N; idim++) {
815
816 // Set up binning for current axis
817 std::vector<PowerSpectrumBinning<N>> Pell_current = Pell_all[idim];
818 for (size_t ell = 0; ell < Pell_current.size(); ell++)
819 Pell_current[ell].reset();
820
821 // Make line of sight direction unit vector
822 std::vector<double> line_of_sight_direction(N, 0.0);
823 line_of_sight_direction[idim] = 1.0;
824
825 // Transform to redshift-space
826 FML::COSMOLOGY::particles_to_redshiftspace(part, line_of_sight_direction, velocity_to_displacement);
827
828 // Bin particles to grid
829 density_k.set_grid_status_real(true);
830 if (interlacing) {
831
832 // Bin to grid and interlaced grid and deconvolve window function
834 part.get_npart(),
835 part.get_npart_total(),
836 density_k,
837 density_assignment_method);
838 deconvolve_window_function_fourier<N>(density_k, density_assignment_method);
839
840 } else {
841
842 // Bin to grid, fourier transform and deconvolve window function
843 particles_to_grid<N, T>(part.get_particles_ptr(),
844 part.get_npart(),
845 part.get_npart_total(),
846 density_k,
847 density_assignment_method);
848 density_k.fftw_r2c();
849 deconvolve_window_function_fourier<N>(density_k, density_assignment_method);
850 }
851
852 // Compute power-spectrum multipoles
853 compute_power_spectrum_multipoles_fourier(density_k, Pell_current, line_of_sight_direction);
854
855 // Copy over binning we computed
856 Pell_all[idim] = Pell_current;
857
858 // Transform particles back to real-space (we don't want to ruin the particles)
859 // Ideally we should have taken a copy, but this is fine
860 FML::COSMOLOGY::particles_to_redshiftspace(part, line_of_sight_direction, -velocity_to_displacement);
861 }
862
863 // Normalize
864 for (size_t ell = 0; ell < Pell.size(); ell++) {
865 for (int idim = 0; idim < N; idim++) {
866 Pell[ell] += Pell_all[idim][ell];
867 }
868 }
869 for (size_t ell = 0; ell < Pell.size(); ell++) {
870 for (int i = 0; i < Pell[ell].n; i++) {
871 Pell[ell].pofk[i] /= double(N);
872 Pell[ell].count[i] /= double(N);
873 Pell[ell].kbin[i] /= double(N);
874 }
875 }
876
877 // XXX Compute variance of pofk
878
879 // Subtract shotnoise for monopole
880 if (Pell[0].subtract_shotnoise)
881 for (int i = 0; i < Pell[0].n; i++) {
882 Pell[0].pofk[i] -= 1.0 / double(part.get_npart_total());
883 }
884 }
size_t get_npart_total() const
Total number of active particles across all tasks.
Definition: MPIParticles.h:441
T * get_particles_ptr()
Get a pointer to the first particle.
Definition: MPIParticles.h:451
size_t get_npart() const
Number of active particles on the local task.
Definition: MPIParticles.h:436
void compute_power_spectrum_multipoles_fourier(const FFTWGrid< N > &fourier_grid, std::vector< PowerSpectrumBinning< N > > &Pell, std::vector< double > line_of_sight_direction)
Compute power-spectrum multipoles from a Fourier grid where we have a fixed line_of_sight_direction (...
void particles_to_redshiftspace(FML::PARTICLE::MPIParticles< T > &part, std::vector< double > line_of_sight_direction, double velocity_to_displacement)
This takes a set of particles and displace them from realspace to redshiftspace Using a fixed line of...

◆ compute_power_spectrum_multipoles_fourier()

template<int N>
void FML::CORRELATIONFUNCTIONS::compute_power_spectrum_multipoles_fourier ( const FFTWGrid< N > &  fourier_grid,
std::vector< PowerSpectrumBinning< N > > &  Pell,
std::vector< double >  line_of_sight_direction 
)

Compute power-spectrum multipoles from a Fourier grid where we have a fixed line_of_sight_direction (typical coordinate axes like \( (0,0,1) \)).

Pell contains \( P_0,P_1,P_2,\ldots,P_{\ell_{\rm max}} \) where \( \ell_{\rm max} = n-1 \) is the size of Pell

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]fourier_gridThe fourier grid to compute the multipoles from.
[out]PellVector of power-spectrum binnings. The size of Pell is the maximum ell to compute. All binnings has to have nbins, kmin and kmax set. At the end Pell[ ell ] is a binning of P_ell(k).
[in]line_of_sight_directionThe line of sight direction, e.g. \( (0,0,1) \) for the z-axis.

Definition at line 425 of file ComputePowerSpectrum.h.

427 {
428
429 assert_mpi(
430 line_of_sight_direction.size() == N,
431 "[compute_power_spectrum_multipoles_fourier] Line of sight direction has wrong number of dimensions\n");
432 assert_mpi(Pell.size() > 0, "[compute_power_spectrum_multipoles_fourier] Pell must have size > 0\n");
433 assert_mpi(fourier_grid.get_nmesh() > 0,
434 "[compute_power_spectrum_multipoles_fourier] grid must have Nmesh > 0\n");
435
436 int Nmesh = fourier_grid.get_nmesh();
437 auto Local_nx = fourier_grid.get_local_nx();
438 auto Local_x_start = fourier_grid.get_local_x_start();
439
440 // Norm of LOS vector
441 double rmag = 0.0;
442 for (int idim = 0; idim < N; idim++)
443 rmag += line_of_sight_direction[idim] * line_of_sight_direction[idim];
444 rmag = std::sqrt(rmag);
445 assert_mpi(rmag > 0.0,
446 "[compute_power_spectrum_multipoles_fourier] Line of sight vector has zero length\n");
447
448 // Initialize binning just in case
449 for (size_t ell = 0; ell < Pell.size(); ell++)
450 Pell[ell].reset();
451
452 // Bin up mu^k |delta|^2
453#ifdef USE_OMP
454#pragma omp parallel for
455#endif
456 for (int islice = 0; islice < Local_nx; islice++) {
457 [[maybe_unused]] double kmag;
458 [[maybe_unused]] std::array<double, N> kvec;
459 for (auto && fourier_index : fourier_grid.get_fourier_range(islice, islice + 1)) {
460 if (Local_x_start == 0 and fourier_index == 0)
461 continue; // DC mode( k=0)
462
463 // Special treatment of k = 0 plane
464 auto last_coord = fourier_index % (Nmesh / 2 + 1);
465 double weight = last_coord > 0 and last_coord < Nmesh / 2 ? 2.0 : 1.0;
466
467 // Compute kvec, |kvec| and |delta|^2
468 fourier_grid.get_fourier_wavevector_and_norm_by_index(fourier_index, kvec, kmag);
469 auto delta = fourier_grid.get_fourier_from_index(fourier_index);
470 double power = std::norm(delta);
471
472 // Compute mu = k_vec*r_vec
473 double mu2 = 0.0;
474 for (int idim = 0; idim < N; idim++)
475 mu2 += kvec[idim] * line_of_sight_direction[idim];
476 mu2 /= (kmag * rmag);
477 mu2 = mu2 * mu2;
478
479 // Add to bin |delta|^2, |delta|^2mu^2, |delta^2|mu^4, ...
480 double mutotwoell = 1.0;
481 for (size_t ell = 0; ell < Pell.size(); ell += 2) {
482 Pell[ell].add_to_bin(kmag, power * mutotwoell, weight);
483 mutotwoell *= mu2;
484 }
485 }
486 }
487
488#ifdef USE_MPI
489 for (size_t ell = 0; ell < Pell.size(); ell++) {
490 MPI_Allreduce(MPI_IN_PLACE, Pell[ell].pofk.data(), Pell[ell].n, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
491 MPI_Allreduce(MPI_IN_PLACE, Pell[ell].count.data(), Pell[ell].n, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
492 MPI_Allreduce(MPI_IN_PLACE, Pell[ell].kbin.data(), Pell[ell].n, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
493 }
494#endif
495
496 // Normalize
497 for (size_t ell = 0; ell < Pell.size(); ell++)
498 Pell[ell].normalize();
499
500 // Binomial coefficient
501 auto binomial = [](int n, int k) -> double {
502 double res = 1.0;
503 for (int i = 0; i < k; i++) {
504 res *= double(n - i) / double(k - i);
505 }
506 return res;
507 };
508
509 // P_ell(x) = Sum_{k=0}^{ell/2} summand_legendre_polynomial * x^(ell - 2k)
510 auto summand_legendre_polynomial = [&](int k, int ell) -> double {
511 double sign = (k % 2) == 0 ? 1.0 : -1.0;
512 return sign * binomial(ell, k) * binomial(2 * ell - 2 * k, ell) / std::pow(2.0, ell);
513 };
514
515 // Go from <mu^k |delta|^2> to (2ell+1) <L_ell(mu) |delta|^2>
516 std::vector<std::vector<double>> temp;
517 for (int ell = 0; ell < int(Pell.size()); ell++) {
518 std::vector<double> sum(Pell[0].pofk.size(), 0.0);
519 for (int k = 0; k <= ell / 2; k++) {
520 std::vector<double> & mu_power = Pell[ell - 2 * k].pofk;
521 for (size_t i = 0; i < sum.size(); i++)
522 sum[i] += mu_power[i] * summand_legendre_polynomial(k, ell) * double(2 * ell + 1);
523 }
524 temp.push_back(sum);
525 }
526
527 // Copy over data. We now have P0,P1,... in Pell
528 for (size_t ell = 0; ell < Pell.size(); ell++) {
529 Pell[ell].pofk = temp[ell];
530 }
531 }

◆ compute_trispectrum() [1/3]

template<int N>
void FML::CORRELATIONFUNCTIONS::compute_trispectrum ( const FFTWGrid< N > &  fourier_grid,
PolyspectrumBinning< N, 4 > &  tofk 
)

Definition at line 1275 of file ComputePowerSpectrum.h.

1275 {
1276 compute_polyspectrum<N, 4>(fourier_grid, tofk);
1277 }

◆ compute_trispectrum() [2/3]

template<int N>
void FML::CORRELATIONFUNCTIONS::compute_trispectrum ( const FFTWGrid< N > &  fourier_grid,
TrispectrumBinning< N > &  tofk 
)

Computes the trispectrum \( T(k_1,k_2,k_3,k_4) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\delta(k_4)\right> \) from a fourier grid.

This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{4} \) integrals. This is just an alias for compute_polyspectrum<N, 4>

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]fourier_gridA fourier grid.
[out]tofkThe binned trispectrum. We required it to be initialized with the number of bins, kmin and kmax.

◆ compute_trispectrum() [3/3]

template<int N, class T >
void FML::CORRELATIONFUNCTIONS::compute_trispectrum ( int  Ngrid,
T *  part,
size_t  NumPart,
size_t  NumPartTotal,
TrispectrumBinning< N > &  tofk,
std::string  density_assignment_method,
bool  interlacing 
)

Computes the truspectrum \( T(k_1,k_2,k_3,k_4) = \left<\delta(k_1)\delta(k_2)\delta(k_3)\delta(k_4)\right> \) from particles.

Note that with interlacing we change the particle positions, but when returning they should be in the same state as when we started. This method allocates nbins FFTWGrids at the same time and performs \( 2{\rm nbins} \) fourier transforms and does \( {\rm nbins}^{4} \) integrals. This is just an alias for compute_polyspectrum<N, 4> The shot-noise term (something like (B(123) + cyc)/n + (P1P2 + cyc)/n^2 + 1/n^3) is not subtracted as it requires the bispectrum which we currently don't bin up together with the trispectrum

Template Parameters
NThe dimension of the particles.
TThe particle class. Must have a get_pos() method.
Parameters
[in]NgridSize of the grid to use.
[in]partPointer to the first particle.
[in]NumPartNumber of particles on the local task.
[in]NumPartTotalTotal number of particles on all tasks.
[out]tofkThe binned trispectrum. We required it to be initialized with the number of bins, kmin and kmax.
[in]density_assignment_methodThe density assignment method (NGP, CIC, TSC, PCS or PQS)
[in]interlacingUse interlaced grids when computing density field for alias reduction

Definition at line 1057 of file ComputePowerSpectrum.h.

1063 {
1064 compute_polyspectrum<N, T, 3>(
1065 Ngrid, part, NumPart, NumPartTotal, tofk, density_assignment_method, interlacing);
1066 }