$4\pi$ beam convolution¶

TOAST provides an interface, OpSimConviqt, to the spherical harmonic convolution library, libconviqt. It was developed by Gary Prezeau and Martin Reinecke and described in

G. Prézeau and M. Reinecke:
Algorithm for the Evaluation of Reduced Wigner Matrices,
APJS 190 (2010) 267

arXiv:1002.1050. This particular implementation of the algorithm is available at https://github.com/hpc4cmb/libconviqt.

# Load common tools for all lessons
import sys
sys.path.insert(0, "..")
from lesson_tools import (
    fake_focalplane
)

# Capture C++ output in the jupyter cells
%reload_ext wurlitzer

Method¶

libconviqt takes in spherical harmonic expansions of the beam and the sky and then synthesizes TOD samples at sample positions in the proper orientation. For efficiency, the sky is distributed as isolatitude rings and then each process gets the detector samples that fall on their rings. The calculation itself has two steps, first conviqt builds a 3D interpolator of the beam-convolved sky on a grid of $(\theta, \phi, \psi)$ and then the detector samples are interpolated from the grid. Finally the samples are communited back to the processes that own them.

Typically the interpolation step dominates but if there are few detector samples and the sky and beam expansion orders are high, it is possible that building the interpolator is more expensive.

Example¶

In this section we create a TOAST data object with simulated signal and noise and process the data into hit maps, pixels noise matrices and signal maps.

import toast
import toast.todmap
import toast.pipeline_tools
from toast.mpi import MPI

import numpy as np
import matplotlib.pyplot as plt

mpiworld, procs, rank = toast.mpi.get_world()
comm = toast.mpi.Comm(mpiworld)

# A pipeline would create the args object with argparse

class args:
    sample_rate = 10  # Hz
    hwp_rpm = None
    hwp_step_deg = None
    hwp_step_time_s = None
    spin_period_min = 1 # 10
    spin_angle_deg = 20 # 30
    prec_period_min = 100 # 50
    prec_angle_deg = 30 # 65
    coord = "E"
    nside = 64
    nnz = 3
    outdir = "maps"
    sky_file = "slm.fits"
    beam_file = "blm.fits"

# Create a fake focalplane, we could also load one from file.
# The Focalplane class interprets the focalplane dictionary
# created by fake_focalplane() but it can also load the information
# from file.

focalplane = fake_focalplane(samplerate=args.sample_rate, fknee=0.1, alpha=2)
detectors = sorted(focalplane.keys())
detquats = {}
for d in detectors:
    detquats[d] = focalplane[d]["quat"]
    
nsample = 100000
start_sample = 0
start_time = 0
iobs = 0
    
tod = toast.todmap.TODSatellite(
    comm.comm_group,
    detquats,
    nsample,
    coord=args.coord,
    firstsamp=start_sample,
    firsttime=start_time,
    rate=args.sample_rate,
    spinperiod=args.spin_period_min,
    spinangle=args.spin_angle_deg,
    precperiod=args.prec_period_min,
    precangle=args.prec_angle_deg,
    detranks=comm.group_size,
    hwprpm=args.hwp_rpm,
    hwpstep=args.hwp_step_deg,
    hwpsteptime=args.hwp_step_time_s,
)

# Constantly slewing precession axis                                                                                                                                             
precquat = np.empty(4 * tod.local_samples[1], dtype=np.float64).reshape((-1, 4))
toast.todmap.slew_precession_axis(
    precquat,
    firstsamp=start_sample + tod.local_samples[0],
    samplerate=args.sample_rate,
    degday=360.0 / 365.25,
)
tod.set_prec_axis(qprec=precquat)

noise = toast.pipeline_tools.get_analytic_noise(args, comm, focalplane)

obs = {}
obs["name"] = "science_{:05d}".format(iobs)
obs["tod"] = tod
obs["intervals"] = None
obs["baselines"] = None
obs["noise"] = noise
obs["id"] = iobs

