rough implementation of the mapmaker pipeline tools. Includes mappin…

src/toast/pipeline_tools/CMakeLists.txt

@@ -12,6 +12,7 @@ install(FILES
   export.py
   filters.py
   gain.py
+  mapmaker.py
   madam.py
   noise.py
   pointing.py

src/toast/pipeline_tools/__init__.py

@@ -26,6 +26,7 @@
     apply_groundfilter,
 )
 from .gain import add_gainscrambler_args, scramble_gains
+from .mapmaker import add_mapmaker_args, apply_mapmaker
 from .madam import add_madam_args, setup_madam, apply_madam
 from .noise import add_noise_args, simulate_noise, get_analytic_noise
 from .pointing import add_pointing_args, expand_pointing

src/toast/pipeline_tools/madam.py

@@ -20,21 +20,27 @@ def add_madam_args(parser):
     """
 
     parser.add_argument(
-        "--madam-prefix", required=False, default="toast", help="Output map prefix"
+        "--madam-prefix",
+        required=False,
+        default="toast",
+        help="Output map prefix",
+        dest="mapmaker_prefix",
     )
     parser.add_argument(
         "--madam-iter-max",
         required=False,
         default=1000,
         type=np.int,
         help="Maximum number of CG iterations in Madam",
+        dest="mapmaker_iter_max",
     )
     parser.add_argument(
         "--madam-precond-width",
         required=False,
         default=100,
         type=np.int,
         help="Width of the Madam band preconditioner",
+        dest="mapmaker_precond_width",
     )
     parser.add_argument(
         "--madam-precond-width-min",
@@ -54,6 +60,7 @@ def add_madam_args(parser):
         default=10000.0,
         type=np.float,
         help="Destriping baseline length (seconds)",
+        dest="mapmaker_baseline_length",
     )
     parser.add_argument(
         "--madam-baseline-order",
@@ -68,6 +75,7 @@ def add_madam_args(parser):
         default=False,
         action="store_true",
         help="Destripe with the noise filter enabled",
+        dest="mapmaker_noisefilter",
     )
     parser.add_argument(
         "--madam-parfile", required=False, default=None, help="Madam parameter file"
@@ -308,25 +316,25 @@ def setup_madam(args):
                         key, value = result.group(1), result.group(2)
                         pars[key] = value
 
-    pars["base_first"] = args.madam_baseline_length
+    pars["base_first"] = args.mapmaker_baseline_length
     pars["basis_order"] = args.madam_baseline_order
     # Adaptive preconditioner width
     width_min = args.madam_precond_width_min
     width_max = args.madam_precond_width_max
     if width_min is None:
         # madam-precond-width has a default value
-        width_min = args.madam_precond_width
+        width_min = args.mapmaker_precond_width
     if width_max is None:
         # madam-precond-width has a default value
-        width_max = args.madam_precond_width
+        width_max = args.mapmaker_precond_width
     if width_min > width_max:
         # it is not an error for these two to match
         width_min = width_max
     pars["precond_width_min"] = width_min
     pars["precond_width_max"] = width_max
     #
     pars["nside_map"] = args.nside
-    if args.madam_noisefilter:
+    if args.mapmaker_noisefilter:
         if args.madam_baseline_order != 0:
             raise RuntimeError(
                 "Madam cannot build a noise filter when baseline"
@@ -336,8 +344,8 @@ def setup_madam(args):
     else:
         pars["kfilter"] = False
     pars["fsample"] = args.sample_rate
-    pars["iter_max"] = args.madam_iter_max
-    pars["file_root"] = args.madam_prefix
+    pars["iter_max"] = args.mapmaker_iter_max
+    pars["file_root"] = args.mapmaker_prefix
 
     # Translate boolean values.  Madam knows how to do this but it
     # simplifies pipeline_tools/madam.py
@@ -431,7 +439,7 @@ def apply_madam(
             log.info("No Madam outputs requested.  Skipping.")
         return
 
-    if args.madam_noisefilter or not pars["kfirst"]:
+    if args.mapmaker_noisefilter or not pars["kfirst"]:
         # With the noise filter enabled, we want to enforce continuity
         # across the Observation.  Otherwise we fit each interval
         # separately.

src/toast/pipeline_tools/mapmaker.py

@@ -0,0 +1,230 @@
+# Copyright (c) 2019-2020 by the parties listed in the AUTHORS file.
+# All rights reserved.  Use of this source code is governed by
+# a BSD-style license that can be found in the LICENSE file.
+
+import argparse
+import copy
+import os
+import re
+
+import numpy as np
+
+from ..timing import function_timer, Timer
+from ..utils import Logger, Environment
+
+from ..todmap import OpMapMaker
+
+
+def add_mapmaker_args(parser):
+    """ Add mapmaker arguments
+    """
+    parser.add_argument(
+        "--mapmaker-prefix",
+        required=False,
+        default="toast",
+        help="Output map prefix",
+        dest="mapmaker_prefix",
+    )
+    parser.add_argument(
+        "--mapmaker-iter-max",
+        required=False,
+        default=1000,
+        type=np.int,
+        help="Maximum number of CG iterations",
+        dest="mapmaker_iter_max",
+    )
+    parser.add_argument(
+        "--mapmaker-precond-width",
+        required=False,
+        default=100,
+        type=np.int,
+        help="Width of the Madam band preconditioner",
+        dest="mapmaker_precond_width",
+    )
+    parser.add_argument(
+        "--mapmaker-baseline-length",
+        required=False,
+        default=10000.0,
+        type=np.float,
+        help="Destriping baseline length (seconds)",
+        dest="mapmaker_baseline_length",
+    )
+    parser.add_argument(
+        "--mapmaker-noisefilter",
+        required=False,
+        default=False,
+        action="store_true",
+        help="Destripe with the noise filter enabled",
+        dest="mapmaker_noisefilter",
+    )
+    try:
+        parser.add_argument(
+            "--binmap",
+            required=False,
+            action="store_true",
+            help="Write binned maps [default]",
+            dest="write_binmap",
+        )
+        parser.add_argument(
+            "--no-binmap",
+            required=False,
+            action="store_false",
+            help="Do not write binned maps",
+            dest="write_binmap",
+        )
+        parser.set_defaults(write_binmap=True)
+    except argparse.ArgumentError:
+        pass
+
+    try:
+        parser.add_argument(
+            "--hits",
+            required=False,
+            action="store_true",
+            help="Write hit maps [default]",
+            dest="write_hits",
+        )
+        parser.add_argument(
+            "--no-hits",
+            required=False,
+            action="store_false",
+            help="Do not write hit maps",
+            dest="write_hits",
+        )
+        parser.set_defaults(write_hits=True)
+    except argparse.ArgumentError:
+        pass
+
+    try:
+        parser.add_argument(
+            "--wcov",
+            required=False,
+            action="store_true",
+            help="Write white noise covariance [default]",
+            dest="write_wcov",
+        )
+        parser.add_argument(
+            "--no-wcov",
+            required=False,
+            action="store_false",
+            help="Do not write white noise covariance",
+            dest="write_wcov",
+        )
+        parser.set_defaults(write_wcov=True)
+    except argparse.ArgumentError:
+        pass
+
+    try:
+        parser.add_argument(
+            "--wcov-inv",
+            required=False,
+            action="store_true",
+            help="Write inverse white noise covariance [default]",
+            dest="write_wcov_inv",
+        )
+        parser.add_argument(
+            "--no-wcov-inv",
+            required=False,
+            action="store_false",
+            help="Do not write inverse white noise covariance",
+            dest="write_wcov_inv",
+        )
+        parser.set_defaults(write_wcov_inv=True)
+    except argparse.ArgumentError:
+        pass
+
+    # `nside` may already be added
+    try:
+        parser.add_argument(
+            "--nside", required=False, default=512, type=np.int, help="Healpix NSIDE"
+        )
+    except argparse.ArgumentError:
+        pass
+    # Common flag mask may already be added
+    try:
+        parser.add_argument(
+            "--common-flag-mask",
+            required=False,
+            default=1,
+            type=np.uint8,
+            help="Common flag mask",
+        )
+    except argparse.ArgumentError:
+        pass
+
+    try:
+        parser.add_argument(
+            "--zip",
+            required=False,
+            action="store_true",
+            help="Compress the map outputs",
+            dest="zip_maps",
+        )
+        parser.add_argument(
+            "--no-zip",
+            required=False,
+            action="store_false",
+            help="Do not compress the map outputs",
+            dest="zip_maps",
+        )
+        parser.set_defaults(zip_maps=True)
+    except argparse.ArgumentError:
+        pass
+
+    return
+
+
+@function_timer
+def apply_mapmaker(
+    args,
+    comm,
+    data,
+    outpath,
+    cache_name,
+    time_comms=None,
+    telescope_data=None,
+    first_call=True,
+    extra_prefix=None,
+    verbose=True,
+    bin_only=False,
+):
+    log = Logger.get()
+    timer = Timer()
+
+    if outpath is None:
+        outpath = args.out
+
+    mapmaker = OpMapMaker(
+        nside=args.nside,
+        nnz=3,
+        name=cache_name,
+        outdir=outpath,
+        outprefix=args.mapmaker_prefix,
+        write_hits=args.write_hits,
+        zip_maps=args.zip_maps,
+        write_wcov_inv=args.write_wcov_inv,
+        write_wcov=args.write_wcov,
+        write_binned=args.write_binmap,
+        write_destriped=True,
+        write_rcond=True,
+        rcond_limit=1e-3,
+        baseline_length=args.mapmaker_baseline_length,
+        maskfile=None,
+        weightmapfile=None,
+        common_flag_mask=args.common_flag_mask,
+        flag_mask=1,
+        intervals="intervals",
+        subharmonic_order=None,
+        iter_min=3,
+        iter_max=args.mapmaker_iter_max,
+        use_noise_prior=args.mapmaker_noisefilter,
+        precond_width=args.mapmaker_precond_width,
+        pixels="pixels",
+    )
+
+    mapmaker.exec(data)
+
+    if comm.world_rank == 0 and verbose:
+        timer.report_clear("  OpMapMaker")
+
+    return

src/toast/pybind11/pybind11/__init__.py

@@ -7,17 +7,17 @@ def get_include(user=False):
     import sys
 
     # Are we running in a virtual environment?
-    virtualenv = hasattr(sys, 'real_prefix') or \
-        sys.prefix != getattr(sys, "base_prefix", sys.prefix)
+    virtualenv = hasattr(sys, "real_prefix") or sys.prefix != getattr(
+        sys, "base_prefix", sys.prefix
+    )
 
     if virtualenv:
-        return os.path.join(sys.prefix, 'include', 'site',
-                            'python' + sys.version[:3])
+        return os.path.join(sys.prefix, "include", "site", "python" + sys.version[:3])
     else:
-        dist = Distribution({'name': 'pybind11'})
+        dist = Distribution({"name": "pybind11"})
         dist.parse_config_files()
 
-        dist_cobj = dist.get_command_obj('install', create=True)
+        dist_cobj = dist.get_command_obj("install", create=True)
 
         # Search for packages in user's home directory?
         if user: