BGS Archetypes

The goal of this notebook is to derive a set of spectral archetypes from the BGS template set using Guangtun Zhu's SetCoverPy algorithm.

Preliminaries.


In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt

In [2]:
from desispec.io.util import write_bintable, makepath
from desisim.io import write_templates
from desisim.archetypes import compute_chi2, ArcheTypes

In [3]:
import multiprocessing
nproc = multiprocessing.cpu_count() // 2

In [4]:
plt.style.use('seaborn-talk')
%matplotlib inline

Initialize the random seed so the results are reproducible, below.


In [5]:
seed = 123
rand = np.random.RandomState(seed)

Output path and filenames.


In [6]:
version = 'v1.0'

In [7]:
outdir = os.path.join(os.getenv('DESI_ROOT'), 'spectro', 'templates', 'archetypes', 'bgs', version)
print('Setting output directory to {}'.format(outdir))
os.makedirs(outdir, exist_ok=True)
chi2file = os.path.join(outdir, 'bgs_archetypes_chi2_{}.fits'.format(version))
archfile = os.path.join(outdir, 'bgs_archetypes_{}.fits'.format(version))


Setting output directory to /Users/ioannis/work/desi/spectro/templates/archetypes/bgs/v1.0

Read the BGS basis templates.

Read both a set of lower-resolution (1 A/pix) templates sampled over a restricted wavelength range (roughly 3500-7000 A) and the same set at higher resolution (0.2 A/pix) and over the wavelength range 0.12-2 micron. The lower-resolution templates will be used to determine the archetypes (since speed is an issue) while the full-resolution templates is what we actually write out.

In both cases we (arbitrarily) normalize every template to r=18 and adopt a nominal velocity dispersion of 100 km/s.


In [8]:
def _build_templates(args):
    """Filler function for the multiprocessing."""
    build_templates(*args)

In [9]:
def build_templates(bgs, input_meta, verbose=False):
    flux, _, meta = bgs.make_templates(input_meta=input_meta, novdisp=True, 
                                          nocolorcuts=True, verbose=verbose)
    return [flux.astype('f4'), meta]

In [10]:
def read_and_normalize(verbose=False, nproc=1, minwave=1200, maxwave=2e4, 
                       cdelt=0.2, nominal_rmag=18.0, nominal_vdisp=1000, 
                       subset=False):
    """Read and normalize the full set of basis templates.
    
    """
    from astropy.table import vstack
    from desisim.templates import BGS
    from desisim.io import empty_metatable

    bgs = BGS(minwave=minwave, maxwave=maxwave, cdelt=cdelt)
    bgs.normline = None # no emission line
    
    nspec = len(bgs.basemeta)
    if subset:
        nspec = 1000
        these = rand.choice(len(bgs.basemeta), nspec)
        print('Selecting a subset of {} / {} templates!'.format(nspec, len(bgs.basemeta)))
    else:
        these = np.arange(nspec)
    
    input_meta = empty_metatable(nmodel=nspec, objtype='BGS')
    input_meta['TEMPLATEID'] = these
    input_meta['REDSHIFT'] = 0.0
    input_meta['MAG'] = nominal_rmag
    input_meta['VDISP'] = nominal_vdisp
    input_meta['SEED'] = rand.randint(2**32, size=nspec)

    # Not sure why multiprocessing isn't working in this case.
    if nproc > 1:
        chunk = np.array_split(these, nproc)
    
        tempargs = list()
        for ii in range(nproc):
            tempargs.append((bgs, input_meta[chunk[ii]], verbose))
        
        pool = multiprocessing.Pool(nproc)
        out = pool.map(_build_templates, tempargs)
        flux = np.vstack(out[0])
        meta = vstack(out[1])
    else:
        flux, meta = build_templates(bgs, input_meta, verbose)

    nspec, npix = flux.shape
    print('Generated {} rest-frame BGS spectra with {} pixels.'.format(nspec, npix))
    
    return flux, bgs.wave, meta, bgs.basemeta[these]

In [11]:
%time flux, wave, meta, basemeta = read_and_normalize(nproc=1, cdelt=1.0, minwave=3600, maxwave=7000, subset=True)
nspec, npix = flux.shape


INFO:io.py:967:read_basis_templates: Reading /Users/ioannis/work/desi/spectro/templates/basis_templates/v2.3/bgs_templates_v2.1.fits
Selecting a subset of 1000 / 7636 templates!
Generated 1000 rest-frame BGS spectra with 3401 pixels.
CPU times: user 10.9 s, sys: 390 ms, total: 11.3 s
Wall time: 12.1 s

In [12]:
%time hiresflux, hireswave, _, _ = read_and_normalize(nproc=1, cdelt=0.2, minwave=1200, maxwave=2e4, subset=True)
_, hiresnpix = hiresflux.shape


INFO:io.py:967:read_basis_templates: Reading /Users/ioannis/work/desi/spectro/templates/basis_templates/v2.3/bgs_templates_v2.1.fits
Selecting a subset of 1000 / 7636 templates!
Generated 1000 rest-frame BGS spectra with 94001 pixels.
CPU times: user 41.9 s, sys: 4.51 s, total: 46.4 s
Wall time: 48.8 s

In [13]:
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 = rand.choice(nspec, nplot, replace=False)
    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, :])
        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)

In [14]:
plot_subset()


Compute the NxN chi2 matrix.

We use chi2 as the "distance" matrix for the Set Cover problem.

Then, we need to determine what threshold chi2 value differentiates "different" templates.

Note that the threshold chi^2 value can be tuned until the desired number of archetypes is achieved. However, If we want the archetypes to describe each spectrum in the parent sample to a precision of prec=0.1 (10%) then we we should set chi2min to be approximately npix*prec^2.


In [15]:
def write_chi2(chi2):
    from astropy.io import fits
    print('Writing {}'.format(chi2file))
    hdu = fits.PrimaryHDU(chi2)
    hdu.writeto(chi2file, overwrite=True)

In [16]:
%time chi2, amp = compute_chi2(flux)


INFO:archetypes.py:49:compute_chi2: Computing chi2 matrix for spectra 0-0 out of 1000.
INFO:archetypes.py:49:compute_chi2: Computing chi2 matrix for spectra 0-0 out of 1000.
CPU times: user 1min, sys: 17.7 s, total: 1min 17s
Wall time: 1min 19s

In [17]:
write_chi2(chi2)


Writing /Users/ioannis/work/desi/spectro/templates/archetypes/bgs/v1.0/bgs_archetypes_chi2_v1.0.fits

In [18]:
prec = 0.1
chi2min_nominal = npix*prec**2
print(chi2min_nominal, np.log10(chi2min_nominal)) # seems high...


34.010000000000005 1.53160663193

In [19]:
with np.errstate(divide='ignore'):
    logchi2 = np.log10(chi2)
logchi2[chi2 == 0] = -1
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
im = ax[0].imshow(logchi2, origin='lower', interpolation='nearest', 
                  vmin=-1.0, cmap='viridis')
ax[0].set_xlabel('Spectrum Number')
ax[0].set_ylabel('Spectrum Number')
plt.colorbar(im, label='$log_{10}(\chi^{2})$', ax=ax[0])

_ = ax[1].hist(logchi2.reshape(nspec * nspec), bins=30, range=(-1.2, np.max(logchi2)))
ax[1].set_ylabel('Number')
ax[1].set_xlabel('$log_{10}(\chi^{2}$)')


Out[19]:
<matplotlib.text.Text at 0x1130bd748>

Compute and plot the number of archetypes vs chi2 threshold.


In [20]:
def narch_vs_chi2min(Arch):
    """Determine the number of archtypes vs chi2 threshold.
    
    """
    cost = np.ones(nspec) # uniform cost
    chi2min = np.logspace(1, 5, 10)
    print(chi2min)
    narch = np.zeros_like(chi2min)
    
    for ii, cmin in enumerate(chi2min):
        iarch = Arch.get_archetypes(chi2_thresh=cmin)
        narch[ii] = len(iarch)
        
    return narch, chi2min

