import time
tstart = time.time()
!hostname

bigmem1.rtp.pvt

import os
os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE'
import h5py
import hdf5plugin  # REQUIRED to have the compression plugins available
import numpy as np
from scipy import constants, interpolate
import copy
import glob
import re
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
pd.set_option('display.max_rows', 1000)
from uvtools.plot import plot_antpos, plot_antclass
from hera_qm import ant_metrics, ant_class, xrfi
from hera_cal import io, utils, redcal, apply_cal, datacontainer, abscal
from hera_filters import dspec
from hera_notebook_templates.data import DATA_PATH as HNBT_DATA
from IPython.display import display, HTML
import linsolve
display(HTML("<style>.container { width:100% !important; }</style>"))
_ = np.seterr(all='ignore')  # get rid of red warnings
%config InlineBackend.figure_format = 'retina'

# this enables better memory management on linux
import ctypes
def malloc_trim():
    try:
        ctypes.CDLL('libc.so.6').malloc_trim(0) 
    except OSError:
        pass

# figure out whether to save results
SAVE_RESULTS = os.environ.get("SAVE_RESULTS", "TRUE").upper() == "TRUE"
SAVE_OMNIVIS_FILE = os.environ.get("SAVE_OMNIVIS_FILE", "FALSE").upper() == "TRUE"

# get infile names
SUM_FILE = os.environ.get("SUM_FILE", None)
# SUM_FILE = '/lustre/aoc/projects/hera/h6c-analysis/IDR2/2459867/zen.2459867.46002.sum.uvh5' # If sum_file is not defined in the environment variables, define it here.
USE_DIFF = os.environ.get("USE_DIFF", "TRUE").upper() == "TRUE"
if USE_DIFF:
    DIFF_FILE = SUM_FILE.replace('sum', 'diff')
else:
    DIFF_FILE = None
    RTP_ANTCLASS = SUM_FILE.replace('.uvh5', '.rtp_ant_class.csv')
VALIDATION = os.environ.get("VALIDATION", "FALSE").upper() == "TRUE"

# get outfilenames
AM_FILE = (SUM_FILE.replace('.uvh5', '.ant_metrics.hdf5') if SAVE_RESULTS else None)
ANTCLASS_FILE = (SUM_FILE.replace('.uvh5', '.ant_class.csv') if SAVE_RESULTS else None)
OMNICAL_FILE = (SUM_FILE.replace('.uvh5', '.omni.calfits') if SAVE_RESULTS else None)
OMNIVIS_FILE = (SUM_FILE.replace('.uvh5', '.omni_vis.uvh5') if SAVE_RESULTS else None)

for fname in ['SUM_FILE', 'DIFF_FILE', 'AM_FILE', 'ANTCLASS_FILE', 'OMNICAL_FILE', 'OMNIVIS_FILE', 
              'SAVE_RESULTS', 'SAVE_OMNIVIS_FILE', 'USE_DIFF', 'VALIDATION']:
    print(f"{fname} = '{eval(fname)}'")

SUM_FILE = '/mnt/sn1/data1/2460743/zen.2460743.45941.sum.uvh5'
DIFF_FILE = '/mnt/sn1/data1/2460743/zen.2460743.45941.diff.uvh5'
AM_FILE = '/mnt/sn1/data1/2460743/zen.2460743.45941.sum.ant_metrics.hdf5'
ANTCLASS_FILE = '/mnt/sn1/data1/2460743/zen.2460743.45941.sum.ant_class.csv'
OMNICAL_FILE = '/mnt/sn1/data1/2460743/zen.2460743.45941.sum.omni.calfits'
OMNIVIS_FILE = '/mnt/sn1/data1/2460743/zen.2460743.45941.sum.omni_vis.uvh5'
SAVE_RESULTS = 'True'
SAVE_OMNIVIS_FILE = 'False'
USE_DIFF = 'True'
VALIDATION = 'False'

# parse plotting settings
PLOT = os.environ.get("PLOT", "TRUE").upper() == "TRUE"
if PLOT:
    %matplotlib inline

# parse omnical settings
OC_MAX_DIMS = int(os.environ.get("OC_MAX_DIMS", 4))
OC_MIN_DIM_SIZE = int(os.environ.get("OC_MIN_DIM_SIZE", 8))
OC_SKIP_OUTRIGGERS = os.environ.get("OC_SKIP_OUTRIGGERS", "FALSE").upper() == "TRUE"
OC_MIN_BL_LEN = float(os.environ.get("OC_MIN_BL_LEN", 1))
OC_MAX_BL_LEN = float(os.environ.get("OC_MAX_BL_LEN", 1e100))
OC_MAXITER = int(os.environ.get("OC_MAXITER", 50))
OC_MAX_RERUN = int(os.environ.get("OC_MAX_RERUN", 4))
OC_RERUN_MAXITER = int(os.environ.get("OC_MAXITER", 25))
OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE = float(os.environ.get("OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE", 1))
OC_USE_PRIOR_SOL = os.environ.get("OC_USE_PRIOR_SOL", "FALSE").upper() == "TRUE"
OC_PRIOR_SOL_FLAG_THRESH = float(os.environ.get("OC_PRIOR_SOL_FLAG_THRESH", .95))
OC_USE_GPU = os.environ.get("SAVE_RESULTS", "FALSE").upper() == "TRUE"

# parse RFI settings
RFI_DPSS_HALFWIDTH = float(os.environ.get("RFI_DPSS_HALFWIDTH", 300e-9))
RFI_NSIG = float(os.environ.get("RFI_NSIG", 4))

# parse abscal settings
ABSCAL_MODEL_FILES_GLOB = os.environ.get("ABSCAL_MODEL_FILES_GLOB", None)
ABSCAL_MIN_BL_LEN = float(os.environ.get("ABSCAL_MIN_BL_LEN", 1.0))
ABSCAL_MAX_BL_LEN = float(os.environ.get("ABSCAL_MAX_BL_LEN", 140.0))
CALIBRATE_CROSS_POLS = os.environ.get("CALIBRATE_CROSS_POLS", "FALSE").upper() == "TRUE"

# print settings
for setting in ['PLOT', 'SAVE_RESULTS', 'OC_MAX_DIMS', 'OC_MIN_DIM_SIZE', 'OC_SKIP_OUTRIGGERS', 
                'OC_MIN_BL_LEN', 'OC_MAX_BL_LEN', 'OC_MAXITER', 'OC_MAX_RERUN', 'OC_RERUN_MAXITER', 
                'OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE', 'OC_USE_PRIOR_SOL', 'OC_PRIOR_SOL_FLAG_THRESH', 
                'OC_USE_GPU', 'RFI_DPSS_HALFWIDTH', 'RFI_NSIG', 'ABSCAL_MODEL_FILES_GLOB', 
                'ABSCAL_MIN_BL_LEN', 'ABSCAL_MAX_BL_LEN', "CALIBRATE_CROSS_POLS"]:
    print(f'{setting} = {eval(setting)}')

PLOT = True
SAVE_RESULTS = True
OC_MAX_DIMS = 4
OC_MIN_DIM_SIZE = 8
OC_SKIP_OUTRIGGERS = True
OC_MIN_BL_LEN = 1.0
OC_MAX_BL_LEN = 140.0
OC_MAXITER = 50
OC_MAX_RERUN = 10
OC_RERUN_MAXITER = 50
OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE = 1.5
OC_USE_PRIOR_SOL = True
OC_PRIOR_SOL_FLAG_THRESH = 0.95
OC_USE_GPU = False
RFI_DPSS_HALFWIDTH = 3e-07
RFI_NSIG = 4.0
ABSCAL_MODEL_FILES_GLOB = None
ABSCAL_MIN_BL_LEN = 1.0
ABSCAL_MAX_BL_LEN = 140.0
CALIBRATE_CROSS_POLS = True

# ant_metrics bounds for low correlation / dead antennas
am_corr_bad = (0, float(os.environ.get("AM_CORR_BAD", 0.3)))
am_corr_suspect = (float(os.environ.get("AM_CORR_BAD", 0.3)), float(os.environ.get("AM_CORR_SUSPECT", 0.5)))

# ant_metrics bounds for cross-polarized antennas
am_xpol_bad = (-1, float(os.environ.get("AM_XPOL_BAD", -0.1)))
am_xpol_suspect = (float(os.environ.get("AM_XPOL_BAD", -0.1)), float(os.environ.get("AM_XPOL_SUSPECT", 0)))

# bounds on solar altitude (in degrees)
good_solar_altitude = (-90, float(os.environ.get("SUSPECT_SOLAR_ALTITUDE", 0)))
suspect_solar_altitude = (float(os.environ.get("SUSPECT_SOLAR_ALTITUDE", 0)), 90)

# 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", 1.5)))
auto_rfi_suspect = (0, float(os.environ.get("AUTO_RFI_SUSPECT", 2)))

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

# bound on per-xengine non-noiselike power in diff
bad_xengine_zcut = float(os.environ.get("BAD_XENGINE_ZCUT", 10.0))

# bounds on chi^2 per antenna in omnical
oc_cspa_good = (0, float(os.environ.get("OC_CSPA_GOOD", 2)))
oc_cspa_suspect = (0, float(os.environ.get("OC_CSPA_SUSPECT", 3)))

# print bounds
for bound in ['am_corr_bad', 'am_corr_suspect', 'am_xpol_bad', 'am_xpol_suspect', 
              'good_solar_altitude', 'suspect_solar_altitude',
              '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', 'oc_cspa_good', 'oc_cspa_suspect']:
    print(f'{bound} = {eval(bound)}')

am_corr_bad = (0, 0.2)
am_corr_suspect = (0.2, 0.4)
am_xpol_bad = (-1, -0.1)
am_xpol_suspect = (-0.1, 0.0)
good_solar_altitude = (-90, 0.0)
suspect_solar_altitude = (0.0, 90)
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, 1.5)
auto_rfi_suspect = (0, 2.0)
auto_shape_good = (0, 0.1)
auto_shape_suspect = (0, 0.2)
bad_xengine_zcut = 10.0
oc_cspa_good = (0, 2.0)
oc_cspa_suspect = (0, 3.0)

read_start = time.time()
hd = io.HERADataFastReader(SUM_FILE)
data, _, _ = hd.read(read_flags=False, read_nsamples=False)
if USE_DIFF:
    hd_diff = io.HERADataFastReader(DIFF_FILE)
    diff_data, _, _ = hd_diff.read(read_flags=False, read_nsamples=False, dtype=np.complex64, fix_autos_func=np.real)
print(f'Finished loading data in {(time.time() - read_start) / 60:.2f} minutes.')

Finished loading data in 0.60 minutes.

ants = sorted(set([ant for bl in hd.bls for ant in utils.split_bl(bl)]))
auto_bls = [bl for bl in data if (bl[0] == bl[1]) and (utils.split_pol(bl[2])[0] == utils.split_pol(bl[2])[1])]
antpols = sorted(set([ant[1] for ant in ants]))

# print basic information about the file
print(f'File: {SUM_FILE}')
print(f'JDs: {hd.times} ({np.median(np.diff(hd.times)) * 24 * 3600:.5f} s integrations)')
print(f'LSTS: {hd.lsts * 12 / np.pi } hours')
print(f'Frequencies: {len(hd.freqs)} {np.median(np.diff(hd.freqs)) / 1e6:.5f} MHz channels from {hd.freqs[0] / 1e6:.5f} to {hd.freqs[-1] / 1e6:.5f} MHz')
print(f'Antennas: {len(hd.data_ants)}')
print(f'Polarizations: {hd.pols}')

File: /mnt/sn1/data1/2460743/zen.2460743.45941.sum.uvh5
JDs: [2460743.45935588 2460743.45946773] (9.66368 s integrations)
LSTS: [11.57965659 11.58234829] hours
Frequencies: 1536 0.12207 MHz channels from 46.92078 to 234.29871 MHz
Antennas: 264
Polarizations: ['nn', 'ee', 'ne', 'en']

if not USE_DIFF:
    def read_antenna_classification(df, category):
        ac = ant_class.AntennaClassification()
        ac._data = {}
        for antname, class_data, antclass in zip(df['Antenna'], df[category], df[f'{category} Class']):
            try:        
                class_data = float(class_data)
            except:
                pass
            if isinstance(class_data, str) or np.isfinite(class_data):
                ant = (int(antname[:-1]), utils._comply_antpol(antname[-1]))
                ac[ant] = antclass
                ac._data[ant] = class_data
        return ac

    df = pd.read_csv(RTP_ANTCLASS)
    am_totally_dead = read_antenna_classification(df, 'Dead?')
    am_corr = read_antenna_classification(df, 'Low Correlation')
    am_xpol = read_antenna_classification(df, 'Cross-Polarized')
    zeros_class = read_antenna_classification(df, 'Even/Odd Zeros')
    xengine_diff_class = read_antenna_classification(df, 'Bad Diff X-Engines')

if USE_DIFF:
    am = ant_metrics.AntennaMetrics(SUM_FILE, DIFF_FILE, sum_data=data, diff_data=diff_data)
    am.iterative_antenna_metrics_and_flagging(crossCut=am_xpol_bad[1], deadCut=am_corr_bad[1])
    am.all_metrics = {}  # this saves time and disk by getting rid of per-iteration information we never use
    if SAVE_RESULTS:
        am.save_antenna_metrics(AM_FILE, overwrite=True)

if USE_DIFF:
    # Turn ant metrics into classifications
    totally_dead_ants = [ant for ant, i in am.xants.items() if i == -1]
    am_totally_dead = ant_class.AntennaClassification(good=[ant for ant in ants if ant not in totally_dead_ants], bad=totally_dead_ants)
    am_corr = ant_class.antenna_bounds_checker(am.final_metrics['corr'], bad=[am_corr_bad], suspect=[am_corr_suspect], good=[(0, 1)])
    am_xpol = ant_class.antenna_bounds_checker(am.final_metrics['corrXPol'], bad=[am_xpol_bad], suspect=[am_xpol_suspect], good=[(-1, 1)])
ant_metrics_class = am_totally_dead + am_corr + am_xpol
if np.all([ant_metrics_class[utils.split_bl(bl)[0]] == 'bad' for bl in auto_bls]):
    raise ValueError('All antennas are flagged for ant_metrics.')

min_sun_alt = np.min(utils.get_sun_alt(hd.times))
solar_class = ant_class.antenna_bounds_checker({ant: min_sun_alt for ant in ants}, good=[good_solar_altitude], suspect=[suspect_solar_altitude])

if USE_DIFF:
    zeros_class = ant_class.even_odd_zeros_checker(data, diff_data, good=good_zeros_per_eo_spectrum, suspect=suspect_zeros_per_eo_spectrum)
if np.all([zeros_class[utils.split_bl(bl)[0]] == 'bad' for bl in auto_bls]):
    raise ValueError('All antennas are flagged for too many even/odd zeros.')

auto_power_class = ant_class.auto_power_checker(data, good=auto_power_good, suspect=auto_power_suspect)
auto_slope_class = ant_class.auto_slope_checker(data, good=auto_slope_good, suspect=auto_slope_suspect, edge_cut=100, filt_size=17)
if np.all([(auto_power_class + auto_slope_class)[utils.split_bl(bl)[0]] == 'bad' for bl in auto_bls]):
    raise ValueError('All antennas are flagged for bad autocorrelation power/slope.')
overall_class = auto_power_class + auto_slope_class + zeros_class + ant_metrics_class + solar_class

antenna_flags, array_flags = xrfi.flag_autos(data, flag_method="channel_diff_flagger", nsig=RFI_NSIG * 5, 
                                             antenna_class=overall_class, flag_broadcast_thresh=.5)
for key in antenna_flags:
    antenna_flags[key] = array_flags
cache = {}
_, array_flags = xrfi.flag_autos(data, freqs=data.freqs, flag_method="dpss_flagger",
                                 nsig=RFI_NSIG, antenna_class=overall_class,
                                 filter_centers=[0], filter_half_widths=[RFI_DPSS_HALFWIDTH],
                                 eigenval_cutoff=[1e-9], flags=antenna_flags, mode='dpss_matrix', 
                                 cache=cache, flag_broadcast_thresh=.5)

if USE_DIFF:
    xengine_diff_class = ant_class.non_noiselike_diff_by_xengine_checker(data, diff_data, flag_waterfall=array_flags, 
                                                                         antenna_class=overall_class, 
                                                                         xengine_chans=96, bad_xengine_zcut=bad_xengine_zcut)
    
    if np.all([overall_class[utils.split_bl(bl)[0]] == 'bad' for bl in auto_bls]):
        raise ValueError('All antennas are flagged after flagging non-noiselike diffs.')
overall_class += xengine_diff_class

def auto_bl_zscores(data, flag_array, cache={}):
    '''This function computes z-score arrays for each delay-filtered autocorrelation, normalized by the expected noise. 
    Flagged times/channels for the whole array are given 0 weight in filtering and are np.nan in the z-score.'''
    zscores = {}
    for bl in auto_bls:
        wgts = np.array(np.logical_not(flag_array), dtype=np.float64)
        model, _, _ = dspec.fourier_filter(hd.freqs, data[bl], wgts, filter_centers=[0], filter_half_widths=[RFI_DPSS_HALFWIDTH], mode='dpss_solve',
                                            suppression_factors=[1e-9], eigenval_cutoff=[1e-9], cache=cache)
        res = data[bl] - model
        int_time = 24 * 3600 * np.median(np.diff(data.times))
        chan_res = np.median(np.diff(data.freqs))
        int_count = int(int_time * chan_res)
        sigma = np.abs(model) / np.sqrt(int_count / 2)
        zscores[bl] = res / sigma    
        zscores[bl][flag_array] = np.nan

    return zscores

def rfi_from_avg_autos(data, auto_bls_to_use, prior_flags=None, nsig=RFI_NSIG):
    '''Average together all baselines in auto_bls_to_use, then find an RFI mask by looking for outliers after DPSS filtering.'''
    
    # If there are no good autos, return 100% flagged
    if len(auto_bls_to_use) == 0:
        return np.ones(data[next(iter(data))].shape, dtype=bool)
    
    # Compute int_count for all unflagged autocorrelations averaged together
    int_time = 24 * 3600 * np.median(np.diff(data.times_by_bl[auto_bls[0][0:2]]))
    chan_res = np.median(np.diff(data.freqs))
    int_count = int(int_time * chan_res) * len(auto_bls_to_use)
    avg_auto = {(-1, -1, 'ee'): np.mean([data[bl] for bl in auto_bls_to_use], axis=0)}
    
    # Flag RFI first with channel differences and then with DPSS
    antenna_flags, _ = xrfi.flag_autos(avg_auto, int_count=int_count, nsig=(nsig * 5))
    if prior_flags is not None:
        antenna_flags[(-1, -1, 'ee')] = prior_flags
    _, rfi_flags = xrfi.flag_autos(avg_auto, int_count=int_count, flag_method='dpss_flagger',
                                   flags=antenna_flags, freqs=data.freqs, filter_centers=[0],
                                   filter_half_widths=[RFI_DPSS_HALFWIDTH], eigenval_cutoff=[1e-9], nsig=nsig)

    return rfi_flags

