<commit_before>#!/usr/bin/env python

import os
import itertools
import numpy as np
from astropy.io import fits
from astropy.table import Table,join,vstack
from astropy.stats import sigma_clip
from astropy.wcs import WCS
import fitsio

import matplotlib.pyplot as plt
from matplotlib import ticker

from bokpipe.bokphot import aper_phot_image

from bass import ampNums,get_amp_index
import ps1cal

test_exptimes = np.array([100.]*6 + [25.,50.,100.,200.,400.])

def srcor(ra1,dec1,ra2,dec2,sep,return_sep=False):
	from astropy.coordinates import SkyCoord,match_coordinates_sky
	from astropy import units as u
	c1 = SkyCoord(ra1,dec1,unit=(u.degree,u.degree))
	c2 = SkyCoord(ra2,dec2,unit=(u.degree,u.degree))
	idx,d2d,d3c = match_coordinates_sky(c1,c2)
	ii = np.where(d2d.arcsec < sep)[0]
	if return_sep:
		return ii,idx[ii],d2d.arcsec[ii]
	else:
		return ii,idx[ii]

def s82_ra_slice(s82,ra_range):
	return np.where((s82['ra']>ra_range[0]) & (s82['ra']<ra_range[1]))[0]

def load_stripe82_truth(ra_range=None):
	s82 = fits.getdata(os.path.join(os.environ['DESITARGTRUTH'],
	                                'stripe82-dr12-stars.fits.gz'))
	if ra_range is not None:
		s82 = s82[s82_ra_slice(s82,ra_range)]
	return s82

def has_coverage(truth,objs):
	return ( (truth['ra']>objs['ALPHA_J2000'].min()) &
	         (truth['ra']<objs['ALPHA_J2000'].max()) &
	         (truth['dec']>objs['DELTA_J2000'].min()) &
	         (truth['dec']<objs['DELTA_J2000'].max()) )

def load_obsdb():
	return Table.read('../basschute/config/nov2015_mod.fits')

def stripe82_catalogs(bokdir,ccdsfile='bass-ccds-idmnov2015.fits'):
	from boketc import k_ext
	alltabs = []
	mrows = []
	fieldNum = 1
	ccdlist = Table.read(ccdsfile)
	# start with the Stripe 82 truth table as a base catalog
	s82truth = load_stripe82_truth()
	# iterate over the two S82 fields observed (Deep2F3 and a random spot)
	for field,ra_range in zip(['s82cal%s_ra334_%s','deep2%s_ra352_%s'],
	                          [(333.5,336.0),(351.0,353.5)]):
		ii = s82_ra_slice(s82truth,ra_range)
		s82 = s82truth[ii]
		# initialize the table with Stripe 82 fields and do some renaming
		# for convenience
		t = Table(s82truth[ii],masked=True)
		for cn in t.colnames:
			cn_new = cn.replace('psfmagerr','err')
			cn_new = cn_new.replace('psfmag','mag')
			t.rename_column(cn,'sdss_'+cn_new)
		t['sdss_idx'] = ii
		# match the Bok catalogs
		for k,dtyp in [('x','f4'),('y','f4'),('ra','f8'),('dec','f8'),
		               ('idx','i4'),('mag','f4'),('err','f4'),
		               ('field','i4'),('ccd','i4'),('amp','i4'),
		               ('kx','f4'),('apcorr','f4')]:
			for b in 'gr':
				t['bok_'+k+'_'+b] = np.zeros((len(ii),11),dtype=dtyp)
		for b in 'gr':
			t['bok_mag_'+b].mask = True
		for n,b in itertools.product(range(11),'gr'):
			sfx = '123456abcde'[n]
			fieldName = field % (b,sfx)
			catf = os.path.join(bokdir,fieldName+'.cat.fits')
			try:
				cat = fits.open(catf)
			except IOError:
				continue
			for ccdNum in range(1,5):
				objs = cat[ccdNum].data
				ii = np.where(has_coverage(s82,objs))[0]
				m1,m2,d = srcor(objs['ALPHA_J2000'],objs['DELTA_J2000'],
				                s82['ra'][ii],s82['dec'][ii],
				                2.0,return_sep=True)
				print ' %20s[%d] matched %4d out of %4d (%.2f)' % \
				       (fieldName,ccdNum,len(m1),len(ii),np.median(d))
				t['bok_idx_'+b][ii[m2],n] = objs['NUMBER'][m1]
				t['bok_field_'+b][ii[m2],n] = fieldNum
				t['bok_ccd_'+b][ii[m2],n] = ccdNum
				ampIndex = get_amp_index(objs['X_IMAGE'][m1],
				                         objs['Y_IMAGE'][m1])
				t['bok_amp_'+b][ii[m2],n] = 4*(ccdNum-1) + ampIndex + 1
				for k1,k2 in [('x','X_IMAGE'),('y','Y_IMAGE'),
				              ('ra','ALPHA_J2000'),('dec','DELTA_J2000')]:
					t['bok_'+k1+'_'+b][ii[m2],n] = objs[k2][m1]
				flux = np.ma.array(objs['FLUX_APER'][m1,1],
				                   mask=((objs['FLUX_APER'][m1,1]<=0) |
				                         (objs['FLAGS'][m1]>0)))
				fluxerr = objs['FLUXERR_APER'][m1,1]
				mag = -2.5*np.ma.log10(flux/test_exptimes[n])
				t['bok_mag_'+b][ii[m2],n] = mag
				t['bok_mag_'+b][ii[m2],n].mask = mag.mask
				t['bok_err_'+b][ii[m2],n] = 1.0856*np.ma.divide(fluxerr,flux)
				# meta-data about the image
				totflux = objs['FLUX_APER'][m1,-1]
				apcorr = objs['FLUX_APER'][m1] / totflux[:,np.newaxis]
				psfcorr = objs['FLUX_PSF'][m1] / totflux
				fwhm = np.median(objs['FWHM_IMAGE'][m1])
				apcorr = sigma_clip(apcorr[:,1]).mean()
				ccdent = np.where((ccdlist['image_filename'] ==
			                                 fieldName+'.fits') &
				                        (ccdlist['ccdnum']==ccdNum))[0][0]
				airmass = ccdlist['airmass'][ccdent]
				mrows.append((fieldName,fieldNum,b,ccdNum,fwhm,apcorr,
				              airmass,ccdlist['avsky'][ccdent],
				              ccdlist['arawgain'][ccdent]))
				# copy corrections into photo tab
				t['bok_kx_'+b][ii[m2],n] = k_ext[b]*(airmass-1)
				t['bok_apcorr_'+b][ii[m2],n] = 2.5*np.log10(apcorr)
			if b=='r':
				fieldNum += 1
		# eliminate entries without Bok coverage in both bands
		ii = np.where(~(t['bok_mag_g'].mask.all(axis=1) &
		                t['bok_mag_r'].mask.all(axis=1)))[0]
		t = t[ii]
		# now bring in PS1
		for k,dtyp in [('ra','f8'),('dec','f8'),('idx','i4')]:
			t['ps1_'+k] = np.zeros(len(t),dtype=dtyp)
		for j,b in enumerate('grizy'):
			t['ps1_mag_'+b] = np.zeros(len(t),dtype='f4')
		ps1objs = ps1cal.get_ps1_stars(t['sdss_ra'],t['sdss_dec'])
		m1,m2 = srcor(ps1objs['RA'],ps1objs['DEC'],
		              t['sdss_ra'],t['sdss_dec'],2.0)
		for k1,k2 in [('ra','RA'),('dec','DEC'),('idx','OBJ_ID')]:
			t['ps1_'+k1][m2] = ps1objs[k2][m1]
		for j,b in enumerate('grizy'):
			t['ps1_mag_'+b][m2] = ps1objs['MEDIAN'][m1,j]
		# 
		alltabs.append(t)
	# combine tables & calculate residual per-amplifier offsets
	tab = vstack(alltabs)
	for b in 'gr':
		mag = tab['bok_mag_'+b] - tab['bok_kx_'+b] + tab['bok_apcorr_'+b]
		mivar = np.ma.power(tab['bok_err_'+b],-2)
		mean_mag = np.ma.average(sigma_clip(mag,axis=1),weights=mivar,axis=1)
		tab['bok_meanmag_'+b] = mean_mag
		tab['bok_ampcorr_'+b] = np.zeros_like(tab['bok_mag_'+b])
		print '\n%s band amp correction (millimag): ' % b
		for amp in range(1,17):
			amp_mags = np.ma.array(mag,mask=(mag.mask|tab['bok_amp_'+b]!=amp))
			amp_dmag = sigma_clip(mean_mag[:,np.newaxis] - amp_mags)
			ii = np.where(tab['bok_amp_'+b] == amp)
			tab['bok_ampcorr_'+b][ii] = amp_dmag.mean()
			tab.meta['AMPCOR%02d'%amp] = amp_dmag.mean()
			print '%4d +/- %2d ' % (1e3*amp_dmag.mean(),1e3*amp_dmag.std()),
			if (amp%4)==0: print
		print
		mag += tab['bok_ampcorr_'+b]
		tab['bok_magcorr_'+b] = mag 
		mean_mag = np.ma.average(sigma_clip(mag,axis=1),weights=mivar,axis=1)
		tab['bok_meanmagcorr_'+b] = mean_mag 
	# something weird here
	for k in ['mag','magcorr','ampcorr']:
		for b in 'gr':
			tab['bok_'+k+'_'+b].fill_value = np.nan
	tab.write('stripe82bok_nov2015stars.fits',overwrite=True)
	#
	metaTab = Table(rows=mrows,
	                names=('field','fieldNum','filter','ccdNum',
	                       'fwhm','apCorr','airmass','skyADU','gain'))
	metaTab.write('stripe82bok_nov2015sum.fits',overwrite=True)

def load_nov2015_data():
	t = Table.read('stripe82bok_nov2015stars.fits')
	# multi-dim Table columns written to FITS lose their masks...
	for k in ['mag','magcorr','meanmag','meanmagcorr']:
		for b in 'gr':
			t['bok_%s_%s'%(k,b)].mask = np.isnan(t['bok_%s_%s'%(k,b)])
	for b in 'grizy':
		t['ps1_mag_'+b].mask = t['ps1_mag_'+b]==0
	return t

def calc_color_terms(doplots=False,savefit=False):<commit_msg>fix keyword<commit_after>#!/usr/bin/env python

import os
import itertools
import numpy as np
from astropy.io import fits
from astropy.table import Table,join,vstack
from astropy.stats import sigma_clip
from astropy.wcs import WCS
import fitsio

import matplotlib.pyplot as plt
from matplotlib import ticker

from bokpipe.bokphot import aper_phot_image

from bass import ampNums,get_amp_index
import ps1cal

test_exptimes = np.array([100.]*6 + [25.,50.,100.,200.,400.])

def srcor(ra1,dec1,ra2,dec2,sep,return_sep=False):
	from astropy.coordinates import SkyCoord,match_coordinates_sky
	from astropy import units as u
	c1 = SkyCoord(ra1,dec1,unit=(u.degree,u.degree))
	c2 = SkyCoord(ra2,dec2,unit=(u.degree,u.degree))
	idx,d2d,d3c = match_coordinates_sky(c1,c2)
	ii = np.where(d2d.arcsec < sep)[0]
	if return_sep:
		return ii,idx[ii],d2d.arcsec[ii]
	else:
		return ii,idx[ii]

def s82_ra_slice(s82,ra_range):
	return np.where((s82['ra']>ra_range[0]) & (s82['ra']<ra_range[1]))[0]

def load_stripe82_truth(ra_range=None):
	s82 = fits.getdata(os.path.join(os.environ['DESITARGTRUTH'],
	                                'stripe82-dr12-stars.fits.gz'))
	if ra_range is not None:
		s82 = s82[s82_ra_slice(s82,ra_range)]
	return s82

def has_coverage(truth,objs):
	return ( (truth['ra']>objs['ALPHA_J2000'].min()) &
	         (truth['ra']<objs['ALPHA_J2000'].max()) &
	         (truth['dec']>objs['DELTA_J2000'].min()) &
	         (truth['dec']<objs['DELTA_J2000'].max()) )

def load_obsdb():
	return Table.read('../basschute/config/nov2015_mod.fits')

