Merge pull request #16 from CMB-S4/85-tube_LAT

Draft pre-baseline design (85-tube LATs, SAT tube config)

.gitignore

@@ -0,0 +1,63 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+env/
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+*.egg-info/
+.installed.cfg
+*.egg
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*,cover
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# atmospheric cache
+atm_cache/
+
+# output directories
+out*/

pipelines/pickle_hardware.py

@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+
+# Copyright (c) 2020-2020 CMB-S4 Collaboration..
+# Full license can be found in the top level "LICENSE" file.
+
+# This script will convert TOML hardware maps into Python pickle
+# format for efficient access
+
+import os
+import pickle
+import sys
+
+from s4sim import hardware
+
+
+if len(sys.argv) < 2:
+    print("Usage: pickle_hardware.py <TOML file1> [<TOML file2>] ...")
+else:
+    for fname_in in sys.argv[1:]:
+        print(f"Loading {fname_in}")
+        hw = hardware.Hardware(fname_in)
+
+        fname_out = fname_in.replace(".toml", "").replace(".gz", "") + ".pkl"
+        print(f"Writing {fname_out}")
+        with open(fname_out, "wb") as fout:
+            pickle.dump(hw, fout)

pipelines/s4_hardware_to_wafers.py

@@ -0,0 +1,80 @@
+#!/usr/bin/env/python
+
+# Copyright (c) 2020-2020 CMB-S4 Collaboration.
+# Full license can be found in the top level "LICENSE" file.
+"""Trim a hardware model to include only some detectors.
+"""
+
+from collections import OrderedDict
+import sys
+import argparse
+
+import numpy as np
+import toast.qarray as qa
+
+from s4sim.hardware import Hardware
+
+
+XAXIS, YAXIS, ZAXIS = np.eye(3)
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="This program reads a hardware model from disk "
+        "and writes out per-wafer geometry.",
+        usage="s4_hardware_trim [options] (use --help for details)",
+    )
+
+    parser.add_argument(
+        "--hardware", required=True, default=None, help="Input hardware file"
+    )
+
+    parser.add_argument(
+        "--out",
+        required=False,
+        default="trimmed",
+        help="Name (without extensions) of the output hardware file",
+    )
+
+    args = parser.parse_args()
+
+    print("Loading hardware from {}...".format(args.hardware), flush=True)
+    hw = Hardware(args.hardware)
+
+    wafers = {}
+
+    for det_name, det_data in hw.data["detectors"].items():
+        wafer = det_data["wafer"]
+        pol = det_data["pol"]
+        quat = qa.norm(det_data["quat"])
+        if wafer not in wafers:
+            wafers[wafer] = {}
+        wx, wy, wz = qa.rotate(quat, ZAXIS)
+        wdir = np.array([wx, wy, wz])
+        posrot = qa.norm(qa.from_vectors(ZAXIS, wdir))
+        angrot = qa.norm(qa.mult(qa.inv(posrot), quat))
+        psi = np.degrees(qa.to_angles(angrot)[2]) % 180
+        wafers[wafer][det_name] = np.array([wx, wy, psi])
+
+    for wafer_name in sorted(wafers):
+        wafer_data = wafers[wafer_name]
+        ndet = len(wafer_data)
+        det_data = np.empty([ndet, 3])
+        for idet, det_name in enumerate(wafer_data):
+            det_data[idet] = wafer_data[det_name]
+        xoffset, yoffset, psioffset = np.mean(det_data, 0)
+        fname = f"wafer_{wafer_name}.txt"
+        with open(fname, "w") as fout:
+            for det_name in sorted(wafer_data):
+                x, y, psi = wafer_data[det_name]
+                x -= xoffset
+                y -= yoffset
+                phi = np.degrees(np.arcsin(x))
+                theta = np.degrees(np.arcsin(y))
+                fout.write(f"{det_name} {phi} {theta} {psi}\n")
+
+    return
+
+
+if __name__ == "__main__":
+    main()

pipelines/toast_s4_sim.py

@@ -30,6 +30,7 @@
 import re
 import sys
 import traceback
+from time import time
 
 import argparse
 import dateutil.parser
@@ -76,6 +77,7 @@ def parse_arguments(comm):
     toast_tools.add_todground_args(parser)
     toast_tools.add_pointing_args(parser)
     toast_tools.add_polyfilter_args(parser)
+    toast_tools.add_polyfilter2D_args(parser)
     toast_tools.add_groundfilter_args(parser)
     toast_tools.add_atmosphere_args(parser)
     toast_tools.add_noise_args(parser)
@@ -107,6 +109,14 @@ def parse_arguments(comm):
         help="Skip the first Madam call.",
     )
 
+    parser.add_argument(
+        "--pairdiff",
+        required=False,
+        default=False,
+        action="store_true",
+        help="Pair-difference TOD and pointing.",
+    )
+
     parser.add_argument(
         "--outdir", required=False, default="out", help="Output directory"
     )
@@ -138,11 +148,7 @@ def parse_arguments(comm):
         comm = Comm(groupsize=args.group_size)
 
     if comm.world_rank == 0:
-        if not os.path.isdir(args.outdir):
-            try:
-                os.makedirs(args.outdir)
-            except FileExistsError:
-                pass
+        os.makedirs(args.outdir, exist_ok=True)
         timer.report_clear("Parse arguments")
 
     return args, comm
@@ -151,11 +157,7 @@ def parse_arguments(comm):
 def setup_output(args, comm, mc):
     outpath = "{}/{:08}".format(args.outdir, mc)
     if comm.world_rank == 0:
-        if not os.path.isdir(outpath):
-            try:
-                os.makedirs(outpath)
-            except FileExistsError:
-                pass
+        os.makedirs(outpath, exist_ok=True)
     return outpath
 
 
@@ -170,27 +172,84 @@ def outputs_exist(args, comm, outpath):
                     outpath, args.mapmaker_prefix + "_telescope_all_time_all_bmap.fits"
                 )
                 there = os.path.isfile(fname)
+                if there:
+                    print(f"{fname} exists", flush=True)
+                else:
+                    print(f"{fname} does not exist", flush=True)
             if there and args.destripe:
                 fname = os.path.join(
                     outpath, args.mapmaker_prefix + "_telescope_all_time_all_map.fits"
                 )
                 there = os.path.isfile(fname)
+                if there:
+                    print(f"{fname} exists", flush=True)
+                else:
+                    print(f"{fname} does not exist", flush=True)
         if there and (args.apply_polyfilter or args.apply_groundfilter):
             fname = os.path.join(
                 outpath,
                 args.mapmaker_prefix + "_filtered" + "_telescope_all_time_all_bmap.fits",
             )
             there = os.path.isfile(fname)
+            if there:
+                print(f"{fname} exists", flush=True)
+            else:
+                print(f"{fname} does not exist", flush=True)
         if there and (args.filterbin_ground_order or args.filterbin_poly_order):
             fname = os.path.join(
                 outpath,
                 args.filterbin_prefix + "_telescope_all_time_all_filtered.fits",
             )
-            there = os.path.isfile(fname)
+            there = os.path.isfile(fname) or os.path.isfile(fname + ".gz")
+            if there:
+                print(f"{fname} exists", flush=True)
+            else:
+                print(f"{fname} does not exist", flush=True)
     there = comm.comm_world.bcast(there)
     return there
 
 
+def pairdiff(data, args, comm, name, do_pointing):
+    if not args.pairdiff:
+        return
+    t1 = time()
+    if comm.comm_world.rank == 0:
+        print("Pair differencing data", flush=True)
+
+    for obs in data.obs:
+        tod = obs["tod"]
+        for det in tod.local_dets:
+            signal = tod.local_signal(det, name)
+            if det.endswith("A"):
+                pairdet = det[:-1] + "B"
+                if pairdet not in tod.local_dets:
+                    raise RuntimeError(
+                        "Detector pair not available ({}, {})".format(det, pairdet)
+                    )
+            else:
+                continue
+            # signal
+            pairsignal = tod.local_signal(pairdet, name)
+            signal[:], pairsignal[:] = [
+                0.5 * (signal + pairsignal), 0.5 * (signal - pairsignal)
+            ]
+            if do_pointing:
+                # flags
+                flags = tod.local_flags(det)
+                pairflags = tod.local_flags(pairdet)
+                flags |= pairflags
+                pairflags[:] = flags
+                # pointing weights
+                weights = tod.cache.reference("weights_" + det)
+                pairweights = tod.cache.reference("weights_" + pairdet)
+                weights[:], pairweights[:] = [
+                    0.5 * (weights + pairweights), 0.5 * (weights - pairweights)
+                ]
+
+    if comm.comm_world.rank == 0:
+        print("Pair differenced in {:.1f} s".format(time() - t1), flush=True)
+    return
+
 def main():
     log = Logger.get()
     gt = GlobalTimers.get()
@@ -327,6 +386,8 @@ def main():
 
         # Bin and destripe maps
 
+        pairdiff(data, args, comm, totalname, mc == firstmc)
+
         if not args.skip_madam:
             toast_tools.apply_madam(
                 args,
@@ -357,10 +418,12 @@ def main():
                 first_call=(mc == firstmc),
             )
 
-        if args.apply_polyfilter or args.apply_groundfilter:
+        if args.apply_polyfilter or args.apply_groundfilter or args.apply_polyfilter2D:
 
             # Filter signal
 
+            toast_tools.apply_polyfilter2D(args, comm, data, totalname)
+
             toast_tools.apply_polyfilter(args, comm, data, totalname)
 
             toast_tools.apply_groundfilter(args, comm, data, totalname)