Merge pull request #213 from desihub/wdsubtypes

support DA and DB white dwarf subtypes

.travis.yml

@@ -44,10 +44,10 @@ env:
         - ASTROPY_VERSION=1.2.1
         # - SPHINX_VERSION=1.3
         # - DESIUTIL_VERSION=1.8.0
-        - DESIUTIL_VERSION=master
+        - DESIUTIL_VERSION=1.9.3
         - DESIMODEL_VERSION=0.5.0
         - DESIMODEL_DATA=branches/test-0.5.0
-        - DESISIM_TESTDATA_VERSION=0.3.2
+        - DESISIM_TESTDATA_VERSION=0.3.3
         - SPECLITE_VERSION=0.5
         - SPECSIM_VERSION=v0.6
         - SPECTER_VERSION=0.6.0

py/desisim/io.py

@@ -9,6 +9,7 @@
 from glob import glob
 
 from astropy.io import fits
+from astropy.table import Table
 from astropy.stats import sigma_clipped_stats
 import numpy as np
 import multiprocessing
@@ -559,13 +560,19 @@ def _qso_format_version(filename):
         else:
             raise IOError('Unknown QSO basis template format '+filename)
 
-def read_basis_templates(objtype, outwave=None, nspec=None, infile=None, onlymeta=False):
+def read_basis_templates(objtype, subtype='', outwave=None, nspec=None,
+                         infile=None, onlymeta=False):
     """Return the basis (continuum) templates for a given object type.  Optionally
     returns a randomly selected subset of nspec spectra sampled at
     wavelengths outwave.
 
     Args:
-        objtype (str): object type to read (e.g., ELG, LRG, QSO, STAR, FSTD, WD, MWS_STAR, BGS).
+    
+        objtype (str): object type to read (e.g., ELG, LRG, QSO, STAR, FSTD, WD,
+          MWS_STAR, BGS).
+        subtype (str, optional): template subtype, currently only for white
+            dwarfs.  The choices are DA and DB and the default is to read both
+            types.
         outwave (numpy.array, optional): array of wavelength at which to sample
             the spectra.
         nspec (int, optional): number of templates to return
@@ -585,11 +592,8 @@ def read_basis_templates(objtype, outwave=None, nspec=None, infile=None, onlymet
         EnvironmentError: If the required $DESI_BASIS_TEMPLATES environment
             variable is not set.
         IOError: If the basis template file is not found.
-    """
-    from glob import glob
-    from astropy.io import fits
-    from astropy.table import Table
 
+    """
     ltype = objtype.lower()
     if objtype == 'FSTD':
         ltype = 'star'
@@ -601,7 +605,16 @@ def read_basis_templates(objtype, outwave=None, nspec=None, infile=None, onlymet
 
     if onlymeta:
         log.info('Reading {} metadata.'.format(infile))
-        return Table(fits.getdata(infile, 1))
+        meta = Table(fits.getdata(infile, 1))
+
+        if (objtype.upper() == 'WD') and (subtype != ''):
+            keep = np.where(meta['WDTYPE'] == subtype.upper())[0]
+            if len(keep) == 0:
+                log.warning('Unrecognized white dwarf subtype {}!'.format(subtype))
+            else:
+                meta = meta[keep]
+
+        return meta
 
     log.info('Reading {}'.format(infile))
 
@@ -627,6 +640,14 @@ def read_basis_templates(objtype, outwave=None, nspec=None, infile=None, onlymet
         meta = Table(fits.getdata(infile, 1))
         wave = fits.getdata(infile, 2)
 
+        if (objtype.upper() == 'WD') and (subtype != ''):
+            keep = np.where(meta['WDTYPE'] == subtype.upper())[0]
+            if len(keep) == 0:
+                log.warning('Unrecognized white dwarf subtype {}!'.format(subtype))
+            else:
+                meta = meta[keep]
+                flux = flux[keep, :]
+
     # Optionally choose a random subset of spectra. There must be a fast way to
     # do this using fitsio.
     ntemplates = flux.shape[0]
@@ -712,12 +733,13 @@ def _parse_filename(filename):
     elif len(x) == 3:
         return x[0], x[1].lower(), int(x[2])
 
-def empty_metatable(nmodel=1, objtype='ELG', add_SNeIa=None):
+def empty_metatable(nmodel=1, objtype='ELG', subtype='', add_SNeIa=None):
     """Initialize the metadata table for each object type."""
     from astropy.table import Table, Column
 
     meta = Table()
     meta.add_column(Column(name='OBJTYPE', length=nmodel, dtype=(str, 10)))
+    meta.add_column(Column(name='SUBTYPE', length=nmodel, dtype=(str, 10)))
     meta.add_column(Column(name='TEMPLATEID', length=nmodel, dtype='i4',
                            data=np.zeros(nmodel)-1))
     meta.add_column(Column(name='SEED', length=nmodel, dtype='int64',
@@ -772,6 +794,7 @@ def empty_metatable(nmodel=1, objtype='ELG', add_SNeIa=None):
                                data=np.zeros(nmodel)-1, unit='days'))
 
     meta['OBJTYPE'] = objtype.upper()
+    meta['SUBTYPE'] = subtype.upper()
 
     return meta
 

py/desisim/lya_spectra.py

@@ -53,7 +53,6 @@ def get_spectra(lyafile, nqso=None, wave=None, templateid=None, normfilter='sdss
     else:
         templateid = np.array(templateid)
         nqso = len(templateid)
-    print(templateid)
 
     if rand is None:
         rand = np.random.RandomState(seed)
@@ -62,7 +61,6 @@ def get_spectra(lyafile, nqso=None, wave=None, templateid=None, normfilter='sdss
     #heads = [head.read_header() for head in h[templateid + 1]]
     heads = []
     for indx in templateid:
-        print(indx + 1)
         heads.append(h[indx + 1].read_header())
 
     zqso = np.array([head['ZQSO'] for head in heads])

py/desisim/templates.py

@@ -509,9 +509,9 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
             nmodel = len(input_meta)
             alltemplateid_chunk = [input_meta['TEMPLATEID'].data.reshape(nmodel, 1)]
 
-            meta = empty_metatable(nmodel, self.objtype, self.add_SNeIa)
+            meta = empty_metatable(nmodel=nmodel, objtype=self.objtype, add_SNeIa=self.add_SNeIa)
         else:
-            meta = empty_metatable(nmodel, self.objtype, self.add_SNeIa)
+            meta = empty_metatable(nmodel=nmodel, objtype=self.objtype, add_SNeIa=self.add_SNeIa)
 
             # Initialize the random seed.
             rand = np.random.RandomState(seed)
@@ -958,7 +958,7 @@ def make_templates(self, nmodel=100, zrange=(0.5, 1.0), zmagrange=(19.0, 20.5),
 class SUPERSTAR(object):
     """Base class for generating Monte Carlo spectra of the various flavors of stars."""
 
-    def __init__(self, objtype='STAR', minwave=3600.0, maxwave=10000.0, cdelt=0.2,
+    def __init__(self, objtype='STAR', subtype='', minwave=3600.0, maxwave=10000.0, cdelt=0.2,
                  wave=None, colorcuts_function=None, normfilter='decam2014-r',
                  baseflux=None, basewave=None, basemeta=None):
         """Read the appropriate basis continuum templates, filter profiles and
