Map-making#

The primary aim of the LiteBIRD Simulation Framework is to create synthetic timestreams as if the real LiteBIRD acquired them. The process of creating maps out of these timelines is called map-making and it is strictly considered a data-analysis task, not a simulation task. However, since most of the assessments on the quality of a timestream can only be done on maps, the LiteBIRD Simulation Framework provides some facilities to produce maps out of timestreams. These maps are created using the Healpix pixelization scheme and saved in FITS files.

The framework provides the following solutions:

A binner, i.e., a simple map-maker that assumes that only uncorrelated noise is present in the timelines.
A pair-differencing map-maker that solves for Q/U only by

differencing T/B detector pairs on the same wafer pixel.

A destriper, i.e., a more advanced map-maker that can remove the effect of correlated instrumental noise from the timelines before producing a map. We usually refer to correlated noise as 1/f noise, and the purpose of the destriper is to estimate its contribution and remove it from the timelines; then, a classical binner is run over the cleaned timelines.
You can also use save_simulation_for_madam() to save TODs and pointing information to disk and then manually call the Madam map-maker.
An interface to use the BrahMap <https://github.com/anand-avinash/BrahMap> GLS mapmaker. (To use this, you must ensure that the package is installed, thus see BrahMap documentation <https://anand-avinash.github.io/BrahMap/>.)

In this chapter, we assume you have already created the timelines to be used as input to the destriper. Here is a sample code that creates a simple timeline containing white noise for two detectors:

import numpy as np
import astropy.units as u
from numpy.random import MT19937, RandomState, SeedSequence

import litebird_sim as lbs

sim = lbs.Simulation(
    base_path="destriper_output",
    start_time=0,
    duration_s=86400.0,
    random_seed=12345,
)

sim.set_scanning_strategy(
    scanning_strategy=lbs.SpinningScanningStrategy(
        spin_sun_angle_rad=np.deg2rad(30),  # CORE-specific parameter
        spin_rate_hz=0.5 / 60,  # Ditto
        # We use astropy to convert the period (4 days) in
        # seconds
        precession_rate_hz=1.0 / (4 * u.day).to("s").value,
    )
)
instr = lbs.InstrumentInfo(
    name="core",
    spin_boresight_angle_rad=np.deg2rad(65),
)

# We create two detectors, whose polarization angles are separated by π/2
sim.create_observations(
    detectors=[
        lbs.DetectorInfo(name="0A", sampling_rate_hz=10),
        lbs.DetectorInfo(
            name="0B", sampling_rate_hz=10, quat=lbs.quat_rotation_z(np.pi / 2)
        ),
    ],
    tod_dtype=np.float64,
    n_blocks_time=lbs.MPI_COMM_WORLD.size,
    split_list_over_processes=False,
)

# Generate some white noise
rs = RandomState(MT19937(SeedSequence(123456789)))
for curobs in sim.observations:
    curobs.tod *= 0.0
    curobs.tod += rs.randn(*curobs.tod.shape)

The paper [6] explains the theory of binners and destripers. Our implementation closely follows the terminology used in the paper. In this chapter, we will refer to the paper as KS2009.

Binner #

Once you have generated a set of observations, either on a single process or distributed over several mpi processes, you can create a simple binned map with the function make_binned_map(). This function takes one or more Observation objects and the Healpix resolution of the output map (nside), and it produces a coadded map. The algorithm assumes white noise and each detector gets weighted by \(1 / \text{NET}^2\). If the pointing information is not provided in the observation, it can be passed through the optional argument pointings, with a syntax similar to scan_map_in_observations(). The output map is in Galactic coordinates, but you can specify the coordinate system you want via the parameter output_coordinates. This is how it should be called:

result = lbs.make_binned_map(nside=128, observations=observations)
healpy.mollview(result.binned_map[0])

The function obtains pointing information from the Observation (either using the precomputed pointing or computed on the fly). As an alternative, you can provide pointings as a list of numpy arrays. The result is an instance of the class BinnerResult and contains both the I/Q/U maps and the covariance matrix.

Function make_binned_map() has a high-level interface in the class Simulation that bins the content of the observations into maps. The syntax is identical to make_binned_map():

result = sim.make_binned_map(nside=nside)
healpy.mollview(result.binned_map[0])

The BinnerResult class contains the field binned_map, which is a Healpix map containing the binned values of the samples.

Using a simple example, we now explain how the function calculates the binned map. Suppose we measured just two pixels in the map with only seven measurements. (These numbers are ridiculous, but in this way, the matrices we will write will be more manageable!)

Each sample is associated with a direction in the sky (the pointing information), but from the point of view of the binner, the only quantity that matters is the pixel index in the map associated with the pointing direction. In our example, the seven samples will measure just two pixels in the sky with different attack (polarization) angles. The following figure and table show how the seven TOD samples observe the two pixels. The bars refer to the polarization angle of the detector.

_images/destriper-tod-angles.svg — Toy-example of a TOD containing 7 samples.#

The figure shows how the seven samples observe each of the two pixels in the map. The thick bars represent the direction of the attack angle ψ, which coincides with the polarization angle of the detector projected on the center of the pixel.

Samples in the TOD, their polarization angle ψ, and the pixel they hit#
#	ψ	Pixel
1	90°	(1)
2	30°	(1)
3	90°	(2)
4	15°	(2)
5	30°	(2)
6	45°	(2)
7	60°	(1)

A fundamental quantity in KS2009 is the pointing matrix matrix \(P\): it is a \((N_t, 3N_p)\) matrix where \(N_t\) is the number of samples in the TOD (7 in our case) and \(N_p\) is the number of pixels observed by the TOD (2 in our case). Thus, for our simple example \(P\) is

\[\begin{split}P = \begin{pmatrix} \textcolor{#002683}{1}& \textcolor{#002683}{\cos90^\circ}& \textcolor{#002683}{\sin90^\circ}& 0& 0& 0\\ \textcolor{#002683}{1}& \textcolor{#002683}{\cos30^\circ}& \textcolor{#002683}{\sin30^\circ}& 0& 0& 0\\ 0& 0& 0& \textcolor{#268300}{1}& \textcolor{#268300}{\cos90^\circ}& \textcolor{#268300}{\sin90^\circ}\\ 0& 0& 0& \textcolor{#268300}{1}& \textcolor{#268300}{\cos15^\circ}& \textcolor{#268300}{\sin15^\circ}\\ 0& 0& 0& \textcolor{#268300}{1}& \textcolor{#268300}{\cos30^\circ}& \textcolor{#268300}{\sin30^\circ}\\ 0& 0& 0& \textcolor{#268300}{1}& \textcolor{#268300}{\cos45^\circ}& \textcolor{#268300}{\sin45^\circ}\\ \textcolor{#002683}{1}& \textcolor{#002683}{\cos60^\circ}& \textcolor{#002683}{\sin60^\circ}& 0& 0& 0 \end{pmatrix},\end{split}\]

where the number of rows matches the number of samples in the TOD (7), and for each pixel there are three columns corresponding to I, Q, and U.

Enough samples must cover each pixel: there must be at least three non-degenerate polarization angles to recover the three Stokes parameters associated with it (I, Q, and U). KS2009 describes how to compute them through the matrix \(M = P^T\cdot C_w^{-1}\cdot P\), where \(C_w\) is a diagonal matrix with shape \((N_t, N_t)\), where each element along the diagonal is \(\sigma^2\), the white-noise variance of the sample:

\[\begin{split}C_w = \begin{pmatrix} \sigma^2& 0& 0& 0& 0& 0& 0\\ 0& \sigma^2& 0& 0& 0& 0& 0\\ 0& 0& \sigma^2& 0& 0& 0& 0\\ 0& 0& 0& \sigma^2& 0& 0& 0\\ 0& 0& 0& 0& \sigma^2& 0& 0\\ 0& 0& 0& 0& 0& \sigma^2& 0\\ 0& 0& 0& 0& 0& 0& \sigma^2 \end{pmatrix}\end{split}\]

The square matrix \(M = P^T\cdot C_w^{-1}\cdot P\) has shape \((3 N_p, 3 N_p)\): each pixel in the map is associated to a 3×3 block along the diagonal:

\[\begin{split}M = \begin{pmatrix} \textcolor{#002683}{\frac3{\sigma^2}}& \textcolor{#002683}{-\frac1{\sigma^2}}& \textcolor{#002683}{\frac{\sqrt3}{\sigma^2}}& 0& 0& 0\\ \textcolor{#002683}{-\frac1{\sigma^2}}& \textcolor{#002683}{\frac3{2\sigma^2}}& \textcolor{#002683}{0}& 0& 0& 0\\ \textcolor{#002683}{\frac{\sqrt3}{\sigma^2}}& \textcolor{#002683}{0}& \textcolor{#002683}{\frac{3}{2\sigma^2}}& 0& 0& 0\\ 0& 0& 0& \textcolor{#268300}{\frac4{\sigma^2}}& \textcolor{#268300}{\frac{\sqrt3 - 1}{2\sigma^2}}& \textcolor{#268300}{\frac{\sqrt3 + 3}{2\sigma^2}}\\ 0& 0& 0& \textcolor{#268300}{\frac{\sqrt3 - 1}{2\sigma^2}}& \textcolor{#268300}{\frac2{\sigma^2}}& \textcolor{#268300}{\frac{\sqrt3}{2\sigma^2}}\\ 0& 0& 0& \textcolor{#268300}{\frac{\sqrt3 + 3}{2\sigma^2}}& \textcolor{#268300}{\frac{\sqrt 3}{2\sigma^2}}& \textcolor{#268300}{\frac2{\sigma^2}} \end{pmatrix}\end{split}\]

If we assigned different values for \(\sigma\) to the seven elements along the diagonal of \(C_w\), then the 3×3 subblocks of the \(M\) matrix would have changed to reflect the fact that those samples with smaller \(\sigma\) have a greater weight in determining the value of the Stokes parameter.

Once we have the matrices \(P\), \(C_w\), and \(M\), determining the value of the three Stokes parameters I, Q, U for each pixel is trivial. The signal \(s\) measured by a detector is a function of \(I\), \(Q\), \(U\), and the polarization angle \(\psi\):

\[s = I + Q \cos2\psi + U\sin2\psi\]

and thus we can exploit the redundancy of the measurements per each pixel (remember: at least three measurements per pixel with non-degenerate angles!) by running a simple \(\chi^2\) minimization assuming Gaussian noise. The result is the following:

\[m = M^{-1}\cdot P^T\cdot C_w^{-1}\cdot y,\]

where \(y\) is a 7-element vector containing the TOD samples \(s_i\) and \(m\) is a 6-element vector containing the estimates of the Stokes parameters I, Q, and U for the two pixels in the map according to the following order:

\[m = \begin{pmatrix}I_1& Q_1& U_1& I_2& Q_2& U_2\end{pmatrix}\]

The presence of the matrix \(C_w^{-1}\) scales the value of each sample so that the ones affected by larger noise will be made smaller; the \(P^T\) factor sums all the samples in the TOD that fall within the same pixel, and finally \(M^{-1}\) “solves” the linear system for the three parameters I, Q, and U per each pixel.

Once the call to make_binned_map() ends, the field invnpp of the BinnerResult object returned by the function contains an array with shape \((N_p, 3, 3)\), where each 3×3 block is the inverse of the sub-matrix \(M_i\) for the i-th pixel.

Data splits #

The function make_binned_map() is also able to provide data splits both in time and in detector space. The time split is passed through a string in the parameter time_split; the same applies for detector_split. Currently, the supported time splits are:

"full" (default): no split.
"first_half": use the first half of the TOD.
"second_half": use the second half of the TOD.
"odd": use the odd samples of the TOD.
"even": use the even samples of the TOD.
"yearX": use the TOD samples relative to the X-th year of observation. The X-th part must be an integer between 1 and 3.
"surveyX": assuming that a complete survey is performed in 6 months, use the X-th survey. X must be an integer between 1 and 6.

There are three supported options for the detector split:

"full" (default): no split is performed.
"waferXXX": use the detectors in the wafer XXX. The XXX part must specify the wafer (e.g. “L00”, “M01”, or “H02”).

The final data split will correspond to the combination of the two splits.

The function check_valid_splits() will check whether the requested split is part of the list above. If the split is not valid, the function will raise a ValueError. In addition, it will check whether the requested split is compatible with the observation duration and the detector list. Thus, for example, if the observation lasts one year, the split year2 will raise an AssertionError. Similarly, if some detector in the L00 wafer does the observation, waferL03 will also raise an AssertionError.

Pair Differencing (Q/U only)#

The function make_pair_differenced_map() implements QU map-making from detector pairs. For each pair of detectors in the same wafer pixel, the map-maker differences the timelines (T - B) and solves a 2x2 system for Q and U in each sky pixel.

If \(\psi_T\) and \(\psi_B\) are the effective polarization angles for the two detectors in the pair, the model used in map-making is

\[d_T - d_B = Q\,[\cos(2\psi_T) - \cos(2\psi_B)] + U\,[\sin(2\psi_T) - \sin(2\psi_B)].\]

The pair weight is the average of the detector weights:

\[w_{\mathrm{pair}} = \frac{w_T + w_B}{2}.\]

This is the low-level usage:

result = lbs.make_pair_differenced_map(nside=128, observations=observations)
healpy.mollview(result.binned_map[0])  # Q
healpy.mollview(result.binned_map[1])  # U

The PairDifferencingResult output contains:

binned_map with shape (2, N_pix) (Q/U only),
invnpp with shape (N_pix, 2, 2).

The corresponding high-level interface is Simulation.make_pair_differenced_map():

result = sim.make_pair_differenced_map(nside=nside)

Constraints for this map-maker:

each local observation must provide detector attributes pol, wafer, and pixel;
detectors must form valid T/B pairs within each (wafer, pixel) group;
paired detectors must point to the same sky pixel for the selected samples.

Time and detector splits are supported with the same semantics of make_binned_map().

Destriper #

If you know that your simulation contains 1/f noise, you should avoid using the binner and instead employ a better map-making algorithm. With 1/f noise, the hypothesis that the noise be white drops and the binning equation is no longer valid.

The destriping algorithm is implemented by the make_destriped_map(), functionally equivalent to make_binned_map(). However, as the algorithm it implements is more complex, you must provide more information when calling it. Specifically, you should instantiate an instance of the class DestriperParameters:

params = DestriperParameters(
    ...
)

result = lbs.make_destriped_map(
    nside=nside,
    observations=observations,
    params=params,
)
healpy.mollview(result.destriped_map)

Pointing information is handled identically to make_binned_map(). The result is an instance of the class DestriperResult, which is similar to BinnerResult but it contains much more information.

Function make_destriped_map() has a high level interface in the class Simulation that applies the destriper algorithm to all the observations in the simulation. The syntax is identical to make_destriped_map():

result = sim.make_destriped_map(nside=nside)
healpy.mollview(result.destriped_map)

We will now explain how a destriper works and what the meaning of each parameter in the classes DestriperParameters and DestriperResult is. Apart from KS2009, another source of information is the file test/test_destriper.py: it compares the results of the destriper with the analytical solution of the simple 7-sample model we discuss here.

The idea of a destriper is to group consecutive samples in the same TOD into different baselines. Each baseline must contain a number of samples such that the noise within it will be roughly white; thus, the baseline should not contain too many samples! A good rule of thumb is to make the time span covered by one baseline shorter than the inverse of the knee frequency for that detector. The destriper works by assuming that the effect of 1/f noise on all the samples within the baselines is a purely additive factor; this is an approximation, but for experiments like Planck, it has worked well.

_images/destriper-baselines.svg — Baselines in a destriper#

A destriper seeing the 13 samples in the figure could decide to group them into two baselines: 7 within the first baseline, and 6 within the second baseline. The level \(a_i\) of each baseline is represented by a dashed horizontal line.

The purpose of the destriper is to model 1/f noise using the baselines and produce a cleaned map. It does so by doing the following steps:

It estimates the value of the baselines by minimizing a \(\chi^2\) function (more on this later);
It removes the value of the baselines from each sample in the TOD;
It applies a binning operation on the cleaned TOD, using the same equations explained in the section Binner.

To provide a concrete example of how the destriper works, we will resume the toy TOD containing seven samples that we discussed in the section Binner. We will group these samples into two baselines of 4 and 3 samples, respectively.

It’s easier to understand what the algorithm does if we play following the rules of the destriper: instead of thinking of proper 1/f noise, let’s pretend that what we have in the TOD are the baselines: constant values that last for some time and are added to the samples in the TOD. How can the destriper understand the level of each baseline in the figure above? The trick is easy to understand if we recall Table Samples in the TOD, their polarization angle ψ, and the pixel they hit. Let’s consider the first pixel: it was measured by the first two samples, then the detector moved to another position in the sky (the second pixel), but it moved back to the first pixel just while measuring the last sample in the TOD (#7). The point is that the first two samples #1 and #2 belong to the first baseline, since it lasts 4 pixels, while sample #7 belongs to the second one. Yet, once we take 1/f noise into account, the estimated I/Q/U parameter for this pixel must agree. We can rewrite this problem as a \(\chi^2\) minimization problem: the destriper looks for the values of the two baselines that minimize the discrepancy in the values of the pixels as estimated by every sample in the TOD.

Apart from the standard matrices \(C_w\), \(P\), and \(M\), the destriper requires a new matrix \(F,\) whose shape is \((N_t, N_b)\):

\[\begin{split}F = \begin{pmatrix} 1& 0\\ 1& 0\\ 1& 0\\ 1& 0\\ 0& 1\\ 0& 1\\ 0& 1 \end{pmatrix}\end{split}\]

Don’t be fooled by the fact that the number of columns is the same as the number of pixels in the map: in this case, the number of columns corresponds to the number of baselines in the TOD! The operator \(F\) tells which baseline “owns” the samples in the TOD, and, as you can see, we are assigning the first four samples to the first baseline and the last three samples to the last baseline. Applying \(F\) to a set of \(N_t\) samples (the TOD) produces a vector containing \(N_b\) samples (the baselines). As it is always the case that \(N_b \ll N_t,\) this means that \(F\) is a very tall and narrow matrix!

The purpose of \(F\) is to “project” the values of the two baselines into a 7-element TOD:

\[\begin{split}F \cdot \begin{pmatrix}a_1\\a_2\end{pmatrix} = \begin{pmatrix} a_1\\ a_1\\ a_1\\ a_1\\ a_2\\ a_2\\ a_2 \end{pmatrix}.\end{split}\]

The transpose \(F^T\) represents a linear operator that takes a 7-element TOD and sums up all the elements belonging to the same baseline:

\[\begin{split}F^T \cdot \begin{pmatrix} y_1\\ y_2\\ y_3\\ y_4\\ y_5\\ y_6\\ y_7 \end{pmatrix} = \begin{pmatrix}y_1 + y_2 + y_3 + y_4\\ y_5 + y_6 + y_7\end{pmatrix}.\end{split}\]

At the basis of the destriper there is the operator \(Z,\) which is defined as follows:

\[Z = I - P\cdot M^{-1}\cdot P^T\cdot C_w^{-1}\]

(\(I\) is the identity operator), and it is applied to a TOD: \(Z\cdot y\). The purpose of the operator is to “clean up” the TOD from all the components that are not white noise. It does so by creating a map (note the presence on the right of the map-binning operator \(M^{-1}\cdot P^T\cdot C_w^{-1}\) we discussed before) and then scanning the map back into a TOD (the last operator \(P\) on the left of the map-binning operator). Because of the difference between the identity operator \(I\) and this matrix, the result of this “map-and-scan” operation is subtracted from the TOD sample.

The destriper estimates the vector of baselines \(a\) by solving iteratively the following equation, which is the solution of a minimization problem on the \(\chi^2\) we qualitatively discussed before:

\[\left(F^T\cdot C_w^{-1}\cdot Z\cdot F\right) a = F^T\cdot C_w^{-1}\cdot Z\cdot y\]

This equation is of the form \(Ax = b,\) where \(A\) and \(b\) are known and \(x\) is the quantity to determine. (The symbol \(A\) is commonly found in linear algebra textbooks, so we will switch to this notation even if KS2009 uses \(D\).) The solution would just be \(x = A^{-1}b\), but there are two problems with this formula:

Matrix \(A\) is too large to be inverted;
Matrix \(A\) is not invertible!

The second point is alarming, and it might sound like an irrecoverable problem. The fact that \(\det A = 0\) stems from the fact that the solution to the destriping equation is not unique, as if \(a\) is a solution, then \(a + K\) is still a solution for any scalar constant \(K.\) (Remember, the purpose of the destriper is to make the many measurements for I/Q/U within each pixel consistent once it subtracts the baselines and the consistency is the same if it shifts all the baselines by the same amount!) We are looking for a solution \(a\) that is orthogonal to the null space of \(A\).

We can solve both problems by employing the fact that \(A\) is a symmetric semi-definite matrix and thus employing the so-called Conjugate Gradient (CG) method, which is able to produce a solution for the problem \(Ax = b\) even if \(A\) is singular, because it simply minimizes the quantity \(Ax - b\) without trying to invert \(A\). If you are curious about the details of the algorithm, a good pedagogical reference is [2].

The way the algorithm works is to start from a guess for \(x\) (the set of baselines), which must be orthogonal to the null-space of \(A\); in our case, it is enough to require that the mean value of the first guess of the baselines is zero, i.e., \(\sum_i x_i = 0\). Starting from guess \(x^{(0)},\) the algorithm produces a new guess \(x^{(1)}\) such that \(r^{(1)} = Ax^{(1)} - b\) is closer to zero than \(r^{(0)} = Ax^{(0)} - b\). The procedure keeps going until the residual \(r^{(n)} = Ax^{(n)} - b\) is small enough, always satisfying the fact that \(\sum_i x_i^{(n)} = 0\) for any \(n.\) Thus, the final solution will have the property that the average value of the baselines will be zero.

The CG algorithm requires iterations to continue until the residuals \(r^{(n)}\) are “small enough”. But how can we tell this? Vector \(r^{(n)}\) contains \(N_b\) elements, corresponding to the residual for each baseline, and it might be that some of the residuals are small and some are not. The most common way to deal with this is to compute some norm over \(r^{(n)}\) to produce a single scalar and then to check whether this value is below some threshold. There are several possible definitions for the norm:

The standard \(L_2\) norm (Euclidean): \(\left\|r^{(n)}\right\|^2\);
The \(L_\infty\) norm: \(\max_i\left|r_i^{(n)}\right|\).

Function make_destriped_map() adopts \(L_\infty\): it is stricter than \(L_2\), because it does not wash out the presence of a few large residuals even if all the others are negligible.

The convergence of the Conjugated Gradient can be controlled by the following fields in the DestriperParameters class:

iter_max: the maximum number of iterations, i.e., the upper limit for \(n\).
threshold: the value to be reached by the \(L_\infty\) norm when applied to the residuals \(r^{(n)}\): if the norm is smaller than threshold, then the CG algorithm stops.
samples_per_baseline: this can either be an integer, in which case it will be used for all the baselines, or a list of 1D arrays, each containing the length of each baseline for each observation passed through the parameter observations. Note that if you provide an integer, it might be that not all baselines will have exactly that length: it depends on whether the number \(N_t\) of samples in the TOD is evenly divisible by samples_per_baseline or not. The algorithm tries to make the number of elements per baseline as evenly as possible.
use_preconditioner: if this flag is True, the preconditioning matrix \(F^T\cdot C_w^{-1}\cdot F\) will be used with the CG algorithm. This might speed up the convergence.

Once the destriper has completed the process, the following fields in the DestriperResult object can be inspected to check how the CG iterations went:

converged: a Boolean flag telling whether the threshold was reached (True) or not (False)
history_of_stopping_factors: a list of values of the \(L_\infty\) norm applied to the residuals \(r^{(n)},\) one per each iteration of the CG algorithm. Inspecting those values might help in understanding if the destriper was not able to converge because of a too small value of iter_max.

The baselines are saved in the field baselines of the DestriperResult class; this is a list of 2D arrays, where each element in the list is associated with one of the observations passed in the parameter observations. The shape of each 2D arrays is \((N_d, N_b),\) where \(N_d\) is the number of detectors for the observation and \(N_b\) is the number of baselines. A visual representation of the memory layout of the field baselines is shown in the following figure.

_images/destriper-baselines-memory-layout.svg — Memory layout for baselines.#

The image assumes the baselines contain data for two detectors, A and B. In the upper image, the baselines are shown as horizontal lines that span the time range covered by the two (consecutive) observations. We use continuous lines for the baselines of detector A and dashed lines for B. Note that the number of baselines in the first Observation object (5) differs from the number in the second object (3). The lower image offers a visual representation of the layout of the baselines field in the class DestriperResult: a Python list of two elements, each containing 2D arrays with shape \((N_d, N_b)\). The i-th row of each 2D array contains the baselines for the i-th detector.

The field baseline_errors has the same memory layout as baselines and contains a rough estimate of the error per each baseline, assuming that the noise in the baselines is not correlated. (Unrealistic!) Finally, the field baseline_lengths is a list of 1D integer arrays of \(N_b\) elements containing the number of samples in each baseline; it should match the value provided in the field samples_per_baseline in the class DestriperParameters.

Once the destriper has computed a solution for the destriping equation, you can ask to remove the baselines from the TOD using one of the functions remove_baselines_from_tod() or remove_destriper_baselines_from_tod().

You can save the results of the destriper using the function save_destriper_results() and load them back again with the load_destriper_results(). Note that if you are running your code using MPI, you should call both functions on all the MPI processes, and the number of processes should be the same between the two calls.

In addition, func:.make_destriped_map can accept a custom set of baselines as input. The destriper will skip the CG iterations and proceed directly to the map-making step if these have the correct dimensions.

How the N_obs matrix is stored #

The destriper uses a different method to store the matrix \(M\) in memory. As the 3×3 sub-blocks of this matrix need to be inverted often during the CG process, the function make_destriped_map() decomposes each block using Cholesky decomposition: each 3×3 matrix \(M_i\) associated with the pixels in the map is decomposed into

\[M_i = L_i L_i^T,\]

where \(L\) is a lower-triangular matrix, and only the nonzero coefficients of the matrix are saved. The advantages are twofold:

Less memory is required;
Calculating the inverse \(M^{-1}\) is faster and more accurate.

The fields are saved in an instance of the class NobsMatrix, which contains the following fields:

nobs_matrix is a 2D array of shape \((N_p, 6)\), containing the nonzero components of the \(N_p\) matrices;
valid_pixel is a 1D Boolean array containing \(N_p\) elements: each of them is True if the corresponding matrix \(M_i\) was invertible (i.e., it was observed with at least three non-degenerate attack angles), False otherwise;
is_cholesky is a flag that tells whether nobs_matrix contains the nonzero coefficients of the many \(L_i\) matrices or the lower triangular part of \(M_i\). After a successful call to make_destriped_map(), this should always be True: it is set to False during the execution of the destriper but updated to True before the CG iteration starts.

Given that it is often useful to have access to matrix \(M^{-1}\), the method NobsMatrix.get_invnpp() computes this inverse as a 3D array with shape \((N_pix, 3, 3)\), where the first index runs over all the pixels and the last two dimensions are used to store the value of \(M_i^{-1}\).

Data splits #

Similarly to the function make_binned_map(), make_destriped_map() too can provide data splits both in time and detector space. Given that the splits are implemented in the same way, refer to the documentation of make_binned_map() for more details.

The difference is that each split is applied to the TOD before the destriper runs. Therefore, the destriper will only use the samples within the split to compute the baselines.

BrahMap GLS mapmaker #

If you have installed the BrahMap, you can use it seemlessly within the LiteBIRD Simulation Framework through make_brahmap_gls_map(). This provides a high-level interface with BrahMap to perform optimal mapmaking with different noise covariance operators. These operators can be defined using BrahMap and can be used in the interface directly.

The MPI communicator used in map-making must be the one that contains exclusively all the data needed for map-making. In case of litebird_sim, the communicator lbs.MPI_COMM_WORLD. Therefore, it is a suitable communicator to be used in map-making. Hence, communicator used by BrahMap must be updated as following before using any other BrahMap function with litebird_sim data: brahmap.MPI_UTILS.update_communicator(comm=lbs.MPI_COMM_WORLD)

# Creating the inverse white noise covariance operator
inv_cov = brahmap.InvNoiseCovLO_Uncorrelated(
    diag=[...],       # Diagonal elements of the inverse of white noise covariance
                      # matrix

    dtype=np.float64, # Numerical precision of the operator
)

# Performing the GLS map-making
gls_result = brahmap.LBSim_compute_GLS_maps(
    nside=nside,                    # Nside parameter for the output healpix map
    observations=sim.observations,  # List of observations from litebird_sim
    components="tod",                # TOD components to be used in map-making
    inv_noise_cov_operator=inv_cov, # Inverse noise covariance operator
    dtype_float=np.float64,         # Numerical precision to be used in map-making
)

gls_result obtained above is an instance of the class LBSimGLSResult. The output maps can be accessed from this object with gls_result.GLS_maps. For further details refer to BrahMap documentation <https://anand-avinash.github.io/BrahMap/>.

High level interface and data-splits #

The functions Simulation.make_binned_map(), Simulation.make_pair_differenced_map(), Simulation.make_destriped_map(), and Simulation.make_brahmap_gls_map() provide high-level access to the main map-making pipelines offered in litebird_sim. These interfaces simplify the process of generating maps from simulated data and are well-integrated into the overall simulation flow via the Simulation class.

In addition to single-map outputs, these methods support automatic splitting of the data along both time and detector axes, enabling efficient generation of subsets for validation, null tests, and systematics control.

Since the class Simulation is interfaced to litebird_sim.make_binned_map(), litebird_sim.make_pair_differenced_map(), and litebird_sim.make_destriped_map(), it is able to provide data splits both in time and detector space (see Map-making for more details on the splitting and the available options). In addition, the class contains Simulation.make_binned_map_splits(), which is a wrapper around litebird_sim.make_binned_map(), Simulation.make_pair_differenced_map_splits(), which is a wrapper around litebird_sim.make_pair_differenced_map(), and Simulation.make_destriped_map_splits(), which is a wrapper around litebird_sim.make_destriped_map(). These allow performing the mapmaking on multiple choices of splits at once (passed as a list of strings). Indeed, the functions will loop over the cartesian product of the time and detector splits. The default behavior is to perform the mapmaking in each combination and save the result to disk, to avoid memory issues. In particular, in the case of Simulation.make_binned_map_splits() and Simulation.make_pair_differenced_map_splits(), only the maps (I/Q/U for the binner; Q/U for pair-differencing) are saved to disk, unless you set include_inv_covariance to True. This saves the elements of the inverse covariance as extra fields in the output FITS file. This default behavior can be changed by setting the write_to_disk parameter to False. Then, the function returns a dictionary containing the full results of the mapmaking for each split. The keys are strings that describe the split in the format "{Dsplit}_{Tsplit}", such as "waferL00_year1". On the other hand, the default behavior of Simulation.make_destriped_map_splits() is to save to disk the complete mapmaking.DestriperResult class for each split as a FITS file. To avoid this and get a dictionary similar to the one returned by Simulation.make_binned_map_splits(), you can set the write_to_disk parameter to False.

The method Simulation.make_destriped_map_splits() also offers the possibility to recycle the baseline computed from the largest split available. Indeed, if the flag recycle_baselines is set to True, that method enforces the computation of the "full_full" split and then reuses the baselines computed for that split for all the other splits.

Before performing the computation, the function Simulation.check_valid_splits() will check whether the requested split is valid (see Map-making). This is a wrapper around litebird_sim.check_valid_splits(). If the split is not valid, the function will raise a ValueError. In addition, it will check whether the requested split is compatible with the duration of the observation and with the detector list. Thus, for example, if the observation lasts 1 year, the split "year2" will raise an AssertionError. Similarly, if the observation is performed with some detector contained in the L00 wafer, the split "waferL03" will also raise an AssertionError.

Saving files for Madam #

The function save_simulation_for_madam() takes a Simulation object, a list of detectors and a output path (the default is a subfolder of the output path of the simulation) and saves a set of files in it:

Pointing information, saved as FITS files;
TOD samples, saved as FITS files;
A so-called «simulation file», named madam.sim;
A so-called «parameter file», named madam.par.

These files are ready to be used with the Madam map-maker; you just need to pass the parameter file madam.par to one of the executables provided by Madam (the other ones are referenced by the parameter file). For instance, the following command will compute the amount of memory needed to run Madam:

$ inputcheck madam.par

The following command will run Madam:

$ madam madam.par

Of course, in a realistic situation you want to run madam using MPI, so you should call mpiexec, mpirun, or something similar.

Creating several maps with Madam #

There are cases where you want to create several maps out of one simulation. A common case is when you simulate several components to be put in the TOD and store them in different fields within each Observation object:

sim.create_observations(params)

for cur_obs in sim.observations:
    # We'll include several components in the signal:
    # CMB, white noise, 1/f noise, and dipole
    cur_obs.wn_tod = np.zeros_like(cur_obs.tod)
    cur_obs.oof_tod = np.zeros_like(cur_obs.tod)
    cur_obs.cmb_tod = np.zeros_like(cur_obs.tod)
    cur_obs.dip_tod = np.zeros_like(cur_obs.tod)

    # Now fill each of the *_tod fields appropriately

In cases like this, it is often useful to generate several sets of maps that include different subsets of the components:

A map including just the CMB and white noise;
A map including CMB, white noise and 1/f, but not the dipole;
A map including all the components.

You could of course call save_simulation_for_madam() three times, but this would be a waste of space because you would end up with three identical copies of the FITS file containing the pointings, and the last set of FITS files would contain the same components that were saved for the first two maps. (Remember, save_simulation_for_madam() saves both the pointings and the TODs!)

A trick to avoid wasting so much space is to save the FITS files only once, including all the TOD components, and then call save_simulation_for_madam() again for each other map using save_pointings=False and save_tods=False:

# This code will generate *three* sets of Madam files:
# - One including the CMB and white noise
# - One including 1/f as well
# - The last one will include the dipole too
save_files = True

# Iterate over all the maps we want to produce. For each of
# them we specify the name of the subfolder where the Madam
# files will be saved and the list of components to include
for (subfolder, components_to_bin) in [
    ("cmb+wn", ["wn_tod", "cmb_tod"]),
    ("cmb+wn+1f", ["wn_tod", "oof_tod", "cmb_tod"]),
    ("cmb+wn+1f+dip", ["wn_tod", "oof_tod", "cmb_tod", "dip_tod"]),
]:
    save_simulation_for_madam(
        sim=sim,
        params=params,
        madam_subfolder_name=subfolder,
        components=["cmb_tod", "wn_tod", "oof_tod", "dipole_tod"],
        components_to_bin=components_to_bin,
        save_pointings=save_files,
        save_tods=save_files,
     )

     # Set this to False for all the following iterations (maps)
     save_files = False

It is important that the components parameter in the call to save_simulation_for_madam() list all the components, even if you ask Madam not to use them in the first and second maps. The reason is that the function uses this parameter to create a «map» of the components as they are supposed to appear in the FITS files; for example, the cmb_tod field is the third in each TOD file, but this would not be apparent while producing the first map, where it is the second in the list of components that must be used. The .par file will list the components that need to be used to create the map, so it will not be confused if the TOD FITS files contain more components than needed. (This is a neat feature of Madam.)

Note

To understand how this kind of stuff works, it is useful to recap how Madam works, as the possibility to reuse TODs for different maps is linked to the fact that Madam requires two files to be run: the parameter file and the simulation file.

The simulation file represents a «map» of the content of a set of FITS files. A simulation file includes no information about the map-making process; it just tells how many FITS files to read and what is inside each.

The parameter file tells Madam how to create the maps. You can ask Madam to skip parts of the TODs in the parameter file. For example, you do not want to include the dipole in the output map.

When you call save_simulation_for_madam(), the components parameter is used to build the simulation file: thus, if you plan to build more than one map out of the same set of components, you want to have the very same simulation files because they «describe» what’s in the FITS files. This is the reason why we passed the same value to components every time we called save_simulation_for_madam().

However, when we create the three parameter files, each of them differs in the list of components that need to be included. If you inspect the three files cmb+wn/madam.par, cmb+wn+1f/madam.par, and cmb+wn+1f+dip/madam.par, you will see that they only differ for the following lines:

# cmb+wn/madam.par
tod_1 = wn_tod
tod_2 = cmb_tod

# cmb+wn+1f/madam.par
tod_1 = wn_tod
tod_2 = oof_tod
tod_3 = cmb_tod

# cmb+wn+1f+dip/madam.par
tod_1 = wn_tod
tod_2 = oof_tod
tod_3 = cmb_tod
tod_4 = dip_tod

That’s it. The lines with tod_* are enough to do all the magic to build the three maps.

Of course, once the three directories cmb+wn, cmb+wn+1f, and cmb+wn+1f+dip are created, Madam will run successfully only in the first one, cmb+wn. The reason is that only that directory includes the pointing and TOD FITS files! But if you are saving data on a filesystem that supports symbolic links, you can use them to make the files appear in the other directories too. For instance, the following commands will create them directly from a Unix shell (Bash, Sh, Zsh):

ln -srv cmb+wn/*.fits cmb+wn+1f
ln -srv cmb+wn/*.fits cmb+wn+1f+dip

(The -s flag asks to create soft links, the -r flag requires paths to be relative, and -v makes ln be more verbose.)

If you want a fully authomatic procedure, you can create the symbolic links in Python, taking advantage of the fact that save_simulation_for_madam() returns a dictionary containing the information needed to do this programmatically:

# First map
params1 = save_simulation_for_madam(sim, params)

# Second map
params2 = save_simulation_for_madam(
    sim,
    params,
    madam_subfolder_name="madam2",
    save_pointings=False,
    save_tods=False,
)

# Caution: you want to do this only within the first MPI process!
if litebird_sim.MPI_COMM_WORLD.rank == 0:
    for source_file, dest_file in zip(
        params1["tod_files"] + params1["pointing_files"],
        params2["tod_files"] + params2["pointing_files"],
    ):
        # This is a Path object associated with the symlink that we
        # want to create
        source_file_path = source_file["file_name"]
        dest_file_path = dest_file["file_name"]

        dest_file_path.symlink_to(source_file_path)

You can include this snippet of code in the script that calls save_simulation_for_madam(), so that the procedure will be 100% automated.

API reference #

class litebird_sim.mapmaking.BinnerResult(binned_map: Any = None, invnpp: Any = None, coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic, components: list | None = None, detector_split: str = 'full', time_split: str = 'full')#

Bases: object

Result of a call to the make_binned_map() function

This dataclass has the following fields:

binned_map: Healpix map containing the binned value for each pixel
invnpp: inverse of the covariance matrix element for each pixel in the map. It is an array of shape (12 * nside * nside, 3, 3)
coordinate_system: the coordinate system of the output maps (a CoordinateSistem object)
components: list of components included in the map, by default only the field tod is used
detector_split: detector split of the binned map
time_split: time split of the binned map

binned_map: Any = None#

components: list | None = None#

coordinate_system: CoordinateSystem = 1#

detector_split: str = 'full'#

invnpp: Any = None#

time_split: str = 'full'#

class litebird_sim.mapmaking.DestriperParameters(output_coordinate_system: CoordinateSystem = CoordinateSystem.Galactic, samples_per_baseline: int | list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | None = 100, iter_max: int = 100, threshold: float = 1e-07, use_preconditioner: bool = True)#

Bases: object

Parameters used by the function make_destriped_map()

The class is used to tell the destriper how to solve the map-making equation. The fields are:

nside (integer, power of two): resolution of the output map
output_coordinate_system (CoordinateSystem): the coordinate system of the output map
samples_per_baseline (integer): how many consecutive samples in the TOD must be assigned to the same baseline. If None, the destriper algorithm is skipped and a simple binning will be done.
iter_max (integer): maximum number of iterations for the Conjugate Gradient
threshold (float): minimum value of the discrepancy between the estimated baselines and the baselines deduced by the destriped map. The lower this value, the better the map produced by the destriper
use_preconditioner (Boolean): if True, use the preconditioned conjugate gradient algorithm. If False do not use the preconditioner. (The latter is probably useful only to debug the code.)

iter_max: int = 100#

output_coordinate_system: CoordinateSystem = 2#

samples_per_baseline: int | list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | None = 100#

threshold: float = 1e-07#

use_preconditioner: bool = True#

class litebird_sim.mapmaking.DestriperResult(params: DestriperParameters, nside: int, components: list[str], hit_map: ndarray[tuple[Any, ...], dtype[_ScalarT]], binned_map: ndarray[tuple[Any, ...], dtype[_ScalarT]], nobs_matrix_cholesky: NobsMatrix, coordinate_system: CoordinateSystem, baselines: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | ndarray[tuple[Any, ...], dtype[_ScalarT]] | None, baseline_errors: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | ndarray[tuple[Any, ...], dtype[_ScalarT]] | None, baseline_lengths: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | ndarray[tuple[Any, ...], dtype[_ScalarT]] | None, stopping_factor: float | None, history_of_stopping_factors: list[float], destriped_map: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None, converged: bool | str, elapsed_time_s: float, bytes_in_temporary_buffers: int, detector_split: str | None = 'full', time_split: str | None = 'full')#

Bases: object

Result of a call to make_destriped_map()

The fields baselines, baseline_lengths, stopping_factors, destriped_map, and converged are only relevant if you actually used the destriper; otherwise, they will be set to None.

If you are running several MPI processes, keep in mind that the fields hit_map, binned_map, destriped_map, nobs_matrix_cholesky, and bytes_in_temporary_buffers contain the global map, but baselines and baseline_lengths only refer to the TODs within the current MPI process.

List of fields:

params: an instance of the class DestriperParameters
nside: the resolution of the Healpix maps.
hit_map: a map of scalar values, each representing the sum of weights per pixel (normalized over the σ of each detector).
binned_map: the binned map (i.e., without subtracting the 1/f baselines estimated by the destriper).
nobs_matrix_cholesky: an instance of the class NobsMatrix.
coordinate_system: the coordinate system of the maps hit_map, binned_map, and destriped_map.
baselines: a Python list of NumPy arrays (one per each observation passed to the destriper) containing the value of the baselines.
baseline_errors: a Python list of NumPy arrays (one per each observation passed to the destriper) containing an optimistic estimate of the error per each baseline. This error is estimated assuming that there is no correlation between baselines, which is only a rough approximation.
baseline_lengths: a Python list of NumPy arrays (one per each observation passed to the destriper) containing the number of TOD samples per each baseline.
stopping_factor: the maximum residual for the destriping solution stored in the field baselines. It is an assessment of the quality of the solution found by the destriper: the lower, the better.
history_of_stopping_factors: list of stopping factors computed by the iterative algorithm. This list should ideally be monothonically decreasing, as this means that the destriping algorithm was able to converge to the correct solution.
converged: a Boolean flag telling whether the destriper was able to achieve the desired accuracy (the value of params.threshold).
elapsed_time_s: the elapsed time spent by the function make_destriped_map().
bytes_in_temporary_buffers: the number of bytes allocated internally by the destriper for temporary buffers.

baseline_errors: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | ndarray[tuple[Any, ...], dtype[_ScalarT]] | None#

baseline_lengths: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | ndarray[tuple[Any, ...], dtype[_ScalarT]] | None#

baselines: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]] | ndarray[tuple[Any, ...], dtype[_ScalarT]] | None#

binned_map: ndarray[tuple[Any, ...], dtype[_ScalarT]]#

bytes_in_temporary_buffers: int#

components: list[str]#

converged: bool | str#

coordinate_system: CoordinateSystem#

destriped_map: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None#

detector_split: str | None = 'full'#

elapsed_time_s: float#

history_of_stopping_factors: list[float]#

hit_map: ndarray[tuple[Any, ...], dtype[_ScalarT]]#

nobs_matrix_cholesky: NobsMatrix#

nside: int#

params: DestriperParameters#

stopping_factor: float | None#

time_split: str | None = 'full'#

class litebird_sim.mapmaking.ExternalDestriperParameters(nside: int = 512, coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic, nnz: int = 3, baseline_length_s: float = 60.0, iter_max: int = 100, output_file_prefix: str = 'lbs_', return_hit_map: bool = False, return_binned_map: bool = False, return_destriped_map: bool = True, return_npp: bool = False, return_invnpp: bool = False, return_rcond: bool = False)#

Bases: object

Parameters used by the Madam mapmakers to produce a map.

The list of fields in this dataclass is the following:

nside: the NSIDE parameter used to create the maps
coordinate_system: an instance of the CoordinateSystem enum. It specifies if the map must be created in ecliptic (default) or galactic coordinates.
nnz: number of components per pixel. The default is 3 (I/Q/U).
baseline_length_s: length of the baseline for 1/f noise in seconds. The default is 60.0 s.
iter_max: maximum number of iterations. The default is 100
output_file_prefix: prefix to be used for the filenames of the Healpix FITS maps saved in the output directory. The default is lbs_.

The following Boolean flags specify which maps should be returned by the function destripe():

return_hit_map: return the hit map (number of hits per pixel)
return_binned_map: return the binned map (i.e., the map with no baselines removed).
return_destriped_map: return the destriped map. If pure white noise is present in the timelines, this should be the same as the binned map.
return_npp: return the map of the white noise covariance per pixel. It contains the following fields: II, IQ, IU, QQ, QU, and UU (in this order).
return_invnpp: return the map of the inverse covariance per pixel. It contains the following fields: II, IQ, IU, QQ, QU, and UU (in this order).
return_rcond: return the map of condition numbers.

The default is to only return the destriped map.

baseline_length_s: float = 60.0#

coordinate_system: CoordinateSystem = 1#

iter_max: int = 100#

nnz: int = 3#

nside: int = 512#

output_file_prefix: str = 'lbs_'#

return_binned_map: bool = False#

return_destriped_map: bool = True#

return_hit_map: bool = False#

return_invnpp: bool = False#

return_npp: bool = False#

return_rcond: bool = False#

class litebird_sim.mapmaking.HMapsResult(h_maps: dict[str, dict[tuple[int, int], h_map_Re_and_Im]], duration_s: float, sampling_rate_Hz: float, coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic, detector_split: str = 'full', time_split: str = 'full')#

Bases: object

Result of make_h_maps().

Variables:

h_maps – Mapping detector -> {(n, m): h_map_Re_and_Im} containing all computed \(h_{n,m}\) maps.
duration_s – Duration of the simulated observation in seconds.
sampling_rate_Hz – Sampling rate in hertz.
coordinate_system – Coordinate system of the output maps.
detector_split – Detector split used to build the maps.
time_split – Time split used to build the maps.

coordinate_system: CoordinateSystem = 1#

detector_split: str = 'full'#

duration_s: float#

h_maps: dict[str, dict[tuple[int, int], h_map_Re_and_Im]]#

sampling_rate_Hz: float#

time_split: str = 'full'#

class litebird_sim.mapmaking.PairDifferencingResult(binned_map: Any = None, invnpp: Any = None, coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic, components: list | None = None, detector_split: str = 'full', time_split: str = 'full')#

Bases: object

Result of a call to the make_pair_differenced_map() function

This dataclass has the following fields:

binned_map: Healpix map containing the binned value for each pixel

invnpp: inverse of the covariance matrix element for each
pixel in the map. It is an array of shape (12 * nside * nside, 2, 2)

coordinate_system: the coordinate system of the output maps (a CoordinateSistem object)
components: list of components included in the map, by default only the field tod is used
detector_split: detector split of the binned map
time_split: time split of the binned map

binned_map: Any = None#

components: list | None = None#

coordinate_system: CoordinateSystem = 1#

detector_split: str = 'full'#

invnpp: Any = None#

time_split: str = 'full'#

litebird_sim.mapmaking.check_valid_splits(observations: Observation | list[Observation], detector_splits: str | list[str] = 'full', time_splits: str | list[str] = 'full')#

Check if the splits are valid

For each observation in the list, check if the detector and time splits are valid. In particular, the compatibility between the detectors in each observation and the desired split in detector domain is checked. On the other hand, this assess whether the desired time split fits inside the duration of the observation (when this applies). If the splits are not compatible with the input data, an error is raised.

Parameters:

observations (list of Observations) –
observations to be mapped. They are required to have the following attributes as arrays
- pointings: the pointing information (in radians) for each tod
  sample. It must be a tensor with shape (N_d, N_t, 3), with N_d number of detectors and N_t number of samples in the TOD.
- any attribute listed in components (by default, tod) and containing the TOD(s) to be binned together.
If the observations are distributed over some communicator(s), they must share the same group processes. If pointings and psi are not included in the observations, they can be provided through an array (or a list of arrays) of dimension (Ndetectors x Nsamples x 3), containing theta, phi and psi
detector_splits (str | list[str], optional) –
detector-domain splits used to produce maps.
- ”full”: every detector in the observation will be used;
- ”waferXXX”: the mapmaking will be performed on the intersection
  of the detectors specified in the input and the detectors specified in the detector_split. The wafer must be specified in the format “waferXXX”. The valid values for “XXX” are all the 3-digits strings corresponding to the wafers in the LITEBIRD focal plane (e.g. L00, M01, H02).
time_splits (str | list[str], optional) –
time-domain splits used to produce maps. This defaults to “full” indicating that every sample in the observation will be used. In addition, the user can specify a string, or a list of strings, to indicate a subsample of the observation to be used:
- ”full”: every sample in the observation will be used;
- ”first_half” and/or “second_half”: the first and/or second half of the
  observation will be used;
- ”odd” and/or “even”: the odd and/or even samples in the observation
  will be used;
- ”yearX”: the samples in the observation will be
  used according to the year they belong to (relative to the starting time). The valid values for “X” are [“1”, “2”, “3”].
- ”surveyX”: the samples in the observation will be used according
  to the requested survey. In this context, a survey is taken to be complete in 6 months, thus the valid values for “X” are [“1”, “2”, “3”, “4”, “5”, “6”].

litebird_sim.mapmaking.destriper_log_callback(stopping_factor: float, step_number: int, max_steps: int) → None#

A function called by make_destriped_map() for every CG iteration

This function has the purpose to produce a visual feedback while the destriper is running the (potentially long) Conjugated Gradient iteration.

litebird_sim.mapmaking.load_destriper_results(folder: Path, custom_dest_file: str | None = None, custom_base_file: str | None = None) → DestriperResult#

Load the results of a call to make_destriped_map() from disk

This function complements save_destriper_results(). It re-creates an object of type DestriperResult from a set of FITS files. If you are calling this function from multiple MPI processes, you must ensure that it is called by every MPI process at the same time, and that the number of MPI processes is the same that was used when the data were saved. (The function checks for this and halts its execution if the condition is not satisfied.)

Parameters:: folder – The folder containing the FITS files to load
Returns:: A new DestriperResult object

litebird_sim.mapmaking.load_h_maps_from_file(filepath: str) → HMapsResult#

Load \(h_{n,m}\) maps from an HDF5 file.

Parameters:: filepath (str) – Path to an HDF5 file produced by save_hn_maps().
Returns:: Loaded maps and metadata.
Return type:: HnMapResult

litebird_sim.mapmaking.make_binned_map(nside: int, observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], pointings: ~numpy.ndarray | list[~numpy.ndarray] | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, output_coordinate_system: ~litebird_sim.coordinates.CoordinateSystem = CoordinateSystem.Galactic, components: str | list[str] = 'tod', detector_split: str = 'full', time_split: str = 'full', pointings_dtype=<class 'numpy.float64'>, nthreads: int | None = None) → BinnerResult#

Bin Map-maker

Map a list of observations

Parameters:

observations (list of Observations) –
observations to be mapped. They are required to have the following attributes as arrays
- pointings: the pointing information (in radians) for each tod
  sample. It must be a tensor with shape (N_d, N_t, 3), with N_d number of detectors and N_t number of samples in the TOD.
- any attribute listed in components (by default, tod) and containing the TOD(s) to be binned together.
If the observations are distributed over some communicator(s), they must share the same group processes. If pointings and psi are not included in the observations, they can be provided through an array (or a list of arrays) of dimension (Ndetectors x Nsamples x 3), containing theta, phi and psi
nside (int) – HEALPix nside of the output map
pointings (array or list of arrays) – optional, external pointing information, if not included in the observations
hwp (HWP, optional) – An instance of the HWP class (optional)
output_coordinate_system (CoordinateSystem) – the coordinates to use for the output map
components (list[str] or str) – components to include in the map-making. The default is just to use the field tod of each Observation object
detector_split (str) – select the detector split to use in the map-making
time_split (str) – select the time split to use in the map-making.
pointings_dtype (dtype) – data type for pointings generated on the fly. If the pointing is passed or already precomputed this parameter is ineffective. Default is np.float64.
nthreads – int, default=None Number of threads to use in ducc’s Healpix methods. If None, the function reads from the OMP_NUM_THREADS environment variable.

Returns:

An instance of the class MapMakerResult. If the observations are: distributed over MPI Processes, all of them get a copy of the same object.

litebird_sim.mapmaking.make_brahmap_gls_map(nside: int, observations: list | ~typing.Any, pointings: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]] | list[~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]]] | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, components: str | list[str] = 'tod', pointings_flag: ~numpy.ndarray | None = None, inv_noise_cov_operator: brahmap.LBSim_InvNoiseCovLO_UnCorr | None = None, threshold: float = 1e-05, pointings_dtype=<class 'numpy.float64'>, gls_params: brahmap.LBSimGLSParameters | None = None) → brahmap.LBSimGLSResult | tuple[brahmap.LBSimProcessTimeSamples, brahmap.LBSimGLSResult]#

GLS Map-maker using Brahmap.

This function allows the users to do the optimal map-making with BrahMap. Since BrahMap is seemlessly interfaced with LBS, it allows the map-making without storing the TODs on the disk. The function needs an inverse noise covariance operator and an object containing the GLS parameters as arguments.

BrahMap offers a variety of noise covariance opeartors, all compatible with LBS. Noise covariance operator for white noise, for example, can be defined as:

inv_cov = brahmap.LBSim_InvNoiseCovLO_UnCorr(
    obs = ...,
    inverse_noise_variance = ...,
    dtype = ...,
)

The parameters used for GLS can be defined as:

gls_params = brahmap.LBSimGLSParameters(
    output_coordinate_system = lbs.CoordinateSystem.Galactic,
    return_processed_samples = False,
    solver_type = brahmap.SolverType.IQU,
    use_iterative_solver = True,
    isolver_threshold = 1.0e-25,
    isolver_max_iterations = 100,
    return_hit_map = False,
)

For the complete list of available noise covariance operators and advance usage, please refer to the BrahMap documentation: https://anand-avinash.github.io/BrahMap/

Parameters:

nside (int) – Nside of the output map
observations (Observation) – An instance of Observation class, or a list of objects of this kind
pointings (np.ndarray, optional) – A 3×N array containing the values of the θ,φ,ψ angles (in radians), or a list if observations was a list. If no pointings are specified, they will be taken from observations (the most common situation)
hwp (HWP, optional) – An instance of HWP class
components (str or list[str], optional) – The components to be summed for map-making, by default “tod”
pointings_flag (np.ndarray, optional) – An array of flags to include or exclude the individual pointing sample, by default None. A False flag excludes the pointing sample whereas a True flag includes it
inv_noise_cov_operator (optional) – Inverse noise covariance operator, by default None
threshold (float, optional) – Threshold parameter to determine the poorly observed pixels, by default 1.0e-5
pointings_dtype (dtype, optional) – dtype to be used for computing pointings on the fly. The pointing_dtype must be same as the dtype of the TOD; by default np.float64
gls_params (LBSimProcessTimeSamples, optional) – An object that encapsulates the parameter of the GLS, including PCG threshold and max iteration; by default None

Returns:

Returns an LBSimGLSResult object when gls_params.return_processed_samples = False. LBSimGLSResult object encapsulates the output of GLS including the output maps and PCG convergence status. The function returns the tuple object (LBSimProcessTimeSamples, LBSimGLSResult) when gls_params.return_processed_samples = True.

Return type:

LBSimGLSResult | tuple[LBSimProcessTimeSamples, LBSimGLSResult]

Raises:

ImportError – Raises ImportError when the brahmap package couldn’t be imported

litebird_sim.mapmaking.make_destriped_map(nside: int, observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], pointings: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]] | list[~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]]] | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, params: ~litebird_sim.mapmaking.destriper.DestriperParameters = DestriperParameters(output_coordinate_system=<CoordinateSystem.Galactic: 2>, samples_per_baseline=100, iter_max=100, threshold=1e-07, use_preconditioner=True), components: str | list[str] = 'tod', keep_weights: bool = False, keep_pixel_idx: bool = False, keep_pol_angle_rad: bool = False, detector_split: str = 'full', time_split: str = 'full', baselines_list: list[~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]]] | None = None, recycled_convergence: bool = False, callback: ~typing.Any = <function destriper_log_callback>, callback_kwargs: dict[~typing.Any, ~typing.Any] | None = None, pointings_dtype=<class 'numpy.float64'>, nthreads: int | None = None) → DestriperResult#

Applies the destriping algorithm to produce a map out from a TOD

This function takes the samples stored in the list of Observation objects observations and produces a destriped map. Pointings can either be embedded in the observations objects or provided through the parameter pointings.

The params argument is an instance of the class DestriperParameters, and it is used to tune the way the destriper should solve the map-making equation. Have a look at the documentation for this class to get more information.

The samples in the TOD can be saved in several fields within each Observation object. For instance, you could have generated a noise component in observations.noise_tod, the dipole in observations.dipole_tod, etc., and you want to produce a map containing the sum of all these components. In this case, you can pass a list of strings containing the name of the fields as the parameter components, and the destriper will add them together before creating the map.

To show a few clues about the progress of the destriping algorithm, the function accepts the callback argument. This must be a function with the following prototype:

def callback(
    stopping_factor: float, step_number: int, max_steps: int
):
    ...

The parameter stopping_factor is the current stopping factor, i.e., the maximum residual between the current estimate of the baselines and the result of applying the binning to the destriped TOD. The parameter step_number is an integer between 0 and max_steps - 1 and it is increased by 1 each time the callback is called. The callback function can accept more parameters; in this case, their value must be passed using the callback_kwargs parameter. For instance, you might be running the destriper within a graphical program and you would show a progress bar dialog; in this case, one of the additional arguments might be an handle to the progress window:

def my_gui_callback(
    stopping_factor: float, step_number: int, max_steps: int,
    window_handle: GuiWindow,
):
    window_handle.set_progress(step_number / (max_steps - 1))
    window_handle.text_label.set(
        f"Stopping factor is {stopping_factor:.2e}"
    )

To use it, you must call make_destriped_map as follows:

my_window_handle = CreateProgressWindow(...)
make_destriped_map(
    observations=obs_list,
    params=params,
    callback=my_gui_callback,
    callback_kwargs={"window_handle": my_window_handle},
)

Parameters:

observations – an instance of the class Observation, or a list of objects of this kind
pointings – a 3×N array containing the values of the θ,φ,ψ angles (in radians), or a list if observations was a list. If no pointings are specified, they will be taken from observations (the most common situation)
hwp – An instance of the HWP class (optional)
params – an instance of the DestriperParameters class
components – components to extract from the TOD and sum together. The default is to use observations.tod.
keep_weights – the destriper adds a destriper_weights field to each Observation object in observations, and it deletes it once the map has been produced. Setting this parameter to True prevents the field from being deleted. (Useful for debugging.)
keep_pixel_idx – same as keep_weights, but the field to be kept from the Observation classes is destriper_pixel_idx
keep_pol_angle_rad – same as keep_weights, but the field to be kept from the Observation classes is destriper_pol_angle_rad
callback – a function that is called during each iteration of the CG routine. It is meant as a way to provide some visual feedback to the user. The default implementation uses the logging library to print a line. You can provide a custom callback, if you want (see above for an example).
callback_kwargs – additional keyword arguments to be passed to the callback function
pointings_dtype(dtype) – data type for pointings generated on the fly. If the pointing is passed or already precomputed this parameter is ineffective. Default is np.float64.
nthreads – number of threads to use in ducc’s Healpix methods. If None, the function reads from the OMP_NUM_THREADS environment variable.

Returns:

an instance of the DestriperResult containing the destriped map and other useful information

litebird_sim.mapmaking.make_h_maps(observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], nside: int, pointings: ~numpy.ndarray | list[~numpy.ndarray] | None = None, n_m_couples: ~numpy.ndarray = array([[0, 0], [2, 0], [4, 0]]), hwp: ~litebird_sim.hwp.HWP | None = None, output_coordinate_system: ~litebird_sim.coordinates.CoordinateSystem = CoordinateSystem.Galactic, detector_split: str = 'full', time_split: str = 'full', pointings_dtype=<class 'numpy.float64'>, save_to_file: bool = True, output_directory: str = './h_n_maps') → HMapsResult#

Generate complex harmonic maps \(h_{n,m}\) from observations.

The map for (n, m) = (0, 0) contains hit counts per pixel.

Parameters:

observations (Observation | list[Observation]) – Observation object or list of observations used to build the maps.
nside (int) – HEALPix nside of the output maps.
pointings (np.ndarray | list[np.ndarray] | None) – Optional external pointings. Expected shape is (n_detectors, n_samples, 3) per observation.
n_m_couples (np.ndarray) – Array of (n, m) pairs with shape (N, 2).
hwp (HWP | None) – Half-wave plate model.
output_coordinate_system (CoordinateSystem) – Coordinate system of the output maps.
detector_split (str) – Detector split to apply.
time_split (str) – Time split to apply.
pointings_dtype (type) – Data type used when pointings are computed on the fly.
save_to_file (bool) – If True, save maps to HDF5 files.
output_directory (str) – Output directory for generated HDF5 files.

Returns:

Result container with all computed maps and metadata.

Return type:

HnMapResult

litebird_sim.mapmaking.make_pair_differenced_map(nside: int, observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], pointings: ~numpy.ndarray | list[~numpy.ndarray] | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, output_coordinate_system: ~litebird_sim.coordinates.CoordinateSystem = CoordinateSystem.Galactic, components: str | list[str] = 'tod', detector_split: str = 'full', time_split: str = 'full', pointings_dtype=<class 'numpy.float64'>, nthreads: int | None = None) → PairDifferencingResult#

QU pair-differencing map-maker

Map a list of observations

Parameters:

observations (list of Observations) –
observations to be mapped. They are required to have the following attributes as arrays
- pointings: the pointing information (in radians) for each tod
  sample. It must be a tensor with shape (N_d, N_t, 3), with N_d number of detectors and N_t number of samples in the TOD.
  any attribute listed in components (by default, tod) and
  containing the TOD(s) to be differenced and binned together.
  
  The local detector set is required to be arranged in T/B pairs sharing the same wafer and focal-plane pixel. The map-maker differences the TOD of each pair and solves only for the Q/U Stokes components.
If the observations are distributed over some communicator(s), they must share the same group processes. If pointings and psi are not included in the observations, they can be provided through an array (or a list of arrays) of dimension (Ndetectors x Nsamples x 3), containing theta, phi and psi
nside (int) – HEALPix nside of the output map
pointings (array or list of arrays) – optional, external pointing information, if not included in the observations
hwp (HWP, optional) – An instance of the HWP class (optional)
output_coordinate_system (CoordinateSystem) – the coordinates to use for the output map
components (list[str] or str) – components to include in the map-making. The default is just to use the field tod of each Observation object
detector_split (str) – select the detector split to use in the map-making
time_split (str) – select the time split to use in the map-making.
pointings_dtype (dtype) – data type for pointings generated on the fly. If the pointing is passed or already precomputed this parameter is ineffective. Default is np.float64.
nthreads – int, default=None Number of threads to use in ducc’s Healpix methods. If None, the function reads from the OMP_NUM_THREADS environment variable.

Returns:

An instance of the class PairDifferencingResult. If the observations are: distributed over MPI Processes, all of them get a copy of the same object.

litebird_sim.mapmaking.remove_baselines_from_tod(obs_list: list[Observation], baselines: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]], baseline_lengths: list[ndarray[tuple[Any, ...], dtype[_ScalarT]]], component: str = 'tod') → None#

Subtract 1/f baselines from the TODs in obs_list

This functions subtracts the baselines in baselines from the component in the Observation objects listed within the argument obs_list. The two lists baselines and baseline_lengths must have the same number of elements as obs_list.

Parameters:

obs_list – A list of Observation objects containing the TOD to be cleaned from baselines
baselines – A list of NumPy arrays of floats (one per each item in obs_list): each value is a baseline to subtract
baseline_lengths – A list of NumPy arrays of integers (one per each item in obs_list). Each array contains the number of samples within each baseline in the argument baselines; the sum of all the values within each array must be equal to the number of samples in the TOD within the corresponding element in obs_list.
component – The name of the TOD component to clean from the baselines. The default is "tod".

litebird_sim.mapmaking.remove_destriper_baselines_from_tod(obs_list: list[Observation], destriper_result: DestriperResult, component: str = 'tod')#

A wrapper around remove_baselines_from_tod()

This method removes the baselines computed by make_destriped_map() from the TOD in the observations obs_list. The TOD component is specified by the string component, which is tod by default.

Parameters:

obs_list – A list of Observation objects containing the TOD to be cleaned from baselines
destriper_result – The result of a call to make_destriped_map()
component – The name of the TOD component to clean. The default is "tod".

litebird_sim.mapmaking.save_destriper_results(results: DestriperResult, output_folder: Path, custom_dest_file: str | None = None, custom_base_file: str | None = None) → None#

Save the results of a call to make_destriped_map() to disk

The results are saved in a set of FITS files:

A FITS file containing the maps (destriped, binned, hits), the N_obs matrix, and general information about the convergence
A set of FITS files containing the baselines. Each MPI process writes one file containing its baselines.

The only parameter that is not saved is the field results.params.samples_per_baseline: its type is very versatile and would not fit well in the FITS file format. Moreover, the necessary information is already available in the baseline lengths that are saved in the files (see point 2. above).

To load the results from the files, use load_destriper_results().

Parameters:

results – The result of the call to make_destriped_map()
output_folder – The folder where to save all the FITS files

litebird_sim.madam.ensure_parent_dir_exists(file_name: str | Path)#

litebird_sim.madam.save_simulation_for_madam(sim: ~litebird_sim.simulations.Simulation, params: ~litebird_sim.mapmaking.common.ExternalDestriperParameters, detectors: list[~litebird_sim.detectors.DetectorInfo] | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, pointings: ~numpy.ndarray | list[~numpy.ndarray] | None = None, use_gzip: bool = False, output_path: str | ~pathlib.Path | None = None, absolute_paths: bool = True, madam_subfolder_name: str = 'madam', components: list[str] = ['tod'], components_to_bin: list[str] | None = None, pointing_dtype=<class 'numpy.float64'>, output_coordinate_system: ~litebird_sim.coordinates.CoordinateSystem = CoordinateSystem.Ecliptic, save_pointings: bool = True, save_tods: bool = True) → dict[str, Any] | None#

Save the TODs and pointings of a simulation to files suitable to be read by Madam

This function takes all the TOD samples and pointing angles from sim and saves them to the directory specified by output_path (the default is to save them in a sub-folder of the output path of the simulation). The parameter detector must be a list of DetectorInfo objects, and it specifies which detectors will be saved to disk; if it is None, all the detectors in the simulation will be considered. The variable params specifies how Madam should produce the maps; see the documentation for ExternalDestriperParameters for more information.

If use_gzip is true, the TOD and pointing files will be compressed using Gzip (the default is false, as this might slow down I/O). If absolute_paths is True (the default), the parameter and simulation files produced by this routine will be absolute; set it to False if you plan to move the FITS files to some other directory or computer before running Madam.

The parameter madam_subfolder_name is the name of the directory within the output folder of the simulation that will contain the Madam parameter files.

You can use multiple TODs in the map-making process. By default, the code will only dump Observation.tod in the FITS files, but you can specify additional components via the components parameter, which is a list of the fields that must be saved in the FITS files and included in the parameter and simulation files. All these components will be summed in the map-making process.

If you want to create a map using just a subset of the components listed in the components parameter, list them in components_to_bin. (This is usually employed when you pass save_pointings=False and save_tods=False.) If components_to_bin is None (the default), all the elements in components will be used.

If you are using MPI, call this function on all the MPI processes, not just on the one with rank #0.

The flags save_pointings and save_tods are used to tell if you want pointings and TODs to be saved in FITS files or not. If either flag is set to false, the corresponding FITS files will not be produced, but the .sim file for Madam will nevertheless list them as if they were created. This is useful if you plan to reuse files from some other call to save_simulation_for_madam; in this case, a common trick is to create soft links to them in the output directory where the .par and .sim files are saved.

pointings and hwp are the optional parameters. External pointing information, if not included in the observations, must be passed through the pointings parameter. It is assumed that the pointing information is available in ecliptic coordinates. The pointings are therefore as such. To save pointings in other coordinates, the parameter output_coordinate_system can be used. The HWP object should be passed to hwp parameter in order to compute the HWP angles.

When pointings are computed on the fly, they are computed in double precision. It can be modified with the argument pointing_dtype.

The return value is either a dictionary containing all the parameters used to fill Madam files (the parameter file and the simulation file) or None; the dictionary is only returned for the MPI process with rank #0.