In [79]:
try:
    reload(apogee.modelatm.atlas9)
    reload(apogee.modelspec.moog)
    reload(apogee.modelspec.ferre)
except NameError:
    import apogee.modelatm.atlas9
    import apogee.modelspec.moog
    import apogee.modelspec.ferre
import apogee.spec.cannon as apcannon
import apogee.spec.plot as splot
import apogee.spec.window as apwindow
from apogee.tools import paramIndx, elemIndx
from matplotlib import cm
from galpy.util import bovy_plot

A closer look at the ASPCAP analysis of simulated MOOG spectra with just [O/H] variations

Convenience function to generate a model spectrum


In [2]:
def simspec(ofe=0.,teff=4000.,logg=2.5,metals=0.,cfe=0.,nfe=0.,isotopes='solar',solarscaledatm=False):
    # Set up a model atmosphere for this point
    if solarscaledatm:
        atm= apogee.modelatm.atlas9.Atlas9Atmosphere(teff=teff,logg=logg,metals=metals,am=0.,cm=0.)
        abu= ([8,ofe],)
    else:
        atm= apogee.modelatm.atlas9.Atlas9Atmosphere(teff=teff,logg=logg,metals=metals,am=ofe,cm=cfe)
        # Keep the other alpha elements the same
        abu= ([8,0.],[12,-ofe],[14,-ofe],[16,-ofe],[20,-ofe],[22,-ofe])
    synspec= apogee.modelspec.moog.synth(*abu,
                                         modelatm=atm,linelist='moog.201405291900.vac',lsf='combo',
                                         cont='aspcap',vmacro=6.,isotopes=isotopes,
                                         vmicro=2.478-0.325*logg)
    return synspec

Convenience function to analyze a simulated spectrum using the MOOG library


In [3]:
def fit(synspec,initcannon=True,blockOwindows=False):
    # initcannon=True only runs a single FERRE fit initialized close to the best fit
    # blockOwindows=True removes the O windows from the fparam fit
    # Assume SNR=200
    errs= 0.005*numpy.ones_like(synspec)
    if blockOwindows:
        errs[apwindow.tophat('O',apStarWavegrid=True)]= 10.
    fparam= apogee.modelspec.ferre.fit(synspec,errs,initcannon=initcannon,lib='msGK',dr='current')
    felem= apogee.modelspec.ferre.elemfitall(synspec,0.005*numpy.ones_like(synspec),fparam=fparam,lib='msGK',dr='current')
    return (fparam,felem)

In [80]:
def cannon(synspec,blockOwindows=False):
    # blockOwindows=True removes the O windows from the analysis
    # Assume SNR=200
    errs= 0.005*numpy.ones_like(synspec)
    if blockOwindows:
        errs[apwindow.tophat('O',apStarWavegrid=True)]= 10.
    labels= apcannon.polylabels(synspec,errs)
    return labels

A function to plot the results


In [100]:
def plot_results(var,fparam,felem,fparam_ref,labels,var_label=r'$\mathrm{O}$'):
    figsize(18,20) 
    # First do the parameters
    params= ['TEFF','LOGG','METALS','C','N','ALPHA']
    param_labels= [r'$\Delta T_{\mathrm{eff}}\,(\mathrm{K})$',r'$\Delta \log g$',r'$\Delta [\mathrm{M/H}]$',
                   r'$\Delta [\mathrm{C/M}]$',r'$\Delta [\mathrm{N/M}$]',r'$\Delta [\alpha/\mathrm{M}$]']
    cannonlabelIndx= [0,1,2,None,None,3]
    ranges= [[-300.,300.],[-1.,1.],[-0.4,0.4],[-0.4,0.4],[-0.4,0.4],[-0.4,0.4]]
    for ii, param in enumerate(params):
        subplot(7,3,ii+1)
        ty= fparam[:,paramIndx(param)]-fparam_ref[paramIndx(param)]
        tdy= numpy.amax(ty)-numpy.amin(ty)
        yrange= ranges[ii]
        bovy_plot.bovy_plot(var,ty,'ko',xlabel=var_label,ylabel=param_labels[ii],
                            xrange=[var[0]-0.1,var[-1]+0.1],yrange=yrange,gcf=True)
        if not cannonlabelIndx[ii] is None:
            ty= labels[:,cannonlabelIndx[ii]]-fparam_ref[paramIndx(param)]
            bovy_plot.bovy_plot(var,ty,'ro',overplot=True)
    # Now do the elements
    elems= ['C','N','O','Na','Mg','Al','Si','S','K','Ca','Ti','V','Mn','Fe','Ni']
    for ii, elem in enumerate(elems):
        subplot(7,3,ii+7)
        ty= felem[elem]
        if not elem.lower() == 'o':
            ylabel= r'$\Delta \mathrm{%s}$' % elem
        else:
            ylabel= r'$\mathrm{Recovered\ %s}$' % elem
        bovy_plot.bovy_plot(var,ty,'ko',xlabel=var_label,ylabel=ylabel,
                            xrange=[var[0]-0.1,var[-1]+0.1],yrange=[-0.4,0.4],gcf=True)
    #tight_layout()

$T_\mathrm{eff}$ = 4000 K

Let's generate a series of spectra at 4000 K


In [5]:
teff= 4000.
ofes= [-0.3,-0.2,-0.1,0, 0.1, 0.2, 0.3]
synspecs40= []
for ofe in ofes:
    synspecs40.append(simspec(ofe=ofe,teff=teff)[0])


/Library/Python/2.7/site-packages/numpy/core/fromnumeric.py:2507: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`.
  VisibleDeprecationWarning)

/Users/bovy/Repos/apogee/apogee/modelspec/moog.py:179: RuntimeWarning: divide by zero encountered in divide
  out/= cflux

                                                                                

Let's take a look at these spectra


In [7]:
colormap= cm.coolwarm
def color(ofe):
    return colormap((ofe-ofes[0])/(ofes[-1]-ofes[0]))
splot.windows(synspecs40[0],'O1',color=color(ofes[0]))
for ii in range(1,len(synspecs40)): splot.windows(synspecs40[ii],'O1',color=color(ofes[ii]),overplot=True)


Now we can analyze these simulated spectra


In [8]:
results_fparam= []
results_felem= []
for ii in range(len(ofes)):
    tresult= fit(synspecs40[ii])
    results_fparam.append(tresult[0])
    results_felem.append(tresult[1])

In [94]:
results_cannon= []
for ii in range(len(ofes)):
    results_cannon.append(cannon(synspecs40[ii]))

In [90]:
fparam_ref= [4000.,2.5,0.,0.,0.,0.,0.]
results_fparam= numpy.reshape(numpy.array(results_fparam),(len(ofes),7))
newresults_felem= {}
for elem in results_felem[0]: newresults_felem[elem]= numpy.array([results_felem[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam,newresults_felem,fparam_ref,numpy.reshape(numpy.array(results_cannon),(len(ofes),4)),
             var_label=r'$\mathrm{O}$')


Same for a solar-scaled atmosphere (the model Jon ran)


In [10]:
teff= 4000.
ofes= [-0.3,-0.2,-0.1,0, 0.1, 0.2, 0.3]
synspecs40ss= []
for ofe in ofes:
    synspecs40ss.append(simspec(ofe=ofe,teff=teff,solarscaledatm=True)[0])


                                                                                

In [11]:
results_fparam_ss= []
results_felem_ss= []
for ii in range(len(ofes)):
    tresult= fit(synspecs40ss[ii])
    results_fparam_ss.append(tresult[0])
    results_felem_ss.append(tresult[1])

In [92]:
results_cannon_ss= []
for ii in range(len(ofes)):
    results_cannon_ss.append(cannon(synspecs40ss[ii]))

In [93]:
fparam_ref= [4000.,2.5,0.,0.,0.,0.,0.]
results_fparam_ss= numpy.reshape(numpy.array(results_fparam_ss),(len(ofes),7))
newresults_felem_ss= {}
for elem in results_felem_ss[0]: newresults_felem_ss[elem]= numpy.array([results_felem_ss[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam_ss,newresults_felem_ss,fparam_ref,numpy.reshape(numpy.array(results_cannon_ss),(len(ofes),4)),
             var_label=r'$\mathrm{O}$')


$T_\mathrm{eff}$ = 4500 K


In [13]:
teff= 4500.
ofes= [-0.3,-0.2,-0.1,0, 0.1, 0.2, 0.3]
synspecs45= []
for ofe in ofes:
    synspecs45.append(simspec(ofe=ofe,teff=teff)[0])


                                                                                

In [14]:
results_fparam_45= []
results_felem_45= []
for ii in range(len(ofes)):
    tresult= fit(synspecs45[ii])
    results_fparam_45.append(tresult[0])
    results_felem_45.append(tresult[1])

In [97]:
results_cannon_45= []
for ii in range(len(ofes)):
    results_cannon_45.append(cannon(synspecs45[ii]))

In [98]:
fparam_ref= [4500.,2.5,0.,0.,0.,0.,0.]
results_fparam_45= numpy.reshape(numpy.array(results_fparam_45),(len(ofes),7))
newresults_felem_45= {}
for elem in results_felem_45[0]: newresults_felem_45[elem]= numpy.array([results_felem_45[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam_45,newresults_felem_45,fparam_ref,numpy.reshape(numpy.array(results_cannon_45),(len(ofes),4)),
             var_label=r'$\mathrm{O}$')


Block out the [O/H] windows in the fparam fit


In [28]:
results_fparam_45bO= []
results_felem_45bO= []
for ii in range(len(ofes)):
    tresult= fit(synspecs45[ii],blockOwindows=True)
    results_fparam_45bO.append(tresult[0])
    results_felem_45bO.append(tresult[1])

In [95]:
results_cannon_45bO= []
for ii in range(len(ofes)):
    results_cannon_45bO.append(cannon(synspecs45[ii],blockOwindows=True))

In [99]:
fparam_ref= [4500.,2.5,0.,0.,0.,0.,0.]
results_fparam_45bO= numpy.reshape(numpy.array(results_fparam_45bO),(len(ofes),7))
newresults_felem_45bO= {}
for elem in results_felem_45bO[0]: newresults_felem_45bO[elem]= numpy.array([results_felem_45bO[ii][elem] for ii in range(len(ofes))])
plot_results(ofes,results_fparam_45bO,newresults_felem_45bO,fparam_ref,numpy.reshape(numpy.array(results_cannon_45bO),(len(ofes),4)),
             var_label=r'$\mathrm{O}$')


Let's do MCMC on the lowest [O/Fe] synthetic spectrum:


In [26]:
chain= apogee.modelspec.ferre.mcmc(synspecs45[0],0.005*numpy.ones_like(synspecs45[0]),
                                   initsig=[100.,0.25,0.1,0.1,0.1,0.1,0.1],
                                   fparam=numpy.reshape(results_fparam_45[0],(1,7)),nsamples=100)

In [24]:
figsize(8,8)
bovy_plot.scatterplot(chain[:,0],chain[:,-1],'k,')



In [27]:
figsize(8,8)
bovy_plot.scatterplot(chain[:,0],chain[:,-1],'k,')


So, the $T_{\mathrm{eff}} =4600$ K really seems to be the better fit. Let's compare the best-fit spectrum and the input fparam to the data


In [64]:
ispec= apogee.modelspec.ferre.interpolate([4500.,results_fparam_45[0][0]],[2.5,results_fparam_45[0][1]],
                                          [0.,results_fparam_45[0][3]],[-0.2,results_fparam_45[0][6]],
                                          [0.,results_fparam_45[0][5]],[0.,results_fparam_45[0][4]],
                                          lib='msGK',dr='current',apStarWavegrid=True)

In [65]:
print numpy.nanstd(ispec[1]-synspecs45[0])
print numpy.nanstd(ispec[0]-synspecs45[0])


0.00593383625587
0.0131502378388

The first scatter is smaller, so $T_{\mathrm{eff}} = 4600$ K is really the better fit. Below the next two plots I compute the scatter for the well-recovered case of [O/H] = 0, which shows that this level of residual scatter appears typical. A comparison of the spectra and the data follows (blue is the best-fit, green is with the correct $T_{\mathrm{eff}} = 4600$):


In [66]:
splot.waveregions(ispec[1]-synspecs45[0],yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.waveregions(ispec[0]-synspecs45[0],cleanZero=False,overplot=True)



In [67]:
splot.windows(ispec[1]-synspecs45[0],'O1',yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.windows(ispec[0]-synspecs45[0],'O1',cleanZero=False,overplot=True)


For reference, here's what this looks like for a good fit (index 4 is the zero [O/H] offset):


In [68]:
ispecgf= apogee.modelspec.ferre.interpolate([4500.,results_fparam_45[4][0]],[2.5,results_fparam_45[4][1]],
                                          [0.,results_fparam_45[4][3]],[0.,results_fparam_45[4][6]],
                                          [0.,results_fparam_45[4][5]],[0.,results_fparam_45[4][4]],
                                          lib='msGK',dr='current',apStarWavegrid=True)

In [69]:
print numpy.nanstd(ispecgf[1]-synspecs45[4])
print numpy.nanstd(ispecgf[0]-synspecs45[4])


0.00396458072757
0.00603694426985

In [70]:
splot.waveregions(ispecgf[1]-synspecs45[4],yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.waveregions(ispecgf[0]-synspecs45[4],cleanZero=False,overplot=True)



In [71]:
splot.windows(ispecgf[1]-synspecs45[4],'O1',yrange=[-0.05,0.05],cleanZero=False,labelLines=False,ylabel=r'$\Delta f/f_c$')
splot.windows(ispecgf[0]-synspecs45[4],'O1',cleanZero=False,overplot=True)



In [ ]: