%pylab inline
from matplotlib import ticker
import os,sys
sys.path.append('../bok')
from collections import OrderedDict
from astropy.table import Table,join
from astropy.stats import sigma_clip
import bokrmphot
_ = os.environ.setdefault('BOKRMDIR','../bok')

"""
Explanation: Bok Photometric Analysis
This notebook examines the photometric reliability of the Bok data using repeat measurements of SDSS standard stars. For each reference star the mean and rms magnitudes are calculated from all individual measurements. The reference stars are then grouped into bins of magnitude. The distribution of the rms values (the photometric scatter) within the magnitude bins gives an indication of the photometric reliability.
First, some preliminaries:
End of explanation
"""


def load_pstable(psTabFn,photRun=None,season='2014',**kwargs):
    if not os.path.exists(psTabFn):
        psTab = bokrmphot.binned_phot_stats(photRun,season=season,**kwargs)
        psTab.write(psTabFn)
    psTab = Table.read(psTabFn)
    for c in psTab.colnames:
        if c.startswith('sig'):
            psTab[c].format = "{:.4f}"
        elif c=='outlierFrac':
            psTab[c].format = "{:.5f}"
    return psTab

"""
Explanation: Load a photometric statistics table. This table contains the statistics of the photometry within magnitude bins. The "season" keyword restricts the data to a single observing season, so that the statistics can be compared between seasons. The resulting table is saved to a FITS file so that it can be reloaded later.
End of explanation
"""


def plot_compare_scatters(psTabs,labels,colors='rgbcymk',filt='g',ax=None,units='mag',sfx=None):
    if sfx is None:
        sfx = ['','']
    def _append_arr(arr):
        return arr
        # used this for drawstyle=steps-post, but no equiv. for fill_between
        #return np.concatenate([arr,[arr[-1]]])
    if units=='mag':
        mscl = 1.0
    elif units=='mmag':
        mscl = 1e3
    if ax is None:
        figure()
        ax = subplot(111)
    for tab,l,c,_sfx in zip(psTabs,labels,colors,sfx):
        ax.fill_between(tab['mbins'],_append_arr(mscl*tab['sig%s25'%_sfx]),
                        _append_arr(mscl*tab['sig%s75'%_sfx]),
                        edgecolor='none',color=c,alpha=0.5)
        ax.plot(tab['mbins'],_append_arr(mscl*tab['sig%s50'%_sfx]),color=c,lw=1.5,label=l)
    legend(loc='upper left')
    xlabel('%s magnitude' % filt)
    ylabel('per-object std [%s]'%units)
    xlim(17,19)

"""
Explanation: Plot the results from the photometric statistics table. A central line shows the median rms within magnitude bins, and a shaded region shows the inter-quartile range of the rms. I.e., half of all reference stars at a given magnitude had an rms value less than the median line.
End of explanation
"""


psNov2015g = load_pstable('../phot_stats_Nov2015g.fits')
print psNov2015g['mbins','sig50','sig50Clip','outlierFrac']

psJan2017g = load_pstable('../phot_stats_Jan2017g.fits')
print psJan2017g['mbins','sig50','sig50Clip','outlierFrac']

psFeb2017g = load_pstable('../phot_stats_Feb2017g.fits')
print psFeb2017g['mbins','sig50','sig50Clip','outlierFrac']

psJan2017i = load_pstable('../phot_stats_Jan2017i.fits')
print psJan2017i['mbins','sig50','sig50Clip','outlierFrac']

psFeb2017i = load_pstable('../phot_stats_Feb2017i.fits')
print psFeb2017i['mbins','sig50','sig50Clip','outlierFrac']

"""
Explanation: Comparing some early 2014 reductions
The next set of photometric statistics come from some early reductions of the first-year Bok data. Nov2015g was an early run of bokpipe. Note that the first round of difference imaging was based on an IDL  pipeline; the associated catalogs are versioned 20140905 and have an asymptotic scattter of ~2.5% in g band, much worse than the bokpipe runs. Jan2017 and Feb2017 were subsequent runs with incremental improvements; the changes can be recovered from the git repository.
End of explanation
"""


