In [1]:
import copy, glob
import progressbar as PGB
import h5py
import numpy as NP
import numpy.ma as MA
from scipy import interpolate, stats
import matplotlib.pyplot as PLT
import matplotlib.colors as PLTC
import matplotlib.ticker as PLTick
import yaml, argparse, warnings
from astropy.io import ascii
import astropy.units as U
from astropy.stats import sigma_clipped_stats, histogram
from astroutils import DSP_modules as DSP
from astroutils import constants as CNST
from astroutils import mathops as OPS
from astroutils import nonmathops as NMO
from astroutils import lookup_operations as LKP
import astroutils
import prisim
from prisim import interferometry as RI
from prisim import bispectrum_phase as BSP
from prisim import delay_spectrum as DS
from IPython.core.debugger import set_trace
%matplotlib inline
In [2]:
print('AstroUtils git # {0}\nPRISim git # {1}'.format(astroutils.__githash__, prisim.__githash__))
In [44]:
inparmsfile = '/lustre/aoc/users/nthyagar/codes/mine/python/projects/closure/multiday_EQ28XX_data_Fornax_transit_closure_PS_analysis_parms.yaml'
with open(inparmsfile, 'r') as parms_file:
parms = yaml.safe_load(parms_file)
In [39]:
print(parms)
In [40]:
projectdir = parms['dirStruct']['projectdir']
datadir = projectdir + parms['dirStruct']['datadir']
figdir = datadir + parms['dirStruct']['figdir']
modelsdir = parms['dirStruct']['modeldir']
infiles = parms['dirStruct']['infiles']
visfile = parms['dirStruct']['visfile']
visfiletype = parms['dirStruct']['visfiletype']
hdf5_infile = parms['dirStruct']['hdf5_infile']
model_hdf5files = parms['dirStruct']['model_hdf5files']
model_labels = parms['dirStruct']['model_labels']
telescope_parms = parms['telescope']
site_latitude = telescope_parms['latitude']
site_longitude = telescope_parms['longitude']
preprocessinfo = parms['preProcessing']
preprocess = preprocessinfo['action']
flagchans = preprocessinfo['flagchans']
if flagchans is not None:
flagchans = NP.asarray(preprocessinfo['flagchans']).reshape(-1)
flagants = preprocessinfo['flagants']
if flagants is not None:
flagants = NP.asarray(preprocessinfo['flagants']).reshape(-1)
daybinsize = preprocessinfo['daybinsize']
ndaybins = preprocessinfo['ndaybins']
lstbinsize = preprocessinfo['lstbinsize']
band_center = preprocessinfo['band_center']
freq_resolution = preprocessinfo['freq_resolution']
dspecinfo = parms['delaySpectrum']
subbandinfo = dspecinfo['subband']
freq_window_centers = NP.asarray(subbandinfo['freq_center'])
freq_window_bw = NP.asarray(subbandinfo['bw_eff'])
freq_window_shape = subbandinfo['shape']
freq_window_fftpow = subbandinfo['fftpow']
pad = dspecinfo['pad']
apply_flags = dspecinfo['applyflags']
if apply_flags:
applyflags_str = 'Y'
else:
applyflags_str = 'N'
In [6]:
bl = NP.asarray(dspecinfo['bl'])
if bl.shape[0] != 3:
raise ValueError('Input bl must be made of three vectors forming the triad')
bltol = dspecinfo['bltol']
infile = infiles[0]
infile_no_ext = hdf5_infile.split('.hdf5')[0]
# visdata = NP.load(visfile)
if visfile is None:
visinfo = None
else:
if visfiletype == 'hdf5':
visinfo = NMO.load_dict_from_hdf5(visfile+'.hdf5')
blind, blrefind, dbl = LKP.find_1NN(visinfo['baseline']['blvect'], bl, distance_ULIM=bltol, remove_oob=True)
if blrefind.size != 3:
blind_missing = NP.setdiff1d(NP.arange(3), blind, assume_unique=True)
blind_next, blrefind_next, dbl_next = LKP.find_1NN(visinfo['baseline']['blvect'], -1*bl[blind_missing,:], distance_ULIM=bltol, remove_oob=True)
if blind_next.size + blind.size != 3:
raise ValueError('Exactly three baselines were not found in the reference baselines')
else:
blind = NP.append(blind, blind_missing[blind_next])
blrefind = NP.append(blrefind, blrefind_next)
else:
blind_missing = []
vistriad = MA.array(visinfo['vis_real'][blrefind,:,:] + 1j * visinfo['vis_imag'][blrefind,:,:], mask=visinfo['mask'][blrefind,:,:])
if len(blind_missing) > 0:
vistriad[-blrefind_next.size:,:,:] = vistriad[-blrefind_next.size:,:,:].conj()
else:
visinfo = RI.InterferometerArray(None, None, None, init_file=visfile)
In [7]:
print(visinfo.skyvis_freq.shape)
print(visinfo.lst/15.0) # LST in hours
In [8]:
tmpnpzdata = NP.load(datadir+infile)
nchan = tmpnpzdata['flags'].shape[-1]
freqs = band_center + freq_resolution * (NP.arange(nchan) - int(0.5*nchan))
In [9]:
print(tmpnpzdata['last'])
In [10]:
cpObj = BSP.ClosurePhase(datadir+hdf5_infile, freqs, infmt='hdf5')
In [11]:
print(cpObj.cpinfo['raw']['lst'].shape)
print(cpObj.cpinfo['raw']['lst'])
In [12]:
print(cpObj.cpinfo.keys())
In [13]:
print(cpObj.cpinfo['raw'].keys())
print(cpObj.cpinfo['raw']['triads'].shape)
print(cpObj.cpinfo['raw']['days'].shape)
print(cpObj.cpinfo['raw']['lst-day'].shape)
print(cpObj.cpinfo['raw']['lst'].shape)
print(cpObj.cpinfo['raw']['cphase'].shape)
In [14]:
print('day bin size = '+str(daybinsize), 'ndaybins={0:0d}'.format(ndaybins), 'LST bin size = {0:.1f}s'.format(lstbinsize))
cpObj.smooth_in_tbins(daybinsize=daybinsize, ndaybins=ndaybins, lstbinsize=lstbinsize)
In [15]:
print(cpObj.cpinfo['processed'].keys())
In [16]:
print(cpObj.cpinfo['processed']['prelim'].keys())
print(cpObj.cpinfo['processed']['prelim']['cphase'].keys())
print(cpObj.cpinfo['processed']['prelim']['lstbins'].shape)
print(cpObj.cpinfo['processed']['prelim']['lstbins'])
print(cpObj.cpinfo['processed']['prelim']['daybins'].shape)
print(cpObj.cpinfo['processed']['prelim']['daybins'])
print(cpObj.cpinfo['processed']['prelim']['cphase']['median'].shape)
print(cpObj.cpinfo['processed']['prelim']['eicp']['median'].shape)
In [17]:
cpObj.subtract(NP.zeros(1024))
print(cpObj.cpinfo['processed'].keys())
In [18]:
# ndaybins=4
print('ndaybins={0:0d}'.format(4), 'LST bin size = {0:.1f}s'.format(lstbinsize))
cpObj.subsample_differencing(daybinsize=None, ndaybins=4, lstbinsize=lstbinsize)
In [19]:
print(cpObj.cpinfo['errinfo'].keys())
print(cpObj.cpinfo['errinfo']['daybins'].shape)
print(cpObj.cpinfo['errinfo']['lstbins'].shape)
print(len(cpObj.cpinfo['errinfo']['list_of_pair_of_pairs']))
print(cpObj.cpinfo['errinfo']['list_of_pair_of_pairs'])
print(cpObj.cpinfo['errinfo']['eicp_diff'].keys())
print(cpObj.cpinfo['errinfo']['eicp_diff']['0'].keys())
print(cpObj.cpinfo['errinfo']['eicp_diff']['0']['median'].shape)
print(cpObj.cpinfo['errinfo']['eicp_diff']['1']['median'].shape)
In [20]:
cpDSobj = BSP.ClosurePhaseDelaySpectrum(cpObj)
In [21]:
if visinfo is not None:
if visfiletype == 'hdf5':
visscaleinfo = {'vis': vistriad, 'lst': visinfo['header']['LST'], 'smoothinfo': {'op_type': 'interp1d', 'interp_kind': 'linear'}}
else:
visscaleinfo = {'vis': visinfo, 'bltriplet': bl, 'smoothinfo': {'op_type': 'interp1d', 'interp_kind': 'linear'}}
else:
visscaleinfo = None
In [22]:
print('freq window centers = ', freq_window_centers)
print('freq window BW eff = ', freq_window_bw)
print('freq window shape = '+freq_window_shape)
print('freq window fftpow = {0:.1f}'.format(freq_window_fftpow))
print('pad = {0:.1f}'.format(pad))
In [23]:
cpds = cpDSobj.FT(freq_window_bw, freq_center=freq_window_centers, shape=freq_window_shape, fftpow=freq_window_fftpow, pad=pad, datapool='prelim', visscaleinfo=visscaleinfo, method='fft', resample=True, apply_flags=apply_flags)
In [47]:
plot_info = parms['plot']
plots = [key for key in plot_info if plot_info[key]['action']]
print(plots)
In [28]:
if ('1' in plots) or ('1a' in plots) or ('1b' in plots) or ('1c' in plots) or ('1d' in plots):
triads = map(tuple, cpDSobj.cPhase.cpinfo['raw']['triads'])
ntriads = len(triads)
lst = cpDSobj.cPhase.cpinfo['raw']['lst']
ntimes = lst.size
tbins = cpDSobj.cPhase.cpinfo['processed']['prelim']['lstbins']
ntbins = tbins.size
dlst = lst[1] - lst[0]
dtbins = cpDSobj.cPhase.cpinfo['processed']['prelim']['dlstbins']
flags = cpDSobj.cPhase.cpinfo['raw']['flags']
wts_raw = cpDSobj.cPhase.cpinfo['processed']['native']['wts'].data
wts_proc = cpDSobj.cPhase.cpinfo['processed']['prelim']['wts'].data
freq_wts = cpds['freq_wts']
triad = plot_info['1a']['triad']
if triad is not None:
triads_to_plot = tuple(triad)
else:
triads_to_plot = triads
print(triad)
# triad = (0, 1, 12)
# triad_ind = triads.index(triad)
In [29]:
print(triads_to_plot)
# print(triad_ind)
print(lst)
print(dlst*3600) # Difference between adjacent LST (in seconds)
print(tbins)
print(wts_raw.shape)
print(wts_raw[:,0,triad_ind,:].shape)
print(wts_raw[:,0,triad_ind,620].shape)
print(dtbins*3.6e3)
print(wts_proc.shape)
# print(wts_proc[:,0,triad_ind,:].shape)
# print(wts_proc[:,0,triad_ind,:].min(), wts_proc[:,0,triad_ind,:].max())
# print(wts_proc[:,0,triad_ind,620])
In [27]:
if '1a' in plots:
ncol = 3
nrow = NP.ceil(1.0 * len(triads_to_plot) / ncol).astype(int)
for irow in range(nrow):
fig, axs = PLT.subplots(ncols=ncol, figsize=(5*ncol,4))
for icol in range(ncol):
triadnum = irow*ncol + icol
if triadnum < len(triads_to_plot):
triplet = triads_to_plot[triadnum]
triad_ind = triads.index(triplet)
axs[icol].imshow(wts_raw[:,0,triad_ind,:], origin='lower', extent=[1e-6*freqs.min(), 1e-6*freqs.max(), lst.min(), lst.max()+NP.mean(dlst)], vmin=wts_raw.min(), vmax=wts_raw.max(), interpolation='none', cmap='gray')
axs[icol].text(0.5, 0.97, '({0[0]:0d}, {0[1]:0d}, {0[2]:0d})'.format(triplet), transform=axs[icol].transAxes, fontsize=12, weight='semibold', ha='center', va='top', color='red')
axs[icol].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
axs[icol].set_ylim(lst.min(), lst.max()+NP.mean(dlst))
axs[icol].set_aspect('auto')
axs[icol].set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium')
axs[icol].set_ylabel('LST [Hours]', fontsize=12, weight='medium')
else:
axs[icol].spines['top'].set_visible(False)
axs[icol].spines['right'].set_visible(False)
axs[icol].spines['bottom'].set_visible(False)
axs[icol].spines['left'].set_visible(False)
axs[icol].set_xticklabels([])
axs[icol].set_xticks([])
axs[icol].set_yticklabels([])
axs[icol].set_yticks([])
# PLT.savefig(figdir + '{0}_time_frequency_flags_triad_{1[0]:0d}_{1[1]:0d}_{1[2]:0d}.pdf'.format(infile_no_ext, triad), bbox_inches=0)
In [28]:
for dayind in range(ndaybins):
for irow in range(nrow):
fig, axs = PLT.subplots(ncols=ncol, figsize=(5*ncol,4))
for icol in range(ncol):
triadnum = irow*ncol + icol
if triadnum < len(triads_to_plot):
triplet = triads_to_plot[triadnum]
triad_ind = triads.index(triplet)
wtsimg = axs[icol].imshow(wts_proc[:,dayind,triad_ind,:], origin='lower', extent=[1e-6*freqs.min(), 1e-6*freqs.max(), tbins.min(), tbins.max()+NP.mean(dtbins)], vmin=wts_proc.min(), vmax=wts_proc.max(), interpolation='none', cmap='gray')
axs[icol].text(0.5, 0.97, '({0[0]:0d}, {0[1]:0d}, {0[2]:0d})'.format(triplet), transform=axs[icol].transAxes, fontsize=12, weight='semibold', ha='center', va='top', color='red')
axs[icol].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
axs[icol].set_ylim(tbins.min(), tbins.max()+NP.mean(dtbins))
axs[icol].set_aspect('auto')
axs[icol].set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium')
axs[icol].set_ylabel('LST [Hours]', fontsize=12, weight='medium')
else:
axs[icol].spines['top'].set_visible(False)
axs[icol].spines['right'].set_visible(False)
axs[icol].spines['bottom'].set_visible(False)
axs[icol].spines['left'].set_visible(False)
axs[icol].set_xticklabels([])
axs[icol].set_xticks([])
axs[icol].set_yticklabels([])
axs[icol].set_yticks([])
# fig = PLT.figure(figsize=(6,5))
# ax = fig.add_subplot(111)
# wtsimg = ax.imshow(wts_proc[:,0,triad_ind,:], origin='lower', extent=[1e-6*freqs.min(), 1e-6*freqs.max(), tbins.min(), tbins.max()+NP.mean(dtbins)], vmin=wts_proc.min(), vmax=wts_proc.max(), interpolation='none', cmap='gray')
# ax.text(0.5, 0.97, '({0[0]:0d}, {0[1]:0d}, {0[2]:0d})'.format(triad), transform=ax.transAxes, fontsize=12, weight='semibold', ha='center', va='top', color='red')
# ax.set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
# ax.set_ylim(tbins.min(), tbins.max()+NP.mean(dtbins))
# ax.set_aspect('auto')
# ax.set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium')
# ax.set_ylabel('LST [hours]', fontsize=12, weight='medium')
# cbax = fig.add_axes([0.86, 0.2, 0.02, 0.75])
# cbar = fig.colorbar(wtsimg, cax=cbax, orientation='vertical')
# cbax.yaxis.tick_right()
# # cbax.yaxis.set_label_position('right')
# fig.subplots_adjust(top=0.95)
# fig.subplots_adjust(left=0.2)
# fig.subplots_adjust(bottom=0.2)
# fig.subplots_adjust(right=0.85)
# # PLT.savefig(figdir + '{0}_time_frequency_wts_triad_{1[0]:0d}_{1[1]:0d}_{1[2]:0d}.pdf'.format(infile_no_ext, triad), bbox_inches=0)
In [29]:
if '1b' in plots:
triad = tuple(plot_info['1b']['triad'])
triad_ind = triads.index(triad)
net_wts_raw = wts_raw[:,0,triad_ind,:][NP.newaxis,:,:] * freq_wts[:,NP.newaxis,:] # nspw x nlst x nchan
net_wts_proc = wts_proc[:,0,triad_ind,:][NP.newaxis,:,:] * freq_wts[:,NP.newaxis,:] # nspw x nlst x nchan
# net_wts_raw = wts_raw[triad_ind,0,:,:][NP.newaxis,:,:] * freq_wts[:,:,NP.newaxis]
# net_wts_proc = wts_proc[triad_ind,0,:,:][NP.newaxis,:,:] * freq_wts[:,:,NP.newaxis]
nrow = freq_wts.shape[0]
fig, axs = PLT.subplots(nrows=nrow, sharex=True, sharey=True, figsize=(3.5,6))
for axind in range(len(axs)):
# wtsimg = axs[axind].imshow(net_wts_proc[axind,:,:], origin='lower', extent=[1e-6*freqs.min(), 1e-6*freqs.max(), tbins.min(), tbins.max()+NP.mean(dtbins)], norm=PLTC.LogNorm(vmin=1e-6, vmax=net_wts_proc.max()), interpolation='none', cmap='binary')
wtsimg = axs[axind].imshow(net_wts_proc[axind,:,:], origin='lower', extent=[1e-6*freqs.min(), 1e-6*freqs.max(), tbins.min(), tbins.max()+NP.mean(dtbins)], vmin=1e-6, vmax=net_wts_proc.max(), interpolation='none', cmap='binary')
if axind == 0:
axs[axind].text(0.97, 0.97, '({0[0]:0d}, {0[1]:0d}, {0[2]:0d})'.format(triad), transform=axs[axind].transAxes, fontsize=12, weight='semibold', ha='right', va='top', color='red')
axs[axind].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
axs[axind].set_ylim(tbins.min(), tbins.max()+NP.mean(dtbins))
axs[axind].set_aspect('auto')
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(top=0.95)
fig.subplots_adjust(left=0.2)
fig.subplots_adjust(bottom=0.12)
fig.subplots_adjust(right=0.85)
cbax = fig.add_axes([0.86, 0.12, 0.02, 0.3])
cbar = fig.colorbar(wtsimg, cax=cbax, orientation='vertical')
cbax.yaxis.tick_right()
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
# big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium', labelpad=20)
big_ax.set_ylabel('LST [hours]', fontsize=12, weight='medium', labelpad=35)
# PLT.savefig(figdir + '{0}_time_frequency_netwts_triad_{1[0]:0d}_{1[1]:0d}_{1[2]:0d}.pdf'.format(infile_no_ext, triad), bbox_inches=0)
In [30]:
if '1d' in plots:
datastage = plot_info['1d']['datastage']
if datastage.lower() not in ['native', 'prelim']:
raise ValueError('Input datastage value invalid')
elif datastage.lower() == 'native':
cphase = cpObj.cpinfo['processed'][datastage]['cphase']
datastr = '{0}'.format(datastage)
else:
statistic = plot_info['1d']['statistic']
cphase = cpObj.cpinfo['processed'][datastage]['cphase'][statistic]
datastr = '{0}_{1}'.format(datastage, statistic)
mask = cphase.mask
timetriad_selection = plot_info['1d']['selection']
# timetriad_selection = {'lstrange': [3.1, 3.6], 'triads': None, 'dayind': 0}
timetriad_selection = {'lstrange': [3.1, 3.6], 'triads': None, 'dayind': 18}
if timetriad_selection is not None:
dayind = timetriad_selection['dayind']
else:
dayind = 0
for key in timetriad_selection:
if timetriad_selection[key] is not None:
if key == 'triads':
triads = map(tuple, timetriad_selection[key])
elif key == 'lstrange':
lstrange = timetriad_selection[key]
if datastage.lower() == 'native':
lstbins = cpObj.cpinfo['raw']['lst'][:,dayind]
else:
lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
if lstrange is None:
lstinds = NP.arange(lstbins.size)
else:
lstrange = NP.asarray(lstrange)
lstinds = NP.where(NP.logical_and(lstbins >= lstrange.min(), lstbins <= lstrange.max()))[0]
else:
if key == 'triads':
triads = map(tuple, cpDSobj.cPhase.cpinfo['raw']['triads'])
elif key == 'lstrange':
if datastage.lower() == 'native':
lstbins = cpObj.cpinfo['raw']['lst'][:,dayind]
else:
lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
lstinds = NP.arange(lstbins.size)
sparseness = plot_info['1d']['sparseness']
if sparseness < 1.0:
sparseness = 1.0
sparsestr = '{0:.1f}'.format(sparseness)
sparsenum = NP.ceil(freqs.size / sparseness).astype(NP.int)
if sparsenum == freqs.size:
indchan = NP.arange(freqs.size)
applyflags = plot_info['1d']['applyflags']
if applyflags:
flags_str = 'flags'
else:
flags_str = 'noflags'
print(lstbins)
print(lstrange)
print(lstbins[lstinds])
print(lstinds)
print(applyflags)
print(timetriad_selection['dayind'])
print(statistic)
In [31]:
ncol = 3
nrow = min(4, int(NP.ceil(1.0*lstinds.size/ncol)))
npages = int(NP.ceil(1.0 * lstinds.size / (nrow*ncol)))
nlst_remain = lstinds.size
for pagei in range(npages):
if pagei > 0:
nlst_remain = lstinds.size - pagei * nrow * ncol
nrow = min(4, int(NP.ceil(1.0*nlst_remain/ncol)))
fig, axs = PLT.subplots(nrows=nrow, ncols=ncol, sharex=True, sharey=True, figsize=(12,10))
for i in range(nrow):
for j in range(ncol):
lstind = (lstinds.size - nlst_remain) + i*ncol+j
lind = lstinds[lstind]
if lstind < lstinds.size:
for triad in triads:
triad_ind = triads.index(triad)
if sparsenum < freqs.size:
indchan = NP.sort(NP.random.randint(freqs.size, size=sparsenum))
axs[i,j].plot(1e-6*freqs[indchan], cphase[lind,dayind,triad_ind,indchan].data, marker='.', ms=2, ls='none')
if applyflags:
flagind = mask[lind,dayind,triad_ind,:]
axs[i,j].plot(1e-6*freqs[flagind], cphase[lind,dayind,triad_ind,flagind].data, marker='.', ms=1, color='black', ls='none')
axs[i,j].text(0.5, 0.97, '{0:.2f} hrs'.format(lstbins[lind]), transform=axs[i,j].transAxes, fontsize=10, weight='medium', ha='center', va='top', color='black')
else:
axs[i,j].axis('off')
axs[i,j].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
axs[i,j].set_ylim(-3.5,3.5)
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(top=0.95)
fig.subplots_adjust(left=0.1)
fig.subplots_adjust(bottom=0.15)
fig.subplots_adjust(right=0.98)
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
# big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium', labelpad=20)
big_ax.set_ylabel(r'$\phi_\nabla$ [radians]', fontsize=12, weight='medium', labelpad=35)
# PLT.savefig(figdir + '{0}_cp_spectra_{1}_{2}_{3}_triads_day_{4}_{5:.1f}x_sparse_page_{6:03d}_of_{7:0d}.png'.format(infile_no_ext, flags_str, datastr, len(triads), dayind, sparseness, pagei+1, npages), bbox_inches=0)
# PLT.savefig(figdir + '{0}_cp_spectra_{1}_{2}_{3}_triads_day_{4}_{5:.1f}x_sparse_page_{6:03d}_of_{7:0d}.eps'.format(infile_no_ext, flags_str, datastr, len(triads), dayind, sparseness, pagei+1, npages), bbox_inches=0)
In [30]:
if '1h' in plots:
statistic = plot_info['1h']['statistic']
timetriad_selection = plot_info['1h']['selection']
cpdev = cpObj.cpinfo['processed']['prelim']['cphase'][statistic]
if timetriad_selection is not None:
dayind = timetriad_selection['dayind']
else:
dayind = 0
for key in timetriad_selection:
if timetriad_selection[key] is not None:
if key == 'triads':
triads = map(tuple, timetriad_selection[key])
elif key == 'lstrange':
lstrange = timetriad_selection[key]
lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
if lstrange is None:
lstinds = NP.arange(lstbins.size)
else:
lstrange = NP.asarray(lstrange)
lstinds = NP.where(NP.logical_and(lstbins >= lstrange.min(), lstbins <= lstrange.max()))[0]
else:
if key == 'triads':
triads = map(tuple, cpDSobj.cPhase.cpinfo['raw']['triads'])
elif key == 'lstrange':
lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
lstinds = NP.arange(lstbins.size)
print(cpdev.shape)
In [33]:
ncol = 3
nrow = min(4, int(NP.ceil(1.0*lstinds.size/ncol)))
npages = int(NP.ceil(1.0 * lstinds.size / (nrow*ncol)))
nlst_remain = lstinds.size
for pagei in range(npages):
if pagei > 0:
nlst_remain = lstinds.size - pagei * nrow * ncol
nrow = min(4, int(NP.ceil(1.0*nlst_remain/ncol)))
fig, axs = PLT.subplots(nrows=nrow, ncols=ncol, sharex=True, sharey=True, figsize=(10,8))
for i in range(nrow):
for j in range(ncol):
lstind = (lstinds.size - nlst_remain) + i*ncol+j
if lstind < lstinds.size:
lind = lstinds[lstind]
for triad in triads:
triad_ind = triads.index(triad)
axs[i,j].plot(1e-6*freqs, NP.degrees(cpdev[lind,dayind,triad_ind,:]), marker='.', ms=2, ls='none')
axs[i,j].text(0.5, 0.97, '{0:.2f} hrs'.format(lstbins[lind]), transform=axs[i,j].transAxes, fontsize=10, weight='medium', ha='center', va='top', color='black')
else:
axs[i,j].axis('off')
axs[i,j].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
axs[i,j].set_ylim(0,19)
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(top=0.95)
fig.subplots_adjust(left=0.1)
fig.subplots_adjust(bottom=0.15)
fig.subplots_adjust(right=0.98)
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
# big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium', labelpad=20)
big_ax.set_ylabel(r'$\sigma_{\phi_\nabla}$ [degrees]', fontsize=12, weight='medium', labelpad=35)
# PLT.savefig(figdir + '{0}_cp{1}_spectra_{2}_{2}_{4}_triads_day_{4}_{5:.1f}x_sparse_page_{6:03d}_of_{8:0d}.png'.format(statistic, infile_no_ext, flags_str, datastr, len(triads), dayind, sparseness, pagei+1, npages), bbox_inches=0)
# PLT.savefig(figdir + '{0}_cp{1}_spectra_{2}_{2}_{4}_triads_day_{4}_{5:.1f}x_sparse_page_{6:03d}_of_{8:0d}.eps'.format(statistic, infile_no_ext, flags_str, datastr, len(triads), dayind, sparseness, pagei+1, npages), bbox_inches=0)
In [25]:
if ('2' in plots) or ('2a' in plots) or ('2b' in plots) or ('2c' in plots) or ('2c_err' in plots) or ('2d' in plots):
dir_PS = plot_info['2']['PS_dir']
infile_pfx_a = plot_info['2']['infile_pfx_a']
outfile_pfx_a = plot_info['2']['outfile_pfx_a']
infile_pfx_b = plot_info['2']['infile_pfx_b']
outfile_pfx_b = plot_info['2']['outfile_pfx_b']
sampling = plot_info['2']['sampling']
statistic = plot_info['2']['statistic']
cohax = plot_info['2']['cohax']
incohax = plot_info['2']['incohax']
collapseax_a = plot_info['2']['collapseax_a']
collapseax_b = plot_info['2']['collapseax_b']
datapool = plot_info['2']['datapool']
pspec_unit_type = plot_info['2']['units']
ps_errtype = plot_info['2']['errtype']
nsigma = plot_info['2']['nsigma']
beaminfo = plot_info['2']['beaminfo']
xlim = plot_info['2']['xlim']
if infile_pfx_a is not None:
ps_infile_a = datadir + dir_PS + infile_pfx_a + '_' + infile_no_ext + '.hdf5'
pserr_infile_a = datadir + dir_PS + infile_pfx_a + '_' + infile_no_ext + '_errinfo.hdf5'
if outfile_pfx_a is not None:
ps_outfile_a = datadir + dir_PS + outfile_pfx_a + '_' + infile_no_ext + '.hdf5'
pserr_outfile_a = datadir + dir_PS + outfile_pfx_a + '_' + infile_no_ext + '_errinfo.hdf5'
if infile_pfx_b is not None:
ps_infile_b = datadir + dir_PS + infile_pfx_b + '_' + infile_no_ext + '.hdf5'
pserr_infile_b = datadir + dir_PS + infile_pfx_b + '_' + infile_no_ext + '_errinfo.hdf5'
if outfile_pfx_b is not None:
ps_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '.hdf5'
pserr_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '_errinfo.hdf5'
timetriad_selection = plot_info['2']['selection']
if timetriad_selection is not None:
dayind = timetriad_selection['days']
for key in timetriad_selection:
if timetriad_selection[key] is not None:
if key == 'triads':
triads = map(tuple, timetriad_selection[key])
elif key == 'lstrange':
lstrange = timetriad_selection[key]
lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
if lstrange is None:
lstinds = NP.arange(lstbins.size)
else:
lstrange = NP.asarray(lstrange)
lstinds = NP.where(NP.logical_and(lstbins >= lstrange.min(), lstbins <= lstrange.max()))[0]
if lstinds.size == 0:
raise ValueError('No data found in the specified LST range.')
else:
if key == 'triads':
triads = map(tuple, cpDSobj.cPhase.cpinfo['raw']['triads'])
elif key == 'lstrange':
lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
lstinds = NP.arange(lstbins.size)
selection = {'triads': triads, 'lst': lstinds, 'days': dayind}
autoinfo = {'axes': cohax}
xinfo_a = {'axes': incohax, 'avgcov': False, 'collapse_axes': collapseax_a, 'dlst_range': timetriad_selection['dlst_range']}
xinfo_b = {'axes': incohax, 'avgcov': False, 'collapse_axes': collapseax_b, 'dlst_range': timetriad_selection['dlst_range']}
if pspec_unit_type == 'K':
pspec_unit = 'mK2 Mpc3'
else:
pspec_unit = 'Jy2 Mpc'
subselection = plot_info['2']['subselection']
mdl_ndaybins = plot_info['2']['modelinfo']['mdl_ndaybins']
mdl_day = plot_info['2']['modelinfo']['mdl_day']
mdl_cohax = plot_info['2']['modelinfo']['mdl_cohax']
mdl_incohax = plot_info['2']['modelinfo']['mdl_incohax']
mdl_collapseax_a = plot_info['2']['modelinfo']['mdl_collapax_a']
mdl_collapseax_b = plot_info['2']['modelinfo']['mdl_collapax_b']
mdl_dir_PS = plot_info['2']['modelinfo']['PS_dir']
mdl_infile_pfx_a = plot_info['2']['modelinfo']['infile_pfx_a']
mdl_outfile_pfx_a = plot_info['2']['modelinfo']['outfile_pfx_a']
mdl_infile_pfx_b = plot_info['2']['modelinfo']['infile_pfx_b']
mdl_outfile_pfx_b = plot_info['2']['modelinfo']['outfile_pfx_b']
if model_hdf5files is not None:
mdl_autoinfo = [{'axes': mdl_cohax[i]} for i in range(len(model_hdf5files))]
mdl_xinfo_a = [{'axes': mdl_incohax[i], 'avgcov': False, 'collapse_axes': mdl_collapseax_a[i], 'dlst_range': timetriad_selection['dlst_range']} for i in range(len(model_hdf5files))]
mdl_xinfo_b = [{'axes': mdl_incohax[i], 'avgcov': False, 'collapse_axes': mdl_collapseax_b[i], 'dlst_range': timetriad_selection['dlst_range']} for i in range(len(model_hdf5files))]
if statistic is None:
statistic = ['mean', 'median']
else:
statistic = [statistic]
In [26]:
if infile_pfx_a is not None:
xcpdps2_a = BSP.read_CPhase_cross_power_spectrum(ps_infile_a)
xcpdps2_a_errinfo = BSP.read_CPhase_cross_power_spectrum(pserr_infile_a)
else:
xcpdps2_a = cpDSobj.compute_power_spectrum(selection=selection, autoinfo=autoinfo, xinfo=xinfo_a, units=pspec_unit_type, beamparms=beaminfo)
xcpdps2_a_errinfo = cpDSobj.compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_a, units=pspec_unit_type, beamparms=beaminfo)
if outfile_pfx_a is not None:
BSP.save_CPhase_cross_power_spectrum(xcpdps2_a, ps_outfile_a)
BSP.save_CPhase_cross_power_spectrum(xcpdps2_a_errinfo, pserr_outfile_a)
if infile_pfx_b is not None:
xcpdps2_b = BSP.read_CPhase_cross_power_spectrum(ps_infile_b)
xcpdps2_b_errinfo = BSP.read_CPhase_cross_power_spectrum(pserr_infile_b)
else:
xcpdps2_b = cpDSobj.compute_power_spectrum(selection=selection, autoinfo=autoinfo, xinfo=xinfo_b, units=pspec_unit_type, beamparms=beaminfo)
xcpdps2_b_errinfo = cpDSobj.compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_b, units=pspec_unit_type, beamparms=beaminfo)
if outfile_pfx_b is not None:
BSP.save_CPhase_cross_power_spectrum(xcpdps2_b, ps_outfile_b)
BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_errinfo, pserr_outfile_b)
In [27]:
nsamples_incoh = xcpdps2_a[sampling]['whole']['nsamples_incoh']
nsamples_coh = xcpdps2_a[sampling]['whole']['nsamples_coh']
In [28]:
model_cpObjs = []
model_cpDSobjs = []
cpds_models = []
xcpdps2_a_models = []
xcpdps2_a_errinfo_models = []
xcpdps2_b_models = []
xcpdps2_b_errinfo_models = []
if model_hdf5files is not None:
if mdl_infile_pfx_a is not None:
if isinstance(mdl_infile_pfx_a, list):
if (len(mdl_infile_pfx_a) > 0):
if not isinstance(mdl_dir_PS, list):
if isinstance(mdl_dir_PS, str):
mdl_dir_PS = [mdl_dir_PS] * len(model_hdf5files)
else:
raise TypeError('PS directory for models must be a list of strings')
else:
if len(mdl_dir_PS) != len(model_hdf5files):
raise ValueError('Input model PS directories must match the number of models being analyzed.')
else:
raise TypeError('Input model PS infile_a prefixes must be specified as a list of strings')
if mdl_infile_pfx_b is not None:
if isinstance(mdl_infile_pfx_b, list):
if len(mdl_infile_pfx_b) != len(mdl_infile_pfx_b):
raise ValueError('Length of input model PS infile_b prefixes must match the length of input model PS infile_a prefixes')
else:
raise TypeError('Input model PS infile_b prefixes must be specified as a list of strings')
progress = PGB.ProgressBar(widgets=[PGB.Percentage(), PGB.Bar(marker='-', left=' |', right='| '), PGB.Counter(), '/{0:0d} Models '.format(len(model_hdf5files)), PGB.ETA()], maxval=len(model_hdf5files)).start()
for i in range(len(model_hdf5files)):
mdl_infile_no_ext = model_hdf5files[i].split('.hdf5')[0]
mdl_ps_infile_a_provided = False
mdl_pserr_infile_a_provided = False
mdl_ps_infile_b_provided = False
mdl_pserr_infile_b_provided = False
if mdl_infile_pfx_a is not None:
if len(mdl_infile_pfx_a) > 0:
if mdl_infile_pfx_a[i] is not None:
if not isinstance(mdl_infile_pfx_a[i], str):
raise TypeError('Input {0}-th model cross PS file must be a string'.format(i+1))
else:
try:
model_xcpdps2_a = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'_'+mdl_infile_no_ext+'.hdf5')
except IOError as xcption:
mdl_ps_infile_a_provided = False
warnings.warn('Provided model cross-power spectrum infile_a "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'.hdf5'))
else:
mdl_ps_infile_a_provided = True
xcpdps2_a_models += [copy.deepcopy(model_xcpdps2_a)]
try:
model_xcpdps2_a_errinfo = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'_'+mdl_infile_no_ext+'_errinfo.hdf5')
except IOError as xcption:
mdl_pserr_infile_a_provided = False
warnings.warn('Provided model cross-power spectrum infile_a "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'_errinfo.hdf5'))
else:
mdl_pserr_infile_a_provided = True
xcpdps2_a_errinfo_models += [copy.deepcopy(model_xcpdps2_a_errinfo)]
if mdl_infile_pfx_b is not None:
if len(mdl_infile_pfx_b) > 0:
if mdl_infile_pfx_b[i] is not None:
if not isinstance(mdl_infile_pfx_b[i], str):
raise TypeError('Input {0}-th model cross PS file must be a string'.format(i+1))
else:
try:
model_xcpdps2_b = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'_'+mdl_infile_no_ext+'.hdf5')
except IOError as xcption:
mdl_ps_infile_b_provided = False
warnings.warn('Provided model cross-power spectrum infile_b "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'.hdf5'))
else:
mdl_ps_infile_b_provided = True
xcpdps2_b_models += [copy.deepcopy(model_xcpdps2_b)]
try:
model_xcpdps2_b_errinfo = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'_'+mdl_infile_no_ext+'_errinfo.hdf5')
except IOError as xcption:
mdl_pserr_infile_b_provided = False
warnings.warn('Provided model cross-power spectrum infile_b "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'_errinfo.hdf5'))
else:
mdl_pserr_infile_b_provided = True
xcpdps2_b_errinfo_models += [copy.deepcopy(model_xcpdps2_b_errinfo)]
if (not mdl_ps_infile_a_provided) or (not mdl_pserr_infile_a_provided) or (not mdl_ps_infile_b_provided) or (not mdl_pserr_infile_b_provided):
model_cpObj = BSP.ClosurePhase(modelsdir+model_hdf5files[i], freqs, infmt='hdf5')
model_cpObj.smooth_in_tbins(daybinsize=daybinsize, ndaybins=mdl_ndaybins[i], lstbinsize=lstbinsize)
model_cpObj.subsample_differencing(daybinsize=None, ndaybins=4, lstbinsize=lstbinsize)
model_cpObj.subtract(NP.zeros(1024))
model_cpObjs += [copy.deepcopy(model_cpObj)]
model_cpDSobjs += [BSP.ClosurePhaseDelaySpectrum(model_cpObjs[i])]
cpds_models += [model_cpDSobjs[i].FT(freq_window_bw, freq_center=freq_window_centers, shape=freq_window_shape, fftpow=freq_window_fftpow, pad=pad, datapool='prelim', visscaleinfo=visscaleinfo, method='fft', resample=True, apply_flags=apply_flags)]
if not mdl_ps_infile_a_provided:
xcpdps2_a_models += [model_cpDSobjs[i].compute_power_spectrum(selection=selection, autoinfo=mdl_autoinfo[i], xinfo=mdl_xinfo_a[i], units=pspec_unit_type, beamparms=beaminfo)]
if not mdl_pserr_infile_a_provided:
xcpdps2_a_errinfo_models += [model_cpDSobjs[i].compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_a, units=pspec_unit_type, beamparms=beaminfo)]
if not mdl_ps_infile_b_provided:
xcpdps2_b_models += [model_cpDSobjs[i].compute_power_spectrum(selection=selection, autoinfo=mdl_autoinfo[i], xinfo=mdl_xinfo_b[i], units=pspec_unit_type, beamparms=beaminfo)]
if not mdl_pserr_infile_b_provided:
xcpdps2_b_errinfo_models += [model_cpDSobjs[i].compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_b, units=pspec_unit_type, beamparms=beaminfo)]
else:
model_cpObjs += [None]
model_cpDSobjs += [None]
cpds_models += [None]
if mdl_outfile_pfx_a is not None:
if isinstance(mdl_outfile_pfx_a, str):
mdl_outfile_pfx_a = [mdl_outfile_pfx_a] * len(model_hdf5files)
if not isinstance(mdl_outfile_pfx_a, list):
raise TypeError('The model cross-power spectrum outfile prefixes must be specified as a list with item for each model.')
if len(mdl_outfile_pfx_a) != len(mdl_dir_PS):
raise ValueError('Invalid number of model cross-power output files specified')
mdl_ps_outfile_a = mdl_dir_PS[i] + mdl_outfile_pfx_a[i] + '_' + mdl_infile_no_ext + '.hdf5'
mdl_pserr_outfile_a = mdl_dir_PS[i] + mdl_outfile_pfx_a[i] + '_' + mdl_infile_no_ext + '_errinfo.hdf5'
BSP.save_CPhase_cross_power_spectrum(xcpdps2_a_models[-1], mdl_ps_outfile_a)
BSP.save_CPhase_cross_power_spectrum(xcpdps2_a_errinfo_models[-1], mdl_pserr_outfile_a)
if mdl_outfile_pfx_b is not None:
if isinstance(mdl_outfile_pfx_b, str):
mdl_outfile_pfx_b = [mdl_outfile_pfx_b] * len(model_hdf5files)
if not isinstance(mdl_outfile_pfx_b, list):
raise TypeError('The model cross-power spectrum outfile prefixes must be specified as a list with item for each model.')
if len(mdl_outfile_pfx_b) != len(mdl_dir_PS):
raise ValueError('Invalid number of model cross-power output files specified')
mdl_ps_outfile_b = mdl_dir_PS[i] + mdl_outfile_pfx_b[i] + '_' + mdl_infile_no_ext + '.hdf5'
mdl_pserr_outfile_b = mdl_dir_PS[i] + mdl_outfile_pfx_b[i] + '_' + mdl_infile_no_ext + '_errinfo.hdf5'
BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_models[-1], mdl_ps_outfile_b)
BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_errinfo_models[-1], mdl_pserr_outfile_b)
progress.update(i+1)
progress.finish()
In [29]:
spw = subselection['spw']
if spw is None:
spwind = NP.arange(xcpdps2_a[sampling]['z'].size)
else:
spwind = NP.asarray(spw)
lstind = NMO.find_list_in_list(xcpdps2_a[sampling][datapool[0]]['diagoffsets'][1], NP.asarray(subselection['lstdiag']))
dayind = NP.asarray(subselection['day'])
dayind_models = NP.asarray(mdl_day)
triadind = NMO.find_list_in_list(xcpdps2_a[sampling][datapool[0]]['diagoffsets'][3], NP.asarray(subselection['triaddiag']))
errshade = {}
for errtype in ps_errtype:
if errtype.lower() == 'ssdiff':
errshade[errtype] = '0.8'
elif errtype.lower() == 'psdiff':
errshade[errtype] = '0.6'
mdl_colrs = ['red', 'green', 'blue', 'cyan', 'gray', 'orange']
In [30]:
if '2b' in plots:
for stat in statistic:
for zind in spwind:
for lind in lstind:
for di,dind in enumerate(dayind):
maxabsvals = []
minabsvals = []
maxvals = []
minvals = []
fig, axs = PLT.subplots(nrows=1, ncols=len(datapool), sharex=True, sharey=True, figsize=(4.0*len(datapool), 3.6))
if len(datapool) == 1:
axs = [axs]
for dpoolind,dpool in enumerate(datapool):
for trno,trind in enumerate([triadind[0]]):
if model_hdf5files is not None:
for mdlind, mdl in enumerate(model_labels):
if dpool in xcpdps2_a_models[mdlind][sampling]:
psval = (2/3.0) * xcpdps2_a_models[mdlind][sampling][dpool][stat][zind,lind,dayind_models[di][mdlind][0],dayind_models[di][mdlind][1],trind,:].to(pspec_unit).value
# negind = psval.real < 0.0
# posind = NP.logical_not(negind)
maxabsvals += [NP.abs(psval.real).max()]
minabsvals += [NP.abs(psval.real).min()]
maxvals += [psval.real.max()]
minvals += [psval.real.min()]
axs[dpoolind].plot(xcpdps2_a_models[mdlind][sampling]['kprll'][zind,:], psval.real, ls='none', marker='.', ms=3, color=mdl_colrs[mdlind], label='{0}'.format(mdl))
if dpool in xcpdps2_a[sampling]:
psval = (2/3.0) * xcpdps2_a[sampling][dpool][stat][zind,lind,dind[0],dind[1],trind,:].to(pspec_unit).value
psrms = (2/3.0) * NP.nanstd(xcpdps2_a_errinfo[sampling]['errinfo'][stat][zind,lind,:,trind,:], axis=0).to(pspec_unit).value
maxabsvals += [NP.abs(psval.real + psrms).max()]
minabsvals += [NP.abs(psval.real).min()]
maxvals += [(psval.real + psrms).max()]
minvals += [(psval.real - psrms).min()]
# axs[dpoolind].plot(xcpdps2_a[sampling]['kprll'][zind,:], psval.real, ls='none', marker='.', ms=1, color='black', label='FG+N')
axs[dpoolind].errorbar(xcpdps2_a[sampling]['kprll'][zind,:], psval.real, yerr=psrms, xerr=None, ecolor='0.8', ls='none', marker='.', ms=4, color='black', label='FG+N')
legend = axs[dpoolind].legend(loc='center', bbox_to_anchor=(0.5,0.3), shadow=False, fontsize=8)
if trno == 0:
axs[dpoolind].text(0.05, 0.97, 'Real', transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
axs[dpoolind].text(0.95, 0.97, r'$z=$'+' {0:.1f}'.format(xcpdps2_a[sampling]['z'][zind]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='right', va='top', color='black')
axs[dpoolind].text(0.05, 0.92, r'$\Delta$'+'LST = {0:.1f} s'.format(lind*3.6e3*xcpdps2_a['dlst'][0]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
axs[dpoolind].text(0.05, 0.87, 'G{0[0]:0d}{0[1]:0d}'.format(dind), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
axt = axs[dpoolind].twiny()
axt.set_xlim(1e6*xcpdps2_a[sampling]['lags'].min(), 1e6*xcpdps2_a[sampling]['lags'].max())
minvals = NP.asarray(minvals)
maxvals = NP.asarray(maxvals)
minabsvals = NP.asarray(minabsvals)
maxabsvals = NP.asarray(maxabsvals)
if xlim is None:
axs[dpoolind].set_xlim(0.99*xcpdps2_a[sampling]['kprll'][zind,:].min(), 1.01*xcpdps2_a[sampling]['kprll'][zind,:].max())
else:
axs[dpoolind].set_xlim(xlim)
if NP.min(minvals) < 0.0:
axs[dpoolind].set_ylim(1.5*NP.min(minvals), 2*NP.max(maxabsvals))
else:
axs[dpoolind].set_ylim(0.5*NP.min(minvals), 2*NP.max(maxabsvals))
axs[dpoolind].set_yscale('symlog', linthreshy=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))))
tickloc = PLTick.SymmetricalLogLocator(linthresh=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))), base=100.0)
axs[dpoolind].yaxis.set_major_locator(tickloc)
axs[dpoolind].grid(color='0.8', which='both', linestyle=':', lw=1)
fig.subplots_adjust(top=0.85)
fig.subplots_adjust(bottom=0.16)
fig.subplots_adjust(left=0.22)
fig.subplots_adjust(right=0.98)
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
# big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$\kappa_\parallel$'+' [pseudo '+r'$h$'+' Mpc'+r'$^{-1}$'+']', fontsize=12, weight='medium', labelpad=20)
if pspec_unit_type == 'K':
big_ax.set_ylabel(r'$\frac{2}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo mK$^2h^{-3}$ Mpc$^3$]', fontsize=12, weight='medium', labelpad=35)
else:
big_ax.set_ylabel(r'$\frac{2}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo Jy$^2h^{-1}$ Mpc]', fontsize=12, weight='medium', labelpad=35)
big_axt = big_ax.twiny()
big_axt.set_xticks([])
big_axt.set_xlabel(r'$\tau$'+' ['+r'$\mu$'+'s]', fontsize=12, weight='medium', labelpad=20)
PLT.savefig(figdir + '{0}_symlog_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_lstdiag_{5:0d}_day_{6[0]:0d}_{6[1]:0d}_triaddiags_{7:0d}_flags_{8}.pdf'.format(infile_no_ext, xcpdps2_a[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_a['dlst'][0], subselection['lstdiag'][lind], dind, xcpdps2_a[sampling][datapool[0]]['diagoffsets'][3][trind], applyflags_str), bbox_inches=0)
print(figdir + '{0}_symlog_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_lstdiag_{5:0d}_day_{6[0]:0d}_{6[1]:0d}_triaddiags_{7:0d}_flags_{8}.pdf'.format(infile_no_ext, xcpdps2_a[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_a['dlst'][0], subselection['lstdiag'][lind], dind, xcpdps2_a[sampling][datapool[0]]['diagoffsets'][3][trind], applyflags_str))
In [31]:
if ('2c' in plots) or ('2c_err' in plots) or ('2d' in plots):
kprll_min_for_rms = plot_info['2c']['kprll_min']
if kprll_min_for_rms is None:
kprll_min_for_rms = 0.0
avg_incohax_a = plot_info['2c']['incohax_a']
diagoffsets_incohax_a = plot_info['2c']['diagoffsets_a']
diagoffsets_a = []
avg_incohax_b = plot_info['2c']['incohax_b']
diagoffsets_incohax_b = plot_info['2c']['diagoffsets_b']
diagoffsets_b = []
for combi,incax_comb in enumerate(avg_incohax_a):
diagoffsets_a += [{}]
for incaxind,incax in enumerate(incax_comb):
diagoffsets_a[-1][incax] = NP.asarray(diagoffsets_incohax_a[combi][incaxind])
xcpdps2_a_avg, excpdps2_a_avg = BSP.incoherent_cross_power_spectrum_average(xcpdps2_a, excpdps=xcpdps2_a_errinfo, diagoffsets=diagoffsets_a)
avg_xcpdps2_a_models = []
avg_excpdps2_a_models = []
for combi,incax_comb in enumerate(avg_incohax_b):
diagoffsets_b += [{}]
for incaxind,incax in enumerate(incax_comb):
diagoffsets_b[-1][incax] = NP.asarray(diagoffsets_incohax_b[combi][incaxind])
# xcpdps2_b_avg, excpdps2_b_avg = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b, excpdps=None, diagoffsets=diagoffsets_b)
xcpdps2_b_avg, excpdps2_b_avg = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b, excpdps=xcpdps2_b_errinfo, diagoffsets=diagoffsets_b)
avg_xcpdps2_b_models = []
avg_excpdps2_b_models = []
if model_hdf5files is not None:
progress = PGB.ProgressBar(widgets=[PGB.Percentage(), PGB.Bar(marker='-', left=' |', right='| '), PGB.Counter(), '/{0:0d} Models '.format(len(model_hdf5files)), PGB.ETA()], maxval=len(model_hdf5files)).start()
for i in range(len(model_hdf5files)):
avg_xcpdps2_a_model, avg_excpdps2_a_model = BSP.incoherent_cross_power_spectrum_average(xcpdps2_a_models[i], excpdps=xcpdps2_a_errinfo_models[i], diagoffsets=diagoffsets_a)
avg_xcpdps2_a_models += [copy.deepcopy(avg_xcpdps2_a_model)]
avg_excpdps2_a_models += [copy.deepcopy(avg_excpdps2_a_model)]
# avg_xcpdps2_b_model, avg_excpdps2_b_model = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b_models[i], excpdps=None, diagoffsets=diagoffsets_b)
avg_xcpdps2_b_model, avg_excpdps2_b_model = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b_models[i], excpdps=xcpdps2_b_errinfo_models[i], diagoffsets=diagoffsets_b)
avg_xcpdps2_b_models += [copy.deepcopy(avg_xcpdps2_b_model)]
avg_excpdps2_b_models += [copy.deepcopy(avg_excpdps2_b_model)]
progress.update(i+1)
progress.finish()
In [32]:
# Save incoherent cross power average of the main dataset and its uncertainties
xps_avg_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '_incoh_diag_avg.hdf5'
xpserr_avg_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '_errinfo_incoh_diag_avg.hdf5'
BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_avg, xps_avg_outfile_b)
BSP.save_CPhase_cross_power_spectrum(excpdps2_b_avg, xpserr_avg_outfile_b)
In [33]:
print(diagoffsets_b[0])
print(diagoffsets_b[1])
print(diagoffsets_b[2])
print(diagoffsets_b[3])
print(len(xcpdps2_b_avg['resampled']['whole']['median']))
In [35]:
if '2c' in plots:
# ylim = {0: [-1e4, 1e10], 1: [-1e4, 1e9]} # keys are spw indices
# linthreshy = 1e0
ylim = {0: [-1e10, 1e16], 1: [-1e10, 1e15]} # keys are spw indices
linthreshy = 9e5
lstind = [0]
triadind = [0]
dayind = [0]
dayind_models = NP.zeros(len(model_labels), dtype=int).reshape(1,-1)
for stat in statistic:
for zind in spwind:
kprll_ind_for_rms = NP.where(NP.abs(excpdps2_b_avg[sampling]['kprll'][zind,:]) >= kprll_min_for_rms)[0]
kprll_for_rms = excpdps2_b_avg[sampling]['kprll'][zind,kprll_ind_for_rms]
for lind in lstind:
for di,dind in enumerate(dayind):
for combi in range(len(diagoffsets_b)):
maxabsvals = []
minabsvals = []
maxvals = []
minvals = []
fig, axs = PLT.subplots(nrows=1, ncols=len(datapool), sharex=True, sharey=True, figsize=(4.0*len(datapool), 3.6))
if len(datapool) == 1:
axs = [axs]
for dpoolind,dpool in enumerate(datapool):
for trno,trind in enumerate(triadind):
# if model_hdf5files is not None:
# for mdlind, mdl in enumerate(model_labels):
# if dpool in avg_xcpdps2_b_models[mdlind][sampling]:
# psval = (1/3.0) * avg_xcpdps2_b_models[mdlind][sampling][dpool][stat][combi][zind,lind,dayind_models[di][mdlind],trind,:].to(pspec_unit).value
# maxabsvals += [NP.abs(psval.real).max()]
# minabsvals += [NP.abs(psval.real).min()]
# maxvals += [psval.real.max()]
# minvals += [psval.real.min()]
# axs[dpoolind].plot(avg_xcpdps2_b_models[mdlind][sampling]['kprll'][zind,:], psval.real, ls='none', marker='.', ms=3, color=mdl_colrs[mdlind], label='{0}'.format(mdl))
if dpool in xcpdps2_b_avg[sampling]:
psval = (2/3.0) * xcpdps2_b_avg[sampling][dpool][stat][combi][zind,lind,dind,trind,:].to(pspec_unit).value
ps_ssdiff = (2/3.0) * excpdps2_b_avg[sampling]['errinfo'][stat][combi][zind,lind,:,trind,kprll_ind_for_rms].to(pspec_unit).value # Single RMS across all k_prll bins
psrms_ssdiff = (2/3.0) * NP.nanstd(excpdps2_a_avg[sampling]['errinfo'][stat][combi][zind,lind,:,trind,:], axis=0).to(pspec_unit).value # RMS per k_prll bin
if 2 in avg_incohax_b[combi]:
ind_dayax_in_incohax = avg_incohax_b[combi].index(2)
if 0 in diagoffsets_incohax_b[combi][ind_dayax_in_incohax]:
rms_inflation_factor = 2.0 * NP.sqrt(2.0)
else:
rms_inflation_factor = NP.sqrt(2.0)
else:
rms_inflation_factor = NP.sqrt(2.0)
psrms_psdiff = (2/3.0) * (xcpdps2_a_avg[sampling][dpool][stat][combi][zind,lind,1,1,trind,:] - xcpdps2_a_avg[sampling][dpool][stat][combi][zind,lind,0,0,trind,:]).to(pspec_unit).value
psrms_psdiff = NP.abs(psrms_psdiff.real) / rms_inflation_factor
psrms_max = NP.amax(NP.vstack((psrms_ssdiff, psrms_psdiff)), axis=0)
maxabsvals += [NP.abs(psval.real + nsigma*psrms_max).max()]
minabsvals += [NP.abs(psval.real).min()]
maxvals += [(psval.real + nsigma*psrms_max).max()]
minvals += [(psval.real - nsigma*psrms_max).min()]
for errtype in ps_errtype:
if errtype.lower() == 'ssdiff':
axs[dpoolind].errorbar(xcpdps2_b_avg[sampling]['kprll'][zind,:], psval.real, yerr=nsigma*psrms_ssdiff, xerr=None, ecolor=errshade[errtype.lower()], ls='none', marker='.', ms=4, color='black')
psmean_ssdiff_across_kbins, psmedian_ssdiff_across_kbins, psrms_ssdiff_across_kbins = sigma_clipped_stats(NP.vstack((ps_ssdiff.real, ps_ssdiff.imag)))
pos_kprll_ind_for_rms = NP.where(kprll_for_rms >= 0.0)[0]
neg_kprll_ind_for_rms = NP.where(kprll_for_rms <= 0.0)[0]
axs[dpoolind].hlines([-psrms_ssdiff_across_kbins, psrms_ssdiff_across_kbins], kprll_for_rms[pos_kprll_ind_for_rms].min(), kprll_for_rms[pos_kprll_ind_for_rms].max(), linestyles='--', lw=1.5, color='0.25')
axs[dpoolind].hlines([-psrms_ssdiff_across_kbins, psrms_ssdiff_across_kbins], kprll_for_rms[neg_kprll_ind_for_rms].min(), kprll_for_rms[neg_kprll_ind_for_rms].max(), linestyles='--', lw=1.5, color='0.25')
elif errtype.lower() == 'psdiff':
axs[dpoolind].errorbar(xcpdps2_b_avg[sampling]['kprll'][zind,:], psval.real, yerr=nsigma*psrms_psdiff, xerr=None, ecolor=errshade[errtype.lower()], ls='none', marker='.', ms=4, color='black', label='FG+N')
# legend = axs[dpoolind].legend(loc='center', bbox_to_anchor=(0.5,0.3), shadow=False, fontsize=8)
if trno == 0:
# axs[dpoolind].text(0.05, 0.97, 'Real', transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
axs[dpoolind].text(0.95, 0.97, r'$z=$'+' {0:.1f}'.format(xcpdps2_b_avg[sampling]['z'][zind]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='right', va='top', color='black')
# axs[dpoolind].text(0.05, 0.92, r'$\Delta$'+'LST = {0:.1f} s'.format(lind*3.6e3*xcpdps2_a_avg['dlst'][0]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
# axs[dpoolind].text(0.05, 0.87, 'G{0[0]:0d}{0[1]:0d}'.format(dind), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
axt = axs[dpoolind].twiny()
axt.set_xlim(1e6*xcpdps2_b_avg[sampling]['lags'].min(), 1e6*xcpdps2_b_avg[sampling]['lags'].max())
axs[dpoolind].axhline(y=0, xmin=0, xmax=1, ls='-', lw=1, color='black')
minvals = NP.asarray(minvals)
maxvals = NP.asarray(maxvals)
minabsvals = NP.asarray(minabsvals)
maxabsvals = NP.asarray(maxabsvals)
if xlim is None:
axs[dpoolind].set_xlim(0.99*xcpdps2_b_avg[sampling]['kprll'][zind,:].min(), 1.01*xcpdps2_b_avg[sampling]['kprll'][zind,:].max())
else:
axs[dpoolind].set_xlim(xlim)
# if NP.min(minvals) < 0.0:
# axs[dpoolind].set_ylim(1.5*NP.min(minvals), 2*NP.max(maxabsvals))
# else:
# axs[dpoolind].set_ylim(0.5*NP.min(minvals), 2*NP.max(maxabsvals))
axs[dpoolind].set_ylim(ylim[zind])
# axs[dpoolind].set_yscale('symlog', linthreshy=10**NP.ceil(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))))
# tickloc = PLTick.SymmetricalLogLocator(linthresh=10**NP.ceil(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))), base=100.0)
axs[dpoolind].set_yscale('symlog', linthreshy=linthreshy)
tickloc = PLTick.SymmetricalLogLocator(linthresh=linthreshy, base=100.0)
axs[dpoolind].yaxis.set_major_locator(tickloc)
yticklocs = NP.asarray(axs[dpoolind].get_yticks())
ytickspacings = NP.diff(yticklocs)
tickinds_to_remove = NP.setdiff1d(NP.where(NP.abs(yticklocs)<=linthreshy)[0], NP.argmin(NP.abs(yticklocs)))
if tickinds_to_remove.size > 0:
new_yticklocs = yticklocs[NP.setdiff1d(NP.arange(yticklocs.size), tickinds_to_remove)]
axs[dpoolind].yaxis.set_major_locator(PLTick.FixedLocator((new_yticklocs.tolist())))
axs[dpoolind].grid(color='0.8', which='both', linestyle=':', lw=1)
fig.subplots_adjust(top=0.85)
fig.subplots_adjust(bottom=0.16)
fig.subplots_adjust(left=0.22)
fig.subplots_adjust(right=0.98)
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
# big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$\kappa_\parallel$'+' [pseudo '+r'$h$'+' Mpc'+r'$^{-1}$'+']', fontsize=12, weight='medium', labelpad=20)
if pspec_unit_type == 'K':
big_ax.set_ylabel(r'$\frac{2}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo mK$^2h^{-3}$ Mpc$^3$]', fontsize=12, weight='medium', labelpad=35)
else:
big_ax.set_ylabel(r'$\frac{2}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo Jy$^2h^{-1}$ Mpc]', fontsize=12, weight='medium', labelpad=35)
big_axt = big_ax.twiny()
big_axt.set_xticks([])
big_axt.set_xlabel(r'$\tau$'+' ['+r'$\mu$'+'s]', fontsize=12, weight='medium', labelpad=20)
PLT.savefig(figdir + '{0}_symlog_incoh_avg_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_b_avg[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg['dlst'][0], applyflags_str, combi), bbox_inches=0)
print(figdir + '{0}_symlog_incoh_avg_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_b_avg[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg['dlst'][0], applyflags_str, combi))
In [36]:
print(dayind)
In [48]:
if '2c_err' in plots:
current_label = plot_info['2c_err']['current_lbl']
other_label = plot_info['2c_err']['other_label']
other_pserr_file = plot_info['2c_err']['other_pserr_file']
other_excpdps2_b_avg = BSP.read_CPhase_cross_power_spectrum(other_pserr_file)
lstind = [0]
triadind = [0]
dayind = [0]
for stat in statistic:
print(stat)
print('=======================================')
for zind in spwind:
kprll_ind_current_for_rms = NP.where(NP.abs(excpdps2_b_avg[sampling]['kprll'][zind,:]) >= kprll_min_for_rms)[0]
kprll_current_for_rms = excpdps2_b_avg[sampling]['kprll'][zind,kprll_ind_current_for_rms]
kprll_ind_other_for_rms = NP.where(NP.abs(other_excpdps2_b_avg[sampling]['kprll'][zind,:]) >= kprll_min_for_rms)[0]
kprll_other_for_rms = other_excpdps2_b_avg[sampling]['kprll'][zind,kprll_ind_other_for_rms]
for lind in lstind:
for dpoolind,dpool in enumerate(datapool):
for di,dind in enumerate(dayind):
for trno,trind in enumerate(triadind):
for combi in range(len(diagoffsets_b)):
ps_ssdiff_current = (2/3.0) * excpdps2_b_avg[sampling]['errinfo'][stat][combi][zind,lind,:,trind,kprll_ind_current_for_rms].to(pspec_unit).value
ps_ssdiff_other = (2/3.0) * other_excpdps2_b_avg[sampling]['errinfo'][stat][combi][zind,lind,:,trind,kprll_ind_other_for_rms].to(pspec_unit).value
hist_current_real, hist_current_real_bins = histogram(ps_ssdiff_current.real.ravel(), bins='knuth', density=False)
hist_other_real, hist_other_real_bins = histogram(ps_ssdiff_other.real.ravel(), bins='knuth', density=False)
hist_current_imag, hist_current_imag_bins = histogram(ps_ssdiff_current.imag.ravel(), bins='knuth', density=False)
hist_other_imag, hist_other_imag_bins = histogram(ps_ssdiff_other.imag.ravel(), bins='knuth', density=False)
hist_current, hist_current_bins = histogram(NP.vstack((ps_ssdiff_current.real,ps_ssdiff_current.imag)).ravel(), bins='knuth', density=False)
hist_other, hist_other_bins = histogram(NP.vstack((ps_ssdiff_other.real,ps_ssdiff_other.imag)).ravel(), bins='knuth', density=False)
current_real_mean, current_real_median, current_real_std = sigma_clipped_stats(ps_ssdiff_current.real)
current_imag_mean, current_imag_median, current_imag_std = sigma_clipped_stats(ps_ssdiff_current.imag)
current_mean, current_median, current_std = sigma_clipped_stats(NP.vstack((ps_ssdiff_current.real, ps_ssdiff_current.imag)))
other_real_mean, other_real_median, other_real_std = sigma_clipped_stats(ps_ssdiff_other.real)
other_imag_mean, other_imag_median, other_imag_std = sigma_clipped_stats(ps_ssdiff_other.imag)
other_mean, other_median, other_std = sigma_clipped_stats(NP.vstack((ps_ssdiff_other.real, ps_ssdiff_other.imag)))
real_kval, real_pval = stats.ks_2samp(ps_ssdiff_current.real.ravel(), ps_ssdiff_other.real.ravel())
imag_kval, imag_pval = stats.ks_2samp(ps_ssdiff_current.imag.ravel(), ps_ssdiff_other.imag.ravel())
kval, pval = stats.ks_2samp(NP.vstack((ps_ssdiff_other.real, ps_ssdiff_other.imag)).ravel(), NP.vstack((ps_ssdiff_current.real, ps_ssdiff_current.imag)).ravel())
print('===================')
print('SpW: {0}, Diagcomb: {1}'.format(zind, combi))
print('-------------------')
print('Current (real): Mean = {0},\t Median = {1},\t RMS = {2}'.format(current_real_mean, current_real_median, current_real_std))
print('Other (real): Mean = {0},\t Median = {1},\t RMS = {2}'.format(other_real_mean, other_real_median, other_real_std))
print('K-S (real): k-val = {0},\t p-val = {1}'.format(real_kval, real_pval))
print('\nCurrent (imag): Mean = {0},\t Median = {1},\t RMS = {2}'.format(current_imag_mean, current_imag_median, current_imag_std))
print('Other (imag): Mean = {0},\t Median = {1},\t RMS = {2}'.format(other_imag_mean, other_imag_median, other_imag_std))
print('K-S (imag): k-val = {0},\t p-val = {1}'.format(imag_kval, imag_pval))
print('\nCurrent: Mean = {0},\t Median = {1},\t RMS = {2}'.format(current_mean, current_median, current_std))
print('Other: Mean = {0},\t Median = {1},\t RMS = {2}'.format(other_mean, other_median, other_std))
print('K-S: k-val = {0},\t p-val = {1}'.format(kval, pval))
fig, axs = PLT.subplots(nrows=3, sharex=True, sharey=True, figsize=(4.0, 4.5))
hist_current_real, hist_current_real_bins, patches = axs[0].hist(ps_ssdiff_current.real.ravel(), bins=hist_current_real_bins, density=False, histtype='step', ls='-', color='black', lw=2, label='{0} (real)'.format(current_label))
hist_other_real, hist_other_real_bins, patches = axs[0].hist(ps_ssdiff_other.real.ravel(), bins=hist_other_real_bins, density=False, histtype='step', ls='-', color='gray', lw=2, label='{0} (real)'.format(other_label))
hist_current_imag, hist_current_imag_bins, patches = axs[1].hist(ps_ssdiff_current.imag.ravel(), bins=hist_current_imag_bins, density=False, histtype='step', ls='-', color='black', lw=2, label='{0} (imag)'.format(current_label))
hist_other_imag, hist_other_imag_bins, patches = axs[1].hist(ps_ssdiff_other.imag.ravel(), bins=hist_other_imag_bins, density=False, histtype='step', ls='-', color='gray', lw=2, label='{0} (imag)'.format(other_label))
hist_current, hist_current_bins, patches = axs[2].hist(NP.vstack((ps_ssdiff_current.real,ps_ssdiff_current.imag)).ravel(), bins=hist_current_bins, density=False, histtype='step', ls='-', color='black', lw=2, label='{0}'.format(current_label))
hist_other, hist_other_bins, patches = axs[2].hist(NP.vstack((ps_ssdiff_other.real,ps_ssdiff_other.imag)).ravel(), bins=hist_current_bins, density=False, histtype='step', ls='-', color='gray', lw=2, label='{0}'.format(other_label))
axs[0].legend(loc='upper right')
axs[1].legend(loc='upper right')
axs[2].legend(loc='upper right')
axs[2].set_xlim(-5*other_std, 5*other_std)
fig.subplots_adjust(hspace=0, wspace=0)
In [37]:
if '2d' in plots:
kbin_min = plot_info['2d']['kbin_min']
kbin_max = plot_info['2d']['kbin_max']
num_kbins = plot_info['2d']['num_kbins']
kbintype = plot_info['2d']['kbintype']
if (kbin_min is None) or (kbin_max is None):
kbins = None
else:
if num_kbins is None:
raise ValueError('Input num_kbins must be set if kbin range is provided')
if kbintype == 'linear':
kbins = NP.linspace(kbin_min, kbin_max, num=num_kbins, endpoint=True)
elif kbintype == 'log':
if kbin_min > 0.0:
kbins = NP.geomspace(kbin_min, kbin_max, num=num_kbins, endpoint=True)
elif kbin_min == 0.0:
eps_k = 1e-3
kbins = NP.geomspace(kbin_min+eps_k, kbin_max, num=num_kbins, endpoint=True)
else:
eps_k = 1e-3
kbins_pos = NP.geomspace(eps_k, kbin_max, num=num_kbins, endpoint=True)
ind_kbin_thresh = NP.argmin(kbins_pos[kbins_pos >= NP.abs(kbin_min)])
kbins_neg = -1 * kbins_pos[:ind_kbin_thresh+1][::-1]
kbins = NP.hstack((kbins_neg, kbins_pos))
else:
raise ValueError('Input kbintype must be set to "linear" or "log"')
xcpdps2_a_avg_kbin = BSP.incoherent_kbin_averaging(xcpdps2_a_avg, kbins=kbins, kbintype=kbintype)
excpdps2_a_avg_kbin = BSP.incoherent_kbin_averaging(excpdps2_a_avg, kbins=kbins, kbintype=kbintype)
xcpdps2_a_avg_kbin_models = []
excpdps2_a_avg_kbin_models = []
xcpdps2_b_avg_kbin = BSP.incoherent_kbin_averaging(xcpdps2_b_avg, kbins=kbins, kbintype=kbintype)
excpdps2_b_avg_kbin = BSP.incoherent_kbin_averaging(excpdps2_b_avg, kbins=kbins, kbintype=kbintype)
xcpdps2_b_avg_kbin_models = []
excpdps2_b_avg_kbin_models = []
if model_hdf5files is not None:
for i in range(len(model_hdf5files)):
xcpdps2_a_avg_kbin_models += [BSP.incoherent_kbin_averaging(avg_xcpdps2_a_models[i], kbins=kbins, kbintype=kbintype)]
excpdps2_a_avg_kbin_models += [BSP.incoherent_kbin_averaging(avg_excpdps2_a_models[i], kbins=kbins, kbintype=kbintype)]
xcpdps2_b_avg_kbin_models += [BSP.incoherent_kbin_averaging(avg_xcpdps2_b_models[i], kbins=kbins, kbintype=kbintype)]
excpdps2_b_avg_kbin_models += [BSP.incoherent_kbin_averaging(avg_excpdps2_b_models[i], kbins=kbins, kbintype=kbintype)]
In [40]:
lstind = [0]
triadind = [0]
dayind = [0]
dayind_models = NP.zeros(len(model_labels), dtype=int).reshape(1,-1)
for stat in statistic:
for zind in spwind:
for lind in lstind:
for di,dind in enumerate(dayind):
for pstype in ['PS', 'Del2']:
for combi in range(len(diagoffsets_b)):
maxabsvals = []
minabsvals = []
maxvals = []
minvals = []
fig, axs = PLT.subplots(nrows=1, ncols=len(datapool), sharex=True, sharey=True, figsize=(4.0*len(datapool), 3.6))
if len(datapool) == 1:
axs = [axs]
for dpoolind,dpool in enumerate(datapool):
for trno,trind in enumerate(triadind):
# if model_hdf5files is not None:
# for mdlind, mdl in enumerate(model_labels):
# if dpool in xcpdps2_b_avg_kbin_models[mdlind][sampling]:
# if pstype == 'PS':
# psval = (2/3.0) * xcpdps2_b_avg_kbin_models[mdlind][sampling][dpool][stat][pstype][combi][zind,lind,dayind_models[di][mdlind],trind,:].to(pspec_unit).value
# else:
# psval = (2/3.0) * xcpdps2_b_avg_kbin_models[mdlind][sampling][dpool][stat][pstype][combi][zind,lind,dayind_models[di][mdlind],trind,:].to('K2').value
# kval = xcpdps2_b_avg_kbin_models[mdlind][sampling]['kbininfo'][dpool][stat][combi][zind,lind,dayind_models[di][mdlind],trind,:].to('Mpc-1').value
# maxabsvals += [NP.nanmin(NP.abs(psval.real))]
# minabsvals += [NP.nanmin(NP.abs(psval.real))]
# maxvals += [NP.nanmax(psval.real)]
# minvals += [NP.nanmin(psval.real)]
# axs[dpoolind].plot(kval, psval.real, ls='none', marker='.', ms=3, color=mdl_colrs[mdlind], label='{0}'.format(mdl))
if dpool in xcpdps2_b_avg_kbin[sampling]:
if pstype == 'PS':
psval = (2/3.0) * xcpdps2_b_avg_kbin[sampling][dpool][stat][pstype][combi][zind,lind,dind,trind,:].to(pspec_unit).value
psrms_ssdiff = (2/3.0) * NP.nanstd(excpdps2_b_avg_kbin[sampling]['errinfo'][stat][pstype][combi][zind,lind,:,trind,:], axis=0).to(pspec_unit).value
psrms_psdiff = (2/3.0) * (xcpdps2_a_avg_kbin[sampling][dpool][stat][pstype][combi][zind,lind,1,1,trind,:] - xcpdps2_a_avg_kbin[sampling][dpool][stat][pstype][combi][zind,lind,0,0,trind,:]).to(pspec_unit).value
else:
psval = (2/3.0) * xcpdps2_b_avg_kbin[sampling][dpool][stat][pstype][combi][zind,lind,dind,trind,:].to('mK2').value
psrms_ssdiff = (2/3.0) * NP.nanstd(excpdps2_b_avg_kbin[sampling]['errinfo'][stat][pstype][combi][zind,lind,:,trind,:], axis=0).to('mK2').value
psrms_psdiff = (2/3.0) * (xcpdps2_a_avg_kbin[sampling][dpool][stat][pstype][combi][zind,lind,1,1,trind,:] - xcpdps2_a_avg_kbin[sampling][dpool][stat][pstype][combi][zind,lind,0,0,trind,:]).to('mK2').value
if 2 in avg_incohax_b[combi]:
ind_dayax_in_incohax = avg_incohax_b[combi].index(2)
if 0 in diagoffsets_incohax_b[combi][ind_dayax_in_incohax]:
rms_inflation_factor = 2.0 * NP.sqrt(2.0)
else:
rms_inflation_factor = NP.sqrt(2.0)
else:
rms_inflation_factor = NP.sqrt(2.0)
psrms_psdiff = NP.abs(psrms_psdiff.real) / rms_inflation_factor
psrms_max = NP.amax(NP.vstack((psrms_ssdiff, psrms_psdiff)), axis=0)
kval = xcpdps2_b_avg_kbin[sampling]['kbininfo'][dpool][stat][combi][zind,lind,dind,trind,:].to('Mpc-1').value
maxabsvals += [NP.nanmax(NP.abs(psval.real + nsigma*psrms_max.real))]
minabsvals += [NP.nanmin(NP.abs(psval.real))]
maxvals += [NP.nanmax(psval.real + nsigma*psrms_max.real)]
minvals += [NP.nanmin(psval.real - nsigma*psrms_max.real)]
for errtype in ps_errtype:
if errtype.lower() == 'ssdiff':
axs[dpoolind].errorbar(kval, psval.real, yerr=nsigma*psrms_ssdiff, xerr=None, ecolor=errshade[errtype.lower()], ls='none', marker='.', ms=4, color='black')
elif errtype.lower() == 'psdiff':
axs[dpoolind].errorbar(kval, psval.real, yerr=nsigma*psrms_psdiff, xerr=None, ecolor=errshade[errtype.lower()], ls='none', marker='.', ms=4, color='black', label='FG+N')
# legend = axs[dpoolind].legend(loc='center', bbox_to_anchor=(0.5,0.3), shadow=False, fontsize=8)
if trno == 0:
# axs[dpoolind].text(0.05, 0.97, 'Real', transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
axs[dpoolind].text(0.95, 0.97, r'$z=$'+' {0:.1f}'.format(xcpdps2_b_avg_kbin['resampled']['z'][zind]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='right', va='top', color='black')
# axs[dpoolind].text(0.05, 0.92, r'$\Delta$'+'LST = {0:.1f} s'.format(lind*3.6e3*xcpdps2_a_avg_kbin['dlst'][0]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
# axs[dpoolind].text(0.05, 0.87, 'G{0[0]:0d}{0[1]:0d}'.format(dind), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
axs[dpoolind].axhline(y=0, xmin=0, xmax=1, ls='-', lw=1, color='black')
minvals = NP.asarray(minvals)
maxvals = NP.asarray(maxvals)
minabsvals = NP.asarray(minabsvals)
maxabsvals = NP.asarray(maxabsvals)
axs[dpoolind].set_xlim(0.99*NP.nanmin(xcpdps2_b_avg_kbin['resampled']['kbininfo']['kbin_edges'][zind].to('Mpc-1').value), 1.01*NP.nanmax(xcpdps2_b_avg_kbin['resampled']['kbininfo']['kbin_edges'][zind].to('Mpc-1').value))
if NP.min(minvals) < 0.0:
axs[dpoolind].set_ylim(1.5*NP.nanmin(minvals), 2*NP.nanmax(maxabsvals))
else:
axs[dpoolind].set_ylim(0.5*NP.nanmin(minvals), 2*NP.nanmax(maxabsvals))
axs[dpoolind].set_yscale('symlog', linthreshy=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))))
tickloc = PLTick.SymmetricalLogLocator(linthresh=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))), base=100.0)
axs[dpoolind].set_yscale('symlog', linthreshy=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))))
tickloc = PLTick.SymmetricalLogLocator(linthresh=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))), base=100.0)
axs[dpoolind].yaxis.set_major_locator(tickloc)
axs[dpoolind].grid(color='0.8', which='both', linestyle=':', lw=1)
fig.subplots_adjust(top=0.85)
fig.subplots_adjust(bottom=0.16)
fig.subplots_adjust(left=0.22)
fig.subplots_adjust(right=0.98)
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
# big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$\kappa_\parallel$'+' [pseudo '+r'$h$'+' Mpc'+r'$^{-1}$'+']', fontsize=12, weight='medium', labelpad=20)
if pstype == 'PS':
big_ax.set_ylabel(r'$\frac{1}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo mK$^2h^{-3}$ Mpc$^3$]', fontsize=12, weight='medium', labelpad=40)
else:
big_ax.set_ylabel(r'$\frac{1}{3}\, \Delta_\nabla^2(\kappa_\parallel)$ [pseudo mK$^2$]', fontsize=12, weight='medium', labelpad=40)
# big_axt = big_ax.twiny()
# big_axt.set_xticks([])
# big_axt.set_xlabel(r'$\tau$'+' ['+r'$\mu$'+'s]', fontsize=12, weight='medium', labelpad=20)
if pstype == 'PS':
PLT.savefig(figdir + '{0}_symlog_incoh_kbin_avg_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_a_avg_kbin[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg_kbin['dlst'][0], applyflags_str, combi), bbox_inches=0)
print(figdir + '{0}_symlog_incoh_kbin_avg_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_a_avg_kbin[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg_kbin['dlst'][0], applyflags_str, combi))
else:
PLT.savefig(figdir + '{0}_symlog_incoh_kbin_avg_real_cpDel2_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_a_avg_kbin[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg_kbin['dlst'][0], applyflags_str, combi), bbox_inches=0)
figdir + '{0}_symlog_incoh_kbin_avg_real_cpDel2_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_a_avg_kbin[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg_kbin['dlst'][0], applyflags_str, combi)
In [ ]:
if ('3' in plots) or ('3a' in plots) or ('3b' in plots) or ('3c' in plots):
sampling = plot_info['3']['sampling']
HI_PS_dir = plot_info['3']['21cm_PS_dir']
sim_rootdir = plot_info['3']['sim_rootdir']
visdirs = plot_info['3']['visdirs']
simvisdirs = [sim_rootdir+visdir for visdir in visdirs]
simlabels = plot_info['3']['simlabels']
visfile_prefix = plot_info['3']['visfile_prfx']
theory_HI_PS_files = glob.glob(HI_PS_dir+'ps_*')
z_theory_HI_PS_files = NP.asarray([fname.split('/')[-1].split('_')[3].split('z')[1] for fname in theory_HI_PS_files], dtype=NP.float)
h_Planck15 = DS.cosmoPlanck15.h
z_freq_window_centers = CNST.rest_freq_HI / freq_window_centers - 1
psfile_inds = [NP.argmin(NP.abs(z_theory_HI_PS_files - z_freq_window_center)) for z_freq_window_center in z_freq_window_centers]
In [ ]:
simvis_objs = [RI.InterferometerArray(None, None, None, init_file=simvisdir+visfile_prefix) for simvisdir in simvisdirs]
simDS_objs = [DS.DelaySpectrum(interferometer_array=simvis_obj) for simvis_obj in simvis_objs]
simDPS_objs = []
for simind,simlbl in enumerate(simlabels):
dspec = simDS_objs[simind].delay_transform(action='store')
subband_dspec = simDS_objs[simind].subband_delay_transform({key: freq_window_bw for key in ['cc', 'sim']}, freq_center={key: freq_window_centers for key in ['cc', 'sim']}, shape={key: freq_window_shape for key in ['cc', 'sim']}, fftpow={key: freq_window_fftpow for key in ['cc', 'sim']}, pad={key: pad for key in ['cc', 'sim']}, bpcorrect=False, action='return_resampled')
simDPS_objs = []
for simind,simlbl in enumerate(simlabels):
simDPS_objs += [DS.DelayPowerSpectrum(simDS_objs[simind])]
simDPS_objs[simind].compute_power_spectrum()
In [ ]:
In [28]:
testfile = datadir+'test_xcpdps.hdf5'
save_CPhase_cross_power_spectrum(xcpdps2, testfile)
In [29]:
testefile = datadir+'test_excpdps.hdf5'
save_CPhase_cross_power_spectrum(xcpdps2_errinfo, testefile)
In [30]:
xcpdps2_copy = read_CPhase_cross_power_spectrum(testfile)
In [36]:
xcpdps2_copy['resampled']['whole']['mean'].shape
Out[36]:
In [37]:
excpdps2_copy = read_CPhase_cross_power_spectrum(testefile)
In [58]:
if '2c' in plots:
avg_incohax = plot_info['2c']['incohax']
diagoffsets_incohax = plot_info['2c']['diagoffsets']
diagoffsets = []
for combi,incax_comb in enumerate(avg_incohax):
diagoffsets += [{}]
for incaxind,incax in enumerate(incax_comb):
diagoffsets[-1][incax] = NP.asarray(diagoffsets_incohax[combi][incaxind])
print(diagoffsets)
xcpdps2_avg, excpdps2_avg = BSP.incoherent_cross_power_spectrum_average([xcpdps2, xcpdps2_copy], excpdps=[xcpdps2_errinfo, excpdps2_copy], diagoffsets=diagoffsets)
In [59]:
print(len(xcpdps2_avg['resampled']['whole']['mean']))
print(xcpdps2_avg['resampled']['whole']['mean'][0].shape)
In [140]:
kbins = None
kbins = NP.geomspace(1e-3, 3, num=20, endpoint=True)
xcpdps2_avg_kbin = incoherent_kbin_averaging(xcpdps2_avg, kbins=kbins, kbintype='log')
In [100]:
print(xcpdps2_avg_kbin['resampled']['kbininfo']['counts'][0].shape)
print(xcpdps2_avg_kbin['resampled']['kbininfo']['kbin_edges'][0].shape)
print(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][0].shape)
print(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][0][0,0,0,1,0,:])
print(xcpdps2_avg_kbin['resampled']['z'])
print(xcpdps2_avg_kbin.keys())
print(xcpdps2_avg_kbin['resampled'].keys())
In [134]:
fig = PLT.figure(); ax = fig.add_subplot(111)
ax.plot(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][0][0,0,0,1,0,:], xcpdps2_avg_kbin['resampled']['whole']['mean']['PS'][0][0,0,0,1,0,:].to('K2 Mpc3').real / 3, ls='none', marker='.', ms=3)
ax.plot(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][1][0,0,0,1,0,:], xcpdps2_avg_kbin['resampled']['whole']['mean']['PS'][1][0,0,0,1,0,:].to('K2 Mpc3').real / 3, ls='none', marker='.', ms=3)
ax.plot(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][2][0,0,0,1,0,:], xcpdps2_avg_kbin['resampled']['whole']['mean']['PS'][2][0,0,0,1,0,:].to('K2 Mpc3').real / 3, ls='none', marker='.', ms=3)
ax.set_yscale('symlog', linthreshy=1e-1)
In [137]:
fig = PLT.figure(); ax = fig.add_subplot(111)
ax.plot(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][0][0,0,0,1,0,:], xcpdps2_avg_kbin['resampled']['whole']['mean']['Del2'][0][0,0,0,1,0,:].to('K2').real / 3, ls='none', marker='.', ms=3)
ax.plot(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][1][0,0,0,1,0,:], xcpdps2_avg_kbin['resampled']['whole']['mean']['Del2'][1][0,0,0,1,0,:].to('K2').real / 3, ls='none', marker='.', ms=3)
ax.plot(xcpdps2_avg_kbin['resampled']['kbininfo']['whole']['mean'][2][0,0,0,1,0,:], xcpdps2_avg_kbin['resampled']['whole']['mean']['Del2'][2][0,0,0,1,0,:].to('K2').real / 3, ls='none', marker='.', ms=3)
ax.set_yscale('symlog', linthreshy=1e-3)
# ax.set_xscale('log')
In [5]:
bl = NP.asarray(dspecinfo['bl'])
if bl.shape[0] != 3:
raise ValueError('Input bl must be made of three vectors forming the triad')
bltol = dspecinfo['bltol']
cohax = dspecinfo['cohax']
incohax = dspecinfo['incohax']
collapseax = dspecinfo['collapseax']
infile = infiles[0]
infile_no_ext = infile.split('.npz')[0]
In [6]:
visinfo = NMO.load_dict_from_hdf5(visfile)
blind, blrefind, dbl = LKP.find_1NN(visinfo['baseline']['blvect'], bl, distance_ULIM=bltol, remove_oob=True)
if blrefind.size != 3:
raise ValueError('Exactly three baselines were not found in the reference baselines')
vistriad = MA.array(visinfo['vis_real'][blrefind,:,:] + 1j * visinfo['vis_imag'][blrefind,:,:], mask=visinfo['mask'][blrefind,:,:])
In [7]:
print(vistriad.shape)
vistriad.shape = (nbl, nlst, nchan)
In [8]:
tmpnpzdata = NP.load(datadir+infile)
nchan = tmpnpzdata['flags'].shape[-1]
freqs = band_center + freq_resolution * (NP.arange(nchan) - int(0.5*nchan))
infile = 'EQ14.hdf5'
In [9]:
cpObj = BSP.ClosurePhase(datadir+infile, freqs, infmt='hdf5')
cpDSobj = BSP.ClosurePhaseDelaySpectrum(cpObj)
In [10]:
cpObj.cpinfo.keys()
Out[10]:
In [11]:
print(cpObj.cpinfo['raw'].keys())
print(cpObj.cpinfo['raw']['triads'].shape)
print(cpObj.cpinfo['raw']['triads'])
print(cpObj.cpinfo['raw']['cphase'].shape)
print(cpObj.cpinfo['raw']['days'].shape)
print(cpObj.cpinfo['raw']['lst-day'].shape)
print(cpObj.cpinfo['raw']['lst'])
print(cpObj.cpinfo['raw']['flags'].shape)
In [12]:
print(cpObj.cpinfo['processed'].keys())
In [13]:
print(cpObj.cpinfo['processed']['native'].keys())
print(cpObj.cpinfo['processed']['native']['cphase'].shape)
print(cpObj.cpinfo['processed']['native']['wts'].shape)
print(cpObj.cpinfo['processed']['native']['eicp'].shape)
In [14]:
print(cpObj.cpinfo['processed']['prelim'])
In [15]:
print(daybinsize)
print(lstbinsize)
In [16]:
cpObj.smooth_in_tbins(daybinsize=daybinsize, lstbinsize=lstbinsize)
In [56]:
print(cpObj.cpinfo['processed']['prelim'].keys())
print(cpObj.cpinfo['processed']['prelim']['cphase'].keys())
print(cpObj.cpinfo['processed']['prelim']['cphase']['mean'].shape)
print(cpObj.cpinfo['processed']['prelim']['eicp']['median'].shape)
print(cpObj.cpinfo['processed']['prelim']['wts'].shape)
print(cpObj.cpinfo['processed']['prelim']['daybins'])
print(cpObj.cpinfo['processed']['prelim']['lstbins'])
print(cpObj.cpinfo['processed']['prelim']['dlstbins'])
print(apply_flags)
In [36]:
cpObj.subtract(NP.zeros(1024))
visscaleinfo = {'vis': vistriad, 'lst': visinfo['header']['LST'],
'smoothinfo': {'op_type': 'interp1d', 'interp_kind': 'linear'}}
In [37]:
print(cpObj.cpinfo['processed'].keys())
In [38]:
cpds = cpDSobj.FT(freq_window_bw, freq_center=freq_window_centers, shape=freq_window_shape,
fftpow=freq_window_fftpow, pad=pad, datapool='prelim', visscaleinfo=visscaleinfo,
method='fft', resample=True, apply_flags=apply_flags)
In [39]:
plot_info = parms['plot']
plots = [key for key in plot_info if plot_info[key]['action']]
In [43]:
# triad = tuple(plot_info['1b']['triad'])
triad = (0, 1, 12)
triad_ind = triads.index(triad)
net_wts_raw = wts_raw[:,0,triad_ind,:][NP.newaxis,:,:] * freq_wts[:,NP.newaxis,:] # nspw x nlst x nchan
net_wts_proc = wts_proc[:,0,triad_ind,:][NP.newaxis,:,:] * freq_wts[:,NP.newaxis,:] # nspw x nlst x nchan
In [44]:
nrow = freq_wts.shape[0]
fig, axs = PLT.subplots(nrows=nrow, sharex=True, sharey=True)
if nrow == 1:
axs = [axs]
for axind in range(len(axs)):
wtsimg = axs[axind].imshow(net_wts_proc[axind,:,:], origin='lower', extent=[1e-6*freqs.min(), 1e-6*freqs.max(), tbins.min(), tbins.max()+NP.mean(dtbins)], norm=PLTC.LogNorm(vmin=1e-6, vmax=net_wts_proc.max()), interpolation='none', cmap='binary')
if axind == 0:
axs[axind].text(0.97, 0.97, '({0[0]:0d}, {0[1]:0d}, {0[2]:0d})'.format(triad), transform=axs[axind].transAxes, fontsize=12, weight='semibold', ha='right', va='top', color='red')
axs[axind].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
axs[axind].set_ylim(tbins.min(), tbins.max()+NP.mean(dtbins))
axs[axind].set_aspect('auto')
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(top=0.95)
fig.subplots_adjust(left=0.2)
fig.subplots_adjust(bottom=0.12)
fig.subplots_adjust(right=0.85)
cbax = fig.add_axes([0.86, 0.12, 0.02, 0.85])
cbar = fig.colorbar(wtsimg, cax=cbax, orientation='vertical')
cbax.yaxis.tick_right()
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
# big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium', labelpad=20)
big_ax.set_ylabel('LST [seconds]', fontsize=12, weight='medium', labelpad=35)
# PLT.savefig(figdir + '{0}_time_frequency_netwts_triad_{1[0]:0d}_{1[1]:0d}_{1[2]:0d}.png'.format(infile_no_ext, triad),
# bbox_inches=0)
# PLT.savefig(figdir + '{0}_time_frequency_netwts_triad_{1[0]:0d}_{1[1]:0d}_{1[2]:0d}.eps'.format(infile_no_ext, triad),
# bbox_inches=0)
Out[44]:
In [45]:
ncol = 5
nrow = min(6, int(NP.ceil(1.0*ntriads/ncol)))
npages = int(NP.ceil(1.0 * ntriads / (nrow*ncol)))
for pagei in range(npages):
if pagei > 0:
ntriads_remain = ntriads - pagei * nrow * ncol
nrow = min(6, int(NP.ceil(1.0*ntriads_remain/ncol)))
fig, axs = PLT.subplots(nrows=nrow, ncols=ncol, sharex=True, sharey=True, figsize=(8,6.4))
for i in range(nrow):
for j in range(ncol):
if i*ncol+j < ntriads:
axs[i,j].imshow(wts_raw[:,0,i*ncol+j,:], origin='lower',
extent=[1e-6*freqs.min(), 1e-6*freqs.max(), lst.min(), lst.max()+NP.mean(dlst)],
vmin=0, vmax=1, interpolation='none', cmap='gray')
axs[i,j].text(0.5, 0.97, '({0[0]:0d}, {0[1]:0d}, {0[2]:0d})'.format(triads[i*ncol+j]),
transform=axs[i,j].transAxes, fontsize=10, weight='medium', ha='center', va='top',
color='red')
else:
axs[i,j].axis('off')
axs[i,j].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
axs[i,j].set_ylim(lst.min(), lst.max()+NP.mean(dlst))
axs[i,j].set_aspect('auto')
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(top=0.95)
fig.subplots_adjust(left=0.1)
fig.subplots_adjust(bottom=0.15)
fig.subplots_adjust(right=0.98)
big_ax = fig.add_subplot(111)
# big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
big_ax.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium', labelpad=20)
big_ax.set_ylabel('LST [hour]', fontsize=12, weight='medium', labelpad=35)
In [52]:
autoinfo = {'axes': cohax} # Handles operations along coherent axes
xinfo = {'axes': incohax, 'avgcov': False, 'collapse_axes': collapseax} # Handles operations along incoherent axes
In [53]:
print(autoinfo)
print(xinfo)
print(selection)
In [55]:
xcpdps = cpDSobj.compute_power_spectrum(selection=selection, autoinfo=autoinfo, xinfo=xinfo)
In [ ]:
datapool = plot_info['2']['datapool']
statistic = plot_info['2']['statistic']
nsamples_incoh = xcpdps[datapool]['nsamples_incoh']
nsamples_coh = xcpdps[datapool]['nsamples_coh']
avgax = plot_info['2']['avgax']
if avgax is None:
avgax = [1,2,3]
lstind = plot_info['2']['lstind']
dayind = plot_info['2']['dayind']
triadind = plot_info['2']['triadind']
subselection = plot_info['2']['selection']
axkeys = ['spw', 'lst', 'day', 'triad', 'lags']
axnums = range(len(axkeys))
nominal_num_to_axes = {axnum: axk for axnum,axk in zip(axnums, axkeys)}
nominal_axes_to_num = {v: k for k, v in nominal_num_to_axes.items()}
actual_axes_to_num = {}
for axn,axkey in zip(axnums,axkeys):
if axn in xcpdps[datapool]['axesmap']:
actual_axes_to_num[axkey] = xcpdps[datapool]['axesmap'][axn][0]
else:
actual_axes_to_num[axkey] = axn
for key in xcpdps[datapool]['axesmap']:
if axn > key:
actual_axes_to_num[axkey] += len(xcpdps[datapool]['axesmap'][key]) - 1
actual_num_to_axes = {v: k for k, v in actual_axes_to_num.items()}
lstdiagonal = plot_info['2']['lstdiag']
triaddiagonal = plot_info['2']['triaddiag']
daydiagonal = plot_info['2']['daydiag']
if lstdiagonal is not None:
if not isinstance(lstdiagonal, int):
raise TypeError('Input LST diagonal must be an integer')
if nominal_axes_to_num['lst'] in xinfo['collapse_axes']:
lstind = xcpdps[datapool]['diagoffsets'][xinfo['collapse_axes'].tolist().index(nominal_axes_to_num['lst'])].tolist().index(lstdiagonal)
else:
lstind = lstdiagonal
lstind = NP.asarray(lstind).ravel()
else:
lstind = NP.arange(xcpdps['lst'].size)
if triaddiagonal is not None:
if not isinstance(triaddiagonal, int):
raise TypeError('Input triad diagonal must be an integer')
if nominal_axes_to_num['triad'] in xinfo['collapse_axes']:
triadind = xcpdps[datapool]['diagoffsets'][xinfo['collapse_axes'].tolist().index(nominal_axes_to_num['triad'])].tolist().index(triaddiagonal)
else:
triadind = triaddiagonal
triadind = NP.asarray(triadind).ravel()
else:
triadind = NP.arange(xcpdps['triads'].shape[0])
if daydiagonal is not None:
if not isinstance(daydiagonal, int):
raise TypeError('Input day diagonal must be an integer')
if nominal_axes_to_num['day'] in xinfo['collapse_axes']:
dayind = xcpdps[datapool]['diagoffsets'][xinfo['collapse_axes'].tolist().index(nominal_axes_to_num['day'])].tolist().index(daydiagonal)
else:
dayind = daydiagonal
dayind = NP.asarray(dayind).ravel()
else:
dayind = NP.arange(xcpdps['days'].size)
spw_ind = NP.arange(xcpdps[datapool]['z'].size)
lst_ind = NP.asarray(lstind)
day_ind = NP.asarray(dayind)
triad_ind = NP.asarray(triadind)
lags_ind = NP.arange(xcpdps[datapool]['lags'].size)
axes_to_avg = NP.arange(1,xcpdps[datapool]['mean'].ndim-1)
if statistic is None:
statistic = ['mean', 'median']
In [ ]:
colrs = ['red', 'green', 'blue']
for stat in statistic:
psval = NP.mean(xcpdps[datapool][stat], axis=tuple(axes_to_avg))
fig = PLT.figure(figsize=(3.5,3.5))
ax = fig.add_subplot(111)
for zind,z in enumerate(xcpdps[datapool]['z']):
negind = psval[zind,:] < 0.0
posind = NP.logical_not(negind)
ax.plot(xcpdps[datapool]['kprll'][zind,posind], psval[zind,posind], ls='none', marker='.', ms=1, color=colrs[zind], label=r'$z$={0:.1f}'.format(z))
ax.plot(xcpdps[datapool]['kprll'][zind,negind], NP.abs(psval[zind,negind]), ls='none', marker='|', ms=1, color=colrs[zind])
ax.set_yscale('log')
ax.set_xlim(0.99*xcpdps[datapool]['kprll'][zind,:].min(), 1.01*xcpdps[datapool]['kprll'][zind,:].max())
ax.set_ylim(1e-3, 1e8)
ax.set_xlabel(r'$k_\parallel$'+' ['+r'$h$'+' Mpc'+r'$^{-1}$'+']', fontsize=12, weight='medium')
ax.set_ylabel(r'$P_\nabla(k_\parallel)$ [Jy$^2h^{-1}$ Mpc]', fontsize=12, weight='medium', labelpad=0)
legend = ax.legend(loc='upper right', shadow=False, fontsize=10)
axt = ax.twiny()
axt.set_xlim(1e6*xcpdps[datapool]['lags'].min(), 1e6*xcpdps[datapool]['lags'].max())
axt.set_xlabel(r'$\tau$'+' ['+r'$\mu$'+'s]', fontsize=12, weight='medium')
fig.subplots_adjust(top=0.85)
fig.subplots_adjust(bottom=0.16)
fig.subplots_adjust(left=0.2)
fig.subplots_adjust(right=0.98)
# PLT.savefig(figdir + '{0}_closure_phase_delay_power_spectra_{1}_{2}_triads_{3}x{4:.1f}sx{5:.1f}d_{6}_statistic_nsamples_incoh_{7}_flags_{8}.png'.format(infile_no_ext, datapool, xcpdps['triads_ind'].size, xcpdps['lst'].size, 3.6e3*xcpdps['dlst'][0], xcpdps['dday'][0], stat, nsamples_incoh, applyflags_str), bbox_inches=0)
# PLT.savefig(figdir + '{0}_closure_phase_delay_power_spectra_{1}_{2}_triads_{3}x{4:.1f}sx{5:.1f}d_{6}_statistic_nsamples_incoh_{7}_flags_{8}.eps'.format(infile_no_ext, datapool, xcpdps['triads_ind'].size, xcpdps['lst'].size, 3.6e3*xcpdps['dlst'][0], xcpdps['dday'][0], stat, nsamples_incoh, applyflags_str), bbox_inches=0)
In [ ]:
print(xcpdps[datapool][stat].shape)
In [ ]: