Add Müller convolver

[ziotom78]

This PR is a WIP to implement proper HWP+beam convolution, following what has been discussed in issue #13. We move the discussion here so that it's easier to comment on specific parts of the code.

[mreineck]

I have added a testcase which (at least I assume) shows that something is still wrong with the convolver. Please have a look at https://github.com/litebird/litebird_sim/pull/170/files#diff-389ac96ea9e52593dace164e986e21012662fcfb74deb8a0780e44ea01c0b3f7R75 and let me know if you have any idea what's happening (e.g. that my assumptions are wrong)!

[mreineck]

During the last telecon with the beamconv developers, we established that the expressions for computing the b_lm of the combined detector+HWP system (which I adopted from beamconv) are approximations, and they become worse for small m. It may be that this is the reason for the failing test case.

[mreineck]

I have an idea on which I'd like to hear your opinion:

As far as I know extensive simulations of the beam shapes are done with GRASP, which form the basis of the beam a_lm we use in the simulations. Are there plans to also simulate the beam+HWP system with GRASP? As far as I understand, getting realistic beams for the full system at a small selection of HWP angles (i*pi/5 for i in [0;4]) should be sufficient for an accurate reconstruction at all other angles.

Beam a_lm are small, so this would not be a storage problem, but I expect the GRASP simulations might be expensive.
It would be simple to adjust my code to work with this kind of input.

[ziotom78]

As far as I know extensive simulations of the beam shapes are done with GRASP, which form the basis of the beam a_lm we use in the simulations. Are there plans to also simulate the beam+HWP system with GRASP? As far as I understand, getting realistic beams for the full system at a small selection of HWP angles (i*pi/5 for i in [0;4]) should be sufficient for an accurate reconstruction at all other angles.

Hi Martin, I checked with Cristian Franceschet, who is in charge of running GRASP simulations for LiteBIRD, and he told me that there is no way to simulate an half-wave plate in GRASP. The closest option would be a polarizing grid, which is however not the same, and Cristian is not sure that would be representative enough.

[mreineck]

Oh, that's a pity... I'll continue investigating alternative approaches then.

[marcobortolami]

With the previous commit I fixed the matrix for horizontal linear polarization (there is a 1/2 in the slides here)

[marcobortolami]

I believe that the test was actually working. If you print the convolved TODs they are actually equal. The order of the objects to compare in np.testing.assert_array_less was inverted! Here they say that the first parameter should be "The smaller object to check", while the second parameter should be "The larger object to check". So now np.testing.assert_array_less(np.max(sig) - np.min(sig),1e-4) should be correct as we would like that the difference between the max and the min of the TODs is smaller than the threshold, i.e. 1e-4. Please, tell me if I misunderstood something:)

[marcobortolami]

Yuya Nagano (@okayamayu8) and I tried to compare what is implemented in LBsim, in ducc and in Condor, an external code developed by Yuya Nagano and Yusuke Takase (@yusuke-takase). We tried to do so in 4 different cases:

1. gaussian and symmetric beam, no HWP
2. asymmetric beam, no HWP
3. gaussian and symmetric beam, with HWP
4. asymmetric beam, with HWP

For cases 1 and 3, i.e. gaussian and symmetric beam, we use the following sky and beam:

For cases 2 and 4, i.e. asymmetric beam, we use the following sky and beam:

For cases 1 and 2, i.e. no HWP, we use mueller = np.identity(4) and HWP angles equal to 0.

For cases 3 and 4, i.e. with HWP, we use mueller = np.diag((1.,1.,-1.,-1.)) and random HWP angles.

Below the results (convolved TODs and residuals).

Everything works fine a part for case 4. We are trying to understand why. Yuya tried to convolve in 2 different ways, that are matching between them but not with LBsim. I tried case 4 without the V component, i.e. mueller = np.diag((1.,1.,-1.,0.)), as Condor uses only T,Q,U and the result is the same.

[mreineck]

Thank you very much for these tests! I'm glad that we have good agreement in the simple cases, and I hope that we should be able to find and fix the discrepancies in case 4.

I just did a very quick comparison of case 4 between the Mueller convolver and beamconv, and I see the differences there too. (For a more thorough comparison betwen Mueller convolver and beamconv, see also AdriJD/beamconv#20.)

[marcobortolami]

Hi:) some news about case 4 above, the only one for which we had strange results. We (mostly @okayamayu8) managed to get a nice comparison!
Yuya took the equations of Adri's notes (adri_notes.pdf) and implemented those equations on Condor. He tried running his program and it is consistent with Mueller Convolver (see plots below)!

Possible conclusion: Mueller convolver is actually right and maybe there is a bug on the beamconv side. What do you think?

[okayamayu8]

Hi! I have a question about the situation assumed by this convolver.
If we adopt the unit matrix as the mueller matrix, I think it represents the case where there are no optical elements (w/o HWP).
And to introduce an ideal HWP, I think that we should use MII=1, MQQ=1 and MUU=-1 as the mueller matrix. (I suspect that I am doing something wrong even here).
We then use these mueller matrix for convolution.

input parameters
alm: anything is OK.
blm: symmetric one.
Pointing: rundom (anything is OK)