In [21]:
def qa_narch_vs_chi2min():
    fig, ax = plt.subplots()
    ax.scatter(np.log10(chi2min), narch)
    ax.set_xlabel('$log_{10}(\chi^{2})$ Threshold')
    ax.set_ylabel('Number of Archetypes')
    ax.axvline(x=np.log10(chi2min_nominal), color='red', ls='-')
    ax.grid(True)

In [22]:
Arch = ArcheTypes(chi2)

In [23]:
narch, chi2min = narch_vs_chi2min(Arch)


[  1.00000000e+01   2.78255940e+01   7.74263683e+01   2.15443469e+02
   5.99484250e+02   1.66810054e+03   4.64158883e+03   1.29154967e+04
   3.59381366e+04   1.00000000e+05]
This Best solution: UB=549.0, LB=548.6211221073318, UB1=549.0, LB1=548.6211221073318
Current Best Solution: UB=549.0, LB=548.6211221073318, change=0.06901236660623726% @ niters=0
Current Best Solution: UB=549.0, LB=548.6211221073318, change=0.06901236660623726% @ niters=1
Final Best solution: 549.0
Took 0.152 minutes to reach current solution.
This Best solution: UB=378.0, LB=-79.8381236225398, UB1=378.0, LB1=-79.8381236225398
Current Best Solution: UB=378.0, LB=-79.8381236225398, change=121.12119672553962% @ niters=0
Current Best Solution: UB=369.0, LB=368.8757069389536, change=0.03368375638113476% @ niters=4
Final Best solution: 369.0
Took 0.714 minutes to reach current solution.
This Best solution: UB=195.0, LB=-1920.8756516870749, UB1=195.0, LB1=-1920.8756516870749
Current Best Solution: UB=195.0, LB=-1920.8756516870749, change=1085.0644367626026% @ niters=0
This Best solution: UB=181.0, LB=178.9145390495732, UB1=182.0, LB1=178.9171497319192
Current Best Solution: UB=181.0, LB=178.9145390495732, change=1.1521883704015508% @ niters=5
This Best solution: UB=199.0, LB=-2066.622305356709, UB1=199.0, LB1=-2066.622305356709
Current Best Solution: UB=181.0, LB=178.6151919824015, change=1.3175734903859049% @ niters=10
This Best solution: UB=181.0, LB=179.33174335992732, UB1=182.0, LB1=179.35450824732777
Current Best Solution: UB=181.0, LB=179.33174335992732, change=0.9216887514213725% @ niters=15
Iteration in re-initialization reaches maximum number = 20
Current Best Solution: UB=181.0, LB=179.62648182456036, change=0.7588498206848852% @ niters=20
Final Best solution: 181.0
Took 4.038 minutes to reach current solution.
This Best solution: UB=78.0, LB=-5946.206911745324, UB1=78.0, LB1=-5946.206911745324
Current Best Solution: UB=78.0, LB=-5946.206911745324, change=7723.342194545287% @ niters=0
This Best solution: UB=73.0, LB=63.02899988431683, UB1=74.0, LB1=63.03623564924859
Current Best Solution: UB=73.0, LB=63.02899988431683, change=13.658904268059135% @ niters=5
This Best solution: UB=82.0, LB=-25308.86290810633, UB1=82.0, LB1=-25308.86290810633
Current Best Solution: UB=73.0, LB=63.02899988431683, change=13.658904268059135% @ niters=10
This Best solution: UB=72.0, LB=63.594207173933995, UB1=75.0, LB1=63.60909739899008
Current Best Solution: UB=72.0, LB=63.594207173933995, change=11.674712258425007% @ niters=15
Iteration in re-initialization reaches maximum number = 20
Current Best Solution: UB=72.0, LB=63.594207173933995, change=11.674712258425007% @ niters=20
Final Best solution: 72.0
Took 3.750 minutes to reach current solution.
This Best solution: UB=30.0, LB=-12337.295807791095, UB1=30.0, LB1=-12337.295807791095
Current Best Solution: UB=30.0, LB=-12337.295807791095, change=41224.31935930365% @ niters=0
This Best solution: UB=26.0, LB=20.501681939875308, UB1=27.0, LB1=20.501885095896117
Current Best Solution: UB=26.0, LB=20.501681939875308, change=21.14737715432574% @ niters=5
This Best solution: UB=31.0, LB=-64212.44939139554, UB1=31.0, LB1=-64212.44939139554
Current Best Solution: UB=26.0, LB=21.40154577554768, change=17.686362401739686% @ niters=10
This Best solution: UB=26.0, LB=21.041681913046915, UB1=28.0, LB1=21.066646098719687
Current Best Solution: UB=26.0, LB=21.041681913046915, change=19.07045418058879% @ niters=15
Iteration in re-initialization reaches maximum number = 20
Current Best Solution: UB=26.0, LB=21.041681913046915, change=19.07045418058879% @ niters=20
Final Best solution: 26.0
Took 4.259 minutes to reach current solution.
This Best solution: UB=11.0, LB=-19092.78154069383, UB1=11.0, LB1=-19092.78154069383
Current Best Solution: UB=11.0, LB=-19092.78154069383, change=173670.74127903485% @ niters=0
This Best solution: UB=11.0, LB=5.491339834476656, UB1=12.0, LB1=5.491387848997409
Current Best Solution: UB=10.0, LB=5.813155268675934, change=41.86844731324066% @ niters=5
This Best solution: UB=11.0, LB=-125888.82595050837, UB1=11.0, LB1=-125888.82595050837
Current Best Solution: UB=9.0, LB=5.7046339304176295, change=36.61517855091523% @ niters=10
This Best solution: UB=11.0, LB=5.754843611435234, UB1=12.0, LB1=5.755400723684248
Current Best Solution: UB=9.0, LB=5.7516454300213535, change=36.0928285553183% @ niters=15
Iteration in re-initialization reaches maximum number = 20
Current Best Solution: UB=9.0, LB=5.9015883837751835, change=34.42679573583129% @ niters=20
Final Best solution: 9.0
Took 5.476 minutes to reach current solution.
This Best solution: UB=6.0, LB=-27802.981698290594, UB1=6.0, LB1=-27802.981698290594
Current Best Solution: UB=6.0, LB=-27802.981698290594, change=463483.0283048433% @ niters=0
This Best solution: UB=5.0, LB=3.091012136349569, UB1=5.0, LB1=3.091012136349569
Current Best Solution: UB=5.0, LB=3.091012136349569, change=38.17975727300862% @ niters=5
This Best solution: UB=5.0, LB=-198563.36816751453, UB1=5.0, LB1=-198563.36816751453
Current Best Solution: UB=5.0, LB=3.123981255035511, change=37.52037489928978% @ niters=10
This Best solution: UB=7.0, LB=3.093196011464177, UB1=7.0, LB1=3.093196011464177
Current Best Solution: UB=5.0, LB=3.0982251650618444, change=38.03549669876311% @ niters=15
Iteration in re-initialization reaches maximum number = 20
Current Best Solution: UB=5.0, LB=3.0716195917041293, change=38.56760816591741% @ niters=20
Final Best solution: 5.0
Took 6.887 minutes to reach current solution.
This Best solution: UB=3.0, LB=-39471.58636443817, UB1=3.0, LB1=-39471.58636443817
Current Best Solution: UB=3.0, LB=-39471.58636443817, change=1315819.5454812723% @ niters=0
Current Best Solution: UB=2.0, LB=1.9994048018126314, change=0.02975990936843198% @ niters=2
Final Best solution: 2.0
Took 0.942 minutes to reach current solution.
This Best solution: UB=2.0, LB=-49218.913574261984, UB1=2.0, LB1=-49218.913574261984
Current Best Solution: UB=2.0, LB=-49218.913574261984, change=2461045.678713099% @ niters=0
This Best solution: UB=3.0, LB=1.0020002429014234, UB1=3.0, LB1=1.0020002429014234
Current Best Solution: UB=2.0, LB=1.001911896844928, change=49.904405157753594% @ niters=5
This Best solution: UB=2.0, LB=-437386.9532060084, UB1=2.0, LB1=-437386.9532060084
Current Best Solution: UB=2.0, LB=1.001911896844928, change=49.904405157753594% @ niters=10
This Best solution: UB=2.0, LB=0.8014834423995294, UB1=2.0, LB1=0.8014834423995294
Current Best Solution: UB=2.0, LB=0.9526893668837277, change=52.36553165581361% @ niters=15
Iteration in re-initialization reaches maximum number = 20
Current Best Solution: UB=2.0, LB=1.0020666721243243, change=49.89666639378378% @ niters=20
Final Best solution: 2.0
Took 3.189 minutes to reach current solution.
This Best solution: UB=1.0, LB=-52387.79684054009, UB1=1.0, LB1=-52387.79684054009
Current Best Solution: UB=1.0, LB=-52387.79684054009, change=5238879.684054009% @ niters=0
This Best solution: UB=1.0, LB=0.999251422271413, UB1=1.0, LB1=0.999251422271413
Current Best Solution: UB=1.0, LB=0.999251422271413, change=0.07485777285870521% @ niters=5
Current Best Solution: UB=1.0, LB=0.999251422271413, change=0.07485777285870521% @ niters=6
Final Best solution: 1.0
Took 0.849 minutes to reach current solution.

