diff --git a/schedview/compute/__init__.py b/schedview/compute/__init__.py index 1b02e05d..185ebcba 100644 --- a/schedview/compute/__init__.py +++ b/schedview/compute/__init__.py @@ -13,6 +13,9 @@ "compute_maps", "compute_metric_by_visit", "compute_hpix_metric_in_bands", + "compute_scalar_metric_at_one_mjd", + "compute_scalar_metric_at_mjds", + "compute_mixed_scalar_metric", "often_repeated_fields", "count_visits_by_sim", "match_visits_across_sims", @@ -40,7 +43,13 @@ from .survey import compute_maps, make_survey_reward_df try: - from .maf import compute_hpix_metric_in_bands, compute_metric_by_visit + from .maf import ( + compute_hpix_metric_in_bands, + compute_metric_by_visit, + compute_mixed_scalar_metric, + compute_scalar_metric_at_mjds, + compute_scalar_metric_at_one_mjd, + ) except ModuleNotFoundError as e: if not e.args == ("No module named 'rubin_sim'",): raise e diff --git a/schedview/compute/maf.py b/schedview/compute/maf.py index b57c1ef9..c9b3cb4a 100644 --- a/schedview/compute/maf.py +++ b/schedview/compute/maf.py @@ -1,14 +1,24 @@ +import datetime import sqlite3 +from collections.abc import Mapping from pathlib import Path from tempfile import TemporaryDirectory +from typing import Any, Callable, Sequence +import numpy as np import pandas as pd from rubin_scheduler.scheduler.utils import SchemaConverter from rubin_sim import maf +from schedview.dayobs import DayObs from schedview.util import band_column -__all__ = ["compute_metric_by_visit"] +__all__ = [ + "compute_metric_by_visit", + "compute_scalar_metric_at_one_mjd", + "compute_mixed_scalar_metric", + "make_metric_progress_df", +] def _visits_to_opsim(visits, opsim): @@ -42,7 +52,7 @@ def _visits_to_opsim(visits, opsim): merged_opsim.to_sql("observations", con) -def compute_metric(visits, metric_bundle): +def compute_metric(visits, metric_bundle, sqlite=True): """Compute metrics with MAF. Parameters @@ -52,6 +62,10 @@ def compute_metric(visits, metric_bundle): database). metric_bundle : `maf.MetricBundle`, `dict`, or `list` of `maf.MetricBundle` The metric bundle(s) to run. + sqlite : `bool` + Write visits to an sqlite3 database and then read the visits back + from it when computing the metrics. Needed when metric bundles have + constraints. Returns ------- @@ -62,11 +76,15 @@ def compute_metric(visits, metric_bundle): metric_bundles = [metric_bundle] if passed_one_bundle else metric_bundle with TemporaryDirectory() as working_dir: - visits_db = Path(working_dir).joinpath("visits.db").as_posix() - _visits_to_opsim(visits, visits_db) + if sqlite: + visits_db = Path(working_dir).joinpath("visits.db").as_posix() + _visits_to_opsim(visits, visits_db) - bundle_group = maf.MetricBundleGroup(metric_bundles, visits_db, out_dir=working_dir) - bundle_group.run_all() + bundle_group = maf.MetricBundleGroup(metric_bundles, visits_db, out_dir=working_dir) + bundle_group.run_all() + else: + bundle_group = maf.MetricBundleGroup(metric_bundles, None, out_dir=working_dir) + bundle_group.run_current(None, visits.to_records(index=False)) return metric_bundle @@ -81,7 +99,7 @@ def compute_metric_by_visit(visits, metric, constraint=""): database). metric : `rubin_sim.maf.metrics.BaseMetric` The metric to compute. - constraint : `str` + constraint : `str` or `None` The SQL query to filter visits to be used. Returns @@ -95,7 +113,10 @@ def compute_metric_by_visit(visits, metric, constraint=""): slicer = maf.OneDSlicer("observationId", bin_size=1) metric_bundle = maf.MetricBundle(slicer=slicer, metric=metric, constraint=constraint) - compute_metric(visits, metric_bundle) + # If the constraint is set, we need to use sqlite + use_sqlite = constraint is None or len(constraint) > 0 + + compute_metric(visits, metric_bundle, sqlite=use_sqlite) result = pd.Series(metric_bundle.metric_values, index=slicer.slice_points["bins"][:-1].astype(int)) result.index.name = "observationId" return result @@ -125,15 +146,334 @@ def compute_hpix_metric_in_bands(visits, metric, constraint="", nside=32): # Do only the filters we actually used used_bands = visits[band_column(visits)].unique() + # If the constraint is set, we need to use sqlite + use_sqlite = len(constraint) > 0 + bundles = {} for this_band in used_bands: - this_constraint = f"{band_column(visits)} == '{this_band}'" - if len(constraint) > 0: - this_constraint += f" AND {constraint}" + band_visits = visits.query(f"{band_column(visits)} == '{this_band}'") slicer = maf.HealpixSlicer(nside=nside, verbose=False) - bundles[this_band] = maf.MetricBundle(metric, slicer, this_constraint) + bundles[this_band] = maf.MetricBundle(metric, slicer, constraint) + compute_metric(band_visits, bundles[this_band], sqlite=use_sqlite) - compute_metric(visits, bundles) metric_values = {b: bundles[b].metric_values for b in bundles if bundles[b].metric_values is not None} return metric_values + + +def compute_scalar_metric_at_one_mjd( + mjd: float, + visits: pd.DataFrame, + slicer: maf.BaseSlicer, + metric: maf.BaseMetric, + summary_metric: maf.BaseMetric | None = None, + run_name: str | None = None, + mjd_column: str = "observationStartMJD", +) -> dict[str, float]: + """Compute a scalar MAF metric at a specific modified Julian date. + + Parameters + ---------- + mjd : `float` + The modified Julian date to use as the cutoff for visits. + visits : `pandas.DataFrame` + The DataFrame of visits. + slicer : `rubin_sim.maf.slicers.BaseSlicer` + The slicer to use for computing the metric. + metric : `rubin_sim.maf.metrics.BaseMetric` + The metric to compute. + summary_metric : `rubin_sim.maf.metrics.BaseMetric` or `None`, optional + The summary metric to compute. + run_name : `str` or `None`, optional + The name to use for the run. + mjd_column : `str`, optional + The name of the column containing the MJD values. + Default is "observationStartMJD". + + Returns + ------- + metric_values : `dict[str, float]` + A dictionary with the metric name as key and the computed scalar value + as value. + + Raises + ------ + ValueError + If no visits are found before the specified MJD. + """ + visits = visits.loc[visits[mjd_column] < mjd, :] + + if len(visits) == 0: + raise ValueError("No visits") + + if run_name is None: + run_name = "Run" + datetime.datetime.now().isoformat() + + if summary_metric is not None: + bundle = maf.MetricBundle( + metric, + slicer, + summary_metrics=[summary_metric], + run_name=run_name, + ) + else: + bundle = maf.MetricBundle( + metric, + slicer, + run_name=run_name, + ) + + compute_metric(visits, bundle, sqlite=False) + + if summary_metric is None: + metric_name = metric.name + metric_values = bundle.metric_values + else: + bundle.compute_summary_stats() + metric_name = summary_metric.name + metric_values = ( + tuple(bundle.summary_values.values()) + if isinstance(bundle.summary_values, Mapping) + else bundle.summary_values + ) + + assert len(metric_values) == 1 + return {metric_name: metric_values[0]} + + +def compute_scalar_metric_at_mjds( + mjds: Sequence[float], + *args: Any, + **kwargs: Any, +) -> pd.Series: + """Compute a scalar MAF metric at multiple modified Julian dates. + + Parameters + ---------- + mjds : `sequence` of `float` + The MJDs to use as cutoffs for visits. + *args : `Any` + Positional arguments to pass to `compute_scalar_metric_at_one_mjd`. + **kwargs : `Any` + Keyword arguments to pass to `compute_scalar_metric_at_one_mjd`. + + Returns + ------- + metric_values : `pandas.Series` + A Series with the computed metric values, indexed by the MJDs used. + + See Also + -------- + compute_scalar_metric_at_one_mjd + The function used to compute the metric at each individual MJD. + """ + metric_values = [] + name = None + mjds_with_data = [] + for mjd in mjds: + try: + metric_value_dict = compute_scalar_metric_at_one_mjd(mjd, *args, **kwargs) + except ValueError: + continue + these_keys = tuple(metric_value_dict.keys()) + assert len(these_keys) == 1 + if name is None: + name = these_keys[0] + else: + assert these_keys[0] == name + mjds_with_data.append(mjd) + metric_values.append(metric_value_dict[name]) + + metric_values = pd.Series(metric_values, index=mjds_with_data, name=name) + return metric_values + + +def compute_mixed_scalar_metric( + start_visits: pd.DataFrame, + end_visits: pd.DataFrame, + transition_mjds: Sequence[float], + *args, + mjd_column: str = "observationStartMJD", + **kwargs, +) -> pd.Series: + """Compute a scalar MAF metric at multiple transition MJDs using a + visits from `start_visits` up to each transition MJD, and visits + from `end_visits` after + + Parameters + ---------- + start_visits : `pandas.DataFrame` + DataFrame of visits to use before each transition MJD. + end_visits : `pandas.DataFrame` + DataFrame of visits to use after each transition MJD. + transition_mjds : `sequence` of `float` + MJDs that define the transition points for combining visits. + mjd_column : `str`, optional + The name of the column containing MJD values. + Default is "observationStartMJD". + *args : `Any` + Positional arguments to pass to `compute_scalar_metric_at_one_mjd`. + **kwargs : `Any` + Keyword arguments to pass to `compute_scalar_metric_at_one_mjd`. + + Returns + ------- + metric_values : `pandas.Series` + A Series with the computed metric values, indexed by the transition + MJDs used. + """ + mjds_with_data = [] + metric_values = [] + name = None + for transition_mjd in transition_mjds: + visits = pd.concat( + ( + start_visits.loc[start_visits[mjd_column] <= transition_mjd, :].dropna( + axis="columns", how="all" + ), + end_visits.loc[transition_mjd < end_visits[mjd_column], :].dropna(axis="columns", how="all"), + ) + ) + try: + metric_value_dict = compute_scalar_metric_at_one_mjd(*args, visits=visits, **kwargs) + except ValueError: + continue + + these_keys = tuple(metric_value_dict.keys()) + assert len(these_keys) == 1 + if name is None: + name = these_keys[0] + else: + assert these_keys[0] == name + mjds_with_data.append(transition_mjd) + metric_values.append(metric_value_dict[name]) + + metric_values = pd.Series(metric_values, index=mjds_with_data, name=name) + return metric_values + + +def make_metric_progress_df( + completed_visits: pd.DataFrame, + baseline_visits: pd.DataFrame, + start_dayobs: datetime.date | int | str | DayObs, + last_completed_dayobs: datetime.date | int | str | DayObs, + extrapolation_dayobs: datetime.date | int | str | DayObs, + slicer_factory: Callable[[], maf.BaseSlicer], + metric_factory: Callable[[], maf.BaseMetric], + summary_metric_factory: Callable[[], maf.BaseMetric] | None = None, + freq: dict | str = "MS", +) -> pd.DataFrame: + """Compute metric values over time for completed, baseline, and mixed + (chimera) sets of visits. + + Parameters + ---------- + completed_visits : `pandas.DataFrame` + Completed visits. + baseline_visits : `pandas.DataFrame` + Visits from the baseline survey strategy. + start_dayobs : `datetime.date`, `int`, `str`, or `DayObs` + The start date for the time range to compute metrics. + last_completed_dayobs : `datetime.date`, `int`, `str`, or `DayObs` + The last dayobs completed. + extrapolation_dayobs : `datetime.date`, `int`, `str`, or `DayObs` + The end date for the time range to compute metrics. + slicer_factory : `callable` + A callable that returns a `maf.BaseSlicer` instance to be used for + computing metrics. + metric_factory : `callable` + A callable that returns a `maf.BaseMetric` instance to be used for + computing metrics. + summary_metric_factory : `callable` or `None`, optional + A callable that returns a `maf.BaseMetric` instance to be used as a + summary metric. + If `None`, no summary metric is computed. + freq : `dict` or `str`, optional + Frequency for computing metrics, as used in the ``freq`` argument + to `pandas.date_range`. If a string, it's used for both + completed and future dates. If a dict, it should have 'completed' + and 'future' keys. Default is "MS" (month start). + + Returns + ------- + metric_values : `pandas.DataFrame` + A DataFrame with columns: + - 'date': The date for each time point + - 'snapshot': Metric values computed using completed visits + - 'baseline': Metric values computed using baseline visits + - 'chimera': Metric values computed using a combination + The index is the MJD values used for computation. + """ + # Be flexible in how we accept the start and end dates. + start_dayobs = DayObs.from_date(start_dayobs) + last_completed_dayobs = DayObs.from_date(last_completed_dayobs) + end_dayobs = DayObs.from_date(extrapolation_dayobs) + # The asserts make type checkers happy. + assert isinstance(start_dayobs, DayObs) + assert isinstance(last_completed_dayobs, DayObs) + assert isinstance(end_dayobs, DayObs) + + # If we only get one value for freq, use it for both completed + # and future dates. + if isinstance(freq, str): + freq = {"completed": freq, "future": freq} + + timeframes = set(["completed", "future"]) + assert isinstance(freq, dict) + assert set(freq.keys()) == timeframes, "freq must be a dict with 'completed' and 'future' keys" + + timestamps = {} + dayobs_end_mjds = {} + dayobs_start_mjds = {} + for timeframe in timeframes: + timestamps[timeframe] = pd.date_range( + start=start_dayobs.date, end=end_dayobs.date, freq=freq[timeframe] + ) + dayobs_end_mjds[timeframe] = np.array( + [DayObs.from_date(d).end.mjd for d in timestamps[timeframe].date] + ) + dayobs_start_mjds[timeframe] = np.array( + [DayObs.from_date(d).start.mjd for d in timestamps[timeframe].date] + ) + + max_completed_mjd = last_completed_dayobs.end.mjd + dayobs_end_mjds["completed"] = dayobs_end_mjds["completed"][ + dayobs_start_mjds["completed"] < max_completed_mjd + ] + dayobs_end_mjds["future"] = dayobs_end_mjds["future"][dayobs_start_mjds["future"] > max_completed_mjd] + dayobs_end_mjds["all"] = np.concatenate([dayobs_end_mjds["completed"], dayobs_end_mjds["future"]]) + + dates = pd.to_datetime([DayObs.from_time(end_mjd - 0.1).date for end_mjd in dayobs_end_mjds["all"]]) + jds = np.array([DayObs(d).jd for d in dates]) + metric_values = pd.DataFrame( + {"date": dates, "jd": jds}, index=pd.Index(dayobs_end_mjds["all"], name="mjd") + ) + + # create a helper function to make it easier to avoid trying to call + # summary_metric_factory when there isn't one. + def maf_args(): + args = {"slicer": slicer_factory(), "metric": metric_factory()} + if summary_metric_factory is not None: + assert isinstance(summary_metric_factory, Callable) + args["summary_metric"] = summary_metric_factory() + return args + + metric_values["snapshot"] = compute_scalar_metric_at_mjds( + dayobs_end_mjds["completed"], visits=completed_visits, **maf_args() + ) + + metric_values["baseline"] = compute_scalar_metric_at_mjds( + dayobs_end_mjds["all"], visits=baseline_visits, **maf_args() + ) + + metric_values["chimera"] = compute_mixed_scalar_metric( + completed_visits, + baseline_visits, + transition_mjds=dayobs_end_mjds["completed"], + mjd=end_dayobs.mjd, + **maf_args(), + ) + + metric_values.set_index("jd", inplace=True) + + return metric_values diff --git a/schedview/compute/visits.py b/schedview/compute/visits.py index dd55c878..11ff9c09 100644 --- a/schedview/compute/visits.py +++ b/schedview/compute/visits.py @@ -1,7 +1,10 @@ import datetime import warnings +import astropy.units as u import numpy as np +import pandas as pd +from astropy.coordinates import SkyCoord, search_around_sky from astropy.time import Time from rubin_scheduler.site_models import SeeingModel @@ -342,3 +345,71 @@ def accum_teff_by_night(visits): """ stats_by_night = accum_stats_by_target_band_night(visits) return stats_by_night.loc[:, "t_eff"] + + +def match_visits_to_pointings( + visits: pd.DataFrame, + pointings: dict, + ra_col: str = "s_ra", + decl_col: str = "s_dec", + name_col: str = "pointing_name", + match_radius: float = 1.75, +) -> pd.DataFrame: + """Match visits to pointings based on coordinates. + + Parameters + ---------- + visits : `pd.DataFrame` + DataFrame containing visit data with equatorial coordinates. + pointings : `dict` + Dictionary of pointings where keys are pointing names and values are + coordinate tuples (ra, dec) in degrees. + ra_col : `str`, optional + Name of the column containing right ascension values in the visits + DataFrame. Default is "s_ra". + decl_col : `str`, optional + Name of the column containing declination values in the visits + DataFrame. Default is "s_dec". + name_col : `str`, optional + Name of the column to be added to the output DataFrame to identify + which pointing each visit matches to. Default is "pointing_name". + match_radius : `float`, optional + Matching radius in degrees for associating visits with pointings. + Default is 1.75 degrees. + + Returns + ------- + pointing_visits : `pd.DataFrame` + DataFrame containing the original visits DataFrame with an additional + column (named by ``name_col``) identifying which pointing each visit + matches to. The returned DataFrame maintains the original visit data + but is filtered to only include visits that match to at least one + pointing. The DataFrame will contain one row for every pointing/visit + match, so selecting on a given pointing will result in a DataFrame + with all visits within match_radius of that pointing. If a visits + covers multiple pointings, it can appear multiple times in the + DataFrame. + + """ + + pointings_df = pd.DataFrame(pointings).T + pointings_df.columns = pd.Index(["ra", "decl"]) + pointing_coords = SkyCoord( + ra=pointings_df.ra.values * u.deg, dec=pointings_df.decl.values * u.deg, frame="icrs" + ) + + visit_centers = SkyCoord( + ra=visits[ra_col].values * u.deg, dec=visits[decl_col].values * u.deg, frame="icrs" + ) + + pointing_matches = search_around_sky(pointing_coords, visit_centers, match_radius * u.deg) + + visit_match_dfs = {} + for pointing_idx, pointing_name in enumerate(pointings): + visit_idx = pointing_matches[1][pointing_matches[0] == pointing_idx] + visit_match_dfs[pointing_name] = visits.iloc[visit_idx, :].copy() + visit_match_dfs[pointing_name][name_col] = pointing_name + + pointing_visits = pd.concat([visit_match_dfs[n] for n in visit_match_dfs]) + + return pointing_visits diff --git a/schedview/dayobs.py b/schedview/dayobs.py index 382691f1..7196174d 100644 --- a/schedview/dayobs.py +++ b/schedview/dayobs.py @@ -199,6 +199,12 @@ def from_time( dayobs_datetime: datetime.datetime = maybe_datetime else: raise ValueError(f"{cls.__name__} currently only accepts scalars.") + case np.datetime64(): + maybe_datetime = Time(arg).datetime + if isinstance(maybe_datetime, datetime.datetime): + dayobs_datetime: datetime.datetime = maybe_datetime + else: + raise ValueError(f"{cls.__name__} currently only accepts scalars.") case Time(): maybe_datetime = arg.datetime if isinstance(maybe_datetime, datetime.datetime): diff --git a/schedview/plot/__init__.py b/schedview/plot/__init__.py index 99cb0de4..60025125 100644 --- a/schedview/plot/__init__.py +++ b/schedview/plot/__init__.py @@ -27,6 +27,8 @@ "create_overhead_histogram", "plot_overhead_vs_slew_distance", "PLOT_BAND_COLORS", + "LIGHT_PLOT_BAND_COLORS", + "LIGHT_EXTRA_COLORS", "make_band_cmap", "create_cadence_plot", "create_visit_table", @@ -35,10 +37,11 @@ "make_timeline_scatterplots", "make_html_table_of_sim_archive_metadata", "mpl_fig_to_html", + "make_metric_line_plots", ] from .cadence import create_cadence_plot -from .colors import PLOT_BAND_COLORS, make_band_cmap +from .colors import LIGHT_EXTRA_COLORS, LIGHT_PLOT_BAND_COLORS, PLOT_BAND_COLORS, make_band_cmap from .multisim import generate_sim_indicators, overplot_kernel_density_estimates from .nightbf import plot_infeasible, plot_rewards from .nightly import plot_airmass_vs_time, plot_alt_vs_time, plot_polar_alt_az @@ -48,6 +51,7 @@ create_survey_visit_summary_table, plot_overhead_vs_slew_distance, ) +from .progress import make_metric_line_plots from .rewards import ( area_timeline_for_tier, create_survey_reward_plot, diff --git a/schedview/plot/colors.py b/schedview/plot/colors.py index 222999e5..fac615b2 100644 --- a/schedview/plot/colors.py +++ b/schedview/plot/colors.py @@ -1,73 +1,106 @@ +from typing import Any, Callable, Dict + import bokeh.core.property import bokeh.transform -# Follow RTN-045 for mapping filters to plot colors - -PLOT_BAND_COLORS = { - "u": "#1600ea", - "g": "#31de1f", - "r": "#b52626", - "i": "#370201", - "z": "#ba52ff", - "y": "#61a2b3", -} - -# Use -# glasbey.extend_palette( -# PLOT_BAND_COLORS.values(), -# palette_size=48, -# colorblind_safe=True, -# grid_size=256 -# ) -# to generate extra colors, attempting to get distinguishable ones. -EXTRA_COLORS = [ - "#007000", - "#ff8c72", - "#006169", - "#000073", - "#9461ff", - "#002e00", - "#d8a2ff", - "#26af00", - "#dc1a25", - "#36164c", - "#10d9aa", - "#48936f", - "#0113ff", - "#84010e", - "#9d7cb5", - "#01593c", - "#7c3ca8", - "#e2b200", - "#015cff", - "#924e47", - "#00b2fd", - "#81c9d3", - "#ff5701", - "#2b8683", - "#00457c", - "#c87168", - "#0ed767", - "#0080aa", - "#7c7f31", - "#3702a6", - "#364801", - "#983e00", - "#5dae9d", - "#540005", - "#79a762", - "#c073f8", - "#0d7849", - "#280162", - "#73a9c9", - "#013ed6", - "#5a3c71", - "#9804f1", +try: + from glasbey import extend_palette +except ModuleNotFoundError: + + def extend_palette(*args, **kwargs) -> Any: + raise NotImplementedError + + +extend_palette: Callable[..., Any] + +try: + from lsst.utils.plotting import get_multiband_plot_colors +except ModuleNotFoundError: + + def get_multiband_plot_colors(*arg, **kwargs) -> Dict[str, str]: + raise NotImplementedError + + +get_multiband_plot_colors: Callable[..., Any] + +BANDS: tuple[str] = ("u", "g", "r", "i", "z", "y") + +# Palettes from RTN-045 as of 2026-02-17, supplemented so we +# have extra colors for non-standard "bands" (e.g. pinhole) using +# glasbey.extend_palette with kwargs: +# palette_size=12, colorblind_safe=True, grid_size=256 +# +DEFAULT_DARK_COLORS: list[str] = [ + "#0c71ff", + "#49be61", + "#c61c00", + "#ffc200", + "#f341a2", + "#5d0000", + "#0000b6", + "#04e9ff", + "#006200", + "#750065", + "#e74eff", + "#ff5f00", +] +DEFAULT_BAND_COLORS: Dict[str, str] = {b: DEFAULT_DARK_COLORS[i] for i, b in enumerate(BANDS)} + +DEFAULT_LIGHT_COLORS: list[str] = [ + "#3eb7ff", + "#30c39f", + "#ff7e00", + "#2af5ff", + "#a7f9c1", + "#fdc900", + "#85b972", + "#12c5c8", + "#f2e585", + "#fba96c", + "#23ffdf", + "#7ccb2b", ] +DEFAULT_LIGHT_BAND_COLORS: Dict[str, str] = {b: DEFAULT_LIGHT_COLORS[i] for i, b in enumerate(BANDS)} + +PLOT_BAND_COLORS: Dict[str, str] +try: + PLOT_BAND_COLORS = get_multiband_plot_colors() +except NotImplementedError: + PLOT_BAND_COLORS = DEFAULT_BAND_COLORS + +LIGHT_PLOT_BAND_COLORS: Dict[str, str] +try: + LIGHT_PLOT_BAND_COLORS = get_multiband_plot_colors(dark_background=True) +except NotImplementedError: + LIGHT_PLOT_BAND_COLORS = DEFAULT_LIGHT_BAND_COLORS + +EXTRA_COLORS: list[str] +if DEFAULT_DARK_COLORS[: len(PLOT_BAND_COLORS)] == list(PLOT_BAND_COLORS.values()): + EXTRA_COLORS = DEFAULT_DARK_COLORS +else: + try: + EXTRA_COLORS = extend_palette( + PLOT_BAND_COLORS.values(), palette_size=12, colorblind_safe=True, grid_size=256 + ) + except NotImplementedError: + EXTRA_COLORS = list(PLOT_BAND_COLORS.values()) + DEFAULT_DARK_COLORS[len(PLOT_BAND_COLORS) :] + +LIGHT_EXTRA_COLORS: list[str] +if DEFAULT_LIGHT_COLORS[: len(LIGHT_PLOT_BAND_COLORS)] == list(LIGHT_PLOT_BAND_COLORS.values()): + LIGHT_EXTRA_COLORS = DEFAULT_LIGHT_COLORS +else: + try: + LIGHT_EXTRA_COLORS = extend_palette( + LIGHT_PLOT_BAND_COLORS.values(), palette_size=12, colorblind_safe=True, grid_size=256 + ) + except NotImplementedError: + LIGHT_EXTRA_COLORS = ( + list(LIGHT_PLOT_BAND_COLORS.values()) + DEFAULT_LIGHT_COLORS[len(LIGHT_PLOT_BAND_COLORS) :] + ) def make_band_cmap( - field_name="band", bands=("u", "g", "r", "i", "z", "y") + field_name="band", bands=("u", "g", "r", "i", "z", "y"), light=False ) -> bokeh.core.property.vectorization.Field: """Make a bokeh cmap transform for bands @@ -77,20 +110,25 @@ def make_band_cmap( Name of field with the band value. bands : `list` A full list of bands that need colors. + light : `bool` + Use light palette Returns ------- cmap : `bokeh.core.property.vectorization.Field` The bokeh color map. """ + band_colors = LIGHT_PLOT_BAND_COLORS if light else PLOT_BAND_COLORS + extra_colors = LIGHT_EXTRA_COLORS if light else EXTRA_COLORS + # Always assign standard colors to standard bands - assigned_colors, assigned_bands = list(PLOT_BAND_COLORS.values()), list(PLOT_BAND_COLORS.keys()) + assigned_colors, assigned_bands = list(band_colors.values()), list(band_colors.keys()) # Go through any we miss and assign extra colors to extra bands. for extra_index, band in enumerate(set(bands)): if band not in assigned_bands: assigned_bands.append(band) - assigned_colors.append(EXTRA_COLORS[extra_index % len(EXTRA_COLORS)]) + assigned_colors.append(extra_colors[extra_index % len(extra_colors)]) cmap = bokeh.transform.factor_cmap(field_name, tuple(assigned_colors), tuple(assigned_bands)) return cmap diff --git a/schedview/plot/monthlyhg.py b/schedview/plot/monthlyhg.py new file mode 100644 index 00000000..a53a478f --- /dev/null +++ b/schedview/plot/monthlyhg.py @@ -0,0 +1,66 @@ +import calendar +import datetime +from typing import Callable + +import matplotlib as mpl +import matplotlib.pyplot as plt +import pandas as pd +from rubin_sim import maf + +from .util import mpl_fig_to_html + + +def plot_collapsed_monthly_hourglass_metric( + metric_bundle: maf.MetricBundle, + name: str, + first_date: datetime.date | str, + last_date: datetime.date | str, + plotter_factory: Callable[[int, int], maf.BasePlotter] = maf.plots.MonthHourglassPlot, +) -> str: + """Generate a collapsible HTML string of hourglass plots. + + Parameters + ---------- + metric_bundle : Any + The ``MetricBundle`` to run. + name : str + Base name used in the plot titles. + first_date : `datetime.date` or `str` + Starting date for the range. + last_date : `datetime.date` or `str` + End date for the range, inclusive. + plotter_factory : Callable[[int, int], maf.BasePlotter], optional + Factory that creates a MAF plotter for a given month and year. The + default creates a ``MonthHourglassPlot``. + + Returns + ------- + html_figures : `str` + HTML string containing the collapsed hourglass figures. + """ + + # Use plt.ioff to stop the jupyter notebook from showing the fig + # natively, so we can use the collapseable html instead of rather + # than in addition to it. + with plt.ioff(): + # Set savefig.bbox to tight to avoid clipping off the legend. + with mpl.rc_context({"savefig.bbox": "tight"}): + figs_html = [] + + # If we pass a list of plot_funcs to the MetricBundle, they all + # get the same plot_dict, an so same title. We want different + # titles for each month, so iterate here: + for month_timestamp in pd.date_range( + pd.Timestamp(first_date).replace(day=1), last_date, freq="MS" + ): + year, month = month_timestamp.year, month_timestamp.month + title = f"{name} for {calendar.month_name[month]}, {year}" + metric_bundle.plot_dict["title"] = title + metric_bundle.plot_funcs = [plotter_factory(month, year)] + these_figs = metric_bundle.plot() + assert len(these_figs) == 1 + this_fig = next(iter(these_figs.values())) + figs_html.append(mpl_fig_to_html(this_fig, summary=f"Hourglass of {title}")) + + all_figs_html = "\n".join(figs_html) + return all_figs_html diff --git a/schedview/plot/progress.py b/schedview/plot/progress.py new file mode 100644 index 00000000..31fad747 --- /dev/null +++ b/schedview/plot/progress.py @@ -0,0 +1,90 @@ +from collections.abc import Mapping +from types import MappingProxyType + +import bokeh.layouts +import bokeh.models +import bokeh.plotting +import pandas as pd + +__all__ = [ + "make_metric_line_plots", +] + + +SNAPSHOT_SCALAR_METRIC_FIGURE_NAME = "snapshot_scalar_metric_figure" +EXTRAPOLATED_SCALAR_METRIC_FIGURE_NAME = "extrapolated_scalar_metric_figure" + + +def make_metric_line_plots( + metric_name: str, + metric_values: pd.DataFrame, + snapshot_line_kwargs: Mapping = MappingProxyType({"color": "black"}), + snapshot_scatter_kwargs: Mapping = MappingProxyType({"color": "black"}), + baseline_line_kwargs: Mapping = MappingProxyType({"color": "lightgreen", "width": 5}), + baseline_ray_kwargs: Mapping = MappingProxyType({"color": "lightgreen", "width": 5}), + chimera_line_kwargs: Mapping = MappingProxyType({"color": "black"}), + chimera_scatter_kwargs: Mapping = MappingProxyType({"color": "black"}), +) -> bokeh.layouts.Row: + """Create Bokeh line and scatter plots for a scalar metric. + + Parameters + ---------- + metric_name : str + Used in the titles of the generated figures. + metric_values : pd.DataFrame + DataFrame containing a datetime column named ``date`` and metric + columns ``baseline``, ``snapshot`` and ``chimera``. The index is + expected to be a ``pd.DatetimeIndex``. + snapshot_line_kwargs : `Mapping` + Keyword args passed to ``bokeh.figure.line`` for the snapshot plot. + snapshot_scatter_kwargs : `Mapping` + Keyword args passed to ``bokeh.figure.scatter`` for the snapshot plot. + baseline_line_kwargs : `Mapping` + Keyword args passed to ``bokeh.figure.line`` for the baseline line. + baseline_ray_kwargs : `Mapping` + Keyword args passod to ``bokeh.figure.ray`` for the baseline + extrpolated metric value. + chimera_line_kwargs : `Mapping` + Keyword args passed to ``bokeh.figure.line`` for the chimera metric. + chimera_scatter_kwargs : `Mapping` + Keyword args passed to ``bokeh.figure.scatter`` for the chimera metric. + + + Returns + ------- + showable : `bokeh.layouts.Row` + A Bokeh layout row containing the two figures. + """ + + metric_values_ds = bokeh.models.ColumnDataSource(metric_values) + + # Build the left plot showing instantaneous metric values. + metric_at_date_fig = bokeh.plotting.figure( + title=f"{metric_name} at date", + x_axis_label="Date", + y_axis_label="Metric value", + x_axis_type="datetime", + name=SNAPSHOT_SCALAR_METRIC_FIGURE_NAME, + ) + metric_at_date_fig.line(x="date", y="baseline", source=metric_values_ds, **baseline_line_kwargs) + metric_at_date_fig.line(x="date", y="snapshot", source=metric_values_ds, **snapshot_line_kwargs) + metric_at_date_fig.scatter(x="date", y="snapshot", source=metric_values_ds, **snapshot_scatter_kwargs) + + # Build the right plot showing the metric values extrapolated to the + # target date. + baseline_final_depth = metric_values.loc[metric_values.index.max(), "baseline"] + chimera_metric_fig = bokeh.plotting.figure( + title=f"{metric_name} extrapolated from date to end with baseline", + x_axis_label="Date", + y_axis_label="Metric value", + x_axis_type="datetime", + name=EXTRAPOLATED_SCALAR_METRIC_FIGURE_NAME, + ) + chimera_metric_fig.ray( + x=metric_values["date"].min(), y=baseline_final_depth, length=0, angle=0, **baseline_ray_kwargs + ) + chimera_metric_fig.line(x="date", y="chimera", source=metric_values_ds, **chimera_line_kwargs) + chimera_metric_fig.scatter(x="date", y="chimera", source=metric_values_ds, **chimera_scatter_kwargs) + + metric_figs = bokeh.layouts.row(metric_at_date_fig, chimera_metric_fig) + return metric_figs diff --git a/schedview/plot/util.py b/schedview/plot/util.py index 7938d68d..9a092d2f 100644 --- a/schedview/plot/util.py +++ b/schedview/plot/util.py @@ -7,7 +7,9 @@ def mpl_fig_to_html( - fig: mpl.figure.Figure, template: str = '' + fig: mpl.figure.Figure, + template: str = '', + summary: str = "", ) -> str: """Convert a Matplotlib figure to an HTML ```` tag. @@ -19,6 +21,9 @@ def mpl_fig_to_html( A format string that will receive the base64‑encoded PNG data. The default template embeds the image in a data URI inside an ```` element. + summary: `str`, option + Embed the figure in a collapseable html element with this summary, + if not the empty string. Defaults to the empty string. Returns ------- @@ -30,4 +35,8 @@ def mpl_fig_to_html( byte_buffer.seek(0) png_base64 = base64.b64encode(byte_buffer.read()).decode("utf-8") fig_html = template.format(png_base64=png_base64) + + if len(summary) > 0: + fig_html = f"
{summary}{fig_html}
" + return fig_html diff --git a/schedview/reports.py b/schedview/reports.py index d677aac3..f1a883f4 100644 --- a/schedview/reports.py +++ b/schedview/reports.py @@ -11,7 +11,7 @@ def find_reports( report_dir: str = "/sdf/data/rubin/shared/scheduler/reports", url_base: str = "https://usdf-rsp-int.slac.stanford.edu/schedview-static-pages", ) -> pd.DataFrame: - """Find staticic schedview reports by walking the report directory. + """Find static schedview reports by walking the report directory. Parameters ---------- @@ -82,7 +82,8 @@ def find_reports( def make_report_link_table( - reports: pd.DataFrame, report_columns=("prenight", "multiprenight", "nightsum", "compareprenight") + reports: pd.DataFrame, + report_columns=("prenight", "multiprenight", "nightsum", "compareprenight", "preprogress"), ) -> str: """Generate an html table of links to reports. diff --git a/tests/test_compute_maf.py b/tests/test_compute_maf.py index 12cb2ab3..9bd4f414 100644 --- a/tests/test_compute_maf.py +++ b/tests/test_compute_maf.py @@ -1,27 +1,47 @@ import unittest +from functools import partial import numpy as np -from rubin_scheduler.utils import SURVEY_START_MJD +import pandas as pd +from rubin_sim import maf from rubin_sim.data import get_baseline +from schedview import DayObs from schedview.collect import read_opsim +from schedview.compute.maf import make_metric_progress_df try: from rubin_sim import maf HAVE_MAF = True - from schedview.compute import compute_hpix_metric_in_bands, compute_metric_by_visit + from schedview.compute import ( + compute_hpix_metric_in_bands, + compute_metric_by_visit, + compute_mixed_scalar_metric, + compute_scalar_metric_at_mjds, + compute_scalar_metric_at_one_mjd, + ) except ModuleNotFoundError: HAVE_MAF = False class TestComputeMAF(unittest.TestCase): + @classmethod + def setUpClass(cls): + # If maf is not installed, we are not doing any tests, so + # do not bother reading cls.visits. + if HAVE_MAF: + cls.visits = read_opsim(get_baseline()) + cls.start_dayobs = DayObs.from_time(cls.visits.observationStartMJD.min()) + cls.start_mjd = cls.start_dayobs.start.mjd + else: + cls.visits = None + @unittest.skipUnless(HAVE_MAF, "No maf installation") def test_compute_metric_by_visit(self): - visits = read_opsim(get_baseline()) - mjd_start = SURVEY_START_MJD - constraint = f"observationStartMjd BETWEEN {mjd_start + 0.5} AND {mjd_start+1.5}" + visits = self.visits + constraint = f"observationStartMjd BETWEEN {self.start_mjd} AND {self.start_mjd+1}" metric = maf.SumMetric(col="t_eff", metric_name="Total Teff") values = compute_metric_by_visit(visits, metric, constraint=constraint) self.assertGreater(len(values), 10) @@ -30,14 +50,170 @@ def test_compute_metric_by_visit(self): @unittest.skipUnless(HAVE_MAF, "No maf installation") def test_compute_hpix_metric_in_bands(self): - visits = read_opsim(get_baseline()) - mjd_start = SURVEY_START_MJD - constraint = f"observationStartMjd BETWEEN {mjd_start+0.5} AND {mjd_start+1.5}" + visits = self.visits.query(f"observationStartMJD < {self.start_mjd + 1}") metric = maf.SumMetric(col="t_eff", metric_name="Total Teff") - values = compute_hpix_metric_in_bands(visits, metric, constraint=constraint) + values = compute_hpix_metric_in_bands(visits, metric) self.assertGreater(len(values.keys()), 1) for band in values.keys(): self.assertTrue(band in "urgizy") if len(values[band]) > 0: self.assertGreater(np.min(values[band]), 0.0) self.assertLess(np.max(values[band]), 300) + + @unittest.skipUnless(HAVE_MAF, "No maf installation") + def test_compute_scalar_metric_at_one_mjd(self): + visits = self.visits + # Use a date that should have visits + mjd = self.start_mjd + 5.0 + metric = maf.SumMetric(col="t_eff", metric_name="Total Teff") + slicer = maf.UniSlicer() + + result = compute_scalar_metric_at_one_mjd(mjd, visits, slicer, metric) + + self.assertIsInstance(result, dict) + self.assertIn("Total Teff", result) + self.assertIsInstance(result["Total Teff"], (int, float)) + self.assertGreaterEqual(result["Total Teff"], 0.0) + + @unittest.skipUnless(HAVE_MAF, "No maf installation") + def test_compute_scalar_metric_at_one_mjd_with_summary_metric(self): + visits = self.visits + # Use a date that should have visits + mjd = self.start_mjd + 5.0 + metric = maf.CountMetric(col="t_eff", metric_name="Count Teff") + slicer = maf.UniSlicer() + summary_metric = maf.SumMetric(col="t_eff", metric_name="Sum Teff") + + result = compute_scalar_metric_at_one_mjd(mjd, visits, slicer, metric, summary_metric=summary_metric) + + self.assertIsInstance(result, dict) + self.assertIn("Sum Teff", result) + self.assertIsInstance(result["Sum Teff"], (int, float)) + self.assertGreaterEqual(result["Sum Teff"], 0.0) + + @unittest.skipUnless(HAVE_MAF, "No maf installation") + def test_compute_scalar_metric_at_one_mjd_no_visits(self): + visits = self.visits + # Use a date that should have no visits + mjd = self.start_mjd - 10.0 + metric = maf.SumMetric(col="t_eff", metric_name="Total Teff") + slicer = maf.UniSlicer() + + with self.assertRaises(ValueError): + compute_scalar_metric_at_one_mjd(mjd, visits, slicer, metric) + + @unittest.skipUnless(HAVE_MAF, "No maf installation") + def test_compute_scalar_metric_at_mjds(self): + visits = self.visits + # Use dates where some have visits and some don't + mjds = self.start_mjd + np.arange(5) + metric = maf.SumMetric(col="t_eff", metric_name="Total Teff") + slicer = maf.UniSlicer() + + result = compute_scalar_metric_at_mjds(mjds, visits, slicer, metric) + + self.assertIsInstance(result, pd.Series) + self.assertLessEqual(len(result), len(mjds[mjds > 0])) + self.assertEqual(result.name, "Total Teff") + + # All values should be numeric and non-negative + for value in result.values: + self.assertIsInstance(value, (int, float)) + self.assertGreaterEqual(value, 0.0) + + @unittest.skipUnless(HAVE_MAF, "No maf installation") + def test_compute_mixed_scalar_metric(self): + + # Use the first month or so of the baseline as the test end visits + end_visits = self.visits[self.visits.observationStartMJD < self.start_mjd + 30] + + # Modify the end visits to create a sample start visits. + start_visits = end_visits.sample(n=100, random_state=42).copy() + + # Shift observationStartDate by random amounts less than 0.0001 days, + # which is just under 9 seconds. + np.random.seed(42) + start_visits["observationStartMJD"] = start_visits["observationStartMJD"] + np.random.uniform( + -0.0001, 0.0001, len(start_visits) + ) + + transition_mjds = self.start_mjd + np.arange(1, 6) * 2 + + metric = maf.SumMetric(col="t_eff", metric_name="Total Teff") + slicer = maf.UniSlicer() + + result = compute_mixed_scalar_metric( + start_visits, end_visits, transition_mjds, mjd=self.start_mjd + 30, slicer=slicer, metric=metric + ) + + self.assertIsInstance(result, pd.Series) + self.assertEqual(result.name, "Total Teff") + + # Check that result index is a subset of transition_mjds + # (some may fail due to no visits) + self.assertTrue(set(result.index).issubset(set(transition_mjds))) + + # All values should be numeric and non-negative + for value in result.values: + self.assertIsInstance(value, (int, float)) + self.assertGreaterEqual(value, 0.0) + + @unittest.skipUnless(HAVE_MAF, "No maf installation") + def test_make_metric_progress_df(self): + # Create some sample data for completed and baseline visits + # Using a subset of the baseline visits for both + current_mjd = self.start_mjd + 5.0 + current_dayobs = DayObs.from_time(current_mjd) + completed_visits = self.visits[self.visits.observationStartMJD < current_mjd] + extrapolation_dayobs = DayObs.from_time(current_mjd + 8) + + # Test with default frequency + result = make_metric_progress_df( + completed_visits=completed_visits, + baseline_visits=self.visits, + start_dayobs=self.start_dayobs, + last_completed_dayobs=current_dayobs, + extrapolation_dayobs=extrapolation_dayobs, + slicer_factory=maf.UniSlicer, + metric_factory=partial(maf.CountMetric, col="observationStartMJD", metric_name="count"), + freq="D", + ) + assert isinstance(result, pd.DataFrame) + assert result.index.name == "jd" + assert len(result) > 0 + assert "date" in result.columns + + assert np.all(result.loc[:, "baseline"].values >= 0) + assert np.all(result.loc[: current_dayobs.end.jd, "snapshot"] >= 0) + assert np.all(np.isnan(result.loc[current_dayobs.end.jd :, "snapshot"])) + assert np.all(result.loc[: current_dayobs.end.jd, "chimera"] >= 0) + assert np.all(np.isnan(result.loc[current_dayobs.end.jd :, "chimera"])) + + result = make_metric_progress_df( + completed_visits=completed_visits, + baseline_visits=self.visits, + start_dayobs=self.start_dayobs, + last_completed_dayobs=current_dayobs, + extrapolation_dayobs=extrapolation_dayobs, + slicer_factory=partial(maf.HealpixSlicer, nside=8, verbose=False), + metric_factory=partial(maf.CountExplimMetric, metric_name="fO"), + summary_metric_factory=partial( + maf.FOArea, + nside=8, + norm=False, + metric_name="fOArea", + asky=18000, + n_visit=5, + badval=0, + ), + freq="D", + ) + assert isinstance(result, pd.DataFrame) + assert result.index.name == "jd" + assert len(result) > 0 + assert "date" in result.columns + assert np.all(result.loc[:, "baseline"].values >= 0) + assert np.all(result.loc[: current_dayobs.end.jd, "snapshot"] >= 0) + assert np.all(np.isnan(result.loc[current_dayobs.end.jd :, "snapshot"])) + assert np.all(result.loc[: current_dayobs.end.jd, "chimera"] >= 0) + assert np.all(np.isnan(result.loc[current_dayobs.end.jd :, "chimera"])) diff --git a/tests/test_compute_visits.py b/tests/test_compute_visits.py index dd729641..6c985256 100644 --- a/tests/test_compute_visits.py +++ b/tests/test_compute_visits.py @@ -1,6 +1,7 @@ import unittest import numpy as np +import pandas as pd from astropy.time import Time from rubin_scheduler.utils import SURVEY_START_MJD from rubin_sim.data import get_baseline @@ -61,3 +62,38 @@ def test_accum_stats_by_target_band_night(self): self.assertEqual(night_teff.index.names[1], "day_obs_iso8601") for col_name in night_teff.columns: self.assertTrue(col_name in "ugrizy") + + def test_match_visits_to_pointings_basic(self): + """Test basic matching functionality""" + # Create sample visit data + # Use RA=0 everywhere so separation is just differenc in declination + visits_data = {"s_ra": [0.0, 0.0, 0.0], "s_dec": [0.0, 1.0, 2.0], "filter": ["g", "r", "i"]} + visits = pd.DataFrame(visits_data) + + # Create pointings dictionary + pointings = { + "pointing1": (0.0, 0.0), + "pointing2": (0.0, 1.0), + "pointing3": (0.0, 3.0), + "pointing4": (0.0, 90.0), + } + + # Test matching + result = schedview.compute.visits.match_visits_to_pointings(visits, pointings, match_radius=1.5) + + # Check that we got results + self.assertIsInstance(result, pd.DataFrame) + self.assertIn("pointing_name", result.columns) + + # Check that we have the right number of matching visits + # pointing1 should match visits 0 and 1 + self.assertEqual(len(result[result["pointing_name"] == "pointing1"]), 2) + + # pointing2 should match visits 0, 1, and 2 + self.assertEqual(len(result[result["pointing_name"] == "pointing2"]), 3) + + # pointing3 should match only visit 2 + self.assertEqual(len(result[result["pointing_name"] == "pointing3"]), 1) + + # pointing4 should have no matches + self.assertEqual(len(result[result["pointing_name"] == "pointing4"]), 0) diff --git a/tests/test_plot_monthlyhg.py b/tests/test_plot_monthlyhg.py new file mode 100644 index 00000000..84026a48 --- /dev/null +++ b/tests/test_plot_monthlyhg.py @@ -0,0 +1,56 @@ +import datetime +import unittest +from unittest.mock import MagicMock +from xml.etree import ElementTree as ET + +from matplotlib import pyplot as plt + +from schedview.plot.monthlyhg import plot_collapsed_monthly_hourglass_metric + + +class TestPlotCollapsedMonthlyHourglassMetric(unittest.TestCase): + + def setUp(self): + """Create a mock MetricBundle to support testing without maf.""" + self.bundle = MagicMock() + self.bundle.plot_dict = {} + self.bundle.plot_funcs = [] + self.bundle.plot.return_value = {"dummy": plt.figure()} + + def test_html_output_for_two_months(self): + """Run the function over a two‑month span and sanity check result.""" + + first_date = datetime.date.fromisoformat("2027-01-01") + last_date = datetime.date.fromisoformat("2027-02-01") + + result = plot_collapsed_monthly_hourglass_metric( + metric_bundle=self.bundle, + name="TestMetric", + first_date=first_date, + last_date=last_date, + ) + + assert isinstance(result, str) + assert len(result) > 0 + + # Parse the HTML using ElementTree + root = ET.fromstring(f"{result}") + + # Check we have all the children we expect + details_elements = root.findall(".//details") + assert len(details_elements) == 2 + + summary_elements = root.findall(".//summary") + assert len(summary_elements) == 2 + + img_elements = root.findall(".//img") + assert len(img_elements) == 2 + + # Verify that each img has a src attribute with base64 data + for img in img_elements: + src = img.get("src", "") + assert src.startswith("data:image/png;base64,") + + +if __name__ == "__main__": + unittest.main()