meerklass / meerpower

Code repo for power spectrum analysis pipeline
MIT License
0 stars 0 forks source link

Optimise map re-gridding and model the effects #5

Open stevecunnington opened 1 year ago

stevecunnington commented 1 year ago

Re-gridding map observations from sky coordinates (RA,Dec,z) onto comoving (Mpc/h) grid, appears to be causing a lot of damped and distorted power. The attached plot demonstrates the problem using lognormal mocks. Generating a mock straight onto the Cartesian grid and then measuring the power recovers the input as expected (blue points). But re-gridding the mock onto the sky, then sampling onto a Cartesian grid causes a loss and distortion of power. I currently use Nearest Grid Point mass assignment in the re-gridding and sample onto a grid with half the pixels in each dimension as the sky map, but enclosed on the footprint. Increasing the number of test particles per pixel in the re-gridding does not make this power damping go away. The errors here are from the variance over generation of 5 lognormal mocks. image

This might have something to do with aliasing effects which are more significant for us due to large RA,Dec pixels? There is literature on this for galaxies (e.g. Jing+05 https://arxiv.org/pdf/astro-ph/0409240.pdf, Blake+10 https://arxiv.org/pdf/1003.5721.pdf and Sefusatti https://arxiv.org/pdf/1512.07295.pdf) but I haven't found any detailed study for intensity mapping and pixelised maps (Blake19 https://arxiv.org/pdf/1902.07439.pdf discusses applying it to IM but no isolated study of the effects).

Implementing aliasing corrections seems to help (green points). The damped model (grey) includes convolving with the re-gridded window and modelling of the finite size of frequency channels and pixel size. This seems to provide good agreement but a more detailed investigation to understand this fully and devise optimal re-gridding strategies is needed.

Alkistis commented 1 year ago

@stevecunnington this is very odd indeed ... I have never seen damping of power because of aliasing/interlacing, only the increase from k = 0.2 onwards looks familiar (what is the Nyquist frequency?). Is it possible that we are looking at 2 effects here, can it be there is a normalisation factor associated with the re-gridding missing somewhere?

stevecunnington commented 1 year ago

I'm sure there are many things going on here yes. The maps in this test undergo 2 coordinate transforms, cartesian->sky, then back to cartesian, with re-gridding at each step. So there's probably something missing. It's just strange that the aliasing correction and simple damping model appears to be handling it! But we should definitely understand why

The current cartesian grid I choose is quite low resolution to avoid empty cells. It seems to work well from visual inspection but does result in a very low Nyquist frequency for the angular pixels: k_nyq_x = 0.1983 k_nyq_y = 0.1722 k_nyq_z = 1.3451

So we measure right up to and beyond the angular nyquist frequency! It surprises me that this even works but I suppose the higher k_nyq_z allows power to be resolved at k>0.2.

stevecunnington commented 1 year ago

I feel like there should be a dedicated paper on this stuff. If not someone should write one...

Alkistis commented 1 year ago

I'm sure there are many things going on here yes. The maps in this test undergo 2 coordinate transforms, cartesian->sky, then back to cartesian, with re-gridding at each step. So there's probably something missing. It's just strange that the aliasing correction and simple damping model appears to be handling it! But we should definitely understand why

The current cartesian grid I choose is quite low resolution to avoid empty cells. It seems to work well from visual inspection but does result in a very low Nyquist frequency for the angular pixels: k_nyq_x = 0.1983 k_nyq_y = 0.1722 k_nyq_z = 1.3451

So we measure right up to and beyond the angular nyquist frequency! It surprises me that this even works but I suppose the higher k_nyq_z allows power to be resolved at k>0.2.

Thanks for explaining this! I'm sure it will be sorted :-)

Alkistis commented 1 year ago

I feel like there should be a dedicated paper on this stuff. If not someone should write one...

Maybe worth talking to Chris B about this by the way!

laurawolz commented 1 year ago

Very interesting. So do you think, measuring so closely to the angular frequency is part of the problem? Can you point me to a paper/equations of the green points, your dampening corrections? Is it the correction from the Blake10 paper?

For our data, we would only have one transformation, from sky to cartesian, right? Would it be worth generating mocks in sky coordinates to isolate the effect of the second transformation only?

stevecunnington commented 1 year ago

Yes, we're right at the Nyquist frequency for the angular scales. But I really don't know if having a good resolution along the LoS is enough to avoid severe aliasing. These preliminary tests suggest not(?) I think Chris' 2010 paper covers it best: eq14 in https://arxiv.org/pdf/1003.5721.pdf. It's essentially summing over mass-assignment functions (eq15) each time shifted by integer Nyquist frequencies. The Jing05 paper was the original dedicated paper on this. I think galaxy surveys may now prefer interlacing techniques as Alkistis mentioned but I don't really understand these nor know if they're applicable to pixelised intensity maps.

You're right on this second point. I'm undergoing two transformations when in reality it's one. Almost impossible to completely avoid this in sims though as they're always generated in physical space, so would just need to have a sim with a coordinate transformation and mass assignment scheme that's maximally "lossless", then test our mass assignment regridding back to comoving methods on this sim.

laurawolz commented 1 year ago

Naive question. Why can you not keep the original pixelisation of the sky map in RA/DEC and only transform the dimensions to comoving units?

stevecunnington commented 1 year ago

Because even for our relatively small cube, the volume footprint is slightly "conical" in comoving space, so just estimating the L_x,L_y,L_z of the (RA,Dec,z) cube and FFTing, would only be an approximation. I have some maps handy from the demo script to demonstrate exactly this point...

image