diff --git a/docs/fbs-output-schema.rst b/docs/fbs-output-schema.rst index 23548547..dd008eec 100644 --- a/docs/fbs-output-schema.rst +++ b/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 | +-----------------------+-------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+ diff --git a/rubin_scheduler/scheduler/basis_functions/basis_functions.py b/rubin_scheduler/scheduler/basis_functions/basis_functions.py index 37c32208..ee79cdaf 100644 --- a/rubin_scheduler/scheduler/basis_functions/basis_functions.py +++ b/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 diff --git a/rubin_scheduler/scheduler/basis_functions/feasibility_funcs.py b/rubin_scheduler/scheduler/basis_functions/feasibility_funcs.py index b7474aaa..dda6044f 100644 --- a/rubin_scheduler/scheduler/basis_functions/feasibility_funcs.py +++ b/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,14 +510,18 @@ 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) @@ -482,13 +529,13 @@ def __init__( 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 diff --git a/rubin_scheduler/scheduler/detailers/dither_detailer.py b/rubin_scheduler/scheduler/detailers/dither_detailer.py index b07cdaae..04ff3542 100644 --- a/rubin_scheduler/scheduler/detailers/dither_detailer.py +++ b/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 diff --git a/rubin_scheduler/scheduler/example/__init__.py b/rubin_scheduler/scheduler/example/__init__.py index 133d7100..eb803c6e 100644 --- a/rubin_scheduler/scheduler/example/__init__.py +++ b/rubin_scheduler/scheduler/example/__init__.py @@ -1 +1,2 @@ +from .comcam_sv_scheduler_example import * from .example_scheduler import * diff --git a/rubin_scheduler/scheduler/example/comcam_sv_scheduler_example.py b/rubin_scheduler/scheduler/example/comcam_sv_scheduler_example.py new file mode 100644 index 00000000..ec0053ce --- /dev/null +++ b/rubin_scheduler/scheduler/example/comcam_sv_scheduler_example.py @@ -0,0 +1,950 @@ +__all__ = ( + "get_model_observatory", + "update_model_observatory_sunset", + "standard_masks", + "simple_rewards", + "simple_pairs_survey", + "simple_greedy_survey", + "get_basis_functions_field_survey", + "get_field_survey", + "get_sv_fields", + "prioritize_fields", + "get_comcam_sv_schedulers", +) + +import copy + +import numpy as np +from astropy.time import Time + +import rubin_scheduler.scheduler.basis_functions as basis_functions +import rubin_scheduler.scheduler.detailers as detailers +from rubin_scheduler.scheduler.detailers import CameraSmallRotPerObservationListDetailer +from rubin_scheduler.scheduler.model_observatory import ( + KinemModel, + ModelObservatory, + rotator_movement, + tma_movement, +) +from rubin_scheduler.scheduler.schedulers import ComCamFilterSched, CoreScheduler, FilterSwapScheduler +from rubin_scheduler.scheduler.surveys import BlobSurvey, FieldSurvey, GreedySurvey +from rubin_scheduler.scheduler.utils import Footprint, get_current_footprint +from rubin_scheduler.site_models import Almanac, ConstantSeeingData, ConstantWindData +from rubin_scheduler.utils import ddf_locations, survey_start_mjd + + +def get_model_observatory( + dayobs: str = "2024-09-09", + fwhm_500: float = 2.0, + wind_speed: float = 5.0, + wind_direction: float = 340, + tma_percent: float = 10, + rotator_percent: float = 100, + survey_start: Time = Time("2024-09-09T12:00:00", format="isot", scale="utc").mjd, +) -> ModelObservatory: + """Set up a model observatory with constant seeing and wind speed, + running on 'dayobs'. + + Track a constant 'survey_start' over different 'dayobs' so that the + 'night' value in the simulation outputs is consistent, and any + basis functions which might use 'season' values from the conditions + have appropriate values. + + Parameters + ---------- + dayobs : `str` + DAYOBS formatted str (YYYY-MM-DD) for the evening to start + up the observatory. + fwhm_500 : `float`, optional + The value to set for atmospheric component of seeing, + constant seeing throughout the night (arcseconds). + Ad-hoc default for start of comcam on-sky operations about 2.0". + wind_speed : `float`, optional + Set a (constant) wind speed for the night, (m/s). + Default of 5.0 is minimal but noticeable. + wind_direction : `float`, optional + Set a (constant) wind direction for the night (deg). + Default of 340 is midrange between typical wind directions for + the site (see SITCOMTN-126). + tma_percent : `float`, optional + Set a percent of full-performance for the telescope TMA (0-100). + Default of 10(%) is likely for start of comcam on-sky SV surveys. + rotator_percent : `float`, optional + Set a percent of full-performance for the rotator. + Default of 100% is likely for the start of comcam on-sky SV surveys. + survey_start : `float`, optional + MJD of the day of the survey start of operations. + This should be kept constant for a given set of simulations, + so that the "night" value in the output is consistent. + For surveys which use basis functions depending on season, this + is also important to be constant. + + Returns + ------- + observatory : `~.scheduler.model_observatory.ModelObservatory` + A ModelObservatory set up to start operations in the evening + of DAYOBS. + + Notes + ----- + The time for the model observatory will be advanced to the time + of `sunset_start_key` (default -12 degree sunset) in the model + observatory. The filters may not be correct however; use + `update_model_observatory_sunset` to get filters in place. + """ + # Set up a fresh model observatory + mjd_now = Time(f"{dayobs}T12:00:00", format="isot", scale="utc").mjd + + kinematic_model = KinemModel(mjd0=mjd_now) + rot = rotator_movement(rotator_percent) + kinematic_model.setup_camera(**rot) + tma = tma_movement(tma_percent) + kinematic_model.setup_telescope(**tma) + + # Some weather telemetry that might be useful + seeing_data = ConstantSeeingData(fwhm_500=fwhm_500) + wind_data = ConstantWindData(wind_direction=wind_direction, wind_speed=wind_speed) + + # Set up the model observatory + observatory = ModelObservatory( + mjd=mjd_now, + mjd_start=survey_start, + kinem_model=kinematic_model, # Modified kinematics + cloud_data="ideal", # No clouds + seeing_data=seeing_data, # Modified seeing + wind_data=wind_data, # Add some wind + downtimes="ideal", # No downtime + lax_dome=True, # dome crawl? + init_load_length=1, # size of skybrightness files to load first + ) + return observatory + + +def update_model_observatory_sunset( + observatory: ModelObservatory, filter_scheduler: FilterSwapScheduler, twilight: int | float = -12 +) -> ModelObservatory: + """Move model observatory to twilight and ensure correct filters are in + place according to the filter_scheduler. + + Parameters + ---------- + observatory : `~.scheduler.model_observatory.ModelObservatory` + The ModelObservatory simulating the observatory. + filter_scheduler : `~.scheduler.schedulers.FilterScheduler` + The filter scheduler providing appropriate information on + the filters that should be in place on the current observatory day. + twilight : `int` or `float` + If twilight is -12 or -18, the Almanac -12 or -18 degree twilight + times are used to set the current observatory time. + If any other value is provided, it is assumed to be a specific + MJD to start operating the observatory. + Filter choices are based on the time after advancing to twilight. + + Returns + ------- + observatory : `~.scheduler.model_observatory.ModelObservatory` + The ModelObservatory simulating the observatory, updated to the + time of 'twilight' and with mounted_filters matching + the filters chosen by the filter_scheduler for the current time + at 'twilight'. + """ + # Move to *next* possible sunset + if twilight in (-12, -18): + time_key = f"sun_n{twilight*-1 :0d}_setting" + # Use the observatory almanac to find the exact time to forward to + alm_indx = np.searchsorted(observatory.almanac.sunsets[time_key], observatory.mjd, side="right") - 1 + new_mjd = observatory.almanac.sunsets[time_key][alm_indx] + if new_mjd < observatory.mjd: + alm_indx += 1 + new_mjd = observatory.almanac.sunsets[time_key][alm_indx] + observatory.mjd = new_mjd + else: + # Assume twilight was a particular (MJD) time + observatory.mjd = twilight + + # Make sure correct filters are mounted + conditions = observatory.return_conditions() + filters_needed = filter_scheduler(conditions) + observatory.observatory.mount_filters(filters_needed) + return observatory + + +def standard_masks( + nside: int, + moon_distance: float = 30.0, + wind_speed_maximum: float = 20.0, + min_alt: float = 20, + max_alt: float = 86.5, + min_az: float = 0, + max_az: float = 360, + shadow_minutes: float = 30, +) -> list[basis_functions.BaseBasisFunction]: + """A set of standard mask functions. + + Avoids the moon, bright planets, high wind, and + areas on the sky out of bounds, using + the MoonAvoidanceBasisFunction, PlanetMaskBasisFunction, + AvoidDirectWindBasisFunction, and the AltAzShadowMaskBasisFunction. + + Parameters + ---------- + nside : `int` or None + The healpix nside to use. + Default of None uses rubin_scheduler.utils.get_default_nside. + moon_distance : `float`, optional + Moon avoidance distance, in degrees. + wind_speed_maximum : `float`, optional + Wind speed maximum to apply to the wind avoidance basis function, + in m/s. + min_alt : `float`, optional + Minimum altitude (in degrees) to observe. + max_alt : `float`, optional + Maximum altitude (in degrees) to observe. + min_az : `float`, optional + Minimum azimuth angle (in degrees) to observe. + max_az : `float`, optional + Maximum azimuth angle (in degrees) to observe. + shadow_minutes : `float`, optional + Avoid inaccessible alt/az regions, as well as parts of the sky + which will move into those regions within `shadow_minutes` (minutes). + + Returns + ------- + mask_basis_functions : `list` [`BaseBasisFunction`] + Mask basis functions should always be used with a weight of 0. + The masked (np.nan or -np.inf) regions will remain masked, + but the basis function values won't influence the reward. + """ + masks = [] + # Add the Moon avoidance mask + masks.append(basis_functions.MoonAvoidanceBasisFunction(nside=nside, moon_distance=moon_distance)) + # Add a mask around bright planets + masks.append(basis_functions.PlanetMaskBasisFunction(nside=nside)) + # Add the wind avoidance mask + masks.append(basis_functions.AvoidDirectWind(nside=nside, wind_speed_maximum=wind_speed_maximum)) + # Avoid inaccessible parts of the sky, as well as places that will + # move into those places within shadow_minutes. + masks.append( + basis_functions.AltAzShadowMaskBasisFunction( + nside=nside, + min_alt=min_alt, + max_alt=max_alt, + # min_az=min_az, + # max_az=max_az, + shadow_minutes=shadow_minutes, + ) + ) + return masks + + +def simple_rewards( + footprints: Footprint, + filtername: str, + nside: int = 32, + m5_weight: float = 6.0, + footprint_weight: float = 1.5, + slewtime_weight: float = 3.0, + stayfilter_weight: float = 3.0, + repeat_weight: float = -20, +) -> list[basis_functions.BaseBasisFunction]: + """A simple set of rewards for area-based surveys. + + Parameters + ---------- + footprints : `Footprint` + A Footprint class, which takes a target map and adds a + time-dependent weighting on top (such as seasonal weighting). + filtername : `str` + The filtername active for these rewards. + nside : `int`, optional + The nside for the rewards. + m5_weight : `float`, optional + The weight to give the M5Diff basis function. + footprint_weight : `float`, optional + The weight to give the footprint basis function. + slewtime_weight : `float`, optional + The weight to give the slewtime basis function. + stayfilter_weight : `float`, optional + The weight to give the FilterChange basis function. + repeat_weight : `float`, optional + The weight to give the VisitRepeat basis function. + This is negative by default, to avoid revisiting the same part + of the sky within `gap_max` (3 hours) in any filter (except in the + pairs survey itself). + + Returns + ------- + reward_functions : `list` [(`BaseBasisFunction`, `float`)] + List of tuples, each tuple is a reward function followed by + its respective weight. + + Notes + ----- + Increasing the m5_weight moves visits toward darker skies. + Increasing the footprint weight distributes visits more evenly. + Increasing the slewtime weight acquires visits with shorter slewtime. + The balance in the defaults here has worked reasonably for pair + surveys in simulations. + """ + reward_functions = [] + # Add M5 basis function (rewards dark sky) + reward_functions.append( + (basis_functions.M5DiffBasisFunction(filtername=filtername, nside=nside), m5_weight) + ) + # Add a footprint basis function + # (rewards sky with fewer visits) + reward_functions.append( + ( + basis_functions.FootprintBasisFunction( + filtername=filtername, + footprint=footprints, + out_of_bounds_val=np.nan, + nside=nside, + ), + footprint_weight, + ) + ) + # Add a reward function for small slewtimes. + reward_functions.append( + ( + basis_functions.SlewtimeBasisFunction(filtername=filtername, nside=nside), + slewtime_weight, + ) + ) + # Add a reward to stay in the same filter as much as possible. + reward_functions.append( + (basis_functions.FilterChangeBasisFunction(filtername=filtername), stayfilter_weight) + ) + # And this is technically a mask, to avoid asking for observations + # which are not possible. However, it depends on the filters + # requested, compared to the filters available in the camera. + reward_functions.append((basis_functions.FilterLoadedBasisFunction(filternames=filtername), 0)) + # And add a basis function to avoid repeating the same pointing + # (VisitRepeat can be used to either encourage or discourage repeats, + # depending on the repeat_weight value). + reward_functions.append( + ( + basis_functions.VisitRepeatBasisFunction( + gap_min=0, gap_max=3 * 60.0, filtername=None, nside=nside, npairs=20 + ), + repeat_weight, + ) + ) + return reward_functions + + +def simple_pairs_survey( + nside: int = 32, + filtername: str = "g", + filtername2: str | None = None, + mask_basis_functions: list[basis_functions.BaseBasisFunction] | None = None, + reward_basis_functions: list[basis_functions.BaseBasisFunction] | None = None, + reward_basis_functions_weights: list[float] | None = None, + survey_start: float | None = None, + footprints_hp: np.ndarray | None = None, + camera_rot_limits: list[float] = [-80.0, 80.0], + pair_time: float = 30.0, + exptime: float = 30.0, + nexp: int = 1, +) -> BlobSurvey: + """Set up a simple blob survey to acquire pairs of visits. + + Parameters + ---------- + nside : `int`, optional + Nside for the surveys. + filtername : `str`, optional + Filtername for the first visit of the pair. + filtername2 : `str` or None, optional + Filtername for the second visit of the pair. If None, the + first filter will be used for both visits. + mask_basis_functions : `list` [`BaseBasisFunction`] or None + List of basis functions to use as masks (with implied weight 0). + If None, `standard_masks` is used with default parameters. + reward_basis_functions : `list` [`BaseBasisFunction`] or None + List of basis functions to use as rewards. + If None, a basic set of basis functions will be used. + reward_basis_functions_weights : `list` [`float`] or None + List of values to use as weights for the reward basis functions. + If None, default values for the basic set will be used. + survey_start : `float` or None + The start of the survey, in MJD. + If None, `survey_start_mjd()` is used. + This should be the start of the survey, not the current time. + footprints_hp : `np.ndarray` (N,) or None + An array of healpix maps with the target survey area, with dtype + like [(filtername, ' 0] + [ + (i[0], i[1]) for i in rf2 if i[1] > 0 + ] + # Then put back in the FilterLoadedBasisFunction with both filters. + filternames = [fn for fn in [filtername, filtername2] if fn is not None] + reward_functions.append((basis_functions.FilterLoadedBasisFunction(filternames=filternames), 0)) + + # unpack the basis functions and weights + reward_basis_functions_weights = [val[1] for val in reward_functions] + reward_basis_functions = [val[0] for val in reward_functions] + + # Set up blob surveys. + if filtername2 is None: + scheduler_note = "pair_%i, %s" % (pair_time, filtername) + else: + scheduler_note = "pair_%i, %s%s" % (pair_time, filtername, filtername2) + + # Set up detailers for each requested observation. + detailer_list = [] + # Avoid camera rotator limits. + detailer_list.append( + detailers.CameraRotDetailer(min_rot=np.min(camera_rot_limits), max_rot=np.max(camera_rot_limits)) + ) + # Convert rotTelPos to rotSkyPos_desired + detailer_list.append(detailers.Rottep2RotspDesiredDetailer(telescope="rubin")) + # Reorder visits in a blob so that closest to current altitude is first. + detailer_list.append(detailers.CloseAltDetailer()) + # Sets a flush-by date to avoid running into prescheduled visits. + detailer_list.append(detailers.FlushForSchedDetailer()) + # Add a detailer to label visits as either first or second of the pair. + if filtername2 is not None: + detailer_list.append(detailers.TakeAsPairsDetailer(filtername=filtername2)) + + # Set up the survey. + ignore_obs = ["DD"] + + BlobSurvey_params = { + "slew_approx": 7.5, + "filter_change_approx": 140.0, + "read_approx": 2.4, + "search_radius": 30.0, + "flush_time": pair_time * 3, + "smoothing_kernel": None, + "nside": nside, + "seed": 42, + "dither": True, + "twilight_scale": False, + } + + pair_survey = BlobSurvey( + reward_basis_functions + mask_basis_functions, + reward_basis_functions_weights + mask_basis_functions_weights, + filtername1=filtername, + filtername2=filtername2, + exptime=exptime, + ideal_pair_time=pair_time, + scheduler_note=scheduler_note, + ignore_obs=ignore_obs, + nexp=nexp, + detailers=detailer_list, + **BlobSurvey_params, + ) + + return pair_survey + + +def simple_greedy_survey( + nside: int = 32, + filtername: str = "r", + mask_basis_functions: list[basis_functions.BaseBasisFunction] | None = None, + reward_basis_functions: list[basis_functions.BaseBasisFunction] | None = None, + reward_basis_functions_weights: list[float] | None = None, + survey_start: float | None = None, + footprints_hp: np.ndarray | None = None, + camera_rot_limits: list[float] = [-80.0, 80.0], + exptime: float = 30.0, + nexp: int = 1, +) -> GreedySurvey: + """Set up a simple greedy survey to just observe single visits. + + Parameters + ---------- + nside : `int`, optional + Nside for the surveys. + filtername : `str`, optional + Filtername for the visits. + mask_basis_functions : `list` [`BaseBasisFunction`] or None + List of basis functions to use as masks (with implied weight 0). + If None, `standard_masks` is used with default parameters. + reward_basis_functions : `list` [`BaseBasisFunction`] or None + List of basis functions to use as rewards. + If None, a basic set of basis functions will be used. + reward_basis_functions_weights : `list` [`float`] or None + List of values to use as weights for the reward basis functions. + If None, default values for the basic set will be used. + survey_start : `float` or None + The start of the survey, in MJD. + If None, `survey_start_mjd()` is used. + This should be the start of the survey, not the current time. + footprints_hp : `np.ndarray` (N,) or None + An array of healpix maps with the target survey area, with dtype + like [(filtername, ' list[basis_functions.BaseBasisFunction]: + """Get the basis functions for a comcam SV field survey. + + Parameters + ---------- + nside : `int` + The nside value for the healpix grid. + wind_speed_maximum : `float` + Maximum wind speed tolerated for the observations of the survey, + in m/s. + moon_distance : `float`, optional + Moon avoidance distance, in degrees. + min_alt : `float`, optional + Minimum altitude (in degrees) to observe. + max_alt : `float`, optional + Maximum altitude (in degrees) to observe. + min_az : `float`, optional + Minimum azimuth angle (in degrees) to observe. + max_az : `float`, optional + Maximum azimuth angle (in degrees) to observe. + shadow_minutes : `float`, optional + Avoid inaccessible alt/az regions, as well as parts of the sky + which will move into those regions within `shadow_minutes` (minutes). + For the FieldSurvey, this should probably be the length of time + required by the sequence in the FieldSurvey, to avoid tracking into + inaccessible areas of sky. + + Returns + ------- + bfs : `list` of `~.scheduler.basis_functions.BaseBasisFunction` + """ + sun_alt_limit = -12.0 + moon_distance = 30 + + bfs = standard_masks( + nside=nside, + moon_distance=moon_distance, + wind_speed_maximum=wind_speed_maximum, + min_alt=min_alt, + max_alt=max_alt, + min_az=min_az, + max_az=max_az, + shadow_minutes=shadow_minutes, + ) + + bfs += [ + basis_functions.NotTwilightBasisFunction(sun_alt_limit=sun_alt_limit), + # Avoid revisits within 30 minutes + basis_functions.AvoidFastRevisitsBasisFunction(nside=nside, filtername=None, gap_min=30.0), + # reward fields which are rising, but don't mask out after zenith + basis_functions.RewardRisingBasisFunction(nside=nside, slope=0.1, penalty_val=0), + # Reward parts of the sky which are darker -- + # note that this is only for r band, so relying on skymap in r band. + # if there isn't a strong reason to go with the darkest pointing, + # it might be reasonable to just drop this basis function + basis_functions.M5DiffBasisFunction(filtername="r", nside=nside), + ] + return bfs + + +def get_field_survey( + field_ra_deg: float, + field_dec_deg: float, + field_name: str, + basis_functions: list[basis_functions.BaseBasisFunction], + detailers: list[detailers.BaseDetailer], + nside: int = 32, +) -> FieldSurvey: + """Set up a comcam SV field survey. + + Parameters + ---------- + field_ra_deg : `float` + The RA (in degrees) of the field. + field_dec_deg : `float` + The Dec (in degrees) of the field. + field_name : `str` + The name of the field. This is used for the survey_name and + transferred to the 'target' information in the output observation. + Also used in 'scheduler_note', which is important for the FieldSurvey + to know whether to count particular observations for the Survey. + basis_functions : `list` of [`~.scheduler.basis_function` objects] + Basis functions for the field survey. + A default set can be obtained from `get_basis_functions_field_survey`. + detailers : `list` of [`~.scheduler.detailer` objects] + Detailers for the survey. + Detailers can add information to output observations, including + specifying rotator or dither positions. + nside : `int`, optional + Nside for the survey. Default 32. + + Returns + ------- + field_survey : `~.scheduler.surveys.FieldSurvey` + The configured FieldSurvey. + + Notes + ----- + The sequences for a given field survey can be set via kwargs, + not necessarily easily accessible here. Only portions of the sequence + which correspond to mounted filters will be requested by the FieldSurvey. + + field_survey.extra_features['ObsRecord'] tracks how many observations + have been accepted by the Survey (and can be useful for diagnostics). + """ + field_survey = FieldSurvey( + basis_functions, + field_ra_deg, + field_dec_deg, + sequence="ugrizy", + nvisits={"u": 20, "g": 20, "r": 20, "i": 20, "z": 20, "y": 20}, + exptimes={"u": 38, "g": 30, "r": 30, "i": 30, "z": 30, "y": 30}, + nexps={"u": 1, "g": 2, "r": 2, "i": 2, "z": 2, "y": 2}, + ignore_obs=None, + accept_obs=[field_name], + survey_name=field_name, + scheduler_note=None, + readtime=2.4, + filter_change_time=120.0, + nside=nside, + flush_pad=30.0, + detailers=detailers, + ) + return field_survey + + +def get_sv_fields() -> dict[str, dict[str, float]]: + """Default potential fields for the SV surveys. + + Returns + ------- + fields_dict : `dict` {`str` : {'RA' : `float`, 'Dec' : `float`}} + A dictionary keyed by field_name, containing RA and Dec (in degrees) + for each field. + """ + fields = ( + ("Rubin_SV_095_-25", 95.0, -25.0), # High stellar densty, low extinction + ("Rubin_SV_125_-15", 125.0, -15.0), # High stellar densty, low extinction + ("DESI_SV3_R1", 179.60, 0.000), # DESI, GAMA, HSC DR2, KiDS-N + ("Rubin_SV_225_-40", 225.0, -40.0), # 225 High stellar densty, low extinction + ("DEEP_A0", 216, -12.5), # DEEP Solar Systen + ("Rubin_SV_250_2", 250.0, 2.0), # 250 High stellar densty, low extinction + ("Rubin_SV_300_-41", 300.0, -41.0), # High stellar densty, low extinction + ("Rubin_SV_280_-48", 280.0, -48.0), # High stellar densty, low extinction + ("DEEP_B0", 310, -19), # DEEP Solar System + ("ELAIS_S1", 9.45, -44.0), # ELAIS-S1 LSST DDF + ("XMM_LSS", 35.708333, -4.75), # LSST DDF + ("ECDFS", 53.125, -28.1), # ECDFS + ("COSMOS", 150.1, 2.1819444444444445), # COSMOS + ("EDFS_A", 58.9, -49.315), # EDFS_a + ("EDFS_B", 63.6, -47.6), # EDFS_b + ) + + fields_dict = dict(zip([f[0] for f in fields], [{"RA": f[1], "Dec": f[2]} for f in fields])) + + # Update ddf_locations to + survey_ddfs = ddf_locations() + ddf_radec = {} + for k in survey_ddfs: + kk = k.lower().replace("_", "").replace("-", "") + ddf_radec[kk] = survey_ddfs[k] + for field in fields_dict: + # see if the field names match any ddf field name, allowing for + # upper/lower case and -vs_ differences + ff = field.lower().replace("_", "").replace("-", "") + if ff in ddf_radec: + fields_dict[field]["RA"] = ddf_radec[ff][0] + fields_dict[field]["Dec"] = ddf_radec[ff][1] + + return fields_dict + + +def prioritize_fields( + priority_fields: list[str] | None = None, field_dict: dict[str, dict[str, float]] | None = None +) -> list[list[FieldSurvey]]: + """Add the remaining field names in field_dict into the last + tier of 'priority_fields' field names, creating a complete + survey tier list of lists. + + Parameters + ---------- + priority_fields : `list` [`list`] + A list of lists, where each final list corresponds to a 'tier' + of FieldSurveys, and contains those field survey names. + These names must be present in field_dict. + field_dict : `dict` {`str` : {'RA' : `float`, 'Dec' : `float`}} or None + Dictionary containing field information for the FieldSurveys. + Default None will fetch the SV fields from 'get_sv_fields'. + + Returns + ------- + tiers : `list` [`list`] + The tiers to pass to the core scheduler, after including the + non-prioritized fields from field_dict. + """ + if field_dict is None: + field_dict = get_sv_fields() + else: + field_dict = copy.deepcopy(field_dict) + tiers = [] + if priority_fields is not None: + for tier in priority_fields: + tiers.append(tier) + for field in tier: + del field_dict[field] + remaining_fields = list(field_dict.keys()) + tiers.append(remaining_fields) + return tiers + + +def get_comcam_sv_schedulers( + starting_tier: int = 0, + tiers: list[list[str]] | None = None, + field_dict: dict[str, dict[str, float]] | None = None, + nside: int = 32, +) -> (CoreScheduler, ComCamFilterSched): + """Set up a CoreScheduler and FilterScheduler generally + appropriate for ComCam SV observing. + + Parameters + ---------- + starting_tier : `int`, optional + Starting to tier to place the surveys coming from the 'tiers' + specified here. + Default 0, to start at first tier. If an additional + survey will be added at highest tier after (such as cwfs), then + set starting tier to 1+ and add these surveys as a list to + scheduler.survey_lists[tier] etc. + tiers : `list` [`str`] or None + Field names for each of the field surveys in tiers. + Should be a list of lists - [tier][surveys_in_tier] + [[field1, field2],[field3, field4, field5]. + Fields should be present in the 'field_dict'. + Default None will use all fields in field_dict. + field_dict : `dict` {`str` : {'RA' : `float`, 'Dec' : `float`}} or None + Dictionary containing field information for the FieldSurveys. + Default None will fetch the SV fields from 'get_sv_fields'. + nside : `int` + Nside for the scheduler. Default 32. + Generally, don't change without serious consideration. + + Returns + ------- + scheduler, filter_scheduler : `~.scheduler.schedulers.CoreScheduler`, + `~.scheduler.schedulers.ComCamFilterSched` + A CoreScheduler and FilterScheduler that are generally + appropriate for ComCam. + """ + if field_dict is None: + field_dict = get_sv_fields() + + if tiers is None: + tiers = [list(field_dict.keys())] + + surveys = [] + + i = 0 + for t in tiers: + if len(t) == 0: + continue + j = i + starting_tier + i += 1 + surveys.append([]) + for kk, fieldname in enumerate(t): + bfs = get_basis_functions_field_survey() + detailer = CameraSmallRotPerObservationListDetailer(per_visit_rot=0.5) + surveys[j].append( + get_field_survey( + field_dict[fieldname]["RA"], field_dict[fieldname]["Dec"], fieldname, bfs, [detailer] + ) + ) + + scheduler = CoreScheduler(surveys, nside=nside) + filter_scheduler = ComCamFilterSched() + return scheduler, filter_scheduler diff --git a/rubin_scheduler/scheduler/example/example_scheduler.py b/rubin_scheduler/scheduler/example/example_scheduler.py index 50f36be8..07f5a943 100644 --- a/rubin_scheduler/scheduler/example/example_scheduler.py +++ b/rubin_scheduler/scheduler/example/example_scheduler.py @@ -1,4 +1,15 @@ -__all__ = ("example_scheduler", "sched_argparser", "set_run_info", "run_sched") +__all__ = ( + "example_scheduler", + "sched_argparser", + "set_run_info", + "run_sched", + "gen_long_gaps_survey", + "gen_greedy_surveys", + "generate_blobs", + "generate_twi_blobs", + "generate_twilight_near_sun", + "standard_bf", +) import argparse import os @@ -33,7 +44,7 @@ iers.conf.auto_download = False -def example_scheduler(nside=32, mjd_start=survey_start_mjd()): +def example_scheduler(nside: int = 32, mjd_start: float = survey_start_mjd()) -> CoreScheduler: """Provide an example baseline survey-strategy scheduler. Parameters @@ -370,7 +381,6 @@ def blob_for_long( "slew_approx": 7.5, "filter_change_approx": 140.0, "read_approx": 2.0, - "min_pair_time": 15.0, "search_radius": 30.0, "alt_max": 85.0, "az_range": None, @@ -751,7 +761,6 @@ def generate_blobs( "slew_approx": 7.5, "filter_change_approx": 140.0, "read_approx": 2.0, - "min_pair_time": 15.0, "search_radius": 30.0, "alt_max": 85.0, "az_range": None, @@ -992,7 +1001,6 @@ def generate_twi_blobs( "slew_approx": 7.5, "filter_change_approx": 140.0, "read_approx": 2.0, - "min_pair_time": 10.0, "search_radius": 30.0, "alt_max": 85.0, "az_range": None, diff --git a/rubin_scheduler/scheduler/features/features.py b/rubin_scheduler/scheduler/features/features.py index 754fcb5c..efd8b5c8 100644 --- a/rubin_scheduler/scheduler/features/features.py +++ b/rubin_scheduler/scheduler/features/features.py @@ -121,8 +121,6 @@ class SurveyInNight(BaseSurveyFeature): ---------- survey_str : `str`, optional String to search for in observation `scheduler_note`. - String does not have to match `scheduler_note` exactly, - just be contained in `scheduler_note`. Default of "" means any observation will match. """ @@ -191,8 +189,8 @@ class NObsCount(BaseSurveyFeature): Parameters ---------- - filtername : `str` (None) - The filter to count (if None, all filters counted) + filtername : `str` or None + The filter to count. Default None, all filters counted. """ @@ -319,15 +317,21 @@ def add_observation(self, observation, indx=None): class NObsSurvey(BaseSurveyFeature): """Count the number of observations, whole sky (not per pixel). + Because this feature will belong to a survey, it would count all + observations that are counted for that survey. + Parameters ---------- - note : `str` (None) - Only count observations that contain str in their note field + note : `str` or None + Count observations that match `str` in their scheduler_note field. + Note can be a substring of scheduler_note, and will still match. """ def __init__(self, note=None): self.feature = 0 self.note = note + if self.note == "": + self.note = None def add_observations_array(self, observations_array, observations_hpid): if self.note is None: @@ -352,29 +356,34 @@ class LastObservation(BaseSurveyFeature): Parameters ---------- - survey_name : `str` (None) - Only records if the survey name matches (or survey_name set to None) + scheduler_note : `str` or None, optional + Value of the scheduler_note to match, if not None. + survey_name : `str` or None, optional + Backwards compatible version of scheduler_note. Deprecated. """ - def __init__(self, survey_name=None): - # Will this work for observations read from a database??? - # "note" is definitely NOT guaranteed to match the survey_name. - self.survey_name = survey_name + def __init__(self, scheduler_note=None, survey_name=None): + if scheduler_note is None and survey_name is not None: + self.scheduler_note = survey_name + else: + self.scheduler_note = scheduler_note # Start out with an empty observation self.feature = utils.empty_observation() def add_observations_array(self, observations_array, observations_hpid): - if self.survey_name is not None: - good = np.where(observations_array["scheduler_note"] == self.survey_name)[0] - if np.size(good) < 0: - self.feature = observations_array[good[-1]] + if self.scheduler_note is not None: + valid_indx = np.ones(observations_array.size, dtype=bool) + tmp = [self.scheduler_note in name for name in observations_array["scheduler_note"]] + valid_indx = valid_indx * np.array(tmp) + if len(tmp) > 0: + self.feature = observations_array[valid_indx][-1] else: if len(observations_array) > 0: self.feature = observations_array[-1] def add_observation(self, observation, indx=None): - if self.survey_name is not None: - if self.survey_name in observation["scheduler_note"]: + if self.scheduler_note is not None: + if self.scheduler_note in observation["scheduler_note"]: self.feature = observation else: self.feature = observation @@ -388,6 +397,7 @@ def __init__(self, sequence_ids=""): # observations... # Start out with an empty observation self.feature = utils.empty_observation() + send_unused_deprecation_warning(self.__class__.__name__) def add_observation(self, observation, indx=None): if observation["survey_id"] in self.sequence_ids: @@ -399,6 +409,7 @@ class LastFilterChange(BaseSurveyFeature): def __init__(self): self.feature = {"mjd": 0.0, "previous_filter": None, "current_filter": None} + send_unused_deprecation_warning(self.__class__.__name__) def add_observation(self, observation, indx=None): if self.feature["current_filter"] is None: @@ -417,28 +428,46 @@ class NObservations(BaseSurveyFeature): Parameters ---------- - filtername : `str` ('r') + filtername : `str` or `list` [`str`] or None String or list that has all the filters that can count. - nside : `int` (32) - The nside of the healpixel map to use - + Default None counts all filters. + nside : `int` + The nside of the healpixel map to use. + Default None uses scheduler default. + scheduler_note : `str` or None, optional + The scheduler_note to match. + Scheduler_note values which match this OR which contain this value + as a subset of their string will match. + survey_name : `str` or None + The scheduler_note value to match. + Deprecated in favor of scheduler_note, but provided for backward + compatibility. Will be removed in the future. """ - def __init__(self, filtername=None, nside=None, survey_name=None): + def __init__(self, filtername=None, nside=None, scheduler_note=None, survey_name=None): if nside is None: nside = utils.set_default_nside() self.feature = np.zeros(hp.nside2npix(nside), dtype=float) self.filtername = filtername - 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 + # This feature is not used with "scheduler_note" in the baseline + # survey. Should it really match scheduler_notes which contain + # self.scheduler_note or should it be vice versa? + # (i.e. this works as: self.scheduler_note = "survey" matches only + # scheduler_note == "survey", but self.scheduler_note = "survey a" + # will match both "survey" and "survey a". self.bins = np.arange(hp.nside2npix(nside) + 1) - 0.5 def add_observations_array(self, observations_array, observations_hpid): valid_indx = np.ones(observations_hpid.size, dtype=bool) if self.filtername is not None: valid_indx[np.where(observations_hpid["filter"] != self.filtername)[0]] = False - if self.survey_name is not None: - tmp = [name in self.survey_name for name in observations_hpid["scheduler_note"]] + if self.scheduler_note is not None: + tmp = [name in self.scheduler_note for name in observations_hpid["scheduler_note"]] valid_indx = valid_indx * np.array(tmp) data = observations_hpid[valid_indx] if np.size(data) > 0: @@ -449,7 +478,7 @@ def add_observations_array(self, observations_array, observations_hpid): def add_observation(self, observation, indx=None): if self.filtername is None or observation["filter"][0] in self.filtername: - if self.survey_name is None or observation["scheduler_note"] in self.survey_name: + if self.scheduler_note is None or observation["scheduler_note"] in self.scheduler_note: self.feature[indx] += 1 @@ -879,14 +908,13 @@ def add_observation(self, observation, indx=None): class NObsNight(BaseSurveyFeature): """ - Track how many times a healpixel has been observed in a night - (Note, even if there are two, it might not be a good pair.) + Track how many times a healpixel has been observed in a night. Parameters ---------- filtername : `str` or None Filter to track. None tracks observations in any filter. - nside : `int` or NOne + nside : `int` or None Scale of the healpix map. Default of None uses the scheduler default nside. """ diff --git a/rubin_scheduler/scheduler/model_observatory/kinem_model.py b/rubin_scheduler/scheduler/model_observatory/kinem_model.py index 30bc7488..188a835d 100644 --- a/rubin_scheduler/scheduler/model_observatory/kinem_model.py +++ b/rubin_scheduler/scheduler/model_observatory/kinem_model.py @@ -13,10 +13,80 @@ from .jerk import jerk_time -__all__ = ("KinemModel",) +__all__ = ( + "tma_movement", + "rotator_movement", + "KinemModel", +) two_pi = 2.0 * np.pi +def tma_movement(percent=70): + """Get a dictionary of parameters to pass to `setup_telescope` + defining altitude and azimuth speed, acceleration, and jerk + in terms of 'percent' of total performance. + + Parameters + ---------- + percent : `float`, optional + Default performance for the scheduler simulations for operations + has been 70% (70, default). + Expected performance at the start of comcam on-sky + science operations is about 10%. + + Returns + ------- + tma : `dict` {`str`: `float`} + A dictionary which can be passed as kwargs to + KinematicModel.setup_telescope(**tma). + """ + # See https://confluence.lsstcorp.org/display/LSSTCOM/TMA+Motion+Settings + # Expected performance at end of comcam on-sky is probably 10% + if percent > 125: + percent = 125 + print("Cannot exceed 125 percent, by requirements.") + tma = {} + scale = percent / 100.0 + tma["azimuth_maxspeed"] = np.min([10.0 * scale, 7.0]) + tma["azimuth_accel"] = 10.0 * scale + tma["azimuth_jerk"] = np.max([1.0, 40.0 * scale]) + tma["altitude_maxspeed"] = 5.0 * scale + tma["altitude_accel"] = 5.0 * scale + tma["altitude_jerk"] = np.max([1.0, 20.0 * scale]) + tma["settle_time"] = 3.0 + return tma + + +def rotator_movement(percent=100): + """Get a dictionary of parameters to pass to `setup_camera` + defining rotator max speed, acceleration and jerk, + in terms of 'percent' of total performance. + + Parameters + ---------- + percent : `float`, optional + Default performance for the scheduler simulations for operations + has been 100% (100, default). + Expected performance at the start of comcam on-sky + science operations is approximately full performance. + + Returns + ------- + rot : `dict` {`str`: `float`} + A dictionary which can be passed as kwargs to + KinematicModel.setup_camera(**rot). + """ + # Kevin and Brian say these can run 100%, are independent of TMA movement + if percent > 125: + percent = 125 + print("Cannot exceed 125 percent, by requirements.") + rot = {} + rot["maxspeed"] = 3.5 * percent / 100 + rot["accel"] = 1.0 * percent / 100 + rot["jerk"] = 4.0 * percent / 100 + return rot + + class Radec2altazpa: """Class to make it easy to swap in different alt/az conversion if wanted.""" @@ -92,6 +162,11 @@ def mount_filters(self, filter_list): List of the mounted filters. """ self.mounted_filters = filter_list + # Make sure we're using one of the available filters. + if self.current_filter not in self.mounted_filters: + self.current_filter = self.mounted_filters[-1] + if self.start_filter not in self.mounted_filters: + self.start_filter = self.mounted_filters[-1] def setup_camera( self, diff --git a/rubin_scheduler/scheduler/model_observatory/model_observatory.py b/rubin_scheduler/scheduler/model_observatory/model_observatory.py index 93692a07..e212172e 100644 --- a/rubin_scheduler/scheduler/model_observatory/model_observatory.py +++ b/rubin_scheduler/scheduler/model_observatory/model_observatory.py @@ -291,11 +291,12 @@ def __init__( good_mjd, to_set_mjd = self.check_mjd(to_set_mjd) self.mjd = to_set_mjd - sun_moon_info = self.almanac.get_sun_moon_positions(mjd) # Create the map of the season offsets - this map is constant ra, dec = _hpid2_ra_dec(nside, np.arange(hp.nside2npix(self.nside))) ra_deg = np.degrees(ra) self.season_map = calc_season(ra_deg, [self.mjd_start], self.mjd_start).flatten() + # Set the sun_ra_start information, for the rolling footprints + sun_moon_info = self.almanac.get_sun_moon_positions(self.mjd_start) self.sun_ra_start = sun_moon_info["sun_RA"] + 0 # Conditions object to update and return on request # (at present, this is not updated -- recreated, below). @@ -611,15 +612,14 @@ def _update_rot_sky_pos(self, observation): obs_pa = _approx_altaz2pa(alt, az, self.site.latitude_rad) - # If the observation has a rotTelPos set, use it to compute - # rotSkyPos + # If the observation has a rotTelPos set, use it to compute rotSkyPos if np.isfinite(observation["rotTelPos"]): observation["rotSkyPos"] = self.rc._rottelpos2rotskypos(observation["rotTelPos"], obs_pa) observation["rotTelPos"] = np.nan else: - # Fall back to rotSkyPos_desired + # Try to fall back to rotSkyPos_desired possible_rot_tel_pos = self.rc._rotskypos2rottelpos(observation["rotSkyPos_desired"], obs_pa) - + # If in range, use rotSkyPos_desired for rotSkyPos if (possible_rot_tel_pos > rot_limit[0]) | (possible_rot_tel_pos < rot_limit[1]): observation["rotSkyPos"] = observation["rotSkyPos_desired"] observation["rotTelPos"] = np.nan diff --git a/rubin_scheduler/scheduler/schedulers/core_scheduler.py b/rubin_scheduler/scheduler/schedulers/core_scheduler.py index 2722de35..8228903c 100644 --- a/rubin_scheduler/scheduler/schedulers/core_scheduler.py +++ b/rubin_scheduler/scheduler/schedulers/core_scheduler.py @@ -214,8 +214,17 @@ def request_observation(self, mjd=None): Returns ------- - observations : observation object (ra,dec,filter,rotangle) + observation : observation object (ra,dec,filter,rotangle) Returns None if the queue fails to fill + + Notes + ----- + Calling request_observation repeatedly, even without updating + the time or conditions, can return different requested + observations. This is because the internal queue is filled + when it becomes empty, and then subsequent calls to + request_observation will pop successive visits from this queue, + until the queue is empty or is flushed. """ if mjd is None: mjd = self.conditions.mjd @@ -286,8 +295,9 @@ def _fill_queue(self): for ns, surveys in enumerate(self.survey_lists): rewards = np.zeros(len(surveys)) for i, survey in enumerate(surveys): + # For each survey, find the highest reward value. rewards[i] = np.nanmax(survey.calc_reward_function(self.conditions)) - # If we have a good reward, break out of the loop + # If we have a tier with a good reward, break out of the loop if np.nanmax(rewards) > -np.inf: self.survey_index[0] = ns break @@ -305,6 +315,7 @@ def _fill_queue(self): ) self.queue = result + self.queue_filled = self.conditions.mjd if len(self.queue) == 0: self.log.warning(f"Failed to fill queue at time {self.conditions.mjd}") diff --git a/rubin_scheduler/scheduler/schedulers/filter_scheduler.py b/rubin_scheduler/scheduler/schedulers/filter_scheduler.py index cbe41203..58e88a84 100644 --- a/rubin_scheduler/scheduler/schedulers/filter_scheduler.py +++ b/rubin_scheduler/scheduler/schedulers/filter_scheduler.py @@ -1,4 +1,6 @@ -__all__ = ("FilterSwapScheduler", "SimpleFilterSched", "FilterSchedUzy") +__all__ = ("FilterSwapScheduler", "SimpleFilterSched", "ComCamFilterSched", "FilterSchedUzy") + +import numpy as np from rubin_scheduler.scheduler.utils import IntRounded @@ -33,6 +35,55 @@ def __call__(self, conditions): return result +class ComCamFilterSched(FilterSwapScheduler): + """ComCam can only hold 3 filters at a time. + + Pretend we will cycle from ugr, gri, riz, izy + depending on lunar phase. + + Parameters + ---------- + loaded_filter_groups : `tuple` (`tuple`) + Groups of 3 filters, to be loaded at the same time. + Multiple groups can be specified, to be swapped between at the + boundaries of `illum_bins` in lunar phase. + illum_bins : `np.ndarray`, (N,) + Lunar illumination boundaries to define when to swap between + different groups of filters within the `loaded_filter_groups`. + Lunar illumination ranges from 0 to 100. + + Notes + ----- + If illum_bins = np.array([0, 50, 100]), then there should be + two groups of filters to use -- one for use between 0 and 50 percent + illumination, and another for use between 50 and 100 percent illumination. + """ + + def __init__( + self, + loaded_filter_groups=(("u", "g", "r"), ("g", "r", "i"), ("r", "i", "z"), ("i", "z", "y")), + illum_bins=np.arange(0, 100 + 1, 25), + ): + self.loaded_filter_groups = loaded_filter_groups + self.illum_bins = illum_bins + if isinstance(self.illum_bins, list): + self.illum_bins = np.array(illum_bins) + if len(illum_bins) - 1 > len(loaded_filter_groups): + raise ValueError("There are illumination bins with an " "undefined loaded_filter_group") + + def __call__(self, conditions): + moon_at_sunset = conditions.moon_phase_sunset + try: + if len(moon_at_sunset) > 0: + moon_at_sunset = moon_at_sunset[0] + except TypeError: + pass + indx = np.searchsorted(self.illum_bins, moon_at_sunset, side="left") + indx = np.max([0, indx - 1]) + result = list(self.loaded_filter_groups[indx]) + return result + + class FilterSchedUzy(FilterSwapScheduler): """ remove u in bright time. Alternate between removing z and y in diff --git a/rubin_scheduler/scheduler/sim_runner.py b/rubin_scheduler/scheduler/sim_runner.py index ecd1174e..530aa20e 100644 --- a/rubin_scheduler/scheduler/sim_runner.py +++ b/rubin_scheduler/scheduler/sim_runner.py @@ -1,5 +1,6 @@ __all__ = ("sim_runner",) +import copy import sqlite3 import sys import time @@ -86,7 +87,6 @@ def sim_runner( step_none = step_none / 60.0 / 24.0 # to days mjd_run = end_mjd - mjd_start nskip = 0 - new_night = False mjd_last_flush = -1 @@ -141,11 +141,14 @@ def sim_runner( # An observation failed to execute, usually it was outside # the altitude limits. if observatory.mjd == mjd_last_flush: - raise RuntimeError("Scheduler has failed to provide a valid observation multiple times.") + raise RuntimeError( + "Scheduler has failed to provide a valid observation multiple times " + f" at time ({observatory.mjd} from survey {scheduler.survey_index}." + ) # if this is a first offence, might just be that targets set. - # Flush queue and get some new targets. + # Flush queue and try to get some new targets. scheduler.flush_queue() - mjd_last_flush = observatory.mjd + 0 + mjd_last_flush = copy.deepcopy(observatory.mjd) if new_night: # find out what filters we want mounted conditions = observatory.return_conditions() @@ -185,8 +188,8 @@ def sim_runner( ) observations["alt"] = np.radians(alt) observations["az"] = np.radians(az) - observations["psudo_pa"] = np.radians(pa) - observations["rotTelPos"] = rc._rotskypos2rottelpos(observations["rotSkyPos"], observations["psudo_pa"]) + observations["pseudo_pa"] = np.radians(pa) + observations["rotTelPos"] = rc._rotskypos2rottelpos(observations["rotSkyPos"], observations["pseudo_pa"]) # Also include traditional parallactic angle pa = _approx_altaz2pa(observations["alt"], observations["az"], lsst.latitude_rad) diff --git a/rubin_scheduler/scheduler/surveys/base_survey.py b/rubin_scheduler/scheduler/surveys/base_survey.py index 6a4cc63b..4dd45820 100644 --- a/rubin_scheduler/scheduler/surveys/base_survey.py +++ b/rubin_scheduler/scheduler/surveys/base_survey.py @@ -188,7 +188,11 @@ def add_observation(self, observation, **kwargs): def _check_feasibility(self, conditions): """ - Check if the survey is feasable in the current conditions + Check if the survey is feasible in the current conditions + + Returns + ------- + result : `bool` """ result = True for bf in self.basis_functions: diff --git a/rubin_scheduler/scheduler/surveys/field_survey.py b/rubin_scheduler/scheduler/surveys/field_survey.py index c8a97637..d9023019 100644 --- a/rubin_scheduler/scheduler/surveys/field_survey.py +++ b/rubin_scheduler/scheduler/surveys/field_survey.py @@ -8,7 +8,7 @@ from rubin_scheduler.utils import ra_dec2_hpid -from ..features import NObsSurvey +from ..features import LastObservation, NObsSurvey from ..utils import empty_observation from . import BaseSurvey @@ -136,6 +136,9 @@ def __init__( survey_name=self.survey_name, ) self.accept_obs = accept_obs + if isinstance(self.accept_obs, str): + self.accept_obs = [self.accept_obs] + self.indx = ra_dec2_hpid(self.nside, self.ra_deg, self.dec_deg) # Set all basis function equal. self.basis_weights = np.ones(len(basis_functions)) / len(basis_functions) @@ -206,8 +209,9 @@ def __init__( self.indx = ra_dec2_hpid(self.nside, self.ra_deg, self.dec_deg) # Tucking this here so we can look at how many observations - # recorded for this field - self.extra_features["ObsRecord"] = NObsSurvey() + # recorded for this field and what was the last one. + self.extra_features["ObsRecorded"] = NObsSurvey() + self.extra_features["LastObs"] = LastObservation() def _generate_survey_name(self): if self.target_name is not None: @@ -237,7 +241,7 @@ def check_continue(self, observation, conditions): def add_observation(self, observation, **kwargs): """Add observation one at a time.""" - # Check each posible ignore string + # Check each possible ignore string checks = [io not in str(observation["scheduler_note"]) for io in self.ignore_obs] passed_ignore = all(checks) passed_accept = True @@ -287,12 +291,16 @@ def add_observations_array(self, observations_array_in, observations_hpid_in): not_ignore = np.where(np.char.find(observations_hpid["scheduler_note"], ig) == -1)[0] observations_hpid = observations_hpid[not_ignore] - for acc in self.accept_obs: - accept = np.where(np.char.find(observations_array["scheduler_note"], acc) == 1)[0] + if self.accept_obs is not None: + accept_indx = [] + accept_hp_indx = [] + for acc in self.accept_obs: + accept_indx.append(np.where(observations_array["scheduler_note"] == acc)[0]) + accept_hp_indx.append(np.where(observations_hpid["scheduler_note"] == acc)[0]) + accept = np.concatenate(accept_indx) + accept_hp = np.concatenate(accept_hp_indx) observations_array = observations_array[accept] - - accept = np.where(np.char.find(observations_hpid["scheduler_note"], acc) == 1)[0] - observations_hpid = observations_hpid[accept] + observations_hpid = observations_hpid[accept_hp] for feature in self.extra_features: self.extra_features[feature].add_observations_array(observations_array, observations_hpid) @@ -305,6 +313,7 @@ def add_observations_array(self, observations_array_in, observations_hpid_in): self.reward_checked = False def calc_reward_function(self, conditions): + # only calculates reward at the index for the RA/Dec of the field self.reward_checked = True if self._check_feasibility(conditions): self.reward = 0 diff --git a/rubin_scheduler/scheduler/surveys/long_gap_survey.py b/rubin_scheduler/scheduler/surveys/long_gap_survey.py index 3d08fe53..15e487a8 100644 --- a/rubin_scheduler/scheduler/surveys/long_gap_survey.py +++ b/rubin_scheduler/scheduler/surveys/long_gap_survey.py @@ -94,13 +94,13 @@ def __init__( self.mjd_step = hour_step / 24.0 def _generate_survey_name(self): - self.survey_name = ( - f"Long Gap ({self.blob_survey.survey_name} +" f" {self.scripted_survey.survey_name})" - ) + self.blob_survey.survey_name = f"Long Gap {self.blob_survey.survey_name}" + self.scripted_survey.survey_name = "Long Gap Scripted Triplet" + self.survey_name = f"Long Gap ({self.blob_survey.survey_name} + Triplet)" def _schedule_obs(self, observations): """Take incoming observations and decide if they should be added - to the scripted survey to try and be observered again later + to the scripted survey to try and be observed again later """ # Only match if we have completed the second of a pair and are @@ -194,8 +194,10 @@ def _schedule_obs(self, observations): sched_array["target_name"] = "" sched_array["observation_reason"] = "FBS" sched_array["json_block"] = "Imaging" - # See if we need to append things to the scripted survey - # object + # Don't let the desired rotSkyPos block the observation. + sched_array["rotSkyPos_desired"] = sched_array["rotSkyPos"] + sched_array["rotSkyPos"] = np.nan + # See if we need to append things to the scripted survey object if self.scripted_survey.obs_wanted is not None: sched_array = np.concatenate([self.scripted_survey.obs_wanted, sched_array]) diff --git a/rubin_scheduler/scheduler/surveys/scripted_surveys.py b/rubin_scheduler/scheduler/surveys/scripted_surveys.py index 4c1a4463..94d97b56 100644 --- a/rubin_scheduler/scheduler/surveys/scripted_surveys.py +++ b/rubin_scheduler/scheduler/surveys/scripted_surveys.py @@ -312,7 +312,7 @@ def set_script(self, obs_wanted): self.id_start = np.max(self.obs_wanted["scripted_id"]) + 1 self.mjd_start = self.obs_wanted["mjd"] - self.obs_wanted["mjd_tol"] - # Here is the atribute that core scheduler checks to + # Here is the attribute that core scheduler checks to # broadcast scheduled observations in the conditions object. self.scheduled_obs = self.obs_wanted["mjd"] @@ -334,7 +334,7 @@ def generate_observations_rough(self, conditions): if self.before_twi_check: # Note that if detailers are adding lots of exposures, this # calculation has the potential to not be right at all. - # Also assumes slew time is negligable. + # Also assumes slew time is negligible. exptime_needed = np.sum(observations["exptime"]) / 3600.0 / 24.0 # to days filter_change_needed = n_filter_changes * self.filter_change_time tot_time_needed = exptime_needed + filter_change_needed diff --git a/rubin_scheduler/scheduler/surveys/surveys.py b/rubin_scheduler/scheduler/surveys/surveys.py index 5e0534cb..cb934480 100644 --- a/rubin_scheduler/scheduler/surveys/surveys.py +++ b/rubin_scheduler/scheduler/surveys/surveys.py @@ -37,6 +37,7 @@ def __init__( fields=None, ): extra_features = {} + self.filtername = filtername self.block_size = block_size self.nexp = nexp @@ -119,14 +120,18 @@ class BlobSurvey(GreedySurvey): The approximate slewtime between neerby fields (seconds). Used to calculate how many observations can be taken in the desired time block. + filter_change_approx : `float` + The approximate time it takes to change filters (seconds). + read_approx : `float` + The approximate time required to readout the camera (seconds). + exptime : `float` + The total on-sky exposure time per visit. nexp : `int` The number of exposures to take in a visit. exp_dict : `dict` If set, should have keys of filtername and values of ints that are the nuber of exposures to take per visit. For estimating block time, nexp is still used. - filter_change_approx : `float` - The approximate time it takes to change filters (seconds). ideal_pair_time : `float` The ideal time gap wanted between observations to the same pointing (minutes) @@ -187,13 +192,30 @@ def __init__( area_required=None, max_radius_peak=40.0, fields=None, - **kwargs, + survey_note=None, + search_radius=None, + alt_max=-9999, + az_range=-9999, ): + if search_radius is not None: + warnings.warn("search_radius unused, remove kwarg", DeprecationWarning, 2) + if alt_max != -9999: + warnings.warn("alt_max unused, remove kwarg", DeprecationWarning, 2) + if az_range != -9999: + warnings.warn("az_range unused, remove kwarg", DeprecationWarning, 2) + self.filtername1 = filtername1 self.filtername2 = filtername2 + self.ideal_pair_time = ideal_pair_time + if survey_name is None: - survey_name = self._generate_survey_name() + self._generate_survey_name() + + if scheduler_note is None: + if survey_note is not None: + scheduler_note = survey_note + warnings.warn("survey_note is deprecated - use scheduler_note", DeprecationWarning, 2) super(BlobSurvey, self).__init__( basis_functions=basis_functions, @@ -209,7 +231,7 @@ def __init__( camera=camera, area_required=area_required, fields=fields, - survey_name=survey_name, + survey_name=self.survey_name, scheduler_note=scheduler_note, ) self.flush_time = flush_time / 60.0 / 24.0 # convert to days @@ -237,10 +259,14 @@ def __init__( # then repeat.) if filtername2 is not None: self.time_needed = ( - (ideal_pair_time * 60.0 * 2.0 + exptime + read_approx + filter_change_approx) / 24.0 / 3600.0 + (self.ideal_pair_time * 60.0 * 2.0 + self.exptime + self.read_approx + filter_change_approx) + / 24.0 + / 3600.0 ) # Days else: - self.time_needed = (ideal_pair_time * 60.0 + exptime + read_approx) / 24.0 / 3600.0 # Days + self.time_needed = ( + (self.ideal_pair_time * 60.0 + self.exptime + self.read_approx) / 24.0 / 3600.0 + ) # Days self.filter_set = set(filtername1) if filtername2 is None: self.filter2_set = self.filter_set @@ -250,7 +276,6 @@ def __init__( self.ra, self.dec = _hpid2_ra_dec(self.nside, self.hpids) self.counter = 1 # start at 1, because 0 is default in empty obs - self.ideal_pair_time = ideal_pair_time self.pixarea = hp.nside2pixarea(self.nside, degrees=True) @@ -259,7 +284,9 @@ def __init__( self.filtername = self.filtername1 def _generate_survey_name(self): - self.survey_name = f"Blob survey {self.filtername1}" + self.survey_name = "Pairs" + self.survey_name += f" {self.ideal_pair_time :.1f}" + self.survey_name += f" {self.filtername1}" if self.filtername2 is None: self.survey_name += f"_{self.filtername1}" else: diff --git a/rubin_scheduler/scheduler/utils/utils.py b/rubin_scheduler/scheduler/utils/utils.py index ea5a93f3..aca407b2 100644 --- a/rubin_scheduler/scheduler/utils/utils.py +++ b/rubin_scheduler/scheduler/utils/utils.py @@ -401,7 +401,7 @@ def __init__(self): "slewTime": "slewtime", "slewDistance": "slewdist", "paraAngle": "pa", - "psudoParaAngle": "psudo_pa", + "pseudoParaAngle": "pseudo_pa", "rotTelPos": "rotTelPos", "rotTelPos_backup": "rotTelPos_backup", "rotSkyPos": "rotSkyPos", @@ -418,6 +418,9 @@ def __init__(self): "note": "note", "scheduler_note": "scheduler_note", "target_name": "target_name", + "science_program": "science_program", + "observation_reason": "observation_reason", + "json_block": "json_block", } # For backwards compatibility self.backwards = {"target": "target_name"} @@ -432,7 +435,7 @@ def __init__(self): "azimuth", "slewDistance", "paraAngle", - "psudoParaAngle", + "pseudoParaAngle", "rotTelPos", "rotSkyPos", "rotSkyPos_desired", @@ -623,7 +626,7 @@ def empty_observation(n=1): (for scheduler purposes, use `scheduler_note`). This maps to observation_reason in the ConsDB, although could be overwritten in JSON BLOCK. - Most likely this is just "science" when using the FBS. + Most likely this is just "science" or "FBS" when using the FBS. json_block : `str` (optional) The JSON BLOCK id to use to acquire observations. This is for use by the SchedulerCSC. @@ -670,7 +673,7 @@ def empty_observation(n=1): ("alt", float), ("az", float), ("pa", float), - ("psudo_pa", float), + ("pseudo_pa", float), ("clouds", float), ("moonAlt", float), ("sunAlt", float), diff --git a/tests/scheduler/test_basisfuncs.py b/tests/scheduler/test_basisfuncs.py index 04685f74..784c5e1e 100644 --- a/tests/scheduler/test_basisfuncs.py +++ b/tests/scheduler/test_basisfuncs.py @@ -197,7 +197,7 @@ def test_deprecated(self): warnings.simplefilter("always") dep_bf() # Verify deprecation warning - assert len(w) == 1 + assert len(w) >= 1 assert issubclass(w[-1].category, (DeprecationWarning, FutureWarning)) diff --git a/tests/scheduler/test_comcam_surveys.py b/tests/scheduler/test_comcam_surveys.py new file mode 100644 index 00000000..9cd17978 --- /dev/null +++ b/tests/scheduler/test_comcam_surveys.py @@ -0,0 +1,76 @@ +import unittest + +import numpy as np +from astropy.time import Time + +from rubin_scheduler.scheduler import sim_runner +from rubin_scheduler.scheduler.example import ( + get_comcam_sv_schedulers, + get_model_observatory, + get_sv_fields, + update_model_observatory_sunset, +) +from rubin_scheduler.scheduler.schedulers import ComCamFilterSched +from rubin_scheduler.utils import survey_start_mjd + + +class TestComCamSurveys(unittest.TestCase): + + def test_model_observatory_conveniences(self): + """Test the model observatory convenience functions.""" + + # Just check that we can acquire a model observatory + # and it is set up for the date expected. + survey_start = survey_start_mjd() + survey_start = np.floor(survey_start) + 0.5 + dayobs = Time(survey_start, format="mjd", scale="utc").iso[:10] + observatory = get_model_observatory(dayobs=dayobs, survey_start=survey_start) + conditions = observatory.return_conditions() + assert conditions.mjd == observatory.mjd + # The model observatory automatically advanced to -12 deg sunset + assert (conditions.mjd - survey_start) < 1 + sun_ra_start = conditions.sun_ra_start + mjd_start = observatory.mjd_start + + newday = survey_start + 4 + new_dayobs = Time(newday, format="mjd", scale="utc").iso[:10] + newday = Time(f"{new_dayobs}T12:00:00", format="isot", scale="utc").mjd + observatory = get_model_observatory(dayobs=new_dayobs, survey_start=survey_start) + conditions = observatory.return_conditions() + assert (conditions.mjd - newday) < 1 + # Check that advancing the day did not change the expected location + # of the sun at the *start* of the survey + assert conditions.mjd_start == mjd_start + assert conditions.sun_ra_start == sun_ra_start + + # And update observatory to sunset, using a filter scheduler + # that only has 'g' available + filter_sched = ComCamFilterSched(illum_bins=np.arange(0, 101, 100), loaded_filter_groups=(("g",))) + observatory = update_model_observatory_sunset(observatory, filter_sched, twilight=-18) + assert observatory.observatory.current_filter == "g" + assert observatory.conditions.sun_alt < np.radians(18) + + def test_comcam_sv_sched(self): + """Test the comcam sv survey scheduler setup.""" + # This is likely to change as we go into commissioning, + # so mostly I'm just going to test that the end result is + # a usable scheduler + survey_start = survey_start_mjd() + survey_start = np.floor(survey_start) + 0.5 + dayobs = Time(survey_start, format="mjd", scale="utc").iso[:10] + scheduler, filter_scheduler = get_comcam_sv_schedulers() + observatory = get_model_observatory(dayobs=dayobs, survey_start=survey_start) + observatory = update_model_observatory_sunset(observatory, filter_scheduler) + + observatory, scheduler, observations = sim_runner( + observatory, scheduler, filter_scheduler, survey_length=30 + ) + assert len(observations) > 24000 + sv_fields = set(np.unique(observations["scheduler_note"])) + all_sv_fields = set(list(get_sv_fields().keys())) + # Probably won't observe all of the fields, but can do some. + assert sv_fields.intersection(all_sv_fields) == sv_fields + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/scheduler/test_features.py b/tests/scheduler/test_features.py index 98204df5..4bbb5a6d 100644 --- a/tests/scheduler/test_features.py +++ b/tests/scheduler/test_features.py @@ -10,6 +10,52 @@ from rubin_scheduler.utils import survey_start_mjd +def make_observations_list(nobs=1): + observations_list = [] + for i in range(0, nobs): + observation = empty_observation() + observation["mjd"] = survey_start_mjd() + i * 30 / 60 / 60 / 24 + observation["RA"] = np.radians(30) + observation["dec"] = np.radians(-20) + observation["filter"] = "r" + observation["scheduler_note"] = "test" + observations_list.append(observation) + return observations_list + + +def make_observations_arrays(observations_list, nside=32): + # Turn list of observations (that should already have useful info) + # into observations_array plus observations_hpids_array. + observations_array = np.empty(len(observations_list), dtype=observations_list[0].dtype) + for i, obs in enumerate(observations_list): + observations_array[i] = obs + # Build observations_hpids_array. + # Find list of lists of healpixels + # (should match [indxs, indxs, indxs2] from above) + pointing2indx = HpInLsstFov(nside=nside) + list_of_hpids = pointing2indx(observations_array["RA"], observations_array["dec"]) + # Unravel list-of-lists (list_of_hpids) to match against observations + hpids = [] + big_array_indx = [] + for i, indxs in enumerate(list_of_hpids): + for indx in indxs: + hpids.append(indx) + big_array_indx.append(i) + hpids = np.array(hpids, dtype=[("hpid", int)]) + # Set up format / dtype for observations_hpids_array + names = list(observations_array.dtype.names) + types = [observations_array[name].dtype for name in names] + names.append(hpids.dtype.names[0]) + types.append(hpids["hpid"].dtype) + ndt = list(zip(names, types)) + observations_hpids_array = np.empty(hpids.size, dtype=ndt) + # Populate observations_hpid_array - big_array_indx points + # between index in observations_array and index in hpid + observations_hpids_array[list(observations_array.dtype.names)] = observations_array[big_array_indx] + observations_hpids_array[hpids.dtype.names[0]] = hpids + return observations_array, observations_hpids_array + + class TestFeatures(unittest.TestCase): def test_pair_in_night(self): pin = features.PairInNight(gap_min=25.0, gap_max=45.0) @@ -283,6 +329,156 @@ def test_NObservationsCurrentSeason(self): # in these cases with added requirements .. but will leave it # to the "restore" test in test_utils.py. + def test_NObsSurvey(self): + # Make some observations to count + observations_list = make_observations_list(2) + observations_list[0]["scheduler_note"] = "survey a" + observations_list[1]["scheduler_note"] = "survey b" + observations_array, observations_hpid_array = make_observations_arrays(observations_list) + # Count the observations matching any note + count_feature = features.NObsSurvey(note=None) + # ... it matters significantly that we pass obs[0] and not obs. + for obs in observations_list: + count_feature.add_observation(obs[0]) + self.assertTrue(count_feature.feature == 2) + # and count again using add_observations_array + count_feature = features.NObsSurvey(note=None) + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature == 2) + # Count using a note to match + # Count the observations matching specific note + count_feature = features.NObsSurvey(note="survey a") + for obs in observations_list: + count_feature.add_observation(obs[0]) + self.assertTrue(count_feature.feature == 1) + # and count again using add_observations_array + count_feature = features.NObsSurvey(note="survey a") + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature == 1) + # Count the observations matching subset of note + count_feature = features.NObsSurvey(note="survey") + for obs in observations_list: + count_feature.add_observation(obs[0]) + self.assertTrue(count_feature.feature == 2) + # and count again using add_observations_array + count_feature = features.NObsSurvey(note="survey") + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature == 2) + + def test_LastObservation(self): + # Make some observations to count + observations_list = make_observations_list(4) + observations_list[0]["mjd"] = survey_start_mjd() + observations_list[1]["mjd"] = survey_start_mjd() + 10 + observations_list[2]["mjd"] = survey_start_mjd() + 20 + observations_list[3]["mjd"] = survey_start_mjd() + 30 + observations_list[0]["scheduler_note"] = "survey a" + observations_list[1]["scheduler_note"] = "survey" + observations_list[2]["scheduler_note"] = "survey a" + observations_list[3]["scheduler_note"] = "survey b" + observations_array, observations_hpid_array = make_observations_arrays(observations_list) + # Observations matching any note + count_feature = features.LastObservation(scheduler_note=None) + for obs in observations_list: + count_feature.add_observation(obs[0]) + self.assertTrue(count_feature.feature["mjd"] == observations_list[-1]["mjd"]) + # and count again using add_observations_array + count_feature = features.LastObservation(scheduler_note=None) + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature["mjd"] == observations_list[-1]["mjd"]) + # Observations matching a specific note. + count_feature = features.LastObservation(scheduler_note="survey a") + for obs in observations_list: + count_feature.add_observation(obs[0]) + self.assertTrue(count_feature.feature["mjd"] == observations_list[-2]["mjd"]) + # and count again using add_observations_array + count_feature = features.LastObservation(scheduler_note="survey a") + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature["mjd"] == observations_list[-2]["mjd"]) + # Observations matching a subset of note. + count_feature = features.LastObservation(scheduler_note="survey") + for obs in observations_list: + count_feature.add_observation(obs[0]) + self.assertTrue(count_feature.feature["mjd"] == observations_list[-1]["mjd"]) + # and count again using add_observations_array + count_feature = features.LastObservation(scheduler_note="survey") + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature["mjd"] == observations_list[-1]["mjd"]) + + def test_NObservations(self): + # Make some observations to count + observations_list = make_observations_list(12) + indexes = [] + nside = 32 + pointing2hpindx = HpInLsstFov(nside=nside) + for i, obs in enumerate(observations_list): + obs["mjd"] = survey_start_mjd() + i + obs["rotSkyPos"] = 0 + if i < 6: + obs["filter"] = "r" + else: + obs["filter"] = "g" + if i % 2 == 0: + obs["RA"] = np.radians(30) + else: + obs["RA"] = np.radians(10) + if i % 3 == 0: + obs["scheduler_note"] = "survey a" + elif i % 3 == 1: + obs["scheduler_note"] = "survey b" + else: + obs["scheduler_note"] = "survey" + indexes.append(pointing2hpindx(obs["RA"], obs["dec"], rotSkyPos=obs["rotSkyPos"])) + observations_array, observations_hpid_array = make_observations_arrays(observations_list) + # Observations matching any note or filter + count_feature = features.NObservations(filtername=None, scheduler_note=None) + for obs, indx in zip(observations_list, indexes): + count_feature.add_observation(obs[0], indx) + self.assertTrue(count_feature.feature.max() == 6) + # and count again using add_observations_array + count_feature = features.NObservations(filtername=None, scheduler_note=None) + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature.max() == 6) + # Observations matching a specific note - or are partial matches. + count_feature = features.NObservations(scheduler_note="survey a") + for obs, indx in zip(observations_list, indexes): + count_feature.add_observation(obs[0], indx) + self.assertTrue(count_feature.feature.max() == 4) + # and count again using add_observations_array + # Observations matching a specific note. + count_feature = features.NObservations(scheduler_note="survey a") + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature.max() == 4) + # Observations matching a subset of note. + # It's not obvious that this is what this SHOULD do, and it's not + # used with "note" in the example /baseline scheduler. + count_feature = features.NObservations(scheduler_note="survey") + for obs, indx in zip(observations_list, indexes): + count_feature.add_observation(obs[0], indx) + self.assertTrue(count_feature.feature.max() == 2) + # and count again using add_observations_array + count_feature = features.NObservations(scheduler_note="survey") + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature.max() == 2) + # Observations matching any note but specified filter. + count_feature = features.NObservations(filtername="r", scheduler_note=None) + for obs, indx in zip(observations_list, indexes): + count_feature.add_observation(obs[0], indx) + self.assertTrue(count_feature.feature.max() == 3) + # and count again using add_observations_array + count_feature = features.NObservations(filtername="r", scheduler_note=None) + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature.max() == 3) + # Observations matching specific note and specified filter. + count_feature = features.NObservations(filtername="r", scheduler_note="survey") + for obs, indx in zip(observations_list, indexes): + count_feature.add_observation(obs[0], indx) + self.assertTrue(count_feature.feature.max() == 1) + # and count again using add_observations_array + count_feature = features.NObservations(filtername="r", scheduler_note="survey") + count_feature.add_observations_array(observations_array, observations_hpid=observations_hpid_array) + self.assertTrue(count_feature.feature.max() == 1) + if __name__ == "__main__": unittest.main() diff --git a/tests/scheduler/test_filterschedulers.py b/tests/scheduler/test_filterschedulers.py new file mode 100644 index 00000000..09558939 --- /dev/null +++ b/tests/scheduler/test_filterschedulers.py @@ -0,0 +1,52 @@ +import unittest + +import numpy as np + +from rubin_scheduler.scheduler.features import Conditions +from rubin_scheduler.scheduler.schedulers import ComCamFilterSched, SimpleFilterSched +from rubin_scheduler.utils import survey_start_mjd + + +class TestFilterSchedulers(unittest.TestCase): + + def test_ComCamFilterSched(self): + illum_bins = np.arange(0, 100 + 1, 50) + filter_groups = (("g", "r", "i"), ("i", "z", "y")) + filtersched = ComCamFilterSched(illum_bins=illum_bins, loaded_filter_groups=filter_groups) + mjd = survey_start_mjd() + conditions = Conditions(nside=8, mjd_start=mjd, mjd=mjd) + conditions.moon_phase_sunset = 0 + load_filters = filtersched(conditions) + self.assertTrue(load_filters == ["g", "r", "i"]) + conditions.moon_phase_sunset = 40 + load_filters = filtersched(conditions) + self.assertTrue(load_filters == ["g", "r", "i"]) + conditions.moon_phase_sunset = 60 + load_filters = filtersched(conditions) + self.assertTrue(load_filters == ["i", "z", "y"]) + conditions.moon_phase_sunset = 100 + load_filters = filtersched(conditions) + self.assertTrue(load_filters == ["i", "z", "y"]) + + def test_comcamfiltersched_except(self): + illum_bins = np.arange(0, 100 + 1, 25) + filter_groups = (("g", "r", "i"), ("i", "z", "y")) + with self.assertRaises(ValueError): + ComCamFilterSched(illum_bins=illum_bins, loaded_filter_groups=filter_groups) + + def test_SimpleFilterSched(self): + filtersched = SimpleFilterSched(illum_limit=40) + brightmoon_result = ["g", "r", "i", "z", "y"] + newmoon_result = ["u", "g", "r", "i", "z"] + mjd = survey_start_mjd() + conditions = Conditions(nside=8, mjd_start=mjd, mjd=mjd) + conditions.moon_phase_sunset = 0 + load_filters = filtersched(conditions) + self.assertTrue(load_filters == newmoon_result) + conditions.moon_phase_sunset = 50 + load_filters = filtersched(conditions) + self.assertTrue(load_filters == brightmoon_result) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/scheduler/test_surveys.py b/tests/scheduler/test_surveys.py index 53c4803c..23b5de17 100644 --- a/tests/scheduler/test_surveys.py +++ b/tests/scheduler/test_surveys.py @@ -8,7 +8,54 @@ import rubin_scheduler.scheduler.surveys as surveys from rubin_scheduler.scheduler.basis_functions import SimpleArrayBasisFunction from rubin_scheduler.scheduler.model_observatory import ModelObservatory -from rubin_scheduler.scheduler.utils import empty_observation, set_default_nside +from rubin_scheduler.scheduler.utils import HpInLsstFov, empty_observation, set_default_nside +from rubin_scheduler.utils import survey_start_mjd + + +def make_observations_list(nobs=1): + observations_list = [] + for i in range(0, nobs): + observation = empty_observation() + observation["mjd"] = survey_start_mjd() + i * 30 / 60 / 60 / 24 + observation["RA"] = np.radians(30) + observation["dec"] = np.radians(-20) + observation["filter"] = "r" + observation["scheduler_note"] = "test" + observations_list.append(observation) + return observations_list + + +def make_observations_arrays(observations_list, nside=32): + # Turn list of observations (that should already have useful info) + # into observations_array plus observations_hpids_array. + observations_array = np.empty(len(observations_list), dtype=observations_list[0].dtype) + for i, obs in enumerate(observations_list): + observations_array[i] = obs + # Build observations_hpids_array. + # Find list of lists of healpixels + # (should match [indxs, indxs, indxs2] from above) + pointing2indx = HpInLsstFov(nside=nside) + list_of_hpids = pointing2indx(observations_array["RA"], observations_array["dec"]) + # Unravel list-of-lists (list_of_hpids) to match against observations + hpids = [] + big_array_indx = [] + for i, indxs in enumerate(list_of_hpids): + for indx in indxs: + hpids.append(indx) + big_array_indx.append(i) + hpids = np.array(hpids, dtype=[("hpid", int)]) + # Set up format / dtype for observations_hpids_array + names = list(observations_array.dtype.names) + types = [observations_array[name].dtype for name in names] + names.append(hpids.dtype.names[0]) + types.append(hpids["hpid"].dtype) + ndt = list(zip(names, types)) + observations_hpids_array = np.empty(hpids.size, dtype=ndt) + # Populate observations_hpid_array - big_array_indx points + # between index in observations_array and index in hpid + observations_hpids_array[list(observations_array.dtype.names)] = observations_array[big_array_indx] + observations_hpids_array[hpids.dtype.names[0]] = hpids + return observations_array, observations_hpids_array class TestSurveys(unittest.TestCase): @@ -34,6 +81,64 @@ def test_field_survey(self): self.assertIsInstance(reward_df, pd.DataFrame) reward_df = survey.make_reward_df(conditions, accum=False) + def test_field_survey_add_observations(self): + nside = 32 + # Just need a placeholder for this. + bfs = [basis_functions.M5DiffBasisFunction(nside=nside)] + + observations_list = make_observations_list(10) + indexes = [] + pointing2hpindx = HpInLsstFov(nside=nside) + for i, obs in enumerate(observations_list): + obs["mjd"] = survey_start_mjd() + i + obs["rotSkyPos"] = 0 + if i < 5: + obs["filter"] = "r" + obs["scheduler_note"] = "r band" + else: + obs["filter"] = "g" + obs["scheduler_note"] = "g band" + obs["RA"] = np.radians(90) + obs["dec"] = np.radians(-30) + indexes.append(pointing2hpindx(obs["RA"], obs["dec"], rotSkyPos=obs["rotSkyPos"])) + observations_array, observations_hpid_array = make_observations_arrays(observations_list) + + # Try adding observations to survey one at a time. + survey = surveys.FieldSurvey(bfs, RA=90.0, dec=-30.0, accept_obs=None) + for obs, indx in zip(observations_list, indexes): + survey.add_observation(obs[0], indx=indx) + self.assertTrue(survey.extra_features["ObsRecorded"].feature == len(observations_list)) + self.assertTrue(survey.extra_features["LastObs"].feature["mjd"] == observations_list[-1]["mjd"]) + # Try adding observations to survey in array + survey = surveys.FieldSurvey(bfs, RA=90.0, dec=-30.0, accept_obs=None) + survey.add_observations_array(observations_array, observations_hpid_array) + self.assertTrue(survey.extra_features["ObsRecorded"].feature == len(observations_list)) + self.assertTrue(survey.extra_features["LastObs"].feature["mjd"] == observations_list[-1]["mjd"]) + + # Now with specific requirements on obs to accept. + # Try adding observations to survey one at a time. + survey = surveys.FieldSurvey(bfs, RA=90.0, dec=-30.0, accept_obs=["r band"]) + for obs, indx in zip(observations_list, indexes): + survey.add_observation(obs[0], indx=indx) + self.assertTrue(survey.extra_features["ObsRecorded"].feature == 5) + self.assertTrue(survey.extra_features["LastObs"].feature["mjd"] == observations_list[4]["mjd"]) + # Try adding observations to survey in array + survey = surveys.FieldSurvey(bfs, RA=90.0, dec=-30.0, accept_obs=["r band"]) + survey.add_observations_array(observations_array, observations_hpid_array) + self.assertTrue(survey.extra_features["ObsRecorded"].feature == 5) + self.assertTrue(survey.extra_features["LastObs"].feature["mjd"] == observations_list[4]["mjd"]) + # Try adding observations to survey one at a time. + survey = surveys.FieldSurvey(bfs, RA=90.0, dec=-30.0, accept_obs=["r band", "g band"]) + for obs, indx in zip(observations_list, indexes): + survey.add_observation(obs[0], indx=indx) + self.assertTrue(survey.extra_features["ObsRecorded"].feature == 10) + self.assertTrue(survey.extra_features["LastObs"].feature["mjd"] == observations_list[-1]["mjd"]) + # Try adding observations to survey in array + survey = surveys.FieldSurvey(bfs, RA=90.0, dec=-30.0, accept_obs=["r band", "g band"]) + survey.add_observations_array(observations_array, observations_hpid_array) + self.assertTrue(survey.extra_features["ObsRecorded"].feature == 10) + self.assertTrue(survey.extra_features["LastObs"].feature["mjd"] == observations_list[-1]["mjd"]) + def test_pointings_survey(self): """Test the pointing survey.""" mo = ModelObservatory() diff --git a/tests/utils/test_rotsky.py b/tests/utils/test_rotsky.py index 8d906dee..e92c6e63 100644 --- a/tests/utils/test_rotsky.py +++ b/tests/utils/test_rotsky.py @@ -59,8 +59,8 @@ def test_rotation_converter(self): with self.assertRaises(ValueError): rc = rotation_converter(telescope="not_a_telescope_name") - def test_psudo_pa(self): - # Check that the psudo parallactic angle is + def test_pseudo_pa(self): + # Check that the pseudo parallactic angle is # somewhat close to the approx parallactic angle lsst = Site("LSST") rng = np.random.default_rng(seed=42)