# Iteratively develop RFI mask, excess RFI classification, and autocorrelation shape classification
stage = 1
rfi_flags = np.array(array_flags)
prior_end_states = set()
while True:
    # compute DPSS-filtered z-scores with current array-wide RFI mask
    zscores = auto_bl_zscores(data, rfi_flags)
    rms = {bl: np.nanmean(zscores[bl]**2)**.5 if np.any(np.isfinite(zscores[bl])) else np.inf for bl in zscores}
    
    # figure out which autos to use for finding new set of flags
    candidate_autos = [bl for bl in auto_bls if overall_class[utils.split_bl(bl)[0]] != 'bad']
    if stage == 1:
        # use best half of the unflagged antennas
        med_rms = np.nanmedian([rms[bl] for bl in candidate_autos])
        autos_to_use = [bl for bl in candidate_autos if rms[bl] <= med_rms]
    elif stage == 2:
        # use all unflagged antennas which are auto RFI good, or the best half, whichever is larger
        med_rms = np.nanmedian([rms[bl] for bl in candidate_autos])
        best_half_autos = [bl for bl in candidate_autos if rms[bl] <= med_rms]
        good_autos = [bl for bl in candidate_autos if (overall_class[utils.split_bl(bl)[0]] != 'bad')
                      and (auto_rfi_class[utils.split_bl(bl)[0]] == 'good')]
        autos_to_use = (best_half_autos if len(best_half_autos) > len(good_autos) else good_autos)
    elif stage == 3:
        # use all unflagged antennas which are auto RFI good or suspect
        autos_to_use = [bl for bl in candidate_autos if (overall_class[utils.split_bl(bl)[0]] != 'bad')]

    # compute new RFI flags
    rfi_flags = rfi_from_avg_autos(data, autos_to_use)

    # perform auto shape and RFI classification
    overall_class = auto_power_class + auto_slope_class + zeros_class + ant_metrics_class + solar_class + xengine_diff_class
    auto_rfi_class = ant_class.antenna_bounds_checker(rms, good=auto_rfi_good, suspect=auto_rfi_suspect, bad=(0, np.inf))
    overall_class += auto_rfi_class
    auto_shape_class = ant_class.auto_shape_checker(data, good=auto_shape_good, suspect=auto_shape_suspect,
                                                    flag_spectrum=np.sum(rfi_flags, axis=0).astype(bool), 
                                                    antenna_class=overall_class)
    overall_class += auto_shape_class
    
    # check for convergence by seeing whether we've previously gotten to this number of flagged antennas and channels
    if stage == 3:
        if (len(overall_class.bad_ants), np.sum(rfi_flags)) in prior_end_states:
            break
        prior_end_states.add((len(overall_class.bad_ants), np.sum(rfi_flags)))
    else:
        stage += 1

auto_class = auto_power_class + auto_slope_class + auto_rfi_class + auto_shape_class
if np.all([overall_class[utils.split_bl(bl)[0]] == 'bad' for bl in auto_bls]):
    raise ValueError('All antennas are flagged after flagging for bad autos power/slope/rfi/shape.')

def rfi_plot(cls, flags=rfi_flags):
    avg_auto = {(-1, -1, 'ee'): np.mean([data[bl] for bl in auto_bls if not cls[utils.split_bl(bl)[0]] == 'bad'], axis=0)}
    plt.figure(figsize=(12, 5), dpi=100)
    plt.semilogy(hd.freqs / 1e6, np.where(flags, np.nan, avg_auto[(-1, -1, 'ee')])[0], label = 'Average Good or Suspect Autocorrelation', zorder=100)
    plt.semilogy(hd.freqs / 1e6, np.where(False, np.nan, avg_auto[(-1, -1, 'ee')])[0], 'r', lw=.5, label=f'{np.sum(flags[0])} Channels Flagged for RFI')
    plt.legend()
    plt.xlabel('Frequency (MHz)')
    plt.ylabel('Uncalibrated Autocorrelation')
    plt.tight_layout()

if PLOT: rfi_plot(overall_class)

def autocorr_plot(cls):    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5), dpi=100, sharey=True, gridspec_kw={'wspace': 0})
    labels = []
    colors = ['darkgreen', 'goldenrod', 'maroon']
    for ax, pol in zip(axes, antpols):
        for ant in cls.ants:
            if ant[1] == pol:
                color = colors[cls.quality_classes.index(cls[ant])]
                ax.semilogy(np.mean(data[utils.join_bl(ant, ant)], axis=0), color=color, lw=.5)
        ax.set_xlabel('Channel', fontsize=12)
        ax.set_title(f'{utils.join_pol(pol, pol)}-Polarized Autos')

    axes[0].set_ylabel('Raw Autocorrelation', fontsize=12)
    axes[1].legend([matplotlib.lines.Line2D([0], [0], color=color) for color in colors], 
                   [cl.capitalize() for cl in cls.quality_classes], ncol=1, fontsize=12, loc='upper right', framealpha=1)
    plt.tight_layout()

if PLOT: autocorr_plot(auto_class)

def array_class_plot(cls, extra_label=""):
    outriggers = [ant for ant in hd.data_ants if ant >= 320]

    if len(outriggers) > 0:
        fig, axes = plt.subplots(1, 2, figsize=(14, 6), dpi=100, gridspec_kw={'width_ratios': [2, 1]})
        plot_antclass(hd.antpos, cls, ax=axes[0], ants=[ant for ant in hd.data_ants if ant < 320], legend=False, title=f'HERA Core{extra_label}')
        plot_antclass(hd.antpos, cls, ax=axes[1], ants=outriggers, radius=50, title='Outriggers')
    else:
        fig, axes = plt.subplots(1, 1, figsize=(9, 6), dpi=100)
        plot_antclass(hd.antpos, cls, ax=axes, ants=[ant for ant in hd.data_ants if ant < 320], legend=False, title=f'HERA Core{extra_label}')

if PLOT: array_class_plot(overall_class)

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

# delete diffs to save memory
if USE_DIFF:
    del diff_data, hd_diff
del cache
malloc_trim()

def classify_off_grid(reds, all_ants):
    '''Returns AntennaClassification of all_ants where good ants are in reds while bad ants are not.'''
    ants_in_reds = set([ant for red in reds for bl in red for ant in utils.split_bl(bl)])
    on_grid = [ant for ant in all_ants if ant in ants_in_reds]
    off_grid = [ant for ant in all_ants if ant not in ants_in_reds]
    return ant_class.AntennaClassification(good=on_grid, bad=off_grid)

def per_pol_filter_reds(reds, pols=['nn', 'ee'], **kwargs):
    '''Performs redcal filtering separately on polarizations (which might have different min_dim_size issues).'''
    return [red for pol in pols for red in redcal.filter_reds(copy.deepcopy(reds), pols=[pol], **kwargs)]

def check_if_whole_pol_flagged(redcal_class, pols=['Jee', 'Jnn']):
    '''Checks if an entire polarization is flagged. If it is, returns True and marks all antennas as bad in redcal_class.'''
    if np.logical_or(*[np.all([redcal_class[ant] == 'bad' for ant in redcal_class.ants if ant[1] == pol]) for pol in pols]):
        print('An entire polarization has been flagged. Stopping redcal.')
        for ant in redcal_class:
            redcal_class[ant] = 'bad'
        return True
    return False

def recheck_chisq(cspa, sol, cutoff, avg_alg):
    '''Recompute chisq per ant without apparently bad antennas to see if any antennas get better.'''
    avg_cspa = {ant: avg_alg(np.where(rfi_flags, np.nan, cspa[ant])) for ant in cspa}
    sol2 = redcal.RedSol(sol.reds, gains={ant: sol[ant] for ant in avg_cspa if avg_cspa[ant] <= cutoff}, vis=sol.vis)
    new_chisq_per_ant = {ant: np.array(cspa[ant]) for ant in sol2.gains}
    if len(set([bl[2] for red in per_pol_filter_reds(sol2.reds, ants=sol2.gains.keys(), antpos=hd.data_antpos, **fr_settings) for bl in red])) >= 2:
        redcal.expand_omni_gains(sol2, sol2.reds, data, chisq_per_ant=new_chisq_per_ant)
    for ant in avg_cspa:
        if ant in new_chisq_per_ant:
            if np.any(np.isfinite(new_chisq_per_ant[ant])):
                if not np.all(np.isclose(new_chisq_per_ant[ant], 0)):
                    new_avg_cspa = avg_alg(np.where(rfi_flags, np.nan, cspa[ant]))
                    if new_avg_cspa > 0:
                        avg_cspa[ant] = np.min([avg_cspa[ant], new_avg_cspa])
    return avg_cspa

# figure out and filter reds and classify antennas based on whether or not they are on the main grid
fr_settings = {'max_dims': OC_MAX_DIMS, 'min_dim_size': OC_MIN_DIM_SIZE, 'min_bl_cut': OC_MIN_BL_LEN, 'max_bl_cut': OC_MAX_BL_LEN}
reds = redcal.get_reds(hd.data_antpos, pols=['ee', 'nn'], pol_mode='2pol', bl_error_tol=2.0)
reds = per_pol_filter_reds(reds, ex_ants=overall_class.bad_ants, antpos=hd.data_antpos, **fr_settings)
if OC_SKIP_OUTRIGGERS:
    reds = redcal.filter_reds(reds, ex_ants=[ant for ant in ants if ant[0] >= 320])