@@ -968,7 +968,9 @@ def __init__(self, objtype='STAR', minwave=3600.0, maxwave=10000.0, cdelt=0.2,
           Only a linearly-spaced output wavelength array is currently supported.
 
         Args:
-          objtype (str): type of object to simulate (default STAR)
+          objtype (str): type of object to simulate (default STAR).
+          subtype (str, optional): stellar subtype, currently only for white
+            dwarfs.  The choices are DA and DB and the default is DA.
           minwave (float, optional): minimum value of the output wavelength
             array (default 3600 Angstrom).
           maxwave (float, optional): minimum value of the output wavelength
@@ -997,6 +999,7 @@ def __init__(self, objtype='STAR', minwave=3600.0, maxwave=10000.0, cdelt=0.2,
         from speclite import filters
 
         self.objtype = objtype.upper()
+        self.subtype = subtype.upper()
         self.colorcuts_function = colorcuts_function
         self.normfilter = normfilter
 
@@ -1010,7 +1013,8 @@ def __init__(self, objtype='STAR', minwave=3600.0, maxwave=10000.0, cdelt=0.2,
         # Read the rest-frame continuum basis spectra, if not specified.
         if baseflux is None or basewave is None or basemeta is None:
             from desisim.io import read_basis_templates
-            baseflux, basewave, basemeta = read_basis_templates(objtype=self.objtype)
+            baseflux, basewave, basemeta = read_basis_templates(objtype=self.objtype,
+                                                                subtype=self.subtype)
         self.baseflux = baseflux
         self.basewave = basewave
         self.basemeta = basemeta
@@ -1105,7 +1109,7 @@ def make_star_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
         # Optionally unpack a metadata table.
         if input_meta is not None:
             nmodel = len(input_meta)
-            meta = empty_metatable(nmodel, self.objtype)
+            meta = empty_metatable(nmodel=nmodel, objtype=self.objtype, subtype=self.subtype)
 
             templateseed = input_meta['SEED'].data
             redshift = input_meta['REDSHIFT'].data
@@ -1166,7 +1170,7 @@ def make_star_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0),
                     mag = rand.uniform(magrange[0], magrange[1], nmodel).astype('f4')
 
             # Initialize the metadata table.
-            meta = empty_metatable(nmodel, self.objtype)
+            meta = empty_metatable(nmodel=nmodel, objtype=self.objtype, subtype=self.subtype)
 
         # Basic error checking and some preliminaries.
         if redshift is not None:
@@ -1449,7 +1453,7 @@ class WD(SUPERSTAR):
     """Generate Monte Carlo spectra of white dwarfs."""
 
     def __init__(self, minwave=3600.0, maxwave=10000.0, cdelt=0.2, wave=None,
-                 normfilter='decam2014-g', colorcuts_function=None,
+                 subtype='DA', normfilter='decam2014-g', colorcuts_function=None,
                  baseflux=None, basewave=None, basemeta=None):
         """Initialize the WD class.  See the SUPERSTAR.__init__ method for documentation
         on the arguments plus the inherited attributes.
@@ -1466,7 +1470,7 @@ def __init__(self, minwave=3600.0, maxwave=10000.0, cdelt=0.2, wave=None,
 
         """
 
-        super(WD, self).__init__(objtype='WD', minwave=minwave, maxwave=maxwave,
+        super(WD, self).__init__(objtype='WD', subtype=subtype, minwave=minwave, maxwave=maxwave,
                                  cdelt=cdelt, wave=wave, colorcuts_function=colorcuts_function,
                                  normfilter=normfilter, baseflux=baseflux, basewave=basewave,
                                  basemeta=basemeta)
@@ -1491,8 +1495,17 @@ def make_templates(self, nmodel=100, vrad_meansig=(0.0, 200.0), gmagrange=(16.0,
           meta (astropy.Table): Table of meta-data [nmodel] for each output spectrum.
 
         Raises:
+          ValueError: If the INPUT_META or STAR_PROPERTIES table contains
+            different values of SUBTYPE.
 
         """
+        for intable in (input_meta, star_properties):
+            if intable is not None:
+                if 'SUBTYPE' in intable.dtype.names:
+                    if (self.subtype != '') and ~np.all(intable['SUBTYPE'] == self.subtype):
+                        log.warning('WD Class initialized with subtype {}, which does not match input table.'.format(self.subtype))
+                        raise ValueError
+
         outflux, wave, meta = self.make_star_templates(nmodel=nmodel, vrad_meansig=vrad_meansig,
                                                        magrange=gmagrange, seed=seed, redshift=redshift,
                                                        mag=mag, input_meta=input_meta,
@@ -1649,10 +1662,10 @@ def make_templates(self, nmodel=100, zrange=(0.5, 4.0), rmagrange=(21.0, 23.0),
             redshift = input_meta['REDSHIFT'].data
             mag = input_meta['MAG'].data
 
-            meta = empty_metatable(nmodel, self.objtype)
+            meta = empty_metatable(nmodel=nmodel, objtype=self.objtype)
 
         else:
-            meta = empty_metatable(nmodel, self.objtype)
+            meta = empty_metatable(nmodel=nmodel, objtype=self.objtype)
 
             # Initialize the random seed.
             rand = np.random.RandomState(seed)

py/desisim/test/test_templates.py

@@ -114,19 +114,28 @@ def test_input_redshift(self):
             self.assertTrue(np.allclose(redshift, meta['REDSHIFT']))
 
     @unittest.skipUnless(desi_basis_templates_available, '$DESI_BASIS_TEMPLATES was not detected.')
-    def test_sne(self):
-        '''Test options for adding in SNeIa spectra'''
-        print('In function test_sne, seed = {}'.format(self.seed))
-        for T in [BGS]:
-            template_factory = T(wave=self.wave, add_SNeIa=True)
-            flux, wave, meta = template_factory.make_templates(self.nspec, seed=self.seed,
-                                                               nocolorcuts=True, 
-                                                               sne_rfluxratiorange=(0.5, 0.7))
-            self._check_output_size(flux, wave, meta)
-            self.assertTrue('SNE_TEMPLATEID' in meta.dtype.names)
-            self.assertTrue('SNE_RFLUXRATIO' in meta.dtype.names)
-            self.assertTrue('SNE_EPOCH' in meta.dtype.names)
-
+    def test_wd_subtype(self):
+        '''Test option of specifying the white dwarf subtype.'''
+        print('In function test_wd_subtype, seed = {}'.format(self.seed))
+        wd = WD(wave=self.wave, subtype='DA')
+        flux, wave, meta = wd.make_templates(self.nspec, seed=self.seed, nocolorcuts=True)
+        self._check_output_size(flux, wave, meta)
+        np.all(meta['SUBTYPE'] == 'DA')
+
+        wd = WD(wave=self.wave, subtype='DB')
+        flux, wave, meta = wd.make_templates(self.nspec, seed=self.seed, nocolorcuts=True)
+        np.all(meta['SUBTYPE'] == 'DB')
+
+    @unittest.skipUnless(desi_basis_templates_available, '$DESI_BASIS_TEMPLATES was not detected.')
+    @unittest.expectedFailure
+    def test_wd_subtype_failure(self):
+        '''Test a known failure of specifying the white dwarf subtype.'''
+        print('In function test_wd_subtype_failure, seed = {}'.format(self.seed))
+        wd = WD(wave=self.wave, subtype='DA')
+        flux1, wave1, meta1 = wd.make_templates(self.nspec, seed=self.seed, nocolorcuts=True)
+        meta1['SUBTYPE'][0] = 'DB'
+        flux2, wave2, meta2 = wd.make_templates(input_meta=meta1)
+
     @unittest.skipUnless(desi_basis_templates_available, '$DESI_BASIS_TEMPLATES was not detected.')
     def test_input_meta(self):
         '''Test that input meta table option works.'''