Observations¶

The Observation class, is the container for the data acquired by the telescope during a scaning period (and the relevant information about it).

Serial applications¶

In a serial code Observation is equivalent to an empty class in which you can put anything:

import litebird_sim as lbs
import numpy as np

obs = lbs.Observation(
    detectors=2,
    start_time=0.0,
    sampling_rate_hz=5.0,
    n_samples=5,
)

obs.my_new_attr = 'value'
setattr(obs, 'another_new_attr', 'another value')

Across the framework, the coherence in the names and content of the attributes is guaranteed by convention (no check is done by the Observation class).

DO WE NEED A LIST OF THEM?

There are two more conventions. First, obs.tod is a 2-D array. obs.tod[0] is the time stream of the first detector. obs.tod[:, 0] are the first time sample of all the detectors. Second, detector-specific attributes are collected in arrays: obs.calibration_factors is a 1-D array with one entry per detector.

Thus, typical operations look like this:

# Collect detector properties in arrays
obs.calibration_factors = np.array([1.1, 1.2])
obs.wn_levels = np.array([2.1, 2.2])

# Apply to each detector its own calibration factor
obs.tod *=  obs.calibration_factors[:, None]

# Add white noise at a different level for each detector
obs.tod += (np.random.normal(size=obs.tod.shape)
            * obs.wn_levels[:, None])

Parallel applications¶

The only work that the .Observation class actually does is handlying parallelism. obs.tod can be distributed over a n_blocks_det by n_blocks_time grid of MPI ranks. The blocks can be changed at run-time.

The coherence between the serial and parallel operations is achieved by distributing also the arrays of detector properties: each rank keep in memory only an interval of calibration_factors or wn_levels, the same detector interval of tod.

The main advantage is that the example operations in the Serial section are achieved with the same lines of code. The price to pay is that you have to set detector properties with special methods.

import litebird_sim as lbs
from mpi4py import MPI

comm = MPI.COMM_WORLD

obs = lbs.Observation(
    detectors=2,
    start_time=0.0,
    sampling_rate_hz=5.0,
    n_samples=5,
    n_blocks_det=2, # Split the detector axis in 2
    comm=comm  # across the processes of this communicator
)

# Add detector properties either with a global array that contains
# all the detectors (which will be scattered across the processor grid)
obs.setattr_det_global('calibration_factors', np.array([1.1, 1.2]))

# Or with the local array, if somehow you have it already
if comm.rank == 0:
    wn_level_local = np.array([2.1])
elif comm.rank == 1:
    wn_level_local = np.array([2.2])
else:
    wn_level_local = np.array([])
obs.setattr_det('wn_levels', wn_level_local)

# Operate on the local portion of the data just like in serial code

# Apply to each detector its own calibration factor
obs.tod *=  obs.calibration_factors[:, None]

# Add white noise at a different level for each detector
obs.tod += (np.random.normal(size=obs.tod.shape)
            * obs.wn_levels[:, None])

# Change the data distribution
obs.set_blocks(n_blocks_det=1, n_blocks_time=1)
# Now the rank 0 has exactly the data of the serial obs object

For clarity, here is a visualization of how data (a detector attribute and the TOD) gets distributed.

_images/observation_data_distribution.png

When n_blocks_det != 1, keep in mind that obs.tod[0] or obs.wn_levels[0] are quantities of the first local detector, not global. This should not be a problem as the only thing that matter is that the two quantities refer to the same detector. If you need the global detector index, you can get it with obs.det_idx[0], which is create at at construction time.

Other notable functionalities¶

Time starting time can be represented either as floating-point values (appropriate in 99% of the cases) or MJD; in the latter case, it is handled through the AstroPy library.

Instead of adding detector attributes after construction, you can pass a list of dictionaries (one entry for each detectors).

import litebird_sim as lbs
from astropy.time import Time

# Second case: use MJD to track the time
obs_mjd = lbs.Observation(
    detectors=[{"name": "A"}, {"name": "B"}]
    start_time=Time("2020-02-20", format="iso"),
    sampling_rate_hz=5.0,
    nsamples=5,
)

API reference¶

class litebird_sim.observations.Observation(detectors: Union[int, List[dict]], n_samples: int, start_time: Union[float, astropy.time.core.Time], sampling_rate_hz: float, allocate_tod=True, dtype_tod=<class 'numpy.float32'>, n_blocks_det=1, n_blocks_time=1, comm=None, root=0)¶

Bases: object

An observation made by one or multiple detectors over some time window

This class encodes the data acquired by one or multiple detectors over a finite amount of time, and it is the fundamental block of a simulation. The characteristics of the detector are assumed to be stationary over the time span covered by a simulation; these can include:

Noise parameters
Gain

To access the TOD, use one of the following methods:

get_times() returns the array of time values (one per each sample in the TOD)
tod() returns the array of samples

Parameters

detectors (int/list of dict) – Either the number of detectors or a list of dictionaries with one entry for each detector. The keys of the dictionaries will become attributes of the observation. If a detector is missing a key it will be set to nan. If an MPI communicator is passed to comm, detectors is relevant only for the root process
n_samples (int) – The number of samples in this observation.
start_time – Start time of the observation. It can either be a astropy.time.Time type or a floating-point number. In the latter case, it must be expressed in seconds.
sampling_rate_hz (float) – The sampling frequency, in Hertz.
dtype_tod (dtype) – Data type of the TOD array
n_blocks_det (int) – divide the detector axis of the tod in n_blocks_det blocks
n_blocks_time (int) – divide the time axis of the tod in n_blocks_time blocks
comm – either None (do not use MPI) or a MPI communicator object, like mpi4py.MPI.COMM_WORLD. Its size is required to be at least n_blocks_det times n_blocks_time
root (int) – rank of the process receiving the detector list, if detectors is a list of dictionaries, otherwise it is ignored.

detector_global_info(name, info, root=0)¶

Piece of information on the detectors

Variant of detector_info() to be used when the information comes from a single MPI rank (root). In particular,

In the root process, info is required to have n_detectors elements (not self.tod.shape[1]). For other processes info is irrelevant

info is scattered from the root rank to all the processes in self.comm so that self.name will have self.tod.shape[0] elements and self.name[i] is a property of self.tod[i]

When changing n_blocks_det, set_n_blocks() is aware of name and will redistribute info in such a way that self.name[i] is a property of self.tod[i] in the new block distribution

Parameters

name (str) – Name of the information
info (array) – Array containing n_detectors entries. Relevant only for thr root process
root (int) – Rank of the root process

detector_info(name, info)¶

Piece of information on the detectors

Store info as the attribute name of the observation. The difference with respect to self.name = info, relevant only in MPI programs, are

info is assumed to contain a number of elements equal to self.tod.shape[0] and to have the same order (info[i] is a property of self.tod[i])

When changing n_blocks_det, set_n_blocks() is aware of name and will redistribute info in such a way that self.name[i] is a property of self.tod[i] in the new block distribution

Parameters

name (str) – Name of the information
info (array) – Information to be stored in the attribute name. The array must contain one entry for each local detector.

get_det2ecl_quaternions(spin2ecliptic_quats: litebird_sim.scanning.Spin2EclipticQuaternions, detector_quats, bore2spin_quat)¶

Return the detector-to-Ecliptic quaternions

This function returns a (D, N, 4) tensor containing the quaternions that convert a vector in detector’s coordinates into the frame of reference of the Ecliptic. The number of quaternions is equal to the number of samples hold in this observation.

Given that the z axis in the frame of reference of a detector points along the main beam axis, this means that if you use these quaternions to rotate the vector z = [0, 0, 1], you will end up with the sequence of vectors pointing towards the points in the sky (in Ecliptic coordinates) that are observed by the detector.

