Single File Calibration¶

by Josh Dillon, Aaron Parsons, Tyler Cox, and Zachary Martinot, last updated December 15, 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
bigmem1.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/2461081/zen.2461081.45950.sum.uvh5'
DIFF_FILE = '/mnt/sn1/data1/2461081/zen.2461081.45950.diff.uvh5'
AM_FILE = '/mnt/sn1/data1/2461081/zen.2461081.45950.sum.ant_metrics.hdf5'
ANTCLASS_FILE = '/mnt/sn1/data1/2461081/zen.2461081.45950.sum.ant_class.csv'
OMNICAL_FILE = '/mnt/sn1/data1/2461081/zen.2461081.45950.sum.omni.calfits'
OMNIVIS_FILE = '/mnt/sn1/data1/2461081/zen.2461081.45950.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_RERUN_MAXITER", 250))
OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE = float(os.environ.get("OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE", 1.5))
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_RESTART_FROM_FC_EVERY_ITER = os.environ.get("OC_RESTART_FROM_FC_EVERY_ITER", "FALSE").upper() == "TRUE"
OC_USE_GPU = os.environ.get("OC_USE_GPU", "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", "TRUE").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', 'OC_RESTART_FROM_FC_EVERY_ITER', '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 = 25
OC_MAX_RERUN = 10
OC_RERUN_MAXITER = 125
OC_MAX_CHISQ_FLAGGING_DYNAMIC_RANGE = 1.5
OC_USE_PRIOR_SOL = True
OC_PRIOR_SOL_FLAG_THRESH = 0.95
OC_USE_GPU = False
OC_RESTART_FROM_FC_EVERY_ITER = 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 0.76 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/2461081/zen.2461081.45950.sum.uvh5
JDs: [2461081.45944446 2461081.4595563 ] (9.66368 s integrations)
LSTS: [9.79180822 9.79449993] 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.')
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)
All antennas are flagged, so this cell is being skipped.

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
All antennas are flagged, so this cell is being skipped.

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
All antennas are flagged, so this cell is being skipped.
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'], thresh=.75):
    '''Checks if nearly an entire polarization is flagged (depending on thresh). 
    If it is, returns True and marks all antennas as bad in redcal_class.'''
    flag_fracs = np.array([np.mean([redcal_class[ant] == 'bad' for ant in redcal_class.ants if ant[1] == pol]) for pol in pols])
    if np.any(flag_fracs > thresh):
        for pol, frac in zip(pols, flag_fracs):
            if frac > thresh:
                print(f'Polarization {pol} is {frac:.3%} flagged > {thresh:.3%} threshold. 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:
            sol = None
            
        malloc_trim()
All antennas are flagged, so this cell is being skipped.
In [37]:
if not all_flagged():
    all_high_chisq_found = False
    metric = 'median'
    
    # iteratively rerun redundant calibration
    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 because an entire pol is flagged
        if check_if_whole_pol_flagged(redcal_class):
            break
    
        # change settings for final run
        if all_high_chisq_found or (i == OC_MAX_RERUN):
            sol = None  # start from scratch
            rc_settings['oc_maxiter'] = rc_settings['check_after'] = OC_RERUN_MAXITER
        
        # optionally redo firstcal and start from there every iteration
        if OC_RESTART_FROM_FC_EVERY_ITER:
            sol = None

        # re-run redundant calibration using previous solution (if it's not the final run and if not OC_RESTART_FROM_FC_EVERY_ITER)
        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
        avg_cspa = recheck_chisq(meta['chisq_per_ant'], sol, oc_cspa_suspect[1], (np.nanmean if metric == 'mean' else np.nanmedian))

        # flag bad antennas
        cspa_class = ant_class.antenna_bounds_checker(avg_cspa, good=oc_cspa_good, suspect=oc_cspa_suspect, bad=[(-np.inf, np.inf)])
        for ant in cspa_class.bad_ants:
            if avg_cspa[ant] < np.max(list(avg_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

        if len(cspa_class.bad_ants) > 0:
            print(f'Removing {cspa_class.bad_ants} for high {metric} unflagged chi^2.')
            for ant in cspa_class.bad_ants:
                print(f'\t{ant}: {avg_cspa[ant]:.3f}')
    
        if check_if_whole_pol_flagged(redcal_class) or all_high_chisq_found or (i == OC_MAX_RERUN):
            break

        if len(cspa_class.bad_ants) == 0:
            if metric == 'median':
                metric = 'mean'  # switch to mean if no new median offenders found
            else:
                all_high_chisq_found = True  # no new antennas to flag, the next iteration will be the final one
    
    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
            max_dly = np.max(np.abs(list(meta['fc_meta']['dlys'].values())))
            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.52 0.31 -45 1.9E+02 13 -0.0047
3n No 0.52 0.31 -45 1.9E+02 18 -0.13
4e No 0.35 0.15 -45 1.9E+02 14 0.041
4n No 0.33 0.15 -45 1.9E+02 16 0.2
5e No 0.57 0.34 -45 1.9E+02 17 -0.1
5n No 0.56 0.34 -45 1.9E+02 21 -0.17
7e No 0.56 0.33 -45 1.9E+02 23 -0.049
7n No 0.56 0.33 -45 1.9E+02 23 -0.13
8e No 0.56 0.34 -45 1.9E+02 26 -0.075
8n No 0.56 0.34 -45 1.9E+02 12 -0.055
9e No 0.56 0.34 -45 1.9E+02 20 -0.038
9n No 0.56 0.34 -45 1.9E+02 21 -0.11
10e No 0.56 0.34 -45 1.9E+02 16 -0.029
10n No 0.55 0.34 -45 1.9E+02 16 -0.051
15e No 0.56 0.34 -45 1.9E+02 33 -0.1
15n No 0.55 0.34 -45 1.9E+02 36 -0.19
16e No 0.58 0.35 -45 1.9E+02 19 -0.054
16n No 0.57 0.35 -45 1.9E+02 18 -0.035
17e No 0.58 0.34 -45 1.9E+02 19 -0.081
17n No 0.57 0.34 -45 1.9E+02 21 -0.096
18e No 0.029 0.0091 -45 1.9E+02 0.64 0.33
18n No 0.036 0.0091 -45 1.9E+02 0.57 0.4
19e No 0.58 0.35 -45 1.9E+02 21 -0.069
19n No 0.57 0.35 -45 1.9E+02 20 -0.14
20e No 0.56 0.33 -45 1.9E+02 22 0.021
20n No 0.56 0.33 -45 1.9E+02 23 -0.14
21e No 0.54 0.32 -45 1.9E+02 8.9 0.031
21n No 0.54 0.32 -45 1.9E+02 9.3 -0.027
27e No 0.15 0.0066 -45 1.9E+02 11 0.025
27n No 0.1 0.0066 -45 1.9E+02 7.2 0.37
28e No 0.12 0.05 -45 1.9E+02 6.3 0.42
28n No 0.21 0.05 -45 1.9E+02 2.7 0.054
29e No 0.58 0.33 -45 1.9E+02 22 -0.043
29n No 0.51 0.33 -45 1.9E+02 18 0.57
30e No 0.57 0.34 -45 1.9E+02 18 -0.031
30n No 0.52 0.34 -45 1.9E+02 2.6 -0.043
31e No 0.59 0.35 -45 1.9E+02 18 -0.11
31n No 0.58 0.35 -45 1.9E+02 18 -0.085
32e No 0.55 0.32 -45 1.9E+02 5.9 0.0037
32n No 0.44 0.32 -45 1.9E+02 12 0.86
33e No 0.57 0.39 -45 1.9E+02 19 -0.021
33n No 0.37 0.39 -45 1.9E+02 24 -0.14
36e No 0.56 0.34 -45 1.9E+02 27 -0.37
36n No 0.54 0.34 -45 1.9E+02 23 -0.41
37e No 0.57 0.35 -45 1.9E+02 13 -0.055
37n No 0.56 0.35 -45 1.9E+02 27 -0.13
38e No 0.56 0.33 -45 1.9E+02 22 -0.011
38n No 0.55 0.33 -45 1.9E+02 32 -0.14
40e No 0.58 0.35 -45 1.9E+02 13 -0.053
40n No 0.54 0.35 -45 1.9E+02 3.5 -0.084
41e No 0.59 0.34 -45 1.9E+02 22 -0.14
41n No 0.58 0.34 -45 1.9E+02 29 -0.16
42e No 0.59 0.37 -45 1.9E+02 22 -0.087
42n No 0.52 0.37 -45 1.9E+02 1.8 -0.047
43e No 0.59 0.34 -45 1.9E+02 20 -0.092
43n No 0.58 0.34 -45 1.9E+02 18 -0.053
44e No 0.15 0.024 -45 1.9E+02 0.97 0.67
44n No 0.13 0.024 -45 1.9E+02 0.78 0.62
45e No 0.58 0.33 -45 1.9E+02 20 -0.083
45n No 0.57 0.33 -45 1.9E+02 22 -0.059
46e No 0.56 0.34 -45 1.9E+02 11 2.1
46n No 0.57 0.34 -45 1.9E+02 19 -0.12
50e No 0.57 0.36 -45 1.9E+02 20 -0.063
50n No 0.57 0.36 -45 1.9E+02 19 -0.17
51e No 0.58 0.35 -45 1.9E+02 20 -0.062
51n No 0.57 0.35 -45 1.9E+02 18 -0.18
52e No 0.59 0.34 -45 1.9E+02 22 -0.32
52n No 0.58 0.34 -45 1.9E+02 20 -0.38
53e No 0.59 0.35 -45 1.9E+02 16 -0.12
53n No 0.59 0.35 -45 1.9E+02 25 -0.26
54e No 0.59 0.33 -45 1.9E+02 18 -0.037
54n No 0.58 0.33 -45 1.9E+02 20 -0.11
55e No 0.58 0.33 -45 1.9E+02 19 -0.054
55n No 0.57 0.33 -45 1.9E+02 20 -0.14
56e No 0.6 0.34 -45 1.9E+02 25 -0.1
56n No 0.59 0.34 -45 1.9E+02 24 -0.076
57e No 0.6 0.46 -45 1.9E+02 23 -0.083
57n No 0.049 0.46 -45 1.9E+02 0.63 0.37
58e No 0.6 0.34 -45 1.9E+02 17 0.018
58n No 0.59 0.34 -45 1.9E+02 18 -0.019
59e No 0.59 0.34 -45 1.9E+02 29 -0.035
59n No 0.59 0.34 -45 1.9E+02 23 -0.053
60e No 0.59 0.35 -45 1.9E+02 15 0.013
60n No 0.57 0.35 -45 1.9E+02 6.8 -0.0065
65e No 0.58 0.36 -45 1.9E+02 19 -0.088
65n No 0.57 0.36 -45 1.9E+02 19 -0.15
66e No 0.6 0.36 -45 1.9E+02 23 -0.098
66n No 0.59 0.36 -45 1.9E+02 24 -0.15
67e No 0.6 0.35 -45 1.9E+02 22 -0.09
67n No 0.59 0.35 -45 1.9E+02 30 -0.14
68e No 0.6 0.35 -45 1.9E+02 17 -0.098
68n No 0.59 0.35 -45 1.9E+02 16 -0.18
69e No 0.6 0.34 -45 1.9E+02 19 -0.045
69n No 0.59 0.34 -45 1.9E+02 23 -0.11
70e No 0.6 0.34 -45 1.9E+02 18 -0.057
70n No 0.59 0.34 -45 1.9E+02 25 -0.069
71e No 0.61 0.34 -45 1.9E+02 16 -0.083
71n No 0.59 0.34 -45 1.9E+02 9.6 -0.073
72e No 0.61 0.34 -45 1.9E+02 19 0.071
72n No 0.59 0.34 -45 1.9E+02 14 -0.13
73e No 0.6 0.34 -45 1.9E+02 20 0.0095
73n No 0.59 0.34 -45 1.9E+02 17 -0.066
74e No 0.6 0.34 -45 1.9E+02 14 -0.047
74n No 0.59 0.34 -45 1.9E+02 14 -0.11
75e No 0.22 -0.13 -45 1.9E+02 1.6 0.097
75n No 0.15 -0.13 -45 1.9E+02 0.73 0.23
76e No 0.032 -0.00022 -45 1.9E+02 2.9 0.34
76n No 0.033 -0.00022 -45 1.9E+02 2.6 0.35
79e No 0.56 0.33 -45 1.9E+02 20 0.053
79n No 0.56 0.33 -45 1.9E+02 28 -0.063
80e No 0.51 0.31 -45 1.9E+02 12 0.13
80n No 0.52 0.31 -45 1.9E+02 14 0.087
81e No 0.52 0.36 -45 1.9E+02 0.74 -0.37
81n No 0.57 0.36 -45 1.9E+02 3 -0.51
82e No 0.58 0.34 -45 1.9E+02 17 0.07
82n No 0.58 0.34 -45 1.9E+02 16 -0.11
83e No 0.57 0.33 -45 1.9E+02 6 0.014
83n No 0.55 0.33 -45 1.9E+02 17 0.39
84e No 0.59 0.33 -45 1.9E+02 19 0.069
84n No 0.59 0.33 -45 1.9E+02 20 -0.068
85e No 0.61 0.34 -45 1.9E+02 16 -0.15
85n No 0.58 0.34 -45 1.9E+02 21 0.13
86e No 0.27 -0.3 -45 1.9E+02 21 -0.02
86n No 0.27 -0.3 -45 1.9E+02 29 -0.16
87e No 0.48 0.33 -45 1.9E+02 13 0.94
87n No 0.6 0.33 -45 1.9E+02 23 -0.15
88e Yes -45 1.5E+03 0 0
88n Yes -45 1.5E+03 0 0
89e No 0.048 0.00026 -45 1.9E+02 0.93 0.32
89n No 0.046 0.00026 -45 1.9E+02 0.87 0.35
90e Yes -45 1.5E+03 0 0
90n Yes -45 1.5E+03 0 0
91e No 0.6 0.35 -45 1.9E+02 19 -0.12
91n No 0.58 0.35 -45 1.9E+02 13 -0.13
92e No 0.59 0.34 -45 1.9E+02 13 -0.034
92n No 0.59 0.34 -45 1.9E+02 15 -0.077
93e No 0.54 0.36 -45 1.9E+02 2.3 -0.018
93n No 0.59 0.36 -45 1.9E+02 18 -0.099
94e No 0.58 0.33 -45 1.9E+02 26 -0.049
94n No 0.57 0.33 -45 1.9E+02 33 -0.1
95e No 0.55 0.33 -45 1.9E+02 19 0.082
95n No 0.56 0.33 -45 1.9E+02 28 -0.049
96e No 0.55 0.33 -45 1.9E+02 16 0.077
96n No 0.55 0.33 -45 1.9E+02 17 -0.0024
97e No 0.56 0.35 -45 1.9E+02 19 -0.034
97n No 0.51 0.35 -45 1.9E+02 8.3 0.13
98e No 0.57 0.35 -45 1.9E+02 14 -0.082
98n No 0.56 0.35 -45 1.9E+02 13 -0.18
99e No 0.18 0.41 -45 1.9E+02 1 0.39
99n No 0.57 0.41 -45 1.9E+02 62 -0.2
100e No 0.59 0.35 -45 1.9E+02 25 0.055
100n No 0.6 0.35 -45 1.9E+02 21 -0.083
101e No 0.61 0.34 -45 1.9E+02 19 -0.13
101n No 0.6 0.34 -45 1.9E+02 17 -0.19
102e No 0.61 0.35 -45 1.9E+02 21 -0.089
102n No 0.6 0.35 -45 1.9E+02 21 -0.069
103e No 0.57 0.34 -45 1.9E+02 4.5 0.04
103n No 0.59 0.34 -45 1.9E+02 25 -0.081
104e No 0.61 0.45 -45 1.9E+02 14 0.0061
104n No 0.025 0.45 -45 1.9E+02 0.046 0.14
105e No 0.61 0.34 -45 1.9E+02 16 -0.1
105n No 0.6 0.34 -45 1.9E+02 17 -0.17
106e No 0.61 0.35 -45 1.9E+02 19 -0.19
106n No 0.6 0.35 -45 1.9E+02 16 -0.28
107e Yes -45 1.5E+03 0 0
107n Yes -45 1.5E+03 0 0
108e No 0.6 0.34 -45 1.9E+02 17 0.036
108n No 0.59 0.34 -45 1.9E+02 14 0.02
109e No 0.6 0.34 -45 1.9E+02 16 -0.017
109n No 0.59 0.34 -45 1.9E+02 18 -0.087
110e No 0.6 0.34 -45 1.9E+02 17 0.016
110n No 0.59 0.34 -45 1.9E+02 18 -0.052
111e No 0.59 0.34 -45 1.9E+02 21 -0.078
111n No 0.59 0.34 -45 1.9E+02 20 -0.12
112e No 0.6 0.34 -45 1.9E+02 26 -0.089
112n No 0.59 0.34 -45 1.9E+02 18 -0.15
113e No 0.57 0.33 -45 1.9E+02 23 0.006
113n No 0.52 0.33 -45 1.9E+02 20 0.49
114e No 0.57 0.33 -45 1.9E+02 21 0.044
114n No 0.56 0.33 -45 1.9E+02 21 -0.068
115e No 0.54 0.33 -45 1.9E+02 21 -0.028
115n No 0.53 0.33 -45 1.9E+02 14 0.016
116e No 0.56 0.35 -45 1.9E+02 11 0.023
116n No 0.56 0.35 -45 1.9E+02 15 -0.024
117e No 0.58 0.35 -45 1.9E+02 19 -0.017
117n No 0.59 0.35 -45 1.9E+02 18 -0.094
118e No 0.6 0.36 -45 1.9E+02 21 -0.041
118n No 0.59 0.36 -45 1.9E+02 43 -0.15
119e No 0.6 0.34 -45 1.9E+02 28 -0.13
119n No 0.6 0.34 -45 1.9E+02 26 -0.14
120e No 0.61 0.35 -45 1.9E+02 30 -0.17
120n No 0.59 0.35 -45 1.9E+02 50 -0.18
121e No 0.6 0.34 -45 1.9E+02 14 -0.016
121n No 0.61 0.34 -45 1.9E+02 17 -0.19
122e No 0.61 0.34 -45 1.9E+02 18 -0.15
122n No 0.6 0.34 -45 1.9E+02 18 -0.16
123e No 0.61 0.33 -45 1.9E+02 17 -0.13
123n No 0.6 0.33 -45 1.9E+02 18 -0.19
124e No 0.61 0.34 -45 1.9E+02 16 -0.012
124n No 0.6 0.34 -45 1.9E+02 18 0.0017
125e No 0.6 0.34 -45 1.9E+02 22 -0.13
125n No 0.59 0.34 -45 1.9E+02 24 -0.27
126e No 0.6 0.34 -45 1.9E+02 17 -0.12
126n No 0.6 0.34 -45 1.9E+02 18 -0.22
127e No 0.6 0.34 -45 1.9E+02 14 -0.056
127n No 0.6 0.34 -45 1.9E+02 15 -0.11
128e No 0.61 0.35 -45 1.9E+02 19 -0.066
128n No 0.59 0.35 -45 1.9E+02 9 -0.01
129e No 0.61 0.34 -45 1.9E+02 18 -0.052
129n No 0.61 0.34 -45 1.9E+02 18 -0.13
130e No 0.6 0.34 -45 1.9E+02 25 -0.081
130n No 0.59 0.34 -45 1.9E+02 20 -0.048
131e No 0.57 0.33 -45 1.9E+02 15 0.031
131n No 0.56 0.33 -45 1.9E+02 14 0.02
132e No 0.57 0.34 -45 1.9E+02 22 -0.021
132n No 0.57 0.34 -45 1.9E+02 19 -0.042
133e No 0.56 0.34 -45 1.9E+02 20 -0.026
133n No 0.56 0.34 -45 1.9E+02 19 -0.035
134e No 0.56 0.35 -45 1.9E+02 33 -0.087
134n No 0.54 0.35 -45 1.9E+02 19 0.0025
135e No 0.052 0.41 -45 1.9E+02 7.3 2.2
135n No 0.57 0.41 -45 1.9E+02 18 -0.071
136e No 0.59 0.36 -45 1.9E+02 21 -0.13
136n No 0.57 0.36 -45 1.9E+02 17 0.17
137e No 0.58 0.34 -45 1.9E+02 16 -0.016
137n No 0.58 0.34 -45 1.9E+02 35 -0.13
138e No 0.59 0.34 -45 1.9E+02 22 -0.093
138n No 0.59 0.34 -45 1.9E+02 18 -0.07
139e No 0.59 0.33 -45 1.9E+02 33 -0.076
139n No 0.58 0.33 -45 1.9E+02 28 -0.08
140e No 0.6 0.34 -45 1.9E+02 18 -0.081
140n No 0.59 0.34 -45 1.9E+02 20 0.039
141e No 0.61 0.34 -45 1.9E+02 20 -0.081
141n No 0.6 0.34 -45 1.9E+02 20 -0.12
142e No 0.6 0.34 -45 1.9E+02 22 0.001
142n No 0.56 0.34 -45 1.9E+02 14 1
143e No 0.61 0.34 -45 1.9E+02 27 -0.085
143n No 0.61 0.34 -45 1.9E+02 30 -0.11
144e No 0.61 0.34 -45 1.9E+02 17 -0.091
144n No 0.61 0.34 -45 1.9E+02 15 -0.13
145e No 0.62 0.35 -45 1.9E+02 19 -0.12
145n No 0.61 0.35 -45 1.9E+02 22 -0.099
146e No 0.58 0.32 -45 1.9E+02 21 0.049
146n No 0.57 0.32 -45 1.9E+02 16 0.0066
147e No 0.61 0.34 -45 1.9E+02 17 -0.044
147n No 0.61 0.34 -45 1.9E+02 19 -0.11
148e No 0.59 0.33 -45 1.9E+02 23 0.01
148n No 0.59 0.33 -45 1.9E+02 21 -0.067
149e No 0.61 0.34 -45 1.9E+02 20 -0.047
149n No 0.6 0.34 -45 1.9E+02 21 -0.13
150e No 0.61 0.35 -45 1.9E+02 22 -0.1
150n No 0.61 0.35 -45 1.9E+02 18 -0.1
151e Yes -45 1.5E+03 0 0
151n Yes -45 1.5E+03 0 0
152e No 0.55 0.32 -45 1.9E+02 18 0.0073
152n No 0.56 0.32 -45 1.9E+02 19 -0.042
153e No 0.56 0.33 -45 1.9E+02 18 -0.0048
153n No 0.55 0.33 -45 1.9E+02 19 0.038
154e No 0.54 0.34 -45 1.9E+02 20 -0.026
154n No 0.55 0.34 -45 1.9E+02 20 -0.053
155e No 0.57 0.35 -45 1.9E+02 18 -0.066
155n No 0.57 0.35 -45 1.9E+02 21 -0.12
156e No 0.31 0.11 -45 1.9E+02 16 0.54
156n No 0.27 0.11 -45 1.9E+02 12 0.83
157e No 0.58 0.34 -45 1.9E+02 18 -0.089
157n No 0.58 0.34 -45 1.9E+02 19 -0.099
158e No 0.6 0.34 -45 1.9E+02 18 0.029
158n No 0.59 0.34 -45 1.9E+02 19 -0.076
159e No 0.52 0.32 -45 1.9E+02 10 0.093
159n No 0.56 0.32 -45 1.9E+02 21 -0.0074
160e No 0.6 0.33 -45 1.9E+02 18 -0.1
160n No 0.6 0.33 -45 1.9E+02 19 -0.091
161e No 0.6 0.32 -45 1.9E+02 21 0.04
161n No 0.46 0.32 -45 1.9E+02 11 0.95
162e No 0.61 0.33 -45 1.9E+02 21 -0.0025
162n No 0.6 0.33 -45 1.9E+02 9.8 -0.038
163e No 0.61 0.34 -45 1.9E+02 11 -0.11
163n No 0.6 0.34 -45 1.9E+02 20 -0.14
164e No 0.61 0.34 -45 1.9E+02 13 -0.021
164n No 0.61 0.34 -45 1.9E+02 14 -0.13
165e No 0.61 0.34 -45 1.9E+02 13 -0.13
165n No 0.6 0.34 -45 1.9E+02 9.2 -0.15
166e No 0.61 0.33 -45 1.9E+02 26 -0.12
166n No 0.6 0.33 -45 1.9E+02 32 -0.17
167e No 0.62 0.35 -45 1.9E+02 18 -0.084
167n No 0.61 0.35 -45 1.9E+02 43 -0.15
168e No 0.61 0.35 -45 1.9E+02 26 -0.07
168n No 0.61 0.35 -45 1.9E+02 20 -0.12
169e No 0.61 0.34 -45 1.9E+02 14 -0.055
169n No 0.61 0.34 -45 1.9E+02 17 -0.091
170e No 0.052 0.36 -45 1.9E+02 0.58 0.34
170n No 0.54 0.36 -45 1.9E+02 35 -0.097
171e Yes -45 1.5E+03 0 0
171n Yes -45 1.5E+03 0 0
172e Yes -45 1.5E+03 0 0
172n Yes -45 1.5E+03 0 0
173e No 0.55 0.33 -45 1.9E+02 17 0.012
173n No 0.55 0.33 -45 1.9E+02 19 0.034
174e No 0.55 0.34 -45 1.9E+02 21 -0.0046
174n No 0.56 0.34 -45 1.9E+02 34 -0.12
175e No 0.53 0.34 -45 1.9E+02 17 0.035
175n No 0.49 0.34 -45 1.9E+02 8.8 0.091
176e No 0.56 0.35 -45 1.9E+02 48 -0.16
176n No 0.57 0.35 -45 1.9E+02 16 -0.12
177e No 0.59 0.35 -45 1.9E+02 13 -0.099
177n No 0.59 0.35 -45 1.9E+02 18 -0.14
178e No 0.59 0.34 -45 1.9E+02 21 -0.049
178n No 0.59 0.34 -45 1.9E+02 42 -0.19
179e No 0.6 0.35 -45 1.9E+02 18 -0.0057
179n No 0.6 0.35 -45 1.9E+02 20 -0.063
180e No 0.61 0.34 -45 1.9E+02 20 -0.0044
180n No 0.45 0.34 -45 1.9E+02 14 0.76
181e No 0.61 0.34 -45 1.9E+02 20 -0.068
181n No 0.6 0.34 -45 1.9E+02 19 -0.037
182e No 0.61 0.34 -45 1.9E+02 14 0.038
182n No 0.61 0.34 -45 1.9E+02 16 -0.12
183e No 0.61 0.34 -45 1.9E+02 27 -0.095
183n No 0.6 0.34 -45 1.9E+02 36 -0.12
184e No 0.6 0.33 -45 1.9E+02 9.6 0.05
184n No 0.6 0.33 -45 1.9E+02 18 -0.04
185e No 0.61 0.33 -45 1.9E+02 13 0.0029
185n No 0.6 0.33 -45 1.9E+02 15 -0.041
186e No 0.6 0.33 -45 1.9E+02 33 -0.1
186n No 0.59 0.33 -45 1.9E+02 22 -0.12
187e No 0.61 0.34 -45 1.9E+02 24 -0.083
187n No 0.6 0.34 -45 1.9E+02 15 -0.045
188e No 0.59 0.34 -45 1.9E+02 20 -0.041
188n No 0.55 0.34 -45 1.9E+02 37 0.15
189e No 0.59 0.33 -45 1.9E+02 19 -0.033
189n No 0.59 0.33 -45 1.9E+02 24 -0.14
190e No 0.61 0.35 -45 1.9E+02 19 -0.019
190n No 0.6 0.35 -45 1.9E+02 18 -0.075
191e No 0.61 0.35 -45 1.9E+02 21 -0.076
191n No 0.6 0.35 -45 1.9E+02 30 -0.036
192e No 0.57 0.34 -45 1.9E+02 23 -0.0082
192n No 0.58 0.34 -45 1.9E+02 29 -0.085
193e No 0.56 0.33 -45 1.9E+02 21 0.047
193n No 0.55 0.33 -45 1.9E+02 18 0.001
194e No 0.54 0.32 -45 1.9E+02 35 -0.12
194n No 0.53 0.32 -45 1.9E+02 29 -0.11
195e No 0.55 0.34 -45 1.9E+02 24 -0.047
195n No 0.54 0.34 -45 1.9E+02 24 -0.037
196e No 0.039 0.42 -45 1.9E+02 2.9 0.36
196n No 0.56 0.42 -45 1.9E+02 21 -0.079
197e No 0.48 0.26 -45 1.9E+02 15 0.67
197n No 0.48 0.26 -45 1.9E+02 19 0.59
198e No 0.58 0.34 -45 1.9E+02 30 -0.094
198n No 0.57 0.34 -45 1.9E+02 36 -0.15
200e No 0.048 0.44 -45 1.9E+02 2.8 0.36
200n No 0.58 0.44 -45 1.9E+02 28 -0.086
201e No 0.59 0.34 -45 1.9E+02 55 -0.13
201n No 0.59 0.34 -45 1.9E+02 45 -0.15
202e No 0.6 0.34 -45 1.9E+02 34 -0.046
202n No 0.59 0.34 -45 1.9E+02 21 -0.027
203e No 0.59 0.33 -45 1.9E+02 19 -0.03
203n No 0.59 0.33 -45 1.9E+02 18 -0.13
204e No 0.6 0.33 -45 1.9E+02 12 -0.32
204n No 0.59 0.33 -45 1.9E+02 15 -0.46
205e No 0.59 0.33 -45 1.9E+02 25 -0.07
205n No 0.59 0.33 -45 1.9E+02 24 -0.1
206e No 0.59 0.34 -45 1.9E+02 33 -0.072
206n No 0.58 0.34 -45 1.9E+02 22 -0.094
207e No 0.59 0.42 -45 1.9E+02 26 -0.005
207n No 0.046 0.42 -45 1.9E+02 2.7 0.38
208e No 0.59 0.34 -45 1.9E+02 33 -0.13
208n No 0.57 0.34 -45 1.9E+02 44 -0.024
209e No 0.57 0.33 -45 1.9E+02 28 -0.064
209n No 0.56 0.33 -45 1.9E+02 32 -0.0047
210e No 0.6 0.35 -45 1.9E+02 23 -0.36
210n No 0.57 0.35 -45 1.9E+02 16 -0.16
211e No 0.58 0.35 -45 1.9E+02 31 -0.046
211n No 0.58 0.35 -45 1.9E+02 22 -0.05
212e No 0.56 0.34 -45 1.9E+02 27 -0.029
212n No 0.53 0.34 -45 1.9E+02 14 -0.014
213e No 0.54 0.34 -45 1.9E+02 19 0.015
213n No 0.56 0.34 -45 1.9E+02 21 -0.079
214e No 0.55 0.35 -45 1.9E+02 24 -0.0091
214n No 0.55 0.35 -45 1.9E+02 21 -0.034
215e No 0.56 0.35 -45 1.9E+02 37 -0.094
215n No 0.54 0.35 -45 1.9E+02 18 0.091
216e No 0.57 0.34 -45 1.9E+02 22 -0.0033
216n No 0.55 0.34 -45 1.9E+02 15 0.068
217e No 0.57 0.34 -45 1.9E+02 28 0.0046
217n No 0.56 0.34 -45 1.9E+02 20 -0.018
218e No 0.57 0.41 -45 1.9E+02 31 -0.064
218n No 0.21 0.41 -45 1.9E+02 0.65 0.14
219e No 0.58 0.34 -45 1.9E+02 23 -0.07
219n No 0.57 0.34 -45 1.9E+02 19 -0.09
220e No 0.58 0.33 -45 1.9E+02 21 -0.00068
220n No 0.57 0.33 -45 1.9E+02 19 -0.063
221e No 0.58 0.33 -45 1.9E+02 22 0.012
221n No 0.57 0.33 -45 1.9E+02 18 0.019
222e No 0.58 0.33 -45 1.9E+02 19 0.057
222n No 0.59 0.33 -45 1.9E+02 28 -0.084
223e No 0.57 0.31 -45 1.9E+02 19 -0.033
223n No 0.57 0.31 -45 1.9E+02 18 -0.091
224e No 0.58 0.33 -45 1.9E+02 25 0.0013
224n No 0.58 0.33 -45 1.9E+02 23 -0.064
225e No 0.59 0.34 -45 1.9E+02 27 -0.034
225n No 0.58 0.34 -45 1.9E+02 23 -0.087
226e No 0.59 0.34 -45 1.9E+02 26 -0.09
226n No 0.52 0.34 -45 1.9E+02 31 0.32
227e No 0.58 0.34 -45 1.9E+02 21 -0.026
227n No 0.56 0.34 -45 1.9E+02 16 -0.00085
228e No 0.56 0.32 -45 1.9E+02 20 0.039
228n No 0.56 0.32 -45 1.9E+02 19 -0.029
229e No 0.56 0.33 -45 1.9E+02 19 0.036
229n No 0.56 0.33 -45 1.9E+02 19 -0.063
231e No 0.56 0.34 -45 1.9E+02 40 -0.13
231n No 0.56 0.34 -45 1.9E+02 33 -0.15
232e No 0.55 0.35 -45 1.9E+02 49 -0.14
232n No 0.53 0.35 -45 1.9E+02 16 -0.013
233e No 0.55 0.33 -45 1.9E+02 18 -0.025
233n No 0.47 0.33 -45 1.9E+02 31 0.46
234e No 0.56 0.35 -45 1.9E+02 30 -0.015
234n No 0.5 0.35 -45 1.9E+02 8.5 0.11
237e No 0.56 0.32 -45 1.9E+02 19 0.0096
237n No 0.57 0.32 -45 1.9E+02 22 -0.022
238e No 0.58 0.33 -45 1.9E+02 38 -0.093
238n No 0.57 0.33 -45 1.9E+02 24 -0.0055
239e No 0.58 0.33 -45 1.9E+02 24 -0.037
239n No 0.58 0.33 -45 1.9E+02 31 -0.11
240e No 0.53 0.27 -45 1.9E+02 20 0.14
240n No 0.49 0.27 -45 1.9E+02 26 0.8
241e No 0.58 0.34 -45 1.9E+02 21 -0.06
241n No 0.57 0.34 -45 1.9E+02 19 -0.065
242e No 0.58 0.34 -45 1.9E+02 27 -0.098
242n No 0.57 0.34 -45 1.9E+02 26 -0.095
243e No 0.58 0.34 -45 1.9E+02 36 -0.043
243n No 0.57 0.34 -45 1.9E+02 20 -0.057
244e No 0.57 0.33 -45 1.9E+02 19 -0.057
244n No 0.57 0.33 -45 1.9E+02 19 -0.052
245e No 0.58 0.34 -45 1.9E+02 29 -0.04
245n No 0.57 0.34 -45 1.9E+02 22 0.016
246e No 0.56 0.33 -45 1.9E+02 20 0.13
246n No 0.56 0.33 -45 1.9E+02 28 -0.035
250e No 0.54 0.33 -45 1.9E+02 26 -0.046
250n No 0.54 0.33 -45 1.9E+02 21 -0.027
251e No 0.074 0.097 -45 1.9E+02 5 0.44
251n No 0.22 0.097 -45 1.9E+02 23 -0.061
252e No 0.55 0.32 -45 1.9E+02 19 0.028
252n No 0.54 0.32 -45 1.9E+02 19 -0.046
253e No 0.55 0.33 -45 1.9E+02 17 0.042
253n No 0.53 0.33 -45 1.9E+02 12 0.069
254e No 0.58 0.33 -45 1.9E+02 36 -0.066
254n No 0.57 0.33 -45 1.9E+02 31 -0.061
255e No 0.54 0.33 -45 1.9E+02 11 0.15
255n No 0.56 0.33 -45 1.9E+02 21 -0.0067
256e No 0.58 0.34 -45 1.9E+02 28 -0.054
256n No 0.54 0.34 -45 1.9E+02 13 0.026
257e No 0.52 0.33 -45 1.9E+02 35 0.15
257n No 0.57 0.33 -45 1.9E+02 28 -0.16
261e No 0.54 0.32 -45 1.9E+02 12 0.045
261n No 0.52 0.32 -45 1.9E+02 9.2 0.055
262e No 0.58 0.34 -45 1.9E+02 15 -0.26
262n No 0.58 0.34 -45 1.9E+02 20 -0.43
266e No 0.56 0.34 -45 1.9E+02 10 -0.33
266n No 0.56 0.34 -45 1.9E+02 15 -0.47
267e No 0.53 0.33 -45 1.9E+02 14 -0.033
267n No 0.55 0.33 -45 1.9E+02 16 -0.085
268e No 0.53 0.34 -45 1.9E+02 11 0.067
268n No 0.56 0.34 -45 1.9E+02 21 -0.048
269e No 0.56 0.35 -45 1.9E+02 22 -0.03
269n No 0.48 0.35 -45 1.9E+02 7.3 0.092
270e No 0.57 0.33 -45 1.9E+02 27 -0.006
270n No 0.57 0.33 -45 1.9E+02 33 -0.077
271e No 0.048 0.46 -45 1.9E+02 3 0.34
271n No 0.55 0.46 -45 1.9E+02 20 -0.075
272e No 0.54 0.32 -45 1.9E+02 19 0.097
272n No 0.55 0.32 -45 1.9E+02 28 -0.049
273e No 0.22 -0.3 -45 1.9E+02 79 -0.2
273n No 0.24 -0.3 -45 1.9E+02 31 -0.12
277e No 0.42 0.34 -45 1.9E+02 11 0.67
277n No 0.57 0.34 -45 1.9E+02 27 -0.058
278e No 0.29 0.35 -45 1.9E+02 4.4 0.31
278n No 0.49 0.35 -45 1.9E+02 8.4 0.089
281e No 0.54 0.33 -45 1.9E+02 19 0.02
281n No 0.54 0.33 -45 1.9E+02 18 -0.03
282e No 0.55 0.34 -45 1.9E+02 21 -0.0067
282n No 0.56 0.34 -45 1.9E+02 37 -0.13
283e No 0.56 0.34 -45 1.9E+02 24 -0.08
283n No 0.54 0.34 -45 1.9E+02 16 -0.033
284e No 0.57 0.34 -45 1.9E+02 30 -0.061
284n No 0.56 0.34 -45 1.9E+02 35 -0.12
285e No 0.052 0.0028 -45 1.9E+02 2.9 0.37
285n No 0.055 0.0028 -45 1.9E+02 2.5 0.41
286e No 0.23 -0.3 -45 1.9E+02 33 -0.05
286n No 0.22 -0.3 -45 1.9E+02 78 -0.2
287e No 0.55 0.3 -45 1.9E+02 43 -0.1
287n No 0.25 0.3 -45 1.9E+02 29 -0.17
290e No 0.51 0.3 -45 46 76 -0.18
290n No 0.5 0.3 -45 1.9E+02 75 -0.19
291e No 0.53 0.43 -45 1.9E+02 72 -0.18
291n No 0.07 0.43 -45 1.9E+02 0.65 -0.25
292e No 0.55 0.34 -45 1.9E+02 19 0.0057
292n No 0.55 0.34 -45 1.9E+02 39 -0.11
293e No 0.42 0.34 -45 1.9E+02 23 0.66
293n No 0.54 0.34 -45 1.9E+02 16 -0.0083
294e No 0.47 0.34 -45 1.9E+02 24 0.37
294n No 0.54 0.34 -45 1.9E+02 23 -0.12
295e No 0.55 0.44 -45 1.9E+02 5.3 -0.47
295n No 0.047 0.44 -45 1.9E+02 0.8 -0.14
299e No 0.55 0.34 -45 1.9E+02 26 -0.004
299n No 0.53 0.34 -45 1.9E+02 18 -0.051
300e No 0.56 0.34 -45 1.9E+02 28 -0.095
300n No 0.42 0.34 -45 1.9E+02 22 0.69
301e No 0.56 0.35 -45 1.9E+02 32 -0.094
301n No 0.52 0.35 -45 1.9E+02 14 -0.039
302e No 0.43 0.37 -45 1.9E+02 5.7 0.19
302n No 0.55 0.37 -45 1.9E+02 22 -0.076
306e No 0.54 0.34 -45 1.9E+02 15 -0.015
306n No 0.42 0.34 -45 1.9E+02 18 0.61
307e No 0.4 0.33 -45 1.9E+02 17 0.67
307n No 0.53 0.33 -45 1.9E+02 19 -0.088
311e No 0.51 0.35 -45 1.9E+02 17 0.049
311n No 0.27 0.35 -45 1.9E+02 4.3 0.57
312e No 0.21 -0.15 -45 1.9E+02 17 0.71
312n No 0.19 -0.15 -45 1.9E+02 8 0.07
313e No 0.4 0.34 -45 1.9E+02 11 0.69
313n No 0.54 0.34 -45 1.9E+02 46 -0.14
314e No 0.55 0.34 -45 1.9E+02 41 -0.088
314n No 0.54 0.34 -45 1.9E+02 36 -0.091
315e No 0.53 0.33 -45 1.9E+02 21 -0.048
315n No 0.51 0.33 -45 1.9E+02 18 0.052
316e No 0.53 0.34 -45 1.9E+02 22 -0.025
316n No 0.54 0.34 -45 1.9E+02 25 -0.12
317e No 0.55 0.35 -45 1.9E+02 32 -0.13
317n No 0.53 0.35 -45 1.9E+02 15 -0.052
318e No 0.53 0.35 -45 1.9E+02 65 -0.21
318n No 0.51 0.35 -45 1.9E+02 79 -0.22
319e No 0.55 0.36 -45 1.9E+02 26 -0.078
319n No 0.54 0.36 -45 1.9E+02 31 -0.1
320e No 0.061 0.0055 -45 1.9E+02 1.4 0.45
320n No 0.056 0.0055 -45 1.9E+02 1.4 0.37
321e No 0.43 0.33 -45 1.9E+02 8.3 -0.09
321n No 0.47 0.33 -45 1.9E+02 12 -0.22
322e No 0.04 0.005 -45 1.9E+02 3 0.34
322n No 0.033 0.005 -45 1.9E+02 2.8 0.38
323e No 0.045 -0.0011 -45 1.9E+02 2.8 0.39
323n No 0.041 -0.0011 -45 1.9E+02 2.8 0.39
324e No 0.46 0.31 -45 1.9E+02 30 -0.076
324n No 0.47 0.31 -45 1.9E+02 32 -0.14
325e No 0.49 0.33 -45 1.9E+02 28 -0.071
325n No 0.49 0.33 -45 1.9E+02 22 -0.076
326e No 0.045 0.0014 -45 1.9E+02 3.1 0.35
326n No 0.049 0.0014 -45 1.9E+02 2.8 0.36
327e No 0.47 0.3 -45 1.9E+02 31 -0.062
327n No 0.48 0.3 -45 1.9E+02 32 -0.13
328e No 0.054 0.36 -45 1.9E+02 3 0.34
328n No 0.48 0.36 -45 1.9E+02 22 -0.1
329e No 0.36 0.24 -45 1.9E+02 16 0.72
329n No 0.049 0.24 -45 1.9E+02 2.8 0.38
331e No 0.37 0.23 -45 1.9E+02 23 0.71
331n No 0.41 0.23 -45 1.9E+02 18 0.64
332e No 0.45 0.3 -45 1.9E+02 28 -0.085
332n No 0.45 0.3 -45 1.9E+02 15 0.09
333e No 0.43 0.32 -45 1.9E+02 9.3 0.02
333n No 0.47 0.32 -45 1.9E+02 14 -0.16
336e No 0.15 -0.28 -45 1.9E+02 52 -0.14
336n No 0.16 -0.28 -45 1.9E+02 20 -0.0006
339e Yes -45 1.5E+03 0 0
339n Yes -45 1.5E+03 0 0
340e No 0.48 0.33 -45 1.9E+02 31 -0.2
340n No 0.49 0.33 -45 1.9E+02 22 -0.22
342e No 0.48 0.34 -45 1.9E+02 16 0.083
342n No 0.51 0.34 -45 1.9E+02 28 -0.071
343e No 0.05 0.4 -45 1.9E+02 3.1 0.37
343n No 0.5 0.4 -45 1.9E+02 15 0.022
344e No 0.038 0.0053 -45 1.9E+02 2.8 0.37
344n No 0.042 0.0053 -45 1.9E+02 2.5 0.38
345e Yes -45 1.5E+03 0 0
345n Yes -45 1.5E+03 0 0
346e No 0.33 0.23 -45 1.9E+02 5 0.22
346n No 0.17 0.23 -45 1.9E+02 2.7 0.34
347e No 0.034 0.0037 -45 1.9E+02 2.9 0.33
347n No 0.038 0.0037 -45 1.9E+02 2.5 0.36
348e No 0.047 0.09 -45 1.9E+02 3.1 0.4
348n No 0.15 0.09 -45 1.9E+02 2.9 0.36
349e No 0.47 0.33 -45 1.9E+02 59 -0.18
349n No 0.45 0.33 -45 1.9E+02 75 -0.21
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/2461081/zen.2461081.45950.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.dev97+gc2668d3f7
hera_filters: 0.1.7
hera_qm: 2.2.1.dev4+gf6d02113b
hera_notebook_templates: 0.0.1.dev1313+g92178a09c
In [66]:
print(f'Finished execution in {(time.time() - tstart) / 60:.2f} minutes.')
Finished execution in 1.98 minutes.