In [1]:
%matplotlib inline
In [2]:
from __future__ import print_function
import pylab as pl
import numpy as np
import rtpipe
import rtlib_cython as rtlib
from astropy import coordinates
from scipy.optimize import curve_fit
import scipy.stats
import emcee, corner
In [3]:
# SNR from realfast using telcal gain solutions
snrdet = {'57623': 35.2, '57633_scan7': 119.3, '57633_scan13': 16.3, '57638': 9.9, '57643': 69.9,
'57645': 11.1, '57646': 16.0, '57648': 17.6, '57649': 25.2}
# flux scale from Bryan
calstring = """2.0520 2.89698 0.00279
2.1800 ******* *******
2.3080 ******* *******
2.4360 3.53585 0.00377
2.5640 3.69554 0.00376
2.6920 3.85507 0.00423
2.8200 4.00438 0.00486
2.9480 4.11069 0.00562
3.0520 4.20375 0.00631
3.1800 4.29385 0.00662
3.3080 4.36557 0.00715
3.4360 4.43684 0.00786
3.5640 4.46937 0.00850
3.6920 4.52488 0.00860
3.8200 4.53571 0.00969
3.9480 4.54625 0.00859"""
# parse flux scale
freq = []
flux = []
eflux = []
for line in calstring.split('\n'):
if '*' not in line:
result = line.split()
freq.append(float(result[0]))
flux.append(float(result[1]))
eflux.append(float(result[2]))
calfreq = np.array(freq)
calflux = np.array(flux)
print(calfreq, calflux)
In [4]:
def getscannum(sdmfile):
sdm = rtpipe.parsesdm.getsdm(sdmfile)
for scan in sdm.scans():
try:
print('Scan {0} binary data file: {1}'.format(scan.idx, scan.bdf.fname))
bdfscan = int(scan.idx)
except IOError:
pass
return bdfscan
def read_cut(sdmfile, scan, segment, dm=558., dt=1, gainfile=None, **kwargs):
if not gainfile:
gainfile = '.'.join(sdmfile.split('.')[:-1] + ['GN'])
st = rtpipe.RT.set_pipeline(sdmfile, scan, dmarr=[dm], dtarr=[dt], flaglist=[('badap', 3., 0.2)],
uvoversample=1.5, gainfile=gainfile, flagantsol=True,
timesub='mean', logfile=False, savecands=False,
savenoise=False, **kwargs)
data = rtpipe.RT.pipeline_reproduce(st, candloc=[segment,0,0,0,0], product='data')
u, v, w = rtpipe.parsesdm.get_uvw_segment(st, segment)
return st, data, u, v, w
def imloc(st, data, u, v, w, integ, verbose=0, weight='uniform'):
if weight == 'uniform':
weightarr = np.ones_like(data)
elif weight == 'uvdist':
uvd = np.sqrt(u**2+v**2)
weightarr = np.ones_like(data)
for bl in range(st['nbl']):
weightarr[:,bl] *= uvd[bl]/uvd.max()
im = rtpipe.RT.sample_image(st, data*weightarr, u, v, w, i=integ)
l,m = rtpipe.RT.calc_lm(st, im)
ra_center, dec_center = st['radec']
ra, dec = rtpipe.reproduce.source_location(ra_center, dec_center, l, m)
ras = ra.split()
decs = dec.split()
loc = coordinates.SkyCoord('{0}h{1}m{2}s'.format(ras[0], ras[1], ras[2]),
'{0}d{1}m{2}s'.format(decs[0], decs[1], decs[2]), frame='icrs')
if verbose:
print('Peak location:')
print(l, m)
print(loc.ra.hms, loc.dec.dms)
return im, loc
def correctdata(st, data, u, v, w, corr='ph,dm', lm = (-3.835e-04,5.406e-04)):
""" lm gives (ra, dec) = (5 31 58.703708986 33 8 52.5067634154)
as quoted in Chatterjee et al (2017)
"""
data2 = data.copy()
if 'ph' in corr:
l1, m1 = lm
rtlib.phaseshift_threaded(data2, st, l1, m1, u, v)
if 'dm' in corr:
rtlib.dedisperse_par(data2, st['freq'], st['inttime'], st['dmarr'][0], [0, st['nbl']])
return data2
# get array2 for bin in array near value
def find_nearest(array, array2, value):
idx = (np.abs(array-value)).argmin()
return array2[idx]
def getscale(st):
# get flux scaling at nearest frequency
scale = []
for i in range(len(st['freq'])):
freq = st['freq'][i]
scale.append(find_nearest(calfreq, calflux, freq))
# print(i, st['freq'][i], scale)
scale = np.array(scale, dtype='complex64')[None,None,:,None]
return scale
def correct_all(st, data, u, v, w):
scale = getscale(st)
dataph = correctdata(st, data*scale, u, v, w, corr='ph')
dataphdm = correctdata(st, dataph, u, v, w, corr='dm')
return dataphdm
def correct_and_plot(st, data, u, v, w, integ):
scale = getscale(st)
dataph = correctdata(st, data*scale, u, v, w, corr='ph')
dataphdm = correctdata(st, dataph, u, v, w, corr='dm')
fig = pl.figure(figsize=(15,8))
ax = fig.add_subplot(131)
sgram = dataphdm[integ-10:integ+10].mean(axis=3).mean(axis=1).real.transpose()
pl.imshow(sgram, interpolation='nearest', origin='lower', cmap='viridis')
pl.xlabel('Integration')
pl.ylabel('Channel')
ax = fig.add_subplot(132)
spectra = dataphdm[integ-10:integ+10].mean(axis=2).mean(axis=1).real
pl.plot(spectra[:,0], 'c.', label='RR')
pl.plot(spectra[:,1], 'm.', label='LL')
pl.legend()
pl.xlabel('Integration')
pl.ylabel('Amplitude')
ax = fig.add_subplot(133)
datadm = correctdata(st, data*scale, u, v, w, corr='dm')
im, loc = imloc(st, datadm, u, v, w, integ, verbose=1)
pl.imshow(im.transpose(), interpolation='nearest', cmap='viridis', aspect='equal')
print('Image SNR {0}'.format(im.max()/im.std()))
pl.xlabel('RA offset (pix)', fontsize=20)
pl.ylabel('Dec offset (pix)', fontsize=20)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=3, color='w')
ax.yaxis.set_tick_params(width=3, color='w')
fig = pl.figure(figsize=(3,6))
ax = fig.add_subplot(111)
sgram = dataph[integ-10:integ+50].mean(axis=3).mean(axis=1).real.transpose()
extent = (0, 60*st['inttime'], st['freq'][0], st['freq'][-1])
pl.imshow(sgram, interpolation='nearest', origin='lower', cmap='viridis',
vmax=sgram.max()*0.8, extent=extent)
pl.xlabel('Time (s, relative)', fontsize=14)
pl.ylabel('Frequency (GHz)', fontsize=14)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=2, color='k')
ax.yaxis.set_tick_params(width=2, color='k')
pl.xticks(np.linspace(0, 60*st['inttime'], 4))
pl.colorbar()
fig = pl.figure(figsize=(6,3))
ax = fig.add_subplot(111)
spectra = dataphdm[integ].mean(axis=0).real
print('(RR-LL)/(RR+LL): {0}'.format(np.mean((spectra[:,0]-spectra[:,1])).real/np.mean((spectra[:,0]+spectra[:,1])).real))
print('(RR+LL): {0}'.format(np.mean((spectra[:,0]+spectra[:,1])).real))
nint, nbl, nch, npol = data.shape
pl.scatter(range(0, nch), spectra[:,0], marker='.', c='c', label='RR', edgecolors='none')
pl.scatter(range(0, nch), spectra[:,1], marker='.', c='m', label='LL', edgecolors='none')
resamp = 16
minch = 0
if nch < 256:
minch = resamp - (256-nch)
print(minch, nch, resamp, spectra[minch:,0].shape)
pl.scatter(range(minch, nch, resamp), spectra[minch:,0].reshape((256-minch)/resamp, resamp).mean(axis=1),
marker='_', c='c', s=500, edgecolors='c')
pl.scatter(range(minch, nch, resamp), spectra[minch:,1].reshape((256-minch)/resamp, resamp).mean(axis=1),
marker='_', c='m', s=500, edgecolors='m')
off = np.ma.masked_equal(dataphdm.mean(axis=1).real[range(0,integ-2)+range(integ+3,len(dataphdm))], 0)
noise = off.std(axis=0)
noiseRR = np.median(noise[:,0])
noiseLL = np.median(noise[:,1])
pl.errorbar(230, spectra.max(), yerr=noiseRR, fmt='c.', ecolor='c')
pl.text(230, spectra.max(), s='RR', verticalalignment='center', horizontalalignment='left')
pl.errorbar(240, spectra.max(), yerr=noiseLL, fmt='m.', ecolor='m')
pl.text(240, spectra.max(), s='LL', verticalalignment='center', horizontalalignment='left')
# pl.errorbar(range(nch), spectra[:,0], yerr=noise[:,0], fmt='c.', ecolor='c')
# pl.errorbar(range(nch), spectra[:,1], yerr=noise[:,1], fmt='m.', ecolor='m')
# pl.legend()
pl.xlabel('Frequency (GHz)', fontsize=14)
pl.ylabel('Flux density (Jy)', fontsize=14)
pl.xlim(-5, nch+5)
pl.xticks(range(0, nch+1, 32), [str(np.round(st['freq'][0]+ch*0.004, 2)) for ch in range(0, nch+1, 32)])
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=2, color='k')
ax.yaxis.set_tick_params(width=2, color='k')
return spectra
def getwidespectrum(st, data, u, v, w, i0, i1):
scale = getscale(st)
dataph = correctdata(st, data*scale, u, v, w, corr='ph')
dataphdm = correctdata(st, dataph, u, v, w, corr='dm')
spectrum = dataphdm[i0:i1].mean(axis=3).mean(axis=1).real
return spectrum
def find_dm(st, data, u, v, w):
scale = getscale(st)
dmorig = st['dmarr']
snrs = []
for dm in range(530, 600, 2):
st['dmarr'] = [dm]
datadm = correctdata(st, data*scale, u, v, w, corr='dm')
datadmph = correctdata(st, datadm, u, v, w, corr='ph')
lc = datadmph.mean(axis=3).mean(axis=2).mean(axis=1).real
integ = np.where(lc.max() == lc)[0][0]
im, loc = imloc(st, datadm, u, v, w, integ)
imc = im.transpose()[-1000:, 1000:]
mad = np.median(np.abs(imc - np.median(imc)))
std = mad/0.67
snr = im.max()/std
snrs.append((dm, integ, snr))
st['dmarr'] = dmorig
snrs = np.array(snrs)
snrmax_ind = np.where(snrs[:,2] == snrs[:,2].max())[0][0]
print('SNR vs DM: {0}'.format(snrs))
print('Max SNR of {0} at DM={1}'.format(snrs[snrmax_ind, 2], snrs[snrmax_ind, 0]))
return snrs[snrmax_ind]
sc = lambda sig: np.sqrt(2*sig**2*np.pi)
lum_s = lambda ld, s: s*1e-23 * (4*np.pi*ld**2) * 5e-3 * 1.024e9 # lum in erg, s in Jy, ld in cm
lum_ext = lambda ld, po: po[0]/sc(po[2])*1e-23 * (4*np.pi*ld**2) * 5e-3 * 2.355*fbin*4*po[2]*1e6 # lum in erg, s is norm fit peak in Jy, ld in cm
In [5]:
read = {}
spectrum = {}
spectrumwide = {}
dmmax = {}
In [6]:
sdmfile = '16A-459_TEST_1hr.57623.72670021991.cut'
key = '57623'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 2, npix_max=4096)
In [7]:
key = '57623'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
In [8]:
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [9]:
sdmfile = '16A-459_TEST_1hr_000.57633.66130137732.scan7.cut'
key = '57633_scan7'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 2, gainfile='16A-459_TEST_1hr_000.57633.66130137732.GN',
npix_max=4096)
In [10]:
key = '57633_scan7'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [11]:
sdmfile = '16A-459_TEST_1hr_000.57633.66130137732.scan13.cut'
key = '57633_scan13'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 2, gainfile='16A-459_TEST_1hr_000.57633.66130137732.GN', npix_max=4096)
In [12]:
key = '57633_scan13'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [13]:
sdmfile = '16A-496_sb32698778_1_02h00m.57638.42695471065.cut'
key = '57638'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 4, npix_max=5500, chans=range(2,256))
In [14]:
key = '57638'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [15]:
sdmfile = '16A-496_sb32698778_1_02h00m_001.57643.38562630787.cut/'
key = '57643'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 4, npix_max=5500)
In [16]:
key = '57643'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [17]:
sdmfile = '16A-496_sb32698778_1_02h00m.57645.38915079861.cut'
key = '57645'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 6, npix_max=6500, chans=range(4,256))
In [18]:
key = '57645'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [19]:
sdmfile = '16A-496_sb32698778_1_02h00m_000.57646.38643644676.cut'
key = '57646'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 6, npix_max=7000, chans=range(3,256))
In [20]:
key = '57646'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [21]:
sdmfile = '16A-496_sb32698778_1_02h00m_000.57648.37452900463.cut/'
key = '57648'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 7, npix_max=7400, chans=range(2,256))
In [22]:
key = '57648'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [23]:
sdmfile = '16A-496_sb32698778_1_02h00m_001.57649.37461215278.cut/'
key = '57649'
scannum = getscannum(sdmfile)
read[key] = read_cut(sdmfile, scannum, 7, npix_max=7400)
In [24]:
key = '57649'
st, data, u, v, w = read[key]
snrmax = find_dm(st, data, u, v, w)
dmmax[key] = snrmax
dm = snrmax[0]
integ = int(snrmax[1])
st['dmarr'] = [dm]
spectrum[key] = correct_and_plot(st, data, u, v, w, integ)
st['dmarr'] = [560.5]
_ = correct_and_plot(st, data, u, v, w, integ)
#i0 = integ - 1
#i1 = integ + 2
#spectrumwide[key] = getwidespectrum(st, data, u, v, w, i0, i1)
print('Image SNR (orig): {0}'.format(snrdet[key]))
In [25]:
# full-res spectral fits to norm function
mul = 1
poptg = {'57623': [ mul*33.18660225 , 68.6786758 , 31.97989185],
'57633_scan7': [ mul*261.01468348, 163.37613653 , 54.94315806],
'57633_scan13': [ mul*25.91234327 , 0. , 36.83456147],
'57643': [ mul*54.70339027 , 68.70560927, 55.53711055],
'57649': [ mul*17.31563655, 104.88902813 , 93.62882441],
'57638': [ mul*8.03104856, 159.12796452, 43.9087101 ],
'57646': [ mul*1.67436907e+01, 0., 4.26219468e+01],
'57645': [ mul*3.63342113, 85.9990355, 22.58927418],
'57648': [ mul*12.25215549, 87.65698285, 44.20621731]}
# same, but bin spectra by 4
poptg4 = {'57633_scan7': [ 65.34364137, 40.53096697, 13.78088794],
'57649': [ 4.3092322, 25.87960046, 23.20757215],
'57643': [ 13.77221702, 16.65750292, 13.9484172 ],
'57638': [ 2.02851581, 39.29549388, 11.28088088],
'57646': [ 3.98608102e+00, 0, 1.00213376e+01],
'57633_scan13': [ 6.40875573e+00, 0, 9.26735728e+00],
'57645': [ 1.17831665, 17.02810437, 9.24786317],
'57648': [ 3.05757414, 21.54651779, 11.21822984],
'57623': [ 8.30028279, 16.80016577 , 8.00799541]}
In [26]:
spec_ci = {}
thetas = {}
post = {}
In [27]:
samples = {}
In [28]:
#######
# 2d modeling with emcee
######
def normdm(freq, dx, dy, a, loc, scale, i0, dm, di, index=2, inttime=0.005):
model = np.zeros((dy, dx), dtype='float')
chans = np.arange(dy)
normarr = a*scipy.stats.norm(loc, scale).pdf(chans)
# old way
# delay = (4.1488e-3 * dm * (1./freq**index - 1./freq[-1]**index))/inttime
# i = np.round(i0 + delay, 0).astype(int) # or np.floor?
# new way
i = i0 + (4.1488e-3 * dm * (1./freq**index - 1./freq[-1]**index))/inttime
i_f = np.floor(i).astype(int)
imax = np.ceil(i + di).astype(int)
imin = i_f
i_r = imax - imin
# print(i_r)
if np.any(i_r == 1):
ir1 = np.where(i_r == 1)
# print(ir1)
model[chans[ir1], i_f[ir1]] += normarr[ir1]
if np.any(i_r == 2):
ir2 = np.where(i_r == 2)
i_c = np.ceil(i).astype(int)
f1 = (di - (i_c - i))/di
f0 = 1 - f1
# print(np.vstack((ir2, f0[ir2], f1[ir2])).transpose())
model[chans[ir2], i_f[ir2]] += f0[ir2]*normarr[ir2]
model[chans[ir2], i_f[ir2]+1] += f1[ir2]*normarr[ir2]
if np.any(i_r == 3):
ir3 = np.where(i_r == 3)
f2 = (i + di - (imax - 1))/di
f0 = ((i_f + 1) - i)/di
f1 = 1 - f2 - f0
# print(np.vstack((ir3, f0[ir3], f1[ir3], f2[ir3])).transpose())
model[chans[ir3], i_f[ir3]] += f0[ir3]*normarr[ir3]
model[chans[ir3], i_f[ir3]+1] += f1[ir3]*normarr[ir3]
model[chans[ir3], i_f[ir3]+2] += f2[ir3]*normarr[ir3]
return model
def lnlike(theta, freq, z, zerr):
a, loc, scale, i0, dm, di = theta
dy, dx = z.shape
model = normdm(freq, dx, dy, a, loc, scale, i0, dm, di)
# return -0.5*(np.sum(((z-model)/np.where(model > zerr, model, zerr))**2))
return -0.5*(np.sum(((z-model)/zerr)**2))
def lnprior(theta, imin=0, imax=150, sigma_im=8):
index = 2
a, loc, scale, i0, dm, di = theta
if ((np.sqrt(di/1)*a/sc(scale)/0.005 > sigma_im) and (1 < a < 1000) and (0 < loc < 256) and (1 < scale < 256) and
(imin < i0 < imax) and (530 < dm < 600) and (0.01 < di < 2) and (1.5 < index < 2.5)):
return 0.
return -np.inf
def lnprob(theta, freq, z, zerr, sigma_im):
imin = 0
imax = z.shape[1]-40 # sweep for dm~560
lp = lnprior(theta, imin=imin, imax=imax, sigma_im=sigma_im)
if not np.isfinite(lp):
return -np.inf
return lp + lnlike(theta, freq, z, zerr)
def madstd(sgram, imax):
return np.ma.median(np.ma.abs(sgram[:,:imax] - np.ma.median(sgram[:,:imax])))/0.67
delay = lambda dm, index, freq: (4.1488e-3 * dm * (1./freq**index - 1./freq[-1]**index))
ndim, nwalkers, nsteps = 6, 100, 700
sigma_im = 8
pct = [16, 50, 84]
#pct = [2.5, 50, 97.5]
di = 0.5
In [29]:
key = '57623'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(15,8))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
pl.tight_layout()
In [30]:
key = '57623'
pos = [np.array(thetas[key]) + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(12,7))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [31]:
samples[key] = sampler.chain[:, 250:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["a", "loc", "scale", "i0", "dm", "di"], show_titles=True)
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [32]:
#key = '57623'
#a, loc, scale, i0, dm, di = thetas[key]
a, loc, scale, i0, dm, di = samples[key][100]
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, a, loc, scale, i0, dm, di)
zerr = 0.066
err = np.where(model <= zerr, zerr, model)
to = 1
weight = ((sgram-model)/np.where(model > to*zerr, model, zerr))**2
fig = pl.figure(figsize=(12,12))
ax = fig.add_subplot(141)
pl.imshow(model[:,95:135], origin='bottom', aspect='auto')
ax = fig.add_subplot(142)
pl.imshow(sgram[:,95:135], origin='bottom', aspect='auto')
ax = fig.add_subplot(143)
pl.imshow(err[:,95:135], origin='bottom', aspect='auto')
ax = fig.add_subplot(144)
pl.imshow(weight[:,95:135], origin='bottom', aspect='auto')
pl.tight_layout()
print(-0.5*np.sum(weight))
In [33]:
key = '57633_scan7'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[33]:
In [34]:
pos = [np.array(thetas[key]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [35]:
samples[key] = sampler.chain[:, 380:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["a", "loc", "scale", "i0", "dm", "di"], show_titles=True)
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [36]:
key = '57633_scan13'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[36]:
In [37]:
pos = [np.array(thetas[key]) + 1e-3*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [38]:
samples[key] = sampler.chain[:, 300:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["a", "loc", "scale", "i0", "dm", "di"], show_titles=True)
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [39]:
key = '57638'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[39]:
In [40]:
pos = [np.array(thetas[key]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [41]:
samples[key] = sampler.chain[:, 150:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["a", "loc", "scale", "i0", "dm", "di"], show_titles=True)
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [42]:
key = '57643'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[42]:
In [43]:
pos = [np.array(thetas[key]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [44]:
samples[key] = sampler.chain[:, 300:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["a", "loc", "scale", "i0", "dm", "di"], show_titles=True)
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [45]:
key = '57645'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = (10, 90, 40) # poptg[key] # force to wider value than nominal Gaussian fit
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[45]:
In [46]:
pos = [np.array(thetas[key]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [47]:
samples[key] = sampler.chain[:, 250:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["a", "loc", "scale", "i0", "dm", "di"], show_titles=True)
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [48]:
key = '57646'
di = 0.5
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[48]:
In [49]:
key = '57646'
pos = [np.array(thetas[key]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [50]:
key = '57646'
samples[key] = sampler.chain[:, 250:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["amp", "center", "$\sigma_{ch}$", "i0", "DM", "W$_{int}$"],
label_kwargs={"fontsize": 16})
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
fig.savefig('corner57646.pdf', format='pdf')
In [51]:
a, loc, scale, i0, dm, di = samples[key][100]
print(a, loc, scale, i0, dm, di)
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, a, loc, scale, i0, dm, di)
zerr = std
err = np.where(model <= zerr, zerr, model)
weight = ((sgram-model)/zerr)**2
fig = pl.figure(figsize=(12,12))
ax = fig.add_subplot(141)
pl.imshow(model, origin='bottom', aspect='auto')
ax = fig.add_subplot(142)
pl.imshow(sgram, origin='bottom', aspect='auto')
ax = fig.add_subplot(143)
pl.imshow(err, origin='bottom', aspect='auto')
ax = fig.add_subplot(144)
pl.imshow(weight, origin='bottom', aspect='auto')
pl.tight_layout()
print(-0.5*np.sum(weight))
In [52]:
key = '57648'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[52]:
In [53]:
pos = [np.array(thetas[key]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [54]:
samples[key] = sampler.chain[:, 200:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["a", "loc", "scale", "i0", "dm", "di"], show_titles=True)
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [55]:
a, loc, scale, i0, dm, di = samples[key][100]
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, a, loc, scale, i0, dm, di)
zerr = std
err = np.where(model <= zerr, zerr, model)
to = 1
weight = ((sgram-model)/np.where(model > to*zerr, model, zerr))**2
fig = pl.figure(figsize=(12,12))
ax = fig.add_subplot(141)
pl.imshow(model, origin='bottom', aspect='auto')
ax = fig.add_subplot(142)
pl.imshow(sgram, origin='bottom', aspect='auto')
ax = fig.add_subplot(143)
pl.imshow(err, origin='bottom', aspect='auto')
ax = fig.add_subplot(144)
pl.imshow(weight, origin='bottom', aspect='auto')
pl.tight_layout()
print(-0.5*np.sum(weight))
In [56]:
key = '57649'
st, data, u, v, w = read[key]
scale = getscale(st)
sgram = correctdata(st, data*scale, u, v, w, corr='ph').mean(axis=3).mean(axis=1).real.transpose()
dm, i0, snr = dmmax[key]
a, loc, scale = poptg[key]
thetas[key] = a, loc, scale, i0, dm, di
dy, dx = sgram.shape
model = normdm(st['freq'], dx, dy, *thetas[key])
fig = pl.figure(figsize=(13,5))
ax = fig.add_subplot(121)
pl.imshow(model, interpolation='nearest', origin='bottom')
ax = fig.add_subplot(122)
pl.imshow(sgram, interpolation='nearest', origin='bottom')
Out[56]:
In [57]:
pos = [np.array(thetas[key]) + 1e-2*np.random.randn(ndim) for i in range(nwalkers)]
std = madstd(sgram, 10)
print('Noise from MAD: {0}'.format(std))
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(st['freq'], sgram, std, sigma_im))
end = sampler.run_mcmc(pos, nsteps)
fig = pl.figure(figsize=(10,5))
ax = fig.add_subplot(611)
plot = ax.plot(sampler.chain[::10,::10,0].transpose(), ',k')
ax = fig.add_subplot(612)
plot = ax.plot(sampler.chain[::10,::10,1].transpose(), ',k')
ax = fig.add_subplot(613)
plot = ax.plot(sampler.chain[::10,::10,2].transpose(), ',k')
ax = fig.add_subplot(614)
plot = ax.plot(sampler.chain[::10,::10,3].transpose(), ',k')
ax = fig.add_subplot(615)
plot = ax.plot(sampler.chain[::10,::10,4].transpose(), ',k')
ax = fig.add_subplot(616)
plot = ax.plot(sampler.chain[::10,::10,5].transpose(), ',k')
In [58]:
samples[key] = sampler.chain[:, 300:, :].reshape((-1, ndim))
fig = corner.corner(samples[key], labels=["Amp", "Center", "$\sigma_ch$", "int", "DM", "dint"],
show_titles=True, label_kwargs={"fontsize": 16})
print(sampler.acceptance_fraction.mean())
post[key] = zip(*np.percentile(samples[key], pct, axis=0))
In [59]:
dts = np.array([(5*post[key][5][1], 5*(post[key][5][1]-post[key][5][0]),
5*(post[key][5][2]-post[key][5][1])) for key in sorted(post)])
dms = np.array([(post[key][4][1], post[key][4][1]-post[key][4][0],
post[key][4][2]-post[key][4][1]) for key in sorted(post)])
print(dms, sorted(post.keys()))
print(dts, sorted(post.keys()))
[(key, post[key][5]) for key in post]
Out[59]:
In [60]:
fig= pl.figure(figsize=(10,10))
ax = fig.add_subplot(211)
i=0
keys = ['57623', '57633_scan7', '57633_scan13', '57638', '57643', '57645', '57646', '57648', '57649']
times = [57623, 57633.6, 57633.8, 57638, 57643, 57645, 57646, 57648, 57649]
labels = dict(zip(keys, times))
#dms = []
#for key in keys:
# a, loc, scale, i0, dm, di = post[key]
# dms.append(dm)
#dms = np.array(dms)
#dmerr = np.array([(dm[1]-dm[0], dm[2]-dm[1]) for dm in dms]).transpose()
pl.errorbar(dms[:,0], times, xerr=dms[:,1:].transpose(), fmt='kx')
pl.plot([560.5, 560.5], [times[0], times[-1]], 'k--')
#ax.axhspan(556, 562, alpha=0.3, color='k')
#pl.xlim(-0.3,8.8)
ax.set_xlim(555, 575)
#pl.xlabel('DM (pc cm$^{-3}$)', fontsize=14)
pl.ylabel('Burst Epoch (MJD)', fontsize=16)
xt = pl.setp(ax.get_xticklabels(), fontsize=0)
yt = pl.setp(ax.get_yticklabels(), fontsize=16)
ax.xaxis.set_tick_params(width=3, color='k')
ax.yaxis.set_tick_params(width=3, color='k')
ax2 = fig.add_subplot(212, sharex=ax)
ax2.errorbar(dms[:,0], dts[:,0], yerr=dts[:,1:].transpose(), xerr=dms[:,1:].transpose(), fmt='k.')
ymin, ymax = ax2.get_ylim()
pl.plot([560.5, 560.5], [ymin, ymax], 'k--')
pl.ylabel('Temporal width (ms)', fontsize=16)
pl.xlabel('DM (pc/cm$^3$)', fontsize=16)
xt = pl.setp(ax2.get_xticklabels(), fontsize=16)
yt = pl.setp(ax2.get_yticklabels(), fontsize=16)
ax2.xaxis.set_tick_params(width=3, color='k')
ax2.yaxis.set_tick_params(width=3, color='k')
pl.tight_layout()
fig.savefig('burst_dmdt.pdf', format='pdf')
In [61]:
dme = dms[:,1:].mean(axis=1)
print( (dms[:,0]/dme**2).sum()/(1./dme**2).sum(), np.sqrt(1./(1./dme**2).sum()))
In [62]:
((st['freq'][-1]-st['freq'][0])/delay(5, 2, st['freq'])[0])**(-1)
Out[62]:
In [63]:
def norm(x, a, loc, scale):
return a*scipy.stats.norm(loc, scale).pdf(x)
def fitandplot(spectrum, f0, popt_fixed=None, axis=(), label=None, noise=None):
# create binned spectrum
nch0 = len(spectrum)
spec = []
for ch in range(0, nch0, fbin):
spec.append(spectrum[ch:ch+fbin].mean()) # last bin may have fewer than 4 channels. go with it.
spec = np.array(spec)
nch = len(spec)
print('Binned from {0} to {1} channels'.format(nch0, nch))
if not popt_fixed:
bounds = ([1, 0, 1], [2000, 256/fbin, 50])
popt, pcov = curve_fit(norm, np.arange(nch), spec, bounds=bounds)
else:
popt = popt_fixed
normsol = norm(np.linspace(0, nch, nch0), *popt)
fig = pl.figure(figsize=(15,7))
ax = fig.add_subplot(111)
pl.plot(range(nch0), spectrum, 'k.')
pl.plot(np.arange(nch0), normsol, 'k')
if label:
pl.text(0.85, 0.89, label, horizontalalignment='left', fontsize=24,
verticalalignment='center', transform=ax.transAxes)
if noise:
pl.errorbar(0.85*nch0, 0.93*spectrum.max(), yerr=noise, fmt='k.', ecolor='k')
if len(axis):
pl.axis(axis)
pl.xlabel('Frequency (GHz)', fontsize=18)
pl.ylabel('Flux density (Jy)', fontsize=18)
pl.xlim(-4, nch0+4)
pl.xticks(range(0, nch0+1, 32), [str(np.round(f0+ch*0.004, 2)) for ch in range(0, nch0+1, 32)])
xt = pl.setp(ax.get_xticklabels(), fontsize=18)
yt = pl.setp(ax.get_yticklabels(), fontsize=18)
ax.xaxis.set_tick_params(width=4, color='k')
ax.yaxis.set_tick_params(width=4, color='k')
fig.savefig('spec_{0}.pdf'.format(label), format='pdf')
return popt
In [64]:
fbin = 1
noise = 0.07 # typical 1 sigma noise per channel per integration in Jy
for key in sorted(snrdet.keys()):
print(key)
st, data, u, v, w = read[key]
spec = spectrum[key].mean(axis=1)
a, loc, scale, i0, dm, di = post[key]
fitandplot(spec, st['freq'][0], popt_fixed=(a[1], loc[1], scale[1]), noise=noise, label=labels[key])
print('Already fit {0} (E, DM, peak, center, fwhm): {1:2.1f} & {2:.1f} & {3:3.0f} & {4:1.1f} & {5:.0f}'.format(
key,
np.round(lum_ext(972e6*3.1e18, (a[1], loc[1], scale[1]))/1e38, 1),
dm[1],
np.round(1e3*a[1]/sc(scale[1]), -1),
2.5+fbin*0.004*loc[1],
np.round(2.355*fbin*4*scale[1], -1)))
In [65]:
def fitanddiff(spectrum, popt, noise=0.07):
# create binned spectrum
nch = len(spectrum)
spec = spectrum.mean(axis=1)
normsol = norm(np.linspace(0, nch, nch), *popt)
# diff = (spec-normsol)/(normsol+noise)
# diff = (spec-normsol)/normsol
diff = spec/normsol
return diff[np.where(normsol > 5*noise)]
In [66]:
fluxes = []
diffs = []
for key in snrdet.keys():
a, loc, scale, i0, dm, di = post[key]
popt = (a[1], loc[1], scale[1])
diff = fitanddiff(spectrum[key], popt)
diffs.append(diff)
# fluxes.append(spectrum[key].mean(axis=1))
diffs = np.concatenate(diffs)
#fluxes = np.concatenate(fluxes)
In [67]:
counts, bins = pl.histogram(diffs, bins=50)
binc = [(bins[i+1]+bins[i])/2 for i in range(len(bins)-1)]
pl.plot(binc, counts.astype(float)/counts.sum(), 'k.')
Out[67]:
In [68]:
# https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution
from scipy.special import erfc
emnd = lambda x, mu, sig, ll: (ll/2.)*np.exp((ll/2.)*(2*mu+ll*sig**2 - 2.*x))*erfc((mu+ll*sig**2-x)/(np.sqrt(2)*sig))
In [69]:
counts, bins = pl.histogram(diffs, bins=50)
binc = [(bins[i+1]+bins[i])/2 for i in range(len(bins)-1)]
pl.plot(binc, counts.astype(float)/counts.sum(), 'k.')
model = emnd(np.array(binc), 0.4, 0.3, 2)
pl.plot(binc, model/model.sum(), 'k')
model2 = emnd(np.array(binc), 1, 1, 20) # sig=1?
pl.plot(binc, model2/model2.sum(), 'r')
Out[69]:
In [70]:
# how much of Gaussian model is actually measured within S band?
for key in sorted(snrdet.keys()):
a, loc, scale, i0, dm, di = post[key]
nch = len(spectrum[key])
popt = (a[1], loc[1], scale[1])
inband = norm(np.linspace(0, nch, nch), *popt)
total = norm(np.linspace(-100, nch+100, nch+200), *popt)
print(labels[key], inband.sum()/total.sum())
In [71]:
def plot_acf(st, data, u, v, w, integ, popt=None):
fig = pl.figure(figsize=(8,5))
ax = fig.add_subplot(111)
dataphdm = correct_all(st, data, u, v, w)
spec = dataphdm[integ].mean(axis=2).mean(axis=0).real
spec_off = dataphdm[integ+3].mean(axis=2).mean(axis=0).real
nch = len(spec)
if popt:
spec_mod = norm(np.arange(nch), *popt)
ac = np.correlate(spec, spec, mode='same')
if popt:
ac_mod = np.correlate(spec_mod, spec_mod, mode='same')
ac_off = np.correlate(spec_off, spec_off, mode='same')
pl.plot(ac, 'k.', label="Burst {0}".format(labels[key]))
pl.plot(ac_off, 'k--', label="Off burst")
if popt:
pl.plot(ac_mod, 'k.', label="Model burst", ms=2)
pl.xlabel('Delay (MHz)', fontsize=14)
pl.ylabel('Autocorrelation power (Jy$^2$)', fontsize=14)
xt = pl.setp(ax.get_xticklabels(), fontsize=14)
yt = pl.setp(ax.get_yticklabels(), fontsize=14)
ax.xaxis.set_tick_params(width=2, color='k')
ax.yaxis.set_tick_params(width=2, color='k')
pl.xticks(np.linspace(0, 256, 8), [str(int(4.0*(ch-128))) for ch in np.linspace(0, 256, 8)])
pl.xlim(127.,256)
pl.legend()
In [72]:
for key in snrdet.keys():
print(key)
integ = int(dmmax[key][1])
st, data, u, v, w = read[key]
a, loc, scale, i0, dm, di = post[key]
popt = (a[1], loc[1], scale[1])
plot_acf(st, data, u, v, w, integ, popt=popt)
break
In [73]:
win = 16
specmod = lambda da: np.sqrt((np.mean(da**2) - np.mean(da)**2)/(np.mean(da)**2))
modi = lambda da: da.mean()/da.std()
In [74]:
sms = {}
mis = {}
mns = {}
keys = sorted(spectrum.keys())
for key in keys:
spec = spectrum[key].mean(axis=1)
sm = np.array([specmod(spec[i*win:(i+1)*win]) for i in range(256/win)])
mi = np.array([modi(spec[i*win:(i+1)*win]) for i in range(256/win)])
mn = np.array([spec[i*win:(i+1)*win].mean() for i in range(256/win)])
sms[key] = sm
mis[key] = mi
mns[key] = mn
In [75]:
fig = pl.figure(figsize=(12,12))
for i in range(len(keys)):
ax = fig.add_subplot(3,3,i+1)
if i in [0, 3, 6]:
pl.ylabel('Flux (black) and log10(specmod) (blue))')
if i in [6, 7, 8]:
pl.xlabel('Freq bin')
ax.plot(spectrum[keys[i]].mean(axis=1), 'k')
ax.plot(8+np.arange(win)*16, np.log10(sms[keys[i]]), 'b*')
ax.plot(8+np.arange(win)*16, np.log10(mis[keys[i]]), 'r*')
ax.set_title(labels[keys[i]])
In [76]:
for key in keys:
print(key, mis[key][np.where(mns[key] > 0.1)].mean())