plot_compare_scatters([psNov2015g,psJan2017g],['Nov15g','Jan17g'],'gb')

"""
Explanation: The improvement from the Nov2015g to Jan2017g versions is significant.
End of explanation
"""


plot_compare_scatters([psJan2017g,psFeb2017g],['Jan17g','Feb17g'],'gb')

"""
Explanation: Not much improvement in the Feb2017g version.
End of explanation
"""


plot_compare_scatters([psJan2017i,psFeb2017i],['Jan17i','Feb17i'],'gb')

for psTab,l in zip([psNov2015g,psJan2017g,psFeb2017g],['Nov15g','Jan17g','Feb17g']):
    ii = where(psTab['outlierFrac']>0)[0]
    plot(psTab['mbins'][ii],log10(psTab['outlierFrac'][ii]),label=l)
legend()

for psTab,l in zip([psJan2017i,psFeb2017i],['Jan17i','Feb17i']):
    ii = where(psTab['outlierFrac']>0)[0]
    plot(psTab['mbins'][ii],log10(psTab['outlierFrac'][ii]),label=l)
legend()

"""
Explanation: The scatter increases significantly if fainter stars are included in the zeropoint. Going from mag$<19.5$ to mag$<20.5$ more than doubles the scatter. Do not understand why this is...
End of explanation
"""


procJul2017 = OrderedDict()
for season in ['2014','2015','2016']:
    procJul2017[season] = load_pstable('../archive/summer_2017/phot_stats_%sg_v2.fits'%season,
                                       'cleanstars',season,catdir='../archive/summer_2017/')
    print season
    print procJul2017[season]['mbins','sig50','outlierFrac']

"""
Explanation: 2015/6 data
Jul2017 is the first processing run to include 2015 and 2016 data. The 2014 data was also reprocessed to include all changes to that point. Mainly the improvements were in handling the illumination correction.
End of explanation
"""


plot_compare_scatters(procJul2017.values(),procJul2017.keys(),'gbr')

"""
Explanation: Comparing the results for the three observing seasons. The 2014 data has less scatter, possibly because there are more data and the sky flats are more accurate.
End of explanation
"""


procOct2017 = OrderedDict()
for season in ['2014','2015','2016','2017']:
    procOct2017[season] = load_pstable('../archive/october_2017/phot_stats_%sg_v2.fits'%season,
                                       'cleanstars',season,catdir='../archive/october_2017/')
    print season
    print procOct2017[season]['mbins','sig50','outlierFrac']

"""
Explanation: bokpipe_v0.3.0
The processing version used for the second round of difference imaging is bokpipe_v0.3.0, which was completed in Oct2017. The results from this run are shown below.
End of explanation
"""


figure(figsize=(14,4))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
for pnum,season in enumerate(['2014','2015','2016'],start=1):
    ax = subplot(1,3,pnum)
    plot_compare_scatters([procJul2017[season],procOct2017[season]],
                          [season+'-Jul2017',season+'-Oct2017'],'gb',ax=ax)
    ax.set_ylim(0,0.04)
    ax.grid()

"""
Explanation: This compares the Oct2017 processing to the Jul2017 processing. The main improvement was in cleaning up some gain correction errors, and fixing a long-standing bug where the zeropoints were based on an aperture (15 pixel) that subsequently had an aperture correction applied. A change was made to apply the correction to the zeropoint. Now the 2014 data reach sub-1% repeatability at $g<17.7$!
End of explanation
"""


plot_compare_scatters(procOct2017.values(),procOct2017.keys(),'gbrm')

"""
Explanation: Again comparing the four seasons. The 2017 data mainly have 300s exposure times, whereas previous years mainly had 150s exposures.
End of explanation
"""


procOct2017i = OrderedDict()
for season in ['2014','2015','2016','2017']:
    procOct2017i[season] = load_pstable('../archive/october_2017/phot_stats_%si_v2.fits'%season,
                                       'cleanstars',season,band='i',
                                       catdir='../archive/october_2017/')
    print season
    print procOct2017i[season]['mbins','sig50','outlierFrac']

plot_compare_scatters(procOct2017i.values(),procOct2017.keys(),'gbrm')

"""
Explanation: And now the i-band:
End of explanation
"""


class NewBinnedStats(object):
    def __init__(self,tabf):
        self.tab = Table.read(tabf)
    def get(self,season,filt):
        i = where((self.tab['season']==season)&(self.tab['filter']==filt))[0][0]
        rv = Table()
        rv['mbins'] = array(self.tab.meta['MBINS'].split(',')).astype(float)
        for k in self.tab.colnames:
            if k.startswith('sig'):
                rv[k] = self.tab[k][i]
        return rv

oct17noselfcal = NewBinnedStats('../archive/october_2017/phot_stats_oct17noselfcal.fits')
oct17selfcal = NewBinnedStats('../archive/october_2017/phot_stats_oct17.fits')

figure(figsize=(14,7))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
pnum = 1
for b in 'gi':
    for season in ['2014','2015','2016','2017']:
        ax = subplot(2,4,pnum)
        plot_compare_scatters([oct17noselfcal.get(season,b),oct17selfcal.get(season,b)],
                              [season+'-Oct2017orig',season+'-Oct2017selfcal'],
                              'gb',filt=b,ax=ax,units='mmag')
        if b=='g':
            ax.set_ylim(0,40)
        else:
            ax.set_ylim(0,80)
        ax.grid()
        pnum += 1

"""
Explanation: So the 2016 i-band data has some issues...
self-calibration
As a shortcut to solving the full ubercalibration problem (even with a large number of repeats the Bok data have been found to be somewhat noisy for recovering ubercal parameters), a selfcal() routine was implemented that assumes the extinction parameters from SDSS and uses them to derive per-image zeropoints from the mean magnitude offsets. The above results were all based on zeropoints obtained from external calibration using SDSS; the following results come from the selfcal procedure which provides zeropoints on an internal Bok system.
End of explanation
"""


oct17selfcalphoto = NewBinnedStats('../archive/october_2017/phot_stats_oct17_photo.fits')

figure(figsize=(14,7))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
pnum = 1
for b in 'gi':
    for season in ['2014','2015','2016','2017']:
        ax = subplot(2,4,pnum)
        plot_compare_scatters([oct17selfcal.get(season,b),oct17selfcalphoto.get(season,b)],
                              [season+'-Oct2017selfcal',season+'-Oct2017selfcalphoto'],
                              'br',filt=b,ax=ax,units='mmag')
        if b=='g':
            ax.set_ylim(0,40)
        else:
            ax.set_ylim(0,80)
        ax.grid()
        pnum += 1

"""
Explanation: selfcal has a small effect on g-band photometry (maybe slight improvement in 2015). But huge effect on i-band, and in particular addresses the weird scatter in the 2016 data. Now asymptotically reaching 10 mmag scatter at g,i=17.
Check the improvement by restricting to only photometric frames:
End of explanation
"""


figure(figsize=(14,7))
subplots_adjust(0.05,0.1,0.99,0.99,0.22)
pnum = 1
for b in 'gi':
    for season in ['2014','2015','2016','2017']:
        ax = subplot(2,4,pnum)
        plot_compare_scatters([oct17selfcal.get(season,b),oct17selfcal.get(season,b)],
                              [season+'-Oct2017selfcal-int',season+'-Oct2017selfcal-ext'],
                              ['C0','C1'],filt=b,ax=ax,units='mmag',sfx=['','Ext'])
        if b=='g':
            ax.set_ylim(0,40)
        else:
            ax.set_ylim(0,80)
        ax.grid()
        pnum += 1

"""
Explanation: Finally, compare the internal calib with the external rms (Bok-SDSS):
End of explanation
"""