redcal_class = classify_off_grid(reds, ants)

if OC_USE_PRIOR_SOL:
    # Find closest omnical file
    omnical_files = sorted(glob.glob('.'.join(OMNICAL_FILE.split('.')[:-5]) + '.*.' + '.'.join(OMNICAL_FILE.split('.')[-3:])))
    if len(omnical_files) == 0:
        OC_USE_PRIOR_SOL = False
    else:
        omnical_jds = np.array([float(re.findall("\d+\.\d+", ocf)[-1]) for ocf in omnical_files])
        closest_omnical = omnical_files[np.argmin(np.abs(omnical_jds - data.times[0]))]

        # Load closest omnical file and use it if the antenna flagging is not too dissimilar
        hc = io.HERACal(closest_omnical)
        prior_gains, prior_flags, _, _ = hc.read()
        not_bad_not_prior_flagged = [ant for ant in overall_class if not ant in redcal_class.bad_ants and not np.all(prior_flags[ant])]
        if (len(redcal_class.bad_ants) == len(redcal_class.ants)):
            OC_USE_PRIOR_SOL = False  # all antennas flagged
        elif (len(not_bad_not_prior_flagged) / (len(redcal_class.ants) - len(redcal_class.bad_ants))) < OC_PRIOR_SOL_FLAG_THRESH:
            OC_USE_PRIOR_SOL = False  # too many antennas unflaged that were flagged in the prior sol
        else:
            print(f'Using {closest_omnical} as a starting point for redcal.')

redcal_start = time.time()
rc_settings = {'max_dims': OC_MAX_DIMS, 'oc_conv_crit': 1e-10, 'gain': 0.4, 'run_logcal': False,
               'oc_maxiter': OC_MAXITER, 'check_after': OC_MAXITER, 'use_gpu': OC_USE_GPU}

if check_if_whole_pol_flagged(redcal_class):
    # skip redcal, initialize empty sol and meta 
    sol = redcal.RedSol(reds)
    meta = {'chisq': None, 'chisq_per_ant': None}
else:    
    if OC_USE_PRIOR_SOL:
        # use prior unflagged gains and data to create starting point for next step
        ants_in_reds = set([ant for red in reds for bl in red for ant in utils.split_bl(bl)])
        sol = redcal.RedSol(reds=reds, gains={ant: prior_gains[ant] for ant in not_bad_not_prior_flagged})
        reds_to_update = [[bl for bl in red if (utils.split_bl(bl)[0] in sol.gains) and (utils.split_bl(bl)[1] in sol.gains)] for red in reds]
        reds_to_update = [red for red in reds_to_update if len(red) > 0]
        sol.update_vis_from_data(data, reds_to_update=reds_to_update)
        redcal.expand_omni_gains(sol, reds, data)
        sol.update_vis_from_data(data)
    else:
        # perform first stage of redundant calibration 
        meta, sol = redcal.redundantly_calibrate(data, reds, **rc_settings)
        max_dly = np.max(np.abs(list(meta['fc_meta']['dlys'].values())))  # Needed for RFI delay-slope cal
        median_cspa = recheck_chisq(meta['chisq_per_ant'], sol, oc_cspa_suspect[1] * 5, np.nanmedian)
         # remove particularly bad antennas (5x the bound on median, not mean)
        cspa_class = ant_class.antenna_bounds_checker(median_cspa, good=(oc_cspa_good[0], oc_cspa_suspect[1] * 5), bad=[(-np.inf, np.inf)])
        redcal_class += cspa_class
        print(f'Removing {cspa_class.bad_ants} for >5x high median chi^2.')
        for ant in cspa_class.bad_ants:
            print(f'\t{ant}: {median_cspa[ant]:.3f}')
        
    malloc_trim()

Removing {(np.int32(224), 'Jee'), (np.int32(146), 'Jnn'), (np.int32(105), 'Jee'), (np.int32(85), 'Jee'), (np.int32(194), 'Jnn'), (np.int32(151), 'Jee'), (np.int32(62), 'Jnn'), (np.int32(73), 'Jnn'), (np.int32(95), 'Jnn'), (np.int32(174), 'Jnn'), (np.int32(31), 'Jnn'), (np.int32(228), 'Jee'), (np.int32(65), 'Jee'), (np.int32(261), 'Jee'), (np.int32(282), 'Jnn'), (np.int32(231), 'Jnn'), (np.int32(205), 'Jnn'), (np.int32(217), 'Jee'), (np.int32(87), 'Jee'), (np.int32(272), 'Jee'), (np.int32(242), 'Jnn'), (np.int32(66), 'Jnn'), (np.int32(243), 'Jee'), (np.int32(189), 'Jnn'), (np.int32(222), 'Jnn'), (np.int32(35), 'Jnn'), (np.int32(102), 'Jee'), (np.int32(90), 'Jnn'), (np.int32(212), 'Jee'), (np.int32(245), 'Jee'), (np.int32(136), 'Jnn'), (np.int32(178), 'Jnn'), (np.int32(255), 'Jnn'), (np.int32(104), 'Jee'), (np.int32(59), 'Jnn'), (np.int32(62), 'Jee'), (np.int32(125), 'Jnn'), (np.int32(139), 'Jee'), (np.int32(17), 'Jnn'), (np.int32(61), 'Jnn'), (np.int32(226), 'Jnn'), (np.int32(172), 'Jee'), (np.int32(238), 'Jee'), (np.int32(270), 'Jnn'), (np.int32(207), 'Jee'), (np.int32(77), 'Jee'), (np.int32(110), 'Jee'), (np.int32(231), 'Jee'), (np.int32(205), 'Jee'), (np.int32(63), 'Jnn'), (np.int32(198), 'Jee'), (np.int32(242), 'Jee'), (np.int32(188), 'Jnn'), (np.int32(45), 'Jnn'), (np.int32(78), 'Jnn'), (np.int32(124), 'Jnn'), (np.int32(189), 'Jee'), (np.int32(3), 'Jnn'), (np.int32(211), 'Jee'), (np.int32(126), 'Jnn'), (np.int32(148), 'Jnn'), (np.int32(202), 'Jee'), (np.int32(125), 'Jee')} for >5x high median chi^2.
	(np.int32(224), 'Jee'): 32.945
	(np.int32(146), 'Jnn'): 126.421
	(np.int32(105), 'Jee'): 23.175
	(np.int32(85), 'Jee'): 119.461
	(np.int32(194), 'Jnn'): 68.040
	(np.int32(151), 'Jee'): 32.368
	(np.int32(62), 'Jnn'): 139.447
	(np.int32(73), 'Jnn'): 108.661
	(np.int32(95), 'Jnn'): 63.253
	(np.int32(174), 'Jnn'): 80.038
	(np.int32(31), 'Jnn'): 110.060
	(np.int32(228), 'Jee'): 47.592
	(np.int32(65), 'Jee'): 119.693
	(np.int32(261), 'Jee'): 97.766
	(np.int32(282), 'Jnn'): 35.522
	(np.int32(231), 'Jnn'): 108.455
	(np.int32(205), 'Jnn'): 77.968
	(np.int32(217), 'Jee'): 237.451
	(np.int32(87), 'Jee'): 39.336
	(np.int32(272), 'Jee'): 78.615
	(np.int32(242), 'Jnn'): 79.487
	(np.int32(66), 'Jnn'): 81.151
	(np.int32(243), 'Jee'): 37.530
	(np.int32(189), 'Jnn'): 102.528
	(np.int32(222), 'Jnn'): 72.934
	(np.int32(35), 'Jnn'): 88.349
	(np.int32(102), 'Jee'): 70.402
	(np.int32(90), 'Jnn'): 101.608
	(np.int32(212), 'Jee'): 71.890
	(np.int32(245), 'Jee'): 108.173
	(np.int32(136), 'Jnn'): 73.836
	(np.int32(178), 'Jnn'): 45.401
	(np.int32(255), 'Jnn'): 62.767
	(np.int32(104), 'Jee'): 46.178
	(np.int32(59), 'Jnn'): 90.964
	(np.int32(62), 'Jee'): 43.710
	(np.int32(125), 'Jnn'): 135.802
	(np.int32(139), 'Jee'): 33.842
	(np.int32(17), 'Jnn'): 97.454
	(np.int32(61), 'Jnn'): 96.502
	(np.int32(226), 'Jnn'): 60.012
	(np.int32(172), 'Jee'): 55.456
	(np.int32(238), 'Jee'): 225.457
	(np.int32(270), 'Jnn'): 94.804
	(np.int32(207), 'Jee'): 90.707
	(np.int32(77), 'Jee'): 44.182
	(np.int32(110), 'Jee'): 100.641
	(np.int32(231), 'Jee'): 56.830
	(np.int32(205), 'Jee'): 34.632
	(np.int32(63), 'Jnn'): 89.794
	(np.int32(198), 'Jee'): 261.650
	(np.int32(242), 'Jee'): 55.595
	(np.int32(188), 'Jnn'): 126.116
	(np.int32(45), 'Jnn'): 71.406
	(np.int32(78), 'Jnn'): 82.620
	(np.int32(124), 'Jnn'): 110.708
	(np.int32(189), 'Jee'): 114.262
	(np.int32(3), 'Jnn'): 101.721
	(np.int32(211), 'Jee'): 60.127
	(np.int32(126), 'Jnn'): 133.384
	(np.int32(148), 'Jnn'): 93.250
	(np.int32(202), 'Jee'): 91.919
	(np.int32(125), 'Jee'): 58.007