# Conviqt requires at least minimal focal plane information to be present in the observation
obs["focalplane"] = toast.pipeline_tools.Focalplane(focalplane)
"""
for det in tod.local_dets:
    obs["focalplane"][det] = {
        "epsilon" : focalplane[det]["epsilon"],
    }
    if det.endswith("A"):
        obs["focalplane"][det]["psi_pol_deg"] = 0,
    elif det.endswith("B"):
        obs["focalplane"][det]["psi_pol_deg"] = 90,
"""


data = toast.Data(comm)
data.obs.append(obs)

TOAST INFO: Creating noise model:  0.00 seconds (1 calls)

Create a high resolution point source map to convolve with the beam¶

import healpy as hp
import numpy as np
nside_high = 1024
npix_high = 12 * nside_high ** 2
pointsource_map = np.zeros([3, npix_high])
coords = []
for lon in np.linspace(0, 360, 9, endpoint=False):
    for lat in np.linspace(-90, 90, 7):
        pix = hp.ang2pix(nside_high, lon, lat, lonlat=True)
        # Add a completely unpolarized source and see if beam asymmetries manufacture polarization
        pointsource_map[0, pix] = 1
        coords.append((lon, lat))
coords = np.vstack(coords).T
hp.mollview(np.zeros(12), title="Input signal", cmap="coolwarm")
hp.projplot(np.pi/2 - np.radians(coords[1]), np.radians(coords[0]), 'o')
lmax_high = nside_high * 2
cl, alm = hp.anafast(pointsource_map, lmax=lmax_high, iter=0, alm=True)
hp.write_map("sim_sources_map.fits", hp.reorder(pointsource_map, r2n=True), nest=True, overwrite=True)
hp.write_alm(args.sky_file, alm, overwrite=True)

Create asymmetric beam¶

beam_map = np.zeros([3, npix_high])
x, y, z = hp.pix2vec(nside_high, np.arange(npix_high))

xvar = .01
yvar = 5 * xvar
beam = np.exp(-(x ** 2 / xvar + y ** 2 / yvar))
beam[z < 0] = 0
hp.mollview(beam, cmap="coolwarm", rot=[0, 90])
beam_map = np.zeros([3, npix_high])
beam_map[0] = beam
beam_map[1] = beam
bl, blm = hp.anafast(beam_map, lmax=lmax_high, iter=0, alm=True)
hp.write_alm(args.beam_file, blm, overwrite=True)

Now simulate sky signal¶

import toast

toast.todmap.OpPointingHpix(nside=args.nside, nest=True, mode="IQU").exec(data)

npix = 12 * args.nside ** 2
hitmap = np.zeros(npix)
tod = data.obs[0]["tod"]
for det in tod.local_dets:
    pixels = tod.cache.reference("pixels_{}".format(det))
    hitmap[pixels] = 1
hitmap[hitmap == 0] = hp.UNSEEN
hp.mollview(hitmap, nest=True, title="all hit pixels", cbar=False)
hp.graticule(22.5, verbose=False)

name = "signal"
toast.tod.OpCacheClear(name).exec(data)

conviqt = toast.todmap.OpSimConviqt(
    comm.comm_rank,
    args.sky_file,
    args.beam_file,
    lmax=512,  # Will use maximum from file
    beammmax=16,  # Will use maximum from file
    pol=True,
    fwhm=0,
    order=13,
    calibrate=True,
    dxx=True,
    out=name,
    quat_name=None,
    flag_name=None,
    flag_mask=255,
    common_flag_name=None,
    common_flag_mask=255,
    apply_flags=False,
    remove_monopole=False,
    remove_dipole=False,
    normalize_beam=True,
    verbosity=1,
)
conviqt.exec(data)

TOAST INFO: initialize sky for detector 0A:  0.22 seconds (1 calls)
TOAST INFO: initialize beam for detector 0A:  0.20 seconds (1 calls)
TOAST INFO: get detector pointing for 0A:  0.00 seconds (1 calls)
TOAST INFO: compute pointing angles for detector 0A:  0.01 seconds (1 calls)
TOAST INFO: pack input array for detector 0A:  0.00 seconds (1 calls)

0 : Convolving : pol = 1 lmax = 512 beammmax = 16 order = 13
Conviqt::convolve timing:
convolve                                     1.94 <     1.94 +-     0.00 <     1.94 (       1 <        1.0 +-        0.0 <        1)
    distribute_colatitudes                   0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    fillingBetaSeg                           0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todRedistribution5cm                     0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todgen                                   1.93 <     1.93 +-     0.00 <     1.93 (       1 <        1.0 +-        0.0 <        1)
        arrFillingcm                         0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
        interpolTOD                          1.92 <     1.92 +-     0.00 <     1.92 (       1 <        1.0 +-        0.0 <        1)
            itheta0SetUp                     0.03 <     0.03 +-     0.00 <     0.03 (       2 <        2.0 +-        0.0 <        2)
                ithetacalc                   0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
            conviqt_hemiscm_pol              1.75 <     1.75 +-     0.00 <     1.75 (       1 <        1.0 +-        0.0 <        1)
                wigner_init                  0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
                wigner_prepare               0.00 <     0.00 +-     0.00 <     0.00 (     816 <      816.0 +-        0.0 <      816)
                todAnnulus                   0.08 <     0.08 +-     0.00 <     0.08 (       2 <        2.0 +-        0.0 <        2)
                    alltoall_datacube        0.12 <     0.12 +-     0.00 <     0.12 (       2 <        2.0 +-        0.0 <        2)
            conviqt_tod_loop                 0.14 <     0.14 +-     0.00 <     0.14 (     312 <      312.0 +-        0.0 <      312)
    alltoall                                 0.00 <     0.00 +-     0.00 <     0.00 (       3 <        3.0 +-        0.0 <        3)
    sort TOD                                 0.01 <     0.01 +-     0.00 <     0.01 (       3 <        3.0 +-        0.0 <        3)

TOAST INFO: convolve detector 0A:  1.94 seconds (1 calls)
TOAST INFO: extract convolved data for 0A:  0.00 seconds (1 calls)
TOAST INFO: cache detector 0A:  0.00 seconds (1 calls)
TOAST INFO: conviqt process detector 0A:  2.38 seconds (1 calls)
TOAST INFO: initialize sky for detector 0B:  0.23 seconds (1 calls)
TOAST INFO: initialize beam for detector 0B:  0.22 seconds (1 calls)
TOAST INFO: get detector pointing for 0B:  0.00 seconds (1 calls)
TOAST INFO: compute pointing angles for detector 0B:  0.01 seconds (1 calls)
TOAST INFO: pack input array for detector 0B:  0.00 seconds (1 calls)

0 : Convolving : pol = 1 lmax = 512 beammmax = 16 order = 13
Conviqt::convolve timing:
convolve                                     1.93 <     1.93 +-     0.00 <     1.93 (       1 <        1.0 +-        0.0 <        1)
    distribute_colatitudes                   0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    fillingBetaSeg                           0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todRedistribution5cm                     0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todgen                                   1.92 <     1.92 +-     0.00 <     1.92 (       1 <        1.0 +-        0.0 <        1)
        arrFillingcm                         0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
        interpolTOD                          1.92 <     1.92 +-     0.00 <     1.92 (       1 <        1.0 +-        0.0 <        1)
            itheta0SetUp                     0.03 <     0.03 +-     0.00 <     0.03 (       2 <        2.0 +-        0.0 <        2)
                ithetacalc                   0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
            conviqt_hemiscm_pol              1.73 <     1.73 +-     0.00 <     1.73 (       1 <        1.0 +-        0.0 <        1)
                wigner_init                  0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
                wigner_prepare               0.00 <     0.00 +-     0.00 <     0.00 (     816 <      816.0 +-        0.0 <      816)
                todAnnulus                   0.06 <     0.06 +-     0.00 <     0.06 (       2 <        2.0 +-        0.0 <        2)
                    alltoall_datacube        0.12 <     0.12 +-     0.00 <     0.12 (       2 <        2.0 +-        0.0 <        2)
            conviqt_tod_loop                 0.15 <     0.15 +-     0.00 <     0.15 (     312 <      312.0 +-        0.0 <      312)
    alltoall                                 0.00 <     0.00 +-     0.00 <     0.00 (       3 <        3.0 +-        0.0 <        3)
    sort TOD                                 0.01 <     0.01 +-     0.00 <     0.01 (       3 <        3.0 +-        0.0 <        3)

TOAST INFO: convolve detector 0B:  1.94 seconds (1 calls)
TOAST INFO: extract convolved data for 0B:  0.00 seconds (1 calls)
TOAST INFO: cache detector 0B:  0.00 seconds (1 calls)
TOAST INFO: conviqt process detector 0B:  2.40 seconds (1 calls)
TOAST INFO: initialize sky for detector 1A:  0.23 seconds (1 calls)
TOAST INFO: initialize beam for detector 1A:  0.22 seconds (1 calls)
TOAST INFO: get detector pointing for 1A:  0.00 seconds (1 calls)
TOAST INFO: compute pointing angles for detector 1A:  0.01 seconds (1 calls)
TOAST INFO: pack input array for detector 1A:  0.00 seconds (1 calls)

0 : Convolving : pol = 1 lmax = 512 beammmax = 16 order = 13
Conviqt::convolve timing:
convolve                                     1.93 <     1.93 +-     0.00 <     1.93 (       1 <        1.0 +-        0.0 <        1)
    distribute_colatitudes                   0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    fillingBetaSeg                           0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todRedistribution5cm                     0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todgen                                   1.92 <     1.92 +-     0.00 <     1.92 (       1 <        1.0 +-        0.0 <        1)
        arrFillingcm                         0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
        interpolTOD                          1.91 <     1.91 +-     0.00 <     1.91 (       1 <        1.0 +-        0.0 <        1)
            itheta0SetUp                     0.03 <     0.03 +-     0.00 <     0.03 (       2 <        2.0 +-        0.0 <        2)
                ithetacalc                   0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
            conviqt_hemiscm_pol              1.77 <     1.77 +-     0.00 <     1.77 (       1 <        1.0 +-        0.0 <        1)
                wigner_init                  0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
                wigner_prepare               0.00 <     0.00 +-     0.00 <     0.00 (     816 <      816.0 +-        0.0 <      816)
                todAnnulus                   0.05 <     0.05 +-     0.00 <     0.05 (       2 <        2.0 +-        0.0 <        2)
                    alltoall_datacube        0.12 <     0.12 +-     0.00 <     0.12 (       2 <        2.0 +-        0.0 <        2)
            conviqt_tod_loop                 0.10 <     0.10 +-     0.00 <     0.10 (     314 <      314.0 +-        0.0 <      314)
    alltoall                                 0.00 <     0.00 +-     0.00 <     0.00 (       3 <        3.0 +-        0.0 <        3)
    sort TOD                                 0.01 <     0.01 +-     0.00 <     0.01 (       3 <        3.0 +-        0.0 <        3)

TOAST INFO: convolve detector 1A:  1.93 seconds (1 calls)
TOAST INFO: extract convolved data for 1A:  0.00 seconds (1 calls)
TOAST INFO: cache detector 1A:  0.00 seconds (1 calls)
TOAST INFO: conviqt process detector 1A:  2.40 seconds (1 calls)
TOAST INFO: initialize sky for detector 1B:  0.23 seconds (1 calls)
TOAST INFO: initialize beam for detector 1B:  0.23 seconds (1 calls)
TOAST INFO: get detector pointing for 1B:  0.00 seconds (1 calls)
TOAST INFO: compute pointing angles for detector 1B:  0.01 seconds (1 calls)
TOAST INFO: pack input array for detector 1B:  0.00 seconds (1 calls)

