Beam convolution

[giuspugl]

Hello,
as discussed am posting the links related to the current implementation of the Convolution Operators in toast3,

https://github.com/hpc4cmb/toast/blob/c964a2586575e88144eb80206ca869d2ba012ab6/src/toast/ops/conviqt.py#L696

and the test routines used to validate the implementation
https://github.com/hpc4cmb/toast/blob/c964a2586575e88144eb80206ca869d2ba012ab6/src/toast/tests/ops_sim_tod_conviqt.py#L460

[mreineck]

After way too much procrastinating (sorry!) I'll try to pick this up.

To make sure I understand things right: the task is to get mueller_convolver integrated into the litebird_sim API, correct?

If that is the case, our best starting point is probably @okayamayu8's work, who has been running HWP simulations with a modified version of litebird_sim recently.

@okayamayu8, would it be OK for you to share your code for future integration, perhaps as part of #170?

[okayamayu8]

Sure.
As for the code I shared, you can use it as you like.

[mreineck]

@okayamayu8 , thank you very much!

I'm sorry that I'm probably missing something, but in the Julia notebook you shared I don't actually see any direct usage of litebird_sim or mueller_convolver. My impression was that you were probably working on a Python code that integrated mueller_convolver with litebird_sim and ran on supercomputer nodes. If that is indeed the case, it would be great to use this for this issue. If not, we probably have to implement this from scratch ... unfortunately I don't really feel qualified for this :-/

[okayamayu8]

Sorry, I didn't understand exactly what you wanted to do.
That note uses random pointing as a validation of the mueller convolver.
So it is not working with litebird_sim.
And for mapmake, which we have been discussing recently, we were using a combination of litebird_sim and other external functions (written in julia).
So I had escaped from Python code.
So I am aware that we need to create from scratch a script that is complete with only litebird_sim and mueller convolver.

[mreineck]

The misunderstanding was on my side, sorry! I had implicity assumed that your mapmaking work and the mueller_convolver tests were parts of the same code.

@ziotom78, I'm a bit at a loss now how to proceed. My own knowledge of the litebird_sim data structures is not sufficient to come up with a solid integration; at the same time I know the mueller_convolver details and intricacies pretty well :) Can you suggest someone I could team up with? My estimate is that it would take maybe two days of dedicated work for the two of us...

[ziotom78]

Hi @mreineck , looking at the TOAST code, it seems that you need pointings, TODs, and synthetic maps; they are described at the links https://litebird-sim.readthedocs.io/en/latest/observations.html and https://litebird-sim.readthedocs.io/en/latest/sky_maps.html. Here are a few notes:

• Maps are created by the Mbs module, which is described here: https://litebird-sim.readthedocs.io/en/latest/sky_maps.html. They are usually stored in a dictionary, where the key is the detector name

• TODs are stored in matrices with shape (n_d, N), where n_d is the number of detectors in the Observation object and N is the number of samples in the timeline. TODs can have arbitrary names, so it would be wise to let your function to accept any name for the TOD:

  def add_convolved_sky_to_observations(
      obs_list: List[Observation],
      …,
      map_dictionary: Dict[str, Any],  # Maps to be convolved
      det_to_map_names: Dict[str, str] | None = None, # Correspondence btw detectors and maps
      component: str = "tod",
  ):
      ptg_list = [ob.pointings for ob in obs]  # θ and φ
      psi_list = [ob.psi for ob in obs] # ψ (polarization angle)
  
      if not det_to_map_names:
          # If `det_to_map_names` was not provided, let's assume that for each
          # detector in `obs_list` (like "000_000_003_QA_040_T") there is a key
          # with the same name in `map_dictionary` containing the I/Q/U maps
          # for that detector
          det_to_map_names = [cur_obs.name for cur_obs in obs_list]
  
      for cur_obs, cur_ptg, cur_psi, cur_det_to_map_names in zip(obs_list, ptg_list, psi_list, det_to_map_names):
          cur_tod = getattr(cur_obs, component)
          # Now you're sure that `cur_tod` is a (n_d, N) matrix and `cur_ptg` is a (n_d, N, 2) matrix
          # containing the θ, φ angles (in radians) of the pointings for each detector
  
          for cur_idx, cur_det_name in enumerate(det_to_map_names):
              cur_det_map = map_dictionary[cur_det_name]  # Here are the three I/Q/U maps
  
              # Implement the convolution and add the convolved TOD to cur_tod[cur_idx, :],
              # using the pointings θ=cur_ptg[cur_idx, :, 0], φ=cur_ptg[cur_idx, :, 1], and
              # ψ=cur_psi[cur_idx, :] (all in radians)
              …

  Perhaps there is some smarter solution where you swap the two loops (the outer and inner fors) so that you can cache the datacube for the convolved maps instead of generating it again for every cur_obs object.

