Merge pull request #82 from lsst/u/yoachim/sr_fix

U/yoachim/sr fix

python/lsst/sims/featureScheduler/basis_functions/basis_functions.py

@@ -7,6 +7,8 @@
 import matplotlib.pylab as plt
 import warnings
 from lsst.sims.utils import _hpid2RaDec
+from astropy.coordinates import SkyCoord
+from astropy import units as u
 
 
 __all__ = ['Base_basis_function', 'Constant_basis_function', 'Target_map_basis_function',
@@ -19,7 +21,8 @@
            'Template_generate_basis_function',
            'Footprint_nvis_basis_function', 'Third_observation_basis_function', 'Season_coverage_basis_function',
            'N_obs_per_year_basis_function', 'Cadence_in_season_basis_function', 'Near_sun_twilight_basis_function',
-           'N_obs_high_am_basis_function', 'Good_seeing_basis_function', 'Observed_twice_basis_function']
+           'N_obs_high_am_basis_function', 'Good_seeing_basis_function', 'Observed_twice_basis_function',
+           'Ecliptic_basis_function']
 
 
 class Base_basis_function(object):
@@ -252,6 +255,24 @@ def _calc_value(self, conditions, indx=None):
         return result
 
 
+class Ecliptic_basis_function(Base_basis_function):
+    """Mark the area around the ecliptic
+    """
+
+    def __init__(self, nside=None, distance_to_eclip=25.):
+        super(Ecliptic_basis_function, self).__init__(nside=nside)
+        self.distance_to_eclip = np.radians(distance_to_eclip)
+        ra, dec = _hpid2RaDec(nside, np.arange(hp.nside2npix(self.nside)))
+        self.result = np.zeros(ra.size)
+        coord = SkyCoord(ra=ra*u.rad, dec=dec*u.rad)
+        eclip_lat = coord.barycentrictrueecliptic.lat.radian
+        good = np.where(np.abs(eclip_lat) < self.distance_to_eclip)
+        self.result[good] += 1
+
+    def __call__(self, conditions, indx=None):
+        return self.result
+
+
 class N_obs_per_year_basis_function(Base_basis_function):
     """Reward areas that have not been observed N-times in the last year
 
@@ -1211,13 +1232,13 @@ def __init__(self, nside=None, filtername='r', footprint=None, FWHMeff_limit=0.8
 
     def _calc_value(self, conditions, **kwargs):
         # Seeing is "bad"
-        if conditions.FWHMeff[self.filtername].min() > self.FWHMeff_limit:
+        if int_rounded(conditions.FWHMeff[self.filtername].min()) > self.FWHMeff_limit:
             return 0
         result = self.result.copy()
 
-        diff = int_rounded(self.survey_features['coadd_depth_all'].feature - self.survey_features['coadd_depth_good'].feature)
+        diff = self.survey_features['coadd_depth_all'].feature - self.survey_features['coadd_depth_good'].feature
         # Where are there things we want to observe?
-        good_pix = np.where((diff > self.mag_diff) &
+        good_pix = np.where((int_rounded(diff) > self.mag_diff) &
                             (int_rounded(conditions.FWHMeff[self.filtername]) <= self.FWHMeff_limit))
         # Hm, should this scale by the mag differences? Probably.
         result[good_pix] = diff[good_pix]

python/lsst/sims/featureScheduler/basis_functions/rolling_funcs.py

@@ -5,20 +5,23 @@
 import matplotlib.pylab as plt
 import warnings
 from lsst.sims.featureScheduler.basis_functions import Base_basis_function
+from lsst.sims.utils import _hpid2RaDec
+
 
 
 __all__ = ["Target_map_modulo_basis_function", "Footprint_basis_function", "Footprint_rolling_basis_function"]
 
 
+
 class Footprint_basis_function(Base_basis_function):
     """Basis function that tries to maintain a uniformly covered footprint
 
     Parameters
     ----------
     filtername : str ('r')
         The filter for this footprint
-    footprint : HEALpix np.array
-        The desired footprint. Assumed normalized.
+    footprint : lsst.sims.featureScheduler.utils.Footprint object
+        The desired footprint.
     all_footprints_sum : float (None)
         If using multiple filters, the sum of all the footprints. Needed to make sure basis functions are
         normalized properly across all fitlers.
@@ -27,33 +30,28 @@ class Footprint_basis_function(Base_basis_function):
         another good value to use)
 
     """
-    def __init__(self, filtername='r', nside=None, footprint=None, all_footprints_sum=None,
-                 out_of_bounds_val=-10.):
+    def __init__(self, filtername='r', nside=None, footprint=None,
+                 out_of_bounds_val=-10., window_size=6.):
 
         super(Footprint_basis_function, self).__init__(nside=nside, filtername=filtername)
         self.footprint = footprint
 
-        if all_footprints_sum is None:
-            # Assume the footprints are similar in weight
-            self.all_footprints_sum = np.sum(footprint)*6
-        else:
-            self.all_footprints_sum = all_footprints_sum
-
-        self.footprint_sum = np.sum(footprint)
-
         self.survey_features = {}
         # All the observations in all filters
         self.survey_features['N_obs_all'] = features.N_observations(nside=nside, filtername=None)
         self.survey_features['N_obs'] = features.N_observations(nside=nside, filtername=filtername)
 
         # should probably actually loop over all the target maps?
-        self.out_of_bounds_area = np.where(footprint <= 0)[0]
+        self.out_of_bounds_area = np.where(footprint.get_footprint(self.filtername) <= 0)[0]
         self.out_of_bounds_val = out_of_bounds_val
 
     def _calc_value(self, conditions, indx=None):
 
+        # Find out what the footprint object thinks we should have been observed
+        desired_footprint_normed = self.footprint(conditions.mjd)[self.filtername]
+
         # Compute how many observations we should have on the sky
-        desired = self.footprint / self.all_footprints_sum * np.sum(self.survey_features['N_obs_all'].feature)
+        desired = desired_footprint_normed * np.sum(self.survey_features['N_obs_all'].feature)
         result = desired - self.survey_features['N_obs'].feature
         result[self.out_of_bounds_area] = self.out_of_bounds_val
         return result
@@ -82,7 +80,7 @@ class Footprint_rolling_basis_function(Base_basis_function):
     """
 
     def __init__(self, filtername='r', nside=None, footprints=None, all_footprints_sum=None, all_rolling_sum=None, out_of_bounds_val=-10,
-                 season_modulo=2, season_length=365.25, max_season=None, day_offset=None):
+                 season_modulo=2, season_length=365.25, max_season=None, day_offset=None, window_size=6.):
         super(Footprint_rolling_basis_function, self).__init__(nside=nside, filtername=filtername)
 
         # OK, going to find the parts of the map that are the same everywhere, and compute the
@@ -103,6 +101,8 @@ def __init__(self, filtername='r', nside=None, footprints=None, all_footprints_s
         self.day_offset = day_offset
         self.footprints = footprints
 
+        self.max_day_offset = np.max(self.day_offset)
+
         self.all_footprints_sum = all_footprints_sum
         self.all_rolling_sum = all_rolling_sum
 
@@ -139,18 +139,40 @@ def _calc_value(self, conditions, indx=None):
         result = self.result.copy()
 
         # Compute what season it is at each pixel
-        seasons = utils.season_calc(conditions.night, offset=self.day_offset,
-                                    modulo=self.season_modulo, max_season=self.max_season,
-                                    season_length=self.season_length)
+        seasons = conditions.season(modulo=self.season_modulo,
+                                    max_season=self.max_season, season_length=self.season_length)
+
+        # Update the coverage of the sky so far
+        # If RA to sun is zero, we are at phase np.pi/2.
+        coverage_map_phase = (conditions.ra - conditions.sun_RA_start+np.pi/2.) % (2.*np.pi)
+        zero = utils.step_line(0., 1., 365.25, phase=coverage_map_phase*365.25/2/np.pi)
+        t_elapsed = conditions.mjd - conditions.mjd_start
+        norm_coverage_raw = utils.step_line(t_elapsed, 1., 365.25, phase=coverage_map_phase*365.25/2/np.pi)
+        norm_coverage_raw -= zero
+
+        norm_coverage = norm_coverage_raw/np.max(norm_coverage_raw)
+        norm_footprint = self.footprints[-1] * norm_coverage
 
         # Compute the constant parts of the footprint like before
-        desired = self.footprints[-1] / self.all_footprints_sum * np.sum(self.survey_features['N_obs_all_-1'].feature)
+        desired = norm_footprint / self.all_footprints_sum * np.sum(self.survey_features['N_obs_all_-1'].feature)
         result[self.constant_footprint_indx] = desired[self.constant_footprint_indx] - self.survey_features['N_obs_-1'].feature[self.constant_footprint_indx]
 
         # Now for the rolling sections
+        norm_coverage = 0
+        ns = np.floor((t_elapsed+self.max_day_offset)/365.25)
+        if np.max(ns) > -1:
+            ns_on = ns/self.season_modulo
+            ns_off = ns - ns_on
+            norm_coverage = norm_coverage_raw - ns_off
+            norm_coverage = norm_coverage/np.max(norm_coverage)
+
         for season in np.unique(seasons[self.rolling_footprint_indx]):
+            if season == -1:
+                nf = norm_footprint
+            else:
+                nf = self.footprints[season] * norm_coverage
             season_indx = np.where(seasons[self.rolling_footprint_indx] == season)[0]
-            desired = self.footprints[season][self.rolling_footprint_indx][season_indx] / self.all_rolling_sum * np.sum(self.survey_features['N_obs_all_%i' % season].feature[self.rolling_footprint_indx])
+            desired = nf[self.rolling_footprint_indx[season_indx]] / self.all_rolling_sum * np.sum(self.survey_features['N_obs_all_%i' % season].feature[self.rolling_footprint_indx])
             result[self.rolling_footprint_indx[season_indx]] = desired - self.survey_features['N_obs_%i' % season].feature[self.rolling_footprint_indx][season_indx]
 
         result[self.out_of_bounds_area] = self.out_of_bounds_val

python/lsst/sims/featureScheduler/features/conditions.py

@@ -1,7 +1,7 @@
 import numpy as np
 from lsst.sims.utils import _approx_RaDec2AltAz, Site, _hpid2RaDec, m5_flat_sed
 import healpy as hp
-from lsst.sims.featureScheduler.utils import set_default_nside, match_hp_resolution, approx_altaz2pa
+from lsst.sims.featureScheduler.utils import set_default_nside, match_hp_resolution, approx_altaz2pa, season_calc
 
 __all__ = ['Conditions']
 
@@ -16,7 +16,8 @@ class Conditions(object):
     Unless otherwise noted, all values are assumed to be valid at the time
     given by self.mjd
     """
-    def __init__(self, nside=None, site='LSST', exptime=30.):
+    def __init__(self, nside=None, site='LSST', exptime=30., mjd_start=59853.5, season_offset=None,
+                 sun_RA_start=None):
         """
         Parameters
         ----------
@@ -27,8 +28,13 @@ def __init__(self, nside=None, site='LSST', exptime=30.):
             observatory paramteres like latitude and longitude.
         expTime : float (30)
             The exposure time to assume when computing the 5-sigma limiting depth
+        mjd_start : float (59853.5)
+            The starting MJD of the survey.
+        season_offset : np.array
+            A HEALpix array that specifies the day offset when computing the season for each HEALpix.
+        sun_RA_start : float (None)
 
-        Atributes (Set on init)
+        Attributes (Set on init)
         -----------
         nside : int
             Healpix resolution. All maps are set to this reslution.
@@ -40,7 +46,7 @@ def __init__(self, nside=None, site='LSST', exptime=30.):
         dec : np.array
             A healpix array with the Dec of each healpixel center (radians). Automatically generated.
 
-        Atributes (to be set by user/telemetry stream)
+        Attributes (to be set by user/telemetry stream)
         -------------------------------------------
         mjd : float
             Modified Julian Date (days).
@@ -151,7 +157,10 @@ def __init__(self, nside=None, site='LSST', exptime=30.):
         self.nside = nside
         self.site = Site(site)
         self.exptime = exptime
+        self.mjd_start = mjd_start
         hpids = np.arange(hp.nside2npix(nside))
+        self.season_offset = season_offset
+        self.sun_RA_start = sun_RA_start
         # Generate an empty map so we can copy when we need a new map
         self.zeros_map = np.zeros(hp.nside2npix(nside), dtype=float)
         self.nan_map = np.zeros(hp.nside2npix(nside), dtype=float)
@@ -223,6 +232,13 @@ def __init__(self, nside=None, site='LSST', exptime=30.):
 
         self.targets_of_opportunity = None
 
+        self._season = None
+
+        self.season_modulo = None
+        self.season_max_season = None
+        self.season_length = 365.25
+        self.season_floor = True
+
     @property
     def lmst(self):
         return self._lmst
@@ -311,6 +327,7 @@ def mjd(self, value):
         self._HA = None
         self._lmst = None
         self._az_to_sun = None
+        self._season = None
 
     @property
     def skybrightness(self):
@@ -361,3 +378,21 @@ def az_to_sun(self):
         if self._az_to_sun is None:
             self.calc_az_to_sun()
         return self._az_to_sun
+
+    # XXX, there's probably an elegant decorator that could do this caching automatically
+    def season(self, modulo=None, max_season=None, season_length=365.25, floor=True):
+        if self.season_offset is not None:
+            kwargs_match = (modulo == self.season_modulo) & (max_season == self.season_max_season) & (season_length == self.season_length) & (floor == self.season_floor)
+            if ~kwargs_match:
+                self.season_modulo = modulo
+                self.season_max_season = max_season
+                self.season_length = season_length
+                self.season_floor = floor
+            if (self._season is None) | (~kwargs_match):
+                self._season = season_calc(self.night, offset=self.season_offset,
+                                           modulo=modulo, max_season=max_season,
+                                           season_length=season_length, floor=floor)
+        else:
+            self._season = None
+
+        return self._season

python/lsst/sims/featureScheduler/features/features.py

@@ -481,3 +481,5 @@ def add_observation(self, observation, indx=None):
         if observation['filter'][0] == self.filtername:
             # I think this is how to broadcast things properly.
             self.feature[indx, :] += np.histogram(observation.rotSkyPos, bins=self.bins)[0]
+
+

python/lsst/sims/featureScheduler/modelObservatory/model_observatory.py

@@ -9,7 +9,7 @@
 from lsst.sims.seeingModel import SeeingData, SeeingModel
 from lsst.sims.cloudModel import CloudData
 from lsst.sims.featureScheduler.features import Conditions
-from lsst.sims.featureScheduler.utils import set_default_nside, approx_altaz2pa
+from lsst.sims.featureScheduler.utils import set_default_nside, approx_altaz2pa, create_season_offset
 from lsst.ts.observatory.model import ObservatoryModel, Target
 from astropy.coordinates import EarthLocation
 from astropy.time import Time
@@ -305,9 +305,6 @@ def __init__(self, nside=None, mjd_start=59853.5, seed=42, quickTest=True,
 
         self.mjd_start = mjd_start
 
-        # Conditions object to update and return on request
-        self.conditions = Conditions(nside=self.nside)
-
         self.sim_ToO = sim_ToO
 
         # Create an astropy location
@@ -381,6 +378,14 @@ def __init__(self, nside=None, mjd_start=59853.5, seed=42, quickTest=True,
             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(self.mjd)
+        season_offset = create_season_offset(self.nside, sun_moon_info['sun_RA'])
+        self.sun_RA_start = sun_moon_info['sun_RA'] + 0
+        # Conditions object to update and return on request
+        self.conditions = Conditions(nside=self.nside, mjd_start=mjd_start,
+                                     season_offset=season_offset, sun_RA_start=self.sun_RA_start)
+
+
         self.obsID_counter = 0
 
     def get_info(self):

python/lsst/sims/featureScheduler/sim_runner.py

@@ -37,7 +37,7 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su
         filter_scheduler = simple_filter_sched()
 
     if mjd_start is None:
-        mjd = observatory.mjd
+        mjd = observatory.mjd + 0
         mjd_start = mjd + 0
     else:
         mjd = mjd_start + 0
@@ -79,7 +79,7 @@ def sim_runner(observatory, scheduler, filter_scheduler=None, mjd_start=None, su
                 # ugh, "swap_filter" means "unmount filter"
                 observatory.observatory.swap_filter(filtername)
 
-        mjd = observatory.mjd
+        mjd = observatory.mjd + 0
         if verbose:
             if (mjd-mjd_track) > step:
                 progress = float(mjd-mjd_start)/mjd_run*100