Merge pull request #211 from desihub/targetspectra

miscellaneous tweaks, mostly to templates.py

.travis.yml

@@ -43,7 +43,8 @@ env:
         # - SCIPY_VERSION=0.16
         - ASTROPY_VERSION=1.2.1
         # - SPHINX_VERSION=1.3
-        - DESIUTIL_VERSION=1.8.0
+        # - DESIUTIL_VERSION=1.8.0
+        - DESIUTIL_VERSION=master
         - DESIMODEL_VERSION=0.5.0
         - DESIMODEL_DATA=branches/test-0.5.0
         - DESISIM_TESTDATA_VERSION=0.3.2

doc/changes.rst

@@ -9,6 +9,7 @@ desisim change log
 * pixsim add new required keywords DOSVER, FEEVER, DETECTOR
 * rewrote lya_spectra to achieve factor of 10 speedup.
 * COSMO (astropy.cosmology setup) is a new optional keyword for qso_desi_templates.
+* minor enhancements to templates code
 
 0.17.1 (2016-12-05)
 -------------------

py/desisim/io.py

@@ -559,7 +559,7 @@ 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):
+def read_basis_templates(objtype, 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.
@@ -572,6 +572,7 @@ def read_basis_templates(objtype, outwave=None, nspec=None, infile=None):
         infile (str, optional): full path to input template file to read,
             over-riding the contents of the $DESI_BASIS_TEMPLATES environment
             variable.
+        onlymeta (Bool, optional): read just the metadata table and return
 
     Returns:
         Tuple of (outflux, outwave, meta) where
@@ -598,6 +599,10 @@ def read_basis_templates(objtype, outwave=None, nspec=None, infile=None):
     if infile is None:
         infile = find_basis_template(ltype)
 
+    if onlymeta:
+        log.info('Reading {} metadata.'.format(infile))
+        return Table(fits.getdata(infile, 1))
+
     log.info('Reading {}'.format(infile))
 
     if objtype.upper() == 'QSO':
@@ -769,3 +774,18 @@ def empty_metatable(nmodel=1, objtype='ELG', add_SNeIa=None):
     meta['OBJTYPE'] = objtype.upper()
 
     return meta
+
+def empty_star_properties(nstar=1):
+    """Initialize a "star_properties" table for desisim.templates."""
+    from astropy.table import Table, Column
+
+    star_properties = Table()
+    star_properties.add_column(Column(name='REDSHIFT', length=nstar, dtype='f4'))
+    star_properties.add_column(Column(name='MAG', length=nstar, dtype='f4'))
+    star_properties.add_column(Column(name='TEFF', length=nstar, dtype='f4'))
+    star_properties.add_column(Column(name='LOGG', length=nstar, dtype='f4'))
+    star_properties.add_column(Column(name='FEH', length=nstar, dtype='f4'))
+    star_properties.add_column(Column(name='SEED', length=nstar, dtype='int64',
+                                      data=np.zeros(nstar)-1))
+
+    return star_properties

py/desisim/targets.py

@@ -14,7 +14,7 @@
 from astropy.table import Table, Column, hstack
 
 from desimodel.focalplane import FocalPlane
-from desisim.io import empty_metatable
+from desisim.io import empty_metatable, empty_star_properties
 import desimodel.io
 from desispec.log import get_logger
 log = get_logger()
@@ -483,3 +483,15 @@ def sample_nz(objtype, n):
     #- Sample that distribution
     x = np.random.uniform(0.0, 1.0, size=n)
     return np.interp(x, cdf, zhi)
+
+
+def _default_wave(wavemin=None, wavemax=None, dw=0.2):
+    '''Construct and return the default wavelength vector.'''
+
+    if wavemin is None:
+        wavemin = desimodel.io.load_throughput('b').wavemin
+    if wavemax is None:
+        wavemax = desimodel.io.load_throughput('z').wavemax
+    wave = np.arange(round(wavemin, 1), wavemax, dw)
+
+    return wave

py/desisim/templates.py

@@ -17,65 +17,6 @@
 
 LIGHT = 2.99792458E5  #- speed of light in km/s
 
-class GaussianMixtureModel(object):
-    """Read and sample from a pre-defined Gaussian mixture model.
-
-    """
-    def __init__(self, weights, means, covars, covtype):
-        self.weights = weights
-        self.means = means
-        self.covars = covars
-        self.covtype = covtype
-        self.n_components, self.n_dimensions = self.means.shape
-
-    @staticmethod
-    def save(model, filename):
-        from astropy.io import fits
-        hdus = fits.HDUList()
-        hdr = fits.Header()
-        hdr['covtype'] = model.covariance_type
-        hdus.append(fits.ImageHDU(model.weights_, name='weights', header=hdr))
-        hdus.append(fits.ImageHDU(model.means_, name='means'))
-        hdus.append(fits.ImageHDU(model.covars_, name='covars'))
-        hdus.writeto(filename, clobber=True)
-
-    @staticmethod
-    def load(filename):
-        from astropy.io import fits
-        hdus = fits.open(filename, memmap=False)
-        hdr = hdus[0].header
-        covtype = hdr['covtype']
-        model = GaussianMixtureModel(
-            hdus['weights'].data, hdus['means'].data, hdus['covars'].data, covtype)
-        hdus.close()
-        return model
-
-    def sample(self, n_samples=1, random_state=None):
-
-        if self.covtype != 'full':
-            return NotImplementedError(
-                'covariance type "{0}" not implemented yet.'.format(self.covtype))
-
-        # Code adapted from sklearn's GMM.sample()
-        if random_state is None:
-            random_state = np.random.RandomState()
-
-        weight_cdf = np.cumsum(self.weights)
-        X = np.empty((n_samples, self.n_dimensions))
-        rand = random_state.rand(n_samples)
-        # decide which component to use for each sample
-        comps = weight_cdf.searchsorted(rand)
-        # for each component, generate all needed samples
-        for comp in range(self.n_components):
-            # occurrences of current component in X
-            comp_in_X = (comp == comps)
-            # number of those occurrences
-            num_comp_in_X = comp_in_X.sum()
-            if num_comp_in_X > 0:
-                X[comp_in_X] = random_state.multivariate_normal(
-                    self.means[comp], self.covars[comp], num_comp_in_X)
-        return X
-
 class EMSpectrum(object):
     """Construct a complete nebular emission-line spectrum.
 
@@ -117,6 +58,7 @@ class EMSpectrum(object):
     def __init__(self, minwave=3650.0, maxwave=7075.0, cdelt_kms=20.0, log10wave=None):
         from pkg_resources import resource_filename
         from astropy.table import Table, Column, vstack
+        from desiutil.sklearn import GaussianMixtureModel
 
         # Build a wavelength array if one is not given.
         if log10wave is None:
@@ -445,7 +387,7 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
                               minlineflux=0.0, sne_rfluxratiorange=(0.01, 0.1),
                               seed=None, redshift=None, mag=None, vdisp=None,
                               input_meta=None, nocolorcuts=False, nocontinuum=False,
-                              agnlike=False, restframe=False):
+                              agnlike=False, novdisp=False, restframe=False):
         """Build Monte Carlo galaxy spectra/templates.
 
         This function chooses random subsets of the basis continuum spectra (for
@@ -516,6 +458,8 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
             the output spectrum (useful for testing; default False).  Note that
             this option automatically sets nocolorcuts to True and add_SNeIa to
             False.
+          novdisp (bool, optional): Do not velocity-blur the spectrum (default
+            False).
           agnlike (bool, optional): Adopt AGN-like emission-line ratios (e.g.,
             for the LRGs and some BGS galaxies) (default False, meaning we adopt
             star-formation-like line-ratios).  Option not yet supported.
@@ -553,7 +497,7 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
 
             vzero = np.where(vdisp <= 0)[0]
             if len(vzero) > 0:
-                log.fatal('Velocity disperion is zero or negative in {} spectra!').format(len(vzero))
+                log.fatal('Velocity dispersion is zero or negative!')
                 raise ValueError
 
             if self.add_SNeIa:
@@ -745,7 +689,7 @@ def make_galaxy_templates(self, nmodel=100, zrange=(0.6, 1.6), magrange=(21.0, 2
                     tempid = templateid[this]
 
                     thisemflux = emflux * normlineflux[templateid[this]]
-                    if nocontinuum:
+                    if nocontinuum or novdisp:
                         blurflux = restflux[this, :] * magnorm[this]
                     else:
                         blurflux = (self.vdispblur((restflux[this, :] - thisemflux), vdisp[ii]) + \
@@ -819,7 +763,7 @@ def make_templates(self, nmodel=100, zrange=(0.6, 1.6), rmagrange=(21.0, 23.4),
                        oiiihbrange=(-0.5, 0.2), logvdisp_meansig=(1.9, 0.15),
                        minoiiflux=0.0, sne_rfluxratiorange=(0.1, 1.0), redshift=None,
                        mag=None, vdisp=None, seed=None, input_meta=None, nocolorcuts=False,
-                       nocontinuum=False, agnlike=False, restframe=False):
+                       nocontinuum=False, agnlike=False, novdisp=False, restframe=False):
         """Build Monte Carlo ELG spectra/templates.
 
         See the GALAXY.make_galaxy_templates function for documentation on the
@@ -853,7 +797,8 @@ def make_templates(self, nmodel=100, zrange=(0.6, 1.6), rmagrange=(21.0, 23.4),
                                                          minlineflux=minoiiflux, redshift=redshift, vdisp=vdisp,
                                                          mag=mag, sne_rfluxratiorange=sne_rfluxratiorange,
                                                          seed=seed, input_meta=input_meta, nocolorcuts=nocolorcuts,
-                                                         nocontinuum=nocontinuum, agnlike=agnlike, restframe=restframe)
+                                                         nocontinuum=nocontinuum, agnlike=agnlike, novdisp=novdisp,
+                                                         restframe=restframe)
 
         return outflux, wave, meta
 
@@ -898,7 +843,7 @@ def make_templates(self, nmodel=100, zrange=(0.01, 0.4), rmagrange=(15.0, 19.5),
                        oiiihbrange=(-1.3, 0.6), logvdisp_meansig=(2.0, 0.17),
                        minhbetaflux=0.0, sne_rfluxratiorange=(0.1, 1.0), redshift=None,
                        mag=None, vdisp=None, seed=None, input_meta=None, nocolorcuts=False,
-                       nocontinuum=False, agnlike=False, restframe=False):
+                       nocontinuum=False, agnlike=False, novdisp=False, restframe=False):
         """Build Monte Carlo BGS spectra/templates.
 
          See the GALAXY.make_galaxy_templates function for documentation on the
@@ -932,8 +877,9 @@ def make_templates(self, nmodel=100, zrange=(0.01, 0.4), rmagrange=(15.0, 19.5),
                                                          minlineflux=minhbetaflux, redshift=redshift, vdisp=vdisp,
                                                          mag=mag, sne_rfluxratiorange=sne_rfluxratiorange,
                                                          seed=seed, input_meta=input_meta, nocolorcuts=nocolorcuts,
-                                                         nocontinuum=nocontinuum, agnlike=agnlike, restframe=restframe)
-
+                                                         nocontinuum=nocontinuum, agnlike=agnlike, novdisp=novdisp,
+                                                         restframe=restframe)
+
         return outflux, wave, meta
 
 class LRG(GALAXY):
@@ -968,7 +914,7 @@ def __init__(self, minwave=3600.0, maxwave=10000.0, cdelt=0.2, wave=None,
     def make_templates(self, nmodel=100, zrange=(0.5, 1.0), zmagrange=(19.0, 20.5),
                        logvdisp_meansig=(2.3, 0.1), sne_rfluxratiorange=(0.1, 1.0),
                        redshift=None, mag=None, vdisp=None, seed=None,
-                       input_meta=None, nocolorcuts=False, agnlike=False,
+                       input_meta=None, nocolorcuts=False, novdisp=False, agnlike=False,
                        restframe=False):
         """Build Monte Carlo BGS spectra/templates.
 
@@ -999,7 +945,7 @@ def make_templates(self, nmodel=100, zrange=(0.5, 1.0), zmagrange=(19.0, 20.5),
                                                          logvdisp_meansig=logvdisp_meansig, redshift=redshift,
                                                          vdisp=vdisp, mag=mag, sne_rfluxratiorange=sne_rfluxratiorange,
                                                          seed=seed, input_meta=input_meta, nocolorcuts=nocolorcuts,
-                                                         agnlike=agnlike, restframe=restframe)
+                                                         novdisp=novdisp, agnlike=agnlike, restframe=restframe)
 
         # Pack into the metadata table some additional information.
         good = np.where(meta['TEMPLATEID'] != -1)[0]

py/desisim/test/test_templates.py

@@ -168,10 +168,11 @@ def test_star_properties(self):
         for T in [STAR]:
             Tx = T(wave=self.wave)
             flux, wave, meta = Tx.make_templates(star_properties=star_properties, seed=self.seed)
+            #import pdb ; pdb.set_trace()
             badkeys = list()
             for key in meta.colnames:
                 if key in star_properties.colnames:
-                    if not np.all(meta[key] == star_properties[key]):
+                    if not np.allclose(meta[key], star_properties[key]):
                         badkeys.append(key)
             self.assertEqual(len(badkeys), 0, 'mismatch for spectral type {} in keys {}'.format(meta['OBJTYPE'][0], badkeys))