notebook.community

Edit and run

Notebook to develop S/N QA (mainly figures) [v2]

v2 -- 18.7 and cframe!
v3 -- Refactor SNRFit




In [52]:



    


# imports
from importlib import reload

import numpy as np
from scipy import optimize

from matplotlib import pyplot as plt

import yaml

from desispec.qa import qalib
from desispec.io import read_frame, read_params

Run QA

desi_qa_frame --frame_file=cframe-b1-00000021.fits

Test file





In [2]:



    


tst_file = '/home/xavier/DESI/DESI_SCRATCH/18.7/spectro/redux/mini/QA/exposures/20200316/00000021/qa-b1-00000021.yaml'



    











In [3]:



    


# Read yaml
with open(tst_file, 'r') as infile:
    qa_data = yaml.load(infile)



    


















    






/home/xavier/Projects/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
  This is separate from the ipykernel package so we can avoid doing imports until

Check it





In [4]:



    


qa_data['20200316'][21]['b1']['S2N']['PARAMS']



    


















    
Out[4]:








{'BIN_SZ': 0.1,
 'B_CONT': [[4000.0, 4500.0], [5250.0, 5550.0]],
 'PCHI_RESID': 0.05,
 'PER_RESID': 95.0,
 'R_CONT': [[5950.0, 6200.0], [6990.0, 7230.0]],
 'SKYCONT_ALARM_RANGE': [50.0, 600.0],
 'SKYCONT_WARN_RANGE': [100.0, 400.0],
 'SKYRESID_ALARM_RANGE': [-10.0, 10.0],
 'SKYRESID_WARN_RANGE': [-5.0, 5.0],
 'Z_CONT': [[8120.0, 8270.0], [9110.0, 9280.0]]}

















In [5]:



    


qa_data['20200316'][21]['b1']['S2N']['METRICS'].keys()



    


















    
Out[5]:








dict_keys(['EXPTIME', 'FIDSNR_TGT', 'FILTERS', 'FITCOEFF_TGT', 'FITCOVAR_TGT', 'FIT_FILTER', 'MAGNITUDES', 'MEDIAN_SNR', 'NUM_NEGATIVE_SNR', 'OBJLIST', 'SCIENCE_FIBERID', 'SNR_MAG_TGT', 'SNR_RESID', 'STAR_FIBERID', 'r2'])

















In [6]:



    


qa_data['20200316'][21]['b1']['S2N']['METRICS']['OBJLIST']



    


















    
Out[6]:








['SCIENCE', 'STD']

Short cut





In [7]:



    


s2n_dict = qa_data['20200316'][21]['b1']['S2N']['METRICS']

Re-shape mags





In [8]:



    


all_mags = np.resize(np.array(s2n_dict['MAGNITUDES']), (500,3))
all_mags.shape



    


















    
Out[8]:








(500, 3)

















In [9]:



    


fidx = np.where(np.array(qa_data['20200316'][21]['b1']['S2N']['METRICS']['FILTERS']) == qa_data['20200316'][21]['b1']['S2N']['METRICS']['FIT_FILTER'])[0]
fidx



    


















    
Out[9]:








array([0])

















In [10]:



    


mags = all_mags[:,fidx].flatten()
mags.shape



    


















    
Out[10]:








(500,)

Good sources (some are Infinite!)





In [11]:



    


gd_mag = np.isfinite(mags)
np.sum(~gd_mag)



    


















    
Out[11]:








40

Plots

Scatter me





In [12]:



    


plt.clf()
ax=plt.gca()
#
ax.scatter(mags[gd_mag], np.array(s2n_dict['MEDIAN_SNR'])[gd_mag])
#
plt.show()

Residual





In [13]:



    


exptime = qa_data['20200316'][21]['b1']['S2N']['METRICS']['EXPTIME']
exptime



    


















    
Out[13]:








1068.323891763128

















In [14]:



    


r2= qa_data['20200316'][21]['b1']['S2N']['METRICS']['r2']
funcMap={"linear":lambda x,a,b:a+b*x,
         "poly":lambda x,a,b,c:a+b*x+c*x**2,
         "astro":lambda x,a,b:(exptime*a*x)/np.sqrt(exptime*(a*x+b)+r2)
        }



    











In [15]:



    


fitfunc = funcMap['astro']



    











In [16]:



    


sci_idx = qa_data['20200316'][21]['b1']['S2N']['METRICS']['OBJLIST'].index('SCIENCE')
sci_idx



    


















    
Out[16]:








