OPSIM-1176 : Add additional tools for comcam sv surveys (#91)

docs/fbs-output-schema.rst

@@ -62,7 +62,7 @@ All values are for the center of the field of view (e.g., airmass, altitude, etc
 +-----------------------+-------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | paraAngle             | degrees           | Paralactic angle of the observation                                                                                                                                                               |
 +-----------------------+-------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
-| psudoParaAngle        | degrees           | Psudo-paralactic angle of the observation, see https://smtn-019.lsst.io/v/DM-44258/index.html                                                                                                                                                              |
+| psudoParaAngle        | degrees           | pseudo-paralactic angle of the observation, see https://smtn-019.lsst.io/v/DM-44258/index.html                                                                                                                                                              |
 +-----------------------+-------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | proposalId            | int               | deprecated                                                                                                                                                                                        |
 +-----------------------+-------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

rubin_scheduler/scheduler/basis_functions/basis_functions.py

@@ -44,6 +44,7 @@
     "BalanceVisits",
     "RewardNObsSequence",
     "FilterDistBasisFunction",
+    "RewardRisingBasisFunction",
 )
 
 import warnings
@@ -943,18 +944,9 @@ def __init__(self, filtername="r", nside=None, gap_min=25.0, penalty_val=np.nan)
 
     def _calc_value(self, conditions, indx=None):
         result = np.ones(hp.nside2npix(self.nside), dtype=float)
-        if indx is None:
-            indx = np.arange(result.size)
-        diff = IntRounded(conditions.mjd - self.survey_features["Last_observed"].feature[indx])
+        diff = IntRounded(conditions.mjd - self.survey_features["Last_observed"].feature)
         bad = np.where(diff < self.gap_min)[0]
-        # If this is used with a FieldSurvey or if indx is single value:
-        if isinstance(indx, np.int64):
-            if diff < self.gap_min:
-                result = self.penalty_val
-            else:
-                result = 0
-        else:
-            result[indx[bad]] = self.penalty_val
+        result[bad] = self.penalty_val
         return result
 
 
@@ -2261,3 +2253,44 @@ def __init__(self, n_obs_survey, note_survey, nside=None):
 
     def _calc_value(self, conditions, indx=None):
         return self.survey_features["n_obs_survey"].feature % self.n_obs_survey
+
+
+class RewardRisingBasisFunction(BaseBasisFunction):
+    """Reward parts of the sky that are rising.
+    Optionally, mask out parts of the sky that are not rising.
+
+    This produces a reward that increases
+    as the field rises toward zenith, then abruptly
+    falls as the field passes zenith.
+    Negative hour angles (or hour angles > 180 degrees)
+    indicate a rising point on the sky.
+
+    Parameters
+    ----------
+    slope : `float`
+        Sets the 'slope' of how fast the basis function
+        value changes with hour angle.
+    penalty_val : `float` or `np.nan`, optional
+        Sets the value for the part of the sky which is
+        not rising (hour angle between 0 and 180).
+        Using a value of np.nan will mask this region of sky,
+        a value of 0 will just make this non-rewarding.
+    nside : `int` or None, optional
+        Nside for the healpix map, default of None uses scheduler default.
+    """
+
+    def __init__(self, slope=0.1, penalty_val=0, nside=None):
+        super().__init__(nside=nside)
+        self.slope = slope
+        self.penalty_val = penalty_val
+
+    # Probably not needed
+    def check_feasibility(self, conditions):
+        return True
+
+    def _calc_value(self, conditions, indx=None):
+        # HA should be available in the conditions object
+        value = self.slope * conditions.HA
+        past_zenith = np.where(conditions.HA < np.pi)
+        value[past_zenith] = self.penalty_val
+        return value

rubin_scheduler/scheduler/basis_functions/feasibility_funcs.py

@@ -32,6 +32,16 @@
 from rubin_scheduler.scheduler.utils import IntRounded
 
 
+def send_unused_deprecation_warning(name):
+    message = (
+        f"The feasibility basis function {name} is not in use, "
+        "may be broken, and will be deprecated shortly. "
+        "Please contact the rubin_scheduler maintainers if "
+        "this is in use elsewhere."
+    )
+    warnings.warn(message, FutureWarning)
+
+
 class FilterLoadedBasisFunction(BaseBasisFunction):
     """Check that the filter(s) needed are loaded.
 
@@ -164,7 +174,8 @@ class LimitObsPnightBasisFunction(BaseBasisFunction):
     def __init__(self, survey_str="", nlimit=100.0):
         super(LimitObsPnightBasisFunction, self).__init__()
         self.nlimit = nlimit
-        self.survey_features["N_in_night"] = features.Survey_in_night(survey_str=survey_str)
+        self.survey_features["N_in_night"] = features.SurveyInNight(survey_str=survey_str)
+        send_unused_deprecation_warning(self.__class__.__name__)
 
     def check_feasibility(self, conditions):
         if self.survey_features["N_in_night"].feature >= self.nlimit:
@@ -309,19 +320,27 @@ def check_feasibility(self, conditions):
 
 
 class ForceDelayBasisFunction(BaseBasisFunction):
-    """Keep a survey from executing to rapidly.
+    """Keep a survey from executing too rapidly.
 
     Parameters
     ----------
-    days_delay : float (2)
+    days_delay : `float`, optional
         The number of days to force a gap on.
+    scheduler_note : `str` or None, optional
+        The value of the scheduler_note to count.
+        Default None will not consider scheduler_note.
+    survey_name : `str` or None, optional
+        Backwards compatible version of scheduler_note.
     """
 
-    def __init__(self, days_delay=2.0, survey_name=None):
+    def __init__(self, days_delay=2.0, scheduler_note=None, survey_name=None):
         super(ForceDelayBasisFunction, self).__init__()
         self.days_delay = days_delay
-        self.survey_name = survey_name
-        self.survey_features["last_obs_self"] = features.LastObservation(survey_name=self.survey_name)
+        if scheduler_note is None and survey_name is not None:
+            self.scheduler_note = survey_name
+        else:
+            self.scheduler_note = scheduler_note
+        self.survey_features["last_obs_self"] = features.LastObservation(scheduler_note=self.scheduler_note)
 
     def check_feasibility(self, conditions):
         result = True
@@ -336,23 +355,35 @@ class SoftDelayBasisFunction(BaseBasisFunction):
 
     Parameters
     ----------
-
+    fractions : `list` [`float`]
+    delays : `list` [`float`]
+    scheduler_note : `str` or None, optional
+        The scheduler_note to identify observations from a given survey or
+        survey mode.
+    survey_name : `str` or None, optional
+        Deprecated version of scheduler_note. Overriden by scheduler_note
+        if not None.
     """
 
-    def __init__(self, fractions=[0.000, 0.009, 0.017], delays=[0.0, 0.5, 1.5], survey_name=None):
+    def __init__(
+        self, fractions=[0.000, 0.009, 0.017], delays=[0.0, 0.5, 1.5], scheduler_note=None, survey_name=None
+    ):
         if len(fractions) != len(delays):
             raise ValueError("fractions and delays must be same length")
         super(SoftDelayBasisFunction, self).__init__()
         self.delays = delays
-        self.survey_name = survey_name
-        self.survey_features["last_obs_self"] = features.LastObservation(survey_name=self.survey_name)
+        if scheduler_note is None and survey_name is not None:
+            self.scheduler_note = survey_name
+        else:
+            self.scheduler_note = scheduler_note
+        self.survey_features["last_obs_self"] = features.LastObservation(scheduler_note=self.scheduler_note)
         self.fractions = fractions
-        self.survey_features["Ntot"] = features.NObsSurvey()
-        self.survey_features["N_survey"] = features.NObsSurvey(note=self.survey_name)
+        self.survey_features["N_total"] = features.NObsSurvey(note=None)
+        self.survey_features["N_note"] = features.NObsSurvey(note=self.scheduler_note)
 
     def check_feasibility(self, conditions):
         result = True
-        current_ratio = self.survey_features["N_survey"].feature / self.survey_features["Ntot"].feature
+        current_ratio = self.survey_features["N_note"].feature / self.survey_features["N_total"].feature
         indx = np.searchsorted(self.fractions, current_ratio)
         if indx == len(self.fractions):
             indx -= 1
@@ -402,33 +433,41 @@ def check_feasibility(self, conditions):
 
 
 class FractionOfObsBasisFunction(BaseBasisFunction):
-    """Limit the fraction of all observations that can be labled a certain
-    survey name. Useful for keeping DDFs from exceeding a given fraction
+    """Limit the fraction of all observations that can be labelled a certain
+    scheduler note.
+
+    Useful for keeping DDFs from exceeding a given fraction
     of the total survey.
 
     Parameters
     ----------
-    frac_total : float
+    frac_total : `float`
         The fraction of total observations that can be of this survey
-    survey_name : str
-        The name of the survey
+    scheduler_note : `str` or None, optional
+        The scheduler_note to identify observations from a given survey or
+        survey mode.
+    survey_name : `str` or None, optional
+        Deprecated version of scheduler_note. Overriden by scheduler_note
+        if scheduler_note not None.
     """
 
-    def __init__(self, frac_total, survey_name=""):
+    def __init__(self, frac_total, scheduler_note=None, survey_name=None):
         super(FractionOfObsBasisFunction, self).__init__()
-        self.survey_name = survey_name
+        if scheduler_note is None and survey_name is not None:
+            self.scheduler_note = survey_name
+        else:
+            self.scheduler_note = scheduler_note
         self.frac_total = frac_total
-        self.survey_features["Ntot"] = features.NObsSurvey()
-        self.survey_features["N_survey"] = features.NObsSurvey(note=self.survey_name)
+        self.survey_features["N_total"] = features.NObsSurvey(note=None)
+        self.survey_features["N_note"] = features.NObsSurvey(note=self.scheduler_note)
 
     def check_feasibility(self, conditions):
         # If nothing has been observed, fine to go
         result = True
-        if self.survey_features["Ntot"].feature == 0:
-            return result
-        ratio = self.survey_features["N_survey"].feature / self.survey_features["Ntot"].feature
-        if ratio > self.frac_total:
-            result = False
+        if self.survey_features["N_total"].feature > 0:
+            ratio = self.survey_features["N_note"].feature / self.survey_features["N_total"].feature
+            if ratio > self.frac_total:
+                result = False
         return result
 
 
@@ -451,8 +490,12 @@ class LookAheadDdfBasisFunction(BaseBasisFunction):
     ha_limits : list of lists (None)
         limits for what hour angles are acceptable (hours). e.g.,
         to give 4 hour window around HA=0, ha_limits=[[22,24], [0,2]]
-    survey_name : str ('')
-        The name of the survey
+    scheduler_note : `str` or None, optional
+        The scheduler_note to identify observations from a given survey or
+        survey mode.
+    survey_name : `str` or None, optional
+        Deprecated version of scheduler_note. Overriden by scheduler_note
+        if scheduler_note not None.
     time_jump : float (44.)
         The amount of time to assume will jump ahead if another survey
         executes (minutes)
@@ -467,28 +510,32 @@ def __init__(
         time_needed=30.0,
         RA=0.0,
         ha_limits=None,
-        survey_name="",
+        scheduler_note=None,
+        survey_name=None,
         time_jump=44.0,
         sun_alt_limit=-18.0,
     ):
         super(LookAheadDdfBasisFunction, self).__init__()
         if aggressive_fraction > frac_total:
             raise ValueError("aggressive_fraction should be less than frac_total")
-        self.survey_name = survey_name
+        if scheduler_note is None and survey_name is not None:
+            self.scheduler_note = survey_name
+        else:
+            self.scheduler_note = scheduler_note
         self.frac_total = frac_total
         self.ra_hours = RA / 360.0 * 24.0
         self.ha_limits = np.array(ha_limits)
         self.sun_alt_limit = str(int(sun_alt_limit)).replace("-", "n")
         self.time_jump = time_jump / 60.0 / 24.0  # To days
         self.time_needed = time_needed / 60.0 / 24.0  # To days
         self.aggressive_fraction = aggressive_fraction
-        self.survey_features["Ntot"] = features.NObsSurvey()
-        self.survey_features["N_survey"] = features.NObsSurvey(note=self.survey_name)
+        self.survey_features["N_total"] = features.NObsSurvey(note=None)
+        self.survey_features["N_note"] = features.NObsSurvey(note=self.scheduler_note)
 
     def check_feasibility(self, conditions):
         result = True
         target_ha = (conditions.lmst - self.ra_hours) % 24
-        ratio = self.survey_features["N_survey"].feature / self.survey_features["Ntot"].feature
+        ratio = self.survey_features["N_note"].feature / self.survey_features["N_total"].feature
         available_time = getattr(conditions, "sun_" + self.sun_alt_limit + "_rising") - conditions.mjd
         # If it's more that self.time_jump to hour angle zero
         # See if there will be enough time to twilight in the future

rubin_scheduler/scheduler/detailers/dither_detailer.py

@@ -1,4 +1,9 @@
-__all__ = ("DitherDetailer", "CameraRotDetailer", "EuclidDitherDetailer")
+__all__ = (
+    "DitherDetailer",
+    "EuclidDitherDetailer",
+    "CameraRotDetailer",
+    "CameraSmallRotPerObservationListDetailer",
+)
 
 import numpy as np
 
@@ -290,3 +295,89 @@ def __call__(self, observation_list, conditions):
             obs["rotTelPos"] = offsets[i]
 
         return observation_list
+
+
+class CameraSmallRotPerObservationListDetailer(BaseDetailer):
+    """
+    Randomly set the camera rotation for each observation list.
+
+    Generates a small sequential offset for sequential visits
+    in the same filter; adds a random offset for each filter change.
+
+    Parameters
+    ----------
+    max_rot : `float`, optional
+        The maximum amount to offset the camera (degrees).
+        Default of 85 allows some padding for camera rotator.
+    min_rot : `float`, optional
+        The minimum to offset the camera (degrees)
+        Default of -85 allows some padding for camera rotator.
+    seed : `int`, optional
+        Seed for random number generation (per night).
+    per_visit_rot : `float`, optional
+        Sequential rotation to add per visit.
+    telescope : `str`, optional
+        Telescope name. Options of "rubin" or "auxtel". Default "rubin".
+        This is used to determine conversions between rotSkyPos and rotTelPos.
+    """
+
+    def __init__(self, max_rot=85.0, min_rot=-85.0, seed=42, per_visit_rot=0.0, telescope="rubin"):
+        self.survey_features = {}
+
+        self.current_night = -1
+        self.max_rot = np.radians(max_rot)
+        self.min_rot = np.radians(min_rot)
+        self.rot_range = self.max_rot - self.min_rot
+        self.seed = seed
+        self.per_visit_rot = np.radians(per_visit_rot)
+        self.offset = None
+        self.rc = rotation_converter(telescope=telescope)
+
+    def _generate_offsets_filter_change(self, filter_list, mjd, initial_offset):
+        """Generate a random camera rotation for each filter change
+        or add a small offset for each sequential observation.
+        """
+        mjd_hash = round(100 * (np.asarray(mjd).item() % 100))
+        rng = np.random.default_rng(mjd_hash * self.seed)
+
+        offsets = np.zeros(len(filter_list))
+
+        # Find the locations of the filter changes
+        filter_changes = np.where(np.array(filter_list[:-1]) != np.array(filter_list[1:]))[0]
+        filter_changes = np.concatenate([np.array([-1]), filter_changes])
+        # But add one because of counting and offsets above.
+        filter_changes += 1
+        # Count visits per filter in the sequence.
+        nvis_per_filter = np.concatenate(
+            [np.diff(filter_changes), np.array([len(filter_list) - 1 - filter_changes[-1]])]
+        )
+        # Set up the random rotator offsets for each filter change
+        # This includes first rotation .. maybe not needed?
+        for fchange_idx, nvis_f in zip(filter_changes, nvis_per_filter):
+            rot_range = self.rot_range - self.per_visit_rot * nvis_f
+            # At the filter change spot, update to random offset
+            offsets[fchange_idx] = rng.random() * rot_range + self.min_rot
+            # After the filter change point, add incremental rotation
+            # (we'll wipe this when we get to next fchange_idx)
+            offsets[fchange_idx:] += self.per_visit_rot * np.arange(len(filter_list) - fchange_idx)
+
+        return offsets
+
+    def __call__(self, observation_list, conditions):
+        # Generate offsets in camera rotator
+        filter_list = [np.asarray(obs["filter"]).item() for obs in observation_list]
+        offsets = self._generate_offsets_filter_change(filter_list, conditions.mjd, conditions.rot_tel_pos)
+
+        for i, obs in enumerate(observation_list):
+            alt, az = _approx_ra_dec2_alt_az(
+                obs["RA"],
+                obs["dec"],
+                conditions.site.latitude_rad,
+                conditions.site.longitude_rad,
+                conditions.mjd,
+            )
+            obs_pa = _approx_altaz2pa(alt, az, conditions.site.latitude_rad)
+            obs["rotSkyPos"] = self.rc._rottelpos2rotskypos(offsets[i], obs_pa)
+            obs["rotTelPos"] = offsets[i]
+
+        return observation_list

rubin_scheduler/scheduler/example/__init__.py

@@ -1 +1,2 @@
+from .comcam_sv_scheduler_example import *
 from .example_scheduler import *