In [24]:
qa_narch_vs_chi2min()


Choose a chi2 threshold value then get the final set of archetypes.


In [25]:
def write_archetypes():
    """ToDo: Write out the responsibility indices for each archetype."""
    from astropy.table import Column
    outmeta = meta[iarch]
    outmeta.add_column(Column(name='RESPONSIBILITY', length=len(iarch), dtype='int8'))
    outmeta['RESPONSIBILITY'] = resp
    print('Writing {}'.format(archfile))
    write_templates(archfile, hiresflux[iarch, :], hireswave, outmeta, objtype='BGS Archetypes')

In [39]:
chi2_thresh = 10**2.5
print('Choosing a log10(chi2) threshold value of {:.1f}.'.format(np.log10(chi2_thresh)))


Choosing a log10(chi2) threshold value of 2.5.

In [40]:
_iarch, _resp, _respindx = Arch.get_archetypes(chi2_thresh=chi2_thresh, responsibility=True)
print('Generated {} archetypes.'.format(len(_iarch)))


This Best solution: UB=52.0, LB=-7982.826243784226, UB1=52.0, LB1=-7982.826243784226
Current Best Solution: UB=52.0, LB=-7982.826243784226, change=15451.588930354279% @ niters=0
This Best solution: UB=48.0, LB=38.59242051548476, UB1=52.0, LB1=38.79008669197965
Current Best Solution: UB=47.0, LB=40.747988950926164, change=13.302151168242204% @ niters=5
This Best solution: UB=52.0, LB=-29668.60852356831, UB1=52.0, LB1=-29668.60852356831
Current Best Solution: UB=47.0, LB=40.520325872377136, change=13.786540697069924% @ niters=10
This Best solution: UB=48.0, LB=41.23414839517082, UB1=51.0, LB1=41.24960219401543
Current Best Solution: UB=47.0, LB=41.83594869617272, change=10.987343199632507% @ niters=15
Iteration in re-initialization reaches maximum number = 20
Current Best Solution: UB=47.0, LB=42.21995243595384, change=10.170313966055653% @ niters=20
Final Best solution: 47.0
Took 3.754 minutes to reach current solution.
Generated 47 archetypes.

