A Notebook to Create Plots for Price-Jones & Bovy 2017


In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator,AutoMinorLocator
import spectralspace.sample.access_spectrum as acs
from spectralspace.analysis.empca_residuals import *
from scipy.interpolate import interp1d
import os, glob
from spectralspace.examples.ncells_calculation import *
from spectralspace.examples.comparative_plots import *
from spectralspace.examples.pc_plotter import *
from apogee.tools import wv2pix,pix2wv,toApStarGrid
%pylab inline


This call to matplotlib.use() has no effect because the backend has already
been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.


Populating the interactive namespace from numpy and matplotlib
pylab import has clobbered these variables: ['colors']
`%matplotlib` prevents importing * from pylab and numpy

This call to matplotlib.use() has no effect because the backend has already
been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.


In [2]:
import matplotlib
font = {'family': 'serif',
        'weight': 'normal',
        'size'  :  20
}

matplotlib.rc('font',**font)

default_cmap = 'plasma'
datadir = './data'
datadir = '/geir_data/scr/price-jones/Data/apogee_dim_reduction'
figdir = './fig/'
if not os.path.exists(datadir):
    os.system('mkdir -p {0}'.format(datadir))
    
if not os.path.exists(figdir):
    os.system('mkdir -p {0}'.format(figdir))

Figures

Figure 1 - Hertzsprung-Russell Diagrams

Red clump sample

To generate data to create this plot, run the next box and fill in the prompts with responses given

Which data release? (Enter for 12): <press Enter> Type done at any prompt when finished Data key: TEFF Default is full range. Match or slice? slice Upper limit (Enter for maximum): <press Enter> Lower limit (Enter for minimum): <press Enter> Found good limits And/or? done

Note that this part can take a few minutes


In [69]:
redclump = empca_residuals('apogee','red_clump',maskFilter,ask=True,degree=2,datadir=datadir)


Which data release? (Enter for 12): 
properties  ['DR', 'EMPCA_wrapper', 'R2compare', '__class__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattribute__', '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_basicStructure', '_dataSource', '_getProperties', '_match', '_sampleInfo', '_sampleType', 'applyMask', 'checkArrays', 'continuumNormalize', 'correctUncertainty', 'data', 'directoryClean', 'filterCopy', 'findCorrection', 'findFit', 'findResiduals', 'fitStatistic', 'func_sort', 'getDirectory', 'imshow', 'initArrays', 'logplot', 'makeArrays', 'makeMatrix', 'multiFit', 'numberStars', 'pixelEMPCA', 'plotHistogram', 'plot_example_fit', 'resizePixelEigvec', 'sample_wrapper', 'samplesplit', 'setDeltaR2', 'setR2', 'setR2noise', 'show_sample_coverage', 'uncorrectUncertainty']
Type done at any prompt when finished
Data key: TEFF
Default is full range. Match or slice? slice
Upper limit (Enter for maximum): 
Lower limit (Enter for minimum): 
Found good limits
And/or? done

In [206]:
fig = plt.figure(figsize=(15,10))
# Create red giant subplot
ax = fig.add_subplot(111)
key = 'MEANFIB'
histogram2d(fig,ax,redclump.data['SIGFIB'][redclump.data[key]!=-9999],redclump.data[key][redclump.data[key]!=-9999],norm='log',vmin=1,vmax=5e2)
plt.axvline(15,color='k',lw=2)
plt.xlabel('sigma fiber')
plt.ylabel('mean fiber number')


Out[206]:
<matplotlib.text.Text at 0x27e21310>

In [155]:
redclump.data.dtype


