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

Set up Cosmological parameters


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

Set up input/output directories


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))


('freq window centers = ', array([1.5e+08]))
('freq window BW eff = ', array([42000000.]))
freq window shape = bhw
freq window fftpow = 2.0
pad = 7.0

Set up cases with Foregrounds and HI


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'}
                    }

Set up cases with Foregrounds only


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'}
               }

Set up cases with HI only


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/'}
               }

Load visibilities of cases with Foregorunds only


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----------')


Loading visibilities for case: 1ps_FG_no_spindex_off_zenith
			 from: /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/
Loaded!
----------

Loading visibilities for case: 1ps_FG_spindex_off_zenith
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex/off_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/
Loaded!
----------

Loading visibilities for case: 1ps_FG_no_spindex_on_zenith
			 from: /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/
	Performing minor astrometric corrections...
Loaded!
----------

Loading visibilities for case: 1ps_FG_spindex_on_zenith
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex/on_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/
	Performing minor astrometric corrections...
Loaded!
----------

Loading visibilities for case: GLEAM
			 from: /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/
Loaded!
----------

Load visibilities of cases with HI only


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----------')


Loading visibilities for case: 1ps_HI_cosine
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_HI_cosine/on_zenith_HERA19/UniformBeam_2x60.0sec_1024x97.7kHz/simdata/
	Performing minor astrometric corrections...
Loaded!
----------

Loading visibilities for case: HI_21cmfast_fiducial
			 from: /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/
Loaded!
----------

Load visibilities of cases with Foregorunds and HI


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----------')


Loading visibilities for case: 1ps_FG_no_spindex_HI_colocated
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex_HI_cosine_colocated/
Loaded!
----------

Loading visibilities for case: 1ps_FG_spindex_HI_displaced
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex_HI_cosine_displaced/
Loaded!
----------

Loading visibilities for case: 1ps_FG_spindex_HI_21cmfast
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex_on_zenith_toy_HI_21cmfast_FaintGalaxies_fiducial/
Loaded!
----------

Loading visibilities for case: 1ps_FG_spindex_HI_colocated
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_spindex_HI_cosine_colocated/
Loaded!
----------

Loading visibilities for case: GLEAM_1ps_HI
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/toy_GLEAM_1ps_HI_cosine/
Loaded!
----------

Loading visibilities for case: 1ps_FG_no_spindex_HI_21cmfast
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex_on_zenith_toy_HI_21cmfast_FaintGalaxies_fiducial/
Loaded!
----------

Loading visibilities for case: 1ps_FG_no_spindex_HI_displaced
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/1ps_FG_no_spindex_HI_cosine_displaced/
Loaded!
----------

Loading visibilities for case: GLEAM_HI_21cmfast
			 from: /lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/toy_GLEAM_HI_21cmfast_FaintGalaxies_fiducial_FoV-30/
Loaded!
----------

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())


(30, 1024, 2)
[[100.+0.j 100.+0.j 100.+0.j ... 100.+0.j 100.+0.j 100.+0.j]
 [100.+0.j 100.+0.j 100.+0.j ... 100.+0.j 100.+0.j 100.+0.j]
 [100.+0.j 100.+0.j 100.+0.j ... 100.+0.j 100.+0.j 100.+0.j]]
0.0
(30, 1024, 2)
[[0.01      +0.j 0.00817585+0.j 0.0033689 +0.j ... 0.00266713+0.j
  0.0033689 +0.j 0.00817585+0.j]
 [0.01      +0.j 0.00817585+0.j 0.0033689 +0.j ... 0.00266713+0.j
  0.0033689 +0.j 0.00817585+0.j]
 [0.01      +0.j 0.00817585+0.j 0.0033689 +0.j ... 0.00266713+0.j
  0.0033689 +0.j 0.00817585+0.j]]
0.0

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))


Antenna Info:
[('0', array([-21.9      ,  12.6439709,   0.       ])),
 ('1', array([-7.3      , 12.6439709,  0.       ])),
 ('2', array([ 7.3      , 12.6439709,  0.       ])),
 ('3', array([21.9      , 12.6439709,  0.       ])),
 ('4', array([-21.9      , -12.6439709,   0.       ])),
 ('5', array([ -7.3      , -12.6439709,   0.       ])),
 ('6', array([  7.3      , -12.6439709,   0.       ])),
 ('7', array([ 21.9      , -12.6439709,   0.       ])),
 ('8', array([-14.6       ,  25.28794179,   0.        ])),
 ('9', array([ 0.        , 25.28794179,  0.        ])),
 ('10', array([14.6       , 25.28794179,  0.        ])),
 ('11', array([-14.6       , -25.28794179,   0.        ])),
 ('12', array([  0.        , -25.28794179,   0.        ])),
 ('13', array([ 14.6       , -25.28794179,   0.        ])),
 ('14', array([-2.92000000e+01, -3.41247498e-16,  0.00000000e+00])),
 ('15', array([-1.46000000e+01, -3.41247498e-16,  0.00000000e+00])),
 ('16', array([ 0.00000000e+00, -3.41247498e-16,  0.00000000e+00])),
 ('17', array([ 1.46000000e+01, -3.41247498e-16,  0.00000000e+00])),
 ('18', array([ 2.92000000e+01, -3.41247498e-16,  0.00000000e+00]))]
--------------
Baseline info:
[(('9', '8'), array([14.6,  0. ,  0. ])),
 (('5', '11'), array([ 7.3      , 12.6439709, -0.       ])),
 (('1', '8'), array([  7.3      , -12.6439709,  -0.       ])),
 (('9', '0'), array([21.9      , 12.6439709,  0.       ])),
 (('0', '4'), array([-0.        , 25.28794179, -0.        ])),
 (('16', '0'), array([ 21.9      , -12.6439709,   0.       ])),
 (('10', '8'), array([29.2,  0. ,  0. ])),
 (('16', '11'), array([14.6       , 25.28794179,  0.        ])),
 (('16', '8'), array([ 14.6       , -25.28794179,   0.        ])),
 (('1', '11'), array([ 7.3       , 37.93191269, -0.        ])),
 (('0', '11'), array([-7.3       , 37.93191269, -0.        ])),
 (('3', '8'), array([ 36.5      , -12.6439709,  -0.       ])),
 (('7', '11'), array([36.5      , 12.6439709, -0.       ])),
 (('2', '4'), array([29.2       , 25.28794179, -0.        ])),
 (('6', '0'), array([ 29.2       , -25.28794179,   0.        ])),
 (('2', '11'), array([21.9       , 37.93191269, -0.        ])),
 (('12', '0'), array([ 21.9       , -37.93191269,   0.        ])),
 (('3', '0'), array([43.8,  0. ,  0. ])),
 (('3', '4'), array([43.8       , 25.28794179, -0.        ])),
 (('8', '11'), array([-0.        , 50.57588358, -0.        ])),
 (('7', '0'), array([ 43.8       , -25.28794179,   0.        ])),
 (('3', '11'), array([36.5       , 37.93191269, -0.        ])),
 (('7', '8'), array([ 36.5       , -37.93191269,  -0.        ])),
 (('9', '11'), array([14.6       , 50.57588358, -0.        ])),
 (('8', '12'), array([-14.6       ,  50.57588358,  -0.        ])),
 (('3', '14'), array([51.1      , 12.6439709, -0.       ])),
 (('18', '0'), array([ 51.1      , -12.6439709,   0.       ])),
 (('18', '14'), array([58.4,  0. ,  0. ])),
 (('10', '11'), array([29.2       , 50.57588358, -0.        ])),
 (('13', '8'), array([ 29.2       , -50.57588358,   0.        ]))]

