Calibration Smoothing¶
by Josh Dillon, last updated July 19, 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
A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.0.2 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.
If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.
Traceback (most recent call last): File "<frozen runpy>", line 198, in _run_module_as_main
File "<frozen runpy>", line 88, in _run_code
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel_launcher.py", line 18, in <module>
app.launch_new_instance()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/traitlets/config/application.py", line 1075, in launch_instance
app.start()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelapp.py", line 739, in start
self.io_loop.start()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/tornado/platform/asyncio.py", line 205, in start
self.asyncio_loop.run_forever()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/base_events.py", line 640, in run_forever
self._run_once()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/base_events.py", line 1992, in _run_once
handle._run()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/events.py", line 88, in _run
self._context.run(self._callback, *self._args)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 545, in dispatch_queue
await self.process_one()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 534, in process_one
await dispatch(*args)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 437, in dispatch_shell
await result
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/ipkernel.py", line 362, in execute_request
await super().execute_request(stream, ident, parent)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 778, in execute_request
reply_content = await reply_content
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/ipkernel.py", line 449, in do_execute
res = shell.run_cell(
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/zmqshell.py", line 549, in run_cell
return super().run_cell(*args, **kwargs)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3075, in run_cell
result = self._run_cell(
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3130, in _run_cell
result = runner(coro)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/async_helpers.py", line 128, in _pseudo_sync_runner
coro.send(None)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3334, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3517, in run_ast_nodes
if await self.run_code(code, result, async_=asy):
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3577, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "/tmp/ipykernel_3228735/530180873.py", line 11, in <module>
from hera_cal import io, utils, smooth_cal
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/__init__.py", line 22, in <module>
from . import utils
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/utils.py", line 35, in <module>
import aipy
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/__init__.py", line 15, in <module>
from . import phs, const, coord, deconv
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/phs.py", line 10, in <module>
from .miriad import ij2bl, bl2ij
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/miriad.py", line 12, in <module>
from . import _miriad
--------------------------------------------------------------------------- ImportError Traceback (most recent call last) File /lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/numpy/core/_multiarray_umath.py:44, in __getattr__(attr_name) 39 # Also print the message (with traceback). This is because old versions 40 # of NumPy unfortunately set up the import to replace (and hide) the 41 # error. The traceback shouldn't be needed, but e.g. pytest plugins 42 # seem to swallow it and we should be failing anyway... 43 sys.stderr.write(msg + tb_msg) ---> 44 raise ImportError(msg) 46 ret = getattr(_multiarray_umath, attr_name, None) 47 if ret is None: ImportError: A module that was compiled using NumPy 1.x cannot be run in NumPy 2.0.2 as it may crash. To support both 1.x and 2.x versions of NumPy, modules must be compiled with NumPy 2.0. Some module may need to rebuild instead e.g. with 'pybind11>=2.12'. If you are a user of the module, the easiest solution will be to downgrade to 'numpy<2' or try to upgrade the affected module. We expect that some modules will need time to support NumPy 2.
A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.0.2 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.
If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.
Traceback (most recent call last): File "<frozen runpy>", line 198, in _run_module_as_main
File "<frozen runpy>", line 88, in _run_code
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel_launcher.py", line 18, in <module>
app.launch_new_instance()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/traitlets/config/application.py", line 1075, in launch_instance
app.start()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelapp.py", line 739, in start
self.io_loop.start()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/tornado/platform/asyncio.py", line 205, in start
self.asyncio_loop.run_forever()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/base_events.py", line 640, in run_forever
self._run_once()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/base_events.py", line 1992, in _run_once
handle._run()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/events.py", line 88, in _run
self._context.run(self._callback, *self._args)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 545, in dispatch_queue
await self.process_one()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 534, in process_one
await dispatch(*args)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 437, in dispatch_shell
await result
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/ipkernel.py", line 362, in execute_request
await super().execute_request(stream, ident, parent)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 778, in execute_request
reply_content = await reply_content
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/ipkernel.py", line 449, in do_execute
res = shell.run_cell(
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/zmqshell.py", line 549, in run_cell
return super().run_cell(*args, **kwargs)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3075, in run_cell
result = self._run_cell(
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3130, in _run_cell
result = runner(coro)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/async_helpers.py", line 128, in _pseudo_sync_runner
coro.send(None)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3334, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3517, in run_ast_nodes
if await self.run_code(code, result, async_=asy):
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3577, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "/tmp/ipykernel_3228735/530180873.py", line 11, in <module>
from hera_cal import io, utils, smooth_cal
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/__init__.py", line 23, in <module>
from . import redcal
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/redcal.py", line 13, in <module>
from .noise import predict_noise_variance_from_autos, infer_dt
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/noise.py", line 12, in <module>
from . import io
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/io.py", line 40, in <module>
import aipy
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/__init__.py", line 15, in <module>
from . import phs, const, coord, deconv
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/phs.py", line 10, in <module>
from .miriad import ij2bl, bl2ij
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/miriad.py", line 12, in <module>
from . import _miriad
--------------------------------------------------------------------------- ImportError Traceback (most recent call last) File /lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/numpy/core/_multiarray_umath.py:44, in __getattr__(attr_name) 39 # Also print the message (with traceback). This is because old versions 40 # of NumPy unfortunately set up the import to replace (and hide) the 41 # error. The traceback shouldn't be needed, but e.g. pytest plugins 42 # seem to swallow it and we should be failing anyway... 43 sys.stderr.write(msg + tb_msg) ---> 44 raise ImportError(msg) 46 ret = getattr(_multiarray_umath, attr_name, None) 47 if ret is None: ImportError: A module that was compiled using NumPy 1.x cannot be run in NumPy 2.0.2 as it may crash. To support both 1.x and 2.x versions of NumPy, modules must be compiled with NumPy 2.0. Some module may need to rebuild instead e.g. with 'pybind11>=2.12'. If you are a user of the module, the easiest solution will be to downgrade to 'numpy<2' or try to upgrade the affected module. We expect that some modules will need time to support NumPy 2.
A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.0.2 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.
If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.
Traceback (most recent call last): File "<frozen runpy>", line 198, in _run_module_as_main
File "<frozen runpy>", line 88, in _run_code
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel_launcher.py", line 18, in <module>
app.launch_new_instance()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/traitlets/config/application.py", line 1075, in launch_instance
app.start()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelapp.py", line 739, in start
self.io_loop.start()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/tornado/platform/asyncio.py", line 205, in start
self.asyncio_loop.run_forever()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/base_events.py", line 640, in run_forever
self._run_once()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/base_events.py", line 1992, in _run_once
handle._run()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/asyncio/events.py", line 88, in _run
self._context.run(self._callback, *self._args)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 545, in dispatch_queue
await self.process_one()
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 534, in process_one
await dispatch(*args)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 437, in dispatch_shell
await result
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/ipkernel.py", line 362, in execute_request
await super().execute_request(stream, ident, parent)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/kernelbase.py", line 778, in execute_request
reply_content = await reply_content
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/ipkernel.py", line 449, in do_execute
res = shell.run_cell(
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/ipykernel/zmqshell.py", line 549, in run_cell
return super().run_cell(*args, **kwargs)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3075, in run_cell
result = self._run_cell(
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3130, in _run_cell
result = runner(coro)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/async_helpers.py", line 128, in _pseudo_sync_runner
coro.send(None)
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3334, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3517, in run_ast_nodes
if await self.run_code(code, result, async_=asy):
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3577, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "/tmp/ipykernel_3228735/530180873.py", line 11, in <module>
from hera_cal import io, utils, smooth_cal
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/__init__.py", line 26, in <module>
from . import abscal
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/abscal.py", line 43, in <module>
from .smooth_cal import pick_reference_antenna, rephase_to_refant
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/hera_cal/smooth_cal.py", line 15, in <module>
import aipy
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/__init__.py", line 15, in <module>
from . import phs, const, coord, deconv
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/phs.py", line 10, in <module>
from .miriad import ij2bl, bl2ij
File "/lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/aipy/miriad.py", line 12, in <module>
from . import _miriad
--------------------------------------------------------------------------- ImportError Traceback (most recent call last) File /lustre/aoc/projects/hera/heramgr/anaconda3/envs/h6c_idr2_validation/lib/python3.12/site-packages/numpy/core/_multiarray_umath.py:44, in __getattr__(attr_name) 39 # Also print the message (with traceback). This is because old versions 40 # of NumPy unfortunately set up the import to replace (and hide) the 41 # error. The traceback shouldn't be needed, but e.g. pytest plugins 42 # seem to swallow it and we should be failing anyway... 43 sys.stderr.write(msg + tb_msg) ---> 44 raise ImportError(msg) 46 ret = getattr(_multiarray_umath, attr_name, None) 47 if ret is None: ImportError: A module that was compiled using NumPy 1.x cannot be run in NumPy 2.0.2 as it may crash. To support both 1.x and 2.x versions of NumPy, modules must be compiled with NumPy 2.0. Some module may need to rebuild instead e.g. with 'pybind11>=2.12'. If you are a user of the module, the easiest solution will be to downgrade to 'numpy<2' or try to upgrade the affected module. We expect that some modules will need time to support NumPy 2.
Parse inputs¶
In [3]:
# get files
SUM_FILE = os.environ.get("SUM_FILE", None)
# SUM_FILE = '/lustre/aoc/projects/hera/h6c-analysis/IDR2/2459867/zen.2459867.43004.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))
PER_POL_REFANT = os.environ.get("PER_POL_REFANT", "False").upper() == "TRUE"
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', 'PER_POL_REFANT']:
print(f'{setting} = {eval(setting)}')
SUM_FILE = /lustre/aoc/projects/hera/Validation/H6C_IDR2/2459861/zen.2459861.25297.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 PER_POL_REFANT = False
Load files and select reference antenna(s)¶
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 1862 *.sum.omni.calfits files starting with /lustre/aoc/projects/hera/Validation/H6C_IDR2/2459861/zen.2459861.25297.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 1862 *.sum.flag_waterfall.h5 files starting with /lustre/aoc/projects/hera/Validation/H6C_IDR2/2459861/zen.2459861.25297.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 1862 *.sum.antenna_flags.h5 files starting with /lustre/aoc/projects/hera/Validation/H6C_IDR2/2459861/zen.2459861.25297.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=False, load_chisq=True, load_cspa=True)
In [8]:
cs.refant = smooth_cal.pick_reference_antenna(cs.gain_grids, cs.flag_grids, cs.freqs, per_pol=True)
for pol in cs.refant:
print(f'Reference antenna {cs.refant[pol][0]} selected for smoothing {pol} gains.')
if not PER_POL_REFANT:
# in this case, rephase both pols separately before smoothing, but also smooth the relative polarization calibration phasor
overall_refant = smooth_cal.pick_reference_antenna({ant: cs.gain_grids[ant] for ant in cs.refant.values()},
{ant: cs.flag_grids[ant] for ant in cs.refant.values()},
cs.freqs, per_pol=False)
print(f'Overall reference antenna {overall_refant} selected.')
other_refant = [ant for ant in cs.refant.values() if ant != overall_refant][0]
relative_pol_phasor = cs.gain_grids[overall_refant] * cs.gain_grids[other_refant].conj() # TODO: is this conjugation right?
relative_pol_phasor /= np.abs(relative_pol_phasor)
Reference antenna 87 selected for smoothing Jnn gains. Reference antenna 123 selected for smoothing Jee gains.
Overall reference antenna (np.int64(87), 'Jnn') selected.
In [9]:
# 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 = {}
for pol in ['Jee', 'Jnn']:
for antnum in ants_to_plot:
refant_here = (cs.refant[pol] if PER_POL_REFANT else overall_refant)
abscal_gains[antnum, pol] = cs.gain_grids[(antnum, pol)] * np.abs(cs.gain_grids[refant_here]) / cs.gain_grids[refant_here]
In [10]:
cs.rephase_to_refant(propagate_refant_flags=True)
Perform smoothing¶
In [11]:
if not PER_POL_REFANT:
# treat the relative_pol_phasor as if it were antenna -1
cs.gain_grids[(-1, other_refant[1])] = relative_pol_phasor
cs.flag_grids[(-1, other_refant[1])] = cs.flag_grids[overall_refant] | cs.flag_grids[other_refant]
In [12]:
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)
In [13]:
if not PER_POL_REFANT:
# put back in the smoothed phasor, ensuring the amplitude is 1 and that data are flagged anywhere either polarization's refant is flagged
smoothed_relative_pol_phasor = cs.gain_grids[(-1, other_refant[-1])] / np.abs(cs.gain_grids[(-1, other_refant[-1])])
for ant in cs.gain_grids:
if ant[0] >= 0 and ant[1] == other_refant[1]:
cs.gain_grids[ant] /= smoothed_relative_pol_phasor
cs.flag_grids[ant] |= (cs.flag_grids[(-1, other_refant[1])])
cs.refant = overall_refant
Plot results¶
In [14]:
lst_grid = utils.JD2LST(cs.time_grid) * 12 / np.pi
lst_grid[lst_grid > lst_grid[-1]] -= 24
In [15]:
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 [16]:
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)
refant = (cs.refant[pol] if isinstance(cs.refant, dict) else cs.refant)
ax.set_title(f'Smoothcal Gain Phase of Ant {ant[0]}{pol[-1]} / Ant {refant[0]}{refant[1][-1]}')
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)
refant = (cs.refant[pol] if isinstance(cs.refant, dict) else cs.refant)
ax.set_title(f'Abscal Gain Phase of Ant {ant[0]}{pol[-1]} / Ant {refant[0]}{refant[1][-1]}')
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)')
refant = (cs.refant[pol] if isinstance(cs.refant, dict) else cs.refant)
ax.set_ylabel(f'Phase of g$_{{{ant[0]}{pol[-1]}}}$ / g$_{{{refant[0]}{refant[1][-1]}}}$')
ax.set_title(f'Median Infrequently-Flagged Gain Phase of Ant {ant[0]}{pol[-1]} / Ant {refant[0]}{refant[1][-1]}')
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)')
refant = (cs.refant[pol] if isinstance(cs.refant, dict) else cs.refant)
ax.set_ylabel(f'Phase of g$_{{{ant[0]}{pol[-1]}}}$ / g$_{{{refant[0]}{refant[1][-1]}}}$')
ax.set_title(f'Mean Infrequently-Flagged Gain Phase of Ant {ant[0]}{pol[-1]} / Ant {refant[0]}{refant[1][-1]}')
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 [17]:
for ant_to_plot in ants_to_plot:
amplitude_plot(ant_to_plot)
Antenna 3 Amplitude Waterfalls
Antenna 239 Amplitude Waterfalls
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 [18]:
for ant_to_plot in ants_to_plot:
phase_plot(ant_to_plot)
Antenna 3 Phase Waterfalls
Antenna 239 Phase Waterfalls
Examine $\chi^2$¶
In [19]:
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']):
refant = (cs.refant[pol] if isinstance(cs.refant, dict) else cs.refant)
im = ax.imshow(np.where(cs.flag_grids[refant], 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 [20]:
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.
In [21]:
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 [22]:
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 [23]:
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 [24]:
cspa_vs_time_plot()
cspa_vs_freq_plot()
Mean of empty slice Passing label as a length 2 sequence when plotting a single dataset is deprecated in Matplotlib 3.9 and will error in 3.11. To keep the current behavior, cast the sequence to string before passing. 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 Passing label as a length 2 sequence when plotting a single dataset is deprecated in Matplotlib 3.9 and will error in 3.11. To keep the current behavior, cast the sequence to string before passing.
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 Passing label as a length 2 sequence when plotting a single dataset is deprecated in Matplotlib 3.9 and will error in 3.11. To keep the current behavior, cast the sequence to string before passing.
Save Results¶
In [25]:
add_to_history = 'Produced by calibration_smoothing notebook with the following environment:\n' + '=' * 65 + '\n' + os.popen('conda env export').read() + '=' * 65
In [26]:
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 [27]:
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.2.dev138+ga1987648 hera_qm: 2.2.0 hera_filters: 0.1.6.dev1+g297dcce
hera_notebook_templates: 0.1.dev996+g3fb11cf pyuvdata: 3.1.3.dev32+ged715863
In [28]:
print(f'Finished execution in {(time.time() - tstart) / 60:.2f} minutes.')
Finished execution in 116.67 minutes.