In [1]:

    

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


    

In [2]:

    

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


    

    




AstroUtils git # 1be75833b358cbc53ff9fd12ec376a5c6335463a
PRISim git # bb083bc1abd3ef44d0b18c094fedfe3dc30ed6ac

In [3]:

    

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


    

In [4]:

    

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, 23, 25]}, '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/models/bispectrum_phase/PS/xCPDPS_collapse_axes_123_EQ28XX_GLEAM_all_RA_1.6_FOV_30_ephemeris_HERA61_noisy_errinfo_incoh_diag_avg.hdf5', 'current_lbl': 'Data', 'other_label': 'Model'}, '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': ['GLEAM_all_RA_1.6_FOV30_ephemeris/RA_1.6_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': [1.6, 2.0], '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': 'Y', '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': [1.5, 2.1], 'triads': None, 'dayind': 0}, 'applyflags': False, 'action': False, 'datastage': 'prelim', 'statistic': 'median', 'sparseness': 10.0}, '2b': {'action': False}, '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': [1.6, 2.0], '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': [511], 'band_center': 150000000.0, 'freq_resolution': 97656.25, 'lstbinsize': 60.0, 'flagants': [50], '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': ['lst_0100_0220_eq28yy.npz'], 'projectdir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/Binned_Data_0100_0220/', 'hdf5_infile': 'EQ28YY_LST_1.6hr.hdf5', 'model_hdf5files': ['EQ28XX_GLEAM_all_RA_1.6_FOV_30_ephemeris_HERA61_noiseless.hdf5', 'EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_GLEAM_all_RA_1.6_FOV_30_ephemeris_HERA61_noiseless.hdf5'], 'datadir': 'EQ28/YY/', 'model_labels': ['FG', 'FG+HI'], 'visfile': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/visibilities/GLEAM_all_RA_1.6_FOV30_ephemeris/RA_1.6_transit_NFbeam_120x10.7sec_1024x97.7kHz/simdata/simvis', 'modeldir': '/lustre/aoc/projects/hera/nthyagar/data/HERA/IDR2.1/ClosurePhase/models/bispectrum_phase/', 'model_npzfiles': ['EQ28XX_GLEAM_all_RA_1.6_FOV_30_ephemeris_HERA61_noiseless.npz', 'EQ28XX_HI_21cmfast_FaintGalaxies_fiducial_GLEAM_all_RA_1.6_FOV_30_ephemeris_HERA61_noiseless.npz'], 'figdir': 'figures/'}, 'telescope': {'latitude': -30.7224, 'longitude': 21.4278}}

In [5]:

    

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

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

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

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


    

In [6]:

    

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

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

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

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


    

In [7]:

    

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


    

    




(1738, 1024, 120)
[1.61665576 1.61963612 1.62261648 1.62559684 1.6285772  1.63155756
 1.63453792 1.63751828 1.64049864 1.643479   1.64645936 1.64943972
 1.65242008 1.65540044 1.6583808  1.66136116 1.66434152 1.66732188
 1.67030224 1.6732826  1.67626296 1.67924332 1.68222368 1.68520404
 1.6881844  1.69116476 1.69414512 1.69712548 1.70010584 1.7030862
 1.70606656 1.70904692 1.71202728 1.71500764 1.717988   1.72096836
 1.72394872 1.72692908 1.72990944 1.7328898  1.73587016 1.73885052
 1.74183088 1.74481124 1.7477916  1.75077196 1.75375232 1.75673268
 1.75971304 1.7626934  1.76567376 1.76865412 1.77163448 1.77461484
 1.7775952  1.78057556 1.78355592 1.78653628 1.78951664 1.792497
 1.79547736 1.79845772 1.80143808 1.80441844 1.8073988  1.81037916
 1.81335952 1.81633988 1.81932024 1.8223006  1.82528096 1.82826132
 1.83124168 1.83422204 1.8372024  1.84018276 1.84316312 1.84614348
 1.84912384 1.8521042  1.85508456 1.85806492 1.86104528 1.86402564
 1.867006   1.86998636 1.87296672 1.87594708 1.87892744 1.8819078
 1.88488816 1.88786852 1.89084888 1.89382924 1.8968096  1.89978996
 1.90277032 1.90575068 1.90873104 1.9117114  1.91469176 1.91767212
 1.92065248 1.92363284 1.9266132  1.92959356 1.93257392 1.93555428
 1.93853464 1.941515   1.94449536 1.94747572 1.95045608 1.95343644
 1.9564168  1.95939716 1.96237752 1.96535788 1.96833823 1.97131859]

In [8]:

    

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


    

In [9]:

    

# 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 [10]:

    

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


    

In [11]:

    

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


    

    




[[1.00281597 1.0026777  1.00257981 ... 1.00426239 1.00424402 1.00424656]
 [1.00580676 1.00566848 1.00557059 ... 1.00725317 1.0072348  1.00723734]
 [1.00879754 1.00865926 1.00856137 ... 1.01024395 1.01022559 1.01022812]
 ...
 [2.19015654 2.19001824 2.18992032 ... 2.19160288 2.19158453 2.1915871 ]
 [2.19314732 2.19300902 2.1929111  ... 2.19459366 2.19457531 2.19457788]
 [2.1961381  2.1959998  2.19590189 ... 2.19758445 2.19756609 2.19756866]]

In [12]:

    

print(cpObj.cpinfo.keys())


    

    




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

In [13]:

    

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


    

    




[u'triads', u'cphase', u'days', u'lst-day', u'lst', u'flags']
(31, 3)
(18,)
(400, 18)
(400, 18)
(400, 18, 31, 1024)

In [14]:

    

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


    

    




('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 [15]:

    

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


    

    




['prelim', 'native']

In [16]:

    

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


    

    




['cphase', 'wts', 'eicp', 'dlstbins', 'daybins', 'lstbins', 'diff_dbins']
['rms', 'median', 'mad', 'mean']
(72,)
[1.01114931 1.02781597 1.04448264 1.06114931 1.07781597 1.09448264
 1.11114931 1.12781597 1.14448264 1.16114931 1.17781597 1.19448264
 1.21114931 1.22781597 1.24448264 1.26114931 1.27781597 1.29448264
 1.31114931 1.32781597 1.34448264 1.36114931 1.37781597 1.39448264
 1.41114931 1.42781597 1.44448264 1.46114931 1.47781597 1.49448264
 1.51114931 1.52781597 1.54448264 1.56114931 1.57781597 1.59448264
 1.61114931 1.62781597 1.64448264 1.66114931 1.67781597 1.69448264
 1.71114931 1.72781597 1.74448264 1.76114931 1.77781597 1.79448264
 1.81114931 1.82781597 1.84448264 1.86114931 1.87781597 1.89448264
 1.91114931 1.92781597 1.94448264 1.96114931 1.97781597 1.99448264
 2.01114931 2.02781597 2.04448264 2.06114931 2.07781597 2.09448264
 2.11114931 2.12781597 2.14448264 2.16114931 2.17781597 2.19448264]
(2,)
[2458102.77777778 2458112.        ]
(72, 2, 31, 1024)
(72, 2, 31, 1024)

In [17]:

    

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


    

    




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

In [18]:

    

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


    

    




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

In [19]:

    

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


    

    




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

In [20]:

    

cpDSobj = BSP.ClosurePhaseDelaySpectrum(cpObj)


    

In [21]:

    

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


    

In [22]:

    

print('freq window centers = ', freq_window_centers)
print('freq window BW eff = ', freq_window_bw)
print('freq window shape = '+freq_window_shape)
print('freq window fftpow = {0:.1f}'.format(freq_window_fftpow))
print('pad = {0:.1f}'.format(pad))


    

    




('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 [23]:

    

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


    

    




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

    

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


    

    




['2d', '2c_err', '3', '2', '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)


    

    




(72, 2, 31, 1024)
[1.6 2. ]
[1.01114931 1.02781597 1.04448264 1.06114931 1.07781597 1.09448264
 1.11114931 1.12781597 1.14448264 1.16114931 1.17781597 1.19448264
 1.21114931 1.22781597 1.24448264 1.26114931 1.27781597 1.29448264
 1.31114931 1.32781597 1.34448264 1.36114931 1.37781597 1.39448264
 1.41114931 1.42781597 1.44448264 1.46114931 1.47781597 1.49448264
 1.51114931 1.52781597 1.54448264 1.56114931 1.57781597 1.59448264
 1.61114931 1.62781597 1.64448264 1.66114931 1.67781597 1.69448264
 1.71114931 1.72781597 1.74448264 1.76114931 1.77781597 1.79448264
 1.81114931 1.82781597 1.84448264 1.86114931 1.87781597 1.89448264
 1.91114931 1.92781597 1.94448264 1.96114931 1.97781597 1.99448264
 2.01114931 2.02781597 2.04448264 2.06114931 2.07781597 2.09448264
 2.11114931 2.12781597 2.14448264 2.16114931 2.17781597 2.19448264]
[36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59]

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

        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)