• Since what is going to be implemented is conceptually similar to the dipole.py component, where you calculate an additional component and sum it into the timelines, you can have a look at litebird_sim/dipole.py.

Here is a sketch of a main.py program to test the convolution:

import numpy as np
import astropy
import litebird_sim as lbs

sim = lbs.Simulation(
    base_path="./example",
    start_time=astropy.time.Time(
        "2030-01-01T00:00:00",
    ),
    duration_s=86_400.0,
    name="My simulation",
    description="My description",
    random_seed=12345,
    imo=lbs.Imo(
        flatfile_location=lbs.PTEP_IMO_LOCATION,
    ),
)

lft_instrument = lbs.InstrumentInfo.from_imo(
    imo=sim.imo, url="/releases/vPTEP/satellite/LFT/instrument_info"
)
sim.set_instrument(lft_instrument)

# Get information about two 40 GHz detectors
det1 = lbs.DetectorInfo.from_imo(
    imo=sim.imo,
    url="/releases/vPTEP/satellite/LFT/L1-040/000_000_003_QA_040_T/detector_info",
)

det2 = lbs.DetectorInfo.from_imo(
    imo=sim.imo,
    url="/releases/vPTEP/satellite/LFT/L1-040/000_000_003_QA_040_B/detector_info",
)

sim.create_observations(
    detectors=[det1, det2],
    tods=[
        # If you want to calculate both the unconvolved and convolved components
        # and keep them separate, here is how you do it:
        lbs.TodDescription(name="unconvolved_sky"),
        lbs.TodDescription(name="convolved_sky"),
    ],
)
# From now on, you have the `Observation` objects in sim.observations,
# and each of them has the fields obs.unconvolved_sky and obs.convolved_sky
# (two 2D matrices with shape n_d, N).

# Let's compute some sky maps using PySM
params = lbs.MbsParameters(
    nside=256,
    make_cmb=True,
    make_fg=True,
    fg_models=["pysm_synch_0", "pysm_freefree_1"],
)
mbs = lbs.Mbs(
    simulation=sim,
    parameters=params,
    detector_list=[det1, det2],
)
input_maps = mbs.run_all()[0]

# And now let's generate the pointings
sim.set_scanning_strategy(
    lbs.SpinningScanningStrategy.from_imo(
        imo=sim.imo, url="/releases/vPTEP/satellite/scanning_parameters"
    )
)

sim.set_hwp(
    lbs.IdealHWP(ang_speed_radpsec=1.0),
)

sim.compute_pointings()
# From now on, each `Observation` object in sim.observations contains
# the two fields `pointings` and `psi`: the former is a (n_d, N, 2) matrix
# containing θ and φ for each detector and sample, while the latter
# contains the polarization angles.

# Let's calculate the component of the unconvolved sky for reference
lbs.scan_map_in_observations(
    sim.observations[0],
    input_maps,
    input_map_in_galactic=False,
    component="unconvolved_sky",  # Store this in `unconvolved_sky`
    interpolation="",
)

# Here is the call to the 4π convolution code! Debugging should happen
# inside this function call
add_convolved_sky_to_observations(
    obs_list=sim.observations,
    map_dictionary=input_maps,
    component="convolved_sky",
)

# Convolution was completed, let's now make a map (not sure if using
# destriping is meaningful here, but the binned map is computed in any case).
result = lbs.make_destriped_map(
    nside=256,
    obs=sim.observations,
    params=lbs.DestriperParameters(),
    components=["convolved_sky"],
)
sim.write_healpix_map("hit_map.fits", result.hit_map)
sim.write_healpix_map("binned_map.fits", result.binned_map)
sim.write_healpix_map("destriped_map.fits", result.destriped_map)

sim.flush()

Perhaps @giuspugl can help you in the implementation. I'm available too: if you need help, let's schedule a Zoom call!

[paganol]

Hi there, I add a couple of comments to what @ziotom78 wrote.

• Keep in mind that sim.set_hwp followed by sim.compute_pointings automatically adds twice the hwp angle to the polarization angle, since MuellerConvolver.signal takes pointing and the hwp angle separatelly, we might think to pass a hwp object to signal
• For handling the maps produced by Mbs I think that the example of lbs.scan_map_in_observations can be very helpful