Sort by Dn(4000).


In [41]:
srt = np.argsort(meta['D4000'][_iarch])
iarch = _iarch[srt]
resp = _resp[srt]
respindx = []
for ss in srt:
    respindx.append(_respindx[ss])

In [42]:
write_archetypes()


Writing /Users/ioannis/work/desi/spectro/templates/archetypes/bgs/v1.0/bgs_archetypes_v1.0.fits
WARNING: The unit 'dex' could not be saved to FITS format [astropy.io.fits.convenience]

Generate some QAplots.


In [43]:
def _markers():
    d4000 = meta['D4000']
    size = 110 * (1+(resp - resp.min()) / resp.ptp())
    shade = (d4000[iarch] - d4000[iarch].min()) / d4000[iarch].ptp()
    col = plt.cm.coolwarm(shade)
    return size, col

In [44]:
def qa_responsibility():
    """Generate a color-color plot with the symbol size scaled by the responsibility.
    
    """
    rz = -2.5 * np.log10(meta['FLUX_R'] / meta['FLUX_Z'])
    d4000 = meta['D4000']

    size, col = _markers()
    
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
    ax1.scatter(rz[iarch], resp, c=col, marker='o', s=size, 
                edgecolor='k')
    ax1.set_xlabel('r - z')
    ax1.set_ylabel('Responsibility')
    ax1.grid(True)

    ax2.scatter(d4000[iarch], resp, c=col, marker='o', s=size, 
               edgecolor='k')
    ax2.set_xlabel('$D_{n}(4000)$')
    ax2.grid(True)

    fig.subplots_adjust(wspace=0.05)

