fixed queb_rotmat, handles all cdelt cases, made tests more general

[zatkins2]

This fixes #288

The fundamental issue is that the previous enmap.queb_rotmat implementation did not respect symmetries required by the DFT of a real map. This induced E<->B mixing at the Nyquist frequency, where these symmetries were being violated by the sine term of the rotmat. This fix respects the symmetries, and also is consistent for all possible variations of the signs of cdeltx and cdelty and well as even/odd map shapes in x and y. I added a new test to catch this.

Handling the possibilities of non-standard cdelt was important to pass test_b_sign in the existing tests -- basically, in the case that a user flips their map in between performing map2harm and harm2map operations. I generalized this test for all possible variations of the signs of cdeltx and cdelty and well as even/odd map shapes in x and y

[codecov]

Codecov Report

Attention: Patch coverage is 89.47368% with 2 lines in your changes missing coverage. Please review.

Project coverage is 39.92%. Comparing base (c0c538d) to head (19230e2).

Files with missing lines	Patch %	Lines
pixell/multimap.py	0.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #289      +/-   ##
==========================================
+ Coverage   39.80%   39.92%   +0.11%     
==========================================
  Files          30       30              
  Lines       10439    10453      +14     
==========================================
+ Hits         4155     4173      +18     
+ Misses       6284     6280       -4     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[zatkins2]

I should mention, I also updated a few things with cases in test_b_sign that were already tested. Namely, I properly handled the input power spectra (they are Dl's starting from l=2 on disk, but were treated as Cl's starting at l=0. This is only important for setting the numerical precision target of the final comparison, but basically left untreated made the tests much harder to pass), and I made the numerical precision target of the final comparison more aggressive. It's kind of arbitrary what to set it to but a relative tolerance of 1e-7 is fine I think.

[zatkins2]

Sorry last thing: I updated the queb_rotmat in multimap to just use the one in enmap rather than copying the same code twice.

[amaurea]

I like this solution to the CDELT sign problem. When noticed it early on when writing enmap, I thought I would have to flip U map's sign any time the pixel ordering was changed, if I wanted things to be consistent, and that would be slow (e.g. suddenly slicing would create a copy). I don't remember why I thought just changing queb_rotmat wouldn't be enough, but it definitely is, and is also the physically reasonable thing to do (U should be defined with respect to coordinate axes, not pixel axes). So good job!

The implementation looks very verbose, though! Does it really have to be that complicated?

I haven't finished thinking about this, but maybe it would be cleaner to have this in laxes instead of queb_rotmat.

[zatkins2]

Feel free to suggest something equivalent and simpler; mathematically I am pretty confident each of these cases must be handled as they are

[zatkins2]

Perhaps we can push and open an issue to consider simplifying in the meantime? I think it's better to have the fix (with tests) pushed, and optimize later.

Fwiw, I think keeping this in queb_rotmat is appropriate -- laxes inherits the conventional numpy behavior which I think should be expected by downstream users. It would be confusing/potentially-breaking to have to respect a different laxes convention. The numpy convention should be handled on a case-by-case basis by laxes consumers as-needed, as in this case.

[amaurea]

I found some time to look at this today. This is surprisingly tricky, so I'm not finished, but I've still found some interesting stuff.

The main thing is that the test in test_queb_rotmat is inconsistent. It uses harm2map without the keep_imag option, but it doesn't ensure that the fourier coefficients it generates represent a real map. For a real map, the f=0 must be real, but additionally the Nyquist frequency must be real for even array lengths. Those cases make up 3/4 of the cases in test. The E/B mxing the test identifies is really just the effect of truncating the complex map in real-space.

I have modified the test into two variants, one that allows the map to be complex, and another that makes sure it really is real. Here they are:

    def test_queb_rotmat_complex(self):
        # Tests that the rotmat respects fft symmetry constraints
        # for real maps -- ie, map2harm, harm2map produces·
        # sensible round-trips. This version tests a map that's
        # complex in real-space
        ishapes = [(10, 10), (10, 11), (11, 10), (11, 11)]
        dtypes = [np.complex64, np.complex128]
        comps = [1, 2] # this is the comp we will set to non-zero

        for ishape, dtype, comp in itertools.product(ishapes, dtypes, comps):
            shape, wcs = enmap.geometry(pos=(0, 0), shape=ishape, res=0.01)

            # define some easy test input to evaluate for mixing
            input_complex_hmap = enmap.zeros((3, *shape), wcs, dtype=dtype)
            input_complex_hmap[comp] += 1. + 1.j
            output_complex_hmap = enmap.zeros((3, *shape), wcs, dtype=dtype)
            output_complex_hmap[comp] += 1. + 1.j

            # do a round-trip and check that we get what we expect.
            # the queb operations mix E and B and could lead to artifacts if
            # we're not careful
            if input_complex_hmap.real.dtype.itemsize == 8:
                atol = 1e-10
            elif input_complex_hmap.real.dtype.itemsize == 4:
                atol = 1e-5
            test_output_complex_hmap = enmap.map2harm(
                enmap.harm2map(input_complex_hmap, keep_imag=True)
                )
            assert np.allclose(
                test_output_complex_hmap, output_complex_hmap, rtol=0, atol=atol
                ), f'{ishape=}, {dtype=}, {comp=}'

    def test_queb_rotmat_real(self):
        # Tests that the rotmat respects fft symmetry constraints
        # for real maps -- ie, map2harm, harm2map produces·
        # sensible round-trips. This version tests a map that's
        # real in real-space. This means that its DC must be real
        # and along even dimensions the nyquist freq must also be real
        ishapes = [(10, 10), (10, 11), (11, 10), (11, 11)]
        dtypes = [np.complex64, np.complex128]
        comps = [1, 2] # this is the comp we will set to non-zero
        np.random.seed(0)

        for ishape, dtype, comp in itertools.product(ishapes, dtypes, comps):
            shape, wcs = enmap.geometry(pos=(0, 0), shape=ishape, res=0.01)

            # define some easy test input to evaluate for mixing
            input_hmap = enmap.map2harm(enmap.rand_gauss((3,)+shape, wcs, dtype=dtype))
            # do a round-trip and check that we get what we expect.
            atol = 1e-10 if dtype == np.complex128 else 1e-5
            output_hmap = enmap.map2harm(enmap.harm2map(input_hmap))
            assert np.allclose(
                input_hmap, output_hmap, rtol=0, atol=atol), f'{ishape=}, {dtype=}, {comp=}'

Both of these tests pass without the modifications to queb_rotmat. As far as I can see what the complicated logic there was doing was basically compensating for this truncation of the imaginary part of the Nyquist element for even-length array that you get without keep_imag when operating on a map that isn't real to begin with.

When it comes to test_b_sign, that's the one I'm currently stuck on. It almost works as is, except for the Nyquist-Nyquist element for doubly even dimensions, where there's a sign flip. The whole rest of Fourier space works. A bit more investigation is needed here.