Out[155]:
dtype([('APSTAR_ID', 'S45'), ('ASPCAP_ID', 'S44'), ('APOGEE_ID', 'S18'), ('TELESCOPE', 'S8'), ('LOCATION_ID', '>i2'), ('FIELD', 'S16'), ('J', '>f4'), ('J_ERR', '>f4'), ('H', '>f4'), ('H_ERR', '>f4'), ('K', '>f4'), ('K_ERR', '>f4'), ('RA', '>f8'), ('DEC', '>f8'), ('GLON', '>f8'), ('GLAT', '>f8'), ('APOGEE_TARGET1', '>i4'), ('APOGEE_TARGET2', '>i4'), ('TARGFLAGS', 'S116'), ('NVISITS', '>i4'), ('COMMISS', '>i2'), ('SNR', '>f4'), ('STARFLAG', '>i4'), ('STARFLAGS', 'S129'), ('ANDFLAG', '>i4'), ('ANDFLAGS', 'S59'), ('VHELIO_AVG', '>f4'), ('VSCATTER', '>f4'), ('VERR', '>f4'), ('VERR_MED', '>f4'), ('STABLERV_CHI2', '>f4', (2,)), ('STABLERV_RCHI2', '>f4', (2,)), ('STABLERV_CHI2_PROB', '>f4', (2,)), ('EXTRATARG', '>i2'), ('PARAM', '>f4', (7,)), ('FPARAM', '>f4', (7,)), ('PARAM_COV', '>f4', (7, 7)), ('FPARAM_COV', '>f4', (7, 7)), ('ELEM', '>f4', (15,)), ('FELEM', '>f4', (15,)), ('ELEM_ERR', '>f4', (15,)), ('FELEM_ERR', '>f4', (15,)), ('TEFF', '>f4'), ('LOGG', '>f4'), ('PARAM_M_H', '>f4'), ('PARAM_ALPHA_M', '>f4'), ('TEFF_ERR', '>f4'), ('LOGG_ERR', '>f4'), ('PARAM_M_H_ERR', '>f4'), ('PARAM_ALPHA_M_ERR', '>f4'), ('ASPCAP_CHI2', '>f4'), ('ASPCAP_CLASS', 'S2'), ('ASPCAPFLAG', '>i4'), ('ASPCAPFLAGS', 'S153'), ('PARAMFLAG', '>i4', (7,)), ('AL_H', '>f4'), ('CA_H', '>f4'), ('C_H', '>f4'), ('FE_H', '>f4'), ('K_H', '>f4'), ('MG_H', '>f4'), ('MN_H', '>f4'), ('NA_H', '>f4'), ('NI_H', '>f4'), ('N_H', '>f4'), ('O_H', '>f4'), ('SI_H', '>f4'), ('S_H', '>f4'), ('TI_H', '>f4'), ('V_H', '>f4'), ('AL_H_ERR', '>f4'), ('CA_H_ERR', '>f4'), ('C_H_ERR', '>f4'), ('FE_H_ERR', '>f4'), ('K_H_ERR', '>f4'), ('MG_H_ERR', '>f4'), ('MN_H_ERR', '>f4'), ('NA_H_ERR', '>f4'), ('NI_H_ERR', '>f4'), ('N_H_ERR', '>f4'), ('O_H_ERR', '>f4'), ('SI_H_ERR', '>f4'), ('S_H_ERR', '>f4'), ('TI_H_ERR', '>f4'), ('V_H_ERR', '>f4'), ('AL_H_FLAG', '>i4'), ('CA_H_FLAG', '>i4'), ('C_H_FLAG', '>i4'), ('FE_H_FLAG', '>i4'), ('K_H_FLAG', '>i4'), ('MG_H_FLAG', '>i4'), ('MN_H_FLAG', '>i4'), ('NA_H_FLAG', '>i4'), ('NI_H_FLAG', '>i4'), ('N_H_FLAG', '>i4'), ('O_H_FLAG', '>i4'), ('SI_H_FLAG', '>i4'), ('S_H_FLAG', '>i4'), ('TI_H_FLAG', '>i4'), ('V_H_FLAG', '>i4'), ('ELEM_CHI2', '>f4', (15,)), ('ELEMFLAG', '>i4', (15,)), ('WASH_M', '>f4'), ('WASH_M_ERR', '>f4'), ('WASH_T2', '>f4'), ('WASH_T2_ERR', '>f4'), ('DDO51', '>f4'), ('DDO51_ERR', '>f4'), ('IRAC_3_6', '>f4'), ('IRAC_3_6_ERR', '>f4'), ('IRAC_4_5', '>f4'), ('IRAC_4_5_ERR', '>f4'), ('IRAC_5_8', '>f4'), ('IRAC_5_8_ERR', '>f4'), ('IRAC_8_0', '>f4'), ('IRAC_8_0_ERR', '>f4'), ('WISE_4_5', '>f4'), ('WISE_4_5_ERR', '>f4'), ('TARG_4_5', '>f4'), ('TARG_4_5_ERR', '>f4'), ('AK_TARG', '>f4'), ('AK_TARG_METHOD', 'S17'), ('WASH_DDO51_GIANT_FLAG', '>i2'), ('WASH_DDO51_STAR_FLAG', '>i2'), ('ALL_VISITS', 'S737'), ('VISITS', 'S665'), ('ALL_VISIT_PK', '>i4', (50,)), ('VISIT_PK', '>i4', (50,)), ('J0', '>f8'), ('H0', '>f8'), ('K0', '>f8'), ('METALS', '>f8'), ('ALPHAFE', '>f8'), ('ADDL_LOGG_CUT', '>i4'), ('RC_DIST', '>f8'), ('RC_DM', '>f8'), ('RC_GALR', '>f8'), ('RC_GALPHI', '>f8'), ('RC_GALZ', '>f8'), ('STAT', '>i4'), ('INVSF', '>f8'), ('PMRA', '>f8'), ('PMDEC', '>f8'), ('PMRA_ERR', '>f8'), ('PMDEC_ERR', '>f8'), ('PMMATCH', '>i4'), ('GALVR', '>f8'), ('GALVT', '>f8'), ('GALVZ', '>f8'), ('PMRA_PPMXL', '>f8'), ('PMDEC_PPMXL', '>f8'), ('PMRA_ERR_PPMXL', '>f8'), ('PMDEC_ERR_PPMXL', '>f8'), ('PMMATCH_PPMXL', '>i4'), ('GALVR_PPMXL', '>f8'), ('GALVT_PPMXL', '>f8'), ('GALVZ_PPMXL', '>f8'), ('MEANFIB', '<f4'), ('SIGFIB', '<f4')])