situation
Case 1: w/o HWP convolution
Case2: w/ HWP but all α=0 convolution (This assumes a stationary HWP at α=0)

I convolve the two cases with the Mueller convolver and my own function (Condor).
As a result, the results of the case1 and case2 match. I got the same result using either tool.
But that seems strange, doesn't it?

Because, by solving the Müller matrix to find the respective tod,
Case 1: d(Ω) = I(Ω) + U(Ω)cos(2ψ) -U(Ω)sin(2ψ)
Case2: d(Ω) = I(Ω) + U(Ω)cos(4α - 2ψ) -U(Ω)sin(4α - 2ψ)
so it is not likely to match at α=0.

[mreineck]

Hi @okayamayu8, your setup certainly looks correct ... I don't immediately know what is going wrong there. Do you have a snippet of code that runs case 1 and case 2 with the Mueller convolver, so that I can have a closer look?

[okayamayu8]

`import ducc0
import numpy as np
import healpy as hp
import mueller_convolver as mc
from litebird_sim import MuellerConvolver
import matplotlib.pyplot as plt
figsize = (10,8)
rng = np.random.default_rng(42)

parameter settings

fct = 1.0 #some factor, I don't now
lmax = 32
mmax = 5
ncomp = 3 #put 1 for T only, 3 for T and P
nptg = 100 #number of pointings to produce
fwhm = 10.0 #deg
nside = 64

epsilon = 1e-4 #desired accuracy for the interpolation; a typical value is 1e-4
ofactor = 1.5 #oversampling factor to be used for the interpolation grids. Should be in the range [1.2; 2], a typical value is 1.5. Increasing this factor makes (adjoint) convolution slower and increases memory consumption, but speeds up interpolation/deinterpolation.
nthreads = 0

alm_size = hp.Alm.getsize(lmax, mmax=mmax)

preparation rundom alm

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

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

#beam alms
blm_gauss = hp.blm_gauss(np.deg2rad(fwhm), lmax=lmax, pol=(ncomp==3)) # generate a Gaussian beam using healpy #random_alm(lmax, kmax, ncomp, rng)
print("blm_gauss.shape: ",blm_gauss.shape)
blm = np.zeros((ncomp,alm_size),dtype=np.complex128)
print("blm.shape: ",blm.shape)
for l in range(lmax+1):
for m in range(ncomp):
index = hp.Alm.getidx(lmax,l,m)
blm[0,index] = blm_gauss[0][index]
if(ncomp==3):
blm[1,index] = blm_gauss[1][index]#np.exp(4jnp.pi/3)
blm[2,index] = blm_gauss[2][index]#np.exp(4jnp.pi/3)

prepare mueller

identity = np.identity(4) * fct
hwp = np.zeros((4,4))
hwp[0,0] = 1.
hwp[1,1] = 1.
hwp[2,2] = -1.
hwp[3,3] = 1.

make rundom pointing

ptg = np.empty((nptg, 3))
ptg[:, 0] = rng.random(nptg) * np.pi
ptg[:, 1] = rng.random(nptg) * 2 * np.pi
ptg[:, 2] = rng.random(nptg) * 2 * np.pi
alpha = rng.random(nptg) * 2 * np.pi

to see result clearly I = 0

slm[0]=0

fullconv = MuellerConvolver(
lmax,
mmax,
slm,
blm,
hwp,
single_precision=False,
epsilon=epsilon,
ofactor=ofactor,
nthreads=nthreads,
)

fullconv_identity = MuellerConvolver(
lmax,
mmax,
slm,
blm,
identity,
single_precision=False,
epsilon=epsilon,
ofactor=ofactor,
nthreads=nthreads,
)

sig = fullconv.signal(ptg, alpha*0)
sig_identity = fullconv_identity.signal(ptg, alpha)

fig = plt.figure(figsize=figsize)
plt.plot(sig, label = "w/ HWP alpha=0")
plt.plot(sig_identity, label = "w/o HWP")
plt.legend()
`

This is the calculation program I used.
Please let me know if I'm using it wrong or something.

[okayamayu8]

If I think about it, this implementation is based on the assumption of HWP, which limits the coefficients of spherical harmonic expansions such as CPP*.
So I would like to understand that w/o HWP cannot be calculated in any way.

[mreineck]

I think that the implementation should produce meaningful results even for identity Mueller matrices.

I'm no longer sure that this is actually a bug ...
If I understand correcly, we have a Gaussian beam that is sensitive to Q polarisation (in its own frame), and in front of it we have an optical element that flips the sign of the U component. So, isn't it expected that the result remains unchanged?

[okayamayu8]

Sorry, in the discussion, I seem to have confused the detector angle with the boresight angle.
Correctly, the two must be distinguished, and I will write below detector angles as xi.
Then the simple scan equations for ① without HWP and ② with HWP are as follows.

① I(Ω) +Q(Ω) cos(2ψ+2xi) - U(Ω) sin(2ψ+2xi)
② I(Ω) +Q(Ω) cos(4α+2ψ-2xi) - U(Ω) sin(4α+2ψ-2xi)

So, if we consider cases such as α=0, xi=0, then ① and ② would have the same result.

The calculation I was doing represents exactly such a case, so I do not think it is a bug in convolver.