FFTPower mu-edges

[adematti]

Hi,

I think there is an issue with current FFTPower implementation (0.3.16dev0);
muedges (now) range from -1 to 1 as set in:

nbodykit/nbodykit/algorithms/fftpower.py

Line 298 in 4aec168

muedges = numpy.linspace(-1, 1, self.attrs['Nmu']+1, endpoint=True)

In most cases one will paint positions and weights on a real mesh

nbodykit/nbodykit/base/catalog.py

Line 787 in 4aec168

def to_mesh(self, Nmesh=None, BoxSize=None, dtype='f4', interlaced=False,

Hence, hermitian symmetry is exploited to use only kz >= 0 modes.
When line-of-sight is z, mu obtained by

nbodykit/nbodykit/algorithms/fftpower.py

Line 620 in 4aec168

mu = slab.mu(los) # defined with respect to specified LOS

will always be > 0. Hence wedges at mu < 0 are empty (NaN power).
I guess one could either (by order of increased code changes):
a) drop Hermitian symmetry by defaulting to complex field as has been done for FKPCatalog, but this is overkill
b) revert to muedges in [0, 1], and take abs(mu) in case of Hermitian symmetry
c) keep muedges in [-1, 1], and create mu < 0 modes on-the-fly in project_to_basis.

Note: about odd poles in FKP estimator, instead of switching to full complex field, one can just use Hermitian antisymmetry in project_to_basis.
Note: in case single precision is used, some (k, mu) modes fall in the wrong bin. Casting slab coordinates to 'f8' around here:

nbodykit/nbodykit/algorithms/fftpower.py

Line 608 in 4aec168

xslab = slab.norm2()

resolves this problem and increases accuracy in the estimated power spectrum at low memory cost.

Another, quite unrelated question: I though that, within pmesh, a particles are "painted" on the mesh nodes,
i.e. not in cell centers. So I do not understand that extra half-cell shift (0.5*pm.BoxSize / pm.Nmesh) here:

nbodykit/nbodykit/algorithms/convpower/fkp.py

Line 457 in 4aec168

offset = self.attrs['BoxCenter'] + 0.5*pm.BoxSize / pm.Nmesh

• I'd have naively taken coordinates at the grid nodes, without this cell-shift. Am I missing something?

Thanks!

[rainwoodman]

Thanks for the reports.

Yes I agree there seems to be a NaN problem. Do you have a simple notebook to reproduce this?

For a fix, perhaps this line:

6ab5a06#diff-bf846022d2da880aa1441037073fa99e26bdde9b7b720c4cc2e15504a3d95aa4R571

shall used y3d.compressed instead of False. If compressed is True then there is no mu < 0 values, and we shall use the hermitian symmetric block. Setting it to False regardless of the input is compressed was incorrect.

The 'symmetric' vs 'anti-symmetric' treatment you are probably referring to this?
6ab5a06#diff-bf846022d2da880aa1441037073fa99e26bdde9b7b720c4cc2e15504a3d95aa4R648

Round off errors(f8): That's reasonable. Could you file a Pull request?

Offset of half of grid: Indeed that shift is weird. I agree with you that eliminating it makes sense. Did you check how much eliminating the shift affects the results?

[adematti]

Sorry for the delay,

The line you quote (6ab5a06#diff-bf846022d2da880aa1441037073fa99e26bdde9b7b720c4cc2e15504a3d95aa4R571) has been reverted back in ede93aa.
If compressed is True then there is no mu < 0 values => that depends on whether the line of sight is the same as the symmetry axis (x, y, z).
The issue appears when line-of-sight is z, in which case there are no mu < 0 values, hence empty mu < 0 wedges.
In other cases (line-of-sights x, y), the result looks about correct (even in this case there is actually some non-negligible difference of power between c16 and f8).

The issue is shown by this code snippet:

from nbodykit.lab import UniformCatalog, FFTPower

catalog = UniformCatalog(nbar=1e-4, BoxSize=1000.)
mesh = catalog.to_mesh(Nmesh=128, resampler='cic', interlaced=True, compensated=True)
power = FFTPower(mesh, mode='2d', Nmu=5, los=[0, 0, 1], dk=0.01)
print(power.power['power'][:,0])

For the 'symmetric' vs 'anti-symmetric' treatment, this is indeed linked to 6ab5a06.
But instead of using a complex mesh, one can stick to a real mesh and use that, given the quantity:
X_ell(k) = F_0(k) F_ell(k)* that is computed here:

nbodykit/nbodykit/algorithms/convpower/fkp.py

Line 635 in 4aec168

A0_1[islab,...] = norm*A0_1[islab]*A0_2[islab].conj()

for odd poles X_ell(-k) = - X_ell(k)*.

Actually, casting coordinates to double precision here:

nbodykit/nbodykit/algorithms/fftpower.py

Line 570 in 4aec168

x3d = y3d.x

solves only ~ half of the wrong bin assignments. I guess that should be done when computing these coordinates, within pmesh itself?

About the shifts in ConvolvedFFTPower, it's a typical ~ 5% effect for the quadrupole, in the test case of https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/tests/test_conv_power.py
(of course, no effect on the monopole, and a bit more for the hexadecapole).

[jegarfa]

Hi, I just arrived to the same issue.
I was trying to reproduce one of the documentation examples, specifically the one that computes the 2D Power Spectrum.
I noticed the same behavior previous commented, the mu_edges are starting from -1 instead of 0.
Here I attached a notebook that reproduce the issue
[example notebook](https://github.com/jegarfa/Example_nbodykit/blob/main/Example_2D_Power_Spectrum.ipynb)

I am not sure if just by editing this line would be enough
6ab5a06#diff-bf846022d2da880aa1441037073fa99e26bdde9b7b720c4cc2e15504a3d95aa4R571