Merge pull request #54 from lsst/feature/SNObjectSED_regions

Feature/sn object sed regions : 

minor BugFix, Speedups and additional feature of SED rectification in sims_catUtils.supernovae.snobject.SNObject + new tests to check feature of SED rectification.

BugFix : A bitwise and changed to and.
Speedup: When the time of observation is outside the temporal range of the model, and uses the modelOutSideTemporalRange 'zero' where the SEDs are 0, speed catsimBandflux calculations by returning 0 rather than assigning zeros to the SED, and doing a bandflux calculation.

SED rectification: At regions where the SALT model is poorly sampled in the training set, the SED can go negative. This forces such negative values to 0, which is important for tools like phosim that will sample the SED to generate photons. Tests are added. Jenkins testbuild https://ci.lsst.codes/job/stack-os-matrix/11442/ passes.

Note: This is a duplicate of a previous mangled pull request #49. That PR has a few more commits and I will delete the branch later.

examples/testSN.py

@@ -11,23 +11,21 @@
     - Band Mag for extincted SED in r Band
 
 SNIaCatalog_tests:
-A Class containing tests to check crictical functionality for SNIaCatalog 
+A Class containing tests to check crictical functionality for SNIaCatalog
 """
 import os
 import sqlite3
 import numpy as np
-import matplotlib.pyplot as plt
 import unittest
 
 # Lsst Sims Dependencies
 import lsst.utils.tests as utilsTests
-from lsst.sims.photUtils import Bandpass
+from lsst.sims.photUtils.PhotometricParameters import PhotometricParameters
 from lsst.sims.photUtils import BandpassDict
 from lsst.sims.utils import ObservationMetaData
 from lsst.sims.utils import spatiallySample_obsmetadata as sample_obsmetadata
 from lsst.sims.catUtils.utils import ObservationMetaDataGenerator
 from lsst.sims.catalogs.generation.db import CatalogDBObject
-from lsst.sims.catalogs.measures.instance import InstanceCatalog
 import eups
 
 # Routines Being Tested
@@ -76,7 +74,7 @@ def setUp(self):
                               x0=1.796112e-06)
 
         self.SNCosmoModel = self.SN_extincted.equivalentSNCosmoModel()
-
+        self.rectify_photParams = PhotometricParameters()
         self.lsstBandPass = BandpassDict.loadTotalBandpassesFromFiles()
         self.SNCosmoBP = sncosmo.Bandpass(wave=self.lsstBandPass['r'].wavelen,
                                           trans=self.lsstBandPass['r'].sb,
@@ -94,7 +92,62 @@ def test_SNstatenotEmpty(self):
         myDict = self.SN_extincted.SNstate
         for key in myDict.keys():
             assert myDict[key] is not None
-
+
+
+    def test_attributeDefaults(self):
+        """
+        Check the defaults and the setter properties for rectifySED and
+        modelOutSideRange
+        """
+        snobj = SNObject(ra=30., dec=-60., source='salt2')
+        assert snobj.rectifySED == True
+        assert snobj.modelOutSideTemporalRange == 'zero'
+
+        snobj.rectifySED = False
+        assert snobj.rectifySED == False
+        assert snobj.modelOutSideTemporalRange == 'zero'
+
+    def test_raisingerror_forunimplementedmodelOutSideRange(self):
+        """
+        check that correct error is raised if the user tries to assign an
+        un-implemented model value to
+        `sims.catUtils.supernovae.SNObject.modelOutSideTemporalRange`
+        """
+        snobj = SNObject(ra=30., dec=-60., source='salt2')
+        assert snobj.modelOutSideTemporalRange == 'zero'
+        with self.assertRaises(ValueError) as context:
+            snobj.modelOutSideTemporalRange = 'False'
+        self.assertEqual('Model not implemented, defaulting to zero method\n',
+                         context.exception.message)
+
+    def test_rectifiedSED(self):
+        """
+        Check for an extreme case that the SN seds are being rectified. This is
+        done by setting up an extreme case where there will be negative seds, and
+        checking that this is indeed the case, and checking that they are not
+        negative if rectified.
+        """
+
+        snobj = SNObject(ra=30., dec=-60., source='salt2')
+        snobj.set(z=0.96, t0=self.mjdobs, x1=-3., x0=1.8e-6)
+        snobj.rectifySED  = False
+        times = np.arange(self.mjdobs - 50., self.mjdobs + 150., 1.)
+        badTimes = []
+        for time in times:
+            sed = snobj.SNObjectSED(time=time,
+                                    bandpass=self.lsstBandPass['r'])
+            if any(sed.flambda < 0.):
+                badTimes.append(time)
+        # Check that there are negative SEDs
+        assert(len(badTimes) > 0)
+        snobj.rectifySED = True
+        for time in badTimes:
+            sed = snobj.SNObjectSED(time=time,
+                                    bandpass=self.lsstBandPass['r'])
+            self.assertGreaterEqual(sed.calcADU(bandpass=self.lsstBandPass['r'],
+                                    photParams=self.rectify_photParams), 0.)
+            self.assertFalse(any(sed.flambda < 0.))
+
     def test_ComparebandFluxes2photUtils(self):
         """
         The SNObject.catsimBandFlux computation uses the sims.photUtils.sed