# iteratively rerun redundant calibration
redcal_done = False
rc_settings['oc_maxiter'] = rc_settings['check_after'] = OC_RERUN_MAXITER
for i in range(OC_MAX_RERUN + 1):
    # refilter reds and update classification to reflect new off-grid ants, if any
    reds = per_pol_filter_reds(reds, ex_ants=(overall_class + redcal_class).bad_ants, antpos=hd.data_antpos, **fr_settings)
    reds = sorted(reds, key=len, reverse=True)
    redcal_class += classify_off_grid(reds, ants)
    ants_in_reds = set([ant for red in reds for bl in red for ant in utils.split_bl(bl)])
    
    # check to see whether we're done
    if check_if_whole_pol_flagged(redcal_class) or redcal_done or (i == OC_MAX_RERUN):
        break

    # re-run redundant calibration using previous solution, updating bad and suspicious antennas
    meta, sol = redcal.redundantly_calibrate(data, reds, sol0=sol, **rc_settings)
    malloc_trim()
    
    # recompute chi^2 for bad antennas without bad antennas to make sure they are actually bad
    mean_cspa = recheck_chisq(meta['chisq_per_ant'], sol, oc_cspa_suspect[1], np.nanmean)
    
    # remove bad antennas
    cspa_class = ant_class.antenna_bounds_checker(mean_cspa, good=oc_cspa_good, suspect=oc_cspa_suspect, bad=[(-np.inf, np.inf)])
    for ant in cspa_class.bad_ants:
        if mean_cspa[ant] < np.max(list(mean_cspa.values())) / OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE:
            cspa_class[ant] = 'suspect'  # reclassify as suspect if they are much better than the worst antennas
    redcal_class += cspa_class
    print(f'Removing {cspa_class.bad_ants} for high mean unflagged chi^2.')
    for ant in cspa_class.bad_ants:
        print(f'\t{ant}: {mean_cspa[ant]:.3f}')

    if len(cspa_class.bad_ants) == 0:
        redcal_done = True  # no new antennas to flag

print(f'Finished redcal in {(time.time() - redcal_start) / 60:.2f} minutes.')

An entire polarization has been flagged. Stopping redcal.
Finished redcal in 0.10 minutes.

overall_class += redcal_class

expanded_reds = redcal.get_reds(hd.data_antpos, pols=['ee', 'nn'], pol_mode='2pol', bl_error_tol=2.0)
expanded_reds = per_pol_filter_reds(expanded_reds, ex_ants=(ant_metrics_class + solar_class + zeros_class + auto_class + xengine_diff_class).bad_ants,
                                    max_dims=OC_MAX_DIMS, min_dim_size=OC_MIN_DIM_SIZE)
if OC_SKIP_OUTRIGGERS:
    expanded_reds = redcal.filter_reds(expanded_reds, ex_ants=[ant for ant in ants if ant[0] >= 320])
if len(sol.gains) > 0:
    redcal.expand_omni_vis(sol, expanded_reds, data, chisq=meta['chisq'], chisq_per_ant=meta['chisq_per_ant'])

# now figure out flags, nsamples etc.
omni_flags = {ant: (~np.isfinite(g)) | (ant in overall_class.bad_ants) for ant, g in sol.gains.items()}
vissol_flags = datacontainer.RedDataContainer({bl: ~np.isfinite(v) for bl, v in sol.vis.items()}, reds=sol.vis.reds)
single_nsamples_array = np.ones((len(hd.times), len(hd.freqs)), dtype=float)
nsamples = datacontainer.DataContainer({bl: single_nsamples_array for bl in data})
vissol_nsamples = redcal.count_redundant_nsamples(nsamples, [red for red in expanded_reds if red[0] in vissol_flags], 
                                                  good_ants=[ant for ant in overall_class if ant not in overall_class.bad_ants])
for bl in vissol_flags:
    vissol_flags[bl][vissol_nsamples[bl] == 0] = True
sol.make_sol_finite()

if not OC_USE_PRIOR_SOL:
    # find channels clearly contaminated by RFI
    not_bad_ants = [ant for ant in overall_class.ants if (overall_class[ant] != 'bad') and (utils.join_bl(ant, ant) in data)]
    if len(not_bad_ants) > 0:
        chan_flags = np.mean([xrfi.detrend_medfilt(data[utils.join_bl(ant, ant)], Kf=8, Kt=2) for ant in not_bad_ants], axis=(0, 1)) > 5

        # hardcoded RFI transmitters and their headings
        # channel: frequency (Hz), heading (rad), chi^2
        phs_sol = {359: ( 90744018.5546875, 0.7853981, 23.3),
                   360: ( 90866088.8671875, 0.7853981, 10.8),
                   385: ( 93917846.6796875, 0.7853981, 27.3),
                   386: ( 94039916.9921875, 0.7853981, 18.1),
                   400: ( 95748901.3671875, 6.0632738, 24.0),
                   441: (100753784.1796875, 0.7853981, 21.7),
                   442: (100875854.4921875, 0.7853981, 19.4),
                   455: (102462768.5546875, 6.0632738, 18.8),
                   456: (102584838.8671875, 6.0632738,  8.8),
                   471: (104415893.5546875, 0.7853981, 13.3),
                   484: (106002807.6171875, 6.0632738, 21.2),
                   485: (106124877.9296875, 6.0632738,  4.0),
                  1181: (191085815.4296875, 0.7853981, 26.3),
                  1182: (191207885.7421875, 0.7853981, 27.0),
                  1183: (191329956.0546875, 0.7853981, 25.6),
                  1448: (223678588.8671875, 2.6075219, 25.7),
                  1449: (223800659.1796875, 2.6075219, 22.6),
                  1450: (223922729.4921875, 2.6075219, 11.6),
                  1451: (224044799.8046875, 2.6075219,  5.9),
                  1452: (224166870.1171875, 2.6075219, 22.6),
                  1510: (231246948.2421875, 0.1068141, 23.9)}

        if not np.isclose(hd.freqs[0], 46920776.3671875, atol=0.001) or len(hd.freqs) != 1536:
            # We have less frequencies than usual (maybe testing)
            phs_sol = {np.argmin(np.abs(hd.freqs - freq)): (freq, heading, chisq) for chan, (freq, heading, chisq) in phs_sol.items() if hd.freqs[0] <= freq <= hd.freqs[-1]}


        rfi_chans = [chan for chan in phs_sol if chan_flags[chan]]
        if len(rfi_chans) >= 2:
            print('Channels used for delay-slope calibration with RFI:', rfi_chans)
            rfi_angles = np.array([phs_sol[chan][1] for chan in rfi_chans])
            rfi_headings = np.array([np.cos(rfi_angles), np.sin(rfi_angles), np.zeros_like(rfi_angles)])
            rfi_chisqs = np.array([phs_sol[chan][2] for chan in rfi_chans])

            # resolve firstcal degeneracy with delay slopes set by RFI transmitters, update cal
            RFI_dly_slope_gains = abscal.RFI_delay_slope_cal([red for red in expanded_reds if red[0] in sol.vis], hd.antpos, sol.vis, hd.freqs, rfi_chans, rfi_headings, rfi_wgts=rfi_chisqs**-1,
                                                             min_tau=-max_dly, max_tau=max_dly, delta_tau=0.1e-9, return_gains=True, gain_ants=sol.gains.keys())
            sol.gains = {ant: g * RFI_dly_slope_gains[ant] for ant, g in sol.gains.items()}
            apply_cal.calibrate_in_place(sol.vis, RFI_dly_slope_gains)
            malloc_trim()
        else:
            print(f"Only {len(rfi_chans)} RFI channels with known headings were flagged for RFI, so RFI-firstcal is being skipped.")

# Load simulated and then downsampled model of autocorrelations that includes receiver noise, then interpolate to upsample
if VALIDATION:
    hd_auto_model = io.HERAData('/lustre/aoc/projects/hera/Validation/H6C_IDR2/sim_data/h6c_validation_autos_for_amp_abscal_with_Trx_100K.uvh5')
else:
    hd_auto_model = io.HERAData(f'{HNBT_DATA}/SSM_autocorrelations_downsampled_sum_pol_convention.uvh5')
model, _, _ = hd_auto_model.read()
per_pol_interpolated_model = {}
for bl in model:
    sorted_lsts, lst_indices = np.unique(model.lsts, return_index=True)
    periodic_model = np.vstack([model[bl][lst_indices, :], model[bl][lst_indices[0], :]])
    periodic_lsts = np.append(sorted_lsts, sorted_lsts[0] + 2 * np.pi)
    lst_interpolated = interpolate.CubicSpline(periodic_lsts, periodic_model, axis=0, bc_type='periodic')(data.lsts)
    per_pol_interpolated_model[bl[2]] = interpolate.CubicSpline(model.freqs, lst_interpolated, axis=1)(data.freqs)
model = {bl: per_pol_interpolated_model[bl[2]] for bl in auto_bls if utils.split_bl(bl)[0] not in overall_class.bad_ants}

# Run abscal and update omnical gains with abscal gains
if len(model) > 0:
    redcaled_autos = {bl: sol.calibrate_bl(bl, data[bl]) for bl in auto_bls if utils.split_bl(bl)[0] not in overall_class.bad_ants}
    g_abscal = abscal.abs_amp_logcal(model, redcaled_autos, verbose=False, return_gains=True, gain_ants=sol.gains)
    sol.gains = {ant: g * g_abscal[ant] for ant, g in sol.gains.items()}
    apply_cal.calibrate_in_place(sol.vis, g_abscal)
    del redcaled_autos, g_abscal