@mreineck, I available for a zoom call as well!

[mreineck]

@ziotom78, @paganol, thank you very much for all the input, this is really valuable! I'll be attending a workshop next week, so might not be able to make much progress, but I'll definitely keep you updated! This feels much more manageable now!

[mreineck]

OK, I'm finally having some time to proceed!

Small comments to the discussion so far:

• the sky to be convolved should not be provided in form of maps, but rather a_lm. Of course I could run any input maps through map2alm, but then there is the question of the band limits, number of iterations etc. that should be used, and in some cases the input even comes naturally in a_lm form, for example the CMB component.
  As far as I see there is no specific a_lm data type in litebird_sim yet, which makes things interesting...
• I also need the description of the beam in a_lm form. I suppose/hope that this will be part of the IMO at some point, but currently it doesn't appear to be. Doing general beam convolution using the Gaussian approximations currently described in the IMO is probably pointless.
  (NOTE: So far I looked into the "public" IMO only. The official IMO page didn't take any of the login data I had to offer; I'll try to get access to that.)

I'm confident that I can make a skeleton code which invokes the necessary convolution functions, but the first iteration will have quite a few holes in it.

[paganol]

Hi @mreineck. A quick comment regarding the first point.
In the current implementation the input maps are usually produced by Mbs, and then passed to scan_map. We could easily add a feature in Mbs for outputting alms if requested.

[mreineck]

Keep in mind that sim.set_hwp followed by sim.compute_pointings automatically adds twice the hwp angle to the polarization angle

I'm not sure I understand the purpose of this, @paganol. In the presence of non-idealities (either non-cicular beams or non-ideal HWP) these two angles cannot be combined to a single value without losing essential information...

I also wonder whether calling psi the "polarization angle" is a choice that describes the quantity well ... even asuming a completely unpolarized sky and detector, this angle is important, if the the beam is not axisymmetric.

[mreineck]

In the current implementation the input maps are usually produced by Mbs, and then passed to scan_map. We could easily add a feature in Mbs for outputting alms if requested.

Thanks for pointing this out! In fact that would be a great feature, which will potentially both save time and improve accuracy.

[mreineck]

Addition to the above: in the Planck simulation pipeline, sky model component "maps" would always be generated in a_lm form first, with a very high band limit. Whenever a real space map was required, the component a_lm would be combined with the appropriate weights and then transformed to the maps at the requested resolution (using a band limit depending on the requested Nside).

[mreineck]

Naive question: is there a reason why Observation.__attr_det_names is not accessible (read-only) from outside? Having a correspondence between detector names and indices in the observation's data arrays could be very convenient.

[ziotom78]

Naive question: is there a reason why Observation.__attr_det_names is not accessible (read-only) from outside? Having a correspondence between detector names and indices in the observation's data arrays could be very convenient.

There is this nice feature implemented by @dpole that any field of the DetectorInfo class associated with a detector is “spread” over each Observation object. Thus, to know the name of the detectors associated with Observation you can use the field name (the same as in DetectorInfo.name):

import litebird_sim as lbs
import astropy

sim = lbs.Simulation(
    base_path="./example",
    start_time=astropy.time.Time(
        "2030-01-01T00:00:00",
    ),
    duration_s=86_400.0,
    name="My simulation",
    description="My description",
    random_seed=12345,
    imo=lbs.Imo(
        flatfile_location=lbs.PTEP_IMO_LOCATION,
    ),
)

lft_instrument = lbs.InstrumentInfo.from_imo(
    imo=sim.imo, url="/releases/vPTEP/satellite/LFT/instrument_info"
)
sim.set_instrument(lft_instrument)

det1 = lbs.DetectorInfo.from_imo(
    imo=sim.imo,
    url="/releases/vPTEP/satellite/LFT/L1-040/" "000_000_003_QA_040_T/detector_info",
)
det2 = lbs.DetectorInfo.from_imo(
    imo=sim.imo,
    url="/releases/vPTEP/satellite/LFT/L1-040/" "000_000_003_QA_040_B/detector_info",
)

import numpy as np

sim.create_observations(
    detectors=[det1, det2],
    tods=[
        lbs.TodDescription(
            name="tod",
            description="TOD",
            dtype=np.float64,
        ),
        lbs.TodDescription(
            name="noise", description="1/f+white noise", dtype=np.float32
        ),
    ],
)

print(sim.observations[0].name)
# Output:
# ['000_000_003_QA_040_T' '000_000_003_QA_040_B']