Merge pull request #15 from lsst/tickets/OPSIM-893

Tickets/opsim 893

rubin_scheduler/scheduler/basis_functions/basis_functions.py

@@ -874,12 +874,15 @@ def __init__(self, nside=None, max_airmass=2.5, penalty=np.nan):
 
     def _calc_value(self, conditions, indx=None):
         result = self.result.copy()
+        valid_airmass = np.isfinite(conditions.airmass)
         good_pix = np.where(
-            (conditions.airmass >= 1.0)
-            & (IntRounded(conditions.airmass) < self.max_airmass)
-            & (IntRounded(np.abs(conditions.az_to_sun)) < IntRounded(np.pi / 2.0))
+            (conditions.airmass[valid_airmass] >= 1.0)
+            & (IntRounded(conditions.airmass[valid_airmass]) < self.max_airmass)
+            & (IntRounded(np.abs(conditions.az_to_sun[valid_airmass])) < IntRounded(np.pi / 2.0))
+        )
+        result[valid_airmass][good_pix] = (
+            conditions.airmass[valid_airmass][good_pix] / self.max_airmass.initial
         )
-        result[good_pix] = conditions.airmass[good_pix] / self.max_airmass.initial
         return result
 
 
@@ -906,12 +909,12 @@ def __init__(self, gap_min=25.0, gap_max=45.0, filtername="r", nside=None, npair
         self.gap_min = IntRounded(gap_min / 60.0 / 24.0)
         self.gap_max = IntRounded(gap_max / 60.0 / 24.0)
         self.npairs = npairs
-
         self.survey_features = {}
         # Track the number of pairs that have been taken in a night
         self.survey_features["Pair_in_night"] = features.PairInNight(
             filtername=filtername, gap_min=gap_min, gap_max=gap_max, nside=nside
         )
+
         # When was it last observed
         # XXX--since this feature is also in Pair_in_night, I should just access that one!
         self.survey_features["Last_observed"] = features.LastObserved(filtername=filtername, nside=nside)
@@ -920,11 +923,16 @@ def _calc_value(self, conditions, indx=None):
         result = np.zeros(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 = conditions.mjd - self.survey_features["Last_observed"].feature[indx]
+        mask = np.isnan(diff)
+        # remove NaNs from diff, but save mask so we exclude those values later.
+        diff[mask] = 0.0
+        ir_diff = IntRounded(diff)
         good = np.where(
-            (diff >= self.gap_min)
-            & (diff <= self.gap_max)
+            (ir_diff >= self.gap_min)
+            & (ir_diff <= self.gap_max)
             & (self.survey_features["Pair_in_night"].feature[indx] < self.npairs)
+            & (~mask)
         )[0]
         result[indx[good]] += 1.0
         return result

rubin_scheduler/scheduler/features/features.py

@@ -265,7 +265,8 @@ def add_observations_array(self, observations_array, observations_hpid):
             if np.size(good) < 0:
                 self.feature = observations_array[good[-1]]
         else:
-            self.feature = observations_array[-1]
+            if len(observations_array) > 0:
+                self.feature = observations_array[-1]
 
     def add_observation(self, observation, indx=None):
         if self.survey_name is not None:
@@ -693,8 +694,8 @@ def __init__(self, filtername="r", nside=None, gap_min=25.0, gap_max=45.0):
         self.feature = np.zeros(hp.nside2npix(nside), dtype=float)
         self.indx = np.arange(self.feature.size)
         self.last_observed = LastObserved(filtername=filtername)
-        self.gap_min = gap_min / (24.0 * 60)  # Days
-        self.gap_max = gap_max / (24.0 * 60)  # Days
+        self.gap_min = IntRounded(gap_min / (24.0 * 60))  # Days
+        self.gap_max = IntRounded(gap_max / (24.0 * 60))  # Days
         self.night = 0
         # Need to keep a full record of times and healpixels observed in a night.
         self.mjd_log = []
@@ -725,21 +726,22 @@ def add_observation(self, observation, indx=None):
                 self.mjd_log = []
                 self.hpid_log = []
 
-            # record the mjds and healpixels that were observed
-            self.mjd_log.extend([np.max(observation["mjd"])] * np.size(indx))
-            self.hpid_log.extend(list(indx))
-
             # Look for the mjds that could possibly pair with observation
-            tmin = observation["mjd"] - self.gap_max
-            tmax = observation["mjd"] - self.gap_min
-            mjd_log = np.array(self.mjd_log)
-            left = np.searchsorted(mjd_log, tmin).max()
-            right = np.searchsorted(mjd_log, tmax, side="right").max()
+            mjd_diff = IntRounded(observation["mjd"] - np.array(self.mjd_log))
+            # normally would use np.searchsorted, but need to use IntRounded
+            # to be sure we are cross-platform repeatable.
+            in_range_indx = np.where((mjd_diff > self.gap_min) & (mjd_diff < self.gap_max))[0]
             # Now check if any of the healpixels taken in the time gap
             # match the healpixels of the observation.
-            matches = np.in1d(indx, self.hpid_log[int(left) : int(right)])
-            # XXX--should think if this is the correct (fastest) order to check things in.
-            self.feature[np.array(indx)[matches]] += 1
+            if in_range_indx.size > 0:
+                left = in_range_indx.min()
+                right = in_range_indx.max() + 1
+                matches = np.in1d(indx, self.hpid_log[left:right])
+                self.feature[np.array(indx)[matches]] += 1
+
+            # record the mjds and healpixels that were observed
+            self.mjd_log.extend([np.max(observation["mjd"])] * np.size(indx))
+            self.hpid_log.extend(list(indx))
 
 
 class RotatorAngle(BaseSurveyFeature):

rubin_scheduler/scheduler/schedulers/core_scheduler.py

@@ -240,12 +240,13 @@ def request_observation(self, mjd=None):
                 )
                 obs_pa = _approx_altaz2pa(alt, az, self.conditions.site.latitude_rad)
                 rot_tel_pos_expected = (obs_pa - observation["rotSkyPos"]) % (2.0 * np.pi)
-                if (IntRounded(rot_tel_pos_expected) > IntRounded(self.rotator_limits[0])) & (
-                    IntRounded(rot_tel_pos_expected) < IntRounded(self.rotator_limits[1])
-                ):
-                    diff = np.abs(self.rotator_limits - rot_tel_pos_expected)
-                    limit_indx = np.min(np.where(diff == np.min(diff))[0])
-                    observation["rotSkyPos"] = (obs_pa - self.rotator_limits[limit_indx]) % (2.0 * np.pi)
+                if np.isfinite(observation["rotSkyPos"]):
+                    if (IntRounded(rot_tel_pos_expected) > IntRounded(self.rotator_limits[0])) & (
+                        IntRounded(rot_tel_pos_expected) < IntRounded(self.rotator_limits[1])
+                    ):
+                        diff = np.abs(self.rotator_limits - rot_tel_pos_expected)
+                        limit_indx = np.min(np.where(diff == np.min(diff))[0])
+                        observation["rotSkyPos"] = (obs_pa - self.rotator_limits[limit_indx]) % (2.0 * np.pi)
             return observation
 
     def _fill_queue(self):

rubin_scheduler/scheduler/utils/utils.py

@@ -72,6 +72,9 @@ class IntRounded:
     Class to help force comparisons be made on scaled up integers,
     preventing machine precision issues cross-platforms
 
+    Note that casting a NaN to an int can have different behaviours
+    cross-platforms, so will throw an error if attempted.
+
     Parameters
     ----------
     inval : number-like thing
@@ -81,6 +84,8 @@ class IntRounded:
     """
 
     def __init__(self, inval, scale=1e5):
+        if np.any(~np.isfinite(inval)):
+            raise ValueError("IntRounded can only take finite values.")
         self.initial = inval
         self.value = np.round(inval * scale).astype(int)
         self.scale = scale