Plot: Annotated antenna positions


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)


Compute closure phases of cases with Foregrounds and HI


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)


Computing closure phases on unique triads for 1ps_FG_no_spindex_HI_colocated...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_spindex_HI_displaced...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_spindex_HI_21cmfast...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_spindex_HI_colocated...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for GLEAM_1ps_HI...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_no_spindex_HI_21cmfast...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_no_spindex_HI_displaced...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for GLEAM_HI_21cmfast...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00

Compute closure phases of cases with Foregrounds only


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)


Computing closure phases on unique triads for 1ps_FG_no_spindex_off_zenith...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_spindex_off_zenith...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_no_spindex_on_zenith...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for 1ps_FG_spindex_on_zenith...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for GLEAM...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00

Compute closure phases of cases with HI only


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)


Computing closure phases on unique triads for 1ps_HI_cosine...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00
Computing closure phases on unique triads for HI_21cmfast_fiducial...
100% |-----------------------------------------| 162/162 Triplets Time: 0:00:00

Predict values for closure phases


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))


Predicted closure phases on unique triads for 1ps_FG_no_spindex_HI_colocated...
Predicted closure phases on unique triads for 1ps_FG_spindex_HI_displaced...
Predicted closure phases on unique triads for 1ps_FG_spindex_HI_21cmfast...
Predicted closure phases on unique triads for 1ps_FG_spindex_HI_colocated...
Predicted closure phases on unique triads for GLEAM_1ps_HI...
Predicted closure phases on unique triads for 1ps_FG_no_spindex_HI_21cmfast...
Predicted closure phases on unique triads for 1ps_FG_no_spindex_HI_displaced...
Predicted closure phases on unique triads for GLEAM_HI_21cmfast...

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'])


(162, 3, 1024, 2)
(162, 1024, 2)
[('0', '1', '2'),
 ('0', '1', '3'),
 ('0', '1', '4'),
 ('0', '1', '5'),
 ('0', '1', '6'),
 ('0', '1', '7'),
 ('0', '1', '8'),
 ('0', '1', '9'),
 ('0', '1', '10'),
 ('0', '1', '11'),
 ('0', '1', '12'),
 ('0', '1', '13'),
 ('0', '1', '16'),
 ('0', '1', '17'),
 ('0', '1', '18'),
 ('0', '2', '4'),
 ('0', '2', '5'),
 ('0', '2', '6'),
 ('0', '2', '7'),
 ('0', '2', '8'),
 ('0', '2', '9'),
 ('0', '2', '10'),
 ('0', '2', '11'),
 ('0', '2', '12'),
 ('0', '2', '13'),
 ('0', '2', '17'),
 ('0', '2', '18'),
 ('0', '3', '4'),
 ('0', '3', '5'),
 ('0', '3', '6'),
 ('0', '3', '7'),
 ('0', '3', '8'),
 ('0', '3', '9'),
 ('0', '3', '10'),
 ('0', '3', '11'),
 ('0', '3', '12'),
 ('0', '3', '13'),
 ('0', '3', '14'),
 ('0', '3', '18'),
 ('0', '4', '8'),
 ('0', '4', '9'),
 ('0', '4', '10'),
 ('0', '4', '11'),
 ('0', '4', '12'),
 ('0', '4', '13'),
 ('0', '4', '14'),
 ('0', '4', '16'),
 ('0', '4', '17'),
 ('0', '4', '18'),
 ('0', '5', '8'),
 ('0', '5', '9'),
 ('0', '5', '10'),
 ('0', '5', '12'),
 ('0', '5', '13'),
 ('0', '5', '14'),
 ('0', '5', '17'),
 ('0', '5', '18'),
 ('0', '6', '8'),
 ('0', '6', '9'),
 ('0', '6', '10'),
 ('0', '6', '13'),
 ('0', '6', '14'),
 ('0', '6', '15'),
 ('0', '6', '18'),
 ('0', '7', '8'),
 ('0', '7', '9'),
 ('0', '7', '10'),
 ('0', '7', '14'),
 ('0', '7', '15'),
 ('0', '7', '16'),
 ('0', '8', '11'),
 ('0', '8', '12'),
 ('0', '8', '13'),
 ('0', '8', '14'),
 ('0', '9', '11'),
 ('0', '9', '12'),
 ('0', '9', '13'),
 ('0', '9', '14'),
 ('0', '9', '15'),
 ('0', '10', '11'),
 ('0', '10', '12'),
 ('0', '10', '13'),
 ('0', '10', '14'),
 ('0', '10', '15'),
 ('0', '10', '16'),
 ('0', '11', '16'),
 ('0', '11', '17'),
 ('0', '11', '18'),
 ('0', '12', '17'),
 ('0', '12', '18'),
 ('0', '13', '18'),
 ('0', '14', '18'),
 ('1', '2', '4'),
 ('1', '2', '11'),
 ('1', '3', '4'),
 ('1', '3', '11'),
 ('1', '3', '14'),
 ('1', '4', '8'),
 ('1', '4', '9'),
 ('1', '4', '10'),
 ('1', '4', '18'),
 ('1', '8', '11'),
 ('1', '8', '12'),
 ('1', '8', '13'),
 ('1', '9', '11'),
 ('1', '10', '11'),
 ('1', '10', '14'),
 ('1', '11', '17'),
 ('1', '11', '18'),
 ('1', '14', '18'),
 ('2', '3', '4'),
 ('2', '3', '11'),
 ('2', '3', '14'),
 ('2', '4', '9'),
 ('2', '4', '10'),
 ('2', '4', '18'),
 ('2', '8', '11'),
 ('2', '8', '12'),
 ('2', '8', '13'),
 ('2', '9', '11'),
 ('2', '10', '11'),
 ('2', '11', '18'),
 ('2', '14', '18'),
 ('3', '4', '10'),
 ('3', '4', '14'),
 ('3', '8', '11'),
 ('3', '8', '12'),
 ('3', '8', '13'),
 ('3', '9', '11'),
 ('3', '10', '11'),
 ('3', '11', '14'),
 ('3', '12', '14'),
 ('3', '13', '14'),
 ('3', '14', '18'),
 ('8', '9', '11'),
 ('8', '9', '12'),
 ('8', '9', '13'),
 ('8', '10', '11'),
 ('8', '10', '12'),
 ('8', '10', '13'),
 ('8', '11', '14'),
 ('8', '11', '15'),
 ('8', '11', '17'),
 ('8', '11', '18'),
 ('8', '12', '14'),
 ('8', '12', '15'),
 ('8', '12', '18'),
 ('8', '13', '14'),
 ('8', '13', '15'),
 ('8', '13', '16'),
 ('8', '14', '18'),
 ('9', '10', '11'),
 ('9', '11', '14'),
 ('9', '11', '15'),
 ('9', '11', '18'),
 ('9', '14', '18'),
 ('10', '11', '14'),
 ('10', '11', '15'),
 ('10', '11', '16'),
 ('10', '14', '18'),
 ('14', '15', '18'),
 ('14', '16', '18')]
[array([[ 14.6,   0. ,   0. ],
       [ 14.6,   0. ,   0. ],
       [-29.2,  -0. ,  -0. ]]),
 array([[ 14.6,   0. ,   0. ],
       [ 29.2,   0. ,   0. ],
       [-43.8,  -0. ,  -0. ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [-14.6       , -25.28794179,  -0.        ],
       [ -0.        ,  25.28794179,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [  0.        , -25.28794179,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [ 14.6       , -25.28794179,   0.        ],
       [-29.2       ,  25.28794179,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [ 29.2       , -25.28794179,   0.        ],
       [-43.8       ,  25.28794179,  -0.        ]]),
 array([[ 14.6      ,   0.       ,   0.       ],
       [ -7.3      ,  12.6439709,   0.       ],
       [ -7.3      , -12.6439709,   0.       ]]),
 array([[ 14.6      ,   0.       ,   0.       ],
       [  7.3      ,  12.6439709,  -0.       ],
       [-21.9      , -12.6439709,  -0.       ]]),
 array([[ 14.6      ,   0.       ,   0.       ],
       [ 21.9      ,  12.6439709,   0.       ],
       [-36.5      , -12.6439709,   0.       ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [ -7.3       , -37.93191269,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [  7.3       , -37.93191269,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [ 21.9       , -37.93191269,   0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[ 14.6      ,   0.       ,   0.       ],
       [  7.3      , -12.6439709,  -0.       ],
       [-21.9      ,  12.6439709,  -0.       ]]),
 array([[ 14.6      ,   0.       ,   0.       ],
       [ 21.9      , -12.6439709,   0.       ],
       [-36.5      ,  12.6439709,   0.       ]]),
 array([[ 14.6      ,   0.       ,   0.       ],
       [ 36.5      , -12.6439709,  -0.       ],
       [-51.1      ,  12.6439709,  -0.       ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [-29.2       , -25.28794179,   0.        ],
       [ -0.        ,  25.28794179,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [-14.6       , -25.28794179,  -0.        ],
       [-14.6       ,  25.28794179,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [  0.        , -25.28794179,   0.        ],
       [-29.2       ,  25.28794179,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [ 14.6       , -25.28794179,   0.        ],
       [-43.8       ,  25.28794179,  -0.        ]]),
 array([[ 29.2      ,   0.       ,   0.       ],
       [-21.9      ,  12.6439709,  -0.       ],
       [ -7.3      , -12.6439709,   0.       ]]),
 array([[ 29.2      ,   0.       ,   0.       ],
       [ -7.3      ,  12.6439709,   0.       ],
       [-21.9      , -12.6439709,  -0.       ]]),
 array([[ 29.2      ,   0.       ,   0.       ],
       [  7.3      ,  12.6439709,  -0.       ],
       [-36.5      , -12.6439709,   0.       ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [-21.9       , -37.93191269,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [ -7.3       , -37.93191269,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [  7.3       , -37.93191269,   0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[ 29.2      ,   0.       ,   0.       ],
       [  7.3      , -12.6439709,  -0.       ],
       [-36.5      ,  12.6439709,   0.       ]]),
 array([[ 29.2      ,   0.       ,   0.       ],
       [ 21.9      , -12.6439709,   0.       ],
       [-51.1      ,  12.6439709,  -0.       ]]),
 array([[ 43.8       ,   0.        ,   0.        ],
       [-43.8       , -25.28794179,   0.        ],
       [ -0.        ,  25.28794179,  -0.        ]]),
 array([[ 43.8       ,   0.        ,   0.        ],
       [-29.2       , -25.28794179,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ]]),
 array([[ 43.8       ,   0.        ,   0.        ],
       [-14.6       , -25.28794179,  -0.        ],
       [-29.2       ,  25.28794179,  -0.        ]]),
 array([[ 43.8       ,   0.        ,   0.        ],
       [  0.        , -25.28794179,   0.        ],
       [-43.8       ,  25.28794179,  -0.        ]]),
 array([[ 43.8      ,   0.       ,   0.       ],
       [-36.5      ,  12.6439709,   0.       ],
       [ -7.3      , -12.6439709,   0.       ]]),
 array([[ 43.8      ,   0.       ,   0.       ],
       [-21.9      ,  12.6439709,  -0.       ],
       [-21.9      , -12.6439709,  -0.       ]]),
 array([[ 43.8      ,   0.       ,   0.       ],
       [ -7.3      ,  12.6439709,   0.       ],
       [-36.5      , -12.6439709,   0.       ]]),
 array([[ 43.8       ,   0.        ,   0.        ],
       [-36.5       , -37.93191269,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[ 43.8       ,   0.        ,   0.        ],
       [-21.9       , -37.93191269,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[ 43.8       ,   0.        ,   0.        ],
       [ -7.3       , -37.93191269,   0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[ 43.8      ,   0.       ,   0.       ],
       [-51.1      , -12.6439709,   0.       ],
       [  7.3      ,  12.6439709,  -0.       ]]),
 array([[ 43.8      ,   0.       ,   0.       ],
       [  7.3      , -12.6439709,  -0.       ],
       [-51.1      ,  12.6439709,  -0.       ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ],
       [ -7.3       , -12.6439709 ,   0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ],
       [-21.9       , -12.6439709 ,  -0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ 36.5       ,  37.93191269,  -0.        ],
       [-36.5       , -12.6439709 ,   0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [  7.3       , -12.6439709 ,  -0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ 21.9       , -12.6439709 ,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ 36.5       , -12.6439709 ,  -0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ -7.3       ,  12.6439709 ,   0.        ],
       [  7.3       ,  12.6439709 ,  -0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ 21.9       ,  12.6439709 ,   0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ 36.5       ,  12.6439709 ,  -0.        ],
       [-36.5       ,  12.6439709 ,   0.        ]]),
 array([[  0.        , -25.28794179,   0.        ],
       [ 51.1       ,  12.6439709 ,  -0.        ],
       [-51.1       ,  12.6439709 ,  -0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ],
       [ -7.3       , -12.6439709 ,   0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ],
       [-21.9       , -12.6439709 ,  -0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ],
       [-36.5       , -12.6439709 ,   0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [  7.3       , -12.6439709 ,  -0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [ 21.9       , -12.6439709 ,   0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ],
       [  7.3       ,  12.6439709 ,  -0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [ 21.9       ,  12.6439709 ,   0.        ],
       [-36.5       ,  12.6439709 ,   0.        ]]),
 array([[ 14.6       , -25.28794179,   0.        ],
       [ 36.5       ,  12.6439709 ,  -0.        ],
       [-51.1       ,  12.6439709 ,  -0.        ]]),
 array([[ 29.2       , -25.28794179,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ],
       [ -7.3       , -12.6439709 ,   0.        ]]),
 array([[ 29.2       , -25.28794179,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ],
       [-21.9       , -12.6439709 ,  -0.        ]]),
 array([[ 29.2       , -25.28794179,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ],
       [-36.5       , -12.6439709 ,   0.        ]]),
 array([[ 29.2       , -25.28794179,   0.        ],
       [  7.3       , -12.6439709 ,  -0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[ 29.2       , -25.28794179,   0.        ],
       [-36.5       ,  12.6439709 ,   0.        ],
       [  7.3       ,  12.6439709 ,  -0.        ]]),
 array([[ 29.2       , -25.28794179,   0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ],
       [ -7.3       ,  12.6439709 ,   0.        ]]),
 array([[ 29.2       , -25.28794179,   0.        ],
       [ 21.9       ,  12.6439709 ,   0.        ],
       [-51.1       ,  12.6439709 ,  -0.        ]]),
 array([[ 43.8       , -25.28794179,   0.        ],
       [-36.5       ,  37.93191269,   0.        ],
       [ -7.3       , -12.6439709 ,   0.        ]]),
 array([[ 43.8       , -25.28794179,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ],
       [-21.9       , -12.6439709 ,  -0.        ]]),
 array([[ 43.8       , -25.28794179,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ],
       [-36.5       , -12.6439709 ,   0.        ]]),
 array([[ 43.8       , -25.28794179,   0.        ],
       [-51.1       ,  12.6439709 ,  -0.        ],
       [  7.3       ,  12.6439709 ,  -0.        ]]),
 array([[ 43.8       , -25.28794179,   0.        ],
       [-36.5       ,  12.6439709 ,   0.        ],
       [ -7.3       ,  12.6439709 ,   0.        ]]),
 array([[ 43.8       , -25.28794179,   0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ]]),
 array([[  7.3       ,  12.6439709 ,  -0.        ],
       [  0.        , -50.57588358,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[  7.3       ,  12.6439709 ,  -0.        ],
       [ 14.6       , -50.57588358,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[  7.3       ,  12.6439709 ,  -0.        ],
       [ 29.2       , -50.57588358,   0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[  7.3       ,  12.6439709 ,  -0.        ],
       [-14.6       , -25.28794179,  -0.        ],
       [  7.3       ,  12.6439709 ,  -0.        ]]),
 array([[ 21.9       ,  12.6439709 ,   0.        ],
       [-14.6       , -50.57588358,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[ 21.9       ,  12.6439709 ,   0.        ],
       [  0.        , -50.57588358,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[ 21.9       ,  12.6439709 ,   0.        ],
       [ 14.6       , -50.57588358,   0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[ 21.9       ,  12.6439709 ,   0.        ],
       [-29.2       , -25.28794179,   0.        ],
       [  7.3       ,  12.6439709 ,  -0.        ]]),
 array([[ 21.9       ,  12.6439709 ,   0.        ],
       [-14.6       , -25.28794179,  -0.        ],
       [ -7.3       ,  12.6439709 ,   0.        ]]),
 array([[ 36.5       ,  12.6439709 ,  -0.        ],
       [-29.2       , -50.57588358,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[ 36.5       ,  12.6439709 ,  -0.        ],
       [-14.6       , -50.57588358,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[ 36.5       ,  12.6439709 ,  -0.        ],
       [  0.        , -50.57588358,   0.        ],
       [-36.5       ,  37.93191269,   0.        ]]),
 array([[ 36.5       ,  12.6439709 ,  -0.        ],
       [-43.8       , -25.28794179,   0.        ],
       [  7.3       ,  12.6439709 ,  -0.        ]]),
 array([[ 36.5       ,  12.6439709 ,  -0.        ],
       [-29.2       , -25.28794179,   0.        ],
       [ -7.3       ,  12.6439709 ,   0.        ]]),
 array([[ 36.5       ,  12.6439709 ,  -0.        ],
       [-14.6       , -25.28794179,  -0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ]]),
 array([[  7.3       , -37.93191269,   0.        ],
       [ 14.6       ,  25.28794179,   0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ]]),
 array([[  7.3       , -37.93191269,   0.        ],
       [ 29.2       ,  25.28794179,  -0.        ],
       [-36.5       ,  12.6439709 ,   0.        ]]),
 array([[  7.3       , -37.93191269,   0.        ],
       [ 43.8       ,  25.28794179,  -0.        ],
       [-51.1       ,  12.6439709 ,  -0.        ]]),
 array([[ 21.9       , -37.93191269,   0.        ],
       [ 14.6       ,  25.28794179,   0.        ],
       [-36.5       ,  12.6439709 ,   0.        ]]),
 array([[ 21.9       , -37.93191269,   0.        ],
       [ 29.2       ,  25.28794179,  -0.        ],
       [-51.1       ,  12.6439709 ,  -0.        ]]),
 array([[ 36.5       , -37.93191269,  -0.        ],
       [ 14.6       ,  25.28794179,   0.        ],
       [-51.1       ,  12.6439709 ,  -0.        ]]),
 array([[ -7.3      , -12.6439709,   0.       ],
       [ 58.4      ,   0.       ,   0.       ],
       [-51.1      ,  12.6439709,  -0.       ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [-29.2       , -25.28794179,   0.        ],
       [ 14.6       ,  25.28794179,   0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [-21.9       , -37.93191269,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [-43.8       , -25.28794179,   0.        ],
       [ 14.6       ,  25.28794179,   0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [-36.5       , -37.93191269,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ]]),
 array([[ 29.2      ,   0.       ,   0.       ],
       [-51.1      , -12.6439709,   0.       ],
       [ 21.9      ,  12.6439709,   0.       ]]),
 array([[-14.6       , -25.28794179,  -0.        ],
       [  7.3       ,  37.93191269,  -0.        ],
       [  7.3       , -12.6439709 ,  -0.        ]]),
 array([[-14.6       , -25.28794179,  -0.        ],
       [ 21.9       ,  37.93191269,  -0.        ],
       [ -7.3       , -12.6439709 ,   0.        ]]),
 array([[-14.6       , -25.28794179,  -0.        ],
       [ 36.5       ,  37.93191269,  -0.        ],
       [-21.9       , -12.6439709 ,  -0.        ]]),
 array([[-14.6       , -25.28794179,  -0.        ],
       [ 51.1       ,  12.6439709 ,  -0.        ],
       [-36.5       ,  12.6439709 ,   0.        ]]),
 array([[ -7.3       ,  12.6439709 ,   0.        ],
       [  0.        , -50.57588358,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ]]),
 array([[ -7.3       ,  12.6439709 ,   0.        ],
       [ 14.6       , -50.57588358,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[ -7.3       ,  12.6439709 ,   0.        ],
       [ 29.2       , -50.57588358,   0.        ],
       [-21.9       ,  37.93191269,  -0.        ]]),
 array([[  7.3       ,  12.6439709 ,  -0.        ],
       [-14.6       , -50.57588358,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ]]),
 array([[ 21.9       ,  12.6439709 ,   0.        ],
       [-29.2       , -50.57588358,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ]]),
 array([[ 21.9       ,  12.6439709 ,   0.        ],
       [-43.8       , -25.28794179,   0.        ],
       [ 21.9       ,  12.6439709 ,   0.        ]]),
 array([[ -7.3       , -37.93191269,   0.        ],
       [ 29.2       ,  25.28794179,  -0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ]]),
 array([[ -7.3       , -37.93191269,   0.        ],
       [ 43.8       ,  25.28794179,  -0.        ],
       [-36.5       ,  12.6439709 ,   0.        ]]),
 array([[-21.9      , -12.6439709,  -0.       ],
       [ 58.4      ,   0.       ,   0.       ],
       [-36.5      ,  12.6439709,   0.       ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [-43.8       , -25.28794179,   0.        ],
       [ 29.2       ,  25.28794179,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [-36.5       , -37.93191269,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ]]),
 array([[ 14.6      ,   0.       ,   0.       ],
       [-51.1      , -12.6439709,   0.       ],
       [ 36.5      ,  12.6439709,  -0.       ]]),
 array([[-29.2       , -25.28794179,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ],
       [  7.3       , -12.6439709 ,  -0.        ]]),
 array([[-29.2       , -25.28794179,   0.        ],
       [ 36.5       ,  37.93191269,  -0.        ],
       [ -7.3       , -12.6439709 ,   0.        ]]),
 array([[-29.2       , -25.28794179,   0.        ],
       [ 51.1       ,  12.6439709 ,  -0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ]]),
 array([[-21.9       ,  12.6439709 ,  -0.        ],
       [  0.        , -50.57588358,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ]]),
 array([[-21.9       ,  12.6439709 ,  -0.        ],
       [ 14.6       , -50.57588358,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ]]),
 array([[-21.9       ,  12.6439709 ,  -0.        ],
       [ 29.2       , -50.57588358,   0.        ],
       [ -7.3       ,  37.93191269,  -0.        ]]),
 array([[ -7.3       ,  12.6439709 ,   0.        ],
       [-14.6       , -50.57588358,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ]]),
 array([[  7.3       ,  12.6439709 ,  -0.        ],
       [-29.2       , -50.57588358,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ]]),
 array([[-21.9       , -37.93191269,   0.        ],
       [ 43.8       ,  25.28794179,  -0.        ],
       [-21.9       ,  12.6439709 ,  -0.        ]]),
 array([[-36.5      , -12.6439709,   0.       ],
       [ 58.4      ,   0.       ,   0.       ],
       [-21.9      ,  12.6439709,  -0.       ]]),
 array([[-43.8       , -25.28794179,   0.        ],
       [ 36.5       ,  37.93191269,  -0.        ],
       [  7.3       , -12.6439709 ,  -0.        ]]),
 array([[-43.8       , -25.28794179,   0.        ],
       [ -7.3       ,  12.6439709 ,   0.        ],
       [ 51.1       ,  12.6439709 ,  -0.        ]]),
 array([[-36.5       ,  12.6439709 ,   0.        ],
       [  0.        , -50.57588358,   0.        ],
       [ 36.5       ,  37.93191269,  -0.        ]]),
 array([[-36.5       ,  12.6439709 ,   0.        ],
       [ 14.6       , -50.57588358,   0.        ],
       [ 21.9       ,  37.93191269,  -0.        ]]),
 array([[-36.5       ,  12.6439709 ,   0.        ],
       [ 29.2       , -50.57588358,   0.        ],
       [  7.3       ,  37.93191269,  -0.        ]]),
 array([[-21.9       ,  12.6439709 ,  -0.        ],
       [-14.6       , -50.57588358,   0.        ],
       [ 36.5       ,  37.93191269,  -0.        ]]),
 array([[ -7.3       ,  12.6439709 ,   0.        ],
       [-29.2       , -50.57588358,   0.        ],
       [ 36.5       ,  37.93191269,  -0.        ]]),
 array([[-36.5       , -37.93191269,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ],
       [ 51.1       ,  12.6439709 ,  -0.        ]]),
 array([[-21.9       , -37.93191269,   0.        ],
       [-29.2       ,  25.28794179,  -0.        ],
       [ 51.1       ,  12.6439709 ,  -0.        ]]),
 array([[ -7.3       , -37.93191269,   0.        ],
       [-43.8       ,  25.28794179,  -0.        ],
       [ 51.1       ,  12.6439709 ,  -0.        ]]),
 array([[-51.1      , -12.6439709,   0.       ],
       [ 58.4      ,   0.       ,   0.       ],
       [ -7.3      ,  12.6439709,   0.       ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [-14.6       , -50.57588358,   0.        ],
       [ -0.        ,  50.57588358,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [  0.        , -50.57588358,   0.        ],
       [-14.6       ,  50.57588358,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [ 14.6       , -50.57588358,   0.        ],
       [-29.2       ,  50.57588358,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [-29.2       , -50.57588358,   0.        ],
       [ -0.        ,  50.57588358,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [-14.6       , -50.57588358,   0.        ],
       [-14.6       ,  50.57588358,  -0.        ]]),
 array([[ 29.2       ,   0.        ,   0.        ],
       [  0.        , -50.57588358,   0.        ],
       [-29.2       ,  50.57588358,  -0.        ]]),
 array([[  0.        , -50.57588358,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ],
       [ 14.6       ,  25.28794179,   0.        ]]),
 array([[  0.        , -50.57588358,   0.        ],
       [ -0.        ,  25.28794179,  -0.        ],
       [ -0.        ,  25.28794179,  -0.        ]]),
 array([[  0.        , -50.57588358,   0.        ],
       [ 29.2       ,  25.28794179,  -0.        ],
       [-29.2       ,  25.28794179,  -0.        ]]),
 array([[  0.        , -50.57588358,   0.        ],
       [ 43.8       ,  25.28794179,  -0.        ],
       [-43.8       ,  25.28794179,  -0.        ]]),
 array([[ 14.6       , -50.57588358,   0.        ],
       [-29.2       ,  25.28794179,  -0.        ],
       [ 14.6       ,  25.28794179,   0.        ]]),
 array([[ 14.6       , -50.57588358,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ],
       [ -0.        ,  25.28794179,  -0.        ]]),
 array([[ 14.6       , -50.57588358,   0.        ],
       [ 29.2       ,  25.28794179,  -0.        ],
       [-43.8       ,  25.28794179,  -0.        ]]),
 array([[ 29.2       , -50.57588358,   0.        ],
       [-43.8       ,  25.28794179,  -0.        ],
       [ 14.6       ,  25.28794179,   0.        ]]),
 array([[ 29.2       , -50.57588358,   0.        ],
       [-29.2       ,  25.28794179,  -0.        ],
       [ -0.        ,  25.28794179,  -0.        ]]),
 array([[ 29.2       , -50.57588358,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ],
       [-14.6       ,  25.28794179,  -0.        ]]),
 array([[-14.6       , -25.28794179,  -0.        ],
       [ 58.4       ,   0.        ,   0.        ],
       [-43.8       ,  25.28794179,  -0.        ]]),
 array([[ 14.6       ,   0.        ,   0.        ],
       [-29.2       , -50.57588358,   0.        ],
       [ 14.6       ,  50.57588358,  -0.        ]]),
 array([[-14.6       , -50.57588358,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ],
       [ 29.2       ,  25.28794179,  -0.        ]]),
 array([[-14.6       , -50.57588358,   0.        ],
       [ -0.        ,  25.28794179,  -0.        ],
       [ 14.6       ,  25.28794179,   0.        ]]),
 array([[-14.6       , -50.57588358,   0.        ],
       [ 43.8       ,  25.28794179,  -0.        ],
       [-29.2       ,  25.28794179,  -0.        ]]),
 array([[-29.2       , -25.28794179,   0.        ],
       [ 58.4       ,   0.        ,   0.        ],
       [-29.2       ,  25.28794179,  -0.        ]]),
 array([[-29.2       , -50.57588358,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ],
       [ 43.8       ,  25.28794179,  -0.        ]]),
 array([[-29.2       , -50.57588358,   0.        ],
       [ -0.        ,  25.28794179,  -0.        ],
       [ 29.2       ,  25.28794179,  -0.        ]]),
 array([[-29.2       , -50.57588358,   0.        ],
       [ 14.6       ,  25.28794179,   0.        ],
       [ 14.6       ,  25.28794179,   0.        ]]),
 array([[-43.8       , -25.28794179,   0.        ],
       [ 58.4       ,   0.        ,   0.        ],
       [-14.6       ,  25.28794179,  -0.        ]]),
 array([[ 14.6,   0. ,   0. ],
       [ 43.8,   0. ,   0. ],
       [-58.4,  -0. ,  -0. ]]),
 array([[ 29.2,   0. ,   0. ],
       [ 29.2,   0. ,   0. ],
       [-58.4,  -0. ,  -0. ]])]

In [20]:
print(cpinfo_case_with_HI.keys())


['1ps_HI_cosine', 'HI_21cmfast_fiducial']

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)


6
[[ 14.6         0.          0.       ]
 [ -7.3        12.6439709   0.       ]
 [ -7.3       -12.6439709   0.       ]]
[0, 2, 1]
[[ 14.6         0.          0.       ]
 [  7.3       -12.6439709  -0.       ]
 [  7.3        12.6439709  -0.       ]]

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)


143
[[  0.         -50.57588358   0.        ]
 [ 43.8         25.28794179  -0.        ]
 [-43.8         25.28794179  -0.        ]]
[19, 18, 20]
[[ -0.          50.57588358  -0.        ]
 [ 43.8         25.28794179  -0.        ]
 [ 43.8        -25.28794179   0.        ]]

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)


6
('0', '1', '8')
(1024,)
(162, 3, 1024, 2)
(162, 1024, 2)
(162, 1024, 2)
(162, 1024, 2)
(162, 1024, 2)

In [24]:
print(figdir + 'closure_phase_comparison_actual_vs_prediction_HERA-19_triad_{0[0]}-{0[1]}-{0[2]}.pdf'.format(triad1))


/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/figures/closure_phase_comparison_actual_vs_prediction_HERA-19_triad_0-1-8.pdf

Plot HI visibility amplitudes


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)


	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)

Plot: Plot Foreground and fluctuations in amplitude due to HI


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)


Plotting 1ps_FG_no_spindex_HI_colocated:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_spindex_HI_displaced:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_spindex_HI_21cmfast:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_spindex_HI_colocated:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting GLEAM_1ps_HI:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_no_spindex_HI_21cmfast:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_no_spindex_HI_displaced:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting GLEAM_HI_21cmfast:
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)

Plot: Compare computed and predicted values for closure phases


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)


Plotting 1ps_FG_no_spindex_HI_colocated
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_spindex_HI_displaced
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_spindex_HI_21cmfast
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_spindex_HI_colocated
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting GLEAM_1ps_HI
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_no_spindex_HI_21cmfast
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting 1ps_FG_no_spindex_HI_displaced
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)
Plotting GLEAM_HI_21cmfast
	Plotting Triad:(0,1,8)
	Plotting Triad:(8,11,18)

Initialize delay spectrum objects for simulated visibilities


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'])

Compute standard delay spectra for simulated visibilities in subbands


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('-----------------')


Foregrounds only
================
Computing standard subband delay transform for case "1ps_FG_no_spindex_off_zenith"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_spindex_off_zenith"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/scipy/signal/signaltools.py:2223: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  Y[sl] = X[sl]
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/scipy/signal/signaltools.py:2225: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  Y[sl] = X[sl]
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/scipy/signal/signaltools.py:2230: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  Y[sl] += X[sl]  # add the component of X at N/2
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_no_spindex_on_zenith"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_spindex_on_zenith"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "GLEAM"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------

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('-----------------')


HI only
=======
Computing standard subband delay transform for case "1ps_HI_cosine"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "HI_21cmfast_fiducial"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------

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('-----------------')


Foregrounds + HI
================
Computing standard subband delay transform for case "1ps_FG_no_spindex_HI_colocated"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_spindex_HI_displaced"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_spindex_HI_21cmfast"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_spindex_HI_colocated"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "GLEAM_1ps_HI"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_no_spindex_HI_21cmfast"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "1ps_FG_no_spindex_HI_displaced"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------
Computing standard subband delay transform for case "GLEAM_HI_21cmfast"-----------------
Preparing to compute delay transform...
	Checking input parameters for compatibility...
	Frequency window weights assigned.
	Input parameters have been verified to be compatible.
	Proceeding to compute delay transform.
	Delay transform computed with padding fraction 1.0
	Delay transform products downsampled by factor of 2.0
delay_transform() completed successfully.
	Renormalized the shaping window to unit power.
	Sub-band(s) delay transform computed
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	Downsampled Sub-band(s) delay transform computed
-----------------

Compute standard delay power spectra for simulated visibilities in subbands


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('-----------------')


Foregrounds only
================
Initializing and computing standard subband delay power spectra for case "1ps_FG_no_spindex_off_zenith"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_spindex_off_zenith"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_no_spindex_on_zenith"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_spindex_on_zenith"
-----------------
Initializing and computing standard subband delay power spectra for case "GLEAM"
-----------------

HI only
=======
Initializing and computing standard subband delay power spectra for case "1ps_HI_cosine"
-----------------
Initializing and computing standard subband delay power spectra for case "HI_21cmfast_fiducial"
-----------------

Foregrounds + HI
================
Initializing and computing standard subband delay power spectra for case "1ps_FG_no_spindex_HI_colocated"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_spindex_HI_displaced"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_spindex_HI_21cmfast"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_spindex_HI_colocated"
-----------------
Initializing and computing standard subband delay power spectra for case "GLEAM_1ps_HI"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_no_spindex_HI_21cmfast"
-----------------
Initializing and computing standard subband delay power spectra for case "1ps_FG_no_spindex_HI_displaced"
-----------------
Initializing and computing standard subband delay power spectra for case "GLEAM_HI_21cmfast"
-----------------

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)


['drz_los', 'kprll', 'kperp', 'jacobian2', 'rz_transverse', 'vis_lag', 'vis_noise_lag', 'dz', 'factor', 'skyvis_lag', 'horizon_kprll_limits', 'Jy2K', 'z', 'rz_los', 'jacobian1']
(2, 1, 30, 2)

Initialize instances of class ClosurePhase (no approximation)


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)


['1ps_FG_spindex_HI_displaced',
 '1ps_FG_no_spindex_HI_displaced',
 'GLEAM_HI_21cmfast',
 '1ps_FG_no_spindex_HI_colocated',
 '1ps_FG_spindex_HI_colocated',
 'GLEAM_1ps_HI',
 '1ps_FG_spindex_HI_21cmfast',
 '1ps_FG_no_spindex_HI_21cmfast']
-------------
['1ps_FG_no_spindex_off_zenith',
 '1ps_FG_spindex_off_zenith',
 '1ps_FG_no_spindex_on_zenith',
 '1ps_FG_spindex_on_zenith',
 'GLEAM']
-------------
[u'raw', 'processed', 'errinfo']
-------------
[u'triads', u'cphase', u'days', u'lst-day', u'lst', u'flags']
(1, 3)
[['8' '11' '18']]
(4,)
(2, 4)
(2, 4)
(2, 4, 1, 1024)

Initialize instances of class ClosurePhase predictions (from first order approximation)


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('-------------')


Working on creating closure phase objects with first-order predictions for case: 1ps_FG_no_spindex_HI_colocated
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------
Working on creating closure phase objects with first-order predictions for case: 1ps_FG_spindex_HI_displaced
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------
Working on creating closure phase objects with first-order predictions for case: 1ps_FG_spindex_HI_21cmfast
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------
Working on creating closure phase objects with first-order predictions for case: 1ps_FG_spindex_HI_colocated
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------
Working on creating closure phase objects with first-order predictions for case: GLEAM_1ps_HI
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------
Working on creating closure phase objects with first-order predictions for case: 1ps_FG_no_spindex_HI_21cmfast
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------
Working on creating closure phase objects with first-order predictions for case: 1ps_FG_no_spindex_HI_displaced
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------
Working on creating closure phase objects with first-order predictions for case: GLEAM_HI_21cmfast
	EQ14
(2, 4, 1, 1024)
	EQ50
(2, 4, 1, 1024)
-------------

Smooth in LST and days


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)


('day bin size = None', 'ndaybins=2', 'LST bin size = 59.9s')
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/prisim/bispectrum_phase.py:1870: UserWarning: LST bin size found to be smaller than the LST resolution in the data. No LST binning/averaging will be performed.
  warnings.warn('LST bin size found to be smaller than the LST resolution in the data. No LST binning/averaging will be performed.')
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/scipy/stats/_binned_statistic.py:607: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  result = result[core]

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)


['prelim', 'native']
--------------
['cphase', 'wts', 'eicp', 'dlstbins', 'daybins', 'diff_dbins', 'lstbins']
['rms', 'median', 'mad', 'mean']
(2,)
[5.55625446 5.57296676]
(2,)
[2458446. 2458448.]
(2, 2, 1, 1024)
(2, 2, 1, 1024)

Subtract a model of Closure Phase (optional)

Creates new keys 'submodel' and 'residual' with a key-value structure similar to 'prelim'


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())


['prelim', 'residual', 'submodel', 'native']

Create subsample differences to keep track of noise from the noisy data


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)


('ndaybins=4', 'LST bin size = 59.9s')
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/prisim/bispectrum_phase.py:2142: UserWarning: LST bin size found to be smaller than the LST resolution in the data. No LST binning/averaging will be performed.
  warnings.warn('LST bin size found to be smaller than the LST resolution in the data. No LST binning/averaging will be performed.')

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)


['wts', 'eicp_diff', 'dlstbins', 'daybins', 'lstbins', 'list_of_pair_of_pairs', 'diff_dbins']
(4,)
(2,)
3
[[0, 1, 2, 3], [0, 2, 1, 3], [0, 3, 1, 2]]
['1', '0']
['median', 'mean']
(2, 3, 1, 1024)
(2, 3, 1, 1024)

Create instances of class ClosurePhaseDelaySpectrum for each of the cases


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])

Prepare to compute delay spectrum of the ClosurePhaseDelaySpectrum instances


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]])


{'EQ14': {'bltriplet': array([[ 14.6      ,   0.       ,   0.       ],
       [ -7.3      ,  12.6439709,   0.       ],
       [ -7.3      , -12.6439709,   0.       ]]),
          'smoothinfo': {'interp_kind': 'linear', 'op_type': 'interp1d'},
          'vis': <prisim.interferometry.InterferometerArray object at 0x7fc2f69a7510>},
 'EQ50': {'bltriplet': array([[  0.        , -50.57588358,   0.        ],
       [ 43.8       ,  25.28794179,  -0.        ],
       [-43.8       ,  25.28794179,  -0.        ]]),
          'smoothinfo': {'interp_kind': 'linear', 'op_type': 'interp1d'},
          'vis': <prisim.interferometry.InterferometerArray object at 0x7fc2f69a7510>}}

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))


('freq window centers = ', array([1.5e+08]))
('freq window BW eff = ', array([42000000.]))
freq window shape = bhw
freq window fftpow = 2.0
pad = 7.0

Compute delay spectrum of the ClosurePhaseDelaySpectrum instance

Foregrounds only


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('-------------')


Foregrounds only
================
Computing delay spectrum for case: 1ps_FG_no_spindex_off_zenith
	EQ14
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/numpy/ma/core.py:2788: ComplexWarning: Casting complex values to real discards the imaginary part
  order=order, subok=True, ndmin=ndmin)
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_spindex_off_zenith
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_no_spindex_on_zenith
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_spindex_on_zenith
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: GLEAM
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------

Foregrounds + HI (no approximation)


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('-------------')


Foregrounds + HI (No approximation)
===================================
Computing delay spectrum for case: 1ps_FG_spindex_HI_displaced
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_no_spindex_HI_displaced
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_no_spindex_HI_21cmfast
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_no_spindex_HI_colocated
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_spindex_HI_colocated
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: GLEAM_1ps_HI
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_spindex_HI_21cmfast
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: GLEAM_HI_21cmfast
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------

Foregrounds + HI (predicted based on linear-order approximation)


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('-------------')


Foregrounds + HI (linear-order prediction)
==========================================
Computing delay spectrum for case: 1ps_FG_spindex_HI_displaced
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_no_spindex_HI_displaced
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: GLEAM_HI_21cmfast
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_no_spindex_HI_colocated
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_spindex_HI_colocated
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: GLEAM_1ps_HI
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_spindex_HI_21cmfast
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------
Computing delay spectrum for case: 1ps_FG_no_spindex_HI_21cmfast
	EQ14
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
	EQ50
	Renormalized the shaping window to unit power.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Determining the interpolating function for downsampling.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
Returning the downsampled data.
-------------

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())


(1, 2, 2, 1, 8192)
1976417599.1726859

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)


{'EQ14': None, 'EQ50': None}

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('-------------')


Foregrounds only
================
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_off_zenith
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_off_zenith
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_on_zenith
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_on_zenith
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: GLEAM
	EQ14
	EQ50
-------------

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('-------------')


Foregrounds + HI (No approximation)
===================================
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_HI_colocated
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_HI_displaced
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_HI_21cmfast
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_HI_colocated
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: GLEAM_1ps_HI
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_HI_21cmfast
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_HI_displaced
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: GLEAM_HI_21cmfast
	EQ14
	EQ50
-------------

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('-------------')


Foregrounds + HI (linear-order approximation)
==============================================
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_HI_colocated
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_HI_displaced
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_HI_21cmfast
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_spindex_HI_colocated
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: GLEAM_1ps_HI
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_HI_21cmfast
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: 1ps_FG_no_spindex_HI_displaced
	EQ14
	EQ50
-------------
Computing closure phase delay power spectrum for case: GLEAM_HI_21cmfast
	EQ14
	EQ50
-------------

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'])


['triads', 'dday', 'triads_ind', 'lst', 'lstXoffsets', 'lst_ind', 'dlst', 'days', 'resampled', 'oversampled', 'day_ind']
[['0' '1' '8']]
['lags', 'residual', 'kprll', 'freq_wts', 'freq_center', 'bw_eff', 'shape', 'lag_corr_length', 'z', 'whole', 'submodel']
['axesmap', 'median', 'nsamples_coh', 'diagoffsets', 'nsamples_incoh', 'diagweights', 'mean']
(1, 2, 3, 1, 8192)
{1: array([0, 1]), 2: array([-1,  0,  1]), 3: array([0])}

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('-----------------')


Plotting deviations between actual and predicted closure phase power spectra
1ps_FG_no_spindex_HI_colocated:
	EQ14
	EQ50
-----------------
1ps_FG_spindex_HI_displaced:
	EQ14
	EQ50
-----------------
1ps_FG_spindex_HI_21cmfast:
	EQ14
	EQ50
-----------------
1ps_FG_spindex_HI_colocated:
	EQ14
	EQ50
-----------------
GLEAM_1ps_HI:
	EQ14
	EQ50
-----------------
1ps_FG_no_spindex_HI_21cmfast:
	EQ14
	EQ50
-----------------
1ps_FG_no_spindex_HI_displaced:
	EQ14
	EQ50
-----------------
GLEAM_HI_21cmfast:
	EQ14
	EQ50
-----------------

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])

Plot standard and BSP delay spectra and for FG + HI cases as a function of $\tau$


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)


Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ14"
Plotting standard delay spectrum for case "GLEAM_1ps_HI" on triad "EQ14"
Plotting closure phase delay spectrum for case "GLEAM_1ps_HI" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ14"
Plotting standard delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ14"
Plotting closure phase delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ50"
Plotting standard delay spectrum for case "GLEAM_1ps_HI" on triad "EQ50"
Plotting closure phase delay spectrum for case "GLEAM_1ps_HI" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ50"
Plotting standard delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ50"
Plotting closure phase delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ50"

Plot standard and BSP delay power spectra and for FG + HI cases as a function of $k_\parallel$


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)


Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ14"
Plotting standard delay spectrum for case "GLEAM_1ps_HI" on triad "EQ14"
Plotting closure phase delay spectrum for case "GLEAM_1ps_HI" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ14"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ14"
Plotting standard delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ14"
Plotting closure phase delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ14"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_colocated" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_displaced" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_21cmfast" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_spindex_HI_colocated" on triad "EQ50"
Plotting standard delay spectrum for case "GLEAM_1ps_HI" on triad "EQ50"
Plotting closure phase delay spectrum for case "GLEAM_1ps_HI" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_21cmfast" on triad "EQ50"
Plotting standard delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ50"
Plotting closure phase delay spectrum for case "1ps_FG_no_spindex_HI_displaced" on triad "EQ50"
Plotting standard delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ50"
Plotting closure phase delay spectrum for case "GLEAM_HI_21cmfast" on triad "EQ50"

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]:
(array([2312, 2313, 2314, 2315, 2316, 2317, 2318, 2319, 2320, 2321, 2322,
        2323, 2324, 2325, 2326, 2327, 2328, 2329, 2330, 2331, 2332, 2333,
        2334, 2335, 2336, 2337, 2338, 2339, 2340, 2341, 2342, 2343, 2344,
        2345, 2346, 2347, 2348, 2349, 2350, 2351, 2352, 2353, 2354, 2355,
        2356, 2357, 2358, 2359, 2360, 2361, 2362, 2363, 2364, 2365, 2366,
        2367, 2368, 2369, 2370, 2371, 2372, 2373, 2374, 2375, 2376, 2377,
        2378, 2379, 2380, 2381, 2382, 2383, 2384, 2385, 2386, 2387, 2388,
        2389, 2390, 2391, 2392, 2393, 2394, 2395, 2396, 2397, 2398, 2399,
        2400, 2401, 2402, 2403, 2404, 2405, 2406, 2407, 2408, 2409, 2410,
        2411, 2412, 2413, 2414, 2415, 2416, 2417, 2418, 2419, 2420, 2421,
        2422, 2423, 2424, 2425, 2426, 2427, 2428, 2429, 2430, 2431, 2432,
        2433, 2434, 2435, 2436, 2437, 2438, 2439, 2440, 2441, 2442, 2443,
        2444, 2445, 2446, 2447, 2448, 2449, 2450, 2451, 2452, 2453, 2454,
        2455, 2456, 2457, 2458, 2459, 2460, 2461, 2462, 2463, 2464, 2465,
        2466, 2467, 2468, 2469, 2470, 2471, 2472, 2473, 2474, 2475, 2476,
        2477, 2478, 2479, 2480, 2481, 2482, 2483, 2484, 2485, 2486, 2487,
        2488, 2489, 2490, 2491, 2492, 2493, 2494, 2495, 2496, 2497, 2498,
        2499, 2500, 2501, 2502, 2503, 2504, 2505, 2506, 2507, 2508, 2509,
        2510, 2511, 2512, 2513, 2514, 2515, 2516, 2517, 2518, 2519, 2520,
        2521, 2522, 2523, 2524, 2525, 2526, 2527, 2528, 2529, 2530, 2531,
        2532, 2533, 2534, 2535, 2536, 2537, 2538, 2539, 2540, 2541, 2542,
        2543, 2544, 2545, 2546, 2547, 2548, 2549, 2550, 2551, 2552, 2553,
        2554, 2555, 2556, 2557, 2558, 2559, 2560, 2561, 2562, 2563, 2564,
        2565, 2566, 2567, 2568, 2569, 2570, 2571, 2572, 2573, 2574, 2575,
        2576, 2577, 2578, 2579, 2580, 2581, 2582, 2583, 2584, 2585, 2586,
        2587, 2588, 2589, 2590, 2591, 2592, 2593, 2594, 2595, 2596, 2597,
        2598, 2599, 2600, 2601, 2602, 2603, 2604, 2605, 2606, 2607, 2608,
        2609]),)

In [ ]: