Brute force convolution test?

mreineck commented 2 years ago

I'm thinking about contributing a test which computes the measured detector signal for a single sample by brute force, according to the following recipe:

given: s_lm and b_lm (T,E,B,V components); M_{HWP} (4x4 real numbers); theta, phi, psi, alpha
rotate b_lm to the proper orientation using healpy.rotator.rotate_alm() (uses b_lm, theta, phi, psi)
convert s_lm to maps (probably Gauss-Legendre to avoid inaccuracies)
apply the appropriate Mueller matrix to the sky maps (uses M_{HWP}, alpha)
convert rotated b_lm to maps (probably Gauss-Legendre to avoid inaccuracies)
multiply sky and beam maps point-wise
integrate product maps over the sphere and return sum over components.

This is of course an extremely expensive way to compute a single sample, but it may be useful to verify the beamconv results for fully general Mueller matrices. The same sort of test (without HWP) is used in https://github.com/mreineck/ducc/blob/135ef6fb68b231cfbd862ab5351f0c9ef9afc44f/python/test/test_totalconvolve.py#L107-L134 and has served me well so far.

Do you think such a test would be a useful addition to beamconv?

AdriJD commented 2 years ago

Hi Martin, I think this would a really valuable test! We have done similar tests before we implemented the non-ideal HWPs (see offset_beam in test.py) but the test you are suggesting sounds more straightforward and would test the HWP details better. Let me know if you would like me to help.

mreineck commented 2 years ago

Thanks! I will work on this next week and ask you for help when (not if) I get stuck!

mreineck commented 2 years ago

Here is a rough first sketch demonstrating what I have in mind. To get exact spherical harmonic transforms I had toi make it depend on ducc0 instead of healpy, but this can be changed if this dependency is too cumbersome.

I'm certain that there are still a lot of mistakes in the details, but do you think this is going in the right direction?

import numpy as np
import ducc0

def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)

def random_alm(rng, lmax, mmax, ncomp):
    res = rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax))) \
     + 1j*rng.uniform(-1., 1., (ncomp, nalm(lmax, mmax)))
    # make a_lm with m==0 real-valued
    res[:, 0:lmax+1].imag = 0.
    return res

def convolve(alm1, alm2, lmax, nthreads=1):
    # Go to a Gauss-Legendre grid of sufficient dimensions
    ntheta, nphi = lmax+1, ducc0.fft.good_size(2*lmax+1, True)
    tmap = ducc0.sht.experimental.synthesis_2d(
        alm=alm1.reshape((1,-1)), ntheta=ntheta, nphi=nphi, lmax=lmax,
        geometry="GL", spin=0, nthreads=nthreads)
    tmap *= ducc0.sht.experimental.synthesis_2d(
        alm=alm2.reshape((1,-1)), ntheta=ntheta, nphi=nphi, lmax=lmax,
        geometry="GL", spin=0, nthreads=nthreads)
    # compute the integral over the resulting map
    res = ducc0.sht.experimental.analysis_2d(
        map=tmap, lmax=0, spin=0, geometry="GL", nthreads=nthreads)
    return np.sqrt(4*np.pi)*res[0,0].real

def apply_mueller_to_sky(slm, lmax, m_hwp, alpha, nthreads=1):
    ncomp = slm.shape[0]

    # Go to a Gauss-Legendre grid of sufficient dimensions
    ntheta, nphi = lmax+1, ducc0.fft.good_size(2*lmax+1, True)
    skymap = np.zeros((ncomp, ntheta, nphi))
    skymap[0:1] = ducc0.sht.experimental.synthesis_2d(
        alm=slm[0:1], ntheta=ntheta, nphi=nphi, lmax=lmax,
        geometry="GL", spin=0, nthreads=nthreads)
    if ncomp >= 3:
        skymap[1:3] = ducc0.sht.experimental.synthesis_2d(
            alm=slm[1:3], ntheta=ntheta, nphi=nphi, lmax=lmax,
            geometry="GL", spin=2, nthreads=nthreads)
    if ncomp == 4:
       skymap[3:4] = ducc0.sht.experimental.synthesis_2d(
           alm=slm[3:4], ntheta=ntheta, nphi=nphi, lmax=lmax,
           geometry="GL", spin=0, nthreads=nthreads)

    # apply Mueller matrix to sky
    m_alpha = np.zeros((4,4))
    m_alpha[0,0] = m_alpha[3,3] = 1
    m_alpha[1,1] = m_alpha[2,2] = np.cos(2*alpha)
    m_alpha[1,2] = np.sin(2*alpha)
    m_alpha[2,1] = -np.sin(2*alpha)
    m_hwp_alpha = m_alpha.T.dot(m_hwp.dot(m_alpha))
    # there must be a better way to do this ...
    for i in range(ntheta):
        for j in range(nphi):
            skymap[:, i, j] = m_hwp_alpha.dot(skymap[:, i, j])

    # go back to spherical harmonics
    res = np.empty_like(slm)
    res[0:1] = ducc0.sht.experimental.analysis_2d(
        map=skymap[0:1], lmax=lmax, spin=0, geometry="GL", nthreads=nthreads)
    if ncomp >= 3:
        res[1:3] = ducc0.sht.experimental.analysis_2d(
            map=skymap[1:3], lmax=lmax, spin=2, geometry="GL", nthreads=nthreads)
    if ncomp == 4:
        res[3:4] = ducc0.sht.experimental.analysis_2d(
            map=skymap[3:4], lmax=lmax, spin=0, geometry="GL", nthreads=nthreads)
    return res