In [186]:
numvisit = np.zeros(len(redclump.data))
for i in range(len(redclump.data)):
    numvisit[i] = len(redclump.data['VISITS'][i].split(','))

In [187]:
fig = plt.figure(figsize=(15,10))
# Create red giant subplot
ax = fig.add_subplot(111)
key = 'TEFF'
histogram2d(fig,ax,redclump.data['SIGFIB'][redclump.data[key]>-9999],redclump.data[key][redclump.data[key]>-9999],norm='log',vmin=1,vmax=5e2)



In [189]:
fig = plt.figure(figsize=(15,10))
# Create red giant subplot
ax = fig.add_subplot(111)
histogram2d(fig,ax,redclump.data['SIGFIB'],numvisit,norm='log',vmin=1,vmax=5e2)



In [53]:
redclump.plotHistogram(redclump.data['SIGFIB'],norm=False,bins=100)
plt.axvline(15)


Out[53]:
<matplotlib.lines.Line2D at 0xa15f950>

In [104]:
fwhminfo = np.load(datadir+'/apogee_dr12_fiberfwhm_atpixel.npy')
fwhms_sample = fwhminfo[(np.round(redclump.matchingData['MEANFIB']-1).astype(int),)]
fwhms = np.median(fwhms_sample,axis=1)

In [133]:
redclump.imshow(fwhminfo)