0 : Convolving : pol = 1 lmax = 512 beammmax = 16 order = 13
Conviqt::convolve timing:
convolve                                     2.05 <     2.05 +-     0.00 <     2.05 (       1 <        1.0 +-        0.0 <        1)
    distribute_colatitudes                   0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    fillingBetaSeg                           0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todRedistribution5cm                     0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todgen                                   2.04 <     2.04 +-     0.00 <     2.04 (       1 <        1.0 +-        0.0 <        1)
        arrFillingcm                         0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
        interpolTOD                          2.04 <     2.04 +-     0.00 <     2.04 (       1 <        1.0 +-        0.0 <        1)
            itheta0SetUp                     0.03 <     0.03 +-     0.00 <     0.03 (       2 <        2.0 +-        0.0 <        2)
                ithetacalc                   0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
            conviqt_hemiscm_pol              1.88 <     1.88 +-     0.00 <     1.88 (       1 <        1.0 +-        0.0 <        1)
                wigner_init                  0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
                wigner_prepare               0.00 <     0.00 +-     0.00 <     0.00 (     816 <      816.0 +-        0.0 <      816)
                todAnnulus                   0.10 <     0.10 +-     0.00 <     0.10 (       2 <        2.0 +-        0.0 <        2)
                    alltoall_datacube        0.14 <     0.14 +-     0.00 <     0.14 (       2 <        2.0 +-        0.0 <        2)
            conviqt_tod_loop                 0.12 <     0.12 +-     0.00 <     0.12 (     314 <      314.0 +-        0.0 <      314)
    alltoall                                 0.00 <     0.00 +-     0.00 <     0.00 (       3 <        3.0 +-        0.0 <        3)
    sort TOD                                 0.01 <     0.01 +-     0.00 <     0.01 (       3 <        3.0 +-        0.0 <        3)

TOAST INFO: convolve detector 1B:  2.05 seconds (1 calls)
TOAST INFO: extract convolved data for 1B:  0.00 seconds (1 calls)
TOAST INFO: cache detector 1B:  0.00 seconds (1 calls)
TOAST INFO: conviqt process detector 1B:  2.52 seconds (1 calls)
TOAST INFO: initialize sky for detector 2A:  0.23 seconds (1 calls)
TOAST INFO: initialize beam for detector 2A:  0.22 seconds (1 calls)
TOAST INFO: get detector pointing for 2A:  0.00 seconds (1 calls)
TOAST INFO: compute pointing angles for detector 2A:  0.01 seconds (1 calls)
TOAST INFO: pack input array for detector 2A:  0.00 seconds (1 calls)

0 : Convolving : pol = 1 lmax = 512 beammmax = 16 order = 13
Conviqt::convolve timing:
convolve                                     1.92 <     1.92 +-     0.00 <     1.92 (       1 <        1.0 +-        0.0 <        1)
    distribute_colatitudes                   0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    fillingBetaSeg                           0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todRedistribution5cm                     0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
    todgen                                   1.91 <     1.91 +-     0.00 <     1.91 (       1 <        1.0 +-        0.0 <        1)
        arrFillingcm                         0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
        interpolTOD                          1.90 <     1.90 +-     0.00 <     1.90 (       1 <        1.0 +-        0.0 <        1)
            itheta0SetUp                     0.03 <     0.03 +-     0.00 <     0.03 (       2 <        2.0 +-        0.0 <        2)
                ithetacalc                   0.00 <     0.00 +-     0.00 <     0.00 (       2 <        2.0 +-        0.0 <        2)
            conviqt_hemiscm_pol              1.75 <     1.75 +-     0.00 <     1.75 (       1 <        1.0 +-        0.0 <        1)
                wigner_init                  0.00 <     0.00 +-     0.00 <     0.00 (       1 <        1.0 +-        0.0 <        1)
                wigner_prepare               0.00 <     0.00 +-     0.00 <     0.00 (     816 <      816.0 +-        0.0 <      816)
                todAnnulus                   0.06 <     0.06 +-     0.00 <     0.06 (       2 <        2.0 +-        0.0 <        2)
                    alltoall_datacube        0.12 <     0.12 +-     0.00 <     0.12 (       2 <        2.0 +-        0.0 <        2)
            conviqt_tod_loop                 0.11 <     0.11 +-     0.00 <     0.11 (     308 <      308.0 +-        0.0 <      308)
    alltoall                                 0.00 <     0.00 +-     0.00 <     0.00 (       3 <        3.0 +-        0.0 <        3)
    sort TOD                                 0.01 <     0.01 +-     0.00 <     0.01 (       3 <        3.0 +-        0.0 <        3)

Destripe the signal and make a map. We use the nascent TOAST mapmaker because it can be run in serial mode without MPI. The TOAST mapmaker is still significantly slower so production runs should used libMadam.

mapmaker = toast.todmap.OpMapMaker(
    nside=args.nside,
    nnz=3,
    name=name,
    outdir=args.outdir,
    outprefix="toast_test_",
    baseline_length=10,
    # maskfile=self.maskfile_binary,
    # weightmapfile=self.maskfile_smooth,
    # subharmonic_order=None,
    iter_max=100,
    use_noise_prior=False,
    # precond_width=30,
)
mapmaker.exec(data)

TOAST INFO: Flag gaps:  0.00 seconds (1 calls)
TOAST INFO: Get detector weights:  0.00 seconds (1 calls)
TOAST INFO: Identify local submaps:  0.00 seconds (1 calls)
TOAST INFO: Accumulate N_pp'^1:  0.03 seconds (1 calls)
TOAST INFO: All reduce N_pp'^1:  0.00 seconds (1 calls)
TOAST INFO: Wrote hits to maps/toast_test_hits.fits
TOAST INFO: Write hits:  0.01 seconds (1 calls)
TOAST INFO: Wrote inverse white noise covariance to maps/toast_test_invnpp.fits
TOAST INFO: Write N_pp'^1:  0.01 seconds (1 calls)
TOAST INFO: Compute reciprocal condition numbers:  0.06 seconds (1 calls)
TOAST INFO: Wrote reciprocal condition numbers to maps/toast_test_rcond.fits
TOAST INFO: Write rcond:  0.01 seconds (1 calls)
TOAST INFO: Invert N_pp'^1:  0.04 seconds (1 calls)
TOAST INFO: Wrote white noise covariance to maps/toast_test_npp.fits
TOAST INFO: Write N_pp':  0.01 seconds (1 calls)
TOAST INFO:   Build noise-weighted map:  0.00 seconds (0 calls)
TOAST INFO:   Apply noise covariance:  0.00 seconds (0 calls)
TOAST INFO:   Write map to maps/toast_test_binned.fits:  0.00 seconds (0 calls)
TOAST INFO: Initializing offset template, step_length = 10
TOAST INFO: Initialize templates:  0.22 seconds (1 calls)
TOAST INFO: Initialize projection matrix:  0.00 seconds (1 calls)
TOAST INFO: Initialize projection matrix:  0.00 seconds (1 calls)
TOAST INFO: Initialize PCG solver:  0.04 seconds (1 calls)
TOAST INFO: Initial residual: 1.4928358557320724e-06
TOAST INFO: Iter =    0 relative residual:   1.3383e-02:  0.33 seconds (1 calls)
TOAST INFO: Iter =    1 relative residual:   5.7900e-04:  0.33 seconds (1 calls)
TOAST INFO: Iter =    2 relative residual:   1.3369e-04:  0.33 seconds (1 calls)
TOAST INFO: Iter =    3 relative residual:   2.8960e-05:  0.28 seconds (1 calls)
TOAST INFO: Iter =    4 relative residual:   8.8193e-06:  0.30 seconds (1 calls)
TOAST INFO: Iter =    5 relative residual:   3.1978e-06:  0.32 seconds (1 calls)
TOAST INFO: Iter =    6 relative residual:   1.1066e-06:  0.31 seconds (1 calls)
TOAST INFO: Iter =    7 relative residual:   3.9732e-07:  0.31 seconds (1 calls)
TOAST INFO: Iter =    8 relative residual:   2.7044e-07:  0.28 seconds (1 calls)
TOAST INFO: Iter =    9 relative residual:   1.5669e-07:  0.28 seconds (1 calls)
TOAST INFO: Iter =   10 relative residual:   7.0114e-08:  0.28 seconds (1 calls)
TOAST INFO: Iter =   11 relative residual:   5.9106e-08:  0.26 seconds (1 calls)
TOAST INFO: Iter =   12 relative residual:   8.2422e-08:  0.26 seconds (1 calls)
TOAST INFO: Iter =   13 relative residual:   4.1546e-08:  0.25 seconds (1 calls)
TOAST INFO: Iter =   14 relative residual:   2.3803e-08:  0.29 seconds (1 calls)
TOAST INFO: Iter =   15 relative residual:   1.9252e-08:  0.29 seconds (1 calls)
TOAST INFO: Iter =   16 relative residual:   1.2358e-08:  0.29 seconds (1 calls)
TOAST INFO: Iter =   17 relative residual:   8.3160e-09:  0.29 seconds (1 calls)
TOAST INFO: Iter =   18 relative residual:   5.6741e-09:  0.28 seconds (1 calls)
TOAST INFO: Iter =   19 relative residual:   4.6278e-09:  0.52 seconds (1 calls)
TOAST INFO: Iter =   20 relative residual:   6.5122e-09:  0.26 seconds (1 calls)
TOAST INFO: Iter =   21 relative residual:   1.0861e-08:  0.28 seconds (1 calls)
TOAST INFO: Iter =   22 relative residual:   1.6557e-08:  0.30 seconds (1 calls)
TOAST INFO: Iter =   23 relative residual:   2.8528e-08:  0.28 seconds (1 calls)
TOAST INFO: Iter =   24 relative residual:   5.8958e-08:  0.26 seconds (1 calls)
TOAST INFO: Iter =   25 relative residual:   5.0146e-08:  0.29 seconds (1 calls)
TOAST INFO: Iter =   26 relative residual:   4.1505e-08:  0.29 seconds (1 calls)
TOAST INFO: Iter =   27 relative residual:   9.8987e-08:  0.29 seconds (1 calls)
TOAST INFO: Iter =   28 relative residual:   3.2317e-07:  0.27 seconds (1 calls)
TOAST INFO: Iter =   29 relative residual:   9.2114e-07:  0.30 seconds (1 calls)
TOAST INFO: Iter =   30 relative residual:   1.3798e-06:  0.30 seconds (1 calls)
TOAST INFO: PCG stalled after 30 iterations:  9.20 seconds (1 calls)
TOAST INFO: 0 : Solution: template amplitudes:
"offset" : 
[ 1.71237528e-06  1.40446578e-06  5.40266906e-07 ...  1.30978221e-06
 -1.71061312e-06  1.11191004e-06]
TOAST INFO: Solve amplitudes:  9.20 seconds (1 calls)
TOAST INFO: Clean TOD:  0.02 seconds (1 calls)
TOAST INFO:   Build noise-weighted map:  0.00 seconds (0 calls)
TOAST INFO:   Apply noise covariance:  0.00 seconds (0 calls)
TOAST INFO:   Write map to maps/toast_test_destriped.fits:  0.00 seconds (0 calls)

Plot a segment of the timelines

plt.figure(figsize=[12, 8])

hitmap = hp.read_map("maps/toast_test_hits.fits")
hitmap[hitmap == 0] = hp.UNSEEN
hp.mollview(hitmap, sub=[2, 2, 1], title="hits")

binmap = hp.read_map("maps/toast_test_binned.fits")
binmap[binmap == 0] = hp.UNSEEN
hp.mollview(binmap, sub=[2, 2, 2], title="binned map", cmap="coolwarm")

destriped = hp.read_map("maps/toast_test_destriped.fits")
destriped[destriped == 0] = hp.UNSEEN
hp.mollview(destriped, sub=[2, 2, 3], title="destriped map", cmap="coolwarm")

inmap = hp.ud_grade(hp.read_map("sim_sources_map.fits"), args.nside)
inmap[hitmap == hp.UNSEEN] = hp.UNSEEN
hp.mollview(inmap, sub=[2, 2, 4], title="input map", cmap="coolwarm")

NSIDE = 64
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING
NSIDE = 64
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING
NSIDE = 64
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING
NSIDE = 1024
ORDERING = NESTED in fits file
INDXSCHM = IMPLICIT
Ordering converted to RING

Exercises¶

Plot the polarization of the simulated signal above
Modify the scan strategy so that the beam elongation is more visible