Calibration Smoothing¶

by Josh Dillon, last updated March 29, 2023

This notebook runs calibration smoothing to the gains coming out of file_calibration notebook. It removes any flags founds on by that notebook and replaces them with flags generated from full_day_rfi and full_day_antenna_flagging. It also plots the results for a couple of antennas.

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

• Figure 1: Full-Day Gain Amplitudes Before and After smooth_cal¶

• Figure 2: Full-Day Gain Phases Before and After smooth_cal¶

• Figure 3: Full-Day $\chi^2$ / DoF Waterfall from Redundant-Baseline Calibration¶

• Figure 4: Average $\chi^2$ per Antenna vs. Time and Frequency¶

In [1]:
import time
tstart = time.time()
In [2]:
import os
os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE'
import h5py
import hdf5plugin  # REQUIRED to have the compression plugins available
import numpy as np
import glob
import copy
import warnings
import matplotlib
import matplotlib.pyplot as plt
from hera_cal import io, utils, smooth_cal
from hera_qm.time_series_metrics import true_stretches
%matplotlib inline
from IPython.display import display, HTML

Parse inputs¶

In [3]:
# get files
SUM_FILE = os.environ.get("SUM_FILE", None)
# SUM_FILE = '/users/jsdillon/lustre/H6C/abscal/2459853/zen.2459853.25518.sum.uvh5'
SUM_SUFFIX = os.environ.get("SUM_SUFFIX", 'sum.uvh5')
CAL_SUFFIX = os.environ.get("CAL_SUFFIX", 'sum.omni.calfits')
SMOOTH_CAL_SUFFIX = os.environ.get("SMOOTH_CAL_SUFFIX", 'sum.smooth.calfits')
ANT_FLAG_SUFFIX = os.environ.get("ANT_FLAG_SUFFIX", 'sum.antenna_flags.h5')
RFI_FLAG_SUFFIX = os.environ.get("RFI_FLAG_SUFFIX", 'sum.flag_waterfall.h5')
FREQ_SMOOTHING_SCALE = float(os.environ.get("FREQ_SMOOTHING_SCALE", 10.0)) # MHz
TIME_SMOOTHING_SCALE = float(os.environ.get("TIME_SMOOTHING_SCALE", 6e5)) # seconds
EIGENVAL_CUTOFF = float(os.environ.get("EIGENVAL_CUTOFF", 1e-12))

for setting in ['SUM_FILE', 'SUM_SUFFIX', 'CAL_SUFFIX', 'SMOOTH_CAL_SUFFIX', 'ANT_FLAG_SUFFIX',
                'RFI_FLAG_SUFFIX', 'FREQ_SMOOTHING_SCALE', 'TIME_SMOOTHING_SCALE', 'EIGENVAL_CUTOFF']:
        print(f'{setting} = {eval(setting)}')
SUM_FILE = /mnt/sn1/data1/2460449/zen.2460449.16885.sum.uvh5
SUM_SUFFIX = sum.uvh5
CAL_SUFFIX = sum.omni.calfits
SMOOTH_CAL_SUFFIX = sum.smooth.calfits
ANT_FLAG_SUFFIX = sum.antenna_flags.h5
RFI_FLAG_SUFFIX = sum.flag_waterfall.h5
FREQ_SMOOTHING_SCALE = 10.0
TIME_SMOOTHING_SCALE = 600000.0
EIGENVAL_CUTOFF = 1e-12

Load files¶

In [4]:
sum_glob = '.'.join(SUM_FILE.split('.')[:-3]) + '.*.' + SUM_SUFFIX
cal_files_glob = sum_glob.replace(SUM_SUFFIX, CAL_SUFFIX)
cal_files = sorted(glob.glob(cal_files_glob))
print(f'Found {len(cal_files)} *.{CAL_SUFFIX} files starting with {cal_files[0]}.')
Found 1572 *.sum.omni.calfits files starting with /mnt/sn1/data1/2460449/zen.2460449.16885.sum.omni.calfits.
In [5]:
rfi_flag_files_glob = sum_glob.replace(SUM_SUFFIX, RFI_FLAG_SUFFIX)
rfi_flag_files = sorted(glob.glob(rfi_flag_files_glob))
print(f'Found {len(rfi_flag_files)} *.{RFI_FLAG_SUFFIX} files starting with {rfi_flag_files[0]}.')
Found 1572 *.sum.flag_waterfall.h5 files starting with /mnt/sn1/data1/2460449/zen.2460449.16885.sum.flag_waterfall.h5.
In [6]:
ant_flag_files_glob = sum_glob.replace(SUM_SUFFIX, ANT_FLAG_SUFFIX)
ant_flag_files = sorted(glob.glob(ant_flag_files_glob))
print(f'Found {len(ant_flag_files)} *.{ANT_FLAG_SUFFIX} files starting with {ant_flag_files[0]}.')
Found 1572 *.sum.antenna_flags.h5 files starting with /mnt/sn1/data1/2460449/zen.2460449.16885.sum.antenna_flags.h5.
In [7]:
cs = smooth_cal.CalibrationSmoother(cal_files, flag_file_list=(ant_flag_files + rfi_flag_files), ignore_calflags=True,
                                    pick_refant=True, propagate_refant_flags=True, load_chisq=True, load_cspa=True)