0

















In [17]:



    


coeff = qa_data['20200316'][21]['b1']['S2N']['METRICS']['FITCOEFF_TGT'][sci_idx]
coeff



    


















    
Out[17]:








[0.02817235631243733, 0.9126296252761991]

















In [18]:



    


x=10**(-0.4*(mags[gd_mag]-22.5))
fit_snr = fitfunc(x,  *coeff)
fit_snr[0:50]



    


















    
Out[18]:








array([ 1.36442528,  0.95790377,  0.49065736,  1.58962013,  1.23940546,
        1.46917792,  0.51455701,  1.36262026,  3.40150129,  1.68522799,
        0.21250464,  1.27527656,  9.59322434,  0.60613509,  0.79580406,
        2.84026592,  0.25478221, 11.27119514,  0.10254784,  1.09194889,
        6.2654671 ,  0.65933352,  0.14025526,  0.37517679,  0.5143374 ,
        6.96764836,  0.78559341,  0.65711137,  0.93635351,  1.63254071,
        0.77759261,  0.25835161,  0.39707204,  0.46172426,  0.81988059,
        3.09893787,  1.16607547,  0.50500013,  0.33444672,  0.35865504,
        0.48262426,  0.88986687,  0.62737275,  2.07023483,  1.08415799,
        0.15792563,  0.46796233,  0.30871621,  0.8111805 ,  0.33784691])

















In [19]:



    


r2



    


















    
Out[19]:








0.0

















In [20]:



    


resid = (np.array(s2n_dict['MEDIAN_SNR'])[gd_mag]-fit_snr)/fit_snr



    











In [21]:



    


plt.clf()
ax=plt.gca()
#
ax.scatter(mags[gd_mag], resid)
#
plt.show()

With fit





In [22]:



    


plt.clf()
ax=plt.gca()
#
ax.scatter(mags[gd_mag], np.array(s2n_dict['MEDIAN_SNR'])[gd_mag])
# Fit
xmag = np.linspace(17., 26., 100)
xval=10**(-0.4*(xmag-22.5))
yval = fitfunc(xval, *coeff)
ax.plot(xmag, yval, '--r')
#
#ax.set_ylim(5., 20.)
#
plt.show()

Note the points at very bright values are likely stars

Refactor

Fitting





In [32]:



    


def fit_astro(flux, A, B):
    return flux*A/np.sqrt(A*flux + B)

Normalize





In [25]:



    


snr_norm = np.array(s2n_dict['MEDIAN_SNR'])[gd_mag] / exptime**(1/2)

Fit





In [29]:



    


flux = 10**(-0.4*(mags[gd_mag]-22.5))



    











In [43]:



    


popt, pcov = optimize.curve_fit(fit_astro, flux, snr_norm, p0=(0.02, 1.))



    


















    






/home/xavier/Projects/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in sqrt
  

















In [44]:



    


popt



    


















    
Out[44]:








array([0.01817394, 0.31367139])

















In [45]:



    


fit_snr = fit_astro(flux, *popt)



    











In [47]:



    


plt.clf()
ax=plt.gca()
#
ax.scatter(mags[gd_mag], np.array(s2n_dict['MEDIAN_SNR'])[gd_mag])
# Fit
xmag = np.linspace(17., 26., 100)
xval=10**(-0.4*(xmag-22.5))
yval = fit_astro(xval, *popt) * exptime**(1/2)
ax.plot(xmag, yval, '--r')
#
ax.set_xlim(18., 25.)
ax.set_ylim(-1., 15.)
#
plt.show()



    


















    
























In [38]:



    


coeff



    


















    
Out[38]:








[0.02817235631243733, 0.9126296252761991]

Test frame

Load





In [91]:



    


cframe_file = '/home/xavier/DESI/DESI_SCRATCH/19.9/spectro/redux/mini/exposures/20200319/00000031/cframe-r6-00000031.fits'



    











In [92]:



    


cframe = read_frame(cframe_file)



    











In [93]:



    


cframe.fibermap[0:5]



    


















    
Out[93]:






Table length=5