def stripe82_catalogs(bokdir,ccdsfile='bass-ccds-idmnov2015.fits'):
	from boketc import k_ext
	alltabs = []
	mrows = []
	fieldNum = 1
	ccdlist = Table.read(ccdsfile)
	# start with the Stripe 82 truth table as a base catalog
	s82truth = load_stripe82_truth()
	# iterate over the two S82 fields observed (Deep2F3 and a random spot)
	for field,ra_range in zip(['s82cal%s_ra334_%s','deep2%s_ra352_%s'],
	                          [(333.5,336.0),(351.0,353.5)]):
		ii = s82_ra_slice(s82truth,ra_range)
		s82 = s82truth[ii]
		# initialize the table with Stripe 82 fields and do some renaming
		# for convenience
		t = Table(s82truth[ii],masked=True)
		for cn in t.colnames:
			cn_new = cn.replace('psfmagerr','err')
			cn_new = cn_new.replace('psfmag','mag')
			t.rename_column(cn,'sdss_'+cn_new)
		t['sdss_idx'] = ii
		# match the Bok catalogs
		for k,dtyp in [('x','f4'),('y','f4'),('ra','f8'),('dec','f8'),
		               ('idx','i4'),('mag','f4'),('err','f4'),
		               ('field','i4'),('ccd','i4'),('amp','i4'),
		               ('kx','f4'),('apcorr','f4')]:
			for b in 'gr':
				t['bok_'+k+'_'+b] = np.zeros((len(ii),11),dtype=dtyp)
		for b in 'gr':
			t['bok_mag_'+b].mask = True
		for n,b in itertools.product(range(11),'gr'):
			sfx = '123456abcde'[n]
			fieldName = field % (b,sfx)
			catf = os.path.join(bokdir,fieldName+'.cat.fits')
			try:
				cat = fits.open(catf)
			except IOError:
				continue
			for ccdNum in range(1,5):
				objs = cat[ccdNum].data
				ii = np.where(has_coverage(s82,objs))[0]
				m1,m2,d = srcor(objs['ALPHA_J2000'],objs['DELTA_J2000'],
				                s82['ra'][ii],s82['dec'][ii],
				                2.0,return_sep=True)
				print ' %20s[%d] matched %4d out of %4d (%.2f)' % \
				       (fieldName,ccdNum,len(m1),len(ii),np.median(d))
				t['bok_idx_'+b][ii[m2],n] = objs['NUMBER'][m1]
				t['bok_field_'+b][ii[m2],n] = fieldNum
				t['bok_ccd_'+b][ii[m2],n] = ccdNum
				ampIndex = get_amp_index(objs['X_IMAGE'][m1],
				                         objs['Y_IMAGE'][m1])
				t['bok_amp_'+b][ii[m2],n] = 4*(ccdNum-1) + ampIndex + 1
				for k1,k2 in [('x','X_IMAGE'),('y','Y_IMAGE'),
				              ('ra','ALPHA_J2000'),('dec','DELTA_J2000')]:
					t['bok_'+k1+'_'+b][ii[m2],n] = objs[k2][m1]
				flux = np.ma.array(objs['FLUX_APER'][m1,1],
				                   mask=((objs['FLUX_APER'][m1,1]<=0) |
				                         (objs['FLAGS'][m1]>0)))
				fluxerr = objs['FLUXERR_APER'][m1,1]
				mag = -2.5*np.ma.log10(flux/test_exptimes[n])
				t['bok_mag_'+b][ii[m2],n] = mag
				t['bok_mag_'+b][ii[m2],n].mask = mag.mask
				t['bok_err_'+b][ii[m2],n] = 1.0856*np.ma.divide(fluxerr,flux)
				# meta-data about the image
				totflux = objs['FLUX_APER'][m1,-1]
				apcorr = objs['FLUX_APER'][m1] / totflux[:,np.newaxis]
				psfcorr = objs['FLUX_PSF'][m1] / totflux
				fwhm = np.median(objs['FWHM_IMAGE'][m1])
				apcorr = sigma_clip(apcorr[:,1]).mean()
				ccdent = np.where((ccdlist['image_filename'] ==
			                                 fieldName+'.fits') &
				                        (ccdlist['ccdnum']==ccdNum))[0][0]
				airmass = ccdlist['airmass'][ccdent]
				mrows.append((fieldName,fieldNum,b,ccdNum,fwhm,apcorr,
				              airmass,ccdlist['avsky'][ccdent],
				              ccdlist['arawgain'][ccdent]))
				# copy corrections into photo tab
				t['bok_kx_'+b][ii[m2],n] = k_ext[b]*(airmass-1)
				t['bok_apcorr_'+b][ii[m2],n] = 2.5*np.log10(apcorr)
			if b=='r':
				fieldNum += 1
		# eliminate entries without Bok coverage in both bands
		ii = np.where(~(t['bok_mag_g'].mask.all(axis=1) &
		                t['bok_mag_r'].mask.all(axis=1)))[0]
		t = t[ii]
		# now bring in PS1
		for k,dtyp in [('ra','f8'),('dec','f8'),('idx','i4')]:
			t['ps1_'+k] = np.zeros(len(t),dtype=dtyp)
		for j,b in enumerate('grizy'):
			t['ps1_mag_'+b] = np.zeros(len(t),dtype='f4')
		ps1objs = ps1cal.get_ps1_stars(t['sdss_ra'],t['sdss_dec'])
		m1,m2 = srcor(ps1objs['RA'],ps1objs['DEC'],
		              t['sdss_ra'],t['sdss_dec'],2.0)
		for k1,k2 in [('ra','RA'),('dec','DEC'),('idx','OBJ_ID')]:
			t['ps1_'+k1][m2] = ps1objs[k2][m1]
		for j,b in enumerate('grizy'):
			t['ps1_mag_'+b][m2] = ps1objs['MEDIAN'][m1,j]
		# 
		alltabs.append(t)
	# combine tables & calculate residual per-amplifier offsets
	tab = vstack(alltabs)
	for b in 'gr':
		mag = tab['bok_mag_'+b] - tab['bok_kx_'+b] + tab['bok_apcorr_'+b]
		mivar = np.ma.power(tab['bok_err_'+b],-2)
		mean_mag = np.ma.average(sigma_clip(mag,axis=1),weights=mivar,axis=1)
		tab['bok_meanmag_'+b] = mean_mag
		tab['bok_ampcorr_'+b] = np.zeros_like(tab['bok_mag_'+b])
		print '\n%s band amp correction (millimag): ' % b
		for amp in range(1,17):
			amp_mags = np.ma.array(mag,mask=(mag.mask|tab['bok_amp_'+b]!=amp))
			amp_dmag = sigma_clip(mean_mag[:,np.newaxis] - amp_mags)
			ii = np.where(tab['bok_amp_'+b] == amp)
			tab['bok_ampcorr_'+b][ii] = amp_dmag.mean()
			tab.meta['AMPC%s%02d'%(b.upper(),amp)] = amp_dmag.mean()
			print '%4d +/- %2d ' % (1e3*amp_dmag.mean(),1e3*amp_dmag.std()),
			if (amp%4)==0: print
		print
		mag += tab['bok_ampcorr_'+b]
		tab['bok_magcorr_'+b] = mag 
		mean_mag = np.ma.average(sigma_clip(mag,axis=1),weights=mivar,axis=1)
		tab['bok_meanmagcorr_'+b] = mean_mag 
	# something weird here
	for k in ['mag','magcorr','ampcorr']:
		for b in 'gr':
			tab['bok_'+k+'_'+b].fill_value = np.nan
	tab.write('stripe82bok_nov2015stars.fits',overwrite=True)
	#
	metaTab = Table(rows=mrows,
	                names=('field','fieldNum','filter','ccdNum',
	                       'fwhm','apCorr','airmass','skyADU','gain'))
	metaTab.write('stripe82bok_nov2015sum.fits',overwrite=True)

def load_nov2015_data():
	t = Table.read('stripe82bok_nov2015stars.fits')
	# multi-dim Table columns written to FITS lose their masks...
	for k in ['mag','magcorr','meanmag','meanmagcorr']:
		for b in 'gr':
			t['bok_%s_%s'%(k,b)].mask = np.isnan(t['bok_%s_%s'%(k,b)])
	for b in 'grizy':
		t['ps1_mag_'+b].mask = t['ps1_mag_'+b]==0
	return t

def calc_color_terms(doplots=False,savefit=False):
{'fraction_non_alphanumeric': 0.2017625231910946, 'fraction_numerical': 0.04267161410018553, 'mean_word_length': 3.8640025781501772, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6305833', 'n_tokens_mistral': 7394, 'n_tokens_neox': 6484, 'n_words': 948}