In [1]:
import pprint, copy
import numpy as NP
import numpy.ma as MA
import matplotlib.pyplot as PLT
import yaml, argparse, warnings
import astropy.units as U
import astropy.cosmology as cosmology
import astropy.constants as FCNST
import astropy.cosmology as cosmology
import scipy.signal as SPS
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]:
cosmoPlanck15 = cosmology.Planck15 # Planck 2015 cosmology
cosmo100 = cosmoPlanck15.clone(name='Modified Planck 2015 cosmology with h=1.0', H0=100.0) # Modified Planck 2015 cosmology with h=1.0, H= 100 km/s/Mpc
In [3]:
bspdir = '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/'
figdir = bspdir + 'figures/'
In [4]:
freq_window_centers = NP.asarray([150e6])
freq_window_bw = NP.asarray([42e6])
freq_window_shape = 'bhw'
freq_window_fftpow = 2
pad = 7.0
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 [5]:
case_with_FG_HI = {'1ps_FG_no_spindex_HI_colocated': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex_HI_cosine_colocated/',
'bspfile_sfx': '_HERA19_1ps_FG_no_spindex_HI_cosine_colocated_UniformBeam_noiseless.hdf5'},
'1ps_FG_spindex_HI_colocated': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex_HI_cosine_colocated/',
'bspfile_sfx': '_HERA19_1ps_FG_spindex_HI_cosine_colocated_UniformBeam_noiseless.hdf5'},
'1ps_FG_no_spindex_HI_displaced': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex_HI_cosine_displaced/',
'bspfile_sfx': '_HERA19_1ps_FG_no_spindex_HI_cosine_displaced_UniformBeam_noiseless.hdf5'},
'1ps_FG_spindex_HI_displaced': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex_HI_cosine_displaced/',
'bspfile_sfx': '_HERA19_1ps_FG_spindex_HI_cosine_displaced_UniformBeam_noiseless.hdf5'},
'GLEAM_1ps_HI': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/toy_GLEAM_1ps_HI_cosine/',
'bspfile_sfx': '_HERA19_GLEAM_FG_FoV-30_HI_cosine_on_zenith_UniformBeam_noiseless.hdf5'},
'GLEAM_HI_21cmfast': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/toy_GLEAM_HI_21cmfast_FaintGalaxies_fiducial_FoV-30/',
'bspfile_sfx': '_HERA19_GLEAM_FG_FoV-30_HI_21cmfast_FaintGalaxies_fiducial_UniformBeam_noiseless.hdf5'},
'1ps_FG_no_spindex_HI_21cmfast': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex_on_zenith_toy_HI_21cmfast_FaintGalaxies_fiducial/',
'bspfile_sfx': '_HERA19_1ps_FG_no_spindex_on_zenith_HI_21cmfast_FaintGalaxies_fiducial_UniformBeam_noiseless.hdf5'},
'1ps_FG_spindex_HI_21cmfast': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex_on_zenith_toy_HI_21cmfast_FaintGalaxies_fiducial/',
'bspfile_sfx': '_HERA19_1ps_FG_no_spindex_on_zenith_HI_21cmfast_FaintGalaxies_fiducial_UniformBeam_noiseless.hdf5'}
}
In [6]:
key_mapping_FG_HI = {'1ps_FG_no_spindex_HI_colocated': {'FG': '1ps_FG_no_spindex_on_zenith',
'HI': '1ps_HI_cosine'},
'1ps_FG_spindex_HI_colocated': {'FG': '1ps_FG_spindex_on_zenith',
'HI': '1ps_HI_cosine'},
'1ps_FG_no_spindex_HI_displaced': {'FG': '1ps_FG_no_spindex_off_zenith',
'HI': '1ps_HI_cosine'},
'1ps_FG_spindex_HI_displaced': {'FG': '1ps_FG_spindex_off_zenith',
'HI': '1ps_HI_cosine'},
'GLEAM_1ps_HI': {'FG': 'GLEAM',
'HI': '1ps_HI_cosine'},
'GLEAM_HI_21cmfast': {'FG': 'GLEAM',
'HI': 'HI_21cmfast_fiducial'},
'1ps_FG_no_spindex_HI_21cmfast': {'FG': '1ps_FG_no_spindex_on_zenith',
'HI': 'HI_21cmfast_fiducial'},
'1ps_FG_spindex_HI_21cmfast': {'FG': '1ps_FG_spindex_on_zenith',
'HI': 'HI_21cmfast_fiducial'}
}
In [7]:
case_with_FG = {'1ps_FG_no_spindex_on_zenith': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex/on_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/',
'bspfile_sfx': '_HERA19_1ps_FG_no_spindex_on_zenith_UniformBeam_noiseless.hdf5'},
'1ps_FG_no_spindex_off_zenith': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex/off_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/',
'bspfile_sfx': '_HERA19_1ps_FG_no_spindex_off_zenith_UniformBeam_noiseless.hdf5'},
'1ps_FG_spindex_on_zenith': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex/on_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/',
'bspfile_sfx': '_HERA19_1ps_FG_spindex_on_zenith_UniformBeam_noiseless.hdf5'},
'1ps_FG_spindex_off_zenith': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex/off_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/',
'bspfile_sfx': '_HERA19_1ps_FG_spindex_off_zenith_UniformBeam_noiseless.hdf5'},
'GLEAM': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/toy_GLEAM_FoV-30/on_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/',
'bspfile_sfx': '_HERA19_GLEAM_FG_FoV-30_UniformBeam_noiseless.hdf5'}
}
In [8]:
case_with_HI = {'1ps_HI_cosine': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_HI_cosine/on_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/'},
'HI_21cmfast_fiducial': {'visdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/toy_HI_21cmfast_FaintGalaxies_fiducial/on_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/'}
}
In [9]:
for case_id in case_with_FG:
print('\nLoading visibilities for case: {0}\n\t\t\t from: {1}'.format(case_id, case_with_FG[case_id]['visdir']))
case_with_FG[case_id]['vis'] = RI.InterferometerArray(None, None, None, init_file=case_with_FG[case_id]['visdir']+'simvis')
# Perform minor astrometric corrections when the point source foreground is supposed to be at zenith but is not
if 'on_zenith' in case_id:
print('\tPerforming minor astrometric corrections...')
case_with_FG[case_id]['vis'].skyvis_freq[:,:,0].real = NP.abs(case_with_FG[case_id]['vis'].skyvis_freq[:,:,0])
case_with_FG[case_id]['vis'].skyvis_freq[:,:,0].imag = NP.zeros_like(case_with_FG[case_id]['vis'].skyvis_freq[:,:,0].real)
case_with_FG[case_id]['vis'].vis_freq = case_with_FG[case_id]['vis'].skyvis_freq + case_with_FG[case_id]['vis'].vis_noise_freq
print('Loaded!\n----------')
In [10]:
for case_id in case_with_HI:
print('\nLoading visibilities for case: {0}\n\t\t\t from: {1}'.format(case_id, case_with_HI[case_id]['visdir']))
case_with_HI[case_id]['vis'] = RI.InterferometerArray(None, None, None, init_file=case_with_HI[case_id]['visdir']+'simvis')
# Perform minor astrometric corrections when the point source foreground is supposed to be at zenith but is not
if '1ps' in case_id:
print('\tPerforming minor astrometric corrections...')
case_with_HI[case_id]['vis'].skyvis_freq[:,:,0].real = NP.abs(case_with_HI[case_id]['vis'].skyvis_freq[:,:,0])
case_with_HI[case_id]['vis'].skyvis_freq[:,:,0].imag = NP.zeros_like(case_with_HI[case_id]['vis'].skyvis_freq[:,:,0].real)
case_with_HI[case_id]['vis'].vis_freq = case_with_HI[case_id]['vis'].skyvis_freq + case_with_HI[case_id]['vis'].vis_noise_freq
print('Loaded!\n----------')
In [11]:
for case_id in case_with_FG_HI:
print('\nLoading visibilities for case: {0}\n\t\t\t from: {1}'.format(case_id, case_with_FG_HI[case_id]['visdir']))
case_with_FG_HI[case_id]['vis'] = RI.InterferometerArray(None, None, None, init_file=case_with_FG_HI[case_id]['visdir']+'simvis')
if ('on_zenith' in key_mapping_FG_HI[case_id]['FG']) or ('1ps' in key_mapping_FG_HI[case_id]['HI']):
case_with_FG_HI[case_id]['vis'].skyvis_freq = case_with_FG[key_mapping_FG_HI[case_id]['FG']]['vis'].skyvis_freq + case_with_HI[key_mapping_FG_HI[case_id]['HI']]['vis'].skyvis_freq
case_with_FG_HI[case_id]['vis'].vis_freq = case_with_FG_HI[case_id]['vis'].skyvis_freq + case_with_FG_HI[case_id]['vis'].vis_noise_freq
print('Loaded!\n----------')
In [12]:
print(case_with_FG['1ps_FG_no_spindex_on_zenith']['vis'].skyvis_freq.shape)
print(case_with_FG['1ps_FG_no_spindex_on_zenith']['vis'].skyvis_freq[[0,1,2],:,0])
print(NP.abs(case_with_FG['1ps_FG_no_spindex_on_zenith']['vis'].skyvis_freq[[0,1,2],:,0].imag).max())
print(case_with_HI['1ps_HI_cosine']['vis'].skyvis_freq.shape)
print(case_with_HI['1ps_HI_cosine']['vis'].skyvis_freq[[0,1,2],:,0])
print(NP.abs(case_with_HI['1ps_HI_cosine']['vis'].skyvis_freq[[0,1,2],:,0].imag).max())
In [13]:
case_id = case_with_HI.keys()[0]
freqs = case_with_HI[case_id]['vis'].channels
antlocs = case_with_HI[case_id]['vis'].layout['positions']
antlabels = case_with_HI[case_id]['vis'].layout['labels']
baselines = case_with_HI[case_id]['vis'].baselines
baseline_labels = case_with_HI[case_id]['vis'].labels
print('Antenna Info:')
pprint.pprint(zip(antlabels, antlocs))
print('--------------')
print('Baseline info:')
pprint.pprint(zip(baseline_labels, baselines))
In [14]:
case_id = case_with_HI.keys()[0]
fig, ax = PLT.subplots(1, 1, figsize=(3.5,3.5))
ax.plot(antlocs[:,0], antlocs[:,1], ls='none', marker='o', mec='black', mfc='none', ms=16, color='black')
for antind,antlabel in enumerate(antlabels):
ax.annotate(antlabel, (antlocs[antind,0], antlocs[antind,1]), xycoords='data', ha='center', va='center', color='black')
ax.set_xlim(NP.min(antlocs[:,0])-5, NP.max(antlocs[:,0])+5)
ax.set_ylim(NP.min(antlocs[:,1])-5, NP.max(antlocs[:,1])+5)
ax.set_aspect('equal')
ax.set_xlabel('X [m]', fontsize=11, weight='medium')
ax.set_ylabel('Y [m]', fontsize=11, weight='medium')
fig.subplots_adjust(left=0.2, right=0.98, bottom=0.12, top=0.98)
PLT.savefig(figdir+'HERA-19_antenna_layout.pdf', bbox_inches=0)
In [15]:
cpinfo_case_with_FG_HI = {}
for case_id in case_with_FG_HI:
print('Computing closure phases on unique triads for {0}...'.format(case_id))
cpinfo_case_with_FG_HI[case_id] = case_with_FG_HI[case_id]['vis'].getClosurePhase(unique=True)
In [16]:
cpinfo_case_with_FG = {}
for case_id in case_with_FG:
print('Computing closure phases on unique triads for {0}...'.format(case_id))
cpinfo_case_with_FG[case_id] = case_with_FG[case_id]['vis'].getClosurePhase(unique=True)
In [17]:
cpinfo_case_with_HI = {}
for case_id in case_with_HI:
print('Computing closure phases on unique triads for {0}...'.format(case_id))
cpinfo_case_with_HI[case_id] = case_with_HI[case_id]['vis'].getClosurePhase(unique=True)
In [18]:
predicted_cpinfo_with_FG_HI = {}
for case_id in case_with_FG_HI:
imag_ratio_triad = (cpinfo_case_with_HI[key_mapping_FG_HI[case_id]['HI']]['skyvis'] / cpinfo_case_with_FG[key_mapping_FG_HI[case_id]['FG']]['skyvis']).imag
predicted_cpinfo_with_FG_HI[case_id] = {'imag_ratio_triad': imag_ratio_triad, 'sum_imag_ratio_triad': NP.sum(imag_ratio_triad, axis=1)}
print('Predicted closure phases on unique triads for {0}...'.format(case_id))
In [19]:
print(imag_ratio_triad.shape)
print(NP.sum(imag_ratio_triad, axis=1).shape)
pprint.pprint(cpinfo_case_with_HI[key_mapping_FG_HI[case_id]['HI']]['antenna_triplets'])
pprint.pprint(cpinfo_case_with_HI[key_mapping_FG_HI[case_id]['HI']]['baseline_triplets'])
In [20]:
print(cpinfo_case_with_HI.keys())
In [21]:
triad1_label = 'EQ14'
triad1 = ('0', '1', '8')
triad1_ind = cpinfo_case_with_HI[cpinfo_case_with_HI.keys()[0]]['antenna_triplets'].index(triad1)
triad1_baselines = cpinfo_case_with_HI[cpinfo_case_with_HI.keys()[0]]['baseline_triplets'][triad1_ind]
blind1 = []
for i in range(3):
blind = NP.where(NP.sum(NP.abs(baselines - triad1_baselines[[i]]), axis=1) <= 0.1)[0]
if blind.size == 1:
blind1 += [blind[0]]
if blind.size == 0:
blind = NP.where(NP.sum(NP.abs(baselines + triad1_baselines[[i]]), axis=1) <= 0.1)[0]
if blind.size == 1:
blind1 += [blind[0]]
else:
raise ValueError('Baseline [{0[0]:.1f}, {0[1]:.1f}, {0[2]:.1f}] nor its conjugate found'.format(triad1_baselines[i]))
bltriplet1_in_sim = NP.asarray([baselines[blind1[i]] for i in range(3)])
print(triad1_ind)
print(triad1_baselines)
print(blind1)
print(bltriplet1_in_sim)
In [22]:
triad2_label = 'EQ50'
triad2 = ('8', '11', '18')
triad2_ind = cpinfo_case_with_HI[cpinfo_case_with_HI.keys()[0]]['antenna_triplets'].index(triad2)
triad2_baselines = cpinfo_case_with_HI[cpinfo_case_with_HI.keys()[0]]['baseline_triplets'][triad2_ind]
blind2 = []
for i in range(3):
blind = NP.where(NP.sum(NP.abs(baselines - triad2_baselines[[i]]), axis=1) <= 0.1)[0]
if blind.size == 1:
blind2 += [blind[0]]
if blind.size == 0:
blind = NP.where(NP.sum(NP.abs(baselines + triad2_baselines[[i]]), axis=1) <= 0.1)[0]
if blind.size == 1:
blind2 += [blind[0]]
else:
raise ValueError('Baseline [{0[0]:.1f}, {0[1]:.1f}, {0[2]:.1f}] nor its conjugate found'.format(triad2_baselines[i]))
bltriplet2_in_sim = NP.asarray([baselines[blind2[i]] for i in range(3)])
print(triad2_ind)
print(triad2_baselines)
print(blind2)
print(bltriplet2_in_sim)
In [23]:
print(triad1_ind)
print(triad1)
print(freqs.shape)
print(predicted_cpinfo_with_FG_HI[case_id]['imag_ratio_triad'].shape)
print(cpinfo_case_with_FG_HI[case_id]['closure_phase_skyvis'].shape)
print(cpinfo_case_with_FG[key_mapping_FG_HI[case_id]['FG']]['closure_phase_skyvis'].shape)
print(predicted_cpinfo_with_FG_HI[case_id]['sum_imag_ratio_triad'].shape)
print((cpinfo_case_with_FG_HI[case_id]['closure_phase_skyvis'] - cpinfo_case_with_FG[key_mapping_FG_HI[case_id]['FG']]['closure_phase_skyvis'] - predicted_cpinfo_with_FG_HI[case_id]['sum_imag_ratio_triad']).shape)
In [24]:
print(figdir + 'closure_phase_comparison_actual_vs_prediction_HERA-19_triad_{0[0]}-{0[1]}-{0[2]}.pdf'.format(triad1))
In [27]:
colrs = ['red', 'blue', 'black']
for case_id in case_with_HI:
for blind, triad_ind, triad, triad_label in zip([blind1, blind2], [triad1_ind, triad2_ind], [triad1, triad2], [triad1_label, triad2_label]):
print('\tPlotting Triad:({0[0]},{0[1]},{0[2]})'.format(triad))
fig = PLT.figure(figsize=(3.5,2))
ax = fig.add_subplot(111)
for i in range(3):
ax.plot(freqs/1e6, NP.abs(case_with_HI[case_id]['vis'].skyvis_freq[blind[i],:,0]), color=colrs[i], ls='-', lw=2, alpha=0.5)
In [25]:
colrs = ['red', 'blue', 'black']
for case_id in case_with_FG_HI:
if (('GLEAM' in case_id) or ('21cmfast' in case_id)) and ('cosine' not in key_mapping_FG_HI[case_id]['HI']):
lwidth = 1
plot_envelope = False
else:
lwidth = 0.5
plot_envelope = True
print('Plotting {0}:'.format(case_id))
for blind, triad_ind, triad, triad_label in zip([blind1, blind2], [triad1_ind, triad2_ind], [triad1, triad2], [triad1_label, triad2_label]):
print('\tPlotting Triad:({0[0]},{0[1]},{0[2]})'.format(triad))
fig, axs = PLT.subplots(nrows=2, sharex=True, squeeze=True, figsize=(3.5,4))
for i in range(3):
axs[0].plot(freqs/1e6, NP.abs(case_with_FG[key_mapping_FG_HI[case_id]['FG']]['vis'].skyvis_freq[blind[i],:,0]), color=colrs[i], ls='-', lw=2, alpha=0.5)
axs[1].plot(freqs/1e6, 1e3 * (NP.abs(case_with_FG_HI[case_id]['vis'].skyvis_freq[blind[i],:,0]) - NP.abs(case_with_FG[key_mapping_FG_HI[case_id]['FG']]['vis'].skyvis_freq[blind[i],:,0])), color=colrs[i], ls='-', lw=lwidth, alpha=0.5)
if plot_envelope:
axs[1].plot(freqs/1e6, 1e3 * (NP.abs(SPS.hilbert(NP.abs(case_with_FG_HI[case_id]['vis'].skyvis_freq[blind[i],:,0]) - NP.abs(case_with_FG[key_mapping_FG_HI[case_id]['FG']]['vis'].skyvis_freq[blind[i],:,0])))), color=colrs[i], ls='-', lw=1)
axs[1].plot(freqs/1e6, -1e3 * (NP.abs(SPS.hilbert(NP.abs(case_with_FG_HI[case_id]['vis'].skyvis_freq[blind[i],:,0]) - NP.abs(case_with_FG[key_mapping_FG_HI[case_id]['FG']]['vis'].skyvis_freq[blind[i],:,0])))), color=colrs[i], ls='-', lw=1)
axs[0].text(0.98, 0.95, r'$\nabla$'+' ({0[0]},{0[1]},{0[2]})'.format(triad), transform=axs[0].transAxes, ha='right', va='center')
axs[0].set_ylabel(r'$|V_p^\mathrm{F}(f)|$ [Jy]', fontsize=11, weight='medium')
axs[1].set_ylabel(r'$\delta\,|V_p(f)|$ [mJy]', fontsize=11, weight='medium')
axs[0].set_xlim(105, 195)
axs[1].set_xlim(105, 195)
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(left=0.2, right=0.98, bottom=0.12, top=0.98)
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none')
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$f$ [MHz]', weight='medium', fontsize=11, labelpad=20)
PLT.savefig(figdir + '{0}_visibility_models_{1}_HERA-19.pdf'.format(triad_label, case_id), bbox_inches=0)
In [25]:
for case_id in case_with_FG_HI:
print('Plotting {0}'.format(case_id))
if ('GLEAM' in case_id) or ('21cmfast' in case_id):
lwidth = 1
plot_envelope = False
else:
lwidth = 0.5
plot_envelope = True
for triad_ind, triad, triad_label in zip([triad1_ind, triad2_ind], [triad1, triad2], [triad1_label, triad2_label]):
print('\tPlotting Triad:({0[0]},{0[1]},{0[2]})'.format(triad))
fig, axs = PLT.subplots(nrows=3, sharex=True, squeeze=True, figsize=(3.5,7.5))
colrs = ['red', 'blue', 'black']
for i in range(3):
axs[0].plot(freqs/1e6, 1e3 * predicted_cpinfo_with_FG_HI[case_id]['imag_ratio_triad'][triad_ind,i,:,0], ls='-', lw=lwidth, color=colrs[i], alpha=0.5)
if plot_envelope:
axs[0].plot(freqs/1e6, 1e3 * NP.abs(SPS.hilbert(predicted_cpinfo_with_FG_HI[case_id]['imag_ratio_triad'][triad_ind,i,:,0])), ls='-', lw=1, color=colrs[i])
axs[0].plot(freqs/1e6, -1e3 * NP.abs(SPS.hilbert(predicted_cpinfo_with_FG_HI[case_id]['imag_ratio_triad'][triad_ind,i,:,0])), ls='-', lw=1, color=colrs[i])
axs[0].text(0.98, 0.95, r'$\nabla$'+' ({0[0]},{0[1]},{0[2]})'.format(triad), transform=axs[0].transAxes, ha='right', va='center')
axs[0].set_ylabel(r'$\delta\phi_p^\mathrm{L}(f)$ [milli-radian]', fontsize=11, weight='medium')
axs[0].set_xlim(105, 195)
# axs[0].set_ylim(-0.15, 0.15)
axs[1].plot(freqs/1e6, 1e3 * (cpinfo_case_with_FG_HI[case_id]['closure_phase_skyvis'][triad_ind,:,0] - cpinfo_case_with_FG[key_mapping_FG_HI[case_id]['FG']]['closure_phase_skyvis'][triad_ind,:,0]), ls='-', lw=lwidth, color='black', alpha=0.5)
axs[1].plot(freqs/1e6, 1e3 * predicted_cpinfo_with_FG_HI[case_id]['sum_imag_ratio_triad'][triad_ind,:,0], ls='-', lw=lwidth, color='gray', alpha=0.5)
if plot_envelope:
axs[1].plot(freqs/1e6, 1e3 * NP.abs(SPS.hilbert(predicted_cpinfo_with_FG_HI[case_id]['sum_imag_ratio_triad'][triad_ind,:,0])), ls='-', lw=1, color='gray')
axs[1].plot(freqs/1e6, -1e3 * NP.abs(SPS.hilbert(predicted_cpinfo_with_FG_HI[case_id]['sum_imag_ratio_triad'][triad_ind,:,0])), ls='-', lw=1, color='gray')
axs[1].set_ylabel(r'$\delta\phi_\nabla^\mathrm{L}(f)$ [milli-radian]', fontsize=11, weight='medium')
axs[1].set_xlim(105, 195)
# axs[1].set_ylim(-0.25, 0.25)
axs[2].plot(freqs/1e6, 1e9 * (cpinfo_case_with_FG_HI[case_id]['closure_phase_skyvis'][triad_ind,:,0] - cpinfo_case_with_FG[key_mapping_FG_HI[case_id]['FG']]['closure_phase_skyvis'][triad_ind,:,0] - predicted_cpinfo_with_FG_HI[case_id]['sum_imag_ratio_triad'][triad_ind,:,0]), ls='-', lw=lwidth, color='gray')
axs[2].set_ylabel(r'err$\left(\delta\phi_\nabla^\mathrm{L}(f)\right)$ [nano-radian]', fontsize=11, weight='medium')
axs[2].set_xlim(105, 195)
# axs[2].set_ylim(-10, 10)
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(left=0.23, right=0.98, bottom=0.07, top=0.98)
big_ax = fig.add_subplot(111)
big_ax.set_facecolor('none')
big_ax.set_xticks([])
big_ax.set_yticks([])
big_ax.set_xlabel(r'$f$ [MHz]', weight='medium', labelpad=20)
PLT.savefig(figdir + '{0}_closure_phase_comparison_actual_vs_prediction_{1}_HERA-19.pdf'.format(triad_label, case_id), bbox_inches=0)
In [26]:
simDSobjs_FG = {}
simDSobjs_HI = {}
simDSobjs_FG_HI = {}
for case_id in case_with_FG:
simDSobjs_FG[case_id] = DS.DelaySpectrum(interferometer_array=case_with_FG[case_id]['vis'])
for case_id in case_with_HI:
simDSobjs_HI[case_id] = DS.DelaySpectrum(interferometer_array=case_with_HI[case_id]['vis'])
for case_id in case_with_FG_HI:
simDSobjs_FG_HI[case_id] = DS.DelaySpectrum(interferometer_array=case_with_FG_HI[case_id]['vis'])
In [27]:
print('Foregrounds only\n================')
for case_id in case_with_FG:
print('Computing standard subband delay transform for case "{0}"-----------------'.format(case_id))
dspec = simDSobjs_FG[case_id].delay_transform(action='store')
simDSobjs_FG[case_id].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=None)
print('-----------------')
In [28]:
print('HI only\n=======')
for case_id in case_with_HI:
print('Computing standard subband delay transform for case "{0}"-----------------'.format(case_id))
dspec = simDSobjs_HI[case_id].delay_transform(action='store')
simDSobjs_HI[case_id].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=None)
print('-----------------')
In [29]:
print('Foregrounds + HI\n================')
for case_id in case_with_FG_HI:
print('Computing standard subband delay transform for case "{0}"-----------------'.format(case_id))
dspec = simDSobjs_FG_HI[case_id].delay_transform(action='store')
simDSobjs_FG_HI[case_id].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=None)
print('-----------------')
In [30]:
simDPSobjs_FG = {}
simDPSobjs_HI = {}
simDPSobjs_FG_HI = {}
print('Foregrounds only\n================')
for case_id in case_with_FG:
print('Initializing and computing standard subband delay power spectra for case "{0}"'.format(case_id))
simDPSobjs_FG[case_id] = DS.DelayPowerSpectrum(simDSobjs_FG[case_id])
simDPSobjs_FG[case_id].compute_power_spectrum()
print('-----------------')
print('\nHI only\n=======')
for case_id in case_with_HI:
print('Initializing and computing standard subband delay power spectra for case "{0}"'.format(case_id))
simDPSobjs_HI[case_id] = DS.DelayPowerSpectrum(simDSobjs_HI[case_id])
simDPSobjs_HI[case_id].compute_power_spectrum()
print('-----------------')
print('\nForegrounds + HI\n================')
for case_id in case_with_FG_HI:
print('Initializing and computing standard subband delay power spectra for case "{0}"'.format(case_id))
simDPSobjs_FG_HI[case_id] = DS.DelayPowerSpectrum(simDSobjs_FG_HI[case_id])
simDPSobjs_FG_HI[case_id].compute_power_spectrum()
print('-----------------')
In [31]:
print(simDPSobjs_FG_HI[simDPSobjs_FG_HI.keys()[0]].subband_delay_power_spectra['sim'].keys())
print(simDPSobjs_FG_HI[simDPSobjs_FG_HI.keys()[0]].subband_delay_power_spectra['sim']['horizon_kprll_limits'].shape)
In [32]:
cpObjs_FG_HI = {}
cpObjs_FG = {}
for case_id in case_with_FG_HI:
cpObjs_FG_HI[case_id] = {}
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG_HI[case_id][triad_label] = BSP.ClosurePhase(bspdir+triad_label+case_with_FG_HI[case_id]['bspfile_sfx'],
freqs, infmt='hdf5')
for case_id in case_with_FG:
cpObjs_FG[case_id] = {}
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG[case_id][triad_label] = BSP.ClosurePhase(bspdir+triad_label+case_with_FG[case_id]['bspfile_sfx'],
freqs, infmt='hdf5')
In [33]:
pprint.pprint(cpObjs_FG_HI.keys())
print('-------------')
pprint.pprint(cpObjs_FG.keys())
print('-------------')
pprint.pprint(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo.keys())
print('-------------')
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['raw'].keys())
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['raw']['triads'].shape)
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['raw']['triads'])
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['raw']['days'].shape)
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['raw']['lst-day'].shape)
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['raw']['lst'].shape)
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['raw']['cphase'].shape)
In [34]:
cpObjs_FG_HI_predicted = copy.deepcopy(cpObjs_FG_HI) # Copy the object for convenience
for case_id in case_with_FG_HI:
print('Working on creating closure phase objects with first-order predictions for case: {0}'.format(case_id))
for triad_ind, triad_label in zip([triad1_ind, triad2_ind], [triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
# processed_dict = cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo.pop('processed', None) # Remove processed portions for recomputing
cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['errinfo'] = {}
# raw_cphase = cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['raw'].pop('cphase', None)
# print(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo.keys())
# print(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['raw'].keys())
cp_fluctuations = predicted_cpinfo_with_FG_HI[case_id]['sum_imag_ratio_triad'][triad_ind,:,:].T[:,NP.newaxis,NP.newaxis,:] # (nchan, nlst) --> (nlst,ndays=1,ntriads=1,nchan)
eicp_multiplicative_perturbations = (1 + 1j*cp_fluctuations) / NP.abs(1 + 1j*cp_fluctuations)
# print(eicp_multiplicative_perturbations.shape)
cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native']['eicp'] = cpObjs_FG[key_mapping_FG_HI[case_id]['FG']][triad_label].cpinfo['processed']['native']['eicp'] * eicp_multiplicative_perturbations
cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native']['cphase'] = MA.angle(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native']['eicp'])
cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['raw']['cphase'] = cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native']['cphase'].data
# cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['raw']['cphase'] = cpObjs_FG[key_mapping_FG_HI[case_id]['FG']][triad_label].cpinfo['raw']['cphase'] + cp_fluctuations # Add predicted fluctuations
# cpObjs_FG_HI_predicted[case_id][triad_label].expicp(force_action=True) # Force the initial computations and population of the dictionary attributes
# print(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo.keys())
# print(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['raw'].keys())
# print(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed'].keys())
# print(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native'].keys())
print(cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native']['cphase'].shape)
print('-------------')
In [35]:
daybinsize = None # days
ndaybins = 2
lstbinsize = 59.9 # seconds
print('day bin size = '+str(daybinsize), 'ndaybins={0:0d}'.format(ndaybins), 'LST bin size = {0:.1f}s'.format(lstbinsize))
for case_id in case_with_FG_HI:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG_HI[case_id][triad_label].smooth_in_tbins(daybinsize=daybinsize, ndaybins=ndaybins,
lstbinsize=lstbinsize)
for case_id in case_with_FG_HI:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG_HI_predicted[case_id][triad_label].smooth_in_tbins(daybinsize=daybinsize, ndaybins=ndaybins,
lstbinsize=lstbinsize)
# Force the closure phases and exp(i*cphase) to conform with predictions
for stat in ['mean', 'median']:
cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['prelim']['cphase'][stat] = cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native']['cphase'][:,[0],:,:] * NP.ones(ndaybins).reshape(1,-1,1,1)
cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['prelim']['eicp'][stat] = cpObjs_FG_HI_predicted[case_id][triad_label].cpinfo['processed']['native']['eicp'][:,[0],:,:] * NP.ones(ndaybins).reshape(1,-1,1,1)
for case_id in case_with_FG:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG[case_id][triad_label].smooth_in_tbins(daybinsize=daybinsize, ndaybins=ndaybins,
lstbinsize=lstbinsize)
In [36]:
pprint.pprint(cpObjs_FG[case_id][triad_label].cpinfo['processed'].keys())
print('--------------')
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim'].keys())
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim']['cphase'].keys())
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim']['lstbins'].shape)
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim']['lstbins'])
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim']['daybins'].shape)
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim']['daybins'])
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim']['cphase']['median'].shape)
print(cpObjs_FG[case_id][triad_label].cpinfo['processed']['prelim']['eicp']['median'].shape)
In [37]:
for case_id in case_with_FG_HI:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG_HI[case_id][triad_label].subtract(NP.zeros(1024))
for case_id in case_with_FG_HI:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG_HI_predicted[case_id][triad_label].subtract(NP.zeros(1024))
for case_id in case_with_FG:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG[case_id][triad_label].subtract(NP.zeros(1024))
print(cpObjs_FG[case_id][triad_label].cpinfo['processed'].keys())
In [38]:
# ndaybins=4
print('ndaybins={0:0d}'.format(4), 'LST bin size = {0:.1f}s'.format(lstbinsize))
for case_id in case_with_FG_HI:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG_HI[case_id][triad_label].subsample_differencing(daybinsize=None, ndaybins=4,
lstbinsize=lstbinsize)
for case_id in case_with_FG_HI:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG_HI_predicted[case_id][triad_label].subsample_differencing(daybinsize=None, ndaybins=4,
lstbinsize=lstbinsize)
for case_id in case_with_FG:
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpObjs_FG[case_id][triad_label].subsample_differencing(daybinsize=None, ndaybins=4,
lstbinsize=lstbinsize)
In [39]:
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo'].keys())
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['daybins'].shape)
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['lstbins'].shape)
print(len(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['list_of_pair_of_pairs']))
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['list_of_pair_of_pairs'])
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['eicp_diff'].keys())
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['eicp_diff']['0'].keys())
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['eicp_diff']['0']['median'].shape)
print(cpObjs_FG[cpObjs_FG.keys()[0]][triad_label].cpinfo['errinfo']['eicp_diff']['1']['median'].shape)
In [40]:
cpDSobjs_FG = {}
cpDSobjs_FG_HI = {}
cpDSobjs_FG_HI_predicted = {}
for case_id in cpObjs_FG:
cpDSobjs_FG[case_id] = {}
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpDSobjs_FG[case_id][triad_label] = BSP.ClosurePhaseDelaySpectrum(cpObjs_FG[case_id][triad_label])
for case_id in cpObjs_FG_HI:
cpDSobjs_FG_HI[case_id] = {}
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpDSobjs_FG_HI[case_id][triad_label] = BSP.ClosurePhaseDelaySpectrum(cpObjs_FG_HI[case_id][triad_label])
for case_id in cpObjs_FG_HI_predicted:
cpDSobjs_FG_HI_predicted[case_id] = {}
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
cpDSobjs_FG_HI_predicted[case_id][triad_label] = BSP.ClosurePhaseDelaySpectrum(cpObjs_FG_HI_predicted[case_id][triad_label])
In [41]:
visscaleinfo_FG = {}
for case_id in cpObjs_FG:
visscaleinfo_FG[case_id] = {}
for triad_baselines, triad_label in zip([triad1_baselines, triad2_baselines], [triad1_label, triad2_label]):
visscaleinfo_FG[case_id][triad_label] = {'vis': case_with_FG[case_id]['vis'],
'bltriplet': triad_baselines,
'smoothinfo': {'op_type': 'interp1d', 'interp_kind': 'linear'}
}
In [42]:
pprint.pprint(visscaleinfo_FG[cpObjs_FG.keys()[0]])
In [43]:
freq_window_centers = NP.asarray([150e6])
freq_window_bw = NP.asarray([42e6])
freq_window_shape = 'bhw'
freq_window_fftpow = 2
pad = 7.0
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 [44]:
cpds_FG = {}
print('Foregrounds only\n================')
for case_id in cpDSobjs_FG:
cpds_FG[case_id] = {}
print('Computing delay spectrum for case: {0}'.format(case_id))
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
cpds_FG[case_id][triad_label] = cpDSobjs_FG[case_id][triad_label].FT(freq_window_bw,
freq_center=freq_window_centers,
shape=freq_window_shape, fftpow=freq_window_fftpow,
pad=pad, datapool='prelim',
visscaleinfo=visscaleinfo_FG[case_id][triad_label],
method='fft', resample=True, apply_flags=False)
print('-------------')
In [45]:
cpds_FG_HI = {}
print('Foregrounds + HI (No approximation)\n===================================')
for case_id in cpDSobjs_FG_HI:
cpds_FG_HI[case_id] = {}
print('Computing delay spectrum for case: {0}'.format(case_id))
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
cpds_FG_HI[case_id][triad_label] = cpDSobjs_FG_HI[case_id][triad_label].FT(freq_window_bw,
freq_center=freq_window_centers,
shape=freq_window_shape,
fftpow=freq_window_fftpow,
pad=pad, datapool='prelim',
visscaleinfo=visscaleinfo_FG[key_mapping_FG_HI[case_id]['FG']][triad_label],
method='fft', resample=True, apply_flags=False)
print('-------------')
In [46]:
cpds_FG_HI_predicted = {}
print('Foregrounds + HI (linear-order prediction)\n==========================================')
for case_id in cpDSobjs_FG_HI_predicted:
cpds_FG_HI_predicted[case_id] = {}
print('Computing delay spectrum for case: {0}'.format(case_id))
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
cpds_FG_HI_predicted[case_id][triad_label] = cpDSobjs_FG_HI_predicted[case_id][triad_label].FT(freq_window_bw,
freq_center=freq_window_centers,
shape=freq_window_shape,
fftpow=freq_window_fftpow,
pad=pad, datapool='prelim',
visscaleinfo=visscaleinfo_FG[key_mapping_FG_HI[case_id]['FG']][triad_label],
method='fft', resample=True, apply_flags=False)
print('-------------')
In [47]:
print(cpDSobjs_FG_HI[cpDSobjs_FG_HI.keys()[0]][triad1_label].cPhaseDS['whole']['dspec']['mean'].shape)
print(NP.abs(cpDSobjs_FG_HI[cpDSobjs_FG_HI.keys()[0]][triad1_label].cPhaseDS['whole']['dspec']['mean']).max())
In [48]:
cohax = None
incohax = [1, 2, 3] # LST, days, triads
collapseax = [1, 2, 3] # Collapse along diagonals of LST, days, triads
# lst_in_sim = case_with_FG[case_with_FG.keys()[0]]['vis'].lst
# lstinds = NP.arange(lst_in_sim.size)
# dayind = None
selection = {}
for triad,triad_label in zip([triad1, triad2], [triad1_label, triad2_label]):
selection[triad_label] = None
# selection[triad_label] = {'triads': triad, 'lst': lstinds, 'days': dayind}
In [49]:
print(selection)
In [50]:
autoinfo = {'axes': cohax}
xinfo = {'axes': incohax, 'avgcov': False, 'collapse_axes': collapseax, 'dlst_range': [1.1]}
In [51]:
beaminfo = {'beamfile': None,
'filepathtype': 'default',
'filefmt': 'UVBeam',
'nside': 128,
'pol': 'X',
'chromatic': True,
'select_freq': 150e6,
'spec_interp': 'cubic',
'telescope': {'shape': 'delta',
'size': 14.0,
'orientation': [90.0, 270.0],
'ocoords': 'altaz',
'phased_array': False,
'ground_plane': None
}
}
In [52]:
xcpdps_FG = {}
print('Foregrounds only\n================')
for case_id in case_with_FG:
xcpdps_FG[case_id] = {}
print('Computing closure phase delay power spectrum for case: {0}'.format(case_id))
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
xcpdps_FG[case_id][triad_label] = cpDSobjs_FG[case_id][triad_label].compute_power_spectrum(selection=selection[triad_label],
autoinfo=autoinfo,
xinfo=xinfo, units='K',
beamparms=beaminfo)
print('-------------')
In [53]:
xcpdps_FG_HI = {}
print('Foregrounds + HI (No approximation)\n===================================')
for case_id in case_with_FG_HI:
xcpdps_FG_HI[case_id] = {}
print('Computing closure phase delay power spectrum for case: {0}'.format(case_id))
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
xcpdps_FG_HI[case_id][triad_label] = cpDSobjs_FG_HI[case_id][triad_label].compute_power_spectrum(selection=selection[triad_label],
autoinfo=autoinfo, xinfo=xinfo,
units='K', beamparms=beaminfo)
print('-------------')
In [54]:
xcpdps_FG_HI_predicted = {}
print('Foregrounds + HI (linear-order approximation)\n==============================================')
for case_id in case_with_FG_HI:
xcpdps_FG_HI_predicted[case_id] = {}
print('Computing closure phase delay power spectrum for case: {0}'.format(case_id))
for triadi, triad_label in enumerate([triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
xcpdps_FG_HI_predicted[case_id][triad_label] = cpDSobjs_FG_HI_predicted[case_id][triad_label].compute_power_spectrum(selection=selection[triad_label],
autoinfo=autoinfo, xinfo=xinfo,
units='K', beamparms=beaminfo)
print('-------------')
In [55]:
print(xcpdps_FG_HI[xcpdps_FG_HI.keys()[0]][triad1_label].keys())
print(xcpdps_FG_HI[xcpdps_FG_HI.keys()[0]][triad1_label]['triads'])
print(xcpdps_FG_HI[xcpdps_FG_HI.keys()[0]][triad1_label]['oversampled'].keys())
print(xcpdps_FG_HI[xcpdps_FG_HI.keys()[0]][triad1_label]['oversampled']['whole'].keys())
print(xcpdps_FG_HI[xcpdps_FG_HI.keys()[0]][triad1_label]['oversampled']['whole']['mean'].shape)
print(xcpdps_FG_HI[xcpdps_FG_HI.keys()[0]][triad1_label]['oversampled']['whole']['diagoffsets'])
In [56]:
print('Plotting deviations between actual and predicted closure phase power spectra')
for case_id in case_with_FG_HI:
print('{0}:'.format(case_id))
for blind_triplet, triad_label in zip([blind1, blind2], [triad1_label, triad2_label]):
print('\t{0}'.format(triad_label))
prediction_deviation = NP.abs(xcpdps_FG_HI_predicted[case_id][triad_label]['oversampled']['whole']['mean'].real.to('mK2 Mpc3').value - xcpdps_FG_HI[case_id][triad_label]['oversampled']['whole']['mean'].real.to('mK2 Mpc3').value)
# print(prediction_deviation.shape)
fig = PLT.figure()
ax = fig.add_subplot(111)
ax.plot(prediction_deviation[0,0,1,0,::8])
ax.set_yscale('log')
print('-----------------')
In [58]:
# print(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind1,0])
# print(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind1,1])
In [67]:
spw = [0]
lstind = 0
trind = 0
dayind = 0
delay_HI = 1e-6 # in seconds
colrs = ['red', 'blue', 'forestgreen']
harmonic_colrs = ['cyan', 'paleturquoise', 'lightcyan']
for spwind in spw:
for blind_triplet, triad, triad_label in zip([blind1, blind2], [triad1, triad2], [triad1_label, triad2_label]):
for case_id in case_with_FG_HI:
print('Plotting standard delay spectrum for case "{0}" on triad "{1}"'.format(case_id, triad_label))
fig, axs = PLT.subplots(nrows=1, ncols=2, squeeze=True, sharex=True, figsize=(6,4))
axs[0].axvspan(1e6*NP.min(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,0]), 1e6*NP.max(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,1]), facecolor='gold')
for blno, blind in enumerate(blind_triplet):
axs[0].plot(1e6*simDSobjs_FG[key_mapping_FG_HI[case_id]['FG']].subband_delay_spectra['sim']['lags'][::16], NP.abs(simDSobjs_FG[key_mapping_FG_HI[case_id]['FG']].subband_delay_spectra['sim']['skyvis_lag'][blind,spwind,::16,lstind]), ls='-', color=colrs[blno], alpha=0.5)
axs[0].plot(1e6*simDSobjs_HI[key_mapping_FG_HI[case_id]['HI']].subband_delay_spectra['sim']['lags'][::16], NP.abs(simDSobjs_HI[key_mapping_FG_HI[case_id]['HI']].subband_delay_spectra['sim']['skyvis_lag'][blind,spwind,::16,lstind]), ls=':', color=colrs[blno], alpha=0.5)
if 'cosine' in key_mapping_FG_HI[case_id]['HI']:
for i in range(1):
axs[0].axvspan(1e6*((i+1)*delay_HI + NP.min(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,0])), 1e6*((i+1)*delay_HI + NP.max(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,1])), facecolor=harmonic_colrs[i])
axs[0].axvspan(1e6*(-(i+1)*delay_HI + NP.min(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,0])), 1e6*(-(i+1)*delay_HI + NP.max(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,1])), facecolor=harmonic_colrs[i])
axs[0].annotate('', xy=(-(i+1), 1e6/(1e3)**i), xycoords='data', xytext=(-(i+1), 1e8/(1e3)**i), textcoords='data', arrowprops=dict(arrowstyle="->", connectionstyle="arc3"),)
axs[0].annotate('', xy=(i+1, 1e6/(1e3)**i), xycoords='data', xytext=(i+1, 1e8/(1e3)**i), textcoords='data', arrowprops=dict(arrowstyle="->", connectionstyle="arc3"),)
axs[0].text(0.98, 0.95, r'$\nabla$'+' ({0[0]},{0[1]},{0[2]})'.format(triad), transform=axs[0].transAxes, ha='right', va='center')
axs[0].set_yscale('log')
axs[0].set_xticks(NP.arange(-5,6))
# axs[0].tick_params(axis='x', rotation=45)
axs[0].set_xlim(-3.25, 3.25)
axs[0].set_ylim(1e-5, 1e10)
axs[0].set_xlabel(r'$\tau$ [$\mu$s]', fontsize=11, weight='medium')
axs[0].set_ylabel(r'$\mathcal{V}_p(\tau)$ [Jy Hz]', fontsize=11, weight='medium')
print('Plotting closure phase delay spectrum for case "{0}" on triad "{1}"'.format(case_id, triad_label))
cpds_lags = cpDSobjs_FG[key_mapping_FG_HI[case_id]['FG']][triad_label].cPhaseDS['lags']
cpdsval_FG = cpDSobjs_FG[key_mapping_FG_HI[case_id]['FG']][triad_label].cPhaseDS['whole']['dspec']['mean'][spwind,lstind,dayind,trind,:]
cpdsval_FG_HI = cpDSobjs_FG_HI[case_id][triad_label].cPhaseDS['whole']['dspec']['mean'][spwind,lstind,dayind,trind,:]
cpdsval_FG_HI_predicted = cpDSobjs_FG_HI_predicted[case_id][triad_label].cPhaseDS['whole']['dspec']['mean'][spwind,lstind,dayind,trind,:]
axs[1].axvspan(1e6*NP.sum(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,0]), 1e6*NP.sum(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,1]), facecolor='gold')
axs[1].plot(1e6*cpDSobjs_FG[key_mapping_FG_HI[case_id]['FG']][triad_label].cPhaseDS['lags'][::16], NP.abs(cpdsval_FG[::16]), ls='-', color='black')
axs[1].plot(1e6*cpDSobjs_FG_HI[case_id][triad_label].cPhaseDS['lags'][::16], NP.abs(cpdsval_FG_HI[::16]), ls=':', color='gray')
if 'cosine' in key_mapping_FG_HI[case_id]['HI']:
for i in range(3):
axs[1].axvspan(1e6*((i+1)*delay_HI + NP.sum(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,0])), 1e6*((i+1)*delay_HI + NP.sum(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,1])), facecolor=harmonic_colrs[i])
axs[1].axvspan(1e6*(-(i+1)*delay_HI + NP.sum(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,0])), 1e6*(-(i+1)*delay_HI + NP.sum(simDSobjs_FG_HI[case_id].horizon_delay_limits[lstind,blind_triplet,1])), facecolor=harmonic_colrs[i])
axs[1].annotate('', xy=(-(i+1)*delay_HI*1e6, 1e6/(1e3)**i), xycoords='data', xytext=(-(i+1)*delay_HI*1e6, 1e8/(1e3)**i), textcoords='data', arrowprops=dict(arrowstyle="->", connectionstyle="arc3"),)
axs[1].annotate('', xy=((i+1)*delay_HI*1e6, 1e6/(1e3)**i), xycoords='data', xytext=((i+1)*delay_HI*1e6, 1e8/(1e3)**i), textcoords='data', arrowprops=dict(arrowstyle="->", connectionstyle="arc3"),)
axs[1].plot(1e6*cpDSobjs_FG_HI[case_id][triad_label].cPhaseDS['lags'][::1], NP.abs(cpdsval_FG_HI[::1] - cpdsval_FG_HI_predicted[::1]), ls=':', color='magenta')
axs[1].yaxis.tick_right()
axs[1].set_yscale('log')
axs[1].set_xticks(NP.arange(-5,6))
# axs[1].tick_params(axis='x', rotation=45)
axs[1].set_xlim(-3.25, 3.25)
axs[1].set_ylim(1e-5, 1e10)
axs[1].set_xlabel(r'$\tau$ [$\mu$s]', fontsize=11, weight='medium')
axs[1].set_ylabel(r'$\mathcal{V}_\nabla(\tau)$ [pseudo Jy Hz]', fontsize=11, weight='medium')
axs[1].yaxis.set_label_position('right')
fig.subplots_adjust(hspace=0, wspace=0)
fig.subplots_adjust(left=0.15, right=0.85, bottom=0.12, top=0.98)
PLT.savefig(figdir + '{0}_standard_vs_BSP_delay_spectra_{1}_HERA-19.pdf'.format(triad_label, case_id), bbox_inches=0)
In [60]:
spw = [0]
lstind = 0
trind = 0
dayind = 1
colrs = ['red', 'blue', 'forestgreen']
for spwind in spw:
for blind_triplet, triad_label in zip([blind1, blind2], [triad1_label, triad2_label]):
for case_id in case_with_FG_HI:
print('Plotting standard delay spectrum for case "{0}" on triad "{1}"'.format(case_id, triad_label))
fig, axs = PLT.subplots(nrows=1, ncols=2, squeeze=True, sharex=True, figsize=(4,4))
axs[0].axvspan(NP.min(simDPSobjs_FG_HI[case_id].subband_delay_power_spectra['sim']['horizon_kprll_limits'][lstind,spwind,blind_triplet,0]), NP.max(simDPSobjs_FG_HI[case_id].subband_delay_power_spectra['sim']['horizon_kprll_limits'][lstind,spwind,blind_triplet,1]), facecolor='gold')
for blno, blind in enumerate(blind_triplet):
axs[0].plot(simDPSobjs_FG[key_mapping_FG_HI[case_id]['FG']].subband_delay_power_spectra['sim']['kprll'][spwind,::16], 1e6*simDPSobjs_FG[key_mapping_FG_HI[case_id]['FG']].subband_delay_power_spectra['sim']['skyvis_lag'][blind,spwind,::16,lstind], ls='-', color=colrs[blno], alpha=0.5)
axs[0].plot(simDPSobjs_HI[key_mapping_FG_HI[case_id]['HI']].subband_delay_power_spectra['sim']['kprll'][spwind,::16], 1e6*simDPSobjs_HI[key_mapping_FG_HI[case_id]['HI']].subband_delay_power_spectra['sim']['skyvis_lag'][blind,spwind,::16,lstind], ls=':', color=colrs[blno], alpha=0.5)
axs[0].set_yscale('log')
axs[0].set_xlim(-2.5, 2.5)
axs[0].set_ylim(1e-15, 1e15)
axs[0].set_xlabel(r'$k_\parallel$ [$h^{-1}$ Mpc]', fontsize=11, weight='medium')
axs[0].set_ylabel(r'$P_p(k_\parallel)$ [mK$^2\,$(Mpc$/h$)$^3$]', fontsize=11, weight='medium')
print('Plotting closure phase delay spectrum for case "{0}" on triad "{1}"'.format(case_id, triad_label))
cpdpsval_FG = (2/3.0) * xcpdps_FG[key_mapping_FG_HI[case_id]['FG']][triad_label]['oversampled']['whole']['mean'][spwind,lstind,dayind,trind,:].real.to('mK2 Mpc3').value
cpdpsval_FG_HI = (2/3.0) * xcpdps_FG_HI[case_id][triad_label]['oversampled']['whole']['mean'][spwind,lstind,dayind,trind,:].real.to('mK2 Mpc3').value
cpdpsval_FG_HI_predicted = (2/3.0) * xcpdps_FG_HI_predicted[case_id][triad_label]['oversampled']['whole']['mean'][spwind,lstind,dayind,trind,:].real.to('mK2 Mpc3').value
axs[1].axvspan(NP.sum(simDPSobjs_FG_HI[case_id].subband_delay_power_spectra['sim']['horizon_kprll_limits'][lstind,spwind,blind_triplet,0]), NP.sum(simDPSobjs_FG_HI[case_id].subband_delay_power_spectra['sim']['horizon_kprll_limits'][lstind,spwind,blind_triplet,1]), facecolor='gold')
axs[1].plot(xcpdps_FG[key_mapping_FG_HI[case_id]['FG']][triad_label]['oversampled']['kprll'][spwind,::16], cpdpsval_FG[::16], ls='-', color='black')
# axs[1].plot(xcpdps_FG_HI_predicted[case_id][triad_label]['oversampled']['kprll'][spwind,::16], cpdpsval_FG_HI_predicted[::16], ls=':', color='gray')
axs[1].plot(xcpdps_FG_HI[case_id][triad_label]['oversampled']['kprll'][spwind,::16], cpdpsval_FG_HI[::16], ls=':', color='gray')
# if 'cosine' in key_mapping_FG_HI[case_id]['HI']:
axs[1].plot(xcpdps_FG_HI[case_id][triad_label]['oversampled']['kprll'][spwind,::32], NP.abs(cpdpsval_FG_HI[::32] - cpdpsval_FG_HI_predicted[::32]), ls=':', color='magenta')
# print(xcpdps_FG_HI[case_id][triad_label]['oversampled']['kprll'][spwind,2000:3000:16])
# print(cpdpsval_FG_HI[2312:2609] - cpdpsval_FG_HI_predicted[2312:2609])
# print(NP.abs(cpdpsval_FG_HI[2312:2609] - cpdpsval_FG_HI_predicted[2312:2609]).max())
# maxind = 2312 + NP.abs(cpdpsval_FG_HI[2312:2609] - cpdpsval_FG_HI_predicted[2312:2609]).argmax()
# print(maxind)
# print(cpdpsval_FG_HI[maxind], cpdpsval_FG_HI_predicted[maxind])
axs[1].yaxis.tick_right()
axs[1].set_yscale('log')
axs[1].set_xlim(-2.5, 2.5)
axs[1].set_ylim(1e-15, 1e15)
axs[1].set_xlabel(r'$\kappa_\parallel$ [pseudo $h^{-1}$ Mpc]', fontsize=11, weight='medium')
axs[1].set_ylabel(r'$P_\nabla(\kappa_\parallel)$ [pseudo mK$^2\,$(Mpc$/h$)$^3$]', fontsize=11, weight='medium')
axs[1].yaxis.set_label_position('right')
fig.subplots_adjust(hspace=0, wspace=0)
In [61]:
NP.where(NP.logical_and(xcpdps_FG_HI[case_id][triad_label]['oversampled']['kprll'][spwind,:]<-1.0, xcpdps_FG_HI[case_id][triad_label]['oversampled']['kprll'][spwind,:]>-1.2))
Out[61]:
In [ ]: