Lens map curved with lenspyx

[zatkins2]

Adds lens_map_curved_lenspyx function to pixell.lensing. This would improve things in 2 ways:

1. I believe it would fix https://github.com/simonsobs/pixell/issues/283
2. I think it would be better to delegate the curved-sky lensing operation to lenspyx, since that is "the" optimized/specialized CMB lensing library, rather than try to maintain our own method in parallel.

One thing to note is that this new function has two (minor) limitations that are not fundamental:

1. It relies on a helper function to perform surgery on the lenspyx output maps, since they don't natively know about the desired geometry. This output function should work in most cases but would fail if https://github.com/simonsobs/pixell/issues/202 is not fixed for edge cases
2. It doesn't handle all possible spin permutations, but these could be added and caught. It's just that for now, the default [0, 2] I think will cover 99% of use-cases.

The lensing output maps aren't identical (see below for T, E, B difference maps for a standard cosmology) but their power spectra are very similar (lmax=2700 for phi_alm and cmb_alm, 4 arcmin pixels on full-sky fejer1 geometry), see below (note np.max(np.abs(lenspyx - pixell)) is ~(20, 14, 12) uK):

[zatkins2]

The test failures are due to an incompatibility of numpy 2 and lenspyx

[carronj]

hopefully the pip version of lenspyx also should work with numpy 2 now.

(I dont know if that's relevant, but parsing the comments, you can also directly call lenspyx for 'non-full-sky' geometries, but adapting the geometry input)

[zatkins2]

Thanks @carronj! I think the most recent version of lenspyx on pypi is 2.0.5 from April 2023, so pip would not capture the numpy 2 fixes currently. I think that's what's causing the test issues, but I agree a new lenspyx version on pypi would probably fix this!

For the non-full-sky geometries, do you mean by using restrict and the pbdGeometry in lenspyx.remapping.utils_geom.py? If it's not too much trouble, could you give a demo of using these?

[carronj]

lenspyx pypi has been updated.

(I believe SO might want only a subset of rings of Fejer-like or similar ? -- if so there is no need to compute all rings (which I thought was what was being done), but only those of interest. That's all what I meant. If this understanding is correct, maybe describe to me a bit more precisely what's needed and I can do a demo)

[zatkins2]

Yeah that's right, SO maps (at least for the LAT) will will be a subset of rings in Fejer 1, so agreed no need to compute all the rings (which I think would be accomplished using the restrict geometry method right?). I only also brought up the pbdGeometry because in principle it would be nice if this pixell function could also handle maps that don't include all the 2pi phis as well -- i.e., a proper "square cutout" of the sky. I'm not sure if the pbdGeometry is the right way of handling that but it looked like it may be?

[carronj]

yes one can use 'restrict' for that for example. pbdGeometry is a residual of earlier attempts where I was actually using proper cutouts, which I abandoned, so it is useless in this public version. The SHT gain is absolutely minimal in this case. You'd gain a little bit more for 'general SHT's' I guess, but not a dramatic gain either (unless you have reasons to think otherwise).

[zatkins2]

That makes sense. In that case no need for a demo, I'll add functionality that uses restrict as well.

[amaurea]

Hm, I was planning on implementing this myself, directly using ducc. I already did this for aberration. When I do that, there won't be any advantage to the lenspyx implementation. But I don't have time to do this now that LAT is looming... So I might accept this as an interim solution.

I feel bad about dragging my feet with this now that you've had to implement it yourself :/

That said:

1. Only import lenspyix in the function that uses it, so we avoid having all the lensing break if lensypix isn't installed
2. Use tabs, but you can't use them to align things, since you don't know how big the tabstop will be in another user's text editor. The way you've done it will look good with your choice of tabstop, but not with mine. Instead of trying to align the second line to start after the ( from the line above, just start the line with an extra indent.

[codecov[bot]]

Codecov Report

:x: Patch coverage is 33.57664% with 91 lines in your changes missing coverage. Please review. :white_check_mark: Project coverage is 37.44%. Comparing base (15ea917) to head (8cc3a30). :warning: Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
pixell/lensing.py	33.57%	91 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #284      +/-   ##
==========================================
- Coverage   37.67%   37.44%   -0.23%     
==========================================
  Files          37       37              
  Lines       11318    11402      +84     
==========================================
+ Hits         4264     4270       +6     
- Misses       7054     7132      +78     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[zatkins2]

Sorry, I am returning to this after a while, but I think it's good to go now.

I addressed your comments and made it work on requested cut skies. I did not use restrict, I just extracted from the full-sky result. I'll add an issue to track that restrict might be faster.

[zatkins2]

@amaurea could you approve and merge?

[msyriac]

@zatkins2 tests seem to be failing with the new edits

[zatkins2]

It's good now, it's because I removed the maplmax argument, but the tests still called it. (The maplmax argument, although in lensing.py isn't used anywhere)

[msyriac]

Just looking at the spectra again, what is going on with BB at high ell? I don't remember seeing such a large discrepancy with the pixell implementation. Have you compared BB between pixell and lenspyx with their default accuracies? (I'm also curious how much faster lenspyx is)

[msyriac]

Ah ok, I see the comparison with pixell now, and understand this depends on the lmax you did your sim at. Merging.