Merge pull request #414 from lsst/u/erykoff/uniformity_metrics

U/erykoff/uniformity metrics

pyproject.toml

@@ -40,6 +40,7 @@ dependencies = [
     "skyproj",
     "pyoorb",
     "rubin-scheduler",
+    "pyarrow",
 ]
 
 [project.optional-dependencies]
@@ -68,6 +69,7 @@ run_moving_calc = "rubin_sim.maf.run_moving_calc:run_moving_calc"
 run_moving_fractions = "rubin_sim.maf.run_moving_fractions:run_moving_fractions"
 run_moving_join = "rubin_sim.maf.run_moving_join:run_moving_join"
 scimaf_dir = "rubin_sim.maf.scimaf_dir:scimaf_dir"
+run_selfcal_metric = "rubin_sim.maf.run_selfcal_metric:run_selfcal_metric"
 show_maf = "rubin_sim.maf.show_maf:show_maf"
 make_lsst_obs = "rubin_sim.moving_objects.make_lsst_obs:make_lsst_obs"
 

rubin_sim/maf/maf_contrib/__init__.py

@@ -14,6 +14,7 @@
 from .periodic_star_metric import *
 from .periodic_star_modulation_metric import *
 from .presto_color_kne_pop_metric import *
+from .selfcal_uniformity_metric import *
 from .star_count_mass_metric import *
 from .star_count_metric import *
 from .static_probes_fom_summary_metric import *

rubin_sim/maf/maf_contrib/selfcal_uniformity_metric.py

@@ -0,0 +1,165 @@
+__all__ = ("PhotometricSelfCalUniformityMetric",)
+
+import os
+import time
+
+import healpy as hp
+import numpy as np
+from astropy.coordinates import SkyCoord
+from astropy.table import Table
+from rubin_scheduler.data import get_data_dir
+from rubin_scheduler.utils import healbin
+from scipy.stats import median_abs_deviation
+
+from rubin_sim.maf.metrics import BaseMetric
+from rubin_sim.selfcal import LsqrSolver, OffsetSNR, generate_catalog
+from rubin_sim.selfcal.offsets import OffsetSys
+
+
+def _match(arr1, arr2):
+    st1 = np.argsort(arr1)
+    sub1 = np.searchsorted(arr1, arr2, sorter=st1)
+    if arr2.max() > arr1.max():
+        bad = sub1 == arr1.size
+        sub1[bad] = arr1.size - 1
+
+    (sub2,) = np.where(arr1[st1[sub1]] == arr2)
+    sub1 = st1[sub1[sub2]]
+
+    return sub1, sub2
+
+
+class PhotometricSelfCalUniformityMetric(BaseMetric):
+    def __init__(self, nside_residual=128, highglat_cut=30.0, outlier_nsig=4.0):
+        cols = [
+            "observationid",
+            "fieldra",
+            "fielddec",
+            "fiveSigmaDepth",
+            "rotSkyPos",
+            "filter",
+        ]
+        super().__init__(col=cols, metric_name="PhotometricSelfCalUniformityMetric")
+
+        filename = os.path.join(get_data_dir(), "maf", "monster_stars_uniformity_i15-18_sampled.parquet")
+        self.stars = Table.read(filename)
+
+        # We have to rename dec to decl for the selfcal code.
+        if "dec" in self.stars.dtype.names:
+            self.stars["decl"] = self.stars["dec"]
+
+        self.stars = self.stars.as_array()
+
+        self.nside_residual = nside_residual
+        self.highglat_cut = highglat_cut
+        self.outlier_nsig = outlier_nsig
+
+        self.units = "mmag"
+
+    def run(self, data_slice, slice_point=None):
+        filter_name = data_slice["filter"][0]
+
+        offsets = [OffsetSys(error_sys=0.03), OffsetSNR(lsst_filter=filter_name)]
+
+        visits = np.zeros(
+            len(data_slice),
+            dtype=[
+                ("observationId", "i8"),
+                ("ra", "f8"),
+                ("dec", "f8"),
+                ("fiveSigmaDepth", "f8"),
+                ("rotSkyPos", "f8"),
+            ],
+        )
+        visits["observationId"] = data_slice["observationId"]
+        visits["ra"] = data_slice["fieldRA"]
+        visits["dec"] = data_slice["fieldDec"]
+        visits["fiveSigmaDepth"] = data_slice["fiveSigmaDepth"]
+        visits["rotSkyPos"] = data_slice["rotSkyPos"]
+
+        good_stars = np.isfinite(self.stars[f"{filter_name}mag"])
+        stars = self.stars[good_stars]
+
+        observed_stars = generate_catalog(
+            visits,
+            stars,
+            offsets=offsets,
+            lsst_filter=filter_name,
+            n_patches=16,
+            verbose=False,
+        )
+
+        solver = LsqrSolver(observed_stars)
+
+        print("Starting solver...")
+        t0 = time.time()
+
+        solver.run()
+        fit_patches, fit_stars = solver.return_solution()
+
+        t1 = time.time()
+        dt = (t1 - t0) / 60.0
+        print("runtime= %.1f min" % dt)
+
+        # Trim stars to the ones with solutions.
+        a, b = _match(stars["id"], fit_stars["id"])
+        stars_trimmed = stars[a]
+        fit_stars_trimmed = fit_stars[b]
+
+        # Residuals after fit, removing floating zeropoint
+        resid = stars_trimmed[f"{filter_name}mag"] - fit_stars_trimmed["fit_mag"]
+        resid = resid - np.median(resid)
+
+        resid_map = healbin(
+            stars_trimmed["ra"],
+            stars_trimmed["dec"],
+            resid,
+            self.nside_residual,
+            reduce_func=np.median,
+        )
+
+        (ipring,) = np.where(resid_map > hp.UNSEEN)
+
+        scatter_full = median_abs_deviation(resid_map[ipring], scale="normal")
+
+        # Fraction of (4) sigma outliers.
+        highscat = np.abs(resid_map[ipring]) > self.outlier_nsig * scatter_full
+        outlier_frac_full = highscat.sum() / len(ipring)
+
+        # Cut to high latitude
+        ra, dec = hp.pix2ang(self.nside_residual, ipring, nest=False, lonlat=True)
+        coords = SkyCoord(ra, dec, frame="icrs", unit="deg")
+        b = coords.galactic.b.value
+
+        high_glat = np.abs(b) > self.highglat_cut
+        scatter_highglat = median_abs_deviation(resid_map[ipring[high_glat]], scale="normal")
+
+        # Fraction of (4) sigma outliers.
+        highscat = np.abs(resid_map[ipring[high_glat]]) > self.outlier_nsig * scatter_highglat
+        outlier_frac_highglat = highscat.sum() / len(ipring)
+
+        # Convert to mmag
+        scatter_full = scatter_full * 1000.0
+        scatter_highglat = scatter_highglat * 1000.0
+
+        result = {
+            "scatter_full": scatter_full,
+            "scatter_highglat": scatter_highglat,
+            "outlier_frac_full": outlier_frac_full,
+            "outlier_frac_highglat": outlier_frac_highglat,
+            "uniformity_map": resid_map,
+        }
+
+        return result
+
+    def reduce_scatter_full(self, metric_value):
+        return metric_value["scatter_full"]
+
+    def reduce_scatter_highglat(self, metric_value):
+        return metric_value["scatter_highglat"]
+
+    def reduce_outlier_frac_full(self, metric_value):
+        return metric_value["outlier_frac_full"]
+
+    def reduce_outlier_frac_highglat(self, metric_value):
+        return metric_value["outlier_frac_highglat"]

rubin_sim/maf/run_selfcal_metric.py

@@ -0,0 +1,84 @@
+__all__ = ("run_selfcal_metric",)
+
+import argparse
+import os
+import shutil
+
+import healpy as hp
+import matplotlib
+import numpy as np
+
+matplotlib.use("Agg")
+
+from .db import ResultsDb
+from .maf_contrib import PhotometricSelfCalUniformityMetric
+from .metric_bundles import MetricBundle, MetricBundleGroup
+from .metrics import IdentityMetric
+from .plots import HealpixHistogram, HealpixSkyMap, PlotHandler
+from .slicers import HealpixSlicer, UniSlicer
+
+
+def run_selfcal_metric():
+    """
+    Run the self-calibration metric on one database.
+    """
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--db", type=str, default=None)
+    parser.add_argument(
+        "--no_clobber",
+        dest="no_clobber",
+        default=False,
+        action="store_true",
+        help="Do not remove existing directory outputs",
+    )
+    args = parser.parse_args()
+
+    opsdb = args.db
+    sim_name = os.path.basename(opsdb).replace(".db", "")
+
+    out_dir = sim_name + "_selfcal"
+    if not args.no_clobber:
+        if os.path.isdir(out_dir):
+            shutil.rmtree(out_dir)
+
+    # Set up the metric bundle.
+    map_nside = 128
+    selfcal_metric = PhotometricSelfCalUniformityMetric(nside_residual=map_nside)
+    slicer = UniSlicer()
+    # Exclude DDF visits
+    sql = "note not like '%DD%'"
+    # And run on only year 1 (?)
+    sql += " and night < 366"
+    bundle = MetricBundle(selfcal_metric, slicer, sql, run_name=sim_name, info_label="year 1 no-DD")
+
+    # Set up the resultsDB
+    results_db = ResultsDb(out_dir=out_dir)
+    # Go and run it
+    group = MetricBundleGroup(
+        {"selfcal": bundle}, opsdb, out_dir=out_dir, results_db=results_db, save_early=True
+    )
+    group.run_all(clear_memory=False, plot_now=True)
+
+    # Make plots of the residuals map
+    map_bundle = MetricBundle(
+        IdentityMetric(metric_name="PhotoCal Uniformity"),
+        HealpixSlicer(map_nside),
+        sql,
+        run_name=sim_name,
+        info_label="year 1 no-DD",
+    )
+    tmp_vals = bundle.metric_values[0]["uniformity_map"]
+    tmp_vals = np.where(tmp_vals == hp.UNSEEN, map_bundle.slicer.badval, tmp_vals)
+    map_bundle.metric_values = np.ma.MaskedArray(
+        data=tmp_vals,
+        mask=np.zeros(map_bundle.slicer.shape, "bool"),
+        fill_value=map_bundle.slicer.badval,
+    )
+    map_bundle.write(out_dir=out_dir, results_db=results_db)
+
+    ph = PlotHandler(results_db=results_db, out_dir=out_dir)
+    ph.set_metric_bundles([map_bundle])
+    _ = ph.plot(HealpixSkyMap(), plot_dicts={"color_min": -0.02, "color_max": 0.02})
+    _ = ph.plot(HealpixHistogram(), plot_dicts={"percentile_clip": 99})
+
+    results_db.close()

rubin_sim/selfcal/offsets.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from rubin_sim.maf import m52snr
+from ..maf.utils.astrometry_utils import m52snr
 
 # from lsst.sims.selfcal.clouds.Arma import ArmaSf, Clouds
 

rubin_sim/selfcal/star_tools.py

@@ -1,8 +1,7 @@
 __all__ = ("stars_project", "assign_patches")
 
 import numpy as np
-
-from rubin_sim.utils import gnomonic_project_toxy
+from rubin_scheduler.utils import gnomonic_project_toxy
 
 
 def stars_project(stars, visit):