if ABSCAL_MODEL_FILES_GLOB is not None:
    abscal_model_files = sorted(glob.glob(ABSCAL_MODEL_FILES_GLOB))
elif VALIDATION:
    abscal_model_files = sorted(glob.glob('/lustre/aoc/projects/hera/Validation/H6C_IDR2/sim_data/foregrounds/zen.LST.*.foregrounds.uvh5'))
else:
    # try to find files on site
    abscal_model_files = sorted(glob.glob('/mnt/sn1/data1/abscal_models/H6C/zen.2458894.?????.uvh5'))
    if len(abscal_model_files) == 0:
        # try to find files at NRAO
        abscal_model_files = sorted(glob.glob('/lustre/aoc/projects/hera/h6c-analysis/abscal_models/h6c_abscal_files_unique_baselines/zen.2458894.?????.uvh5'))
print(f'Found {len(abscal_model_files)} abscal model files{" in " + os.path.dirname(abscal_model_files[0]) if len(abscal_model_files) > 0 else ""}.')

Found 425 abscal model files in /mnt/sn1/data1/abscal_models/H6C.

# Try to perform a full abscal of phase
if len(abscal_model_files) == 0:
    DO_FULL_ABSCAL = False
    print('No model files found... not performing full absolute calibration of phase gradients.')
elif np.all([ant in overall_class.bad_ants for ant in ants]):
    DO_FULL_ABSCAL = False
    print('All antennas classified as bad... skipping absolute calibration of phase gradients.')
else:
    abscal_start = time.time()
    # figure out which model files match the LSTs of the data
    matched_model_files = sorted(set(abscal.match_times(SUM_FILE, abscal_model_files, filetype='uvh5')))
    if len(matched_model_files) == 0:
        DO_FULL_ABSCAL = False
        print(f'No model files found matching the LSTs of this file after searching for {(time.time() - abscal_start) / 60:.2f} minutes. '
              'Not performing full absolute calibration of phase gradients.')
    else:
        DO_FULL_ABSCAL = True
        # figure out appropriate model times to load
        hdm = io.HERAData(matched_model_files)
        all_model_times, all_model_lsts = abscal.get_all_times_and_lsts(hdm, unwrap=True)
        d2m_time_map = abscal.get_d2m_time_map(data.times, np.unwrap(data.lsts), all_model_times, all_model_lsts, extrap_limit=.5)

All antennas classified as bad... skipping absolute calibration of phase gradients.

if DO_FULL_ABSCAL:
    abscal_meta = {}
    for pol in ['ee', 'nn']:
        print(f'Performing absolute phase gradient calibration of {pol}-polarized visibility solutions...')
        
        # load matching times and baselines
        unflagged_data_bls = [bl for bl in vissol_flags if not np.all(vissol_flags[bl]) and bl[2] == pol]
        model_bls = copy.deepcopy(hdm.bls)
        model_antpos = hdm.data_antpos
        if len(matched_model_files) > 1:  # in this case, it's a dictionary
            model_bls = list(set([bl for bls in list(hdm.bls.values()) for bl in bls]))
            model_antpos = {ant: pos for antpos in hdm.data_antpos.values() for ant, pos in antpos.items()}
        data_bls, model_bls, data_to_model_bl_map = abscal.match_baselines(unflagged_data_bls, model_bls, data.antpos, model_antpos=model_antpos, 
                                                                         pols=[pol], data_is_redsol=True, model_is_redundant=True, tol=1.0,
                                                                         min_bl_cut=ABSCAL_MIN_BL_LEN, max_bl_cut=ABSCAL_MAX_BL_LEN, verbose=True)
        model, model_flags, _ = io.partial_time_io(hdm, np.unique([d2m_time_map[time] for time in data.times]), bls=model_bls)
        model_bls = [data_to_model_bl_map[bl] for bl in data_bls]
        
        # rephase model to match in lsts
        model_blvecs = {bl: model.antpos[bl[0]] - model.antpos[bl[1]] for bl in model.keys()}
        utils.lst_rephase(model, model_blvecs, model.freqs, data.lsts - model.lsts,
                          lat=hdm.telescope_location_lat_lon_alt_degrees[0], inplace=True)

        # run abscal and apply 
        abscal_meta[pol], delta_gains = abscal.complex_phase_abscal(sol.vis, model, sol.reds, data_bls, model_bls)
        
        # apply gains
        sol.gains = {antpol : g * delta_gains.get(antpol, 1) for antpol, g in sol.gains.items()}
        apply_cal.calibrate_in_place(sol.vis, delta_gains)            
     
    del model, model_flags, delta_gains
    malloc_trim()    
    
    print(f'Finished absolute calibration of tip-tilt phase slopes in {(time.time() - abscal_start) / 60:.2f} minutes.')

if DO_FULL_ABSCAL and CALIBRATE_CROSS_POLS:
    cross_pol_cal_start = time.time()

    # Compute reds for good antennas 
    cross_reds = redcal.get_reds(hd.data_antpos, pols=['en', 'ne'], bl_error_tol=2.0)        
    cross_reds = redcal.filter_reds(cross_reds, ex_ants=overall_class.bad_ants, pols=['en', 'ne'], antpos=hd.antpos, **fr_settings)    
    unflagged_data_bls = [red[0] for red in cross_reds]

    # Get cross-polarized model visibilities
    model_bls = copy.deepcopy(hdm.bls)
    model_antpos = hdm.data_antpos
    if len(matched_model_files) > 1:  # in this case, it's a dictionary
        model_bls = list(set([bl for bls in list(hdm.bls.values()) for bl in bls]))
        model_antpos = {ant: pos for antpos in hdm.data_antpos.values() for ant, pos in antpos.items()}

    data_bls, model_bls, data_to_model_bl_map = abscal.match_baselines(unflagged_data_bls, model_bls, data.antpos, model_antpos=model_antpos, 
                                                                     pols=['en', 'ne'], data_is_redsol=False, model_is_redundant=True, tol=1.0,
                                                                     min_bl_cut=ABSCAL_MIN_BL_LEN, max_bl_cut=ABSCAL_MAX_BL_LEN, verbose=True)
    
    model, model_flags, _ = io.partial_time_io(hdm, np.unique([d2m_time_map[time] for time in data.times]), 
                                               bls=list(set([bl[0:2] for bl in model_bls])), polarizations=['en', 'ne'])
    model_bls = [data_to_model_bl_map[bl] for bl in data_bls]

    # rephase model to match in lsts
    model_blvecs = {bl: model.antpos[bl[0]] - model.antpos[bl[1]] for bl in model.keys()}
    utils.lst_rephase(model, model_blvecs, model.freqs, data.lsts - model.lsts, lat=hdm.telescope_location_lat_lon_alt_degrees[0], inplace=True)


    wgts_here = {}
    data_here = {}

    
    for red in cross_reds:
        data_bl = red[0]
        if data_bl in data_to_model_bl_map:

            wgts_here[data_bl] = np.sum([
                np.logical_not(omni_flags[utils.split_bl(bl)[0]] | omni_flags[utils.split_bl(bl)[1]])
                for bl in red
            ], axis=0)
            data_here[data_bl] = np.nanmean([
                np.where(
                    omni_flags[utils.split_bl(bl)[0]] | omni_flags[utils.split_bl(bl)[1]],
                    np.nan, sol.calibrate_bl(bl, data[bl])
                ) 
                for bl in red
            ], axis=0)
    
    # Run cross-polarized phase calibration
    delta = abscal.cross_pol_phase_cal(
        model=model, data=data_here, wgts=wgts_here, data_bls=data_bls, model_bls=model_bls, return_gains=False, 
        refpol='Jee', gain_ants=sol.gains.keys()
    )
    delta_gains = {antpol: (np.ones_like(delta) if antpol[1] == 'Jee' else np.exp(1j * delta)) for antpol in sol.gains.keys()}
    
    # apply gains
    # \Delta = \phi_e - \phi_n, where V_{en}^{cal} = V_{en}^{uncal} * e^{i \Delta} 
    # and V_{ne}^{cal} = V_{ne}^{uncal} * e^{-i \Delta}
    sol.gains = {antpol: g * delta_gains[antpol] for antpol, g in sol.gains.items()}
    apply_cal.calibrate_in_place(sol.vis, delta_gains)
    del hdm, model, model_flags, delta_gains
    print(f'Finished relative polarized phase calibration in {(time.time() - cross_pol_cal_start) / 60:.2f} minutes.')

def redundant_group_plot():
    if np.all([ant in overall_class.bad_ants for ant in ants]):
        print('All antennas classified as bad. Nothing to plot.')
        return
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 6), dpi=100, sharex='col', sharey='row', gridspec_kw={'hspace': 0, 'wspace': 0})
    for i, pol in enumerate(['ee', 'nn']):
        reds_here = redcal.get_reds(hd.data_antpos, pols=[pol], pol_mode='1pol', bl_error_tol=2.0)
        red = sorted(redcal.filter_reds(reds_here, ex_ants=overall_class.bad_ants), key=len, reverse=True)[0]
        rc_data = {bl: sol.calibrate_bl(bl, data[bl]) for bl in red}
        for bl in red:
            axes[0, i].plot(hd.freqs/1e6, np.angle(rc_data[bl][0]), alpha=.5, lw=.5)
            axes[1, i].semilogy(hd.freqs/1e6, np.abs(rc_data[bl][0]), alpha=.5, lw=.5)
        axes[0, i].plot(hd.freqs / 1e6, np.angle(sol.vis[red[0]][0]), lw=1, c='k')
        axes[1, i].semilogy(hd.freqs / 1e6, np.abs(sol.vis[red[0]][0]), lw=1, c='k', label=f'Baseline Group:\n{(int(red[0][0]), int(red[0][1]), red[0][2])}')
        axes[1, i].set_xlabel('Frequency (MHz)')
        axes[1, i].legend(loc='upper right')
    axes[0, 0].set_ylabel('Visibility Phase (radians)')
    axes[1, 0].set_ylabel('Visibility Amplitude (Jy)')
    plt.tight_layout()

def abscal_degen_plot():
    if DO_FULL_ABSCAL:
        fig, axes = plt.subplots(3, 1, figsize=(14, 6), dpi=100, sharex=True, gridspec_kw={'hspace': .05})

        for ax, pol in zip(axes[:2], ['ee', 'nn']):
            for kk in range(abscal_meta[pol]['Lambda_sol'].shape[-1]):
                ax.plot(hd.freqs[~rfi_flags[0]] * 1e-6, abscal_meta[pol]['Lambda_sol'][0, ~rfi_flags[0], kk], '.', ms=1, label=f"Component {kk}")

            ax.set_ylim(-np.pi-0.5, np.pi+0.5)
            ax.set_xlabel('Frequency (MHz)')
            ax.set_ylabel('Phase Gradient\nVector Component')
            ax.legend(markerscale=20, title=f'{pol}-polarization', loc='lower right')
            ax.grid()
            
        for pol, color in zip(['ee', 'nn'], ['b', 'r']):
            axes[2].plot(hd.freqs[~rfi_flags[0]]*1e-6, abscal_meta[pol]['Z_sol'].real[0, ~rfi_flags[0]], '.', ms=1, label=pol, color=color)
        axes[2].set_ylim(-.25, 1.05)
        axes[2].set_ylabel('Re[Z($\\nu$)]')
        axes[2].legend(markerscale=20, loc='lower right')
        axes[2].grid()            
        plt.tight_layout()

def polarized_gain_phase_plot():
    if CALIBRATE_CROSS_POLS and DO_FULL_ABSCAL:
        plt.figure(figsize=(14, 4), dpi=100)
        for i, time in enumerate(data.times):
            plt.plot(data.freqs / 1e6, np.where(rfi_flags[i], np.nan, delta[i, :]), '.', ms=1.5, label=f'{time:.6f}')
        plt.ylim(-np.pi-0.5, np.pi+0.5)
        plt.xlabel('Frequency (MHz)')
        plt.ylabel('Relative Phase $\Delta \ (\phi_{ee} - \phi_{nn})$')
        plt.grid()
        plt.legend()

if PLOT: redundant_group_plot()

All antennas classified as bad. Nothing to plot.

if PLOT: abscal_degen_plot()

if PLOT: polarized_gain_phase_plot()

expand_start = time.time()
expanded_reds = redcal.get_reds(hd.data_antpos, pols=['ee', 'nn'], pol_mode='2pol', bl_error_tol=2.0)
sol.vis.build_red_keys(expanded_reds)
redcal.expand_omni_gains(sol, expanded_reds, data, chisq_per_ant=meta['chisq_per_ant'])
if not np.all([ant in overall_class.bad_ants for ant in ants]):
    redcal.expand_omni_vis(sol, expanded_reds, data)

# Replace near-zeros in gains and infs/nans in gains/sols
for ant in sol.gains:
    zeros_in_gains = np.isclose(sol.gains[ant], 0)
    if ant in omni_flags:
        omni_flags[ant][zeros_in_gains] = True
    sol.gains[ant][zeros_in_gains] = 1.0 + 0.0j
sol.make_sol_finite()
malloc_trim()
print(f'Finished expanding gain solution in {(time.time() - expand_start) / 60:.2f} minutes.')

Finished expanding gain solution in 0.78 minutes.

def array_chisq_plot(include_outriggers=True):
    if np.all([ant in overall_class.bad_ants for ant in ants]):
        print('All antennas classified as bad. Nothing to plot.')
        return    
    
    def _chisq_subplot(ants, size=250):
        fig, axes = plt.subplots(1, 2, figsize=(14, 5), dpi=100)
        for ax, pol in zip(axes, ['ee', 'nn']):
            ants_to_plot = set([ant for ant in meta['chisq_per_ant'] if utils.join_pol(ant[1], ant[1]) == pol and (ant[0] in ants)])
            cspas = np.array([np.nanmean(np.where(rfi_flags, np.nan, meta['chisq_per_ant'][ant])) for ant in ants_to_plot])
            xpos = [hd.antpos[ant[0]][0] for ant in ants_to_plot]
            ypos = [hd.antpos[ant[0]][1] for ant in ants_to_plot]
            scatter = ax.scatter(xpos, ypos, s=size, c=cspas, lw=.25, edgecolors=np.where(np.isfinite(cspas) & (cspas > 0), 'none', 'k'), 
                                 norm=matplotlib.colors.LogNorm(vmin=1, vmax=oc_cspa_suspect[1]))
            for ant in ants_to_plot:
                ax.text(hd.antpos[ant[0]][0], hd.antpos[ant[0]][1], ant[0], va='center', ha='center', fontsize=8,
                        c=('r' if ant in overall_class.bad_ants else 'w'))
            plt.colorbar(scatter, ax=ax, extend='both')
            ax.axis('equal')
            ax.set_xlabel('East-West Position (meters)')
            ax.set_ylabel('North-South Position (meters)')
            ax.set_title(f'{pol}-pol $\\chi^2$ / Antenna (Red is Flagged)')
        plt.tight_layout()    
    
    _chisq_subplot([ant for ant in hd.data_ants if ant < 320])
    outriggers = [ant for ant in hd.data_ants if ant >= 320]    
    if include_outriggers & (len(outriggers) > 0):
        _chisq_subplot([ant for ant in hd.data_ants if ant >= 320], size=400)

if PLOT: array_chisq_plot(include_outriggers=(not OC_SKIP_OUTRIGGERS))

All antennas classified as bad. Nothing to plot.

if PLOT: array_class_plot(overall_class, extra_label=", Post-Redcal")

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed y limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

WARNING:matplotlib.axes._base:Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.

to_show = {'Antenna': [f'{ant[0]}{ant[1][-1]}' for ant in ants]}
classes = {'Antenna': [overall_class[ant] if ant in overall_class else '-' for ant in ants]}
to_show['Dead?'] = [{'good': 'No', 'bad': 'Yes'}[am_totally_dead[ant]] if (ant in am_totally_dead) else '' for ant in ants]
classes['Dead?'] = [am_totally_dead[ant] if (ant in am_totally_dead) else '' for ant in ants]
for title, ac in [('Low Correlation', am_corr),
                  ('Cross-Polarized', am_xpol),
                  ('Solar Alt', solar_class),
                  ('Even/Odd Zeros', zeros_class),
                  ('Autocorr Power', auto_power_class),
                  ('Autocorr Slope', auto_slope_class),
                  ('Auto RFI RMS', auto_rfi_class),
                  ('Autocorr Shape', auto_shape_class),
                  ('Bad Diff X-Engines', xengine_diff_class)]:
    to_show[title] = [f'{ac._data[ant]:.2G}' if (ant in ac._data) else '' for ant in ants]
    classes[title] = [ac[ant] if ant in ac else 'bad' for ant in ants]
    
to_show['Redcal chi^2'] = [f'{np.nanmean(np.where(rfi_flags, np.nan, meta["chisq_per_ant"][ant])):.3G}' \
                           if (meta['chisq_per_ant'] is not None and ant in meta['chisq_per_ant']) else '' for ant in ants]
classes['Redcal chi^2'] = [redcal_class[ant] if ant in redcal_class else '' for ant in ants]

df = pd.DataFrame(to_show)
df_classes = pd.DataFrame(classes)
colors = {'good': 'darkgreen', 'suspect': 'goldenrod', 'bad': 'maroon'}
df_colors = df_classes.applymap(lambda x: f'background-color: {colors.get(x, None)}')

table = df.style.hide() \
                .apply(lambda x: pd.DataFrame(df_colors.values, columns=x.columns), axis=None) \
                .set_properties(subset=['Antenna'], **{'font-weight': 'bold', 'border-right': "3pt solid black"}) \
                .set_properties(subset=df.columns[1:], **{'border-left': "1pt solid black"}) \
                .set_properties(**{'text-align': 'center', 'color': 'white'})

HTML(table.to_html())

# Save antenna classification table as a csv
if SAVE_RESULTS:
    for ind, col in zip(np.arange(len(df.columns), 0, -1), df_classes.columns[::-1]):
        df.insert(int(ind), col + ' Class', df_classes[col])
    df.to_csv(ANTCLASS_FILE)

print('Final Ant-Pol Classification:\n\n', overall_class)

Final Ant-Pol Classification:

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


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

# update flags in omnical gains and visibility solutions
for ant in omni_flags:
    omni_flags[ant] |= rfi_flags
for bl in vissol_flags:
    vissol_flags[bl] |= rfi_flags

