Add 4pi beam convolution

[mreineck]

This aims to add 4pi convolution with non-axisymmetric beams (very similar to the "conviqt" functionality in TOAST, but more accurate and hopefully also quicker).
Currently this cannot be fully implemented because the code requires

• a_lm of the sky model maps (probably solved by Alms in Mbs #306, I need to check how I can access them).
• a_lm of the beams. Those will be in the IMO at some point, but have not yet been integrated.

It should be possible to load the beams from files for testing purposes, using healpy.read_alm, but I'm not sure how to provide the corresponding file names in the input parameters and how they would be passed to the beam convolution module.
Unfortunately I'm not very familiar with the simulation framework and the really large set of tools it uses. I'm happy to help tweaking and debugging the underlying algorithms, but I need help from experts with the high-level integration.

[mreineck]

I must say that some of the pre-commit hooks are driving me quite crazy ...
I attempt to add the following lines to litebid_sim/__init__.py:

from .beam_convolution import (
    add_convolved_sky_to_observations,
)

This is what I'm trying to do:

(lbsim) martin@marvin:~/codes/litebird_sim$ git status
On branch beam_comvolution_2
Your branch is up to date with 'origin/beam_comvolution_2'.

Changes not staged for commit:
  (use "git add <file>..." to update what will be committed)
  (use "git restore <file>..." to discard changes in working directory)
	modified:   litebird_sim/__init__.py

no changes added to commit (use "git add" and/or "git commit -a")
(lbsim) martin@marvin:~/codes/litebird_sim$ git diff
diff --git a/litebird_sim/__init__.py b/litebird_sim/__init__.py
index b4a8aa6..1985a0c 100644
--- a/litebird_sim/__init__.py
+++ b/litebird_sim/__init__.py
@@ -150,6 +150,9 @@ from .spacecraft import (
     SpacecraftOrbit,
     SpacecraftPositionAndVelocity,
 )
+from .beam_convolution import (
+    add_convolved_sky_to_observations,
+)
 from .version import __author__, __version__
 
 # Check if the TOAST2 mapmaker is available

So I have the mentioned changes in my tree.
Now I run

(lbsim) martin@marvin:~/codes/litebird_sim$ git commit -a -m "add method to __init__.py"
ruff.....................................................................Failed
- hook id: ruff
- files were modified by this hook

Found 1 error (1 fixed, 0 remaining).

ruff-format..............................................................Passed

If I now look at the status, all changes are gone!

(lbsim) martin@marvin:~/codes/litebird_sim$ git status
On branch beam_comvolution_2
Your branch is up to date with 'origin/beam_comvolution_2'.

nothing to commit, working tree clean

As you can see, the pre-commit hooks somehow remove my changes from __init__.py! What is going wrong here?

[paganol]

Hi @mreineck, I can help in the high level integration. I can work on this branch and prepare a coherent lbs interface, what do you think?
Regarding the pre-commit, I cannot help, I never had this kind of problem

[mreineck]

Thank you so much, @paganol, this is a great relief for me!

Concerning the precommit hook, I managed to get the crucial information of of ruff by ading -v:

(lbsim) martin@marvin:~/codes/litebird_sim$ ruff litebird_sim/__init__.py 
litebird_sim/__init__.py:22:5: F401 [*] `.beam_convolution.add_convolved_sky_to_observations` imported but unused
Found 1 error.
[*] 1 fixable with the `--fix` option.

So I need to add the new function to __all__ further down in the file. Makes a lot of sense, but I really don't feel that silently removing the import was helpful to understand the problem. Maybe we shouldn't use --fix and just let ruff complain about things (and explain them), until everything has been fixed manually.

[mreineck]

I wonder whether we should merge #337 wit this PR, then we can add the signal modulation by the HWP as well (optionally, of course!). The only additional required quantity for this is the time-dependent HWP angle.

[mreineck]

I know that the PR doesn't look like much, but apart from obtaining the missing input data (which, I admit, may be complicated to achieve!), it is pretty much complete. All the numerical details are handled by ducc.

[mreineck]

Thank you, @paganol, your interfacing work looks great to me!
I see that some error messages still refer to scan_map_in_observations; I'm happy to update them, unless there is some deeper connection to scan_map_in_observations that I don't see at the moment.
I'll also copy over mueller_convolver.py from #337 and close #337 afterwards; having these two changes in a single PR makes much more sense!

[paganol]

Hi @mreineck, sorry I forgot to update the error messages! No deeper connection 😅... Of course if we find better names for the lower level functions we can change them. About merging the PRs, I totally agree with you.

[mreineck]

Stupid question: is there a place in litebird_sim where I can put mueller_convolver.py, so that it is not reformatted by ruff? Strictly speaking this piece of code is maintained somewhere else (although it doesn't have an "official" place yet), and reformatting (and especially the re-naming of variables that ruff requests) will be quite a maintenance headache.
Should we have an external/ subdirectory?

[ziotom78]

Stupid question: is there a place in litebird_sim where I can put mueller_convolver.py, so that it is not reformatted by ruff? Strictly speaking this piece of code is maintained somewhere else (although it doesn't have an "official" place yet), and reformatting (and especially the re-naming of variables that ruff requests) will be quite a maintenance headache. Should we have an external/ subdirectory?

Would #fmt: off at the beginning of the file be ok, or do you want the file to be left completely untouched?

[mreineck]

No, that's perfectly fine, thanks for pointing this out!

[paganol]

Hi @mreineck, for matching the latest changes in the pointing API (PR #340) this line

litebird_sim/litebird_sim/beam_convolution.py

Line 65 in b671b68

curr_pointings_det = curr_pointings_det.reshape(-1, 3)

should be removed.

[mreineck]

Thanks for te tip, I'll fix this!
I'm travelling right now, but I hope to get back to this PR next week.

[mreineck]

If this has the blessing of @martamonelli, I'm totally happy with it!
@martamonelli, does this also resolve (some of) the remaining discrepancies with beamconv?

[martamonelli]

Hi all,
the equations (and their derivations) shared by @okayamayu8 make perfect sense to me. Together with @paganol, we came up with a few tests (and minor changes to lbs) that will help us make sure the convolution is doing what it's supposed to. Let's see.

@mreineck: I started working on checking what changes for the comparison with beamconv, but I'm not done yet. I'll let you know as soon as I'm done.

[mreineck]

Sorry for the long delay, @martamonelli, and thanks a lot also for te verification work! Having a tested Mueller convolver is a huge milestone!

[martamonelli]

@mreineck, I've checked whether the changes proposed by Yuya to Mueller convolver take care of the discrepancies with beamconv and they don't. (This does make sense though: Yuya's formulae for the beam convolution assumes the same conventions as scan map, which seem to be different than the ones used in beamconv.)

[mreineck]

Thank you for looking into this!
I'm glad it is a convention issue (although that doesn't make things any less complicated :-) and will port the changes back to my master version of Mueller convolver with additional comments to clarify this.

[paganol]

Good news! I made some tests, after having merged master, and having adjusted the interface. The results are consistent with the previous implementation if the beams are rotated by the proper polarization angle.
This simple piece of code, for a Gaussian beam, is sufficient:

def make_gauss_beam(fwhm_rad,lmax,mmax,psi_pol_rad):
    blm = lbs.gauss_beam_to_alm(lmax,
                                mmax,
                                fwhm_min_rad=fwhm_rad,
                                fwhm_max_rad=fwhm_rad,
                                psi_ell_rad=0,
                                psi_pol_rad=psi_pol_rad,
                                cross_polar_leakage=0)
    return blm

dets = []
blms = {}
for n_det in detlist:
    det = lbs.DetectorInfo.from_imo(
        url=f"/releases/{imo_version}/satellite/{telescope}/{channel}/{n_det}/detector_info",
        imo=imo,
    )
    dets.append(det)
    blms[det.name] = make_gauss_beam(fwhm_rad,lmax,mmax,psi_pol_rad=det.pol_angle_rad)

btw.. potentially we could add a function like this to the code that automatically:

• reads fwhm, ellipticity, angle of the ellipses and polarization angle
• synthesize the elliptical gaussian beam for each detector
• store them in a dictionary.

But we can discuss this later..

Regarding the comparison, here the results for Gaussian beams:
this with scan_map assuming ideal HWP

this for non ideal HWP, assuming for each detector an ideal HWP + a random perturbation

I think we accomplished the first part of the validation. @martamonelli, @okayamayu8, @mreineck, @ziotom78, @patanch, I guess we can move to the asymmetric beam validation and also work a bit on the optimization.

Last thing.. the ~0.2% normalization mismatch in polarization is still there.. whatever I do.. any idea?

[mreineck]

Just curious, do you have an intuition why there is so much ringing in the polarization maps (and not in T) for nonideal HWP?

[paganol]

Hi @mreineck, I also thought it was ringing at first, but I believe it's actually T-to-P leakage caused by the coupling of a non-ideal HWP and the scanning strategy. When I feed (0, Q, U) into the simulation, this is the result I get:

Probably the ringing-like pattern is due to the scanning strategy.. does that sound convincing to you?

[mreineck]

I did think of a purely scanning strategy related issue, i.e. that the equation system that needs to be solved for Q/U was ill-conditioned in some pixels. But that would have affected the ideal HWP as well, I guess. The leakage theory sounds much more convincing to me, thanks!

[mreineck]

I still have absolutely no clue about the normalization issue; this almost feels as if there was a pi=3.14 somewhere in the code ... or a division by Npix-1 instead of Npix.

[paganol]

I still have absolutely no clue about the normalization issue; this almost feels as if there was a pi=3.14 somewhere in the code ... or a division by Npix-1 instead of Npix.

Yeah, I was checking for that too, but I couldn’t find anything.