In [45]:
qa_responsibility()



In [46]:
def qa_colorcolor():
    """Generate a color-color plot with the symbol size scaled by the responsibility.
    
    """
    gr = -2.5 * np.log10(meta['FLUX_G'] / meta['FLUX_R'])
    rz = -2.5 * np.log10(meta['FLUX_R'] / meta['FLUX_Z'])
    d4000 = meta['D4000']

    size, col = _markers()

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
    ax1.scatter(rz, gr, s=30, c='lightgray', edgecolor='k')
    ax1.scatter(rz[iarch], gr[iarch], c=col, marker='o', s=size, 
                alpha=1, edgecolor='k')
    ax1.set_xlabel(r'$r - z$')
    ax1.set_ylabel(r'$g - r$')
    ax1.grid(True)

    ax2.scatter(d4000, rz, s=30, c='lightgray', edgecolor='k')
    ax2.scatter(d4000[iarch], rz[iarch], c=col, marker='o', s=size, 
                alpha=1, edgecolor='k')
    ax2.set_xlabel('$D_{n}(4000)$')
    ax2.set_ylabel(r'$g - r$')
    ax2.grid(True)

    fig.subplots_adjust(wspace=0.3)

In [47]:
qa_colorcolor()



In [62]:
def qa_archetypes(ncol=5, nfilter=11):
    """Plot the archetypes and the spectra for which they're responsible."""
    from scipy.signal import medfilt
    
    _, col = _markers()    
    narch = len(iarch)
    nrow = np.ceil(narch / ncol).astype('int')
    fig, ax = plt.subplots(nrow, ncol, figsize=(2.5*ncol, 2.5*nrow), sharey=True, sharex=True)
    
    ww = (hireswave > 3000) * (hireswave < 1e4)
    for jj, (thisax, indx, rindx, rr) in enumerate(zip(ax.flat, iarch, respindx, resp)):
        if rr > 1:
            for ii in rindx:
                thisax.plot(hireswave[ww], hiresflux[ii, ww], color='lightgrey')
        smoothflux = medfilt(hiresflux[indx, ww], 11)
        thisax.plot(hireswave[ww], smoothflux, color=col[jj])
        thisax.xaxis.set_major_locator(plt.MaxNLocator(2))
        thisax.text(0.95, 0.93, '{:04d}\nResp={}'.format(indx, rr), ha='right', 
             va='top', transform=thisax.transAxes, fontsize=11)
    fig.subplots_adjust(wspace=0.05, hspace=0.05)

In [63]:
qa_archetypes()



In [50]:
def qa_ages_colormag():
    """Generate a color-magnitude plot for the original AGES sample.
    
    """
    Mr = basemeta['SDSS_UGRIZ_ABSMAG_Z01'][:, 2]
    gr = basemeta['SDSS_UGRIZ_ABSMAG_Z01'][:, 1]-basemeta['SDSS_UGRIZ_ABSMAG_Z01'][:, 2]
    size, col = _markers()

    fig, ax = plt.subplots(figsize=(6, 4))
    ax.scatter(Mr, gr, s=30, c='lightgray', edgecolor='k')
    ax.scatter(Mr[iarch], gr[iarch], c=col, marker='o', s=size, 
                alpha=1, edgecolor='k')
    ax.set_xlabel(r'$M_{0.1r}$')
    ax.set_ylabel(r'$^{0.1}(g - r)$')
    ax.set_xlim(-16, -23)
    ax.set_ylim(0, 1.3)
    
    ax.grid(True)

In [51]:
qa_ages_colormag()



In [ ]: