<jupyter_start><jupyter_text>LiteBIRD Simulation ExampleThis notebook does a simple LiteBIRD simulation which you can use as a starting point for testing and customization.  The toast_litebird package is here: https://github.com/hpc4cmb/toast-litebird and the documentation is here: https://hpc4cmb.github.io/toast-litebird/First you must get access to the kernel that has the toast_litebird package.  Open a jupyter terminal and do:```%>  module use /global/common/software/litebird/cori/modulefiles    %>  module load litebird    %>  litebird-jupyter.sh``` Now in this notebook select the the litebird kernel.  You may have to shutdown this notebook and re-open to see the new kernel.<jupyter_code>import os
import sys

import healpy as hp
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np

import toast
from toast.mpi import MPI

from toast.utils import memreport

from toast import pipeline_tools

from toast_litebird import pipeline_tools as lbtools

# Capture C++ output in the jupyter cells
%reload_ext wurlitzer<jupyter_output><empty_output><jupyter_text>Select DetectorsWe can use some command line tools to easily select detectors and dump them to a file for use in a pipeline.  This command creates a full hardware model:<jupyter_code>! lb_hardware_sim --overwrite<jupyter_output><empty_output><jupyter_text>We can look at the details of the hardware file:<jupyter_code>! lb_hardware_info hardware.toml.gz<jupyter_output><empty_output><jupyter_text>Now we can select just some detectors<jupyter_code>! lb_hardware_trim \
--hardware hardware.toml.gz \
--overwrite \
--out selected \
--telescopes LFT \
--match 'band:.*040'
! lb_hardware_info selected.toml.gz<jupyter_output><empty_output><jupyter_text>Plot this<jupyter_code>! lb_hardware_plot --hardware selected.toml.gz --out selected.pdf<jupyter_output><empty_output><jupyter_text>The previous command makes a PDF file.  We can display it:<jupyter_code>from IPython.display import IFrame
IFrame("selected.pdf", width=600, height=300)<jupyter_output><empty_output><jupyter_text>ParametersThese arguments control the entire notebook<jupyter_code>class args:
    # Hardware model
    hardware = "selected.toml.gz"
    bands = "040"
    thinfp = False
    # Observations
    obs_num = 1
    start_time = 0
    sample_rate = 30.0
    obs_time_h = 23.0
    gap_h = 1.0
    # half-wave plate
    hwp_rpm = 91.0
    hwp_step_deg = None
    hwp_step_time_s = None
    # Scanning parameters
    spin_period_min = 10.0
    spin_angle_deg = 50.0      # This is "beta"
    prec_period_min = 96.174
    prec_angle_deg = 45.0      # This is "alpha"
    # Pixelization
    coord = "G"
    nside = 512
    mode = "IQU"
    single_precision_pointing = False
    nside_submap = 16
    # Noise simulation
    common_mode_noise = False
    # Output directory
    outdir = "litebird_out"<jupyter_output><empty_output><jupyter_text>CommunicatorSince this is a serial notebook, this communicator will just have one process.<jupyter_code>mpiworld, procs, rank = toast.mpi.get_world()
comm = toast.mpi.Comm(mpiworld)

memreport("After communicator creation", comm.comm_world)<jupyter_output><empty_output><jupyter_text>FocalplaneLoad the hardware file and create the focalplane.<jupyter_code>hw, telescope = lbtools.get_hardware(args, comm)

focalplane = lbtools.get_focalplane(args, comm, hw)

memreport("After focalplane creation", comm.comm_world)<jupyter_output><empty_output><jupyter_text>Create ObservationsThis uses the parameters at the top of the notebook to simulate regular spaced observations.<jupyter_code>data = lbtools.create_observations(args, comm, focalplane, 1)

memreport("After creating observations", comm.comm_world)<jupyter_output><empty_output><jupyter_text>Pointing matrixHere we translate the boresight quaternions into detector pointing (pixels numbers and Stokes weights).<jupyter_code>pipeline_tools.expand_pointing(args, comm, data)

memreport("After expanding pointing", comm.comm_world)<jupyter_output><empty_output><jupyter_text>Make a boolean hit map for diagnostics<jupyter_code>npix = 12 * args.nside ** 2
hitmap = np.zeros(npix)
for obs in data.obs:
    tod = obs["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)<jupyter_output><empty_output><jupyter_text>Sky signalCreate a synthetic Gaussian map to scan as input signal<jupyter_code>lmax = args.nside * 2
cls = np.zeros([4, lmax + 1])
cls[0] = 1e0
sim_map = hp.synfast(cls, args.nside, lmax=lmax, fwhm=np.radians(15), new=True)

plt.figure(figsize=[12, 8])
for i, m in enumerate(sim_map):
    hp.mollview(sim_map[i], cmap="coolwarm", title="Input signal {}".format("IQU"[i]), sub=[1, 3, 1+i])
hp.write_map("sim_map.fits", hp.reorder(sim_map, r2n=True), nest=True, overwrite=True)<jupyter_output><empty_output><jupyter_text>Scan the sky signal<jupyter_code>all_name = "all_signal"
sky_name = "sky_signal"

# Clear any previous signal from the buffers
toast.tod.OpCacheClear(all_name).exec(data)

distmap = toast.map.DistPixels(
    data,
    nnz=len(args.mode),
    dtype=np.float32,
)
distmap.read_healpix_fits("sim_map.fits")

toast.todmap.OpSimScan(distmap=distmap, out=all_name).exec(data)

# Copy the sky signal, just in case we need it later

toast.tod.OpCacheCopy(input=all_name, output=sky_name, force=True).exec(data)

memreport("After scanning sky signal", comm.comm_world)<jupyter_output><empty_output><jupyter_text>NoiseSimulate noise and make a copy of signal+noise in case we need it later<jupyter_code>copy_name = "copy_signal"

toast.tod.OpSimNoise(out=all_name, realization=0).exec(data)

toast.tod.OpCacheCopy(input=all_name, output=copy_name, force=True).exec(data)

memreport("After simulating noise", comm.comm_world)<jupyter_output><empty_output><jupyter_text>Your own operator hereHere we define an empty operator you can work with<jupyter_code>class MyOperator(toast.Operator):
    def __init__(self, name="signal"):
        """ Arguments:
        name(str) : Cache prefix to operate on
        """
        self._name = name
    
    def exec(self, data):
        # We loop here over all local data but do nothing with it.
        for obs in data.obs:
            tod = obs["tod"]
            for det in tod.local_dets:
                signal = tod.local_signal(det, self._name)
                # Do operations in-place
                signal *= 1.0
                #signal[:] = (some other data)<jupyter_output><empty_output><jupyter_text>Then we apply the operator to the data<jupyter_code>toast.tod.OpCacheCopy(input=copy_name, output=all_name, force=True).exec(data)
MyOperator(name=all_name).exec(data)

memreport("After my operator", comm.comm_world)<jupyter_output><empty_output><jupyter_text>Plot a short segment of the signal before and after the operator<jupyter_code>tod = data.obs[0]["tod"]
times = tod.local_times()

fig = plt.figure(figsize=[12, 8])
for idet, det in enumerate(tod.local_dets):
    cflags = tod.local_common_flags()
    before = tod.local_signal(det, copy_name)
    after = tod.local_signal(det, all_name)

    ind = slice(0, 1000)
    # Flag out turnarounds
    good = cflags[ind] == 0
    ax = fig.add_subplot(8, 8, 1 + idet)
    ax.set_title(det)
    ax.plot(times[ind][good], before[ind][good], '.', label="before")
    ax.plot(times[ind][good], after[ind][good], '.', label="after")
ax.legend(bbox_to_anchor=(1.1, 1.00))
fig.subplots_adjust(hspace=0.6)<jupyter_output><empty_output><jupyter_text>Make a mapDestripe 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`.<jupyter_code># Always begin mapmaking by copying the simulated signal.

destriped_name = "destriped"
toast.tod.OpCacheCopy(
    input=all_name,
    output=destriped_name,
    force=True
).exec(data)

mapmaker = toast.todmap.OpMapMaker(
    nside=args.nside,
    nnz=3,
    name=destriped_name,
    outdir=args.outdir,
    outprefix="toast_test_",
    baseline_length=10,
    iter_max=15,
    use_noise_prior=False,
)
mapmaker.exec(data)

memreport("After map making", comm.comm_world)<jupyter_output><empty_output><jupyter_text>Plot a segment of the timelines<jupyter_code>tod = data.obs[0]["tod"]
times = tod.local_times()

fig = plt.figure(figsize=[12, 8])
for idet, det in enumerate(tod.local_dets):
    sky = tod.local_signal(det, sky_name)
    full = tod.local_signal(det, all_name)
    destriped = tod.local_signal(det, destriped_name)

    ind = slice(0, 1000)
    ax = fig.add_subplot(8, 8, 1 + idet)
    ax.set_title(det)
    ax.plot(times[ind], sky[ind], '.', label="sky", zorder=100)
    ax.plot(times[ind], full[ind] - sky[ind], '.', label="noise")
    ax.plot(times[ind], full[ind] - destriped[ind], '.', label="baselines")
ax.legend(bbox_to_anchor=(1.1, 1.00))
fig.subplots_adjust(hspace=0.6)

fig = plt.figure(figsize=[12, 8])
for idet, det in enumerate(tod.local_dets):
    sky = tod.local_signal(det, sky_name)
    full = tod.local_signal(det, copy_name)
    destriped = tod.local_signal(det, destriped_name)
    ax = fig.add_subplot(8, 8, 1 + idet)
    ax.set_title(det)
    #plt.plot(times[ind], sky[ind], '-', label="signal", zorder=100)
    plt.plot(times, full - sky, '.', label="noise")
    plt.plot(times, full - destriped, '.', label="baselines")
ax.legend(bbox_to_anchor=(1.1, 1.00))
fig.subplots_adjust(hspace=.6)
plt.figure(figsize=[16, 8])

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

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

# Fix the plotting range for input signal and the destriped map
amp = 5.0

destriped = hp.read_map("litebird_out/toast_test_destriped.fits", verbose=False)
destriped[destriped == 0] = hp.UNSEEN
# Remove monopole
good = destriped != hp.UNSEEN
destriped[good] -= np.median(destriped[good])
hp.mollview(destriped, sub=[2, 2, 3], title="destriped map", cmap="coolwarm", min=-amp, max=amp)

inmap = hp.read_map("sim_map.fits", verbose=False)
inmap[hitmap == hp.UNSEEN] = hp.UNSEEN
hp.mollview(inmap, sub=[2, 2, 4], title="input map", cmap="coolwarm", min=-amp, max=amp)

# Plot the white noise covariance

plt.figure(figsize=[12, 8])
wcov = hp.read_map("litebird_out/toast_test_npp.fits", None)
wcov[:, wcov[0] == 0] = hp.UNSEEN
hp.mollview(wcov[0], sub=[3, 3, 1], title="II", cmap="coolwarm")
hp.mollview(wcov[1], sub=[3, 3, 2], title="IQ", cmap="coolwarm")
hp.mollview(wcov[2], sub=[3, 3, 3], title="IU", cmap="coolwarm")
hp.mollview(wcov[3], sub=[3, 3, 5], title="QQ", cmap="coolwarm")
hp.mollview(wcov[4], sub=[3, 3, 6], title="QU", cmap="coolwarm")
hp.mollview(wcov[5], sub=[3, 3, 9], title="UU", cmap="coolwarm")<jupyter_output><empty_output><jupyter_text>Filter & binA filter-and-bin mapmaker is easily created by combining TOAST filter operators and running the mapmaker without destriping:<jupyter_code>filtered_name = "filtered"

toast.tod.OpCacheCopy(input=copy_name, output=filtered_name, force=True).exec(data)

toast.tod.OpPolyFilter(order=3, name=filtered_name).exec(data)

mapmaker = toast.todmap.OpMapMaker(
    nside=args.nside,
    nnz=len(args.mode),
    name=filtered_name,
    outdir=args.outdir,
    outprefix="toast_test_filtered_",
    baseline_length=None,
)
mapmaker.exec(data)

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

binmap = hp.read_map("litebird_out/toast_test_binned.fits", verbose=False)
binmap[binmap == 0] = hp.UNSEEN
hp.mollview(binmap, sub=[1, 3, 1], title="binned map", cmap="coolwarm")

filtered_map = hp.read_map("litebird_out/toast_test_filtered_binned.fits", verbose=False)
filtered_map[filtered_map == 0] = hp.UNSEEN
hp.mollview(filtered_map, sub=[1, 3, 2], title="filtered map", cmap="coolwarm")

inmap = hp.read_map("sim_map.fits", verbose=False)
inmap[binmap == hp.UNSEEN] = hp.UNSEEN
hp.mollview(inmap, sub=[1, 3, 3], title="input map", cmap="coolwarm")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12827417380660955, 'fraction_numerical': 0.016401468788249694, 'mean_word_length': 5.183652875882946, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 88, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '35968', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 4711, 'n_tokens_neox': 4443, 'n_words': 986, 'path': 'hackathon/litebird_sim.ipynb', 'repo_name': 'ickc/toast-workshop', 'size': '12255'}