From 8f1446a105931ed7935394f2e47288b272814a90 Mon Sep 17 00:00:00 2001 From: Peter Yoachim Date: Wed, 25 Nov 2020 15:52:21 -0800 Subject: [PATCH 1/5] adding solar elongation --- .../basis_functions/mask_basis_funcs.py | 31 +++++++++++++++++-- .../featureScheduler/features/conditions.py | 14 ++++++++- 2 files changed, 41 insertions(+), 4 deletions(-) diff --git a/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py b/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py index 5a152b4..eccb8f8 100644 --- a/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py +++ b/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py @@ -8,7 +8,32 @@ __all__ = ['Zenith_mask_basis_function', 'Zenith_shadow_mask_basis_function', 'Moon_avoidance_basis_function', 'Map_cloud_basis_function', - 'Planet_mask_basis_function', 'Mask_azimuth_basis_function'] + 'Planet_mask_basis_function', 'Mask_azimuth_basis_function', 'Solar_elongation_mask_basis_function'] + + +class Solar_elongation_mask_basis_function(Base_basis_function): + """Mask things at various solar elongations + + Parameters + ---------- + min_elong : float (0) + The minimum solar elongation to consider (degrees). + max_elong : float (60.) + The maximum solar elongation to consider (degrees). + """ + + def __init__(self, min_elong=0., max_elong=60., nside=None): + super(Solar_elongation_mask_basis_function, self).__init__(nside=nside) + self.min_elong = np.radians(min_elong) + self.max_elong = np.radians(max_elong) + self.result = np.empty(hp.nside2npix(self.nside), dtype=float).fill(self.penalty) + + def _calc_value(self, conditions, indx=None): + result = self.result.copy() + in_range = np.where((int_rounded(conditions.solar_elongation) >= int_rounded(self.min_elong)) & + (int_rounded(conditions.solar_elongation) <= int_rounded(self.max_elong)))[0] + result[in_range] = 1 + return result class Zenith_mask_basis_function(Base_basis_function): @@ -21,8 +46,8 @@ class Zenith_mask_basis_function(Base_basis_function): max_alt : float (82.) The maximum allowed altitude (degrees) """ - def __init__(self, min_alt=20., max_alt=82.): - super(Zenith_mask_basis_function, self).__init__() + def __init__(self, min_alt=20., max_alt=82., nside=None): + super(Zenith_mask_basis_function, self).__init__(nside=nside) self.update_on_newobs = False self.min_alt = np.radians(min_alt) self.max_alt = np.radians(max_alt) diff --git a/python/lsst/sims/featureScheduler/features/conditions.py b/python/lsst/sims/featureScheduler/features/conditions.py index f05c31e..ccf2b2c 100644 --- a/python/lsst/sims/featureScheduler/features/conditions.py +++ b/python/lsst/sims/featureScheduler/features/conditions.py @@ -1,5 +1,5 @@ import numpy as np -from lsst.sims.utils import _approx_RaDec2AltAz, Site, _hpid2RaDec, m5_flat_sed, _approx_altaz2pa +from lsst.sims.utils import _approx_RaDec2AltAz, Site, _hpid2RaDec, m5_flat_sed, _approx_altaz2pa, _angularSeparation import healpy as hp from lsst.sims.featureScheduler.utils import set_default_nside, match_hp_resolution, season_calc, smallest_signed_angle @@ -151,6 +151,8 @@ def __init__(self, nside=None, site='LSST', exptime=30., mjd_start=59853.5, seas Healpix map of the azimuthal distance to the sun for each healpixel (radians) az_to_anitsun : np.array Healpix map of the azimuthal distance to the anit-sun for each healpixel (radians) + solar_elongation : np.array + Healpix map of the solar elongation (angular distance to the sun) for each healpixel (radians) Attributes (set by the scheduler) ------------------------------- @@ -338,6 +340,7 @@ def mjd(self, value): self._az_to_sun = None self._az_to_antisun = None self._season = None + self._solar_elongation = None @property def skybrightness(self): @@ -378,6 +381,15 @@ def calc_M5Depth(self): self.exptime, self._airmass[good]) + def calc_solar_elongation(self): + self._solar_elongation = _angularSeparation(self.ra, self.dec, self.runRA, self.sunDec) + + @property + def solar_elongation(self): + if self._solar_elongation is None: + self.calc_solar_elongation() + return self._solar_elongation + def calc_az_to_sun(self): self._az_to_sun = smallest_signed_angle(self.ra, self.sunRA) From 2e8a86f12d173a087f42a1f57376605a6f0f17f7 Mon Sep 17 00:00:00 2001 From: Peter Yoachim Date: Wed, 25 Nov 2020 16:20:36 -0800 Subject: [PATCH 2/5] adding area minimum requirements to surveys --- .../featureScheduler/surveys/base_survey.py | 17 +++++++++++++++-- .../sims/featureScheduler/surveys/surveys.py | 16 ++++++++++++---- 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/python/lsst/sims/featureScheduler/surveys/base_survey.py b/python/lsst/sims/featureScheduler/surveys/base_survey.py index e1b5e37..5d254c2 100644 --- a/python/lsst/sims/featureScheduler/surveys/base_survey.py +++ b/python/lsst/sims/featureScheduler/surveys/base_survey.py @@ -174,11 +174,13 @@ class BaseMarkovDF_survey(BaseSurvey): Random number seed, used for randomly orienting sky tessellation. camera : str ('LSST') Should be 'LSST' or 'comcam' + area_required : float (None) + The valid area that should be present in the reward function (square degrees). """ def __init__(self, basis_functions, basis_weights, extra_features=None, smoothing_kernel=None, ignore_obs=None, survey_name='', nside=None, seed=42, - dither=True, detailers=None, camera='LSST'): + dither=True, detailers=None, camera='LSST', area_required=None): super(BaseMarkovDF_survey, self).__init__(basis_functions=basis_functions, extra_features=extra_features, @@ -197,7 +199,7 @@ def __init__(self, basis_functions, basis_weights, extra_features=None, elif self.camera == 'comcam': self.fields_init = comcamTessellate() else: - ValueError('camera %s unknown, should be "LSST" or "comcam"' %camera) + ValueError('camera %s unknown, should be "LSST" or "comcam"' % camera) self.fields = self.fields_init.copy() self.hp2fields = np.array([]) self._hp2fieldsetup(self.fields['RA'], self.fields['dec']) @@ -207,6 +209,11 @@ def __init__(self, basis_functions, basis_weights, extra_features=None, else: self.smoothing_kernel = None + if area_required is None: + self.area_required = area_required + else: + self.area_required = area_required * (np.pi/180.)**2 # To steradians + # Start tracking the night self.night = -1 @@ -307,6 +314,12 @@ def calc_reward_function(self, conditions): if self.smoothing_kernel is not None: self.smooth_reward() + if self.area_required is not None: + good_area = np.where(np.abs(self.reward) >= 0)[0].size * hp.nside2pixarea(self.nside) + if good_area < self.area_required: + self.reward = -np.inf + + return self.reward def generate_observations_rough(self, conditions): diff --git a/python/lsst/sims/featureScheduler/surveys/surveys.py b/python/lsst/sims/featureScheduler/surveys/surveys.py index a88b956..1576861 100644 --- a/python/lsst/sims/featureScheduler/surveys/surveys.py +++ b/python/lsst/sims/featureScheduler/surveys/surveys.py @@ -19,7 +19,7 @@ class Greedy_survey(BaseMarkovDF_survey): def __init__(self, basis_functions, basis_weights, filtername='r', block_size=1, smoothing_kernel=None, nside=None, dither=True, seed=42, ignore_obs=None, survey_name='', - nexp=2, exptime=30., detailers=None, camera='LSST'): + nexp=2, exptime=30., detailers=None, camera='LSST', area_required=None): extra_features = {} @@ -30,7 +30,8 @@ def __init__(self, basis_functions, basis_weights, filtername='r', ignore_obs=ignore_obs, nside=nside, survey_name=survey_name, dither=dither, - detailers=detailers, camera=camera) + detailers=detailers, camera=camera, + area_required=area_required) self.filtername = filtername self.block_size = block_size self.nexp = nexp @@ -132,7 +133,7 @@ def __init__(self, basis_functions, basis_weights, dither=True, seed=42, ignore_obs=None, survey_note='blob', detailers=None, camera='LSST', twilight_scale=True, in_twilight=False, check_scheduled=True, min_area=None, - grow_blob=True): + grow_blob=True, area_required=None): if nside is None: nside = set_default_nside() @@ -142,7 +143,8 @@ def __init__(self, basis_functions, basis_weights, filtername=None, block_size=0, smoothing_kernel=smoothing_kernel, dither=dither, seed=seed, ignore_obs=ignore_obs, - nside=nside, detailers=detailers, camera=camera) + nside=nside, detailers=detailers, camera=camera, + area_required=area_required) self.flush_time = flush_time/60./24. # convert to days self.nexp = nexp self.nexp_dict = nexp_dict @@ -305,6 +307,12 @@ def calc_reward_function(self, conditions): self.reward[az_out] = np.nan else: self.reward = -np.inf + + if self.area_required is not None: + good_area = np.where(np.abs(self.reward) >= 0)[0].size * hp.nside2pixarea(self.nside) + if good_area < self.area_required: + self.reward = -np.inf + self.reward_checked = True return self.reward From 3b50460599b17aeaccef20f426ce74ba50e6d33d Mon Sep 17 00:00:00 2001 From: Peter Yoachim Date: Wed, 9 Dec 2020 16:47:15 -0800 Subject: [PATCH 3/5] working on twilight neo --- .../basis_functions/basis_functions.py | 7 ++++--- .../basis_functions/mask_basis_funcs.py | 6 ++++-- .../featureScheduler/features/conditions.py | 2 +- python/lsst/sims/featureScheduler/sim_runner.py | 2 +- .../featureScheduler/surveys/base_survey.py | 16 ++++++++++++++++ .../sims/featureScheduler/utils/footprints.py | 17 ++++++++++++++++- 6 files changed, 42 insertions(+), 8 deletions(-) diff --git a/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py b/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py index 9330a73..519f51c 100644 --- a/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py +++ b/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py @@ -564,9 +564,10 @@ def __init__(self, nside=None, max_airmass=2.5): def _calc_value(self, conditions, indx=None): result = self.result.copy() - good_pix = np.where((int_rounded(conditions.airmass) < self.max_airmass) & - (int_rounded(conditions.az_to_sun) < int_rounded(np.pi/2.))) - result[good_pix] = conditions.airmass[good_pix] / self.max_airmass.value + good_pix = np.where((conditions.airmass >= 1.) & + (int_rounded(conditions.airmass) < self.max_airmass) & + (int_rounded(np.abs(conditions.az_to_sun)) < int_rounded(np.pi/2.))) + result[good_pix] = conditions.airmass[good_pix] / self.max_airmass.initial return result diff --git a/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py b/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py index eccb8f8..958308c 100644 --- a/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py +++ b/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py @@ -22,11 +22,13 @@ class Solar_elongation_mask_basis_function(Base_basis_function): The maximum solar elongation to consider (degrees). """ - def __init__(self, min_elong=0., max_elong=60., nside=None): + def __init__(self, min_elong=0., max_elong=60., nside=None, penalty=np.nan): super(Solar_elongation_mask_basis_function, self).__init__(nside=nside) self.min_elong = np.radians(min_elong) self.max_elong = np.radians(max_elong) - self.result = np.empty(hp.nside2npix(self.nside), dtype=float).fill(self.penalty) + self.penalty = penalty + self.result = np.empty(hp.nside2npix(self.nside), dtype=float) + self.result.fill(self.penalty) def _calc_value(self, conditions, indx=None): result = self.result.copy() diff --git a/python/lsst/sims/featureScheduler/features/conditions.py b/python/lsst/sims/featureScheduler/features/conditions.py index ccf2b2c..adb123b 100644 --- a/python/lsst/sims/featureScheduler/features/conditions.py +++ b/python/lsst/sims/featureScheduler/features/conditions.py @@ -382,7 +382,7 @@ def calc_M5Depth(self): self._airmass[good]) def calc_solar_elongation(self): - self._solar_elongation = _angularSeparation(self.ra, self.dec, self.runRA, self.sunDec) + self._solar_elongation = _angularSeparation(self.ra, self.dec, self.sunRA, self.sunDec) @property def solar_elongation(self): diff --git a/python/lsst/sims/featureScheduler/sim_runner.py b/python/lsst/sims/featureScheduler/sim_runner.py index 174d37c..2b10f29 100644 --- a/python/lsst/sims/featureScheduler/sim_runner.py +++ b/python/lsst/sims/featureScheduler/sim_runner.py @@ -88,7 +88,7 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su if len(observations) == n_visit_limit: break # XXX--handy place to interupt and debug - # if len(observations) > 3: + #if len(observations) > 25: # import pdb ; pdb.set_trace() runtime = time.time() - t0 print('Skipped %i observations' % nskip) diff --git a/python/lsst/sims/featureScheduler/surveys/base_survey.py b/python/lsst/sims/featureScheduler/surveys/base_survey.py index 5d254c2..a75ba8d 100644 --- a/python/lsst/sims/featureScheduler/surveys/base_survey.py +++ b/python/lsst/sims/featureScheduler/surveys/base_survey.py @@ -221,6 +221,22 @@ def __init__(self, basis_functions, basis_weights, extra_features=None, np.random.seed(seed) self.dither = dither + def _check_feasibility(self, conditions): + """ + Check if the survey is feasable in the current conditions + """ + for bf in self.basis_functions: + result = bf.check_feasibility(conditions) + if not result: + return result + if self.area_required is not None: + reward = self.calc_reward_function(conditions) + good_pix = np.where(np.isfinite(reward) == True)[0] + area = hp.nside2pixarea(self.nside) * np.size(good_pix) + if area < self.area_required: + return False + return result + def _hp2fieldsetup(self, ra, dec, leafsize=100): """Map each healpixel to nearest field. This will only work if healpix resolution is higher than field resolution. diff --git a/python/lsst/sims/featureScheduler/utils/footprints.py b/python/lsst/sims/featureScheduler/utils/footprints.py index e9210c3..272e51a 100644 --- a/python/lsst/sims/featureScheduler/utils/footprints.py +++ b/python/lsst/sims/featureScheduler/utils/footprints.py @@ -22,7 +22,7 @@ 'WFD_healpixels', 'WFD_no_gp_healpixels', 'WFD_bigsky_healpixels', 'WFD_no_dust_healpixels', 'SCP_healpixels', 'NES_healpixels', 'galactic_plane_healpixels', #'low_lat_plane_healpixels', 'bulge_healpixels', - 'magellanic_clouds_healpixels', + 'magellanic_clouds_healpixels', 'Constant_footprint', 'generate_goal_map', 'standard_goals', 'calc_norm_factor', 'filter_count_ratios', 'Step_line', 'Footprints', 'Footprint', 'Step_slopes'] @@ -176,6 +176,21 @@ def __call__(self, mjd, array=False): return self.arr2struc(self.current_footprints) +class Constant_footprint(Footprint): + def __init__(self, nside=32, + filters={'u': 0, 'g': 1, 'r': 2, 'i': 3, 'z': 4, 'y': 5}): + self.nside = nside + self.filters = filters + self.npix = hp.nside2npix(nside) + self.footprints = np.zeros((len(filters), self.npix), dtype=float) + self.out_dtype = list(zip(filters, [float]*len(filters))) + self.to_return = self.arr2struc(self.footprints) + + def __call__(self, mjd, array=False): + return self.to_return + + + class Footprints(Footprint): """An object to combine multiple Footprint objects. """ From 64875b02011fd7f6638fe366499fbd3bfb70475c Mon Sep 17 00:00:00 2001 From: Peter Yoachim Date: Thu, 10 Dec 2020 15:05:33 -0800 Subject: [PATCH 4/5] looks like neo twilight working better --- .../basis_functions/basis_functions.py | 7 ++-- .../basis_functions/mask_basis_funcs.py | 36 ++++++++++++++++++- .../modelObservatory/model_observatory.py | 2 +- .../lsst/sims/featureScheduler/sim_runner.py | 2 ++ .../featureScheduler/surveys/base_survey.py | 1 - .../sims/featureScheduler/surveys/surveys.py | 3 ++ 6 files changed, 45 insertions(+), 6 deletions(-) diff --git a/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py b/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py index 519f51c..1a0b8d3 100644 --- a/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py +++ b/python/lsst/sims/featureScheduler/basis_functions/basis_functions.py @@ -74,7 +74,7 @@ def check_feasibility(self, conditions): """ return True - def _calc_value(self, conditions, **kwarge): + def _calc_value(self, conditions, **kwargs): self.value = 0 # Update the last time we had an mjd self.mjd_last = conditions.mjd + 0 @@ -557,10 +557,11 @@ class Near_sun_twilight_basis_function(Base_basis_function): The maximum airmass to try and observe (unitless) """ - def __init__(self, nside=None, max_airmass=2.5): + def __init__(self, nside=None, max_airmass=2.5, penalty=np.nan): super(Near_sun_twilight_basis_function, self).__init__(nside=nside) self.max_airmass = int_rounded(max_airmass) - self.result = np.zeros(hp.nside2npix(self.nside)) + self.result = np.empty(hp.nside2npix(self.nside)) + self.result.fill(penalty) def _calc_value(self, conditions, indx=None): result = self.result.copy() diff --git a/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py b/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py index 958308c..1cc0384 100644 --- a/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py +++ b/python/lsst/sims/featureScheduler/basis_functions/mask_basis_funcs.py @@ -8,7 +8,41 @@ __all__ = ['Zenith_mask_basis_function', 'Zenith_shadow_mask_basis_function', 'Moon_avoidance_basis_function', 'Map_cloud_basis_function', - 'Planet_mask_basis_function', 'Mask_azimuth_basis_function', 'Solar_elongation_mask_basis_function'] + 'Planet_mask_basis_function', 'Mask_azimuth_basis_function', + 'Solar_elongation_mask_basis_function', 'Area_check_mask_basis_function'] + + + +class Area_check_mask_basis_function(Base_basis_function): + """Take a list of other mask basis functions, and do an additional check for area available + """ + def __init__(self, bf_list, nside=32, min_area=1000.): + super(Area_check_mask_basis_function, self).__init__(nside=nside) + self.bf_list = bf_list + self.result = np.zeros(hp.nside2npix(self.nside), dtype=float) + self.min_area = min_area + + def check_feasibility(self, conditions): + result = True + for bf in self.bf_list: + if not bf.check_feasibility(conditions): + return False + + area_map = self.result.copy() + for bf in self.bf_list: + area_map *= bf(conditions) + + good_pix = np.where(area_map == 0)[0] + if hp.nside2pixarea(self.nside, degrees=True)*good_pix.size < self.min_area: + result = False + return result + + def _calc_value(self, conditions, **kwargs): + result = self.result.copy() + for bf in self.bf_list: + result *= bf(conditions) + return result + class Solar_elongation_mask_basis_function(Base_basis_function): diff --git a/python/lsst/sims/featureScheduler/modelObservatory/model_observatory.py b/python/lsst/sims/featureScheduler/modelObservatory/model_observatory.py index 84fba5a..f589662 100644 --- a/python/lsst/sims/featureScheduler/modelObservatory/model_observatory.py +++ b/python/lsst/sims/featureScheduler/modelObservatory/model_observatory.py @@ -443,7 +443,7 @@ def observe(self, observation): # obsevation['rotSkyPos'] will be ignored. slewtime, visittime = self.observatory.observe(observation, self.mjd, rotTelPos=observation['rotTelPos']) - # inf slewtime means the observation failed (probably outsire alt limits) + # inf slewtime means the observation failed (probably outside alt limits) if ~np.all(np.isfinite(slewtime)): return None, False diff --git a/python/lsst/sims/featureScheduler/sim_runner.py b/python/lsst/sims/featureScheduler/sim_runner.py index 2b10f29..562b13c 100644 --- a/python/lsst/sims/featureScheduler/sim_runner.py +++ b/python/lsst/sims/featureScheduler/sim_runner.py @@ -69,6 +69,8 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su observations.append(completed_obs) filter_scheduler.add_observation(completed_obs[0]) else: + # XXX--note this can result in an infinite loop, where if a survey requests something bad + # say, outside the alt limits. The queue will flush and refill, potentially with bad observations again. scheduler.flush_queue() if new_night: # find out what filters we want mounted diff --git a/python/lsst/sims/featureScheduler/surveys/base_survey.py b/python/lsst/sims/featureScheduler/surveys/base_survey.py index a75ba8d..240f9e7 100644 --- a/python/lsst/sims/featureScheduler/surveys/base_survey.py +++ b/python/lsst/sims/featureScheduler/surveys/base_survey.py @@ -335,7 +335,6 @@ def calc_reward_function(self, conditions): if good_area < self.area_required: self.reward = -np.inf - return self.reward def generate_observations_rough(self, conditions): diff --git a/python/lsst/sims/featureScheduler/surveys/surveys.py b/python/lsst/sims/featureScheduler/surveys/surveys.py index 1576861..1d19a22 100644 --- a/python/lsst/sims/featureScheduler/surveys/surveys.py +++ b/python/lsst/sims/featureScheduler/surveys/surveys.py @@ -313,6 +313,9 @@ def calc_reward_function(self, conditions): if good_area < self.area_required: self.reward = -np.inf + #if ('twi' in self.survey_note) & (np.any(np.isfinite(self.reward))): + # import pdb ; pdb.set_trace() + self.reward_checked = True return self.reward From a096f214bcce0658a9fbbdde547d77fde38f8753 Mon Sep 17 00:00:00 2001 From: Peter Yoachim Date: Fri, 11 Dec 2020 14:50:02 -0800 Subject: [PATCH 5/5] insert break for inf loop --- python/lsst/sims/featureScheduler/sim_runner.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/python/lsst/sims/featureScheduler/sim_runner.py b/python/lsst/sims/featureScheduler/sim_runner.py index 562b13c..830ffd1 100644 --- a/python/lsst/sims/featureScheduler/sim_runner.py +++ b/python/lsst/sims/featureScheduler/sim_runner.py @@ -52,6 +52,8 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su nskip = 0 new_night = False + mjd_last_flush = -1 + while mjd < end_mjd: if not scheduler._check_queue_mjd_only(observatory.mjd): scheduler.update_conditions(observatory.return_conditions()) @@ -69,9 +71,12 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su observations.append(completed_obs) filter_scheduler.add_observation(completed_obs[0]) else: - # XXX--note this can result in an infinite loop, where if a survey requests something bad - # say, outside the alt limits. The queue will flush and refill, potentially with bad observations again. + # 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.") + # if this is a first offence, might just be that targets set. Flush queue and get some new targets. scheduler.flush_queue() + mjd_last_flush = observatory.mjd + 0 if new_night: # find out what filters we want mounted conditions = observatory.return_conditions()