def explicit_convolution(slm, blm, lmax, theta, phi, psi, alpha, m_hwp, nthreads=1):
    # extend blm to full (lmax, lmax) size for rotation
    blm2 = np.zeros((blm.shape[0], nalm(lmax, lmax)), dtype=np.complex128)
    blm2[:, 0:blm.shape[1]] = blm
    # rotate beam to desired orientation
    for i in range(blm2.shape[0]):
        blm2[i] = ducc0.sht.rotate_alm(blm2[i], lmax, psi, theta, phi, nthreads)
    # apply Mueller matrix to sky
    slm2 = apply_mueller_to_sky(slm, lmax, m_hwp, alpha, nthreads)

    # convolve sky and beam component-wise and return the sum
    res = 0.
    for i in range(blm2.shape[0]):
        res += convolve(slm2[i], blm2[i], lmax, nthreads)
    return res

def main():
    lmax = 100
    bmmax = 10
    ncomp = 4
    nthreads = 1

    rng = np.random.default_rng(42)
    slm = random_alm(rng, lmax, lmax, ncomp)
    blm = random_alm(rng, lmax, bmmax, ncomp)

    theta, phi, psi = 0.2, 0.3, 0.4
    alpha = 0.7
    m_hwp = np.array([[1., 0., 0., 0.],
                      [0., 1., 0., 0.],
                      [0., 0.,-1., 0.],
                      [0., 0., 0.,-1.]])
    res = explicit_convolution(slm, blm, lmax, theta, phi, psi, alpha, m_hwp, nthreads)
    print(res)

if __name__ == '__main__':
    main()

mreineck commented 2 years ago

I have added a Python script with the first beginnings of such a test at https://github.com/mreineck/beamconv/tree/low_level_tests. There are still a few problems, most likely because my application of the Mueller matrix is wrong... for example if I change the (2,2) element to -1, the returned values differ. I'm investigating.

Questions:

what do I have to do to enable convolution of the V component? I tried passing input_v=True and beam_v=True to various functions, but so far without success
can I somehow trick beamconv into accepting a set of theta, phi, pa, hwp_angle that I have prepared in numpy arrays? This could make systematic testing a bit easier.

mreineck commented 2 years ago

I have added some documentation to the test code. At the moment I'm a tmy wits' end concerning the problem with the Mueller matrix. Any insights would be more than welcome!

mreineck commented 2 years ago

One observation, which is a bit disconcerting (at least to me):

An ideal HWP with no rotation basically flips the sign of the Stokes U component (let's ignore V for the moment) and leaves Stokes I and Q unchanged, correct?

I just observed that if I take a band-limited Stokes IQU map and flip the sign of U, the resulting map is not band-limited any more. If that's correct, my approach (which involves synthesizing a sky map from the a_lm, applying the Mueller matrix pixel-wise, and then analyzing the resulting map) cannot work ...

AdriJD commented 2 years ago

Hi Martin,

Sorry for reacting so slowly, I should have more time after coming Thursday.

About your last point: I agree that the Mueller matrix of an ideal HWP is given by diag(1, 1, -1, -1). And I think you're also right about the band-limit increasing when you change U to -U. I tested it with GL pixels and with Healpix pixels and it happens in both cases. The extra signal seems to be localized around the poles. I can't really wrap my head around why that is happening. Q -> Q, U -> -U could be interpreted as a parity transformation, but I don't see how that increases the band-limit. Or, you could interpret Q -> Q, U -> -U as a complex conjugation of P = Q + iU, but I still don't see how that increases the band-limit. Let's perhaps discuss this tomorrow when we meet.

I don't think it's very easy to trick beamconv to use theta and phi directly. I'll let you know if I find a way to do it.

Another thing for your test: perhaps it would be interesting to think about applying the rotated HWP matrix to the beam before rotating it to the position on the sky instead of applying it to the sky. I think that might make the comparison to beamconv easier to interpret, since in beamconv the HWP angle is defined with respect to a coordinate frame fixed to the telescope. When using that definition one should rotate the HWP angle with respect to the sky when rotating the beam.

mreineck commented 2 years ago

Concerning the loss of the band limit, this could be a somewhat intuitive example:

Assume that the map is using a pixelization with a ring directly at the North Pole, i.e. where all pixels on the first ring actually fall in the same place. In that case, Q and U in these pixels must follow Q(phi) + i U(phi) = exp(2*i*phi)(Q(0) + iU(0)). Most "realistic" Mueller matrices will destroy this property and produce an inconsistent map, if applied in the way I'm currently doing it. So of course I hope that what I'm doing is wrong :-)

mreineck commented 2 years ago

I have switched the test to the complex-valued variant including the T matrices. Unfortunately this doesn't seem to change anything.

AdriJD / beamconv

Brute force convolution test? #8