In [1]:

    

import copy, glob
import progressbar as PGB
import h5py
import numpy as NP
import numpy.ma as MA
from scipy import interpolate, stats
import matplotlib.pyplot as PLT
import matplotlib.colors as PLTC
import matplotlib.ticker as PLTick
import yaml, argparse, warnings
from astropy.io import ascii
import astropy.units as U
from astropy.stats import sigma_clipped_stats, histogram
from astroutils import DSP_modules as DSP
from astroutils import constants as CNST
from astroutils import mathops as OPS
from astroutils import nonmathops as NMO
from astroutils import lookup_operations as LKP
import astroutils
import prisim
from prisim import interferometry as RI
from prisim import bispectrum_phase as BSP
from prisim import delay_spectrum as DS
from IPython.core.debugger import set_trace
%matplotlib inline


    

In [6]:

    

print('AstroUtils git # {0}\nPRISim git # {1}'.format(astroutils.__githash__, prisim.__githash__))


    

    




AstroUtils git # 1be75833b358cbc53ff9fd12ec376a5c6335463a
PRISim git # bb083bc1abd3ef44d0b18c094fedfe3dc30ed6ac

In [5]:

    

inparmsfile = '/lustre/aoc/users/nthyagar/codes/mine/python/projects/closure/multiday_EQ28_model_Fornax_transit_closure_PS_analysis_parms.yaml'
with open(inparmsfile, 'r') as parms_file:
    parms = yaml.safe_load(parms_file)


    

In [6]:

    

print(parms)


    

    




