New API for pointings

[ziotom78]

I am starting a discussion to decide how to rework the code used to compute the pointing information and the angle of the HWP. First, I'll present the structure of the current code, and then I will present a few ideas to improve the performance.

Basic principles

The code computes the final pointing of each detector according to the transformations explained in the chapter Scanning Strategy of the User's Manual. Each rotation is encoded using a quaternion, and there are several of them:

1. One quaternion to convert from the reference frame of a detector to the reference frame of the boresight of the instrument (LFT, MFT, HFT); this quaternion at the moment is considered to be constant in time.

2. One quaternion to convert from the reference frame of the boresight of the instrument to the reference frame of the spacecraft (sometimes called “spin reference frame”, because the z axis is aligned with the spin axis). This quaternion is constant in time as well.

3. One quaternion to convert from the reference frame of the spacecraft to Ecliptic coordinates. This quaternion does of course depend on time, as it encodes the constant precession of the spin axis around the Sun-Earth axis and the revolution of the Earth around the Sun.

The code multiplies the quaternions to get the overall transformation from the reference frame of a detector to Ecliptic coordinates; these quaternions are sampled with a customizable sample rate, which is usually chosen to be a fraction of 1 Hz because the typical angular speeds are slower than that.

How the code works now

• The full quaternions are combined generated by Simulation.set_scanning_strategy. These are sampled at a lower frequency than the scientific TOD—typically ~1 min—and are saved in the field Simulation.spin2ecliptic_quats, which is a (N_q, 4) matrix, where N_q is the number of quaternions. As the quaternions are related to the spacecraft, one array is enough for a full focal-plane simulation. These quaternions encode the details of the precession of the spin axis around the Sun-Earth axis and model the revolution of the Earth (and thus of the L₂ point where LiteBIRD will stay) around the Sun.

• The real juice is in get_pointings, which works on a single Observation. It takes three quaternions:

  1. The detector→boresight quaternion (just one per detector), grabbed from the Observation object (this is the only field read by the function from this object);
  2. The boresight→spin-axis quaternion (just one);
  3. The spin→ecliptic coordinates (many; they have been calculated by Simulation.set_scanning_strategy, see above).

  The function calls get_det2ecl_quaternions, which employs ducc0 to compute all the quaternions at the nominal sampling frequency of the TOD. There is just one buffer with shape (N, 4) that gets allocated at the beginning, and then when looping over the detectors it keeps overwriting it every time.

  Once the fully sampled quaternions have been computed, the function all_compute_pointing_and_orientation converts them into actual pointings, i.e., θ, φ, and ψ. These pointings are saved by get_pointings into Observation.

• The method Simulation.compute_pointings iterates over a set of Observations to call get_pointings for each of them. The purpose of this function is to save the pointings within each Observation object, but it also returns them to the caller. (This is a reference to the pointings in the Observation objects, so no actual copies of the arrays are made.)

Pointings are used in four modules:

1. scan_map.py;

2. dipole.py;

3. The binning map-maker;

4. The destriper.

What we want from the code

The actual API has served us well, but it has one significant limitation: all the pointings are stored in memory, and this takes more space than the scientific signal itself! (We need just one floating-point for each scientific sample but 4 floating-point numbers for the pointings: θ, φ, ψ, and the HWP angle α, which however is the same for all the detectors.) Moreover, there are a few other limitations:

1. The HWP angles are added to the orientation angles. However, we would like to keep them separate, because this is needed to correctly perform 4π beam convolution.

2. We want to simulate time-dependent effects, in particular:

   • Effects that alter the orientation of the detectors with respect to the boresight of the focal plane;
   • Effects that alter the orientation of the boresight of an instrument with respect to the rest frame of the spacecraft;
   • Non-idealities in the HWP that can affect the pointings. These effects should not be modeled through lowly-sampled quaternions, as the HWP spins much more rapidly than the spacecraft.

A proposal for a new API

The idea we want to discuss here is a new API which can help to save memory. Since all the code that uses pointings so far requires to access pointing information for one detector at a time, the idea is to add a method that returns the pointings for just one detector, instead of keeping the pointings of all the detectors in one huge matrix.

In the following code snippets, I will always include typing hints to make the code clearer. I'll assume that the following imports are valid:

from typing import Optional
import numpy as np
import numpy.typing as npt

How to handle time-dependent effects?

At the moment, there are three families of quaternions:

1. Detector→boresight quaternions, which are assumed to be constant in time;
2. Boresight→spin quaternions, which are constant;
3. Spin→Eliptic quaternions, which are time-dependent.

To ease the burden to interpolate time-dependent quaternions, the code implements the class Spin2EclipticQuaternions that keeps an array of quaternions, each associated with a time, and can interpolate them at any arbitrary sampling frequency. Of course, this class is used only with the quaternions in point 3.

We should make Spin2EclipticQuaternions more general by renaming it to TimeDependentQuaternions and add the ability to multiply two TimeDependentQuaternions instances, returning a new TimeDependentQuaternions. In this way we could make the quaternions in point 1. and 2. above time-dependent and combine the three families as easy as we are doing now. This would enable us to simulate a quite large amount of time-dependent optical systematics caused either by mechanical or thermal fluctuations.

It is conceptually easy to implement the multiplication of two TimeDependentQuaternions objects, but special care must be applied to times: the code must ensure that both quaternions are sampled at the same frequency!

How to handle systematics caused by the HWP?

We already have a HWP class, from which IdealHWP is derived. We should change the API so that the class has a virtual method named apply_hwp_to_pointings with the following signature:

class HWP:
    …

    def apply_hwp_to_pointings(
        self,
        start_time_s: float,
        delta_time_s: float,  # This will be 1/(19 Hz) in most cases
        bore2ecl_quaternions: npt.ArrayLike,  # Boresight→Ecliptic quaternions at 19 Hz
        hwp_angle: npt.ArrayLike,
    ) -> None:
        …

An ideal HWP would just use hwp_angle, storing the value 2πνt into it (with ν being the angular speed of the HWP), and ignore the presence of bore2ecl_quaternions. However, in a derived class there might be the simulation of some systematic effects that introduces some noise in the pointings; in the latter case, the code is free to modify quaternions by multiplying each of them by a new quaternion on the left.

Note that we do not want to use TimeDependentQuaternions here, as that class is optimized for low-sampling quaternions, while a HWP spins quickly! Instead, apply_hwp_to_pointings gets the fully sampled quaternions.

How to properly query pointings for one detector

The amount of information needed to be carried over in the code is larger if we want to recompute all the pointings whenever we need them. (The current status is that the full pointing information is stored in the (N, 3) pointing matrix containing θ, φ, and ψ, and nothing else is needed.) We must always be able to know (1) if there are time-dependent effect to be applied to the detector/boresight orientation, (2) if there are non-idealities in the HWP, (3) which coordinate system are we using.

The cleanest way to do this is to keep all of these information in a new class, which we might call PointingProvider. This class should keep an instance of a HWP class (or any derived class), an instance of TimeDependentQuaternions that enables to convert from the reference frame of the boresight to the reference frame of the Ecliptic/Galactic coordinate system, and the specification of the coordinate system itself.

A copy of PointingProvider could be stored in each Observation object, as having pointings immediately available there has proved to be useful so far. (The fact that NumPy are reference-counted grants that the amount of memory used by the code won't explode.)

The class should implement the following methods:

class PointingProvider:
    def __init__(
        self,
        # Note that we require here *boresight*→Ecliptic instead of *spin*→Ecliptic
        bore2ecliptic_quats: TimeDependentQuaternions,
        hwp: Optional[HWP],
    ):
        …

    def get_pointings(
        self,
        # These can now be time-dependent
        detector_quat: TimeDependentQuaternions,
        # The times at which the pointings must be computed (typically @ 19 Hz)
        times: npt.ArrayLike,
        # If None, a new NumPy array will be allocated
        pointing_buffer: Optional[npt.ArrayLike] = None,
        # Ditto
        psi_buffer: Optional[npt.ArrayLike] = None,
        # If no HWP was specified in __init__ and here we have None,
        # no allocation will be done
        hwp_buffer: Optional[npt.ArrayLike] = None,
    ) -> (np.ndarray, np.ndarray):
        # Compute the pointings following this algorithm:
        # - Resample `self.bore2ecliptic_quats` to 19 Hz
        # - Call self.hwp.apply_hwp_to_pointings to “hack” the quaternions and get the HWP angle
        # - Convert the quaternions into actual pointings

        …

The problem with this design is that PointingProvider.get_pointings() requires to allocate a NumPy array with shape (N, 4) for the N quaternions at 19 Hz. A possible optimization could be to cache this array, as it is very likely that multiple calls to get_pointings will need the same number of samples (it will be typically iterated over several detectors). Here is a possible implementation:

class PointingProvider:
    def __init__(
        …
    ):
        …

        # This will be the buffer used internally to keep the 19 Hz quaternions
        self.quaternion_buffer = None

    def get_pointings(
        …
    ):
        …

        # Figure out the expected shape of the quaternion buffer, which should be (N, 4)
        expected_shape = (…, 4)
        if (self.quaternion_buffer is None) or (self.quaternion_buffer.shape != expected_shape):
            self.quaternion_buffer = np.empty(expected_shape)

        # Ok, now use `self.quaternion_buffer` to store the fully-sampled quaternions

How PointingProvider should be used by callers

Assuming that each object Observation have a member variable called pointing_provider of type PointingProvider, one could get the pointings using the following syntax:

# We ask the function to allocate the memory needed for the output arrays
(pointings, psi, hwp_angle) = obs.pointing_provider.get_pointings(
    detector_quat=obs.detector_quat[0],
    times=obs.get_time(),
)

Of course, we could build a tiny method that wraps this call inside the Observation class:

class Observation:
    # …

    def get_pointings(self, detector_idx: int):
        return self.pointing_provider.get_pointings(
            detector_quat=self.detector_quat[detector_idx],
            times=self.get_time(),
        )

In this way, the previous snipped would become

(pointings, psi, hwp_angle) = obs.get_pointings(0)  # First detector

[mreineck]

I like the design a lot, thanks!

The only thing I'm not sure about is having the the internal self.quaternion_buffer ... this will take memory management out of the hands of the users, and they cannot take it back. My feeling is that it could be enough to allow optional specification of output arrays (pointing_buffer, psi_buffer, and hwp_buffer) as you already do.

[ziotom78]

I see your point. We might keep this approach but add a couple more methods to PointingProvider:

1. A method that returns the expected shape of the quaternions, so that the caller can pre-allocate the buffer on their own;
2. A method that sets the quaternion buffer (e.g., PointingProvider.set_quaternion_buffer()).

In this way, if one does not care about how the buffer is allocated will get the optimal solution, because the buffer will be allocated just once by PointingProvider.get_pointings() (using the same code above), while more advanced users can pre-allocate the buffer with PointingProvider.set_quaternion_buffer() and then let PointingProvider.get_pointings() use it transparently.

[fpiacentini]

Dear Maurizio,
thanks for this proposal.

The MHFT instrument team is looking into an option with a reflecting HWP. The reflection on the HWP may induce a time variation of the detector-boresight quaternion, synchronous with the HWP rotation, if the HWP is tilted with respect to its rotation axis.
A similar effect can be induced in the case of a transmissive HWP, if the HWP surfaces are not parallel to each other (wedge).

Other long-term effects, such as thermo-elastic deformations, may induce a slow evolution of the detector-boresight quaternion.

With your new proposal, these effects can be included using a time variable detector-boresight quaternion. The only issue, is that in case of an HWP synchronous effect, the sampling of the pointing must be as fast as the sampling of the detectors.

Alessandro Novelli and Silvia Micheli have simulated the HWP synchronous variation with a workaround, that is currently quite slow to run. To my knowledge, this effect is already included in falcons.

I understand that you are proposing something in this sentence:

An ideal HWP would just use hwp_angle, storing the value 2πνt into it (with ν being the angular speed of the HWP), and ignore the presence of bore2ecl_quaternions. However, in a derived class there might be the simulation of some systematic effects that introduces some noise in the pointings; in the latter case, the code is free to modify quaternions by multiplying each of them by a new quaternion on the left.

but I am not sure if this can be used to the systematics I have in mind.

[ziotom78]

With your new proposal, these effects can be included using a time variable detector-boresight quaternion. The only issue, is that in case of an HWP synchronous effect, the sampling of the pointing must be as fast as the sampling of the detectors.

Yes, in fact this is what I was saying in this paragraph:

Note that we do not want to use TimeDependentQuaternions here, as that class is optimized for low-sampling quaternions, while a HWP spins quickly! Instead, apply_hwp_to_pointings gets the fully sampled quaternions.

Thus, they should be sampled at 19 Hz, the same sampling frequency of the detectors.

[yusuke-takase]

I think the proposal will work effectively to simulate several pointing systematics, especially for the pointing disturbance due to vibrations.
I'm not sure that I should comment it now, but let me confirm. If we focus on some physical vibration for example due to a cooler, which is mounted on the spacecraft, in other word, the vibration model is characterized by xyz coordinate in payload with specific frequency spectra. So if we want to estimate this effect correctly, we may multiply quaternions which describes vibration to spin quaternion before going to go Ecliptic quaternions. Is it correct understanding of your structure, @ziotom78 ?

If I comment on the HWP wedge effect, I simulated the systematics by Falcons 3 years ago, and way was just I generate true pointing and wrong pointing i.e. with constant offset due to a wedge, and then rotate the wrong pointing with HWP rotation rate around true pointing as rotation axis. Every calculation scheme was done in 19Hz. However, if I think about implementation again now, I will buy @ziotom78 's aproach because it is more coherent.

[ziotom78]

So if we want to estimate this effect correctly, we may multiply quaternions which describes vibration to spin quaternion before going to go Ecliptic quaternions. Is it correct understanding of your structure, @ziotom78 ?

Yes, that was my idea, I believe that in this way it should be reasonably straightforward to implement.

[ziotom78]

PR #319 is a tentative implementation of the ideas in this discussion.