<jupyter_start><jupyter_text>BGS Reference SpectraThe goal of this notebook is to generate a set of ~1000 BGS-like spectra which can serve as a set of reference spectra for the BGS exposure-time model (see, e.g., https://github.com/desihub/desisurvey/issues/77).The basic methodology is to assign BGS templates to galaxies from the MXXL/BGS mock which lie in the green valley and therefore have an intermediate-strength 4000-A break and weak or no emission lines.  For these types of spectra we expect to have a somewhat lower redshift efficiency, and so they should yield conservative S/N estimates for the survey.  Preliminaries<jupyter_code>import os
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import desitarget.mock.io as mockio
import multiprocessing
nproc = multiprocessing.cpu_count() // 2 
%matplotlib inline
sns.set(style='ticks', font_scale=1.4, palette='Set2')
col = sns.color_palette()<jupyter_output><empty_output><jupyter_text>Establish the I/O path, output filename, and random seed.<jupyter_code>version = 'v1.0'
refdir = os.path.join(os.getenv('DESI_ROOT'), 'spectro', 'sim', 'bgs-refspec')
refspecfile = os.path.join(refdir, 'bgs-refspec-{}.fits'.format(version))
simseed = 123
simrand = np.random.RandomState(simseed)
nrefspec = 1000<jupyter_output><empty_output><jupyter_text>Read a handful of healpixels from the BGS/MXXL mock.This mock has the correct redshift distribution and correlation between redshift and apparent magnitude for BGS targets (among many other attributes).  Note that we set the velocity dispersion to a fiducial value of 100 km/s.<jupyter_code>mockfile = os.path.join(os.getenv('DESI_ROOT'), 'mocks', 'bgs', 'MXXL', 'desi_footprint', 'v0.0.4', 'BGS_r20.6.hdf5')
mockdata = mockio.read_durham_mxxl_hdf5(mockfile, rand=simrand, nside=32, nproc=nproc,
                                        healpixels=[3151,3150,3149,3148])
print(mockdata.keys())
mockdata['VDISP'] = np.repeat(100.0, len(mockdata['RA'])) # [km/s]
def qa_radec():
    fig, ax = plt.subplots()
    ax.scatter(mockdata['RA'], mockdata['DEC'], s=1)
    ax.set_xlabel('RA')
    ax.set_ylabel('Dec')    
qa_radec()
def qa_zmag():
    fig, ax = plt.subplots(1, 2, figsize=(12, 4))
    _ = ax[0].hist(mockdata['Z'], bins=100)
    ax[0].set_xlabel('Redshift')
    ax[0].set_ylabel('Number of Galaxies')
    ax[1].scatter(mockdata['Z'], mockdata['MAG'], s=1, alpha=0.75)
    ax[1].axhline(y=20.6, ls='--', color='k')
    ax[1].set_xlabel('Redshift')
    ax[1].set_ylabel(r'$r_{SDSS}$ (AB mag)')
qa_zmag()<jupyter_output><empty_output><jupyter_text>Select galaxies in the green valley.<jupyter_code>magaxis = np.arange(-23, -16, 0.1)
def green_valley(Mr, offset=0.0):
    return np.polyval([-0.04, 0.75], Mr+20) + offset
in_greenvalley = np.where(
    (mockdata['SDSS_01gr'] < green_valley(mockdata['SDSS_absmag_r01'], offset=+0.05)) * 
    (mockdata['SDSS_01gr'] > green_valley(mockdata['SDSS_absmag_r01'], offset=-0.05)))[0]    
print('Selected {} / {} green-valley galaxies.'.format(len(in_greenvalley), len(mockdata['RA'])))
greenvalley_mockdata = dict(Z=mockdata['Z'][in_greenvalley],
                            SEED=mockdata['SEED'][in_greenvalley],
                            MAG=mockdata['MAG'][in_greenvalley],
                            VDISP=mockdata['VDISP'][in_greenvalley],
                            SDSS_absmag_r01=mockdata['SDSS_absmag_r01'][in_greenvalley],
                            SDSS_01gr=mockdata['SDSS_01gr'][in_greenvalley])
def qa_colormag():
    fig, ax = plt.subplots()
    ax.scatter(mockdata['SDSS_absmag_r01'], mockdata['SDSS_01gr'], s=1)
    ax.scatter(mockdata['SDSS_absmag_r01'][in_greenvalley], 
               mockdata['SDSS_01gr'][in_greenvalley], s=1)
    ax.set_xlabel(r'$M_{0.1r}$')
    ax.set_ylabel(r'$^{0.1}(g - r)$')
    ax.set_xlim(-14, -24)
    ax.plot(magaxis, green_valley(magaxis, offset=0.0), color='k', ls='-')
    ax.plot(magaxis, green_valley(magaxis, offset=+0.05), color='k', ls='--')    
    ax.plot(magaxis, green_valley(magaxis, offset=-0.05), color='k', ls='--')        
qa_colormag()<jupyter_output><empty_output><jupyter_text>Use a KD tree to assign a basis template to each mock galaxy.To assign a (continuum) template we use the algorithm developed for the Data Challenge, which uses a simple KD tree constructed from the rest-frame color, magnitude, and redshift of each mock galaxy.  See ```desitarget.mock.spectra``` for more details.<jupyter_code>class BGStree(object):
    """Build a KD Tree."""
    def __init__(self):
        from speclite import filters
        from scipy.spatial import cKDTree as KDTree
        from desisim.io import read_basis_templates

        self.bgs_meta = read_basis_templates(objtype='BGS', onlymeta=True)
        self.bgs_tree = KDTree(self._bgs())

    def _bgs(self):
        """Quantities we care about: redshift (z), M_0.1r, and 0.1(g-r).

        """
        zobj = self.bgs_meta['Z'].data
        mabs = self.bgs_meta['SDSS_UGRIZ_ABSMAG_Z01'].data
        rmabs = mabs[:, 2]
        gr = mabs[:, 1] - mabs[:, 2]

        return np.vstack((zobj, rmabs, gr)).T        
    
    def query(self, objtype, matrix, subtype=''):
        """Return the nearest template number based on the KD Tree.

        Args:
          objtype (str): object type
          matrix (numpy.ndarray): (M,N) array (M=number of properties,
            N=number of objects) in the same format as the corresponding
            function for each object type (e.g., self.bgs).
          subtype (str, optional): subtype (only for white dwarfs)

        Returns:
          dist: distance to nearest template
          indx: index of nearest template
        
        """
        if objtype.upper() == 'BGS':
            dist, indx = self.bgs_tree.query(matrix)
        else:
            print('Unrecognized SUBTYPE {}!'.format(subtype))
            raise ValueError
                
        return dist, indx    
class BGStemplates(object):
    """Generate spectra.  

    """
    def __init__(self, wavemin=None, wavemax=None, dw=0.2,
                 rand=None, verbose=False):

        from desimodel.io import load_throughput

        self.tree = BGStree()

        # Build a default (buffered) wavelength vector.
        if wavemin is None:
            wavemin = load_throughput('b').wavemin - 10.0
        if wavemax is None:
            wavemax = load_throughput('z').wavemax + 10.0
            
        self.wavemin = wavemin
        self.wavemax = wavemax
        self.dw = dw
        self.wave = np.arange(round(wavemin, 1), wavemax, dw)

        self.rand = rand
        self.verbose = verbose

        # Initialize the templates once
        from desisim.templates import BGS
        self.bgs_templates = BGS(wave=self.wave, normfilter='sdss2010-r')
        self.bgs_templates.normline = None # no emission lines!
        
    def bgs(self, data, index=None, mockformat='durham_mxxl_hdf5'):
        """Generate spectra for BGS.

        Currently only the MXXL (durham_mxxl_hdf5) mock is supported.  DATA
        needs to have Z, SDSS_absmag_r01, SDSS_01gr, VDISP, and SEED, which are
        assigned in mock.io.read_durham_mxxl_hdf5.  See also BGSKDTree.bgs().

        """
        from desisim.io import empty_metatable

        objtype = 'BGS'
        if index is None:
            index = np.arange(len(data['Z']))
            
        input_meta = empty_metatable(nmodel=len(index), objtype=objtype)
        for inkey, datakey in zip(('SEED', 'MAG', 'REDSHIFT', 'VDISP'),
                                  ('SEED', 'MAG', 'Z', 'VDISP')):
            input_meta[inkey] = data[datakey][index]

        if mockformat.lower() == 'durham_mxxl_hdf5':
            alldata = np.vstack((data['Z'][index],
                                 data['SDSS_absmag_r01'][index],
                                 data['SDSS_01gr'][index])).T
            _, templateid = self.tree.query(objtype, alldata)
        else:
            raise ValueError('Unrecognized mockformat {}!'.format(mockformat))

        input_meta['TEMPLATEID'] = templateid
        flux, _, meta = self.bgs_templates.make_templates(input_meta=input_meta,
                                                          nocolorcuts=True, novdisp=False,
                                                          verbose=self.verbose)

        return flux, meta<jupyter_output><empty_output><jupyter_text>Generate the templates and write them out.<jupyter_code>def bgs_make_templates():
    """Generate the actual templates.  Iterate until we build the desired 
    number of models after target selection.
    
    """
    from astropy.table import vstack
    from desitarget.cuts import isBGS_bright, isBGS_faint
    
    natatime = np.min( (50, nrefspec) )
    ngood = 0
    flux, meta = [], []
    while ngood < nrefspec:
        these = BGSmaker.rand.choice(len(in_greenvalley), natatime)
        flux1, meta1 = BGSmaker.bgs(greenvalley_mockdata, index=these)
            
        keep = np.logical_or( isBGS_bright(rflux=meta1['FLUX_R'].data), 
                             isBGS_faint(rflux=meta1['FLUX_R'].data) )
        ngood1 = np.count_nonzero(keep)
        if ngood1 > 0:
            ngood += ngood1
            flux.append(flux1[keep, :])
            meta.append(meta1[keep])
            
    meta = vstack(meta)[:nrefspec]
    flux = np.vstack(flux)[:nrefspec, :]

    return flux, BGSmaker.wave, meta
BGSmaker = BGStemplates(wavemin=1500.0, wavemax=2e4, rand=simrand, verbose=False)
%time flux, wave, meta = bgs_make_templates()
def plot_subset(nplot=25, ncol=5):
    """Plot a random sampling of the basis templates."""
    nspec, npix = flux.shape
    nrow = np.ceil(nplot / ncol).astype('int')
    these = simrand.choice(nspec, nplot, replace=False)
    these = np.sort(these)
    ww = (wave > 5500) * (wave < 5550)

    fig, ax = plt.subplots(nrow, ncol, figsize=(2.2*ncol, 2.2*nrow), sharey=True, sharex=True)
    for thisax, indx in zip(ax.flat, these):
        thisax.plot(wave, flux[indx, :] / np.median(flux[indx, ww]))
        thisax.text(0.95, 0.93, '{:0d}'.format(indx), ha='right', 
             va='top', transform=thisax.transAxes, fontsize=11)
        thisax.xaxis.set_major_locator(plt.MaxNLocator(3))
    fig.subplots_adjust(wspace=0.05, hspace=0.05)
plot_subset()
meta
def write_templates(outfile):
    import astropy.units as u
    from astropy.io import fits

    hx = fits.HDUList()
    hdu_wave = fits.PrimaryHDU(wave)
    hdu_wave.header['EXTNAME'] = 'WAVE'
    hdu_wave.header['BUNIT'] = 'Angstrom'
    hdu_wave.header['AIRORVAC']  = ('vac', 'Vacuum wavelengths')
    hx.append(hdu_wave)    
    
    fluxunits = 1e-17 * u.erg / (u.s * u.cm**2 * u.Angstrom)
    hdu_flux = fits.ImageHDU(flux)
    hdu_flux.header['EXTNAME'] = 'FLUX'
    hdu_flux.header['BUNIT'] = str(fluxunits)
    hx.append(hdu_flux)
    
    hdu_meta = fits.table_to_hdu(meta)
    hdu_meta.header['EXTNAME'] = 'METADATA'
    hx.append(hdu_meta)

    print('Writing {}'.format(outfile))
    hx.writeto(outfile, clobber=True)
write_templates(refspecfile)<jupyter_output>WARNING: The unit 'dex' could not be saved to FITS format [astropy.io.fits.convenience]
{'fraction_non_alphanumeric': 0.1297709923664122, 'fraction_numerical': 0.02092501122586439, 'mean_word_length': 2.9183673469387754, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '110182', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3997, 'n_tokens_neox': 3797, 'n_words': 840, 'path': 'doc/nb/bgs-reference-spectra.ipynb', 'repo_name': 'HiramHerrera/desisim', 'size': '11135'}