for pol in cs.refant:
    print(f'Reference antenna {cs.refant[pol][0]} selected for {pol}.')
Mean of empty slice
Reference antenna 222 selected for Jee.
Reference antenna 94 selected for Jnn.
In [8]:
# duplicate a small number of abscal gains for plotting
antnums = set([ant[0] for ant in cs.ants])
flags_per_antnum = [np.sum(cs.flag_grids[ant, 'Jnn']) + np.sum(cs.flag_grids[ant, 'Jee']) for ant in antnums]
refant_nums = [ant[0] for ant in cs.refant.values()]
candidate_ants = [ant for ant, nflags in zip(antnums, flags_per_antnum) if (ant not in refant_nums) and (nflags <= np.percentile(flags_per_antnum, 25))
                  and not np.all(cs.flag_grids[ant, 'Jee']) and not np.all(cs.flag_grids[ant, 'Jnn'])]
ants_to_plot = [func(candidate_ants) for func in (np.min, np.max)]
abscal_gains = {(ant, pol): np.array(cs.gain_grids[(ant, pol)]) for ant in ants_to_plot for pol in ['Jee', 'Jnn']}

Perform smoothing¶

In [9]:
cs.time_freq_2D_filter(freq_scale=FREQ_SMOOTHING_SCALE, time_scale=TIME_SMOOTHING_SCALE, eigenval_cutoff=EIGENVAL_CUTOFF, 
                       method='DPSS', fit_method='lu_solve', fix_phase_flips=True, flag_phase_flip_ints=True)
