Full-Day Autocorrelation Checking¶

by Josh Dillon and Steven Murray, last updated July 13, 2023

This notebook is designed to assess per-day data quality from just autocorrelations, enabling a quick assessment of whether the day is worth pushing through further analysis. In particular, it is designed to look for times that are particularly contaminated by broadband RFI (e.g. from lightning), picking out fraction of days worth analyzing. It's output is a an a priori flag yaml file readable by hera_qm.metrics_io functions read_a_priori_chan_flags(), read_a_priori_int_flags(), and read_a_priori_ant_flags().

Here's a set of links to skip to particular figures:

• Figure 1: Preliminary Array Flag Fraction Summary¶

• Figure 2: z-Score of DPSS-Filtered, Averaged Good Autocorrelation and Initial Flags¶

• Figure 3: Proposed A Priori Time Flags Based on Frequency-Averaged and Convolved z-Score Magnitude¶

In [1]:
import time
tstart = time.time()
In [2]:
import os
os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE'
import h5py
import hdf5plugin  # REQUIRED to have the compression plugins available
import numpy as np
import pandas as pd
import glob
import os
import matplotlib.pyplot as plt
import matplotlib
import copy
from scipy.ndimage import convolve
from hera_cal import io, utils
from hera_cal.smooth_cal import dpss_filters, solve_2D_DPSS
from hera_qm import ant_class, xrfi, metrics_io
from hera_qm.time_series_metrics import true_stretches, impose_max_flag_gap, metric_convolution_flagging
from hera_filters import dspec
import warnings

from IPython.display import display, HTML
%matplotlib inline
display(HTML("<style>.container { width:100% !important; }</style>"))
_ = np.seterr(all='ignore')  # get rid of red warnings
%config InlineBackend.figure_format = 'retina'

Parse input and output files¶

In [3]:
# get filenames
AUTO_FILE = os.environ.get("AUTO_FILE", None)
# AUTO_FILE = '/lustre/aoc/projects/hera/h6c-analysis/IDR2/2459869/zen.2459869.25299.sum.autos.uvh5'

SUM_AUTOS_SUFFIX = os.environ.get("SUM_AUTOS_SUFFIX", 'sum.autos.uvh5')
DIFF_AUTOS_SUFFIX = os.environ.get("DIFF_AUTOS_SUFFIX", 'diff.autos.uvh5')
OUT_YAML_SUFFIX = os.environ.get("OUT_YAML_SUFFIX", '_apriori_flags.yaml')
OUT_YAML_DIR = os.environ.get("OUT_YAML_DIR", None)
# OUT_YAML_DIR = '/lustre/aoc/projects/hera/jsdillon/H6C/'

if OUT_YAML_DIR is None:
    OUT_YAML_DIR = os.path.dirname(AUTO_FILE)

auto_sums_glob = '.'.join(AUTO_FILE.split('.')[:-4]) + '.*.' + SUM_AUTOS_SUFFIX
auto_diffs_glob = auto_sums_glob.replace(SUM_AUTOS_SUFFIX, DIFF_AUTOS_SUFFIX)
out_yaml_file = os.path.join(OUT_YAML_DIR, AUTO_FILE.split('.')[-5] + OUT_YAML_SUFFIX)
In [4]:
auto_sums = sorted(glob.glob(auto_sums_glob))
print(f'Found {len(auto_sums)} *.{SUM_AUTOS_SUFFIX} files starting with {auto_sums[0]}.')

auto_diffs = sorted(glob.glob(auto_diffs_glob))
print(f'Found {len(auto_diffs)} *.{DIFF_AUTOS_SUFFIX} files starting with {auto_diffs[0]}.')

Found 1851 *.sum.autos.uvh5 files starting with /mnt/sn1/2460231/zen.2460231.25447.sum.autos.uvh5.
Found 1851 *.diff.autos.uvh5 files starting with /mnt/sn1/2460231/zen.2460231.25447.diff.autos.uvh5.

Parse settings¶

In [5]:
# bounds on zeros in spectra
good_zeros_per_eo_spectrum = (0, int(os.environ.get("MAX_ZEROS_PER_EO_SPEC_GOOD", 2)))
suspect_zeros_per_eo_spectrum = (0, int(os.environ.get("MAX_ZEROS_PER_EO_SPEC_SUSPECT", 8)))

# bounds on autocorrelation power
auto_power_good = (float(os.environ.get("AUTO_POWER_GOOD_LOW", 5)), float(os.environ.get("AUTO_POWER_GOOD_HIGH", 30)))
auto_power_suspect = (float(os.environ.get("AUTO_POWER_SUSPECT_LOW", 1)), float(os.environ.get("AUTO_POWER_SUSPECT_HIGH", 60)))

# bounds on autocorrelation slope
auto_slope_good = (float(os.environ.get("AUTO_SLOPE_GOOD_LOW", -0.4)), float(os.environ.get("AUTO_SLOPE_GOOD_HIGH", 0.4)))
auto_slope_suspect = (float(os.environ.get("AUTO_SLOPE_SUSPECT_LOW", -0.6)), float(os.environ.get("AUTO_SLOPE_SUSPECT_HIGH", 0.6)))

# bounds on autocorrelation RFI
auto_rfi_good = (0, float(os.environ.get("AUTO_RFI_GOOD", 0.015)))
auto_rfi_suspect = (0, float(os.environ.get("AUTO_RFI_SUSPECT", 0.03)))

# bounds on autocorrelation shape
auto_shape_good = (0, float(os.environ.get("AUTO_SHAPE_GOOD", 0.1)))
auto_shape_suspect = (0, float(os.environ.get("AUTO_SHAPE_SUSPECT", 0.2)))

# parse RFI settings for antenna flagging
RFI_DPSS_HALFWIDTH = float(os.environ.get("RFI_DPSS_HALFWIDTH", 300e-9)) # in s
RFI_NSIG = float(os.environ.get("RFI_NSIG", 6))

# parse settings for identifying mislabeled data by X-engine
BAD_XENGINE_ZCUT = float(os.environ.get("BAD_XENGINE_ZCUT", 10))

# DPSS settings for full-day filtering
FREQ_FILTER_SCALE = float(os.environ.get("FREQ_FILTER_SCALE", 5.0)) # in MHz
TIME_FILTER_SCALE = float(os.environ.get("TIME_FILTER_SCALE", 450.0)) # in s
EIGENVAL_CUTOFF = float(os.environ.get("EIGENVAL_CUTOFF", 1e-12))

# A priori flag settings
FM_LOW_FREQ = float(os.environ.get("FM_LOW_FREQ", 87.5)) # in MHz
FM_HIGH_FREQ = float(os.environ.get("FM_HIGH_FREQ", 108.0)) # in MHz
FM_freq_range = [FM_LOW_FREQ * 1e6, FM_HIGH_FREQ * 1e6]
MAX_SOLAR_ALT = float(os.environ.get("MAX_SOLAR_ALT", 0.0)) # in degrees
SMOOTHED_ABS_Z_THRESH = float(os.environ.get("SMOOTHED_ABS_Z_THRESH", 10))
WHOLE_DAY_FLAG_THRESH = float(os.environ.get("WHOLE_DAY_FLAG_THRESH", 0.5))

for setting in ['good_zeros_per_eo_spectrum', 'suspect_zeros_per_eo_spectrum', 'auto_power_good', 'auto_power_suspect', 
                'auto_slope_good', 'auto_slope_suspect', 'auto_rfi_good', 'auto_rfi_suspect',
                'auto_shape_good', 'auto_shape_suspect', 'BAD_XENGINE_ZCUT', 'RFI_DPSS_HALFWIDTH', 'RFI_NSIG',
                'FREQ_FILTER_SCALE', 'TIME_FILTER_SCALE', 'EIGENVAL_CUTOFF', 'FM_freq_range',
                'MAX_SOLAR_ALT', 'SMOOTHED_ABS_Z_THRESH', 'WHOLE_DAY_FLAG_THRESH']:
        print(f'{setting} = {eval(setting)}')

good_zeros_per_eo_spectrum = (0, 2)
suspect_zeros_per_eo_spectrum = (0, 8)
auto_power_good = (5.0, 30.0)
auto_power_suspect = (1.0, 60.0)
auto_slope_good = (-0.4, 0.4)
auto_slope_suspect = (-0.6, 0.6)
auto_rfi_good = (0, 0.015)
auto_rfi_suspect = (0, 0.03)
auto_shape_good = (0, 0.1)
auto_shape_suspect = (0, 0.2)
BAD_XENGINE_ZCUT = 10.0
RFI_DPSS_HALFWIDTH = 3e-07
RFI_NSIG = 6.0
FREQ_FILTER_SCALE = 5.0
TIME_FILTER_SCALE = 450.0
EIGENVAL_CUTOFF = 1e-12
FM_freq_range = [87500000.0, 108000000.0]
MAX_SOLAR_ALT = 0.0
SMOOTHED_ABS_Z_THRESH = 10.0
WHOLE_DAY_FLAG_THRESH = 0.5

Classify Antennas and Find RFI Per-File¶

In [6]:
cache = {}
In [7]:
def classify_autos_and_preliminary_rfi(auto_sum_file, auto_diff_file):
    
    hd_sum = io.HERADataFastReader(auto_sum_file)
    sum_autos, _, _ = hd_sum.read(read_flags=False, read_nsamples=False)
    hd_diff = io.HERADataFastReader(auto_diff_file)
    diff_autos, _, _ = hd_diff.read(read_flags=False, read_nsamples=False)    
    ants = sorted(set([ant for bl in hd_sum.bls for ant in utils.split_bl(bl)]))
    
    zeros_class = ant_class.even_odd_zeros_checker(sum_autos, diff_autos, good=good_zeros_per_eo_spectrum, suspect=suspect_zeros_per_eo_spectrum)
    auto_power_class = ant_class.auto_power_checker(sum_autos, good=auto_power_good, suspect=auto_power_suspect)
    auto_slope_class = ant_class.auto_slope_checker(sum_autos, good=auto_slope_good, suspect=auto_slope_suspect, edge_cut=100, filt_size=17)
    auto_rfi_class = ant_class.auto_rfi_checker(sum_autos, good=auto_rfi_good, suspect=auto_rfi_suspect, 
                                                filter_half_widths=[RFI_DPSS_HALFWIDTH], nsig=RFI_NSIG, cache=cache)
    auto_class = auto_power_class + auto_slope_class + auto_rfi_class    
    
    final_class = zeros_class + auto_class
    
    if len(final_class.good_ants) + len(final_class.suspect_ants) > 0:
        # Compute int_count for all unflagged autocorrelations averaged together
        int_time = 24 * 3600 * np.median(np.diff(sum_autos.times_by_bl[hd_sum.bls[0][0:2]]))
        chan_res = np.median(np.diff(sum_autos.freqs))
        int_count = int(int_time * chan_res) * (len(final_class.good_ants) + len(final_class.suspect_ants))
    
        avg_auto = {(-1, -1, 'ee'): np.mean([sum_autos[bl] for bl in hd_sum.bls if final_class[utils.split_bl(bl)[0]] != 'bad'], axis=0)}
        
        # Flag RFI first with channel differences and then with DPSS
        antenna_flags, _ = xrfi.flag_autos(avg_auto, int_count=int_count, nsig=(RFI_NSIG * 5))
        _, rfi_flags = xrfi.flag_autos(avg_auto, int_count=int_count, flag_method='dpss_flagger',
                                       flags=antenna_flags, freqs=sum_autos.freqs, filter_centers=[0],
                                       filter_half_widths=[RFI_DPSS_HALFWIDTH], eigenval_cutoff=[1e-9], nsig=RFI_NSIG)
        
        # Classify on auto shapes
        auto_shape_class = ant_class.auto_shape_checker(sum_autos, good=auto_shape_good, suspect=auto_shape_suspect,
                                                flag_spectrum=np.sum(rfi_flags, axis=0).astype(bool), antenna_class=final_class)
        final_class += auto_shape_class
        
        # Classify on excess power in X-engine diffs
        xengine_diff_class = ant_class.non_noiselike_diff_by_xengine_checker(sum_autos, diff_autos, flag_waterfall=rfi_flags, 
                                                                             antenna_class=final_class, 
                                                                             xengine_chans=96, bad_xengine_zcut=BAD_XENGINE_ZCUT)
        final_class += xengine_diff_class

    else:
        rfi_flags = np.ones((len(sum_autos.times), len(sum_autos.freqs)), dtype=bool)

    return final_class, rfi_flags
In [8]:
classifications = {}
preliminary_rfi_flags = {}
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    for auto_sum_file, auto_diff_file in list(zip(auto_sums, auto_diffs)):
        classifications[auto_sum_file], preliminary_rfi_flags[auto_sum_file] = classify_autos_and_preliminary_rfi(auto_sum_file, auto_diff_file)
In [9]:
all_ants = set([ant for auto_sum_file in classifications for ant in classifications[auto_sum_file]])
ant_flag_fracs = {ant: 0 for ant in all_ants}
for classification in classifications.values():
    for ant in classification.bad_ants:
        ant_flag_fracs[ant] += 1
ant_flag_fracs = {ant: ant_flag_fracs[ant] / len(classifications) for ant in all_ants}
In [10]:
def flag_frac_array_plot():
    fig, axes = plt.subplots(1, 2, figsize=(14, 8), dpi=100, gridspec_kw={'width_ratios': [2, 1]})

    def flag_frac_panel(ax, antnums, radius=7, legend=False):

        ang_dict = {'Jee': (225, 405), 'Jnn': (45, 225)}
        hd_sum = io.HERADataFastReader(auto_sums[-1])
        xpos = np.array([hd_sum.antpos[ant[0]][0] for ant in all_ants if ant[0] in antnums])
        ypos = np.array([hd_sum.antpos[ant[0]][1] for ant in all_ants if ant[0] in antnums])
        scatter = ax.scatter(xpos, ypos, c='w', s=0)
        for ant in all_ants:
            antnum, pol = ant
            if antnum in antnums:
                ax.add_artist(matplotlib.patches.Wedge(tuple(hd_sum.antpos[antnum][0:2]), radius, *ang_dict[pol], color='grey'))
                flag_frac = ant_flag_fracs[ant]
                if flag_frac > .05:
                    ax.add_artist(matplotlib.patches.Wedge(tuple(hd_sum.antpos[antnum][0:2]), radius * np.sqrt(flag_frac), *ang_dict[pol], color='r'))
                ax.text(hd_sum.antpos[antnum][0], hd_sum.antpos[antnum][1], str(antnum), color='w',  va='center', ha='center', zorder=100)

        ax.axis('equal')
        ax.set_xlim([np.min(xpos) - radius * 2, np.max(xpos) + radius * 2])
        ax.set_ylim([np.min(ypos) - radius * 2, np.max(ypos) + radius * 2])
        ax.set_xlabel("East-West Position (meters)", size=12)
        ax.set_ylabel("North-South Position (meters)", size=12)

        if legend:
            legend_objs = []
            legend_labels = []

            legend_objs.append(matplotlib.lines.Line2D([0], [0], marker='o', color='w', markeredgecolor='grey', markerfacecolor='grey', markersize=15))
            unflagged_nights = lambda pol: np.sum([1 - ant_flag_fracs[ant] for ant in all_ants if ant[-1] == pol])
            legend_labels.append((' \u2571\n').join([f'{unflagged_nights(pol):.1f} unflagged {pol[-1]}-polarized\nantenna-nights.' for pol in ['Jee', 'Jnn']]))

            legend_objs.append(matplotlib.lines.Line2D([0], [0], marker='o', color='w', markeredgecolor='red', markerfacecolor='red', markersize=15))
            unflagged_nights = lambda pol: np.sum([ant_flag_fracs[ant] for ant in all_ants if ant[-1] == pol])
            legend_labels.append((' \u2571\n').join([f'{unflagged_nights(pol):.1f} flagged {pol[-1]}-polarized\nantenna-nights.' for pol in ['Jee', 'Jnn']]))        
            ax.legend(legend_objs, legend_labels, ncol=1, fontsize=12)

    flag_frac_panel(axes[0], sorted(set([ant[0] for ant in all_ants if ant[0] < 320])), radius=7)
    flag_frac_panel(axes[1], sorted(set([ant[0] for ant in all_ants if ant[0] >= 320])), radius=50, legend=True)
    plt.tight_layout()

Figure 1: Preliminary Array Flag Fraction Summary¶

Per-antenna flagging fraction of data based purely on metrics that only use autocorrelations. This is likely an underestimate of flags, since it ignores low correlation, cross-polarized antennas, and high redcal $\chi^2$, among other factors.

In [11]:
flag_frac_array_plot()

Load and Average Unflagged Autocorrelations¶

In [12]:
min_flag_frac = np.min(list(ant_flag_fracs.values()))
least_flagged_ants = sorted([ant for ant in all_ants if ant_flag_fracs[ant] == min_flag_frac])
least_flagged_autos = [utils.join_bl(ant, ant) for ant in least_flagged_ants]
In [13]:
avg_autos = {}
times = []
for auto_sum_file in auto_sums:
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
    hd_sum = io.HERADataFastReader(auto_sum_file)
    sum_autos, _, _ = hd_sum.read(bls=least_flagged_autos, read_flags=False, read_nsamples=False)
    avg_autos[auto_sum_file] = np.mean([sum_autos[bl] for bl in sum_autos], axis=0)
    times.extend(hd_sum.times)
times = np.array(times)
freqs = hd_sum.freqs
In [14]:
flags = np.vstack([flag for flag in preliminary_rfi_flags.values()])
avg_auto = np.vstack([auto for auto in avg_autos.values()])
In [15]:
def flag_FM(flags, freqs, freq_range=[87.5e6, 108e6]):
    '''Apply flags to all frequencies within freq_range (in Hz).'''
    flags[:, np.logical_and(freqs >= freq_range[0], freqs <= freq_range[1])] = True 
    
def flag_sun(flags, times, max_solar_alt=0):
    '''Apply flags to all times where the solar altitude is greater than max_solar_alt (in degrees).'''
    solar_altitudes_degrees = utils.get_sun_alt(times)
    flags[solar_altitudes_degrees >= max_solar_alt, :] = True
In [16]:
flag_FM(flags, freqs, freq_range=FM_freq_range)
solar_flags = np.zeros_like(flags)
flag_sun(solar_flags, times, max_solar_alt=MAX_SOLAR_ALT)
flags |= solar_flags
all_flagged_times = np.all(flags, axis=1)

DPSS Filter Average Autocorrlations¶

In [17]:
def predict_auto_noise(auto, dt, df, nsamples=1):
    '''Predict noise on an (antenna-averaged) autocorrelation. The product of Delta t and Delta f
    must be unitless. For N autocorrelations averaged together, use nsamples=N.'''
    int_count = int(dt * df) * nsamples
    return np.abs(auto) / np.sqrt(int_count / 2)

# Figure out noise and weights
int_time = 24 * 3600 * np.median(np.diff(times))
chan_res = np.median(np.diff(freqs))
noise = predict_auto_noise(avg_auto, int_time, chan_res, nsamples=len(least_flagged_ants))
wgts = np.where(flags, 0, noise**-2)
In [18]:
time_filters, freq_filters = dpss_filters(freqs=freqs, # Hz
                                          times=times, # JD
                                          freq_scale=FREQ_FILTER_SCALE,
                                          time_scale=TIME_FILTER_SCALE,
                                          eigenval_cutoff=EIGENVAL_CUTOFF)
In [19]:
model, fit_info = solve_2D_DPSS(avg_auto, wgts, time_filters, freq_filters, method='lu_solve')

WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
In [20]:
noise_model = predict_auto_noise(np.abs(model), int_time, chan_res, nsamples=len(least_flagged_ants))
zscore = ((avg_auto - model) / noise_model).real
In [21]:
def plot_z_score(flags, zscore):
    plt.figure(figsize=(14,10), dpi=100)
    extent = [freqs[0]/1e6, freqs[-1]/1e6, times[-1] - int(times[0]), times[0] - int(times[0])]
    plt.imshow(np.where(flags, np.nan, zscore.real), aspect='auto', cmap='bwr', interpolation='none', vmin=-SMOOTHED_ABS_Z_THRESH, vmax=SMOOTHED_ABS_Z_THRESH, extent=extent)
    plt.colorbar(location='top', label='z score', extend='both')
    plt.xlabel('Frequency (MHz)')
    plt.ylabel(f'JD - {int(times[0])}')
    plt.tight_layout()

Figure 2: z-Score of DPSS-Filtered, Averaged Good Autocorrelation and Initial Flags¶

This plot shows the z-score of a DPSS-filtered, deeply averaged autocorrelation, where the noise is inferred from the integration time, channel width, and DPSS model. DPSS was performed using the per-file RFI flagging analogous to that used in the file_calibration notebook, which is generally insensitive to broadband RFI.

In [22]:
plot_z_score(flags, zscore)

Find Bad Time Ranges¶

In [23]:
# Propose flags for night times when there's too much temporal structure due to, e.g. lightning
sigma = TIME_FILTER_SCALE / int_time
metric = np.nanmean(np.where(flags, np.nan, np.abs(zscore)), axis=1)
kernel = np.exp(-np.arange(-len(metric) // 2, len(metric) // 2 + 1)**2 / 2 / sigma**2)
kernel /= np.sum(kernel)
convolved_metric = np.full_like(metric, np.nan)
convolved_metric[~all_flagged_times] = convolve(metric[~all_flagged_times], kernel, mode='reflect')
apriori_time_flags = np.ones_like(metric, dtype=bool)
apriori_time_flags[~all_flagged_times] = metric_convolution_flagging(metric[~all_flagged_times], convolved_metric[~all_flagged_times] >= SMOOTHED_ABS_Z_THRESH, 
                                                               [0, SMOOTHED_ABS_Z_THRESH], sigma=(TIME_FILTER_SCALE / int_time), max_flag_gap=0)

Mean of empty slice
In [24]:
# Flag whole day if too much of the day is flagged
apriori_flag_frac = np.mean(apriori_time_flags[~all_flagged_times])
if apriori_flag_frac > WHOLE_DAY_FLAG_THRESH:
    print(f'A priori time flag fraction of {apriori_flag_frac:.2%} is greater than {WHOLE_DAY_FLAG_THRESH:.2%}... Flagging whole day.')
    apriori_time_flags = np.ones_like(apriori_time_flags)
In [25]:
def apriori_flag_plot():
    plt.figure(figsize=(14, 6))
    plt.semilogy(times - int(times[0]), np.nanmean(np.where(flags, np.nan, np.abs(zscore)), axis=1), label='Average |z| Over Frequency')
    plt.semilogy(times - int(times[0]), convolved_metric, label=f'Convolved on {TIME_FILTER_SCALE}-second timescale')
    plt.axhline(SMOOTHED_ABS_Z_THRESH, color='k', ls='--', label='Threshold on Convolved |z|')
    for i, apf in enumerate(true_stretches(apriori_time_flags)):
        plt.axvspan((times - int(times[0]))[apf.start], (times - int(times[0]))[apf.stop - 1], color='r', zorder=0, alpha=.3, 
                    label=('Proposed A Priori Flags' if i == 0 else None))
    plt.legend()

    plt.xlabel(f'JD - {int(times[0])}')
    plt.ylabel('Frequency Averaged |z-score|')
    plt.tight_layout()

Figure 3: Proposed A Priori Time Flags Based on Frequency-Averaged and Convolved z-Score Magnitude¶

This plot shows the average (over frequency) magnitude of z-scores as a function of time. This metric is smoothed to pick out ranges of times where the DPSS residual reveals persistent temporal structure. Flags due to the sun being above the horizon are also shown. The unflagged range of times is required to be contiguous.

In [26]:
apriori_flag_plot()

Mean of empty slice

Write a priori flags to a yaml¶

Also writing as a priori flags channels that are 100% flagged and antennas that are 100% flagged.

In [27]:
dt = int_time / 3600 / 24 / 2
out_yml_str = 'JD_flags: ' + str([[times[apf][0] - dt, times[apf][-1] + dt] for apf in true_stretches(apriori_time_flags)])
out_yml_str += '\n\nfreq_flags: ' + str([[freqs[apf][0] - chan_res / 2, freqs[apf][-1] + chan_res / 2] for apf in true_stretches(np.all(flags, axis=0))])
out_yml_str += '\n\nex_ants: ' + str(sorted([ant for ant in all_ants if ant_flag_fracs[ant] == 1])).replace("'", "").replace('(', '[').replace(')', ']')
out_yml_str += '\n\nall_ant: ' + str(sorted(all_ants)).replace("'", "").replace('(', '[').replace(')', ']')
out_yml_str += f'\n\njd_range: [{times.min()}, {times.max()}]'
out_yml_str += f'\n\nfreq_range: [{freqs.min()}, {freqs.max()}]'

print(f'Writing the following to {out_yaml_file}\n' + '-' * (25 + len(out_yaml_file)))
print(out_yml_str)
with open(out_yaml_file, 'w') as outfile:
    outfile.writelines(out_yml_str)

Writing the following to /mnt/sn1/2460231/2460231_apriori_flags.yaml
--------------------------------------------------------------------
JD_flags: [[2460231.667527605, 2460231.66842239]]

freq_flags: [[62362670.8984375, 62606811.5234375], [69931030.2734375, 70053100.5859375], [87509155.2734375, 108016967.7734375], [124984741.2109375, 125106811.5234375], [141464233.3984375, 141586303.7109375], [142074584.9609375, 142318725.5859375], [143783569.3359375, 143905639.6484375], [147445678.7109375, 147567749.0234375], [154281616.2109375, 154403686.5234375], [175155639.6484375, 175277709.9609375], [187362670.8984375, 187606811.5234375], [191146850.5859375, 191268920.8984375], [201766967.7734375, 201889038.0859375], [223129272.4609375, 223373413.0859375], [229232788.0859375, 229354858.3984375]]

ex_ants: [[3, Jnn], [15, Jee], [16, Jnn], [18, Jee], [18, Jnn], [27, Jee], [27, Jnn], [28, Jnn], [29, Jee], [30, Jee], [30, Jnn], [32, Jee], [45, Jee], [45, Jnn], [46, Jee], [46, Jnn], [54, Jee], [63, Jee], [63, Jnn], [64, Jee], [64, Jnn], [73, Jee], [73, Jnn], [77, Jee], [78, Jee], [78, Jnn], [84, Jee], [84, Jnn], [85, Jee], [85, Jnn], [86, Jee], [86, Jnn], [87, Jee], [87, Jnn], [93, Jnn], [96, Jnn], [98, Jnn], [99, Jee], [99, Jnn], [101, Jee], [101, Jnn], [102, Jee], [102, Jnn], [103, Jee], [103, Jnn], [104, Jee], [104, Jnn], [106, Jee], [109, Jnn], [110, Jee], [112, Jee], [115, Jee], [119, Jnn], [120, Jee], [120, Jnn], [121, Jee], [121, Jnn], [122, Jee], [122, Jnn], [123, Jee], [123, Jnn], [125, Jee], [159, Jnn], [161, Jnn], [165, Jee], [170, Jee], [174, Jee], [174, Jnn], [176, Jee], [180, Jnn], [182, Jnn], [183, Jnn], [184, Jee], [194, Jee], [194, Jnn], [200, Jee], [203, Jee], [203, Jnn], [208, Jnn], [209, Jee], [209, Jnn], [211, Jnn], [214, Jee], [216, Jee], [216, Jnn], [217, Jee], [218, Jnn], [219, Jee], [219, Jnn], [226, Jnn], [232, Jee], [232, Jnn], [245, Jnn], [246, Jee], [246, Jnn], [251, Jee], [253, Jee], [253, Jnn], [261, Jee], [261, Jnn], [262, Jee], [262, Jnn], [269, Jnn], [282, Jee], [282, Jnn], [283, Jee], [283, Jnn], [295, Jee], [323, Jee], [328, Jnn], [329, Jee], [329, Jnn], [332, Jee], [332, Jnn]]

all_ant: [[3, Jee], [3, Jnn], [4, Jee], [4, Jnn], [5, Jee], [5, Jnn], [7, Jee], [7, Jnn], [8, Jee], [8, Jnn], [9, Jee], [9, Jnn], [10, Jee], [10, Jnn], [15, Jee], [15, Jnn], [16, Jee], [16, Jnn], [17, Jee], [17, Jnn], [18, Jee], [18, Jnn], [19, Jee], [19, Jnn], [20, Jee], [20, Jnn], [21, Jee], [21, Jnn], [22, Jee], [22, Jnn], [27, Jee], [27, Jnn], [28, Jee], [28, Jnn], [29, Jee], [29, Jnn], [30, Jee], [30, Jnn], [31, Jee], [31, Jnn], [32, Jee], [32, Jnn], [33, Jee], [33, Jnn], [34, Jee], [34, Jnn], [35, Jee], [35, Jnn], [36, Jee], [36, Jnn], [37, Jee], [37, Jnn], [38, Jee], [38, Jnn], [40, Jee], [40, Jnn], [41, Jee], [41, Jnn], [42, Jee], [42, Jnn], [45, Jee], [45, Jnn], [46, Jee], [46, Jnn], [47, Jee], [47, Jnn], [48, Jee], [48, Jnn], [49, Jee], [49, Jnn], [50, Jee], [50, Jnn], [51, Jee], [51, Jnn], [52, Jee], [52, Jnn], [53, Jee], [53, Jnn], [54, Jee], [54, Jnn], [55, Jee], [55, Jnn], [56, Jee], [56, Jnn], [57, Jee], [57, Jnn], [61, Jee], [61, Jnn], [62, Jee], [62, Jnn], [63, Jee], [63, Jnn], [64, Jee], [64, Jnn], [65, Jee], [65, Jnn], [66, Jee], [66, Jnn], [67, Jee], [67, Jnn], [68, Jee], [68, Jnn], [69, Jee], [69, Jnn], [70, Jee], [70, Jnn], [71, Jee], [71, Jnn], [72, Jee], [72, Jnn], [73, Jee], [73, Jnn], [77, Jee], [77, Jnn], [78, Jee], [78, Jnn], [79, Jee], [79, Jnn], [80, Jee], [80, Jnn], [81, Jee], [81, Jnn], [82, Jee], [82, Jnn], [83, Jee], [83, Jnn], [84, Jee], [84, Jnn], [85, Jee], [85, Jnn], [86, Jee], [86, Jnn], [87, Jee], [87, Jnn], [88, Jee], [88, Jnn], [89, Jee], [89, Jnn], [90, Jee], [90, Jnn], [91, Jee], [91, Jnn], [92, Jee], [92, Jnn], [93, Jee], [93, Jnn], [94, Jee], [94, Jnn], [95, Jee], [95, Jnn], [96, Jee], [96, Jnn], [97, Jee], [97, Jnn], [98, Jee], [98, Jnn], [99, Jee], [99, Jnn], [100, Jee], [100, Jnn], [101, Jee], [101, Jnn], [102, Jee], [102, Jnn], [103, Jee], [103, Jnn], [104, Jee], [104, Jnn], [105, Jee], [105, Jnn], [106, Jee], [106, Jnn], [107, Jee], [107, Jnn], [108, Jee], [108, Jnn], [109, Jee], [109, Jnn], [110, Jee], [110, Jnn], [111, Jee], [111, Jnn], [112, Jee], [112, Jnn], [113, Jee], [113, Jnn], [114, Jee], [114, Jnn], [115, Jee], [115, Jnn], [116, Jee], [116, Jnn], [117, Jee], [117, Jnn], [118, Jee], [118, Jnn], [119, Jee], [119, Jnn], [120, Jee], [120, Jnn], [121, Jee], [121, Jnn], [122, Jee], [122, Jnn], [123, Jee], [123, Jnn], [124, Jee], [124, Jnn], [125, Jee], [125, Jnn], [126, Jee], [126, Jnn], [127, Jee], [127, Jnn], [128, Jee], [128, Jnn], [129, Jee], [129, Jnn], [130, Jee], [130, Jnn], [131, Jee], [131, Jnn], [132, Jee], [132, Jnn], [133, Jee], [133, Jnn], [134, Jee], [134, Jnn], [135, Jee], [135, Jnn], [136, Jee], [136, Jnn], [137, Jee], [137, Jnn], [138, Jee], [138, Jnn], [139, Jee], [139, Jnn], [140, Jee], [140, Jnn], [141, Jee], [141, Jnn], [142, Jee], [142, Jnn], [143, Jee], [143, Jnn], [144, Jee], [144, Jnn], [145, Jee], [145, Jnn], [146, Jee], [146, Jnn], [147, Jee], [147, Jnn], [148, Jee], [148, Jnn], [149, Jee], [149, Jnn], [150, Jee], [150, Jnn], [151, Jee], [151, Jnn], [152, Jee], [152, Jnn], [153, Jee], [153, Jnn], [154, Jee], [154, Jnn], [155, Jee], [155, Jnn], [156, Jee], [156, Jnn], [157, Jee], [157, Jnn], [158, Jee], [158, Jnn], [159, Jee], [159, Jnn], [160, Jee], [160, Jnn], [161, Jee], [161, Jnn], [162, Jee], [162, Jnn], [163, Jee], [163, Jnn], [164, Jee], [164, Jnn], [165, Jee], [165, Jnn], [166, Jee], [166, Jnn], [167, Jee], [167, Jnn], [168, Jee], [168, Jnn], [169, Jee], [169, Jnn], [170, Jee], [170, Jnn], [171, Jee], [171, Jnn], [172, Jee], [172, Jnn], [173, Jee], [173, Jnn], [174, Jee], [174, Jnn], [175, Jee], [175, Jnn], [176, Jee], [176, Jnn], [177, Jee], [177, Jnn], [178, Jee], [178, Jnn], [179, Jee], [179, Jnn], [180, Jee], [180, Jnn], [181, Jee], [181, Jnn], [182, Jee], [182, Jnn], [183, Jee], [183, Jnn], [184, Jee], [184, Jnn], [185, Jee], [185, Jnn], [186, Jee], [186, Jnn], [187, Jee], [187, Jnn], [189, Jee], [189, Jnn], [190, Jee], [190, Jnn], [191, Jee], [191, Jnn], [192, Jee], [192, Jnn], [193, Jee], [193, Jnn], [194, Jee], [194, Jnn], [195, Jee], [195, Jnn], [196, Jee], [196, Jnn], [197, Jee], [197, Jnn], [198, Jee], [198, Jnn], [199, Jee], [199, Jnn], [200, Jee], [200, Jnn], [201, Jee], [201, Jnn], [202, Jee], [202, Jnn], [203, Jee], [203, Jnn], [204, Jee], [204, Jnn], [205, Jee], [205, Jnn], [206, Jee], [206, Jnn], [207, Jee], [207, Jnn], [208, Jee], [208, Jnn], [209, Jee], [209, Jnn], [210, Jee], [210, Jnn], [211, Jee], [211, Jnn], [212, Jee], [212, Jnn], [213, Jee], [213, Jnn], [214, Jee], [214, Jnn], [215, Jee], [215, Jnn], [216, Jee], [216, Jnn], [217, Jee], [217, Jnn], [218, Jee], [218, Jnn], [219, Jee], [219, Jnn], [220, Jee], [220, Jnn], [221, Jee], [221, Jnn], [222, Jee], [222, Jnn], [223, Jee], [223, Jnn], [224, Jee], [224, Jnn], [225, Jee], [225, Jnn], [226, Jee], [226, Jnn], [227, Jee], [227, Jnn], [228, Jee], [228, Jnn], [229, Jee], [229, Jnn], [231, Jee], [231, Jnn], [232, Jee], [232, Jnn], [233, Jee], [233, Jnn], [234, Jee], [234, Jnn], [235, Jee], [235, Jnn], [237, Jee], [237, Jnn], [238, Jee], [238, Jnn], [239, Jee], [239, Jnn], [240, Jee], [240, Jnn], [241, Jee], [241, Jnn], [242, Jee], [242, Jnn], [243, Jee], [243, Jnn], [244, Jee], [244, Jnn], [245, Jee], [245, Jnn], [246, Jee], [246, Jnn], [250, Jee], [250, Jnn], [251, Jee], [251, Jnn], [252, Jee], [252, Jnn], [253, Jee], [253, Jnn], [261, Jee], [261, Jnn], [262, Jee], [262, Jnn], [266, Jee], [266, Jnn], [267, Jee], [267, Jnn], [268, Jee], [268, Jnn], [269, Jee], [269, Jnn], [281, Jee], [281, Jnn], [282, Jee], [282, Jnn], [283, Jee], [283, Jnn], [295, Jee], [295, Jnn], [320, Jee], [320, Jnn], [321, Jee], [321, Jnn], [323, Jee], [323, Jnn], [324, Jee], [324, Jnn], [325, Jee], [325, Jnn], [326, Jee], [326, Jnn], [327, Jee], [327, Jnn], [328, Jee], [328, Jnn], [329, Jee], [329, Jnn], [331, Jee], [331, Jnn], [332, Jee], [332, Jnn], [333, Jee], [333, Jnn], [336, Jee], [336, Jnn], [340, Jee], [340, Jnn]]

jd_range: [2460231.254416623, 2460231.6683664657]

freq_range: [46920776.3671875, 234298706.0546875]
In [28]:
# check output yaml file
flagged_indices = metrics_io.read_a_priori_int_flags(out_yaml_file, times)
for i, apf in enumerate(apriori_time_flags):
    if i in flagged_indices:
        assert apf
    else:
        assert not apf
flagged_chans = metrics_io.read_a_priori_chan_flags(out_yaml_file, freqs=freqs)
for chan in range(len(freqs)):
    if chan in flagged_chans:
        assert np.all(flags[:, chan])
    else:
        assert not np.all(flags[:, chan])
flagged_ants = metrics_io.read_a_priori_ant_flags(out_yaml_file)
for ant in all_ants:
    if ant_flag_fracs[ant] == 1:
        assert ant in flagged_ants
    else:
        assert ant not in flagged_ants

Metadata¶

In [29]:
for repo in ['hera_cal', 'hera_qm', 'hera_filters', 'hera_notebook_templates', 'pyuvdata']:
    exec(f'from {repo} import __version__')
    print(f'{repo}: {__version__}')

hera_cal: 3.4.0
hera_qm: 2.1.3.dev4+g50af006
hera_filters: 0.1.0
hera_notebook_templates: 0.1.dev486+gfb8560a
pyuvdata: 2.4.0
In [30]:
print(f'Finished execution in {(time.time() - tstart) / 60:.2f} minutes.')

Finished execution in 86.32 minutes.