if SAVE_RESULTS:
    add_to_history = 'Produced by file_calibration notebook with the following environment:\n' + '=' * 65 + '\n' + os.popen('conda env export').read() + '=' * 65    
    
    hd_vissol = io.HERAData(SUM_FILE)
    hc_omni = hd_vissol.init_HERACal(gain_convention='divide', cal_style='redundant')
    hc_omni.pol_convention = hd_auto_model.pol_convention
    hc_omni.gain_scale = hd_auto_model.vis_units
    hc_omni.update(gains=sol.gains, flags=omni_flags, quals=meta['chisq_per_ant'], total_qual=meta['chisq'])
    hc_omni.history += add_to_history
    hc_omni.write_calfits(OMNICAL_FILE, clobber=True)
    del hc_omni
    malloc_trim()
    
    if SAVE_OMNIVIS_FILE:
        # output results, harmonizing keys over polarizations
        all_reds = redcal.get_reds(hd.data_antpos, pols=['ee', 'nn', 'en', 'ne'], pol_mode='4pol', bl_error_tol=2.0)
        bl_to_red_map = {bl: red[0] for red in all_reds for bl in red}
        hd_vissol.read(bls=[bl_to_red_map[bl] for bl in sol.vis], return_data=False)
        hd_vissol.empty_arrays()
        hd_vissol.history += add_to_history
        hd_vissol.update(data={bl_to_red_map[bl]: sol.vis[bl] for bl in sol.vis}, 
                         flags={bl_to_red_map[bl]: vissol_flags[bl] for bl in vissol_flags}, 
                         nsamples={bl_to_red_map[bl]: vissol_nsamples[bl] for bl in vissol_nsamples})
        hd_vissol.pol_convention = hd_auto_model.pol_convention
        hd_vissol.vis_units = hd_auto_model.vis_units
        hd_vissol.write_uvh5(OMNIVIS_FILE, clobber=True)

if SAVE_RESULTS:
    del hd_vissol
    malloc_trim()

if SAVE_RESULTS and not os.path.exists(OMNICAL_FILE):
    print(f'WARNING: No calibration file produced at {OMNICAL_FILE}. Creating a fully-flagged placeholder calibration file.')
    hd_writer = io.HERAData(SUM_FILE)
    io.write_cal(OMNICAL_FILE, freqs=hd_writer.freqs, times=hd_writer.times,
                 gains={ant: np.ones((hd_writer.Ntimes, hd_writer.Nfreqs), dtype=np.complex64) for ant in ants},
                 flags={ant: np.ones((len(data.times), len(data.freqs)), dtype=bool) for ant in ants},
                 quality=None, total_qual=None, outdir='', overwrite=True, history=utils.history_string(add_to_history), 
                 x_orientation=hd_writer.x_orientation, telescope_location=hd_writer.telescope_location, lst_array=np.unique(hd_writer.lsts),
                 antenna_positions=np.array([hd_writer.antenna_positions[hd_writer.antenna_numbers == antnum].flatten() for antnum in set(ant[0] for ant in ants)]),
                 antnums2antnames=dict(zip(hd_writer.antenna_numbers, hd_writer.antenna_names)))

if SAVE_RESULTS and SAVE_OMNIVIS_FILE and not os.path.exists(OMNIVIS_FILE):
    print(f'WARNING: No omnivis file produced at {OMNIVIS_FILE}. Creating an empty visibility solution file.')
    hd_writer = io.HERAData(SUM_FILE)
    hd_writer.initialize_uvh5_file(OMNIVIS_FILE, clobber=True)

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

pyuvdata: 3.1.3
hera_cal: 3.7.1.dev18+g10e9584
hera_filters: 0.1.6.dev9+gf165ec1
hera_qm: 2.2.1.dev2+ga535e9e
hera_notebook_templates: 0.1.dev989+gee0995d

print(f'Finished execution in {(time.time() - tstart) / 60:.2f} minutes.')

Finished execution in 3.32 minutes.

Antenna	Dead?	Low Correlation	Cross-Polarized	Solar Alt	Autocorr Power	Autocorr Slope	Auto RFI RMS	Autocorr Shape	Bad Diff X-Engines	Redcal chi^2
3e	No	0.95	0.089	-55	15	-0.074	1.1	0.072	7	609
3n	No	0.92	0.089	-55	30	-0.22	0.57	0.033	0	321
4e	No	0.86	0.19	-55	0.94	0.24	2.7	0.077		148
4n	No	0.91	0.19	-55	22	-0.15	0.65	0.084	3	2.8E+03
5e	No	0.94	0.09	-55	12	-0.23	0.95	0.062	15	315
5n	No	0.91	0.09	-55	22	-0.31	0.66	0.057	2	1.59E+03
7e	No	0.93	0.1	-55	16	-0.099	1.2	0.036	8	1.02E+03
7n	No	0.9	0.1	-55	26	-0.29	1	0.049	1	747
8e	No	0.9	0.13	-55	62	0.029	10	0.04		30.6
8n	No	0.87	0.13	-55	61	-0.074	0.34	0.033		218
9e	No	0.91	0.11	-55	20	-0.14	0.73	0.038	5	612
9n	No	0.88	0.11	-55	22	-0.26	0.66	0.049	3	548
10e	No	0.78	0.27	-55	12	-0.38	0.79	0.11	23	253
10n	No	0.3	0.27	-55	0.7	0.5	2.3	0.17		3.8
15e	No	0.95	0.083	-55	20	-0.093	0.73	0.037	2	720
15n	No	0.92	0.083	-55	11	-0.63	0.96	0.15		692
16e	No	0.95	0.087	-55	19	-0.11	0.81	0.044	3	715
16n	No	0.92	0.087	-55	19	-0.27	0.73	0.055	4	2.28E+03
17e	No	0.94	0.1	-55	14	-0.18	0.8	0.051	14	101
17n	No	0.9	0.1	-55	46	-0.12	0.46	0.012	0	1.17E+03
18e	No	0.93	0.21	-55	18	-0.12	0.71	0.044	4	663
18n	No	0.83	0.21	-55	20	-0.21	28	0.12	1	419
19e	No	0.93	0.097	-55	18	-0.15	0.8	0.045	7	63.6
19n	No	0.9	0.097	-55	21	-0.31	0.68	0.053	4	902
20e	No	0.82	-0.072	-55	8.7	-0.46	1	0.12	28	789
20n	No	0.78	-0.072	-55	2.7	-0.29	1.7	0.075	32	271
21e	No	0.91	0.12	-55	30	-0.047	10	0.025	0	691
21n	No	0.88	0.12	-55	29	-0.19	0.68	0.036	1	400
22e	No	0.89	0.12	-55	29	-0.031	0.65	0.015	2	223
22n	No	0.87	0.12	-55	35	-0.21	0.53	0.032	1	387
27e	No	0.73	0.3	-55	8.9	0.4	57	0.37	7	30.2
27n	No	0.84	0.3	-55	28	-0.085	20	0.1	0	71.9
28e	No	0.69	0.087	-55	10	0.094	40	0.19	7	116
28n	No	0.87	0.087	-55	13	-0.2	1.1	0.045	12	259
29e	No	0.92	0.066	-55	4.2	0.39	1.7	0.17	32	432
29n	No	0.91	0.066	-55	31	0.38	1.9	0.15	1	1.59E+03
30e	No	0.94	0.087	-55	18	-0.16	0.78	0.047	6	951
30n	No	0.92	0.087	-55	22	-0.3	0.8	0.05	3	1.08E+03
31e	No	0.93	0.11	-55	13	-0.29	0.87	0.079	16	28
31n	No	0.9	0.11	-55	56	-0.086	0.41	0.026	0	390

Single File Calibration¶

Parse inputs and outputs¶

Parse settings¶

Parse bounds¶

Load sum and diff data¶

Classify good, suspect, and bad antpols¶

Load classifications that use diffs if diffs are not available¶

Run ant_metrics¶

Mark sun-up (or high solar altitude) data as suspect¶

Classify antennas responsible for 0s in visibilities as bad:¶

Examine and classify autocorrelation power and slope¶

Find starting set of array flags¶

Classify antennas based on non-noiselike diffs¶

Examine and classify autocorrelation excess RFI and shape, finding consensus RFI mask along the way¶

Figure 1: RFI Flagging¶

Figure 2: Plot of autocorrelations with classifications¶

Summarize antenna classification prior to redundant-baseline calibration¶

Figure 3: Summary of antenna classifications prior to calibration¶

Perform redundant-baseline calibration¶

Perform iterative redcal¶

Expand solution to include calibratable baselines excluded from redcal (e.g. because they were too long)¶

Fix the firstcal delay slope degeneracy using RFI transmitters¶

Perform absolute amplitude calibration using a model of autocorrelations¶

Full absolute calibration of phase gradients¶

Plotting¶

Figure 4: Redundant calibration of a single baseline group¶

Figure 5: Absolute calibration of redcal degeneracies¶

Figure 6: Relative Phase Calibration¶

Attempt to calibrate some flagged antennas¶

Figure 7: chi^2 per antenna across the array¶

Figure 8: Summary of antenna classifications after redundant calibration¶

Table 1: Complete summary of per-antenna classifications¶

Save calibration solutions¶

Output fully flagged calibration file if OMNICAL_FILE is not written¶

Output empty visibility file if OMNIVIS_FILE is not written¶

Metadata¶

Run `ant_metrics`¶

Perform iterative `redcal`¶

Fix the `firstcal` delay slope degeneracy using RFI transmitters¶

Figure 5: Absolute calibration of `redcal` degeneracies¶

Output fully flagged calibration file if `OMNICAL_FILE` is not written¶

Output empty visibility file if `OMNIVIS_FILE` is not written¶