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)