In [114]:
default_cmap='plasma_r'
fig = plt.figure(figsize=(15,6))
# Create red giant subplot
ax = fig.add_subplot(121)
histogram2d(fig,ax,fwhms_sample.flatten(),np.repeat(redclump.teff,7214),bins=200,norm='log',vmin=1,vmax=1e6)



In [115]:
default_cmap='plasma_r'
fig = plt.figure(figsize=(15,6))
# Create red giant subplot
ax = fig.add_subplot(121)
histogram2d(fig,ax,fwhms_sample.flatten(),np.repeat(redclump.logg,7214),bins=200,norm='log',vmin=1,vmax=1e6)



In [117]:
default_cmap='plasma_r'
fig = plt.figure(figsize=(15,6))
# Create red giant subplot
ax = fig.add_subplot(121)
histogram2d(fig,ax,fwhms_sample.flatten()[np.repeat(redclump.fe_h,7214)!=-9999],np.repeat(redclump.fe_h,7214)[np.repeat(redclump.fe_h,7214)!=-9999],bins=200,norm='log',vmin=1,vmax=1e6)



In [105]:
default_cmap='plasma_r'
fig = plt.figure(figsize=(15,6))
# Create red giant subplot
ax = fig.add_subplot(121)
histogram2d(fig,ax,fwhms,redclump.teff,bins=200,vmax=20)



In [106]:
fig = plt.figure(figsize=(15,6))
# Create red giant subplot
ax = fig.add_subplot(121)
histogram2d(fig,ax,fwhms,redclump.logg,bins=200,vmax=20)



In [107]:
fig = plt.figure(figsize=(15,6))
# Create red giant subplot
ax = fig.add_subplot(121)
histogram2d(fig,ax,fwhms[redclump.fe_h!=-9999],redclump.fe_h[redclump.fe_h!=-9999],bins=200,vmax=15)


Red giant sample

Now run the box below and fill in the prompts with the responses given

Which data release? (Enter for 12): <press Enter> Type done at any prompt when finished Data key: TEFF Default is full range. Match or slice? slice Upper limit (Enter for maximum): <press Enter> Lower limit (Enter for minimum): <press Enter> Found good limits And/or? done

Note that this part can take a few minutes


In [151]:
redgiant = empca_residuals('apogee','red_giant',maskFilter,ask=True,degree=2,datadir=datadir)


Which data release? (Enter for 12): 
properties  ['DR', 'EMPCA_wrapper', 'R2compare', '__class__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattribute__', '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_basicStructure', '_dataSource', '_getProperties', '_match', '_sampleInfo', '_sampleType', 'applyMask', 'checkArrays', 'continuumNormalize', 'correctUncertainty', 'data', 'directoryClean', 'filterCopy', 'findCorrection', 'findFit', 'findResiduals', 'fitStatistic', 'func_sort', 'getDirectory', 'imshow', 'initArrays', 'logplot', 'makeArrays', 'makeMatrix', 'multiFit', 'numberStars', 'pixelEMPCA', 'plotHistogram', 'plot_example_fit', 'resizePixelEigvec', 'sample_wrapper', 'samplesplit', 'setDeltaR2', 'setR2', 'setR2noise', 'show_sample_coverage', 'uncorrectUncertainty']
Type done at any prompt when finished
Data key: TEFF
Default is full range. Match or slice? slice
Upper limit (Enter for maximum): 
Lower limit (Enter for minimum): 
Found good limits
And/or? done

In [101]:
from matplotlib.colors import LogNorm
def histogram2d(fig,ax,x,y,bins=100,clabel=False,vmin=0,vmax=110,norm='lin'):
    """
    Create a 2D histogram of data represented by the two dimensions x and y
    
    fig:      Figure to plot in
    ax:       Subplot object to plot in
    x:        Array of data values in 'x'
    y:        Array of data values in 'y'
    bins:     Number of bins in which to divide each axis
    clabel:   Label for the colourbar - no colourbar is plotted if this is not given
    vmin:     Minimum value of the histogram
    vmax:     Maximum value of the histogram
    
    """
    # Create histogram
    H,xedges,yedges = np.histogram2d(x,y,bins=bins)
    # Reorient appropriately
    H = np.rot90(H)
    H = np.flipud(H)
    # Mask where bins are empty
    Hmasked = np.ma.masked_where(H==0,H)
    # Plot histogram
    if norm == 'lin':
        im = ax.pcolormesh(xedges,yedges,Hmasked,
                           cmap = plt.get_cmap(default_cmap),
                           vmin=vmin,vmax=vmax)
    elif norm == 'log':
            im = ax.pcolormesh(xedges,yedges,Hmasked,
                           cmap = plt.get_cmap(default_cmap),
                           norm=LogNorm(vmin=vmin,vmax=vmax))
    # If colourbar is desired, plot and label it
    if clabel:
        cbar=fig.colorbar(im,pad = 0)
        cbar.set_label(label=clabel,fontsize=20)
        cbar.ax.tick_params(labelsize=20)
    elif not clabel:
        cbar=fig.colorbar(im,pad = 0)
        cbar.ax.tick_params(labelsize=20)

In [153]:
# Choose density colormap
default_cmap = 'viridis'
# Create figure of appropriate size
fig = plt.figure(figsize=(15,6))
# Create red giant subplot
ax = fig.add_subplot(121)
histogram2d(fig,ax,redgiant.teff,redgiant.logg,bins=120,vmax=35)
plt.ylim(4,0)
plt.xlim(6000,3500)
# Don't plot on extremes of axis
plt.xticks(np.arange(4000,6000,500)[::-1],fontsize=20)
# Add minor ticks on both axes
xminorlocator = MultipleLocator(250./2)
ax.xaxis.set_minor_locator(xminorlocator)
yminorlocator = MultipleLocator(0.25)
ax.yaxis.set_minor_locator(yminorlocator)
plt.xticks(np.arange(4000,6000,500)[::-1],fontsize=20)
plt.yticks(fontsize=20)
# Adjust tick thickness and length
plt.tick_params(which='both', width=2)
plt.tick_params(which='major',length=7)
plt.tick_params(which='minor',length=4)
# Add axis labels
plt.ylabel(r'$\log g$',fontsize=26)
plt.xlabel(r'$T_{\mathrm{eff}}\,\,(\mathrm{K})$',fontsize=26)
# Shade sample area
plt.axhline(2,color='k',ls='--',lw=2)
plt.axhline(4,color='k',ls='--',lw=2)
plt.fill_between(np.arange(3500,6100,100),4,2,alpha=0.1,color='k')
# Add sample description
plt.text(5800,0.5,'red giant\nstars',va='top')

# Create red clump subplot
ax = fig.add_subplot(122)
histogram2d(fig,ax,redclump.teff,redclump.logg,clabel='number of stars',bins=80,vmax=125)
plt.ylim(4,0)
plt.xlim(6000,3500)
# Don't plot on extremes of x axis
plt.xticks(np.arange(4000,6000,500)[::-1],fontsize=20)
# Add minor ticks to both axes
xminorlocator = MultipleLocator(250./2)
ax.xaxis.set_minor_locator(xminorlocator)
yminorlocator = MultipleLocator(0.25)
ax.yaxis.set_minor_locator(yminorlocator)
plt.yticks(fontsize=20)
# Adjust tick thickness and length
plt.tick_params(which='both', width=2)
plt.tick_params(which='major',length=7)
plt.tick_params(which='minor',length=4)
# Add axis labels
plt.ylabel(r'$\log g$',fontsize=26)
plt.xlabel(r'$T_{\mathrm{eff}}\,\,(\mathrm{K})$',fontsize=26)
# Shade sample region
plt.axvline(4700,color='k',ls='--',lw=2)
plt.axvline(4900,color='k',ls='--',lw=2)
plt.fill_between(np.arange(4700,5000,100),4,0,alpha=0.1,color='k')
# Label sample
plt.text(5800,0.5,'red clump\nstars',va='top')