python/lsst/sims/catUtils/supernovae/snObject.py

@@ -57,6 +57,10 @@ class SNObject(sncosmo.Model):
         will be None, leading to exceptions in extinction calculation.
         Therefore, the value must be set explicitly to 0. to get unextincted
         quantities.
+    rectifySED : Bool, True by Default
+        if the SED is negative at the requested time and wavelength, return 0.
+        instead of the negative value.
+
 
     Methods
     -------
@@ -121,8 +125,19 @@ def __init__(self, ra=None, dec=None, source='salt2-extended'):
         self.ebvofMW = None
         if self._hascoords:
             self.mwEBVfromMaps()
+
+
+        # Behavior of model outside temporal range :
+        # if 'zero' then all fluxes outside the temporal range of the model
+        # are set to 0.
+        self._modelOutSideTemporalRange = 'zero'
+
+        # SED will be rectified to 0. for negative values of SED if this
+        # attribute is set to True
+        self.rectifySED = True
         return
 
+
     @property
     def SNstate(self):
         """
@@ -195,6 +210,21 @@ def fromSNState(cls, snState):
 
         return cls
 
+    @property
+    def modelOutSideTemporalRange(self):
+        """
+        Defines the behavior of the model when sampled at times beyond the model
+        definition.
+        """
+        return self._modelOutSideTemporalRange
+
+    @modelOutSideTemporalRange.setter
+    def modelOutSideTemporalRange(self, value):
+        if value != 'zero':
+            raise ValueError('Model not implemented, defaulting to zero method\n')
+        return self._modelOutSideTemporalRange
+
+
     def equivalentSNCosmoModel(self):
         """
         returns an SNCosmo Model which is equivalent to SNObject
@@ -210,6 +240,8 @@ def equivalentSNCosmoModel(self):
         sncosmoParams['mwebv'] = snState['MWE(B-V)']
         sncosmoModel.set(**sncosmoParams)
         return sncosmoModel
+
+
     @staticmethod
     def equivsncosmoParamDict(SNstate, SNCosmoModel):
         """
@@ -449,15 +481,13 @@ def SNObjectSED(self, time, wavelen=None, bandpass=None,
 
         flambda = np.zeros(len(wavelen))
 
-        # self.mintime() and self.maxtime() are properties describing
-        # the ranges of SNCosmo.Model in time.
-
-        # Set SED to 0 beyond the model phase range, will change this if
-        # SNCosmo includes a more sensible decay later.
-        if (time > self.mintime()) & (time < self.maxtime()):
 
-            # If SNCosmo is requested a SED value beyond the model range
-            # it will crash. Try to prevent that by returning np.nan for
+        # self.mintime() and self.maxtime() are properties describing
+        # the ranges of SNCosmo.Model in time. Behavior beyond this is 
+        # determined by self.modelOutSideTemporalRange
+        if (time > self.mintime()) and (time < self.maxtime()):
+            # If SNCosmo is requested a SED value beyond the wavelength range
+            # of model it will crash. Try to prevent that by returning np.nan for
             # such wavelengths. This will still not help band flux calculations
             # but helps us get past this stage.
 
@@ -475,7 +505,15 @@ def SNObjectSED(self, time, wavelen=None, bandpass=None,
 
             flambda[mask] = self.flux(time=time, wave=wave)
             flambda[mask] = flambda[mask] * 10.0
-
+        else:
+            # use prescription for modelOutSideTemporalRange
+            if self.modelOutSideTemporalRange != 'zero':
+                raise ValueError('Model not implemented, change to zero\n')
+                # Else Do nothing as flambda is already 0.
+                # This takes precedence over being outside wavelength range
+
+        if self.rectifySED:
+            flambda = np.where(flambda > 0., flambda, 0.)
         SEDfromSNcosmo = Sed(wavelen=wavelen, flambda=flambda)
 
         if not applyExtinction:
@@ -521,6 +559,9 @@ def catsimBandFlux(self, time, bandpassobject):
         .. note: If there is an unphysical value of sed in
         the wavelength range, it produces a flux of  `np.nan`
         """
+        # Speedup for cases outside temporal range of model
+        if time <= self.mintime() or time >= self.maxtime() :
+            return 0.
         SEDfromSNcosmo = self.SNObjectSED(time=time,
                                           bandpass=bandpassobject)
         return SEDfromSNcosmo.calcFlux(bandpass=bandpassobject) / 3631.0