FFTPower mu-edges
adematti opened this issue · 3 comments
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
In most cases one will paint positions and weights on a real mesh
nbodykit/nbodykit/base/catalog.py
Line 787 in 4aec168
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
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
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
- I'd have naively taken coordinates at the grid nodes, without this cell-shift. Am I missing something?
Thanks!
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?
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
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
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).
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