# Reduce space between plots and save
plt.subplots_adjust(wspace=0.2)
plt.savefig('{0}/{1}'.format(figdir,'HRdiagram.pdf'))


Figure 2 - Example polynomial fit for NGC 6819 in Teff

Run the box below using the bolded responses to the prompts

Which data release? (Enter for 12): <press Enter> Type done at any prompt when finished Data key: CLUSTER Default is full range. Match or slice? match Match value: N6819 And/or? done


In [154]:
oc = empca_residuals('apogee','clusters',maskFilter,ask=True,degree=2,datadir=datadir)
oc.findResiduals(minStarNum=5,gen=True)


Which data release? (Enter for 12): 
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-154-5c9cb7e9519f> in <module>()
----> 1 oc = empca_residuals('apogee','clusters',maskFilter,ask=True,degree=2,datadir=datadir)
      2 oc.findResiduals(minStarNum=5,gen=True)

/home/price-jones/local/lib/python2.7/site-packages/spectralspace-1.-py2.7.egg/spectralspace/analysis/empca_residuals.pyc in __init__(self, dataSource, sampleType, maskMaker, ask, datadict, datadir, func, badcombpixmask, minSNR, degree, nvecs)
    139         mask.__init__(self,dataSource,sampleType,maskMaker,ask=ask,datadict=datadict,
    140                       minSNR=minSNR,datadir=datadir,func=func,
--> 141                       badcombpixmask=badcombpixmask)
    142         self.degree = degree
    143         self.polynomial = PolynomialFeatures(degree=degree)

/home/price-jones/local/lib/python2.7/site-packages/spectralspace-1.-py2.7.egg/spectralspace/sample/mask_data.pyc in __init__(self, dataSource, sampleType, maskMaker, ask, datadict, datadir, func, badcombpixmask, minSNR)
     93 
     94         """
---> 95         subStarSample.__init__(self,dataSource,sampleType,ask=ask,datadict=datadict,datadir=datadir,func=func)
     96         if isinstance(badcombpixmask,list):
     97             badcombpixmask=np.array(badcombpixmask)

/home/price-jones/local/lib/python2.7/site-packages/spectralspace-1.-py2.7.egg/spectralspace/sample/star_sample.pyc in __init__(self, dataSource, sampleType, ask, datadict, datadir, func)
    543         """
    544         # Create starFilter
--> 545         makeFilter.__init__(self,dataSource,sampleType,ask=ask,datadict=datadict,datadir=datadir,func=func)
    546         import filter_function
    547         reload(filter_function)

/home/price-jones/local/lib/python2.7/site-packages/spectralspace-1.-py2.7.egg/spectralspace/sample/star_sample.pyc in __init__(self, dataSource, sampleType, ask, datadict, datadir, func)
    325 
    326         """
--> 327         starSample.__init__(self,dataSource,sampleType,ask=ask,datadict=datadict)
    328         if ask:
    329             self.done = False

/home/price-jones/local/lib/python2.7/site-packages/spectralspace-1.-py2.7.egg/spectralspace/sample/star_sample.pyc in __init__(self, dataSource, sampleType, ask, datadict)
    138                 os.system('echo RESULTS_VERS $RESULTS_VERS')
    139             self._sampleType = sampleType
--> 140             self._getProperties()
    141 
    142     def _getProperties(self):

/home/price-jones/local/lib/python2.7/site-packages/spectralspace-1.-py2.7.egg/spectralspace/sample/star_sample.pyc in _getProperties(self)
    145         """
    146         if self.DR:
--> 147             self.data = readfn[self._dataSource][self._sampleType]()
    148             if self.DR=='12':
    149                 fib = np.load(os.path.join(os.path.dirname(os.path.realpath(__file__)),'..','data','DR12_supplement','fiberinfo.npy'))

/home/price-jones/local/lib/python2.7/site-packages/spectralspace-1.-py2.7.egg/spectralspace/sample/read_clusterdata.pyc in read_caldata(filename, dr)
     84     data.rename_column('2MASS','ID')
     85     # Now match to allStar to get the location_ids
---> 86     alldata= apread.allStar(raw=True,dr=dr)
     87     locids= numpy.zeros(len(data),dtype='int')-1
     88     hmags= numpy.zeros(len(data),dtype='float')-1

TypeError: allStar() got an unexpected keyword argument 'dr'

In [8]:
oc.plot_example_fit(indep=1,pixel=4313,xlabel='$T_{\mathrm{eff}}$ - med($T_{\mathrm{eff}}$) (K)',figsize=(10,8))
plt.subplots_adjust(left=0.2)
plt.savefig('{0}/example_fit.pdf'.format(figdir))


Figure 3 - Comparing EMPCA model with data

Run the box below using the given responses to the prompts

Which data release? (Enter for 12): <press Enter> Type done at any prompt when finished Data key: TEFF Default is full range. Match or slice? slice Upper limit (Enter for maximum): 4900 Lower limit (Enter for minimum): 4700 Found good limits And/or? and Data key: MEANFIB Default is full range. Match or slice? slice Upper limit (Enter for maximum): 300 Lower limit (Enter for minimum): 100 Found good limits And/or? done

This box will take a few minutes to run


In [6]:
subrc = empca_residuals('apogee','red_clump',maskFilter,ask=True,degree=2,badcombpixmask=7935,datadir = datadir)
subrc.findResiduals(gen=False)


Which data release? (Enter for 12): 12
Type done at any prompt when finished
Data key: TEFF
Default is full range. Match or slice? slice
Upper limit (Enter for maximum): 4900
Lower limit (Enter for minimum): 4700
Found good limits
And/or? and
Data key: MEANFIB
Default is full range. Match or slice? slice
Upper limit (Enter for maximum): 300
Lower limit (Enter for minimum): 100
Found good limits
And/or? done

This box will take a few hours to run


In [ ]:
subrc.pixelEMPCA(nvecs=20,varfunc=meanMed,savename='eig20_minSNR50_corrNone_meanMed.pkl')


       iter        R2             rchi2
EMPCA  1/25       0.01944634      2.22647129
EMPCA  2/25       0.12586238      1.84370218
EMPCA  3/25       0.23497501      1.59595446
EMPCA  4/25       0.25197757      1.56271999
EMPCA  5/25       0.25731060      1.55335551
EMPCA  6/25       0.25893361      1.55134022

In [3]:
subrcmodel = acs.pklread('{0}/red_clump_12_TEFF_up4900.0_lo4700.0_MEANFIB_up300.0_lo100.0/bm7935/eig20_minSNR50_corrNone_meanMed.pkl'.format(datadir))
# Retrieve data
subrcmodel=getarrays(subrcmodel)
# Reconstruct unsaved data from model
subrcmodel=reconstruct_EMPCA_data('{0}/red_clump_12_TEFF_up4900.0_lo4700.0_MEANFIB_up300.0_lo100.0/bm7935'.format(datadir),subrcmodel,minStarNum=5)
# Start with baseline as 
totalapprox = np.tile(np.ma.mean(subrcmodel.residuals,axis=0),(subrcmodel.residuals.shape[0],1))
nvec = 8
for i in range(nvec):
    totalapprox += np.outer(subrcmodel.coeff[:, i], subrcmodel.eigvec[i])

In [4]:
# Choose star
indx=312
# Retrieve plot colours
colors = plt.get_cmap('plasma')(np.linspace(0,0.85,2))
# Set bounds on the apStarGrid
pixup = 8575
pixdown = 0

