Single File Calibration¶

by Josh Dillon, Aaron Parsons, Tyler Cox, and Zachary Martinot, last updated August 11, 2025

This notebook is designed to infer as much information about the array from a single file, including pushing the calibration and RFI mitigation as far as possible. Calibration includes redundant-baseline calibration, RFI-based calibration of delay slopes, model-based calibration of overall amplitudes, and a full per-frequency phase gradient absolute calibration if abscal model files are available.

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

• Figure 1: RFI Flagging¶

• Figure 2: Plot of autocorrelations with classifications¶

• Figure 3: Summary of antenna classifications prior to calibration¶

• Figure 4: Redundant calibration of a single baseline group¶

• Figure 5: Absolute calibration of redcal degeneracies¶

• Figure 6: Relative Phase Calibration¶

• 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¶

In [1]:
import time
tstart = time.time()
!hostname
still3.rtp.pvt
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
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'
In [3]:
# this enables better memory management on linux
import ctypes
def malloc_trim():
    try:
        ctypes.CDLL('libc.so.6').malloc_trim(0) 
    except OSError:
        pass

Parse inputs and outputs¶

To use this notebook interactively, you will have to provide a sum filename path if none exists as an environment variable. All other parameters have reasonable default values.

In [4]:
# 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/2461003/zen.2461003.43882.sum.uvh5'
DIFF_FILE = '/mnt/sn1/data1/2461003/zen.2461003.43882.diff.uvh5'
AM_FILE = '/mnt/sn1/data1/2461003/zen.2461003.43882.sum.ant_metrics.hdf5'
ANTCLASS_FILE = '/mnt/sn1/data1/2461003/zen.2461003.43882.sum.ant_class.csv'
OMNICAL_FILE = '/mnt/sn1/data1/2461003/zen.2461003.43882.sum.omni.calfits'
OMNIVIS_FILE = '/mnt/sn1/data1/2461003/zen.2461003.43882.sum.omni_vis.uvh5'
SAVE_RESULTS = 'True'
SAVE_OMNIVIS_FILE = 'False'
USE_DIFF = 'True'
VALIDATION = 'False'

Parse settings¶

Load settings relating to the operation of the notebook, then print what was loaded (or default).

In [5]:
# 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", "TRUE").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", 10))
OC_RERUN_MAXITER = int(os.environ.get("OC_MAXITER", 50))
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', '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
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

Parse bounds¶

Load settings related to classifying antennas as good, suspect, or bad, then print what was loaded (or default).

In [6]:
# ant_metrics bounds for low correlation / dead antennas
am_corr_bad = (0, float(os.environ.get("AM_CORR_BAD", 0.35)))
am_corr_suspect = (float(os.environ.get("AM_CORR_BAD", 0.35)), float(os.environ.get("AM_CORR_SUSPECT", 0.45)))

# 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.35)
am_corr_suspect = (0.35, 0.45)
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)

Load sum and diff data¶

In [7]:
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 1.15 minutes.
In [8]:
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]))
In [9]:
# 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/2461003/zen.2461003.43882.sum.uvh5
JDs: [2461003.43876498 2461003.43887683] (9.66368 s integrations)
LSTS: [4.16872184 4.17141355] hours
Frequencies: 1536 0.12207 MHz channels from 46.92078 to 234.29871 MHz
Antennas: 290
Polarizations: ['nn', 'ee', 'ne', 'en']

Classify good, suspect, and bad antpols¶

In [10]:
ALL_FLAGGED = False
def all_flagged():
    if ALL_FLAGGED:
        print('All antennas are flagged, so this cell is being skipped.')
    return ALL_FLAGGED

# initialize classes to None to help make Table 1 when everything is flagged
overall_class = None
am_totally_dead = None
am_corr = None
am_xpol = None
solar_class = None
zeros_class = None
auto_power_class = None
auto_slope_class = None
auto_rfi_class = None
auto_shape_class = None
xengine_diff_class = None
meta = None
redcal_class = None 

Load classifications that use diffs if diffs are not available¶

In [11]:
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')

Run ant_metrics¶

This classifies antennas as cross-polarized, low-correlation, or dead. Such antennas are excluded from any calibration.

In [12]:
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)
In [13]:
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]):
    ALL_FLAGGED = True
    print('All antennas are flagged for ant_metrics.')

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

In [14]:
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])

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

This classifier looks for X-engine failure or packet loss specific to an antenna which causes either the even visibilities (or the odd ones, or both) to be 0s.

In [15]:
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]):
    ALL_FLAGGED = True
    print('All antennas are flagged for too many even/odd zeros.')

Examine and classify autocorrelation power and slope¶

These classifiers look for antennas with too high or low power or to steep a slope.

In [16]:
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]):
    ALL_FLAGGED = True
    print('All antennas are flagged for bad autocorrelation power/slope.')
overall_class = auto_power_class + auto_slope_class + zeros_class + ant_metrics_class + solar_class

Find starting set of array flags¶

In [17]:
if not all_flagged():
    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)

Classify antennas based on non-noiselike diffs¶

In [18]:
if not all_flagged():
    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]):
            ALL_FLAGGED = True
            print('All antennas are flagged after flagging non-noiselike diffs.')
    overall_class += xengine_diff_class

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

This classifier iteratively identifies antennas for excess RFI (characterized by RMS of DPSS-filtered autocorrelations after RFI flagging) and bad shape, as determined by a discrepancy with the mean good normalized autocorrelation's shape. Along the way, it iteratively discovers a conensus array-wide RFI mask.

In [19]:
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
In [20]:
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
In [21]:
# Iteratively develop RFI mask, excess RFI classification, and autocorrelation shape classification
if not all_flagged():
    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
In [22]:
auto_class = auto_power_class + auto_slope_class
if auto_rfi_class is not None:
    auto_class += auto_rfi_class
if auto_shape_class is not None:
    auto_class += auto_rfi_class
if np.all([overall_class[utils.split_bl(bl)[0]] == 'bad' for bl in auto_bls]):
    ALL_FLAGGED = True
    print('All antennas are flagged after flagging for bad autos power/slope/rfi/shape.')
All antennas are flagged after flagging for bad autos power/slope/rfi/shape.
In [23]:
if not all_flagged():
    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()
All antennas are flagged, so this cell is being skipped.

Figure 1: RFI Flagging¶

This figure shows RFI identified using the average of all autocorrelations---excluding bad antennas---for the first integration in the file.

In [24]:
if PLOT and not all_flagged(): rfi_plot(overall_class)
All antennas are flagged, so this cell is being skipped.
In [25]:
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()

Figure 2: Plot of autocorrelations with classifications¶

This figure shows a plot of all autocorrelations in the array, split by polarization. Antennas are classified based on their autocorrelations into good, suspect, and bad, by examining power, slope, and RFI-occupancy.

In [26]:
if PLOT and not all_flagged(): autocorr_plot(auto_class)
All antennas are flagged, so this cell is being skipped.

Summarize antenna classification prior to redundant-baseline calibration¶

In [27]:
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}')

Figure 3: Summary of antenna classifications prior to calibration¶

This figure shows the location and classification of all antennas prior to calibration. Antennas are split along the diagonal, with ee-polarized antpols represented by the southeast half of each antenna and nn-polarized antpols represented by the northwest half. Outriggers are split from the core and shown at exaggerated size in the right-hand panel. This classification includes ant_metrics, a count of the zeros in the even or odd visibilities, and autocorrelation power, slope, and RFI occupancy. An antenna classified as bad in any classification will be considered bad. An antenna marked as suspect any in any classification will be considered suspect unless it is also classified as bad elsewhere.

In [28]:
if PLOT and not all_flagged(): array_class_plot(overall_class)
All antennas are flagged, so this cell is being skipped.
In [29]:
# delete diffs to save memory
if USE_DIFF:
    del diff_data, hd_diff
try:
    del cache
except NameError:
    pass
malloc_trim()

Perform redundant-baseline calibration¶

In [30]:
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)
In [31]:
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)]
In [32]:
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
In [33]:
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

Perform iterative redcal¶

In [34]:
# figure out and filter reds and classify antennas based on whether or not they are on the main grid
if not all_flagged():
    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)
All antennas are flagged, so this cell is being skipped.
In [35]:
if OC_USE_PRIOR_SOL and not all_flagged():
    # 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.')
All antennas are flagged, so this cell is being skipped.
In [36]:
if not all_flagged():
    redcal_start = time.time()
    rc_settings = {'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
            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, max_dims=None, **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()
All antennas are flagged, so this cell is being skipped.
In [37]:
if not all_flagged():
    # 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)
        
        # 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, max_dims=None, **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.')
    overall_class += redcal_class
All antennas are flagged, so this cell is being skipped.

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

In [38]:
if not all_flagged():
    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'])
All antennas are flagged, so this cell is being skipped.
In [39]:
if not all_flagged():
    # 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()
All antennas are flagged, so this cell is being skipped.

Fix the firstcal delay slope degeneracy using RFI transmitters¶

In [40]:
if not OC_USE_PRIOR_SOL and not all_flagged():
    # 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.")

Perform absolute amplitude calibration using a model of autocorrelations¶

In [41]:
# 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')
if not all_flagged():
    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}
All antennas are flagged, so this cell is being skipped.
In [42]:
if not all_flagged():
    # 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
All antennas are flagged, so this cell is being skipped.

Full absolute calibration of phase gradients¶

If an ABSCAL_MODEL_FILES_GLOB is provided, try to perform a full absolute calibration of tip-tilt phase gradients across the array using that those model files. Specifically, this step calibrates omnical visbility solutions using unique baselines simulated with a model of the sky and HERA's beam.

In [43]:
if not all_flagged():
    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 ""}.')
All antennas are flagged, so this cell is being skipped.
In [44]:
if not all_flagged():
    # 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 are flagged, so this cell is being skipped.
In [45]:
if not all_flagged():
    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.deg, 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.')
All antennas are flagged, so this cell is being skipped.
In [46]:
if not all_flagged() and 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.deg, 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.')
All antennas are flagged, so this cell is being skipped.

Plotting¶

In [47]:
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()
In [48]:
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()
In [49]:
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()

Figure 4: Redundant calibration of a single baseline group¶

The results of a redundant-baseline calibration of a single integration and a single group, the one with the highest redundancy in each polarization after antenna classification and excision based on the above, plus the removal of antennas with high chi^2 per antenna. The black line is the redundant visibility solution. Each thin colored line is a different baseline group. Phases are shown in the top row, amplitudes in the bottom, ee-polarized visibilities in the left column, and nn-polarized visibilities in the right.

In [50]:
if PLOT and not all_flagged(): redundant_group_plot()
All antennas are flagged, so this cell is being skipped.

Figure 5: Absolute calibration of redcal degeneracies¶

This figure shows the per-frequency phase gradient solutions across the array for both polarizations and all components of the degenerate subspace of redundant-baseline calibraton. While full HERA only has two such tip-tilt degeneracies, a subset of HERA can have up to OC_MAX_DIMS (depending on antenna flagging). In addition to the absolute amplitude, this is the full set of the calibration degrees of freedom not constrainted by redcal. This figure also includes a plot of $Re[Z(\nu)]$, the complex objective function which varies from -1 to 1 and indicates how well the data and the absolute calibration model have been made to agree. Perfect agreement is 1.0 and good agreement is anything above $\sim$0.5 Decorrelation yields values closer to 0, where anything below $\sim$0.3 is suspect.

In [51]:
if PLOT and not all_flagged(): abscal_degen_plot()
All antennas are flagged, so this cell is being skipped.

Figure 6: Relative Phase Calibration¶

This figure shows the relative phase calibration between the ee vs. nn polarizations.

In [52]:
if PLOT and not all_flagged(): polarized_gain_phase_plot()
All antennas are flagged, so this cell is being skipped.

Attempt to calibrate some flagged antennas¶

This attempts to calibrate bad antennas using information from good or suspect antennas without allowing bad antennas to affect their calibration. However, introducing 0s in gains or infs/nans in gains or visibilities can create problems down the line, so those are removed.

In [53]:
if not all_flagged():
    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.')
All antennas are flagged, so this cell is being skipped.
In [54]:
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)

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

This plot shows median (taken over time and frequency) of the normalized chi^2 per antenna. The expectation value for this quantity when the array is perfectly redundant is 1.0. Antennas that are classified as bad for any reason have their numbers shown in red. Some of those antennas were classified as bad during redundant calibration for high chi^2. Some of those antennas were originally excluded from redundant calibration because they were classified as bad earlier for some reason. See here for more details. Note that the color scale saturates at below 1 and above 10.

In [55]:
if PLOT and not all_flagged(): array_chisq_plot(include_outriggers=(not OC_SKIP_OUTRIGGERS))
All antennas are flagged, so this cell is being skipped.

Figure 8: Summary of antenna classifications after redundant calibration¶

This figure is the same as Figure 2, except that it now includes additional suspect or bad antennas based on redundant calibration. This can include antennas with high chi^2, but it can also include antennas classified as "bad" because they would add extra degeneracies to calibration.

In [56]:
if PLOT and not all_flagged(): array_class_plot(overall_class, extra_label=", Post-Redcal")
All antennas are flagged, so this cell is being skipped.
In [57]:
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 (ac is not None and ant in ac._data) else '' for ant in ants]
    classes[title] = [ac[ant] if (ac is not None and 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 is not None and 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 redcal_class is not None and ant in redcal_class else 'bad' 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'})

Table 1: Complete summary of per-antenna classifications¶

This table summarizes the results of the various classifications schemes detailed above. As before, green is good, yellow is suspect, and red is bad. The color for each antenna (first column) is the final summary of all other classifications. Antennas missing from redcal $\chi^2$ were excluded redundant-baseline calibration, either because they were flagged by ant_metrics or the even/odd zeros check, or because they would add unwanted extra degeneracies.

In [58]:
HTML(table.to_html())
Out[58]:
Antenna Dead? Low Correlation Cross-Polarized Solar Alt Even/Odd Zeros Autocorr Power Autocorr Slope Auto RFI RMS Autocorr Shape Bad Diff X-Engines Redcal chi^2
3e No 0.023 0.0016 -39 0 0.69 0.52 INF INF
3n No 0.023 0.0016 -39 0 0.62 0.6 INF INF
4e No 0.54 0.45 -39 0 24 0.096 INF INF 15
4n No 0.03 0.45 -39 0 0.62 0.59 INF INF
5e No 0.52 0.35 -39 0 15 0.12 INF INF 15
5n No 0.52 0.35 -39 0 19 0.2 INF INF 13
7e No 0.28 0.19 -39 0 20 0.15 INF INF
7n No 0.27 0.19 -39 0 19 0.17 INF INF
8e No 0.26 0.17 -39 0 54 0.06 INF INF
8n No 0.25 0.17 -39 0 10 0.26 INF INF
9e No 0.27 0.18 -39 0 17 0.18 INF INF
9n No 0.25 0.18 -39 0 18 0.28 INF INF
10e No 0.36 0.26 -39 0 23 0.14 INF INF 29
10n No 0.35 0.26 -39 0 16 0.29 INF INF 21
15e No 0.57 0.37 -39 0 29 0.11 INF INF 9
15n No 0.56 0.37 -39 0 30 0.11 INF INF 8
16e No 0.57 0.38 -39 0 16 0.16 INF INF 10
16n No 0.56 0.38 -39 0 16 0.3 INF INF 7
17e No 0.57 0.38 -39 0 16 0.13 INF INF 9
17n No 0.55 0.38 -39 0 19 0.24 INF INF 8
18e No 0.5 0.38 -39 0 9.5 0.12 INF INF 19
18n No 0.34 0.38 -39 0 13 0.55 INF INF
19e No 0.38 0.27 -39 0 19 0.15 INF INF 29
19n No 0.37 0.27 -39 0 17 0.2 INF INF 24
20e No 0.38 0.27 -39 0 18 0.16 INF INF 28
20n No 0.37 0.27 -39 0 19 0.2 INF INF 25
21e No 0.33 0.23 -39 0 7.8 0.29 INF INF
21n No 0.32 0.23 -39 0 10 0.35 INF INF
27e No 0.11 -0.014 -39 0 9.9 0.64 INF INF
27n No 0.064 -0.014 -39 0 7.6 1 INF INF
28e No 0.024 0.0027 -39 0 0.72 0.52 INF INF
28n No 0.021 0.0027 -39 1 0.67 0.6 INF INF
29e No 0.5 0.28 -39 0 20 0.92 INF INF
29n No 0.41 0.28 -39 0 14 1.4 INF INF
30e No 0.57 0.39 -39 0 14 0.18 INF INF 7
30n No 0.53 0.39 -39 0 2.3 0.28 INF INF 6
31e No 0.37 0.25 -39 0 16 0.15 INF INF 22
31n No 0.35 0.25 -39 0 15 0.23 INF INF
32e No 0.35 0.23 -39 0 15 0.19 INF INF 16
32n No 0.22 0.23 -39 0 15 1.1 INF INF
33e No 0.41 0.31 -39 0 17 0.21 INF INF 25
33n No 0.25 0.31 -39 0 19 0.35 INF INF
36e No 0.56 0.32 -39 0 23 -0.18 INF INF 11
36n No 0.56 0.32 -39 0 18 -0.1 INF INF 11
37e No 0.61 0.38 -39 0 11 0.17 INF INF 5
37n No 0.61 0.38 -39 0 23 0.19 INF INF 5
38e No 0.61 0.38 -39 0 19 0.16 INF INF 4
38n No 0.61 0.38 -39 0 22 0.18 INF INF 5
40e No 0.56 0.37 -39 0 8.8 0.17 INF INF 15
40n No 0.52 0.37 -39 0 2.8 0.28 INF INF 11
41e No 0.41 0.27 -39 0 18 0.066 INF INF 29
41n No 0.39 0.27 -39 0 19 0.19 INF INF 24
42e No 0.46 0.32 -39 0 19 0.16 INF INF 23
42n No 0.4 0.32 -39 1 1.6 0.27 INF INF 16
43e No 0.5 0.33 -39 0 18 0.13 INF INF 19
43n No 0.48 0.33 -39 0 15 0.3 INF INF 18
44e No 0.024 0.0014 -39 0 0.67 0.52 INF INF
44n No 0.022 0.0014 -39 1 0.61 0.58 INF INF
45e No 0.49 0.34 -39 0 17 0.2 INF INF 24
45n No 0.47 0.34 -39 0 18 0.31 INF INF 22
46e No 0.46 0.33 -39 0 15 1.5 INF INF
46n No 0.46 0.33 -39 0 13 0.24 INF INF 24
50e No 0.54 0.32 -39 0 15 0.16 INF INF 11
50n No 0.54 0.32 -39 0 16 0.13 INF INF 7
51e Yes -39 1.5E+03 0 0 INF INF
51n Yes -39 1.5E+03 0 0 INF INF
52e No 0.55 0.32 -39 0 19 -0.12 INF INF 11
52n No 0.55 0.32 -39 0 16 -0.027 INF INF 12
53e Yes -39 1.5E+03 0 0 INF INF
53n Yes -39 1.5E+03 0 0 INF INF
54e No 0.46 0.29 -39 0 16 0.17 INF INF 23
54n No 0.45 0.29 -39 0 16 0.24 INF INF 23
55e No 0.47 0.28 -39 0 17 0.15 INF INF 15
55n No 0.47 0.28 -39 0 20 0.18 INF INF 16
56e No 0.58 0.36 -39 0 45 0.094 INF INF 15
56n No 0.58 0.36 -39 0 38 0.17 INF INF 12
57e No 0.4 0.28 -39 0 13 0.14 INF INF 26
57n No 0.38 0.28 -39 0 12 0.27 INF INF 26
58e No 0.55 0.37 -39 0 14 0.25 INF INF 18
58n No 0.53 0.37 -39 0 14 0.35 INF INF 16
59e No 0.55 0.38 -39 0 23 0.2 INF INF 23
59n No 0.54 0.38 -39 0 18 0.3 INF INF 19
60e No 0.48 0.33 -39 0 16 0.32 INF INF 18
60n No 0.45 0.33 -39 0 4.9 0.36 INF INF 13
65e No 0.62 0.39 -39 0 17 0.18 INF INF 3
65n No 0.63 0.39 -39 0 16 0.2 INF INF 3
66e No 0.62 0.38 -39 0 22 0.15 INF INF 7
66n No 0.62 0.38 -39 0 20 0.19 INF INF 6
67e No 0.64 0.39 -39 0 21 0.14 INF INF 3
67n No 0.64 0.39 -39 0 25 0.19 INF INF 2
68e No 0.64 0.37 -39 0 14 0.14 INF INF 4
68n No 0.64 0.37 -39 0 12 0.14 INF INF 4
69e No 0.58 0.35 -39 0 15 0.16 INF INF 15
69n No 0.57 0.35 -39 0 20 0.24 INF INF 13
70e No 0.48 0.29 -39 0 26 0.14 INF INF 17
70n No 0.46 0.29 -39 0 53 0.14 INF INF 16
71e No 0.48 0.31 -39 0 15 0.16 INF INF 20
71n No 0.42 0.31 -39 0 1.8 0.36 INF INF 15
72e No 0.46 0.29 -39 0 16 0.29 INF INF 23
72n No 0.44 0.29 -39 0 14 0.25 INF INF 20
73e No 0.51 0.35 -39 0 18 0.24 INF INF 26
73n No 0.49 0.35 -39 0 13 0.24 INF INF 24
74e No 0.5 0.33 -39 0 20 0.19 INF INF 18
74n No 0.49 0.33 -39 0 15 0.2 INF INF 18
75e No 0.45 0.32 -39 0 14 0.3 INF INF 27
75n No 0.44 0.32 -39 0 14 0.32 INF INF 25
76e No 0.023 0.0029 -39 0 3.2 0.56 INF INF
76n No 0.019 0.0029 -39 0 2.9 0.57 INF INF
79e No 0.23 0.15 -39 0 18 0.3 INF INF
79n No 0.23 0.15 -39 0 24 0.27 INF INF
80e No 0.3 0.21 -39 32 19 0.28 INF INF
80n No 0.27 0.21 -39 32 13 0.41 INF INF
81e No 0.59 0.38 -39 1 0.71 -0.17 INF INF
81n No 0.64 0.38 -39 0 2.5 -0.2 INF INF 3
82e No 0.65 0.39 -39 0 18 0.3 INF INF 2
82n No 0.64 0.39 -39 0 17 0.28 INF INF 2
83e No 0.63 0.38 -39 0 5.1 0.24 INF INF 2
83n No 0.61 0.38 -39 0 17 0.6 INF INF
84e No 0.63 0.37 -39 0 18 0.35 INF INF 2
84n No 0.64 0.37 -39 0 17 0.29 INF INF 2
85e No 0.52 0.3 -39 0 13 0.091 INF INF 13
85n No 0.51 0.3 -39 0 19 0.28 INF INF 10
86e No 0.028 0.0013 -39 0 1.2 0.54 INF INF
86n No 0.027 0.0013 -39 0 1.1 0.59 INF INF
87e No 0.42 0.27 -39 0 15 0.98 INF INF
87n No 0.5 0.27 -39 0 20 0.17 INF INF 11
88e Yes -39 1.5E+03 0 0 INF INF
88n Yes -39 1.5E+03 0 0 INF INF
89e No 0.53 0.33 -39 0 6 0.16 INF INF 14
89n No 0.53 0.33 -39 0 15 0.19 INF INF 22
90e Yes -39 1.5E+03 0 0 INF INF
90n Yes -39 1.5E+03 0 0 INF INF
91e No 0.36 0.24 -39 0 16 0.13 INF INF 26
91n No 0.33 0.24 -39 0 10 0.24 INF INF
92e No 0.29 0.19 -39 0 14 0.41 INF INF
92n No 0.29 0.19 -39 0 14 0.26 INF INF
93e No 0.24 0.2 -39 0 2.1 0.24 INF INF
93n No 0.28 0.2 -39 0 14 0.24 INF INF
94e No 0.31 0.21 -39 0 23 0.19 INF INF
94n No 0.3 0.21 -39 0 28 0.25 INF INF
95e No 0.24 0.16 -39 0 18 0.3 INF INF
95n No 0.23 0.16 -39 0 25 0.36 INF INF
96e No 0.28 0.2 -39 0 12 0.37 INF INF
96n No 0.26 0.2 -39 0 11 0.51 INF INF
97e No 0.26 0.19 -39 0 17 0.2 INF INF
97n No 0.17 0.19 -39 0 5.2 0.52 INF INF
98e No 0.63 0.36 -39 0 9.2 0.14 INF INF 2
98n No 0.63 0.36 -39 0 11 0.14 INF INF 4
99e No 0.026 0.37 -39 1 0.071 0.44 INF INF
99n No 0.52 0.37 -39 0 0.32 -0.11 INF INF
100e No 0.61 0.36 -39 0 23 0.2 INF INF 4
100n No 0.61 0.36 -39 0 18 0.29 INF INF 2
101e No 0.51 0.28 -39 0 16 0.12 INF INF 13
101n No 0.5 0.28 -39 0 15 0.17 INF INF 9
102e No 0.52 0.29 -39 0 23 0.18 INF INF 12
102n No 0.51 0.29 -39 0 21 0.3 INF INF 12
103e Yes -39 1.5E+03 0 0 INF INF
103n Yes -39 1.5E+03 0 0 INF INF
104e Yes -39 1.5E+03 0 0 INF INF
104n Yes -39 1.5E+03 0 0 INF INF
105e No 0.37 0.24 -39 0 14 0.19 INF INF 24
105n No 0.35 0.24 -39 0 14 0.19 INF INF 18
106e No 0.55 0.33 -39 0 16 0.15 INF INF 25
106n No 0.54 0.33 -39 0 15 0.17 INF INF 11
107e Yes -39 1.5E+03 0 0 INF INF
107n Yes -39 1.5E+03 0 0 INF INF
108e No 0.37 0.24 -39 0 15 0.28 INF INF 27
108n No 0.34 0.24 -39 0 11 0.37 INF INF
109e No 0.35 0.26 -39 0 18 0.18 INF INF
109n No 0.021 0.26 -39 1 0.67 0.58 INF INF
110e No 0.33 0.23 -39 0 15 0.29 INF INF
110n No 0.32 0.23 -39 0 14 0.31 INF INF
111e No 0.28 0.21 -39 0 19 0.84 INF INF
111n No 0.31 0.21 -39 0 16 0.21 INF INF
112e No 0.3 0.21 -39 0 22 0.19 INF INF
112n No 0.28 0.21 -39 0 14 0.21 INF INF
113e No 0.24 0.16 -39 0 21 0.27 INF INF
113n No 0.21 0.16 -39 0 20 0.56 INF INF
114e No 0.3 0.22 -39 0 20 0.29 INF INF
114n No 0.29 0.22 -39 0 19 0.31 INF INF
115e No 0.24 0.16 -39 0 16 0.2 INF INF
115n No 0.22 0.16 -39 0 12 0.37 INF INF
116e No 0.64 0.35 -39 0 15 0.12 INF INF 5
116n No 0.64 0.35 -39 0 32 0.24 INF INF 2
117e No 0.64 0.38 -39 0 18 0.24 INF INF 4
117n No 0.64 0.38 -39 0 15 0.25 INF INF 4
118e No 0.65 0.38 -39 0 18 0.2 INF INF 2
118n No 0.65 0.38 -39 0 25 0.23 INF INF 1
119e No 0.62 0.37 -39 0 27 0.091 INF INF 4
119n No 0.62 0.37 -39 0 21 0.19 INF INF 4
120e Yes -39 1.5E+03 0 0 INF INF
120n Yes -39 1.5E+03 0 0 INF INF
121e No 0.63 0.37 -39 0 10 0.46 INF INF 4
121n No 0.66 0.37 -39 0 13 0.16 INF INF 3
122e No 0.51 0.29 -39 0 16 0.09 INF INF 13
122n No 0.5 0.29 -39 0 14 0.18 INF INF 10
123e No 0.66 0.39 -39 0 15 0.097 INF INF 4
123n No 0.66 0.39 -39 0 15 0.17 INF INF 2
124e No 0.55 0.33 -39 0 32 0.19 INF INF 25
124n No 0.55 0.33 -39 0 15 0.29 INF INF 23
125e No 0.55 0.34 -39 0 19 0.1 INF INF 16
125n No 0.52 0.34 -39 0 17 0.28 INF INF 13
126e No 0.55 0.34 -39 0 14 0.14 INF INF 14
126n No 0.54 0.34 -39 0 14 0.12 INF INF 10
127e No 0.33 0.23 -39 0 18 0.17 INF INF
127n No 0.3 0.23 -39 0 15 0.23 INF INF
128e No 0.32 0.22 -39 0 17 0.16 INF INF
128n No 0.29 0.22 -39 0 7.6 0.33 INF INF
129e No 0.32 0.23 -39 0 18 0.17 INF INF
129n No 0.31 0.23 -39 0 17 0.21 INF INF
130e No 0.31 0.22 -39 0 21 0.16 INF INF
130n No 0.29 0.22 -39 0 16 0.3 INF INF
131e No 0.25 0.17 -39 0 14 0.28 INF INF
131n No 0.24 0.17 -39 0 12 0.36 INF INF
132e No 0.27 0.19 -39 0 20 0.24 INF INF
132n No 0.25 0.19 -39 0 16 0.32 INF INF
133e No 0.25 0.17 -39 0 17 0.22 INF INF
133n No 0.24 0.17 -39 0 16 0.3 INF INF
134e No 0.26 0.18 -39 0 34 0.15 INF INF
134n No 0.24 0.18 -39 0 19 0.39 INF INF
135e No 0.029 0.39 -39 0 2.9 2.1 INF INF
135n No 0.59 0.39 -39 0 4.5 0.27 INF INF 5
136e No 0.61 0.35 -39 0 18 0.088 INF INF 7
136n No 0.59 0.35 -39 0 15 0.48 INF INF 2
137e No 0.65 0.38 -39 0 25 0.16 INF INF 2
137n No 0.66 0.38 -39 0 33 0.18 INF INF 4
138e No 0.62 0.36 -39 0 18 0.17 INF INF 4
138n No 0.62 0.36 -39 0 14 0.28 INF INF 4
139e No 0.64 0.38 -39 0 32 0.16 INF INF 2
139n No 0.64 0.38 -39 0 24 0.24 INF INF 1
140e No 0.65 0.38 -39 0 15 0.17 INF INF 3
140n No 0.65 0.38 -39 0 16 0.27 INF INF 1
141e No 0.66 0.38 -39 0 17 0.16 INF INF 4
141n No 0.66 0.38 -39 0 16 0.22 INF INF 2
142e No 0.65 0.37 -39 0 19 0.24 INF INF 2
142n No 0.64 0.37 -39 0 16 0.3 INF INF 3
143e No 0.64 0.37 -39 0 24 0.2 INF INF 8
143n No 0.63 0.37 -39 0 26 0.2 INF INF 6
144e No 0.48 0.29 -39 0 14 0.18 INF INF 24
144n No 0.47 0.29 -39 0 14 0.23 INF INF 19
145e No 0.5 0.3 -39 0 17 0.13 INF INF 23
145n No 0.48 0.3 -39 0 17 0.29 INF INF 17
146e No 0.48 0.3 -39 0 18 0.27 INF INF 17
146n No 0.45 0.3 -39 0 13 0.37 INF INF 11
147e Yes -39 1.5E+03 0 0 INF INF
147n Yes -39 1.5E+03 0 0 INF INF
148e Yes -39 1.5E+03 0 0 INF INF
148n Yes -39 1.5E+03 0 0 INF INF
149e Yes -39 1.5E+03 0 0 INF INF
149n Yes -39 1.5E+03 0 0 INF INF
150e No 0.43 0.31 -39 0 17 0.16 INF INF 21
150n No 0.42 0.31 -39 0 15 0.25 INF INF 20
151e No 0.5 0.37 -39 0 29 0.19 INF INF 8
151n No 0.48 0.37 -39 0 16 0.28 INF INF 8
152e No 0.48 0.34 -39 0 16 0.22 INF INF 19
152n No 0.48 0.34 -39 0 17 0.29 INF INF 17
153e No 0.46 0.33 -39 0 16 0.25 INF INF 12
153n No 0.46 0.33 -39 0 17 0.38 INF INF 12
154e No 0.46 0.33 -39 0 18 0.2 INF INF 11
154n No 0.45 0.33 -39 0 17 0.33 INF INF 9
155e No 0.61 0.33 -39 0 16 0.1 INF INF 10
155n No 0.61 0.33 -39 0 17 0.19 INF INF 10
156e No 0.62 0.35 -39 0 18 0.17 INF INF 5
156n No 0.61 0.35 -39 0 14 0.29 INF INF 4
157e No 0.63 0.37 -39 0 15 0.12 INF INF 4
157n No 0.63 0.37 -39 0 15 0.26 INF INF 3
158e No 0.63 0.38 -39 0 16 0.24 INF INF 4
158n No 0.64 0.38 -39 1 16 0.24 INF INF 2
159e No 0.6 0.36 -39 0 12 0.32 INF INF 1
159n No 0.62 0.36 -39 0 18 0.31 INF INF 1
160e No 0.64 0.37 -39 0 15 0.15 INF INF 2
160n No 0.64 0.37 -39 0 17 0.28 INF INF 1
161e No 0.59 0.3 -39 0 18 0.29 INF INF 5
161n No 0.46 0.3 -39 0 14 1.1 INF INF
162e No 0.61 0.35 -39 0 19 0.22 INF INF 6
162n No 0.61 0.35 -39 0 16 0.27 INF INF 5
163e No 0.63 0.38 -39 0 10 0.12 INF INF 7
163n No 0.62 0.38 -39 0 16 0.17 INF INF 8
164e No 0.63 0.38 -39 0 12 0.2 INF INF 7
164n No 0.63 0.38 -39 0 12 0.23 INF INF 6
165e No 0.63 0.39 -39 0 11 0.11 INF INF 8
165n No 0.61 0.39 -39 0 8.1 0.32 INF INF 5
166e No 0.63 0.39 -39 0 11 0.21 INF INF 4
166n No 0.62 0.39 -39 0 16 0.43 INF INF 4
167e No 0.46 0.31 -39 0 14 0.13 INF INF 31
167n No 0.45 0.31 -39 0 37 0.16 INF INF 29
168e No 0.32 0.22 -39 0 23 0.16 INF INF
168n No 0.3 0.22 -39 0 17 0.21 INF INF
169e No 0.29 0.2 -39 0 7.6 0.28 INF INF
169n No 0.3 0.2 -39 0 16 0.2 INF INF
170e No 0.022 0.21 -39 0 0.65 0.57 INF INF
170n No 0.29 0.21 -39 0 12 0.32 INF INF
171e No 0.33 0.36 -39 0 4.7 0.45 INF INF
171n No 0.48 0.36 -39 0 15 0.29 INF INF 8
172e No 0.5 0.37 -39 0 21 0.22 INF INF 10
172n No 0.47 0.37 -39 0 14 0.49 INF INF 4
173e No 0.46 0.32 -39 0 15 0.28 INF INF 18
173n No 0.47 0.32 -39 0 16 0.33 INF INF 18
174e No 0.43 0.31 -39 0 19 0.25 INF INF 28
174n No 0.42 0.31 -39 0 31 0.33 INF INF 22
175e No 0.4 0.29 -39 0 18 0.27 INF INF 16
175n No 0.33 0.29 -39 0 6.3 0.47 INF INF
176e No 0.63 0.38 -39 0 43 0.056 INF INF 3
176n No 0.63 0.38 -39 0 14 0.19 INF INF 3
177e No 0.63 0.39 -39 0 11 0.14 INF INF 2
177n No 0.63 0.39 -39 0 14 0.21 INF INF 1
178e No 0.64 0.38 -39 0 18 0.3 INF INF 4
178n No 0.65 0.38 -39 0 25 0.17 INF INF 1
179e No 0.65 0.39 -39 0 15 0.24 INF INF 3
179n No 0.65 0.39 -39 0 15 0.29 INF INF 1
180e No 0.65 0.37 -39 0 17 0.23 INF INF 2
180n No 0.51 0.37 -39 0 16 0.95 INF INF
181e No 0.66 0.38 -39 0 17 0.15 INF INF 4
181n No 0.66 0.38 -39 0 12 0.3 INF INF 1
182e No 0.051 0.43 -39 1 0.69 0.54 INF INF
182n No 0.61 0.43 -39 0 11 0.23 INF INF 6
183e No 0.66 0.38 -39 0 22 0.17 INF INF 4
183n No 0.66 0.38 -39 0 23 0.22 INF INF 8
184e No 0.64 0.39 -39 0 11 0.15 INF INF 4
184n No 0.62 0.39 -39 0 15 0.34 INF INF 3
185e No 0.63 0.38 -39 0 15 0.19 INF INF 8
185n No 0.62 0.38 -39 0 27 0.16 INF INF 6
186e No 0.49 0.37 -39 0 35 0.99 INF INF
186n No 0.61 0.37 -39 0 22 0.25 INF INF 7
187e No 0.64 0.4 -39 0 19 0.2 INF INF 4
187n No 0.62 0.4 -39 0 11 0.31 INF INF 4
188e No 0.44 0.28 -39 0 17 0.23 INF INF 22
188n No 0.38 0.28 -39 0 36 0.69 INF INF
189e No 0.35 0.23 -39 0 15 0.2 INF INF 21
189n No 0.34 0.23 -39 0 20 0.21 INF INF
190e No 0.36 0.24 -39 0 17 0.27 INF INF 23
190n No 0.35 0.24 -39 0 14 0.23 INF INF 20
191e No 0.36 0.25 -39 0 19 0.18 INF INF 21
191n No 0.35 0.25 -39 0 27 0.27 INF INF 18
192e No 0.48 0.34 -39 0 21 0.23 INF INF 20
192n No 0.48 0.34 -39 0 26 0.31 INF INF 19
193e No 0.46 0.33 -39 0 19 0.39 INF INF 19
193n No 0.46 0.33 -39 0 15 0.34 INF INF 15
194e No 0.41 0.29 -39 0 15 0.22 INF INF 19
194n No 0.42 0.29 -39 0 20 0.23 INF INF 20
195e No 0.41 0.31 -39 0 21 0.19 INF INF 15
195n No 0.4 0.31 -39 0 20 0.35 INF INF 15
196e Yes -39 1.5E+03 0 0 INF INF
196n Yes -39 1.5E+03 0 0 INF INF
197e Yes -39 1.5E+03 0 0 INF INF
197n Yes -39 1.5E+03 0 0 INF INF
198e Yes -39 1.5E+03 0 0 INF INF
198n Yes -39 1.5E+03 0 0 INF INF
200e No 0.031 0.43 -39 0 3.1 0.58 INF INF
200n No 0.6 0.43 -39 0 23 0.25 INF INF 8
201e No 0.58 0.37 -39 0 54 0.24 INF INF 6
201n No 0.61 0.37 -39 0 42 0.15 INF INF 4
202e No 0.62 0.37 -39 0 32 0.2 INF INF 10
202n No 0.6 0.37 -39 0 19 0.33 INF INF 5
203e No 0.63 0.38 -39 0 15 0.21 INF INF 6
203n No 0.63 0.38 -39 0 15 0.21 INF INF 6
204e No 0.62 0.37 -39 0 12 -0.099 INF INF 15
204n No 0.61 0.37 -39 0 15 -0.13 INF INF 14
205e No 0.59 0.37 -39 0 22 0.18 INF INF 17
205n No 0.58 0.37 -39 0 20 0.24 INF INF 10
206e No 0.59 0.38 -39 0 20 0.25 INF INF 17
206n No 0.58 0.38 -39 0 19 0.29 INF INF 10
207e No 0.61 0.41 -39 0 23 0.24 INF INF 17
207n No 0.58 0.41 -39 0 18 0.49 INF INF 12
208e No 0.59 0.4 -39 0 26 0.15 INF INF 17
208n No 0.57 0.4 -39 0 38 0.27 INF INF 13
209e No 0.57 0.39 -39 0 25 0.18 INF INF 12
209n No 0.57 0.39 -39 0 29 0.33 INF INF 11
210e No 0.58 0.4 -39 0 19 -0.048 INF INF 14
210n No 0.57 0.4 -39 0 12 0.0071 INF INF 16
211e No 0.56 0.39 -39 0 29 0.22 INF INF 13
211n No 0.55 0.39 -39 0 20 0.29 INF INF 14
212e No 0.44 0.33 -39 0 24 0.22 INF INF 14
212n No 0.39 0.33 -39 0 9.8 0.39 INF INF 11
213e No 0.42 0.3 -39 0 17 0.24 INF INF 21
213n No 0.43 0.3 -39 0 18 0.27 INF INF 19
214e No 0.42 0.31 -39 0 19 0.23 INF INF 24
214n No 0.41 0.31 -39 0 19 0.33 INF INF 17
215e No 0.61 0.38 -39 0 33 0.14 INF INF 3
215n No 0.57 0.38 -39 0 16 0.46 INF INF 1
216e No 0.6 0.38 -39 0 20 0.24 INF INF 2
216n No 0.59 0.38 -39 0 13 0.37 INF INF 1
217e No 0.62 0.38 -39 0 25 0.26 INF INF 6
217n No 0.61 0.38 -39 0 17 0.32 INF INF 2
218e No 0.61 0.42 -39 0 27 0.17 INF INF 6
218n No 0.16 0.42 -39 1 0.54 0.24 INF INF
219e No 0.63 0.39 -39 0 20 0.17 INF INF 11
219n No 0.61 0.39 -39 0 16 0.3 INF INF 10
220e No 0.62 0.38 -39 0 19 0.24 INF INF 13
220n No 0.6 0.38 -39 0 15 0.27 INF INF 12
221e No 0.63 0.38 -39 0 18 0.28 INF INF 9
221n No 0.63 0.38 -39 0 15 0.34 INF INF 8
222e No 0.63 0.39 -39 0 18 0.3 INF INF 10
222n No 0.64 0.39 -39 0 24 0.26 INF INF 8
223e No 0.61 0.37 -39 0 16 0.2 INF INF 17
223n No 0.62 0.37 -39 0 21 0.2 INF INF 15
224e No 0.6 0.39 -39 0 23 0.4 INF INF 14
224n No 0.62 0.39 -39 0 20 0.28 INF INF 15
225e No 0.61 0.39 -39 0 25 0.19 INF INF 8
225n No 0.59 0.39 -39 0 20 0.29 INF INF 8
226e No 0.6 0.38 -39 0 23 0.18 INF INF 8
226n No 0.59 0.38 -39 0 29 0.28 INF INF 8
227e No 0.55 0.38 -39 0 19 0.24 INF INF 22
227n No 0.52 0.38 -39 0 12 0.39 INF INF 16
228e No 0.52 0.36 -39 0 18 0.32 INF INF 18
228n No 0.52 0.36 -39 0 17 0.41 INF INF 20
229e No 0.54 0.39 -39 0 17 0.26 INF INF 22
229n No 0.54 0.39 -39 0 15 0.27 INF INF 22
231e No 0.44 0.32 -39 0 37 0.12 INF INF 16
231n No 0.44 0.32 -39 0 28 0.16 INF INF 14
232e No 0.43 0.32 -39 0 44 0.11 INF INF 16
232n No 0.4 0.32 -39 0 14 0.35 INF INF 15
233e No 0.59 0.35 -39 0 18 0.2 INF INF 11
233n No 0.47 0.35 -39 0 35 0.92 INF INF
234e No 0.6 0.37 -39 0 26 0.22 INF INF 11
234n No 0.6 0.37 -39 0 28 0.25 INF INF 6
237e No 0.61 0.37 -39 0 18 0.27 INF INF 11
237n No 0.61 0.37 -39 0 19 0.3 INF INF 10
238e No 0.64 0.39 -39 0 28 0.16 INF INF 11
238n No 0.62 0.39 -39 0 18 0.32 INF INF 10
239e No 0.63 0.39 -39 0 22 0.2 INF INF 11
239n No 0.63 0.39 -39 0 27 0.23 INF INF 10
240e No 0.54 0.38 -39 0 7.9 0.43 INF INF 15
240n No 0.59 0.38 -39 0 15 0.38 INF INF 10
241e No 0.61 0.39 -39 0 19 0.19 INF INF 8
241n No 0.6 0.39 -39 0 16 0.35 INF INF 6
242e No 0.62 0.39 -39 32 24 0.17 INF INF
242n No 0.62 0.39 -39 32 28 0.19 INF INF
243e No 0.51 0.38 -39 0 33 0.86 INF INF
243n No 0.6 0.38 -39 0 17 0.28 INF INF 8
244e No 0.56 0.37 -39 0 16 0.19 INF INF 14
244n No 0.56 0.37 -39 0 17 0.32 INF INF 22
245e No 0.57 0.39 -39 0 27 0.21 INF INF 22
245n No 0.55 0.39 -39 0 20 0.36 INF INF 20
246e No 0.55 0.37 -39 0 18 0.31 INF INF 10
246n No 0.55 0.37 -39 0 24 0.32 INF INF 10
250e No 0.6 0.38 -39 0 22 0.19 INF INF 9
250n No 0.6 0.38 -39 0 18 0.3 INF INF 8
251e No 0.053 0.084 -39 0 5.7 0.66 INF INF
251n No 0.2 0.084 -39 0 21 0.21 INF INF
252e No 0.6 0.37 -39 0 17 0.27 INF INF 6
252n No 0.6 0.37 -39 0 16 0.3 INF INF 6
253e No 0.61 0.37 -39 0 18 0.29 INF INF 9
253n No 0.59 0.37 -39 0 11 0.4 INF INF 2
254e No 0.62 0.39 -39 0 33 0.18 INF INF 4
254n No 0.61 0.39 -39 0 28 0.35 INF INF 7
255e No 0.61 0.38 -39 0 15 0.31 INF INF 2
255n No 0.62 0.38 -39 0 42 0.25 INF INF 3
256e No 0.62 0.39 -39 0 26 0.3 INF INF 9
256n No 0.6 0.39 -39 0 11 0.43 INF INF 2
257e No 0.61 0.38 -39 0 38 0.17 INF INF 13
257n No 0.62 0.38 -39 0 23 0.21 INF INF 12
261e No 0.44 0.31 -39 0 5.8 0.4 INF INF 6
261n No 0.44 0.31 -39 0 5.3 0.48 INF INF 4
262e No 0.56 0.38 -39 0 13 0.013 INF INF 13
262n No 0.55 0.38 -39 0 15 0.11 INF INF 14
266e No 0.61 0.4 -39 0 8.9 -0.14 INF INF 7
266n No 0.61 0.4 -39 0 13 -0.19 INF INF 7
267e No 0.55 0.41 -39 0 8 0.32 INF INF 3
267n No 0.61 0.41 -39 0 14 0.25 INF INF 8
268e No 0.59 0.4 -39 0 10 0.31 INF INF 9
268n No 0.62 0.4 -39 0 18 0.31 INF INF 8
269e No 0.61 0.38 -39 0 18 0.27 INF INF 11
269n No 0.6 0.38 -39 0 15 0.36 INF INF 8
270e No 0.21 -0.34 -39 0 25 0.36 INF INF
270n No 0.22 -0.34 -39 0 23 0.17 INF INF
271e No 0.039 0.5 -39 0 3.3 0.56 INF INF
271n No 0.63 0.5 -39 0 45 0.11 INF INF 5
272e No 0.6 0.39 -39 0 17 0.35 INF INF 8
272n No 0.62 0.39 -39 0 25 0.29 INF INF 9
273e No 0.06 0.47 -39 0 3.3 0.58 INF INF
273n No 0.62 0.47 -39 0 33 0.21 INF INF 4
277e No 0.29 0.29 -39 0 21 0.94 INF INF
277n No 0.43 0.29 -39 0 47 0.16 INF INF 18
278e No 0.28 0.25 -39 0 6.5 0.46 INF INF
278n No 0.36 0.25 -39 0 9.6 0.39 INF INF 15
281e No 0.58 0.38 -39 0 17 0.28 INF INF 9
281n No 0.59 0.38 -39 0 15 0.32 INF INF 2
282e No 0.61 0.4 -39 0 18 0.25 INF INF 9
282n No 0.62 0.4 -39 0 32 0.18 INF INF 8
283e No 0.62 0.39 -39 0 21 0.16 INF INF 6
283n No 0.61 0.39 -39 0 18 0.27 INF INF 10
284e No 0.21 -0.34 -39 0 35 0.21 INF INF
284n No 0.22 -0.34 -39 0 28 0.15 INF INF
285e No 0.61 0.38 -39 0 22 0.26 INF INF 9
285n No 0.59 0.38 -39 0 15 0.33 INF INF 8
286e No 0.2 0.12 -39 0 26 0.31 INF INF
286n No 0.033 0.12 -39 0 2.8 0.6 INF INF
287e No 0.64 0.5 -39 0 31 0.076 INF INF 2
287n No 0.039 0.5 -39 0 3.2 0.63 INF INF
290e No 0.55 0.37 -39 0 75 0.059 INF INF
290n No 0.56 0.37 -39 0 73 0.074 INF INF
291e No 0.56 0.46 -39 0 69 0.058 INF INF
291n No 0.049 0.46 -39 0 0.68 -0.0016 INF INF
292e No 0.41 0.3 -39 0 17 0.27 INF INF 18
292n No 0.42 0.3 -39 0 34 0.22 INF INF 18
293e No 0.35 0.32 -39 0 25 0.97 INF INF
293n No 0.47 0.32 -39 0 12 0.31 INF INF 20
294e No 0.4 0.33 -39 0 23 0.66 INF INF
294n No 0.47 0.33 -39 0 19 0.24 INF INF 15
295e No 0.61 0.48 -39 0 4.3 -0.27 INF INF 4
295n No 0.03 0.48 -39 0 0.79 0.089 INF INF
299e No 0.51 0.35 -39 0 23 0.24 INF INF 22
299n No 0.48 0.35 -39 0 11 0.33 INF INF 18
300e No 0.52 0.34 -39 0 24 0.14 INF INF 19
300n No 0.38 0.34 -39 0 22 1.1 INF INF
301e No 0.51 0.35 -39 0 29 0.25 INF INF 18
301n No 0.49 0.35 -39 0 12 0.32 INF INF 21
302e No 0.45 0.45 -39 0 4.8 0.46 INF INF 4
302n No 0.6 0.45 -39 0 18 0.26 INF INF 7
306e No 0.47 0.33 -39 0 13 0.23 INF INF 10
306n No 0.33 0.33 -39 0 18 1 INF INF
307e No 0.3 0.1 -39 0 17 1 INF INF
307n No 0.3 0.1 -39 0 15 1.1 INF INF
311e No 0.58 0.39 -39 0 15 0.24 INF INF 9
311n No 0.55 0.39 -39 0 9.4 0.48 INF INF 6
312e No 0.18 0.13 -39 1 17 0.36 INF INF
312n No 0.028 0.13 -39 0 2.9 0.59 INF INF
313e No 0.44 0.39 -39 0 33 0.79 INF INF
313n No 0.56 0.39 -39 0 40 0.17 INF INF 15
314e No 0.56 0.4 -39 0 65 0.068 INF INF
314n No 0.57 0.4 -39 0 62 0.098 INF INF
315e No 0.53 0.38 -39 0 17 0.21 INF INF 12
315n No 0.52 0.38 -39 0 16 0.36 INF INF 12
316e No 0.54 0.38 -39 0 22 0.21 INF INF 11
316n No 0.54 0.38 -39 0 22 0.24 INF INF 12
317e No 0.54 0.38 -39 0 30 0.16 INF INF 12
317n No 0.51 0.38 -39 0 13 0.32 INF INF 10
318e No 0.58 0.42 -39 0 61 0.023 INF INF
318n No 0.55 0.42 -39 0 76 0.037 INF INF
319e No 0.46 0.34 -39 0 23 0.25 INF INF 13
319n No 0.46 0.34 -39 0 25 0.26 INF INF 12
320e Yes -39 1.5E+03 0 0 INF INF
320n Yes -39 1.5E+03 0 0 INF INF
321e No 0.31 0.24 -39 0 12 0.1 INF INF
321n No 0.32 0.24 -39 0 14 0.15 INF INF
322e No 0.34 0.25 -39 0 51 0.093 INF INF
322n No 0.27 0.25 -39 0 7.5 0.44 INF INF
323e No 0.25 0.21 -39 0 12 0.33 INF INF
323n No 0.29 0.21 -39 0 23 0.17 INF INF
324e No 0.25 0.18 -39 0 28 0.16 INF INF
324n No 0.25 0.18 -39 0 29 0.18 INF INF
325e No 0.47 0.34 -39 0 26 0.16 INF INF 9
325n No 0.46 0.34 -39 0 19 0.29 INF INF 4
326e No 0.021 -0.00024 -39 0 3.5 0.56 INF INF
326n No 0.021 -0.00024 -39 0 3.1 0.58 INF INF
327e No 0.29 0.19 -39 0 25 0.22 INF INF
327n No 0.28 0.19 -39 0 24 0.27 INF INF
328e No 0.021 0.16 -39 1 3.4 0.56 INF INF
328n No 0.21 0.16 -39 0 19 0.27 INF INF
329e No 0.5 0.38 -39 0 21 0.2 INF INF 1
329n No 0.52 0.38 -39 0 26 0.16 INF INF 1
331e No 0.3 0.22 -39 0 29 0.18 INF INF
331n No 0.29 0.22 -39 0 22 0.37 INF INF
332e No 0.32 0.23 -39 0 29 0.17 INF INF
332n No 0.29 0.23 -39 0 13 0.36 INF INF
333e No 0.47 0.38 -39 0 8.9 0.23 INF INF 1
333n No 0.5 0.38 -39 0 12 0.24 INF INF 1
336e No 0.27 0.2 -39 0 23 0.29 INF INF
336n No 0.28 0.2 -39 0 48 0.14 INF INF
339e Yes -39 1.5E+03 0 0 INF INF
339n Yes -39 1.5E+03 0 0 INF INF
340e No 0.31 0.23 -39 0 34 0.042 INF INF
340n No 0.31 0.23 -39 0 24 0.15 INF INF
342e No 0.49 0.38 -39 0 14 0.27 INF INF 8
342n No 0.5 0.38 -39 32 25 0.28 INF INF
343e No 0.53 0.41 -39 0 21 0.25 INF INF 10
343n No 0.51 0.41 -39 0 14 0.42 INF INF 7
344e No 0.11 -0.3 -39 0 69 0.075 INF INF
344n No 0.12 -0.3 -39 0 53 0.08 INF INF
345e Yes -39 1.5E+03 0 0 INF INF
345n Yes -39 1.5E+03 0 0 INF INF
346e No 0.025 0.0017 -39 0 3.3 0.56 INF INF
346n No 0.023 0.0017 -39 0 2.7 0.6 INF INF
347e No 0.49 0.38 -39 0 29 0.2 INF INF 8
347n No 0.5 0.38 -39 0 49 0.16 INF INF 9
348e No 0.035 0.38 -39 0 3.4 0.61 INF INF
348n No 0.48 0.38 -39 0 62 0.1 INF INF
349e No 0.42 0.33 -39 0 55 0.058 INF INF 4
349n No 0.42 0.33 -39 0 73 0.052 INF INF
In [59]:
# 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)    
In [60]:
print('Final Ant-Pol Classification:\n\n', overall_class)
Final Ant-Pol Classification:

 Jee:
----------
bad (290 antpols):
3, 4, 5, 7, 8, 9, 10, 15, 16, 17, 18, 19, 20, 21, 27, 28, 29, 30, 31, 32, 33, 36, 37, 38, 40, 41, 42, 43, 44, 45, 46, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 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, 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, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 250, 251, 252, 253, 254, 255, 256, 257, 261, 262, 266, 267, 268, 269, 270, 271, 272, 273, 277, 278, 281, 282, 283, 284, 285, 286, 287, 290, 291, 292, 293, 294, 295, 299, 300, 301, 302, 306, 307, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 331, 332, 333, 336, 339, 340, 342, 343, 344, 345, 346, 347, 348, 349


Jnn:
----------
bad (290 antpols):
3, 4, 5, 7, 8, 9, 10, 15, 16, 17, 18, 19, 20, 21, 27, 28, 29, 30, 31, 32, 33, 36, 37, 38, 40, 41, 42, 43, 44, 45, 46, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 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, 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, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 250, 251, 252, 253, 254, 255, 256, 257, 261, 262, 266, 267, 268, 269, 270, 271, 272, 273, 277, 278, 281, 282, 283, 284, 285, 286, 287, 290, 291, 292, 293, 294, 295, 299, 300, 301, 302, 306, 307, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 331, 332, 333, 336, 339, 340, 342, 343, 344, 345, 346, 347, 348, 349

Save calibration solutions¶

In [61]:
if not all_flagged():
    # 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
All antennas are flagged, so this cell is being skipped.
In [62]:
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    
    
    if not all_flagged():
        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)
    
        del hd_vissol
        malloc_trim()        
All antennas are flagged, so this cell is being skipped.

Output fully flagged calibration file if OMNICAL_FILE is not written¶

In [63]:
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)
    # create fully flagged unit gains with chi^2 = 0
    hc_omni = hd_writer.init_HERACal(gain_convention='divide', cal_style='redundant')
    hc_omni.history += add_to_history
    hc_omni.pol_convention = hd_auto_model.pol_convention
    hc_omni.gain_scale = hd_auto_model.vis_units
    hc_omni.write_calfits(OMNICAL_FILE, clobber=True)
    del hc_omni
WARNING: No calibration file produced at /mnt/sn1/data1/2461003/zen.2461003.43882.sum.omni.calfits. Creating a fully-flagged placeholder calibration file.

Output empty visibility file if OMNIVIS_FILE is not written¶

In [64]:
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)

Metadata¶

In [65]:
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.2.5.dev1+g5a985ae31
hera_cal: 3.7.7.dev68+g3286222d3
hera_filters: 0.1.7
hera_qm: 2.2.1.dev4+gf6d02113b
hera_notebook_templates: 0.1.dev989+gee0995d
In [66]:
print(f'Finished execution in {(time.time() - tstart) / 60:.2f} minutes.')
Finished execution in 3.57 minutes.