This is a low-level method; you should usually call the method get_pointings(), which wraps this function to compute both the pointing direction and the polarization angle.

See also the method get_ecl2det_quaternions(), which mirrors this one.

get_ecl2det_quaternions(spin2ecliptic_quats: litebird_sim.scanning.Spin2EclipticQuaternions, detector_quats, bore2spin_quat)¶

Return the Ecliptic-to-detector quaternions

This function returns a (D, N, 4) matrix containing the N quaternions of alla the D detectors that convert a vector in Ecliptic coordinates into the frame of reference of the detector itself. The number of quaternions is equal to the number of samples hold in this observation.

This method is useful when you want to simulate how a point source is observed by the detector’s beam: if you know the Ecliptic coordinates of the point sources, you can easily derive the location of the source with respect to the reference frame of the detector’s beam.

get_pointings(spin2ecliptic_quats: litebird_sim.scanning.Spin2EclipticQuaternions, detector_quats, bore2spin_quat)¶

Return the time stream of pointings for the detector

Given a Spin2EclipticQuaternions and a quaternion representing the transformation from the reference frame of a detector to the boresight reference frame, compute a set of pointings for the detector that encompass the time span covered by this observation (i.e., starting from self.start_time and including self.n_samples pointings).

The parameter spin2ecliptic_quats can be easily retrieved by the field spin2ecliptic_quats in a object of Simulation object, once the method Simulation.generate_spin2ecl_quaternions() is called.

The parameter detector_quats is a stack of detector quaternions. For example, it can be

The stack of the field quat of an instance of the class
DetectorInfo
If all you want to do is a simulation using a boresight
direction, you can pass the value np.array([[0., 0., 0., 1.]]), which represents the null rotation.

The parameter bore2spin_quat is calculated through the class Instrument, which has the field bore2spin_quat. If all you have is the angle β between the boresight and the spin axis, just pass quat_rotation_y(β) here.

The return value is a (D x N × 3) tensor: the colatitude (in radians) is stored in column 0 (e.g., result[:, :, 0]), the longitude (ditto) in column 1, and the polarization angle (ditto) in column 2. You can extract the three vectors using the following idiom:

pointings = obs.get_pointings(...)
# Extract the colatitude (theta), longitude (psi), and
# polarization angle (psi) from pointings
theta, phi, psi = [pointings[:, :, i] for i in (0, 1, 2)]

get_times(normalize=False, astropy_times=False)¶

Return a vector containing the time of each sample in the observation

The measure unit of the result depends on the value of astropy_times: if it’s true, times are returned as astropy.time.Time objects, which can be converted to several units (MJD, seconds, etc.); if astropy_times is false (the default), times are expressed in seconds. In the latter case, you should interpret these times as sidereal.

If normalize=True, then the first time is zero. Setting this flag requires that astropy_times=False.

In any case the size of the returned array is self.tod.shape[1].

This can be a costly operation; you should cache this result if you plan to use it in your code, instead of calling this method over and over again.

Note for MPI-parallel codes: the times returned are only those of the local portion of the data. This means that

the size of the returned array is smaller than n_samples whenever there is more than one time-block

self.tod * self.get_times() is a meaningless but always allowed operation

property n_blocks_det¶

property n_blocks_time¶

property n_detectors¶

property n_samples¶: Samples in the whole observation

Note

If you need the time-lenght of the self.tod array, use directly self.tod.shape[1]: n_samples allways refers to the full observation, even when it is distributed over multiple processes and self.tod is just the local block.

property sampling_rate_hz¶

set_n_blocks(n_blocks_det=1, n_blocks_time=1)¶

Change the number of blocks

Parameters

n_blocks_det (int) – new number of blocks in the detector direction
n_blocks_time (int) – new number of blocks in the time direction

Observations¶

Serial applications¶

Parallel applications¶

Other notable functionalities¶

API reference¶

LiteBIRD simulation pipeline

Navigation

Related Topics