# Convert spectra to apStarGrid and reapply mask
fitspec = toApStarGrid(subrcmodel.fitspec[indx])
fitspec = np.ma.masked_array(fitspec,mask=fitspec==0)
modelspec = toApStarGrid(totalapprox[indx])
modelspec = np.ma.masked_array(modelspec,mask=modelspec==0)
compspec = toApStarGrid(subrcmodel.spectra[indx])
compspec = np.ma.masked_array(compspec,mask=(modelspec.mask|fitspec.mask))
resspec = toApStarGrid(subrcmodel.residuals[indx])
resspec = np.ma.masked_array(resspec,mask=compspec.mask)
errspec = toApStarGrid(subrcmodel.spectra_errs[indx])
errspec = np.ma.masked_array(errspec,mask=compspec.mask)
# Create array of wavelength values
wvs = pix2wv(np.arange(pixdown,pixup),apStarWavegrid = True)/1e4
# Specify xtick location
stepsize = 0.05
xticks = np.linspace(wvs[0],wvs[-1]+stepsize,stepsize)

# Initialize figure
plt.figure(figsize=(15,5))

# Begin model subplot
fit = plt.subplot2grid((3,1), (0, 0), rowspan=2)
# Plot spectrum and model
fit.plot(wvs,compspec[pixdown:pixup],lw=0.5,color=colors[0],label='original spectrum')
fit.plot(wvs,(modelspec[pixdown:pixup]+fitspec[pixdown:pixup]),lw=0.5,color=colors[-1],label='model spectrum:\n{0} principal components'.format(nvec))
fit.set_xlim(wvs[0],wvs[-1])
# Add minor ticks and labels
xminorlocator = AutoMinorLocator()
fit.xaxis.set_minor_locator(xminorlocator)
yminorlocator = AutoMinorLocator()
fit.yaxis.set_minor_locator(yminorlocator)
fit.set_xticklabels(['']*len(xticks))
fit.set_ylabel('normalized flux')
ymax=1.1
ymin=0.6
yticks = np.linspace(ymin,ymax,6)
fit.set_yticks(yticks)
fit.set_yticklabels(yticks)
# Create legend
legend = fit.legend(loc='best',fontsize=15)
legend.get_frame().set_linewidth(0.0)
for legobj in legend.legendHandles:
    legobj.set_linewidth(4.0)
    
# Begin residual subplot
res = plt.subplot2grid((3,1), (2, 0))
# Plot model residuals
res.plot(wvs,(compspec[pixdown:pixup] - (modelspec[pixdown:pixup]+fitspec[pixdown:pixup])),lw=0.3,color=colors[-1])
# Mark median measurement uncertainty
res.axhline(-np.ma.median(errspec[pixdown:pixup]),color='k')
res.axhline(np.ma.median(errspec[pixdown:pixup]),color='k')
res.set_xlim(wvs[0],wvs[-1])
# Add minor ticks and labels
ymax=0.02
ymin=-0.02
yticks = np.linspace(ymin,ymax,5)
res.set_ylim(ymin,ymax)
res.set_yticks(yticks)
res.set_yticklabels(yticks)
xminorlocator = AutoMinorLocator()
res.xaxis.set_minor_locator(xminorlocator)
yminorlocator = AutoMinorLocator()
res.yaxis.set_minor_locator(yminorlocator)
res.set_ylabel('residuals')
res.set_xlabel('$\lambda\,\,(\mu m)$',fontsize=26)

# Adjust subplot positions and save figure
plt.subplots_adjust(bottom=0.2)
plt.savefig('{0}/model_comp_star{1}.pdf'.format(figdir,indx))


/usr/lib/python2.7/site-packages/ipykernel/__main__.py:24: DeprecationWarning: object of type <type 'float'> cannot be safely interpreted as an integer.