22 phase flips detected on antenna (240, 'Jnn'). A total of 100 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5183729506.
77 phase flips detected on antenna (283, 'Jnn'). A total of 124 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5203862167.
14 phase flips detected on antenna (142, 'Jnn'). A total of 82 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5116620637.
2 phase flips detected on antenna (189, 'Jnn'). A total of 48 integrations were phase-flipped relative to the 0th integration between 2460449.503161607 and 2460449.5084184683.
20 phase flips detected on antenna (222, 'Jnn'). A total of 40 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5116620637.
70 phase flips detected on antenna (295, 'Jnn'). A total of 107 integrations were phase-flipped relative to the 0th integration between 2460449.490634618 and 2460449.519938824.
23 phase flips detected on antenna (239, 'Jnn'). A total of 97 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
3 phase flips detected on antenna (224, 'Jnn'). A total of 113 integrations were phase-flipped relative to the 0th integration between 2460449.503161607 and 2460449.5203862167.
2 phase flips detected on antenna (124, 'Jnn'). A total of 1 integrations were phase-flipped relative to the 0th integration between 2460449.503161607 and 2460449.503161607.
65 phase flips detected on antenna (285, 'Jnn'). A total of 111 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5203862167.
22 phase flips detected on antenna (85, 'Jnn'). A total of 43 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.511773912.
25 phase flips detected on antenna (51, 'Jnn'). A total of 109 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
12 phase flips detected on antenna (123, 'Jnn'). A total of 9 integrations were phase-flipped relative to the 0th integration between 2460449.5017075813 and 2460449.5084184683.
21 phase flips detected on antenna (53, 'Jnn'). A total of 122 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
22 phase flips detected on antenna (221, 'Jnn'). A total of 70 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5116620637.
4 phase flips detected on antenna (87, 'Jnn'). A total of 5 integrations were phase-flipped relative to the 0th integration between 2460449.502490518 and 2460449.5030497587.
23 phase flips detected on antenna (52, 'Jnn'). A total of 149 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
20 phase flips detected on antenna (141, 'Jnn'). A total of 75 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5184847987.
27 phase flips detected on antenna (208, 'Jnn'). A total of 80 integrations were phase-flipped relative to the 0th integration between 2460449.5017075813 and 2460449.5203862167.
26 phase flips detected on antenna (101, 'Jnn'). A total of 55 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.516807077.
28 phase flips detected on antenna (50, 'Jnn'). A total of 67 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.516807077.
31 phase flips detected on antenna (228, 'Jnn'). A total of 115 integrations were phase-flipped relative to the 0th integration between 2460449.50237867 and 2460449.5203862167.
28 phase flips detected on antenna (66, 'Jnn'). A total of 52 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.516807077.
27 phase flips detected on antenna (179, 'Jnn'). A total of 89 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
2 phase flips detected on antenna (192, 'Jnn'). A total of 63 integrations were phase-flipped relative to the 0th integration between 2460449.5038326955 and 2460449.5107672787.
7 phase flips detected on antenna (20, 'Jnn'). A total of 79 integrations were phase-flipped relative to the 0th integration between 2460449.509313253 and 2460449.5203862167.
22 phase flips detected on antenna (181, 'Jnn'). A total of 72 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5116620637.
43 phase flips detected on antenna (65, 'Jnn'). A total of 135 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5203862167.
16 phase flips detected on antenna (122, 'Jnn'). A total of 78 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5116620637.
20 phase flips detected on antenna (38, 'Jnn'). A total of 43 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5084184683.
9 phase flips detected on antenna (184, 'Jnn'). A total of 153 integrations were phase-flipped relative to the 0th integration between 2460449.502490518 and 2460449.5203862167.
1 phase flips detected on antenna (191, 'Jnn'). A total of 149 integrations were phase-flipped relative to the 0th integration between 2460449.5038326955 and 2460449.5203862167.
22 phase flips detected on antenna (102, 'Jnn'). A total of 12 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.511773912.
14 phase flips detected on antenna (143, 'Jnn'). A total of 74 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5107672787.
25 phase flips detected on antenna (187, 'Jnn'). A total of 54 integrations were phase-flipped relative to the 0th integration between 2460449.5017075813 and 2460449.5203862167.
3 phase flips detected on antenna (185, 'Jnn'). A total of 115 integrations were phase-flipped relative to the 0th integration between 2460449.5028260625 and 2460449.5203862167.
61 phase flips detected on antenna (216, 'Jnn'). A total of 135 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5203862167.
27 phase flips detected on antenna (140, 'Jnn'). A total of 105 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
1 phase flips detected on antenna (56, 'Jnn'). A total of 158 integrations were phase-flipped relative to the 0th integration between 2460449.5028260625 and 2460449.5203862167.
4 phase flips detected on antenna (211, 'Jnn'). A total of 70 integrations were phase-flipped relative to the 0th integration between 2460449.5029379106 and 2460449.5107672787.
4 phase flips detected on antenna (167, 'Jnn'). A total of 4 integrations were phase-flipped relative to the 0th integration between 2460449.5028260625 and 2460449.5084184683.
20 phase flips detected on antenna (121, 'Jnn'). A total of 71 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5116620637.
24 phase flips detected on antenna (67, 'Jnn'). A total of 30 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.511997608.
34 phase flips detected on antenna (139, 'Jnn'). A total of 81 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5193795837.
22 phase flips detected on antenna (160, 'Jnn'). A total of 44 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5116620637.
15 phase flips detected on antenna (162, 'Jnn'). A total of 160 integrations were phase-flipped relative to the 0th integration between 2460449.501372037 and 2460449.5203862167.
3 phase flips detected on antenna (55, 'Jnn'). A total of 156 integrations were phase-flipped relative to the 0th integration between 2460449.5028260625 and 2460449.5203862167.
54 phase flips detected on antenna (217, 'Jnn'). A total of 148 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5193795837.
35 phase flips detected on antenna (237, 'Jnn'). A total of 165 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5203862167.
3 phase flips detected on antenna (225, 'Jnn'). A total of 151 integrations were phase-flipped relative to the 0th integration between 2460449.5029379106 and 2460449.5203862167.
69 phase flips detected on antenna (266, 'Jnn'). A total of 164 integrations were phase-flipped relative to the 0th integration between 2460449.490634618 and 2460449.5203862167.
1 phase flips detected on antenna (193, 'Jnn'). A total of 149 integrations were phase-flipped relative to the 0th integration between 2460449.5038326955 and 2460449.5203862167.
1 phase flips detected on antenna (41, 'Jnn'). A total of 158 integrations were phase-flipped relative to the 0th integration between 2460449.5028260625 and 2460449.5203862167.
62 phase flips detected on antenna (235, 'Jnn'). A total of 97 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5193795837.
21 phase flips detected on antenna (84, 'Jnn'). A total of 153 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
21 phase flips detected on antenna (36, 'Jnn'). A total of 53 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
21 phase flips detected on antenna (159, 'Jnn'). A total of 152 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
6 phase flips detected on antenna (70, 'Jnn'). A total of 56 integrations were phase-flipped relative to the 0th integration between 2460449.5018194295 and 2460449.50830662.
49 phase flips detected on antenna (155, 'Jnn'). A total of 105 integrations were phase-flipped relative to the 0th integration between 2460449.498687682 and 2460449.5203862167.
15 phase flips detected on antenna (182, 'Jnn'). A total of 120 integrations were phase-flipped relative to the 0th integration between 2460449.50237867 and 2460449.5203862167.
22 phase flips detected on antenna (220, 'Jnn'). A total of 68 integrations were phase-flipped relative to the 0th integration between 2460449.5017075813 and 2460449.5116620637.
55 phase flips detected on antenna (281, 'Jnn'). A total of 154 integrations were phase-flipped relative to the 0th integration between 2460449.490634618 and 2460449.5203862167.
21 phase flips detected on antenna (207, 'Jnn'). A total of 126 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.5203862167.
20 phase flips detected on antenna (227, 'Jnn'). A total of 79 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.518596647.
28 phase flips detected on antenna (118, 'Jnn'). A total of 84 integrations were phase-flipped relative to the 0th integration between 2460449.5030497587 and 2460449.519938824.
2 phase flips detected on antenna (205, 'Jnn'). A total of 43 integrations were phase-flipped relative to the 0th integration between 2460449.503161607 and 2460449.5078592277.
1 phase flips detected on antenna (106, 'Jnn'). A total of 108 integrations were phase-flipped relative to the 0th integration between 2460449.5084184683 and 2460449.5203862167.
37 phase flips detected on antenna (98, 'Jnn'). A total of 135 integrations were phase-flipped relative to the 0th integration between 2460449.4942137576 and 2460449.5203862167.
34 phase flips detected on antenna (250, 'Jnn'). A total of 59 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.519938824.
5 phase flips detected on antenna (7, 'Jnn'). A total of 61 integrations were phase-flipped relative to the 0th integration between 2460449.513004241 and 2460449.5203862167.
16 phase flips detected on antenna (238, 'Jnn'). A total of 8 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.504280088.
8 phase flips detected on antenna (244, 'Jnn'). A total of 10 integrations were phase-flipped relative to the 0th integration between 2460449.501260189 and 2460449.50416824.
53 phase flips detected on antenna (135, 'Jnn'). A total of 84 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5203862167.
13 phase flips detected on antenna (8, 'Jnn'). A total of 79 integrations were phase-flipped relative to the 0th integration between 2460449.508977709 and 2460449.5203862167.
30 phase flips detected on antenna (251, 'Jnn'). A total of 167 integrations were phase-flipped relative to the 0th integration between 2460449.4952203906 and 2460449.5201625205.
5 phase flips detected on antenna (99, 'Jnn'). A total of 201 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.5203862167.
4 phase flips detected on antenna (270, 'Jnn'). A total of 18 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.4980165935.
64 phase flips detected on antenna (269, 'Jnn'). A total of 74 integrations were phase-flipped relative to the 0th integration between 2460449.4961151755 and 2460449.5193795837.
5 phase flips detected on antenna (145, 'Jnn'). A total of 110 integrations were phase-flipped relative to the 0th integration between 2460449.5028260625 and 2460449.5203862167.
40 phase flips detected on antenna (46, 'Jnn'). A total of 137 integrations were phase-flipped relative to the 0th integration between 2460449.195579288 and 2460449.47743654.
1 phase flips detected on antenna (144, 'Jnn'). A total of 2930 integrations were phase-flipped relative to the 0th integration between 2460449.192783085 and 2460449.5203862167.
1 phase flips detected on antenna (72, 'Jnn'). A total of 158 integrations were phase-flipped relative to the 0th integration between 2460449.5028260625 and 2460449.5203862167.
28 phase flips detected on antenna (100, 'Jnn'). A total of 124 integrations were phase-flipped relative to the 0th integration between 2460449.4960033274 and 2460449.518820343.

Plot results¶

In [10]:
lst_grid = utils.JD2LST(cs.time_grid) * 12 / np.pi
lst_grid[lst_grid > lst_grid[-1]] -= 24
In [11]:
def amplitude_plot(ant_to_plot):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        # Pick vmax to not saturate 90% of the abscal gains
        vmax = np.max([np.percentile(np.abs(cs.gain_grids[ant_to_plot, pol][~cs.flag_grids[ant_to_plot, pol]]), 99) for pol in ['Jee', 'Jnn']])

        display(HTML(f'<h2>Antenna {ant_to_plot} Amplitude Waterfalls</h2>'))    

        # Plot abscal gain amplitude waterfalls for a single antenna
        fig, axes = plt.subplots(4, 2, figsize=(14,14), gridspec_kw={'height_ratios': [1, 1, .4, .4]})
        for ax, pol in zip(axes[0], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)
            extent=[cs.freqs[0]/1e6, cs.freqs[-1]/1e6, lst_grid[-1], lst_grid[0]]
            im = ax.imshow(np.where(cs.flag_grids[ant], np.nan, np.abs(cs.gain_grids[ant])), aspect='auto', cmap='inferno', 
                           interpolation='nearest', vmin=0, vmax=vmax, extent=extent)
            ax.set_title(f'Smoothcal Gain Amplitude of Antenna {ant[0]}: {pol[-1]}-polarized' )
            ax.set_xlabel('Frequency (MHz)')
            ax.set_ylabel('LST (Hours)')
            ax.set_xlim([cs.freqs[0]/1e6, cs.freqs[-1]/1e6])
            ax.set_yticklabels(ax.get_yticks() % 24)
            plt.colorbar(im, ax=ax,  orientation='horizontal', pad=.15)

        # Now flagged plot abscal waterfall    
        for ax, pol in zip(axes[1], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)
            extent=[cs.freqs[0]/1e6, cs.freqs[-1]/1e6, lst_grid[-1], lst_grid[0]]
            im = ax.imshow(np.where(cs.flag_grids[ant], np.nan, np.abs(abscal_gains[ant])), aspect='auto', cmap='inferno', 
                           interpolation='nearest', vmin=0, vmax=vmax, extent=extent)
            ax.set_title(f'Abscal Gain Amplitude of Antenna {ant[0]}: {pol[-1]}-polarized' )
            ax.set_xlabel('Frequency (MHz)')
            ax.set_ylabel('LST (Hours)')
            ax.set_xlim([cs.freqs[0]/1e6, cs.freqs[-1]/1e6])
            ax.set_yticklabels(ax.get_yticks() % 24)
            plt.colorbar(im, ax=ax,  orientation='horizontal', pad=.15)
            
        # Now plot mean gain spectra 
        for ax, pol in zip(axes[2], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)   
            nflags_spectrum = np.sum(cs.flag_grids[ant], axis=0)
            to_plot = nflags_spectrum <= np.percentile(nflags_spectrum, 75)
            ax.plot(cs.freqs[to_plot] / 1e6, np.nanmean(np.where(cs.flag_grids[ant], np.nan, np.abs(abscal_gains[ant])), axis=0)[to_plot], 'r.', label='Abscal')        
            ax.plot(cs.freqs[to_plot] / 1e6, np.nanmean(np.where(cs.flag_grids[ant], np.nan, np.abs(cs.gain_grids[ant])), axis=0)[to_plot], 'k.', ms=2, label='Smoothed')        
            ax.set_ylim([0, vmax])
            ax.set_xlim([cs.freqs[0]/1e6, cs.freqs[-1]/1e6])    
            ax.set_xlabel('Frequency (MHz)')
            ax.set_ylabel('|g| (unitless)')
            ax.set_title(f'Mean Infrequently-Flagged Gain Amplitude of Antenna {ant[0]}: {pol[-1]}-polarized')
            ax.legend(loc='upper left')

        # Now plot mean gain time series
        for ax, pol in zip(axes[3], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)
            nflags_series = np.sum(cs.flag_grids[ant], axis=1)
            to_plot = nflags_series <= np.percentile(nflags_series, 75)
            ax.plot(lst_grid[to_plot], np.nanmean(np.where(cs.flag_grids[ant], np.nan, np.abs(abscal_gains[ant])), axis=1)[to_plot], 'r.', label='Abscal')        
            ax.plot(lst_grid[to_plot], np.nanmean(np.where(cs.flag_grids[ant], np.nan, np.abs(cs.gain_grids[ant])), axis=1)[to_plot], 'k.', ms=2, label='Smoothed')        
            ax.set_ylim([0, vmax])
            ax.set_xlabel('LST (hours)')
            ax.set_ylabel('|g| (unitless)')
            ax.set_title(f'Mean Infrequently-Flagged Gain Amplitude of Antenna {ant[0]}: {pol[-1]}-polarized')
            ax.set_xticklabels(ax.get_xticks() % 24)
            ax.legend(loc='upper left')

        plt.tight_layout()
        plt.show()    
In [12]:
def phase_plot(ant_to_plot):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")    
        display(HTML(f'<h2>Antenna {ant_to_plot} Phase Waterfalls</h2>'))
        fig, axes = plt.subplots(4, 2, figsize=(14,14), gridspec_kw={'height_ratios': [1, 1, .4, .4]})

        # Plot phase waterfalls for a single antenna    
        for ax, pol in zip(axes[0], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)
            extent=[cs.freqs[0]/1e6, cs.freqs[-1]/1e6, lst_grid[-1], lst_grid[0]]
            im = ax.imshow(np.where(cs.flag_grids[ant], np.nan, np.angle(cs.gain_grids[ant])), aspect='auto', cmap='inferno', 
                           interpolation='nearest', vmin=-np.pi, vmax=np.pi, extent=extent)
            ax.set_title(f'Smoothcal Gain Phase of Ant {ant[0]} / Ant {cs.refant[pol][0]}: {pol[-1]}-polarized')
            ax.set_xlabel('Frequency (MHz)')
            ax.set_ylabel('LST (Hours)')
            ax.set_xlim([cs.freqs[0]/1e6, cs.freqs[-1]/1e6])
            ax.set_yticklabels(ax.get_yticks() % 24)
            plt.colorbar(im, ax=ax,  orientation='horizontal', pad=.15)

        # Now plot abscal phase waterfall    
        for ax, pol in zip(axes[1], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)
            extent=[cs.freqs[0]/1e6, cs.freqs[-1]/1e6, lst_grid[-1], lst_grid[0]]
            im = ax.imshow(np.where(cs.flag_grids[ant], np.nan, np.angle(abscal_gains[ant])), aspect='auto', cmap='inferno', 
                           interpolation='nearest', vmin=-np.pi, vmax=np.pi, extent=extent)
            ax.set_title(f'Abscal Gain Phase of Ant {ant[0]} / Ant {cs.refant[pol][0]}: {pol[-1]}-polarized')
            ax.set_xlabel('Frequency (MHz)')
            ax.set_ylabel('LST (Hours)')
            ax.set_xlim([cs.freqs[0]/1e6, cs.freqs[-1]/1e6])
            ax.set_yticklabels(ax.get_yticks() % 24)
            plt.colorbar(im, ax=ax,  orientation='horizontal', pad=.15)
            
        # Now plot median gain spectra 
        for ax, pol in zip(axes[2], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)   
            nflags_spectrum = np.sum(cs.flag_grids[ant], axis=0)
            to_plot = nflags_spectrum <= np.percentile(nflags_spectrum, 75)
            ax.plot(cs.freqs[to_plot] / 1e6, np.nanmedian(np.where(cs.flag_grids[ant], np.nan, np.angle(abscal_gains[ant])), axis=0)[to_plot], 'r.', label='Abscal')        
            ax.plot(cs.freqs[to_plot] / 1e6, np.nanmedian(np.where(cs.flag_grids[ant], np.nan, np.angle(cs.gain_grids[ant])), axis=0)[to_plot], 'k.', ms=2, label='Smoothed')        
            ax.set_ylim([-np.pi, np.pi])
            ax.set_xlim([cs.freqs[0]/1e6, cs.freqs[-1]/1e6])    
            ax.set_xlabel('Frequency (MHz)')
            ax.set_ylabel(f'Phase of g$_{{{ant[0]}}}$ / g$_{{{cs.refant[pol][0]}}}$')
            ax.set_title(f'Median Infrequently-Flagged Gain Phase of Ant {ant[0]} / Ant {cs.refant[pol][0]}: {pol[-1]}-polarized')
            ax.legend(loc='upper left')

        # # Now plot median gain time series
        for ax, pol in zip(axes[3], ['Jee', 'Jnn']):
            ant = (ant_to_plot, pol)
            nflags_series = np.sum(cs.flag_grids[ant], axis=1)
            to_plot = nflags_series <= np.percentile(nflags_series, 75)
            ax.plot(lst_grid[to_plot], np.nanmean(np.where(cs.flag_grids[ant], np.nan, np.angle(abscal_gains[ant])), axis=1)[to_plot], 'r.', label='Abscal')        
            ax.plot(lst_grid[to_plot], np.nanmean(np.where(cs.flag_grids[ant], np.nan, np.angle(cs.gain_grids[ant])), axis=1)[to_plot], 'k.', ms=2, label='Smoothed')        
            ax.set_ylim([-np.pi, np.pi])    
            ax.set_xlabel('LST (hours)')
            ax.set_ylabel(f'Phase of g$_{{{ant[0]}}}$ / g$_{{{cs.refant[pol][0]}}}$')
            ax.set_title(f'Mean Infrequently-Flagged Gain Phase of Ant {ant[0]} / Ant {cs.refant[pol][0]}: {pol[-1]}-polarized')
            ax.set_xticklabels(ax.get_xticks() % 24)    
            ax.legend(loc='upper left')

        plt.tight_layout()
        plt.show()

Figure 1: Full-Day Gain Amplitudes Before and After smooth_cal¶

Here we plot abscal and smooth_cal gain amplitudes for both of the sample antennas. We also show means across time/frequency, excluding frequencies/times that are frequently flagged.

In [13]:
for ant_to_plot in ants_to_plot:
    amplitude_plot(ant_to_plot)

Antenna 4 Amplitude Waterfalls

No description has been provided for this image

Antenna 295 Amplitude Waterfalls

No description has been provided for this image

Figure 2: Full-Day Gain Phases Before and After smooth_cal¶

Here we plot abscal and smooth_cal phases relative to each polarization's reference antenna for both of the sample antennas. We also show medians across time/frequency, excluding frequencies/times that are frequently flagged.

In [14]:
for ant_to_plot in ants_to_plot:
    phase_plot(ant_to_plot)

Antenna 4 Phase Waterfalls

No description has been provided for this image

Antenna 295 Phase Waterfalls

No description has been provided for this image

Examine $\chi^2$¶

In [15]:
def chisq_plot():
    fig, axes = plt.subplots(1, 2, figsize=(14, 10), sharex=True, sharey=True)
    extent = [cs.freqs[0]/1e6, cs.freqs[-1]/1e6, lst_grid[-1], lst_grid[0]]
    for ax, pol in zip(axes, ['Jee', 'Jnn']):

        im = ax.imshow(np.where(cs.flag_grids[cs.refant[pol]], np.nan, cs.chisq_grids[pol]), vmin=1, vmax=5, 
                       aspect='auto', cmap='turbo', interpolation='none', extent=extent)
        ax.set_yticklabels(ax.get_yticks() % 24)
        ax.set_title(f'{pol[1:]}-Polarized $\\chi^2$ / DoF')
        ax.set_xlabel('Frequency (MHz)')

    axes[0].set_ylabel('LST (hours)')
    plt.tight_layout()
    fig.colorbar(im, ax=axes, pad=.07, label='$\\chi^2$ / DoF', orientation='horizontal', extend='both', aspect=50)

Figure 3: Full-Day $\chi^2$ / DoF Waterfall from Redundant-Baseline Calibration¶

Here we plot $\chi^2$ per degree of freedom from redundant-baseline calibration for both polarizations separately. While this plot is a little out of place, as it was not produced by this notebook, it is a convenient place where all the necessary components are readily available. If the array were perfectly redundant and any non-redundancies in the calibrated visibilities were explicable by thermal noise alone, this waterfall should be all 1.

In [16]:
chisq_plot()
set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
No description has been provided for this image
In [17]:
avg_cspa_vs_time = {ant: np.nanmean(np.where(cs.flag_grids[ant], np.nan, cs.cspa_grids[ant]), axis=1) for ant in cs.ants}
avg_cspa_vs_freq = {ant: np.nanmean(np.where(cs.flag_grids[ant], np.nan, cs.cspa_grids[ant]), axis=0) for ant in cs.ants}
Mean of empty slice
Mean of empty slice
In [18]:
def cspa_vs_time_plot():
    fig, axes = plt.subplots(2, 1, figsize=(14, 8), sharex=True, sharey=True, gridspec_kw={'hspace': 0})
    for ax, pol in zip(axes, ['Jee', 'Jnn']):
        detail_cutoff = np.percentile([np.nanmean(m) for ant, m in avg_cspa_vs_time.items() 
                                       if ant[1] == pol and np.isfinite(np.nanmean(m))], 95)
        for ant in avg_cspa_vs_time:
            if ant[1] == pol and not np.all(cs.flag_grids[ant]):
                if np.nanmean(avg_cspa_vs_time[ant]) > detail_cutoff:
                    ax.plot(lst_grid, avg_cspa_vs_time[ant], label=ant, zorder=100)
                else:
                    ax.plot(lst_grid, avg_cspa_vs_time[ant], c='grey', alpha=.2, lw=.5)
        ax.legend(title=f'{pol[1:]}-Polarized', ncol=2)
        ax.set_ylabel('Mean Unflagged $\\chi^2$ per Antenna')
        ax.set_xlabel('LST (hours)')
        ax.set_xticklabels(ax.get_xticks() % 24)

    plt.ylim([1, 5.4])
    plt.tight_layout()
In [19]:
def cspa_vs_freq_plot():
    fig, axes = plt.subplots(2, 1, figsize=(14, 6), sharex=True, sharey=True, gridspec_kw={'hspace': 0})
    for ax, pol in zip(axes, ['Jee', 'Jnn']):
        detail_cutoff = np.percentile([np.nanmean(m) for ant, m in avg_cspa_vs_freq.items() 
                                       if ant[1] == pol and np.isfinite(np.nanmean(m))], 95)
        for ant in avg_cspa_vs_freq:
            if ant[1] == pol and not np.all(cs.flag_grids[ant]):
                if np.nanmean(avg_cspa_vs_freq[ant]) > detail_cutoff:
                    ax.plot(cs.freqs / 1e6, avg_cspa_vs_freq[ant], label=ant, zorder=100)
                else:
                    ax.plot(cs.freqs / 1e6, avg_cspa_vs_freq[ant], c='grey', alpha=.2, lw=.5)
        ax.legend(title=f'{pol[1:]}-Polarized', ncol=2)
        ax.set_ylabel('Mean Unflagged $\\chi^2$ per Antenna')
        ax.set_xlabel('Frequency (MHz)')

    plt.ylim([1, 5.4])
    plt.tight_layout()

Figure 4: Average $\chi^2$ per Antenna vs. Time and Frequency¶

Here we plot $\chi^2$ per antenna from redundant-baseline calibration, separating polarizations and averaging the unflagged pixels in the waterfalls over frequency or time. The worst 5% of antennas are shown in color and highlighted in the legends, the rest are shown in grey.

In [20]:
cspa_vs_time_plot()
cspa_vs_freq_plot()
Mean of empty slice
set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
Mean of empty slice
set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
Mean of empty slice
No description has been provided for this image
No description has been provided for this image

Save Results¶

In [21]:
add_to_history = 'Produced by calibration_smoothing notebook with the following environment:\n' + '=' * 65 + '\n' + os.popen('conda env export').read() + '=' * 65
In [22]:
cs.write_smoothed_cal(output_replace=(CAL_SUFFIX, SMOOTH_CAL_SUFFIX), add_to_history=add_to_history, clobber=True)
Mean of empty slice

Metadata¶

In [23]:
for repo in ['hera_cal', 'hera_qm', 'hera_filters', 'hera_notebook_templates', 'pyuvdata']:
    exec(f'from {repo} import __version__')
    print(f'{repo}: {__version__}')
hera_cal: 3.6.dev65+ge56a686
hera_qm: 2.1.4
hera_filters: 0.1.5
hera_notebook_templates: 0.1.dev734+g90f16f4
pyuvdata: 2.4.2
In [24]:
print(f'Finished execution in {(time.time() - tstart) / 60:.2f} minutes.')
Finished execution in 18.65 minutes.