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 lensing.lens_map_curved sometimes crashes for single-precision maps #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 pix2sky breaks at RA = +/- 180 due to reliance on floating-point arithmetic #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.