TARGETIDDESI_TARGETBGS_TARGETMWS_TARGETSECONDARY_TARGETTARGET_RATARGET_DECTARGET_RA_IVARTARGET_DEC_IVARBRICKIDBRICK_OBJIDMORPHTYPEPRIORITYSUBPRIORITYREF_IDPMRAPMDECREF_EPOCHPMRA_IVARPMDEC_IVARRELEASEFLUX_GFLUX_RFLUX_ZFLUX_W1FLUX_W2FLUX_IVAR_GFLUX_IVAR_RFLUX_IVAR_ZFLUX_IVAR_W1FLUX_IVAR_W2FIBERFLUX_GFIBERFLUX_RFIBERFLUX_ZFIBERFLUX_W1FIBERFLUX_W2FIBERTOTFLUX_GFIBERTOTFLUX_RFIBERTOTFLUX_ZFIBERTOTFLUX_W1FIBERTOTFLUX_W2MW_TRANSMISSION_GMW_TRANSMISSION_RMW_TRANSMISSION_ZEBVPHOTSYSOBSCONDITIONSNUMOBS_INITPRIORITY_INITNUMOBS_MOREHPXPIXELFIBERPETAL_LOCDEVICE_LOCLOCATIONFIBERSTATUSOBJTYPELAMBDA_REFFIBERASSIGN_XFIBERASSIGN_YFA_TARGETFA_TYPENUMTARGETFIBER_RAFIBER_DECFIBER_RA_IVARFIBER_DEC_IVARPLATEMAKER_XPLATEMAKER_YPLATEMAKER_RAPLATEMAKER_DECNUM_ITERSPECTROIDBRICKNAMELAMBDAREFDELTA_XDELTA_Y


int64int64int64int64int64float64float64float64float64int64int64str4int32float64int64float32float32float32float32float32int16float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32str1int32int64int64int32int64int32int32int32int32int32str3float32float32float32int64uint8int16float64float64float32float32float32float32float32float32int32int32str8float64float64float64


2882303983815257732000153.2290649414062532.367282867431640.00.00000.000.00.00.00.00.000.45957150.577162861.15034563.68708523.43857380.00.00.00.00.00.00.00.00.00.00.00.00.00.00.00.9990.9990.9990.001S000003000622362230TGT5400.0-272.08463294.61902000153.2290649414062532.367282867431640.00.00.00.00.00.0261532p3225400.00.00.0


2882303983815280884000153.1161301120841432.291222420156340.00.00000.000.00.00.00.00.001.25295592.45053486.162202415.6647769.6845270.00.00.00.00.00.00.00.00.00.00.00.00.00.00.00.9990.9990.9990.001S000003001649064900TGT5400.0-246.56375273.41003000153.1161301120841432.291222420156340.00.00.00.00.00.0261530p3225400.00.00.0


2882303983815258582000153.225631713867232.325237274169920.00.00000.000.00.00.00.00.001.01556381.01505322.16479474.52152824.06749630.00.00.00.00.00.00.00.00.00.00.00.00.00.00.00.9990.9990.9990.001S000003002627762770TGT5400.0-271.0014283.4717000153.225631713867232.325237274169920.00.00.00.00.00.0261532p3225400.00.00.0


2882303983899161254000152.5028686523437532.29344558715820.00.00000.000.00.00.00.00.005.95674235.523755.933598511.26574816.1264110.00.00.00.00.00.00.00.00.00.00.00.00.00.00.00.9990.9990.9990.001S000003003633763370TGT5400.0-114.49704269.69635000152.5028686523437532.29344558715820.00.00.00.00.00.0261524p3225400.00.00.0


2882303983815280894000153.1521421218731732.302687438180720.00.00000.000.00.00.00.00.000.69616610.939527632.58050946.94464255.07163670.00.00.00.00.00.00.00.00.00.00.00.00.00.00.00.9990.9990.9990.001S000003004627662760TGT5400.0-254.57202276.7705000153.1521421218731732.302687438180720.00.00.00.00.00.0261532p3225400.00.00.0



















In [94]:



    


params = read_params()

Run





In [95]:



    


reload(qalib)
qadict = qalib.SNRFit(cframe, 'r6', params['qa']['s2n'])



    


















    






Starting SNR Fit
Working on ELG
Working on LRG
Working on QSO
Working on MWS
Working on STAR
End SNR Fit

















In [22]:



    


np.array(s2n_dict['MAGNITUDES']).size



    


















    
Out[22]:








1500

















In [21]:



    


np.array(s2n_dict['MEDIAN_SNR']).size



    


















    
Out[21]:








500

Content source: desihub/desispec

Similar notebooks:

notebook.community | gallery | about