Map-making¶

The primary aim of the LiteBIRD Simulation Framework is to create synthetic timestreams as if they were acquired by the real LiteBIRD Instrument. The process that creates 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 destriper, i.e., a more advanced map-maker that can remove the effect of correlated instrumental noise from the timelines before producing a map. The correlated noise is usually referred 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 ran over the cleaned timelines.
A wrapper that enables you to use the TOAST2 destriper. (To use this, you must ensure that the package toast is installed.)
You can also use save_simulation_for_madam() to save TODs and pointing information to disk and then manually call the Madam mapmaker.

In this chapter, we assume that you have already created the timelines that must be provided 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)
        ),
    ],
    dtype_tod=np.float64,  # Needed if you use the TOAST destriper
    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 theory of binners and destripers is explained in the paper [5], and the 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: a single (or a list of) Observation, the Healpix resolution of the output map (nside) and produces a coadded map. It assumes white noise and each detector gets weighted by \(1 / 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, obs=obs)
healpy.mollview(result.binned_map)

(The pointing information is included in the Observation, alternatively pointings can be provided 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.

The 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)

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

It is now worth explaining how the binned map is calculated, and we will do it using a simple example. Suppose that we have measured just two pixels in the map with only seven measurements. (These numbers are ridiculous, but in this way the matrices we are going to 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 index of the pixel 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 two pixels are observed by the seven TOD samples. 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.

It is important that each pixel be covered by enough samples: to recover the three Stokes parameters associated with it (I, Q, and U), there must be at least three non-degenerate polarization angles. KS2009 describes how to recover this information 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 is performed.
"first_half": the first half of the TOD is used.
"second_half": the second half of the TOD is used.
"odd": the odd samples of the TOD are used.
"even": the even samples of the TOD are used.
"yearX": the TOD samples relative to the X-th year of observation is used. The X-th part must be an integer between 1 and 3.
"surveyX": assuming that a complete survey is performed in 6 months, the X-th survey is used. The X-th part must be an integer between 1 and 6.

As regards the detector split, the supported options are:

"full" (default): no split is performed.
"waferXXX": the detectors in the wafer XXX are used. The XXX part must specify the desired 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 duration of the observation and with the detector list. Thus, for example, if the observation lasts 1 year, the split “year2” will raise an AsserionError. Similarly, if the observation is performed with some detector contained in the L00 wafer, the split “waferL03” will also raise an AsserionError.

Destriper¶

If you know that your simulation contains 1/f noise, you should avoid using the binner and use a better map-making algorithm, as 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(), which is 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, obs=obs, params=params)
healpy.mollview(result.destriped_map)

The result is an instance of the class DestriperResult, which is similar to BinnerResult but it contains much more information.

The 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 is the meaning of each parameter in the classes DestriperParameters and DestriperResult. 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 are discussing 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 has shown to work 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 show a concrete example about how the destriper works, we will resume the toy TOD containing 7 samples that was discussed in section Binner. We will group these samples in 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 true 1/f noise, let’s pretend that what we actually 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 visited 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 1/f noise is taken into account, the estimated I/Q/U parameter for this pixel must agree. This problem can be rewritten 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 actually assigning the first four samples to the first baseline and the last three samples to the last baseline. The application of \(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 quite alarming, and it might sound like it is 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 the baselines are subtracted, and the consistency is the same if all the baselines are shifted by the same amount!) What we are looking for is 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 be such 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 kind of 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 obs. Note that if you provide an integer, it might be that not all baselines will have exactly that length: it depends whether the number \(N_t\) of samples in the TOD is evenly divisibile 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 obs. 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 that 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. The baselines for detector A are shown using continuous lines, while the lines for B are dashed. Note that the number of baselines in the first Observation object (5) is different than 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: it is 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 are 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 a solution for the destriping equation has been computed, 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.

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}\).

TOAST2 Destriper¶

If you install the toast using pip, you can use the TOAST2 destriper within the LiteBIRD Simulation Framework. As TOAST is an optional dependency, you should check if the Framework was able to detect its presence:

import litebird_sim as lbs

if lbs.TOAST_ENABLED:
    # It's ok to use TOAST
    ...
else:
    # TOAST is not present, do something else
    ...

The procedure to use the TOAST2 destriper is similar to the steps required to call the internal destriper: you must create a ExternalDestriperParameters object that specifies which input parameters (apart from the timelines) should be used:

params = lbs.ExternalDestriperParameters(
    nside=16,
    return_hit_map=True,
    return_binned_map=True,
    return_destriped_map=True,
)

The parameters we use here are the resolution of the output map (nside=16), and the kind of results that must be returned: specifically, we are looking here for the hit map (i.e., a map that specifies how many samples were observed while the detector was looking at a specific pixel), the binned map (the same map that would be produced by the binner, see above), and the destriped map.

Note

The TOAST destriper only works with timelines containing 64-bit floating point numbers. As the default data type for timelines created by sim.create_observations is a 32-bit float, if you plan to run the destriper you should pass the flag dtype_tod=np.float64 to sim.create_observations (see the code above), otherwise destripe will create an internal copy of the TOD converted in 64-bit floating-point numbers, which is usually a waste of space.

To run the TOAST2 destriper, you simply call destripe_with_toast2():

result = lbs.destripe_with_toast2(sim, params)

The result is an instance of the class Toast2DestriperResult and contains the three maps we have asked above (hit map, binned map, destriped map).

Note

Unlike the internal destriper, TOAST2 does not return information about the convergence of the CG algorithm, and it is not granted that the norm used to estimate the stopping factor is the same as the one calculated by the internal destriper. Therefore, please avoid comparing the stopping factors of the two destripers!

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.

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 they are not going to be used in the first and second map. The reason is that this parameter is used by the function to create a «map» of the components as they are supposed to be found 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 actually used to create the map, so it will not be confused if the TOD FITS files will 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. No information about the map-making process is included in a simulation file: it just tells how many FITS files should be read and what is inside each of them.

The parameter file is used to tell Madam how you want maps to be created. It’s in the parameter file that you can ask Madam to skip parts of the TODs, for example because 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().

But 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 = None, invnpp: Any | None = None, coordinate_system: CoordinateSystem = CoordinateSystem.Ecliptic, components: List | None = None, detector_split: str | None = None, time_split: str | None = None)¶

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¶

coordinate_system: CoordinateSystem = 1¶

detector_split: str = None¶

invnpp: Any = None¶

time_split: str = None¶

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[_SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]] | None = 100¶

threshold: float = 1e-07¶

use_preconditioner: bool = True¶

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: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None¶

baseline_lengths: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None¶

baselines: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None¶

binned_map: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]¶

bytes_in_temporary_buffers: int¶

components: List[str]¶

converged: bool¶

coordinate_system: CoordinateSystem¶

destriped_map: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes] | None¶

elapsed_time_s: float¶

history_of_stopping_factors: List[float] | None¶

hit_map: _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]¶

nobs_matrix_cholesky: NobsMatrix¶

nside: int¶

params: DestriperParameters¶

stopping_factor: float | None¶

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 TOAST/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¶

litebird_sim.mapmaking.check_valid_splits(obs: 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:

obs (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 (Union[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 (Union[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) → 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.make_binned_map(nside: int, obs: Observation | List[Observation], pointings: ndarray | List[ndarray] | None = None, output_coordinate_system: CoordinateSystem = CoordinateSystem.Galactic, components: List[str] | None = None, detector_split: str = 'full', time_split: str = 'full') → BinnerResult¶

Bin Map-maker

Map a list of observations

Parameters:

obs (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
output_coordinate_system (CoordinateSystem) – the coordinates to use for the output map
components (list[str]) – list of 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.

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.

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

This function takes the samples stored in the list of Observation objects obs and produces a destriped map. Pointings can either be embedded in the obs 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 obs.noise_tod, the dipole in obs.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(
    obs=obs_list,
    params=params,
    callback=my_gui_callback,
    callback_kwargs={"window_handle": my_window_handle},
)

Parameters:

obs – 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 obs was a list. If no pointings are specified, they will be taken from obs (the most common situation)
params – an instance of the DestriperParameters class
components – a list of components to extract from the TOD and sum together. The default is to use obs.tod.
keep_weights – the destriper adds a destriper_weights field to each Observation object in obs, 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

Returns:

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

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) → 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: Simulation, params: ExternalDestriperParameters, detectors: List[DetectorInfo] | None = None, use_gzip: bool = False, output_path: str | Path | None = None, absolute_paths: bool = True, madam_subfolder_name: str = 'madam', components: List[str] = ['tod'], components_to_bin: List[str] | None = None, 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.

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.

Map-making¶

Binner¶

Data splits¶

Destriper¶

TOAST2 Destriper¶

Saving files for Madam¶

API reference¶

LiteBIRD simulation pipeline

Navigation

Related Topics