{'plot': {'1g': {'action': False, 'infile1': 'EQ14new', 'infile2': 'EQ14old', 'avgind': 0, 'obsind': 0}, '1f': {'action': False, 'labels': ['EWs', 'EQ14', 'EQ29'], 'obsind': 0, 'infiles': ['EW14', 'EQ14', 'EQ28'], 'avgind': 0}, '3b': {'action': False}, '1b': {'action': False, 'triad': [0, 1, 12]}, '1': {'action': False}, '2d': {'action': True, 'kbin_min': None, 'kbin_max': None, 'kbintype': 'linear', 'num_kbins': None}, '2c_err': {'action': True, 'other_pserr_file': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/Binned_Data/EQ28/XX/PS/xCPDPS_collapse_axes_123_EQ28XX_errinfo_incoh_diag_avg.hdf5', 'current_lbl': 'Model', 'other_lbl': 'Data'}, '2f': {'collapseax': [1], 'incohax': [1, 2, 3], 'selection': {'triads': None, 'lst': [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52], 'dlst_range': [1.1], 'days': None}, 'datapool': 'whole', 'antloc_file': '/data4/HERA_Simulations/simulations/hera_layouts/antenna_positions_350.dat', 'cohax': None, 'subselection': {'spw': None, 'day': [[0, 0], [0, 1]], 'lstdiag': [0, 1]}, 'statistic': 'median', 'action': False}, '1a': {'action': False, 'triad': None}, '2e': {'action': False, 'subselection': {'triaddiag': [0], 'spw': [0, 1], 'day': [[0, 1]], 'lstdiag': [0, 1]}}, '3': {'21cm_PS_dir': '/lustre/aoc/users/nthyagar/data/EoR_models/21cmFAST/Andrei_Mesinger/Faint_galaxies_fiducial_1024/PS/21cm_PS/', 'simlabels': ['FG', 'HI'], 'sim_rootdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/', 'action': True, 'visdirs': ['Fornax_GLEAM_all_FOV30_ephemeris_HERA61/Fornax_transit_NFbeam_120x10.7sec_1024x97.7kHz/simdata/', 'HI_21cmfast_FaintGalaxies_fiducial/14.0arcmin_v2/Fornax_transit_NFbeam_120x10.7sec_1024x97.7kHz/simdata/'], 'sampling': 'resampled', 'visfile_prfx': 'all-simvis'}, '2': {'incohax': [1, 2, 3], 'selection': {'lstrange': [3.0, 3.6], 'triads': None, 'lst': [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52], 'dlst_range': [1.1], 'days': None}, 'PS_dir': 'PS/', 'xlim': None, 'beaminfo': {'filefmt': 'UVbeam', 'telescope': {'shape': 'dish', 'ocoords': 'altaz', 'orientation': [90.0, 270.0], 'ground_plane': None, 'size': 14.0, 'id': 'hera', 'phased_array': False}, 'nside': 128, 'select_freq': 150000000.0, 'chromatic': True, 'pol': 'X', 'beamfile': 'NF_HERA_antenna_power_pattern_99-201_MHz_nside_128.uvbeam', 'filepathtype': 'default', 'spec_interp': 'cubic'}, 'datapool': ['whole'], 'collapseax_b': [1, 2, 3], 'collapseax_a': [1, 3], 'sampling': 'resampled', 'cohax': None, 'units': 'K', 'infile_pfx_b': 'xCPDPS_collapse_axes_123', 'modelinfo': {'outfile_pfx_b': ['xCPDPS_collapse_axes_123', 'xCPDPS_collapse_axes_123'], 'PS_dir': ['/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/PS/', '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/PS/'], 'mdl_day': [[[0, 0], [0, 0]]], 'mdl_cohax': [None, None], 'mdl_ndaybins': [2, 2], 'infile_pfx_a': ['xCPDPS_collapse_axes_13', 'xCPDPS_collapse_axes_13'], 'outfile_pfx_a': ['xCPDPS_collapse_axes_13', 'xCPDPS_collapse_axes_13'], 'infile_pfx_b': ['xCPDPS_collapse_axes_123', 'xCPDPS_collapse_axes_123'], 'mdl_collapax_a': [[1, 3], [1, 3]], 'mdl_collapax_b': [[1, 2, 3], [1, 2, 3]], 'mdl_incohax': [[1, 2, 3], [1, 2, 3]]}, 'statistic': 'median', 'action': True, 'infile_pfx_a': 'xCPDPS_collapse_axes_13', 'subselection': {'triaddiag': [0], 'spw': [0, 1], 'day': [[0, 1]], 'lstdiag': [0, 1]}, 'outfile_pfx_a': 'xCPDPS_collapse_axes_13', 'errtype': ['ssdiff'], 'outfile_pfx_b': 'xCPDPS_collapse_axes_123', 'nsigma': 1.0}, '1e': {'action': False, 'triplet': [25, 26, 39], 'infile': 'EQ14', 'obsind': 0}, '1d': {'selection': {'lstrange': [2.9, 3.9], 'triads': None, 'dayind': 0}, 'applyflags': False, 'action': False, 'datastage': 'prelim', 'statistic': 'median', 'sparseness': 10.0}, '2b': {'action': True}, '2c': {'kprll_min': 0.85, 'action': True, 'diagoffsets_b': [[[0], [1], [0]], [[0], [1], [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]], [[0], [1], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]], [[0], [-1, 0, 1], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]]], 'diagoffsets_a': [[[0], [0]], [[0], [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]], [[0], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]], [[0], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]]], 'incohax_b': [[1, 2, 3], [1, 2, 3], [1, 2, 3], [1, 2, 3]], 'incohax_a': [[1, 3], [1, 3], [1, 3], [1, 3]]}, '1h': {'action': False, 'selection': {'lstrange': [3.0, 3.7], 'triads': None, 'dayind': 0}, 'freq_bins': [115000000.0, 120000000.0, 125000000.0, 130000000.0, 135000000.0], 'statistic': 'rms'}, '3a': {'action': True, 'bl': [[14.6, 0.0, 0.0], [29.2, 0.0, 0.0]], 'bltol': 0.2, 'lst': None, 'spw': None}, '1c': {'action': False}, '2a': {'action': False}}, 'preProcessing': {'daybinsize': None, 'flagchans': [], 'band_center': 150000000.0, 'freq_resolution': 97656.25, 'lstbinsize': 60.0, 'flagants': [], 'appendaxis': 'obsid', 'ndaybins': 2, 'action': True, 'append': False}, 'delaySpectrum': {'bl': [[29.2, 0.0, 0.0], [14.6, 25.287942, 0.0], [-14.6, 25.287942, 0.0]], 'pad': 1.0, 'bltol': 0.1, 'subband': {'shape': 'bhw', 'fftpow': 2.0, 'freq_center': [126500000.0, 163000000.0], 'bw_eff': [10000000.0, 10000000.0]}, 'applyflags': False}, 'dirStruct': {'visfiletype': 'prisim', 'infiles': ['EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noisy.npz'], 'projectdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/', 'hdf5_infile': 'EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noisy.hdf5', 'model_hdf5files': ['EQ28XX_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_noiseless.hdf5', 'EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noiseless.hdf5'], 'datadir': 'bispectrum_phase/', 'model_labels': ['FG', 'FG+HI'], 'visfile': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/Fornax_GLEAM_all_FOV30_ephemeris_HERA61/Fornax_transit_NFbeam_120x10.7sec_1024x97.7kHz/simdata/simvis-copy', 'modeldir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/', 'model_npzfiles': ['EQ28XX_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_noiseless.npz', 'EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noiseless.npz'], 'figdir': 'figures/'}, 'telescope': {'latitude': -30.7224, 'longitude': 21.4278}}

In [7]:

    

projectdir = parms['dirStruct']['projectdir']
datadir = projectdir + parms['dirStruct']['datadir']
figdir = datadir + parms['dirStruct']['figdir']
modelsdir = parms['dirStruct']['modeldir']
infiles = parms['dirStruct']['infiles']
visfile = parms['dirStruct']['visfile']
visfiletype = parms['dirStruct']['visfiletype']
hdf5_infile = parms['dirStruct']['hdf5_infile']
model_hdf5files = parms['dirStruct']['model_hdf5files']
model_labels = parms['dirStruct']['model_labels']

telescope_parms = parms['telescope']
site_latitude = telescope_parms['latitude']
site_longitude = telescope_parms['longitude']

preprocessinfo = parms['preProcessing']
preprocess = preprocessinfo['action']
flagchans = preprocessinfo['flagchans']
if flagchans is not None:
    flagchans = NP.asarray(preprocessinfo['flagchans']).reshape(-1)
flagants = preprocessinfo['flagants']
if flagants is not None:
    flagants = NP.asarray(preprocessinfo['flagants']).reshape(-1)
daybinsize = preprocessinfo['daybinsize']
ndaybins = preprocessinfo['ndaybins']
lstbinsize = preprocessinfo['lstbinsize']
band_center = preprocessinfo['band_center']
freq_resolution = preprocessinfo['freq_resolution']

dspecinfo = parms['delaySpectrum']
subbandinfo = dspecinfo['subband']
freq_window_centers = NP.asarray(subbandinfo['freq_center'])
freq_window_bw = NP.asarray(subbandinfo['bw_eff'])
freq_window_shape = subbandinfo['shape']
freq_window_fftpow = subbandinfo['fftpow']
pad = dspecinfo['pad']
apply_flags = dspecinfo['applyflags']
if apply_flags:
    applyflags_str = 'Y'
else:
    applyflags_str = 'N'


    

In [8]:

    

bl = NP.asarray(dspecinfo['bl'])
if bl.shape[0] != 3:
    raise ValueError('Input bl must be made of three vectors forming the triad')
bltol = dspecinfo['bltol']

infile = infiles[0]
infile_no_ext = hdf5_infile.split('.hdf5')[0]

# visdata = NP.load(visfile)
if visfile is None:
    visinfo = None
else:
    if visfiletype == 'hdf5':
        visinfo = NMO.load_dict_from_hdf5(visfile+'.hdf5')
        blind, blrefind, dbl = LKP.find_1NN(visinfo['baseline']['blvect'], bl, distance_ULIM=bltol, remove_oob=True)
        if blrefind.size != 3:
            blind_missing = NP.setdiff1d(NP.arange(3), blind, assume_unique=True)
            blind_next, blrefind_next, dbl_next = LKP.find_1NN(visinfo['baseline']['blvect'], -1*bl[blind_missing,:], distance_ULIM=bltol, remove_oob=True)
            if blind_next.size + blind.size != 3:
                raise ValueError('Exactly three baselines were not found in the reference baselines')
            else:
                blind = NP.append(blind, blind_missing[blind_next])
                blrefind = NP.append(blrefind, blrefind_next)
        else:
            blind_missing = []

        vistriad = MA.array(visinfo['vis_real'][blrefind,:,:] + 1j * visinfo['vis_imag'][blrefind,:,:], mask=visinfo['mask'][blrefind,:,:])
        if len(blind_missing) > 0:
            vistriad[-blrefind_next.size:,:,:] = vistriad[-blrefind_next.size:,:,:].conj()
    else:
        visinfo = RI.InterferometerArray(None, None, None, init_file=visfile)


    

In [9]:

    

print(visfile)
print(visinfo.skyvis_freq.shape)
print(visinfo.lst/15.0) # LST in hours


    

    




/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/Fornax_GLEAM_all_FOV30_ephemeris_HERA61/Fornax_transit_NFbeam_120x10.7sec_1024x97.7kHz/simdata/simvis-copy
(1738, 1024, 120)
[3.21168913 3.21466949 3.21764985 3.22063021 3.22361057 3.22659093
 3.22957129 3.23255165 3.23553201 3.23851237 3.24149273 3.24447309
 3.24745345 3.25043381 3.25341417 3.25639453 3.25937489 3.26235525
 3.26533561 3.26831597 3.27129633 3.27427669 3.27725705 3.28023741
 3.28321777 3.28619813 3.28917849 3.29215885 3.29513921 3.29811957
 3.30109993 3.30408029 3.30706065 3.31004101 3.31302137 3.31600173
 3.31898209 3.32196245 3.32494281 3.32792317 3.33090353 3.33388389
 3.33686425 3.33984461 3.34282497 3.34580533 3.34878569 3.35176605
 3.35474641 3.35772677 3.36070713 3.36368749 3.36666785 3.36964821
 3.37262857 3.37560893 3.37858929 3.38156965 3.38455001 3.38753037
 3.39051073 3.39349109 3.39647145 3.39945181 3.40243217 3.40541253
 3.40839289 3.41137325 3.41435361 3.41733397 3.42031433 3.42329469
 3.42627505 3.42925541 3.43223577 3.43521613 3.43819649 3.44117685
 3.44415721 3.44713757 3.45011793 3.45309829 3.45607865 3.45905901
 3.46203937 3.46501973 3.46800009 3.47098045 3.47396081 3.47694117
 3.47992153 3.48290189 3.48588225 3.48886261 3.49184297 3.49482333
 3.49780369 3.50078405 3.50376441 3.50674477 3.50972513 3.51270549
 3.51568585 3.51866621 3.52164657 3.52462693 3.52760729 3.53058765
 3.53356801 3.53654837 3.53952873 3.54250909 3.54548945 3.54846981
 3.55145017 3.55443053 3.55741089 3.56039125 3.56337161 3.56635197]

In [10]:

    

tmpnpzdata = NP.load(datadir+infile)
nchan = tmpnpzdata['flags'].shape[-1]
freqs = band_center + freq_resolution * (NP.arange(nchan) - int(0.5*nchan))


    

In [11]:

    

# eq28yy_npzfile = '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/Binned_Data/EQ28/YY/EQ28YY.hdf5'
# eq28yy_cpObj = BSP.ClosurePhase(eq28yy_npzfile, freqs, infmt='hdf5')
# print(eq28yy_cpObj.cpinfo['raw']['lst'])


    

In [12]:

    

cpObj = BSP.ClosurePhase(datadir+hdf5_infile, freqs, infmt='hdf5')


    

In [13]:

    

print(cpObj.cpinfo['raw']['lst'])


    

    




[[3.21168913 3.21168913 3.21168913 3.21168913]
 [3.21466949 3.21466949 3.21466949 3.21466949]
 [3.21764985 3.21764985 3.21764985 3.21764985]
 [3.22063021 3.22063021 3.22063021 3.22063021]
 [3.22361057 3.22361057 3.22361057 3.22361057]
 [3.22659093 3.22659093 3.22659093 3.22659093]
 [3.22957129 3.22957129 3.22957129 3.22957129]
 [3.23255165 3.23255165 3.23255165 3.23255165]
 [3.23553201 3.23553201 3.23553201 3.23553201]
 [3.23851237 3.23851237 3.23851237 3.23851237]
 [3.24149273 3.24149273 3.24149273 3.24149273]
 [3.24447309 3.24447309 3.24447309 3.24447309]
 [3.24745345 3.24745345 3.24745345 3.24745345]
 [3.25043381 3.25043381 3.25043381 3.25043381]
 [3.25341417 3.25341417 3.25341417 3.25341417]
 [3.25639453 3.25639453 3.25639453 3.25639453]
 [3.25937489 3.25937489 3.25937489 3.25937489]
 [3.26235525 3.26235525 3.26235525 3.26235525]
 [3.26533561 3.26533561 3.26533561 3.26533561]
 [3.26831597 3.26831597 3.26831597 3.26831597]
 [3.27129633 3.27129633 3.27129633 3.27129633]
 [3.27427669 3.27427669 3.27427669 3.27427669]
 [3.27725705 3.27725705 3.27725705 3.27725705]
 [3.28023741 3.28023741 3.28023741 3.28023741]
 [3.28321777 3.28321777 3.28321777 3.28321777]
 [3.28619813 3.28619813 3.28619813 3.28619813]
 [3.28917849 3.28917849 3.28917849 3.28917849]
 [3.29215885 3.29215885 3.29215885 3.29215885]
 [3.29513921 3.29513921 3.29513921 3.29513921]
 [3.29811957 3.29811957 3.29811957 3.29811957]
 [3.30109993 3.30109993 3.30109993 3.30109993]
 [3.30408029 3.30408029 3.30408029 3.30408029]
 [3.30706065 3.30706065 3.30706065 3.30706065]
 [3.31004101 3.31004101 3.31004101 3.31004101]
 [3.31302137 3.31302137 3.31302137 3.31302137]
 [3.31600173 3.31600173 3.31600173 3.31600173]
 [3.31898209 3.31898209 3.31898209 3.31898209]
 [3.32196245 3.32196245 3.32196245 3.32196245]
 [3.32494281 3.32494281 3.32494281 3.32494281]
 [3.32792317 3.32792317 3.32792317 3.32792317]
 [3.33090353 3.33090353 3.33090353 3.33090353]
 [3.33388389 3.33388389 3.33388389 3.33388389]
 [3.33686425 3.33686425 3.33686425 3.33686425]
 [3.33984461 3.33984461 3.33984461 3.33984461]
 [3.34282497 3.34282497 3.34282497 3.34282497]
 [3.34580533 3.34580533 3.34580533 3.34580533]
 [3.34878569 3.34878569 3.34878569 3.34878569]
 [3.35176605 3.35176605 3.35176605 3.35176605]
 [3.35474641 3.35474641 3.35474641 3.35474641]
 [3.35772677 3.35772677 3.35772677 3.35772677]
 [3.36070713 3.36070713 3.36070713 3.36070713]
 [3.36368749 3.36368749 3.36368749 3.36368749]
 [3.36666785 3.36666785 3.36666785 3.36666785]
 [3.36964821 3.36964821 3.36964821 3.36964821]
 [3.37262857 3.37262857 3.37262857 3.37262857]
 [3.37560893 3.37560893 3.37560893 3.37560893]
 [3.37858929 3.37858929 3.37858929 3.37858929]
 [3.38156965 3.38156965 3.38156965 3.38156965]
 [3.38455001 3.38455001 3.38455001 3.38455001]
 [3.38753037 3.38753037 3.38753037 3.38753037]
 [3.39051073 3.39051073 3.39051073 3.39051073]
 [3.39349109 3.39349109 3.39349109 3.39349109]
 [3.39647145 3.39647145 3.39647145 3.39647145]
 [3.39945181 3.39945181 3.39945181 3.39945181]
 [3.40243217 3.40243217 3.40243217 3.40243217]
 [3.40541253 3.40541253 3.40541253 3.40541253]
 [3.40839289 3.40839289 3.40839289 3.40839289]
 [3.41137325 3.41137325 3.41137325 3.41137325]
 [3.41435361 3.41435361 3.41435361 3.41435361]
 [3.41733397 3.41733397 3.41733397 3.41733397]
 [3.42031433 3.42031433 3.42031433 3.42031433]
 [3.42329469 3.42329469 3.42329469 3.42329469]
 [3.42627505 3.42627505 3.42627505 3.42627505]
 [3.42925541 3.42925541 3.42925541 3.42925541]
 [3.43223577 3.43223577 3.43223577 3.43223577]
 [3.43521613 3.43521613 3.43521613 3.43521613]
 [3.43819649 3.43819649 3.43819649 3.43819649]
 [3.44117685 3.44117685 3.44117685 3.44117685]
 [3.44415721 3.44415721 3.44415721 3.44415721]
 [3.44713757 3.44713757 3.44713757 3.44713757]
 [3.45011793 3.45011793 3.45011793 3.45011793]
 [3.45309829 3.45309829 3.45309829 3.45309829]
 [3.45607865 3.45607865 3.45607865 3.45607865]
 [3.45905901 3.45905901 3.45905901 3.45905901]
 [3.46203937 3.46203937 3.46203937 3.46203937]
 [3.46501973 3.46501973 3.46501973 3.46501973]
 [3.46800009 3.46800009 3.46800009 3.46800009]
 [3.47098045 3.47098045 3.47098045 3.47098045]
 [3.47396081 3.47396081 3.47396081 3.47396081]
 [3.47694117 3.47694117 3.47694117 3.47694117]
 [3.47992153 3.47992153 3.47992153 3.47992153]
 [3.48290189 3.48290189 3.48290189 3.48290189]
 [3.48588225 3.48588225 3.48588225 3.48588225]
 [3.48886261 3.48886261 3.48886261 3.48886261]
 [3.49184297 3.49184297 3.49184297 3.49184297]
 [3.49482333 3.49482333 3.49482333 3.49482333]
 [3.49780369 3.49780369 3.49780369 3.49780369]
 [3.50078405 3.50078405 3.50078405 3.50078405]
 [3.50376441 3.50376441 3.50376441 3.50376441]
 [3.50674477 3.50674477 3.50674477 3.50674477]
 [3.50972513 3.50972513 3.50972513 3.50972513]
 [3.51270549 3.51270549 3.51270549 3.51270549]
 [3.51568585 3.51568585 3.51568585 3.51568585]
 [3.51866621 3.51866621 3.51866621 3.51866621]
 [3.52164657 3.52164657 3.52164657 3.52164657]
 [3.52462693 3.52462693 3.52462693 3.52462693]
 [3.52760729 3.52760729 3.52760729 3.52760729]
 [3.53058765 3.53058765 3.53058765 3.53058765]
 [3.53356801 3.53356801 3.53356801 3.53356801]
 [3.53654837 3.53654837 3.53654837 3.53654837]
 [3.53952873 3.53952873 3.53952873 3.53952873]
 [3.54250909 3.54250909 3.54250909 3.54250909]
 [3.54548945 3.54548945 3.54548945 3.54548945]
 [3.54846981 3.54846981 3.54846981 3.54846981]
 [3.55145017 3.55145017 3.55145017 3.55145017]
 [3.55443053 3.55443053 3.55443053 3.55443053]
 [3.55741089 3.55741089 3.55741089 3.55741089]
 [3.56039125 3.56039125 3.56039125 3.56039125]
 [3.56337161 3.56337161 3.56337161 3.56337161]
 [3.56635197 3.56635197 3.56635197 3.56635197]]

In [14]:

    

print(cpObj.cpinfo.keys())


    

    




[u'raw', 'processed', 'errinfo']

In [15]:

    

print(cpObj.cpinfo['raw'].keys())
print(cpObj.cpinfo['raw']['triads'].shape)
print(cpObj.cpinfo['raw']['days'].shape)
print(cpObj.cpinfo['raw']['lst-day'].shape)
print(cpObj.cpinfo['raw']['lst'].shape)
print(cpObj.cpinfo['raw']['cphase'].shape)


    

    




[u'triads', u'cphase', u'days', u'lst-day', u'lst', u'flags']
(37, 3)
(4,)
(120, 4)
(120, 4)
(120, 4, 37, 1024)

In [16]:

    

print('day bin size = '+str(daybinsize), 'ndaybins={0:0d}'.format(ndaybins), 'LST bin size = {0:.1f}s'.format(lstbinsize))
cpObj.smooth_in_tbins(daybinsize=daybinsize, ndaybins=ndaybins, lstbinsize=lstbinsize)


    

    




('day bin size = None', 'ndaybins=2', 'LST bin size = 60.0s')






    




/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 [17]:

    

print(cpObj.cpinfo['processed'].keys())


    

    




['prelim', 'native']

In [18]:

    

print(cpObj.cpinfo['processed']['prelim'].keys())
print(cpObj.cpinfo['processed']['prelim']['cphase'].keys())
print(cpObj.cpinfo['processed']['prelim']['lstbins'].shape)
print(cpObj.cpinfo['processed']['prelim']['lstbins'])
print(cpObj.cpinfo['processed']['prelim']['daybins'].shape)
print(cpObj.cpinfo['processed']['prelim']['daybins'])
print(cpObj.cpinfo['processed']['prelim']['cphase']['median'].shape)
print(cpObj.cpinfo['processed']['prelim']['eicp']['median'].shape)


    

    




['cphase', 'wts', 'eicp', 'dlstbins', 'daybins', 'lstbins', 'diff_dbins']
['rms', 'median', 'mad', 'mean']
(22,)
[3.22002247 3.23668913 3.2533558  3.27002247 3.28668913 3.3033558
 3.32002247 3.33668913 3.3533558  3.37002247 3.38668913 3.4033558
 3.42002247 3.43668913 3.4533558  3.47002247 3.48668913 3.5033558
 3.52002247 3.53668913 3.5533558  3.57002247]
(2,)
[2458445.90257652 2458447.90257652]
(22, 2, 37, 1024)
(22, 2, 37, 1024)

In [19]:

    

cpObj.subtract(NP.zeros(1024))
print(cpObj.cpinfo['processed'].keys())


    

    




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

In [20]:

    

# ndaybins=4
print('ndaybins={0:0d}'.format(4), 'LST bin size = {0:.1f}s'.format(lstbinsize))
cpObj.subsample_differencing(daybinsize=None, ndaybins=4, lstbinsize=lstbinsize)


    

    




('ndaybins=4', 'LST bin size = 60.0s')

In [21]:

    

print(cpObj.cpinfo['errinfo'].keys())
print(cpObj.cpinfo['errinfo']['daybins'].shape)
print(cpObj.cpinfo['errinfo']['lstbins'].shape)
print(len(cpObj.cpinfo['errinfo']['list_of_pair_of_pairs']))
print(cpObj.cpinfo['errinfo']['list_of_pair_of_pairs'])
print(cpObj.cpinfo['errinfo']['eicp_diff'].keys())
print(cpObj.cpinfo['errinfo']['eicp_diff']['0'].keys())
print(cpObj.cpinfo['errinfo']['eicp_diff']['0']['median'].shape)
print(cpObj.cpinfo['errinfo']['eicp_diff']['1']['median'].shape)


    

    




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

In [22]:

    

cpDSobj = BSP.ClosurePhaseDelaySpectrum(cpObj)


    

In [23]:

    

if visinfo is not None:
    if visfiletype == 'hdf5':
        visscaleinfo = {'vis': vistriad, 'lst': visinfo['header']['LST'], 'smoothinfo': {'op_type': 'interp1d', 'interp_kind': 'linear'}}
    else:
        visscaleinfo = {'vis': visinfo, 'bltriplet': bl, 'smoothinfo': {'op_type': 'interp1d', 'interp_kind': 'linear'}}
else:
    visscaleinfo = None


    

In [24]:

    

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.265e+08, 1.630e+08]))
('freq window BW eff = ', array([10000000., 10000000.]))
freq window shape = bhw
freq window fftpow = 2.0
pad = 1.0

In [25]:

    

cpds = cpDSobj.FT(freq_window_bw, freq_center=freq_window_centers, shape=freq_window_shape, fftpow=freq_window_fftpow, pad=pad, datapool='prelim', visscaleinfo=visscaleinfo, method='fft', resample=True, apply_flags=apply_flags)


    

    




/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)
/lustre/aoc/users/nthyagar/src/miniconda2/envs/PRISim-dev/lib/python2.7/site-packages/astroutils/mathops.py:1344: RuntimeWarning: invalid value encountered in less
  mask_out = NP.logical_or(wts_interped.data < eps, NP.isnan(wts_interped.data)) # Mask small, negative, and NaN weights






    




	Renormalized the shaping window to unit power.
	Renormalized the shaping window to unit power.

In [26]:

    

plot_info = parms['plot']
plots = [key for key in plot_info if plot_info[key]['action']]
print(plots)


    

    




['2d', '2c_err', '3', '2', '2b', '2c', '3a']

In [25]:

    

if '1h' in plots:
    statistic = plot_info['1h']['statistic']
    timetriad_selection = plot_info['1h']['selection']
    cpdev = cpObj.cpinfo['processed']['prelim']['cphase'][statistic]
    if timetriad_selection is not None:
        dayind = timetriad_selection['dayind']
    else:
        dayind = 0
    for key in timetriad_selection:
        if timetriad_selection[key] is not None:
            if key == 'triads':
                triads = map(tuple, timetriad_selection[key])
            elif key == 'lstrange':
                lstrange = timetriad_selection[key]
                lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
                if lstrange is None:
                    lstinds = NP.arange(lstbins.size)
                else:
                    lstrange = NP.asarray(lstrange)
                    lstinds = NP.where(NP.logical_and(lstbins >= lstrange.min(), lstbins <= lstrange.max()))[0]
        else:
            if key == 'triads':
                triads = map(tuple, cpDSobj.cPhase.cpinfo['raw']['triads'])
            elif key == 'lstrange':
                lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
                lstinds = NP.arange(lstbins.size)
    print(cpdev.shape)
    print(lstrange)
    print(lstbins)
    print(lstinds)


    

In [26]:

    

ncol = 3
    nrow = min(4, int(NP.ceil(1.0*lstinds.size/ncol)))
    npages = int(NP.ceil(1.0 * lstinds.size / (nrow*ncol)))
    nlst_remain = lstinds.size
    for pagei in range(npages):
        if pagei > 0:
            nlst_remain = lstinds.size - pagei * nrow * ncol
            nrow = min(4, int(NP.ceil(1.0*nlst_remain/ncol)))
        fig, axs = PLT.subplots(nrows=nrow, ncols=ncol, sharex=True, sharey=True, figsize=(10,8))
        for i in range(nrow):
            for j in range(ncol):
                lstind = (lstinds.size - nlst_remain) + i*ncol+j
                if lstind < lstinds.size:
                    lind = lstinds[lstind]
                    for triad in triads:
                        triad_ind = triads.index(triad)
                        axs[i,j].plot(1e-6*freqs, NP.degrees(cpdev[lind,dayind,triad_ind,:]), marker='.', ms=2, ls='none')
                        axs[i,j].text(0.5, 0.97, '{0:.2f} hrs'.format(lstbins[lind]), transform=axs[i,j].transAxes, fontsize=10, weight='medium', ha='center', va='top', color='black')
                else:
                    axs[i,j].axis('off')
                axs[i,j].set_xlim(1e-6*freqs.min(), 1e-6*freqs.max())
                axs[i,j].set_ylim(0,10)

        fig.subplots_adjust(hspace=0, wspace=0)
        fig.subplots_adjust(top=0.95)
        fig.subplots_adjust(left=0.1)
        fig.subplots_adjust(bottom=0.15)
        fig.subplots_adjust(right=0.98)

        big_ax = fig.add_subplot(111)
        big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
        # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
        big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
        big_ax.set_xticks([])
        big_ax.set_yticks([])
        big_ax.set_xlabel(r'$f$ [MHz]', fontsize=12, weight='medium', labelpad=20)
        big_ax.set_ylabel(r'$\sigma_{\phi_\nabla}$ [degrees]', fontsize=12, weight='medium', labelpad=35)

        # PLT.savefig(figdir + '{0}_cp{1}_spectra_{2}_{2}_{4}_triads_day_{4}_{5:.1f}x_sparse_page_{6:03d}_of_{8:0d}.png'.format(statistic, infile_no_ext, flags_str, datastr, len(triads), dayind, sparseness, pagei+1, npages), bbox_inches=0)
        # PLT.savefig(figdir + '{0}_cp{1}_spectra_{2}_{2}_{4}_triads_day_{4}_{5:.1f}x_sparse_page_{6:03d}_of_{8:0d}.eps'.format(statistic, infile_no_ext, flags_str, datastr, len(triads), dayind, sparseness, pagei+1, npages), bbox_inches=0)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-26-3f2e63de8fb2> in <module>()
      1 ncol = 3
----> 2 nrow = min(4, int(NP.ceil(1.0*lstinds.size/ncol)))
      3 npages = int(NP.ceil(1.0 * lstinds.size / (nrow*ncol)))
      4 nlst_remain = lstinds.size
      5 for pagei in range(npages):

NameError: name 'lstinds' is not defined

In [27]:

    

if ('2' in plots) or ('2a' in plots) or ('2b' in plots) or ('2c' in plots) or ('2c_err' in plots) or ('2d' in plots):
    dir_PS = plot_info['2']['PS_dir']
    infile_pfx_a = plot_info['2']['infile_pfx_a']
    outfile_pfx_a = plot_info['2']['outfile_pfx_a']
    infile_pfx_b = plot_info['2']['infile_pfx_b']
    outfile_pfx_b = plot_info['2']['outfile_pfx_b']
    sampling = plot_info['2']['sampling']
    statistic = plot_info['2']['statistic']
    cohax = plot_info['2']['cohax']
    incohax = plot_info['2']['incohax']
    collapseax_a = plot_info['2']['collapseax_a']
    collapseax_b = plot_info['2']['collapseax_b']
    datapool = plot_info['2']['datapool']
    pspec_unit_type = plot_info['2']['units']
    ps_errtype = plot_info['2']['errtype']
    nsigma = plot_info['2']['nsigma']
    beaminfo = plot_info['2']['beaminfo']
    xlim = plot_info['2']['xlim']
    if infile_pfx_a is not None:
        ps_infile_a = datadir + dir_PS + infile_pfx_a + '_' + infile_no_ext + '.hdf5'
        pserr_infile_a = datadir + dir_PS + infile_pfx_a + '_' + infile_no_ext + '_errinfo.hdf5'
    if outfile_pfx_a is not None:
        ps_outfile_a = datadir + dir_PS + outfile_pfx_a + '_' + infile_no_ext + '.hdf5'
        pserr_outfile_a = datadir + dir_PS + outfile_pfx_a + '_' + infile_no_ext + '_errinfo.hdf5'

    if infile_pfx_b is not None:
        ps_infile_b = datadir + dir_PS + infile_pfx_b + '_' + infile_no_ext + '.hdf5'
        pserr_infile_b = datadir + dir_PS + infile_pfx_b + '_' + infile_no_ext + '_errinfo.hdf5'
    if outfile_pfx_b is not None:
        ps_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '.hdf5'
        pserr_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '_errinfo.hdf5'

    timetriad_selection = plot_info['2']['selection']
    if timetriad_selection is not None:
        dayind = timetriad_selection['days']
    for key in timetriad_selection:
        if timetriad_selection[key] is not None:
            if key == 'triads':
                triads = map(tuple, timetriad_selection[key])
            elif key == 'lstrange':
                lstrange = timetriad_selection[key]
                lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
                if lstrange is None:
                    lstinds = NP.arange(lstbins.size)
                else:
                    lstrange = NP.asarray(lstrange)
                    lstinds = NP.where(NP.logical_and(lstbins >= lstrange.min(), lstbins <= lstrange.max()))[0]
                    if lstinds.size == 0:
                        raise ValueError('No data found in the specified LST range.')
        else:
            if key == 'triads':
                triads = map(tuple, cpDSobj.cPhase.cpinfo['raw']['triads'])
            elif key == 'lstrange':
                lstbins = cpObj.cpinfo['processed']['prelim']['lstbins']
                lstinds = NP.arange(lstbins.size)
    selection = {'triads': triads, 'lst': lstinds, 'days': dayind}
    autoinfo = {'axes': cohax}
    xinfo_a = {'axes': incohax, 'avgcov': False, 'collapse_axes': collapseax_a, 'dlst_range': timetriad_selection['dlst_range']}
    xinfo_b = {'axes': incohax, 'avgcov': False, 'collapse_axes': collapseax_b, 'dlst_range': timetriad_selection['dlst_range']}

    if pspec_unit_type == 'K':
        pspec_unit = 'mK2 Mpc3'
    else:
        pspec_unit = 'Jy2 Mpc'

    subselection = plot_info['2']['subselection']
    mdl_ndaybins = plot_info['2']['modelinfo']['mdl_ndaybins']
    mdl_day = plot_info['2']['modelinfo']['mdl_day']
    mdl_cohax = plot_info['2']['modelinfo']['mdl_cohax']
    mdl_incohax = plot_info['2']['modelinfo']['mdl_incohax']
    mdl_collapseax_a = plot_info['2']['modelinfo']['mdl_collapax_a']
    mdl_collapseax_b = plot_info['2']['modelinfo']['mdl_collapax_b']
    mdl_dir_PS = plot_info['2']['modelinfo']['PS_dir']
    mdl_infile_pfx_a = plot_info['2']['modelinfo']['infile_pfx_a']
    mdl_outfile_pfx_a = plot_info['2']['modelinfo']['outfile_pfx_a']
    mdl_infile_pfx_b = plot_info['2']['modelinfo']['infile_pfx_b']
    mdl_outfile_pfx_b = plot_info['2']['modelinfo']['outfile_pfx_b']

    if model_hdf5files is not None:
        mdl_autoinfo = [{'axes': mdl_cohax[i]} for i in range(len(model_hdf5files))]
        mdl_xinfo_a = [{'axes': mdl_incohax[i], 'avgcov': False, 'collapse_axes': mdl_collapseax_a[i], 'dlst_range': timetriad_selection['dlst_range']} for i in range(len(model_hdf5files))]
        mdl_xinfo_b = [{'axes': mdl_incohax[i], 'avgcov': False, 'collapse_axes': mdl_collapseax_b[i], 'dlst_range': timetriad_selection['dlst_range']} for i in range(len(model_hdf5files))]

    if statistic is None:
        statistic = ['mean', 'median']
    else:
        statistic = [statistic]


    

In [28]:

    

if infile_pfx_a is not None:
        xcpdps2_a = BSP.read_CPhase_cross_power_spectrum(ps_infile_a)
        xcpdps2_a_errinfo = BSP.read_CPhase_cross_power_spectrum(pserr_infile_a)
    else:
        xcpdps2_a = cpDSobj.compute_power_spectrum(selection=selection, autoinfo=autoinfo, xinfo=xinfo_a, units=pspec_unit_type, beamparms=beaminfo)
        xcpdps2_a_errinfo = cpDSobj.compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_a, units=pspec_unit_type, beamparms=beaminfo)
    if outfile_pfx_a is not None:
        BSP.save_CPhase_cross_power_spectrum(xcpdps2_a, ps_outfile_a)
        BSP.save_CPhase_cross_power_spectrum(xcpdps2_a_errinfo, pserr_outfile_a)

    if infile_pfx_b is not None:
        xcpdps2_b = BSP.read_CPhase_cross_power_spectrum(ps_infile_b)
        xcpdps2_b_errinfo = BSP.read_CPhase_cross_power_spectrum(pserr_infile_b)
    else:
        xcpdps2_b = cpDSobj.compute_power_spectrum(selection=selection, autoinfo=autoinfo, xinfo=xinfo_b, units=pspec_unit_type, beamparms=beaminfo)
        xcpdps2_b_errinfo = cpDSobj.compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_b, units=pspec_unit_type, beamparms=beaminfo)
    if outfile_pfx_b is not None:
        BSP.save_CPhase_cross_power_spectrum(xcpdps2_b, ps_outfile_b)
        BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_errinfo, pserr_outfile_b)


    

In [29]:

    

nsamples_incoh = xcpdps2_a[sampling]['whole']['nsamples_incoh']
    nsamples_coh = xcpdps2_a[sampling]['whole']['nsamples_coh']


    

In [30]:

    

model_cpObjs = []
    model_cpDSobjs = []
    cpds_models = []
    xcpdps2_a_models = []
    xcpdps2_a_errinfo_models = []
    xcpdps2_b_models = []
    xcpdps2_b_errinfo_models = []

    if model_hdf5files is not None:
        if mdl_infile_pfx_a is not None:
            if isinstance(mdl_infile_pfx_a, list):
                if (len(mdl_infile_pfx_a) > 0):
                    if not isinstance(mdl_dir_PS, list):
                        if isinstance(mdl_dir_PS, str):
                            mdl_dir_PS = [mdl_dir_PS] * len(model_hdf5files)
                        else:
                            raise TypeError('PS directory for models must be a list of strings')
                    else:
                        if len(mdl_dir_PS) != len(model_hdf5files):
                            raise ValueError('Input model PS directories must match the number of models being analyzed.')
            else:
                raise TypeError('Input model PS infile_a prefixes must be specified as a list of strings')

        if mdl_infile_pfx_b is not None:
            if isinstance(mdl_infile_pfx_b, list):
                if len(mdl_infile_pfx_b) != len(mdl_infile_pfx_b):
                    raise ValueError('Length of input model PS infile_b prefixes must match the length of input model PS infile_a prefixes')
            else:
                raise TypeError('Input model PS infile_b prefixes must be specified as a list of strings')

        progress = PGB.ProgressBar(widgets=[PGB.Percentage(), PGB.Bar(marker='-', left=' |', right='| '), PGB.Counter(), '/{0:0d} Models '.format(len(model_hdf5files)), PGB.ETA()], maxval=len(model_hdf5files)).start()

        for i in range(len(model_hdf5files)):
            mdl_infile_no_ext = model_hdf5files[i].split('.hdf5')[0]
            mdl_ps_infile_a_provided = False
            mdl_pserr_infile_a_provided = False
            mdl_ps_infile_b_provided = False
            mdl_pserr_infile_b_provided = False
            if mdl_infile_pfx_a is not None:
                if len(mdl_infile_pfx_a) > 0:
                    if mdl_infile_pfx_a[i] is not None:
                        if not isinstance(mdl_infile_pfx_a[i], str):
                            raise TypeError('Input {0}-th model cross PS file must be a string'.format(i+1))
                        else:
                            try:
                                model_xcpdps2_a = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'_'+mdl_infile_no_ext+'.hdf5')
                            except IOError as xcption:
                                mdl_ps_infile_a_provided = False
                                warnings.warn('Provided model cross-power spectrum infile_a "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'.hdf5'))
                            else:
                                mdl_ps_infile_a_provided = True
                                xcpdps2_a_models += [copy.deepcopy(model_xcpdps2_a)]

                            try:
                                model_xcpdps2_a_errinfo = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'_'+mdl_infile_no_ext+'_errinfo.hdf5')
                            except IOError as xcption:
                                mdl_pserr_infile_a_provided = False
                                warnings.warn('Provided model cross-power spectrum infile_a "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_a[i]+'_errinfo.hdf5'))
                            else:
                                mdl_pserr_infile_a_provided = True
                                xcpdps2_a_errinfo_models += [copy.deepcopy(model_xcpdps2_a_errinfo)]

            if mdl_infile_pfx_b is not None:
                if len(mdl_infile_pfx_b) > 0:
                    if mdl_infile_pfx_b[i] is not None:
                        if not isinstance(mdl_infile_pfx_b[i], str):
                            raise TypeError('Input {0}-th model cross PS file must be a string'.format(i+1))
                        else:
                            try:
                                model_xcpdps2_b = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'_'+mdl_infile_no_ext+'.hdf5')
                            except IOError as xcption:
                                mdl_ps_infile_b_provided = False
                                warnings.warn('Provided model cross-power spectrum infile_b "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'.hdf5'))
                            else:
                                mdl_ps_infile_b_provided = True
                                xcpdps2_b_models += [copy.deepcopy(model_xcpdps2_b)]

                            try:
                                model_xcpdps2_b_errinfo = BSP.read_CPhase_cross_power_spectrum(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'_'+mdl_infile_no_ext+'_errinfo.hdf5')
                            except IOError as xcption:
                                mdl_pserr_infile_b_provided = False
                                warnings.warn('Provided model cross-power spectrum infile_b "{0}" could not be opened. Will proceed with computing of model cross power spectrum based on parameters specified.'.format(mdl_dir_PS[i]+mdl_infile_pfx_b[i]+'_errinfo.hdf5'))
                            else:
                                mdl_pserr_infile_b_provided = True
                                xcpdps2_b_errinfo_models += [copy.deepcopy(model_xcpdps2_b_errinfo)]

            if (not mdl_ps_infile_a_provided) or (not mdl_pserr_infile_a_provided) or (not mdl_ps_infile_b_provided) or (not mdl_pserr_infile_b_provided):
                model_cpObj = BSP.ClosurePhase(modelsdir+model_hdf5files[i], freqs, infmt='hdf5')
                model_cpObj.smooth_in_tbins(daybinsize=daybinsize, ndaybins=mdl_ndaybins[i], lstbinsize=lstbinsize)
                model_cpObj.subsample_differencing(daybinsize=None, ndaybins=4, lstbinsize=lstbinsize)
                model_cpObj.subtract(NP.zeros(1024))
                model_cpObjs += [copy.deepcopy(model_cpObj)]
                model_cpDSobjs += [BSP.ClosurePhaseDelaySpectrum(model_cpObjs[i])]
                cpds_models += [model_cpDSobjs[i].FT(freq_window_bw, freq_center=freq_window_centers, shape=freq_window_shape, fftpow=freq_window_fftpow, pad=pad, datapool='prelim', visscaleinfo=visscaleinfo, method='fft', resample=True, apply_flags=apply_flags)]

                if not mdl_ps_infile_a_provided:
                    xcpdps2_a_models += [model_cpDSobjs[i].compute_power_spectrum(selection=selection, autoinfo=mdl_autoinfo[i], xinfo=mdl_xinfo_a[i], units=pspec_unit_type, beamparms=beaminfo)]
                if not mdl_pserr_infile_a_provided:
                    xcpdps2_a_errinfo_models += [model_cpDSobjs[i].compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_a, units=pspec_unit_type, beamparms=beaminfo)]

                if not mdl_ps_infile_b_provided:
                    xcpdps2_b_models += [model_cpDSobjs[i].compute_power_spectrum(selection=selection, autoinfo=mdl_autoinfo[i], xinfo=mdl_xinfo_b[i], units=pspec_unit_type, beamparms=beaminfo)]
                if not mdl_pserr_infile_b_provided:
                    xcpdps2_b_errinfo_models += [model_cpDSobjs[i].compute_power_spectrum_uncertainty(selection=selection, autoinfo=autoinfo, xinfo=xinfo_b, units=pspec_unit_type, beamparms=beaminfo)]

            else:
                model_cpObjs += [None]
                model_cpDSobjs += [None]
                cpds_models += [None]

            if mdl_outfile_pfx_a is not None:
                if isinstance(mdl_outfile_pfx_a, str):
                    mdl_outfile_pfx_a = [mdl_outfile_pfx_a] * len(model_hdf5files)
                if not isinstance(mdl_outfile_pfx_a, list):
                    raise TypeError('The model cross-power spectrum outfile prefixes must be specified as a list with item for each model.')
                if len(mdl_outfile_pfx_a) != len(mdl_dir_PS):
                    raise ValueError('Invalid number of model cross-power output files specified')
                mdl_ps_outfile_a = mdl_dir_PS[i] + mdl_outfile_pfx_a[i] + '_' + mdl_infile_no_ext + '.hdf5'
                mdl_pserr_outfile_a = mdl_dir_PS[i] + mdl_outfile_pfx_a[i] + '_' + mdl_infile_no_ext + '_errinfo.hdf5'
                BSP.save_CPhase_cross_power_spectrum(xcpdps2_a_models[-1], mdl_ps_outfile_a)
                BSP.save_CPhase_cross_power_spectrum(xcpdps2_a_errinfo_models[-1], mdl_pserr_outfile_a)

            if mdl_outfile_pfx_b is not None:
                if isinstance(mdl_outfile_pfx_b, str):
                    mdl_outfile_pfx_b = [mdl_outfile_pfx_b] * len(model_hdf5files)
                if not isinstance(mdl_outfile_pfx_b, list):
                    raise TypeError('The model cross-power spectrum outfile prefixes must be specified as a list with item for each model.')
                if len(mdl_outfile_pfx_b) != len(mdl_dir_PS):
                    raise ValueError('Invalid number of model cross-power output files specified')
                mdl_ps_outfile_b = mdl_dir_PS[i] + mdl_outfile_pfx_b[i] + '_' + mdl_infile_no_ext + '.hdf5'
                mdl_pserr_outfile_b = mdl_dir_PS[i] + mdl_outfile_pfx_b[i] + '_' + mdl_infile_no_ext + '_errinfo.hdf5'
                BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_models[-1], mdl_ps_outfile_b)
                BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_errinfo_models[-1], mdl_pserr_outfile_b)

            progress.update(i+1)
        progress.finish()


    

    




100% |-----------------------------------------------| 2/2 Models Time: 0:00:08

In [31]:

    

spw = subselection['spw']
    if spw is None:
        spwind = NP.arange(xcpdps2_a[sampling]['z'].size)
    else:
        spwind = NP.asarray(spw)
    lstind = NMO.find_list_in_list(xcpdps2_a[sampling][datapool[0]]['diagoffsets'][1], NP.asarray(subselection['lstdiag']))
    dayind = NP.asarray(subselection['day'])
    dayind_models = NP.asarray(mdl_day)
    triadind = NMO.find_list_in_list(xcpdps2_a[sampling][datapool[0]]['diagoffsets'][3], NP.asarray(subselection['triaddiag']))

    errshade = {}
    for errtype in ps_errtype:
        if errtype.lower() == 'ssdiff':
            errshade[errtype] = '0.8'
        elif errtype.lower() == 'psdiff':
            errshade[errtype] = '0.6'
    
    mdl_colrs = ['red', 'green', 'blue', 'cyan', 'gray', 'orange']


    

In [32]:

    

if '2b' in plots:
        for stat in statistic:
            for zind in spwind:
                for lind in lstind:
                    for di,dind in enumerate(dayind):
                        maxabsvals = []
                        minabsvals = []
                        maxvals = []
                        minvals = []
                        fig, axs = PLT.subplots(nrows=1, ncols=len(datapool), sharex=True, sharey=True, figsize=(4.0*len(datapool), 3.6))
                        if len(datapool) == 1:
                            axs = [axs]
                        for dpoolind,dpool in enumerate(datapool):
                            for trno,trind in enumerate([triadind[0]]):
                                if model_hdf5files is not None:
                                    for mdlind, mdl in enumerate(model_labels):
                                        if dpool in xcpdps2_a_models[mdlind][sampling]:
                                            psval = (2/3.0) * xcpdps2_a_models[mdlind][sampling][dpool][stat][zind,lind,dayind_models[di][mdlind][0],dayind_models[di][mdlind][1],trind,:].to(pspec_unit).value

                                            # negind = psval.real < 0.0
                                            # posind = NP.logical_not(negind)
                                            maxabsvals += [NP.abs(psval.real).max()]
                                            minabsvals += [NP.abs(psval.real).min()]
                                            maxvals += [psval.real.max()]
                                            minvals += [psval.real.min()]
                                            axs[dpoolind].plot(xcpdps2_a_models[mdlind][sampling]['kprll'][zind,:], psval.real, ls='none', marker='.', ms=3, color=mdl_colrs[mdlind], label='{0}'.format(mdl))

                                if dpool in xcpdps2_a[sampling]:
                                    psval = (2/3.0) * xcpdps2_a[sampling][dpool][stat][zind,lind,dind[0],dind[1],trind,:].to(pspec_unit).value
                                    psrms = (2/3.0) * NP.nanstd(xcpdps2_a_errinfo[sampling]['errinfo'][stat][zind,lind,:,trind,:], axis=0).to(pspec_unit).value

                                    maxabsvals += [NP.abs(psval.real + psrms).max()]
                                    minabsvals += [NP.abs(psval.real).min()]
                                    maxvals += [(psval.real + psrms).max()]
                                    minvals += [(psval.real - psrms).min()]

                                    # axs[dpoolind].plot(xcpdps2_a[sampling]['kprll'][zind,:], psval.real, ls='none', marker='.', ms=1, color='black', label='FG+N')
                                    axs[dpoolind].errorbar(xcpdps2_a[sampling]['kprll'][zind,:], psval.real, yerr=psrms, xerr=None, ecolor='0.8', ls='none', marker='.', ms=4, color='black', label='FG+N')

                                legend = axs[dpoolind].legend(loc='center', bbox_to_anchor=(0.5,0.3), shadow=False, fontsize=8)
                                if trno == 0:
                                    axs[dpoolind].text(0.05, 0.97, 'Real', transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
                                    axs[dpoolind].text(0.95, 0.97, r'$z=$'+' {0:.1f}'.format(xcpdps2_a[sampling]['z'][zind]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='right', va='top', color='black')
                                    axs[dpoolind].text(0.05, 0.92, r'$\Delta$'+'LST = {0:.1f} s'.format(lind*3.6e3*xcpdps2_a['dlst'][0]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
                                    axs[dpoolind].text(0.05, 0.87, 'G{0[0]:0d}{0[1]:0d}'.format(dind), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')

                                axt = axs[dpoolind].twiny()
                                axt.set_xlim(1e6*xcpdps2_a[sampling]['lags'].min(), 1e6*xcpdps2_a[sampling]['lags'].max())

                            minvals = NP.asarray(minvals)
                            maxvals = NP.asarray(maxvals)
                            minabsvals = NP.asarray(minabsvals)
                            maxabsvals = NP.asarray(maxabsvals)
                            if xlim is None:
                                axs[dpoolind].set_xlim(0.99*xcpdps2_a[sampling]['kprll'][zind,:].min(), 1.01*xcpdps2_a[sampling]['kprll'][zind,:].max())
                            else:
                                axs[dpoolind].set_xlim(xlim)
                            if NP.min(minvals) < 0.0:
                                axs[dpoolind].set_ylim(1.5*NP.min(minvals), 2*NP.max(maxabsvals))
                            else:
                                axs[dpoolind].set_ylim(0.5*NP.min(minvals), 2*NP.max(maxabsvals))
                            axs[dpoolind].set_yscale('symlog', linthreshy=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))))
                            tickloc = PLTick.SymmetricalLogLocator(linthresh=10**NP.floor(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))), base=100.0)
                            axs[dpoolind].yaxis.set_major_locator(tickloc)
                            axs[dpoolind].grid(color='0.8', which='both', linestyle=':', lw=1)

                        fig.subplots_adjust(top=0.85)
                        fig.subplots_adjust(bottom=0.16)
                        fig.subplots_adjust(left=0.22)
                        fig.subplots_adjust(right=0.98)

                        big_ax = fig.add_subplot(111)
                        big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
                        # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
                        big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
                        big_ax.set_xticks([])
                        big_ax.set_yticks([])
                        big_ax.set_xlabel(r'$\kappa_\parallel$'+' [psuedo '+r'$h$'+' Mpc'+r'$^{-1}$'+']', fontsize=12, weight='medium', labelpad=20)
                        if pspec_unit_type == 'K':
                            big_ax.set_ylabel(r'$\frac{1}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo mK$^2h^{-3}$ Mpc$^3$]', fontsize=12, weight='medium', labelpad=40)
                        else:
                            big_ax.set_ylabel(r'$\frac{1}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo Jy$^2h^{-1}$ Mpc]', fontsize=12, weight='medium', labelpad=40)

                        big_axt = big_ax.twiny()
                        big_axt.set_xticks([])
                        big_axt.set_xlabel(r'$\tau$'+' ['+r'$\mu$'+'s]', fontsize=12, weight='medium', labelpad=20)

                        PLT.savefig(figdir + '{0}_symlog_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_lstdiag_{5:0d}_day_{6[0]:0d}_{6[1]:0d}_triaddiags_{7:0d}_flags_{8}.pdf'.format(infile_no_ext, xcpdps2_a[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_a['dlst'][0], subselection['lstdiag'][lind], dind, xcpdps2_a[sampling][datapool[0]]['diagoffsets'][3][trind], applyflags_str), bbox_inches=0)
                        print(figdir + '{0}_symlog_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_lstdiag_{5:0d}_day_{6[0]:0d}_{6[1]:0d}_triaddiags_{7:0d}_flags_{8}.pdf'.format(infile_no_ext, xcpdps2_a[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_a['dlst'][0], subselection['lstdiag'][lind], dind, xcpdps2_a[sampling][datapool[0]]['diagoffsets'][3][trind], applyflags_str))


    

    




/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/figures/EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noisy_symlog_real_cpdps_z_10.2_median_resampled_dlst_60.0s_lstdiag_0_day_0_1_triaddiags_0_flags_N.pdf
/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/figures/EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noisy_symlog_real_cpdps_z_10.2_median_resampled_dlst_60.0s_lstdiag_1_day_0_1_triaddiags_0_flags_N.pdf
/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/figures/EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noisy_symlog_real_cpdps_z_7.7_median_resampled_dlst_60.0s_lstdiag_0_day_0_1_triaddiags_0_flags_N.pdf
/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/figures/EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noisy_symlog_real_cpdps_z_7.7_median_resampled_dlst_60.0s_lstdiag_1_day_0_1_triaddiags_0_flags_N.pdf






    













    













    













    

In [33]:

    

if ('2c' in plots) or ('2c_err' in plots) or ('2d' in plots):
        kprll_min_for_rms = plot_info['2c']['kprll_min']
        if kprll_min_for_rms is None:
            kprll_min_for_rms = 0.0
        avg_incohax_a = plot_info['2c']['incohax_a']
        diagoffsets_incohax_a = plot_info['2c']['diagoffsets_a']
        diagoffsets_a = []
        avg_incohax_b = plot_info['2c']['incohax_b']
        diagoffsets_incohax_b = plot_info['2c']['diagoffsets_b']
        diagoffsets_b = []
        for combi,incax_comb in enumerate(avg_incohax_a):
            diagoffsets_a += [{}]
            for incaxind,incax in enumerate(incax_comb):
                diagoffsets_a[-1][incax] = NP.asarray(diagoffsets_incohax_a[combi][incaxind])
        xcpdps2_a_avg, excpdps2_a_avg = BSP.incoherent_cross_power_spectrum_average(xcpdps2_a, excpdps=xcpdps2_a_errinfo, diagoffsets=diagoffsets_a)

        avg_xcpdps2_a_models = []
        avg_excpdps2_a_models = []

        for combi,incax_comb in enumerate(avg_incohax_b):
            diagoffsets_b += [{}]
            for incaxind,incax in enumerate(incax_comb):
                diagoffsets_b[-1][incax] = NP.asarray(diagoffsets_incohax_b[combi][incaxind])

        # xcpdps2_b_avg, excpdps2_b_avg = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b, excpdps=None, diagoffsets=diagoffsets_b)
        xcpdps2_b_avg, excpdps2_b_avg = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b, excpdps=xcpdps2_b_errinfo, diagoffsets=diagoffsets_b)
        avg_xcpdps2_b_models = []
        avg_excpdps2_b_models = []

        if model_hdf5files is not None:
            progress = PGB.ProgressBar(widgets=[PGB.Percentage(), PGB.Bar(marker='-', left=' |', right='| '), PGB.Counter(), '/{0:0d} Models '.format(len(model_hdf5files)), PGB.ETA()], maxval=len(model_hdf5files)).start()

            for i in range(len(model_hdf5files)):
                avg_xcpdps2_a_model, avg_excpdps2_a_model = BSP.incoherent_cross_power_spectrum_average(xcpdps2_a_models[i], excpdps=xcpdps2_a_errinfo_models[i], diagoffsets=diagoffsets_a)
                avg_xcpdps2_a_models += [copy.deepcopy(avg_xcpdps2_a_model)]
                avg_excpdps2_a_models += [copy.deepcopy(avg_excpdps2_a_model)]

                # avg_xcpdps2_b_model, avg_excpdps2_b_model = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b_models[i], excpdps=None, diagoffsets=diagoffsets_b)
                avg_xcpdps2_b_model, avg_excpdps2_b_model = BSP.incoherent_cross_power_spectrum_average(xcpdps2_b_models[i], excpdps=xcpdps2_b_errinfo_models[i], diagoffsets=diagoffsets_b)
                avg_xcpdps2_b_models += [copy.deepcopy(avg_xcpdps2_b_model)]
                avg_excpdps2_b_models += [copy.deepcopy(avg_excpdps2_b_model)]

                progress.update(i+1)
            progress.finish()


    

    




100% |-----------------------------------------------| 2/2 Models Time: 0:00:07

In [34]:

    

# Save incoherent cross power average of the main dataset and its uncertainties
    
        xps_avg_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '_incoh_diag_avg.hdf5'
        xpserr_avg_outfile_b = datadir + dir_PS + outfile_pfx_b + '_' + infile_no_ext + '_errinfo_incoh_diag_avg.hdf5'        
        BSP.save_CPhase_cross_power_spectrum(xcpdps2_b_avg, xps_avg_outfile_b)
        BSP.save_CPhase_cross_power_spectrum(excpdps2_b_avg, xpserr_avg_outfile_b)


    

In [68]:

    

print(infile_no_ext)


    

    




EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_Fornax_GLEAM_all_FOV_30_ephemeris_HERA61_14.0_arcmin_v2_noisy

In [66]:

    

if '2c' in plots:

#             ylim = {0: [-1e4, 1e10], 1: [-1e4, 1e9]} # keys are spw indices
#             linthreshy = 1e0
            ylim = {0: [-1e10, 1e16], 1: [-1e10, 1e15]} # keys are spw indices
            linthreshy = 9e5

            lstind = [0]
            triadind = [0]
            dayind = [0]
            dayind_models = NP.zeros(len(model_labels), dtype=int).reshape(1,-1)

            for stat in statistic:
                for zind in spwind:
                    kprll_ind_for_rms = NP.where(NP.abs(excpdps2_b_avg[sampling]['kprll'][zind,:]) >= kprll_min_for_rms)[0]
                    kprll_for_rms = excpdps2_b_avg[sampling]['kprll'][zind,kprll_ind_for_rms]
                    for lind in lstind:
                        for di,dind in enumerate(dayind):
                            for combi in range(len(diagoffsets_b)):
                                maxabsvals = []
                                minabsvals = []
                                maxvals = []
                                minvals = []
                                fig, axs = PLT.subplots(nrows=1, ncols=len(datapool), sharex=True, sharey=True, figsize=(4.0*len(datapool), 3.6))
                                if len(datapool) == 1:
                                    axs = [axs]
                                for dpoolind,dpool in enumerate(datapool):
                                    for trno,trind in enumerate(triadind):
#                                         if model_hdf5files is not None:
#                                             for mdlind, mdl in enumerate(model_labels):
#                                                 if dpool in avg_xcpdps2_b_models[mdlind][sampling]:
#                                                     psval = (1/3.0) * avg_xcpdps2_b_models[mdlind][sampling][dpool][stat][combi][zind,lind,dayind_models[di][mdlind],trind,:].to(pspec_unit).value
#                                                     maxabsvals += [NP.abs(psval.real).max()]
#                                                     minabsvals += [NP.abs(psval.real).min()]
#                                                     maxvals += [psval.real.max()]
#                                                     minvals += [psval.real.min()]
#                                                     axs[dpoolind].plot(avg_xcpdps2_b_models[mdlind][sampling]['kprll'][zind,:], psval.real, ls='none', marker='.', ms=3, color=mdl_colrs[mdlind], label='{0}'.format(mdl))

                                        if dpool in xcpdps2_b_avg[sampling]:
                                            psval = (2/3.0) * xcpdps2_b_avg[sampling][dpool][stat][combi][zind,lind,dind,trind,:].to(pspec_unit).value
                                            ps_ssdiff = (2/3.0) * excpdps2_b_avg[sampling]['errinfo'][stat][combi][zind,lind,:,trind,kprll_ind_for_rms].to(pspec_unit).value # Single RMS across all k_prll bins
                                            psrms_ssdiff = (2/3.0) * NP.nanstd(excpdps2_a_avg[sampling]['errinfo'][stat][combi][zind,lind,:,trind,:], axis=0).to(pspec_unit).value # RMS per k_prll bin
                                            if 2 in avg_incohax_b[combi]:
                                                ind_dayax_in_incohax = avg_incohax_b[combi].index(2)
                                                if 0 in diagoffsets_incohax_b[combi][ind_dayax_in_incohax]:
                                                    rms_inflation_factor = 2.0 * NP.sqrt(2.0)
                                                else:
                                                    rms_inflation_factor = NP.sqrt(2.0)
                                            else:
                                                rms_inflation_factor = NP.sqrt(2.0)
                                            psrms_psdiff = (2/3.0) * (xcpdps2_a_avg[sampling][dpool][stat][combi][zind,lind,1,1,trind,:] - xcpdps2_a_avg[sampling][dpool][stat][combi][zind,lind,0,0,trind,:]).to(pspec_unit).value
                                            psrms_psdiff = NP.abs(psrms_psdiff.real) / rms_inflation_factor

                                            psrms_max = NP.amax(NP.vstack((psrms_ssdiff, psrms_psdiff)), axis=0)

                                            maxabsvals += [NP.abs(psval.real + nsigma*psrms_max).max()]
                                            minabsvals += [NP.abs(psval.real).min()]
                                            maxvals += [(psval.real + nsigma*psrms_max).max()]
                                            minvals += [(psval.real - nsigma*psrms_max).min()]

                                            for errtype in ps_errtype:
                                                if errtype.lower() == 'ssdiff':                                                    
                                                    axs[dpoolind].errorbar(xcpdps2_b_avg[sampling]['kprll'][zind,:], psval.real, yerr=nsigma*psrms_ssdiff, xerr=None, ecolor=errshade[errtype.lower()], ls='none', marker='.', ms=4, color='black')
                                                    psmean_ssdiff_across_kbins, psmedian_ssdiff_across_kbins, psrms_ssdiff_across_kbins = sigma_clipped_stats(NP.vstack((ps_ssdiff.real, ps_ssdiff.imag)))
                                                    pos_kprll_ind_for_rms = NP.where(kprll_for_rms >= 0.0)[0]
                                                    neg_kprll_ind_for_rms = NP.where(kprll_for_rms <= 0.0)[0]
                                                    axs[dpoolind].hlines([-psrms_ssdiff_across_kbins, psrms_ssdiff_across_kbins], kprll_for_rms[pos_kprll_ind_for_rms].min(), kprll_for_rms[pos_kprll_ind_for_rms].max(), linestyles='--', lw=1.5, color='0.25')
                                                    axs[dpoolind].hlines([-psrms_ssdiff_across_kbins, psrms_ssdiff_across_kbins], kprll_for_rms[neg_kprll_ind_for_rms].min(), kprll_for_rms[neg_kprll_ind_for_rms].max(), linestyles='--', lw=1.5, color='0.25')
                                                elif errtype.lower() == 'psdiff':
                                                    axs[dpoolind].errorbar(xcpdps2_b_avg[sampling]['kprll'][zind,:], psval.real, yerr=nsigma*psrms_psdiff, xerr=None, ecolor=errshade[errtype.lower()], ls='none', marker='.', ms=4, color='black', label='FG+N')

#                                         legend = axs[dpoolind].legend(loc='center', bbox_to_anchor=(0.5,0.3), shadow=False, fontsize=8)
                                        if trno == 0:
                                            # axs[dpoolind].text(0.05, 0.97, 'Real', transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
                                            axs[dpoolind].text(0.95, 0.97, r'$z=$'+' {0:.1f}'.format(xcpdps2_b_avg[sampling]['z'][zind]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='right', va='top', color='black')
                                            # axs[dpoolind].text(0.05, 0.92, r'$\Delta$'+'LST = {0:.1f} s'.format(lind*3.6e3*xcpdps2_a_avg['dlst'][0]), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')
                                            # axs[dpoolind].text(0.05, 0.87, 'G{0[0]:0d}{0[1]:0d}'.format(dind), transform=axs[dpoolind].transAxes, fontsize=8, weight='medium', ha='left', va='top', color='black')

                                        axt = axs[dpoolind].twiny()
                                        axt.set_xlim(1e6*xcpdps2_b_avg[sampling]['lags'].min(), 1e6*xcpdps2_b_avg[sampling]['lags'].max())

                                    axs[dpoolind].axhline(y=0, xmin=0, xmax=1, ls='-', lw=1, color='black')

                                    minvals = NP.asarray(minvals)
                                    maxvals = NP.asarray(maxvals)
                                    minabsvals = NP.asarray(minabsvals)
                                    maxabsvals = NP.asarray(maxabsvals)
                                    if xlim is None:
                                        axs[dpoolind].set_xlim(0.99*xcpdps2_b_avg[sampling]['kprll'][zind,:].min(), 1.01*xcpdps2_b_avg[sampling]['kprll'][zind,:].max())
                                    else:
                                        axs[dpoolind].set_xlim(xlim)
#                                     if NP.min(minvals) < 0.0:
#                                         axs[dpoolind].set_ylim(1.5*NP.min(minvals), 2*NP.max(maxabsvals))
#                                     else:
#                                         axs[dpoolind].set_ylim(0.5*NP.min(minvals), 2*NP.max(maxabsvals))
                                    axs[dpoolind].set_ylim(ylim[zind])
                                        
#                                     axs[dpoolind].set_yscale('symlog', linthreshy=10**NP.ceil(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))))
#                                     tickloc = PLTick.SymmetricalLogLocator(linthresh=10**NP.ceil(NP.log10(NP.min(minabsvals[minabsvals > 0.0]))), base=100.0)
                                    axs[dpoolind].set_yscale('symlog', linthreshy=linthreshy)
                                    tickloc = PLTick.SymmetricalLogLocator(linthresh=linthreshy, base=100.0)

                                    axs[dpoolind].yaxis.set_major_locator(tickloc)
                                    yticklocs = NP.asarray(axs[dpoolind].get_yticks())
                                    ytickspacings = NP.diff(yticklocs)
                                    tickinds_to_remove = NP.setdiff1d(NP.where(NP.abs(yticklocs)<=linthreshy)[0], NP.argmin(NP.abs(yticklocs)))
                                    if tickinds_to_remove.size > 0:
                                        new_yticklocs = yticklocs[NP.setdiff1d(NP.arange(yticklocs.size), tickinds_to_remove)]
                                        axs[dpoolind].yaxis.set_major_locator(PLTick.FixedLocator((new_yticklocs.tolist())))
                                  
                                    axs[dpoolind].grid(color='0.8', which='both', linestyle=':', lw=1)

                                fig.subplots_adjust(top=0.85)
                                fig.subplots_adjust(bottom=0.16)
                                fig.subplots_adjust(left=0.22)
                                fig.subplots_adjust(right=0.98)

                                big_ax = fig.add_subplot(111)
                                big_ax.set_facecolor('none') # matplotlib.__version__ >= 2.0.0
                                # big_ax.set_axis_bgcolor('none') # matplotlib.__version__ < 2.0.0
                                big_ax.tick_params(labelcolor='none', top=False, bottom=False, left=False, right=False)
                                big_ax.set_xticks([])
                                big_ax.set_yticks([])
                                big_ax.set_xlabel(r'$\kappa_\parallel$'+' [pseudo '+r'$h$'+' Mpc'+r'$^{-1}$'+']', fontsize=12, weight='medium', labelpad=20)
                                if pspec_unit_type == 'K':
                                    big_ax.set_ylabel(r'$\frac{2}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo mK$^2h^{-3}$ Mpc$^3$]', fontsize=12, weight='medium', labelpad=35)
                                else:
                                    big_ax.set_ylabel(r'$\frac{2}{3}\, P_\nabla(\kappa_\parallel)$ [pseudo Jy$^2h^{-1}$ Mpc]', fontsize=12, weight='medium', labelpad=35)

                                big_axt = big_ax.twiny()
                                big_axt.set_xticks([])
                                big_axt.set_xlabel(r'$\tau$'+' ['+r'$\mu$'+'s]', fontsize=12, weight='medium', labelpad=20)

                                PLT.savefig(figdir + '{0}_symlog_incoh_avg_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_b_avg[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg['dlst'][0], applyflags_str, combi), bbox_inches=0)
                                print(figdir + '{0}_symlog_incoh_avg_real_cpdps_z_{1:.1f}_{2}_{3}_dlst_{4:.1f}s_flags_{5}_comb_{6:0d}.pdf'.format(infile_no_ext, xcpdps2_b_avg[sampling]['z'][zind], stat, sampling, 3.6e3*xcpdps2_b_avg['dlst'][0], applyflags_str, combi))