Guide for hera_sim Defaults and Simulator

This notebook is intended to be a guide for interfacing with the hera_sim.defaults module and using the hera_sim.Simulator class with the run_sim class method.

import matplotlib.pyplot as plt
import numpy as np
import uvtools

import hera_sim
from hera_sim import DATA_PATH, Simulator
from hera_sim.config import CONFIG_PATH

%matplotlib inline

We’ll be using the uvtools.plot.waterfall function a few times throughout the notebook, but with a standardized way of generating plots. So let’s write a wrapper to make these plots nicer:

# let's wrap it so that it'll work on Simulator objects
def waterfall(sim, antpairpol):
    """
    For reference, sim is a Simulator object, and antpairpol is a tuple with the
    form (ant1, ant2, pol).
    """
    freqs = np.unique(sim.data.freq_array) * 1e-9  # GHz
    lsts = np.unique(sim.data.lst_array)  # radians
    vis = sim.data.get_data(antpairpol)

    # extent format is [left, right, bottom, top], vis shape is (NTIMES,NFREQS)
    extent = [freqs.min(), freqs.max(), lsts.max(), lsts.min()]

    fig = plt.figure(figsize=(12, 8))
    axes = fig.subplots(2, 1, sharex=True)
    axes[1].set_xlabel("Frequency [GHz]", fontsize=12)
    for ax in axes:
        ax.set_ylabel("LST [rad]", fontsize=12)

    fig.sca(axes[0])
    cax = uvtools.plot.waterfall(vis, mode="log", extent=extent)
    fig.colorbar(cax, label=r"$\log_{10}(V$/Jy)")

    fig.sca(axes[1])
    cax = uvtools.plot.waterfall(vis, mode="phs", extent=extent)
    fig.colorbar(cax, label="Phase [rad]")

    plt.tight_layout()
    plt.show()

# while we're preparing things, let's make a dictionary of default settings for
# creating a Simulator object

# choose a modest number of frequencies and times to simulate
NFREQ = 128
NTIMES = 32

# choose the channel width so that the bandwidth is 100 MHz
channel_width = 1e8 / NFREQ

# use just two antennas
ants = {0: (20.0, 20.0, 0), 1: (50.0, 50.0, 0)}

# use cross- and auto-correlation baselines
no_autos = False

# actually make the dictionary of initialization parameters
init_params = {
    "n_freq": NFREQ,
    "n_times": NTIMES,
    "antennas": ants,
    "channel_width": channel_width,
    "no_autos": no_autos,
}

# turn off warnings; remove this cell in the future when this isn't a feature
hera_sim.defaults._warn = False

Defaults Basics

Let’s first go over the basics of the defaults module. There are three methods that serve as the primary interface between the user and the defaults module: set, activate, and deactivate. These do what you’d expect them to do, but the set method requires a bit of extra work to interface with if you want to set custom default parameter values. We’ll cover this later; for now, we’ll see how you can switch between defaults characteristic to different observing seasons. Currently, we support default settings for the H1C observing season and the H2C observing season, and these may be accessed using the strings 'h1c' and 'h2c', respectively.

# first, let's note that the defaults are initially deactivated
hera_sim.defaults._override_defaults

False

# so let's activate the season defaults
hera_sim.defaults.activate()
hera_sim.defaults._override_defaults

True

# now that the defaults are activated, let's see what some differences are
fqs = np.linspace(0.1, 0.2, 1024)
lsts = np.linspace(0, 2 * np.pi, 100)

# note that the defaults object is initialized with H1C defaults, but let's
# make sure that it's set to H1C defaults in case this cell is rerun later
hera_sim.defaults.set("h1c")

h1c_beam = hera_sim.defaults("omega_p")(fqs)
h1c_bandpass = hera_sim.defaults("bp_poly")(fqs)

# now change the defaults to the H2C observing season
hera_sim.defaults.set("h2c")

h2c_beam = hera_sim.defaults("omega_p")(fqs)
h2c_bandpass = hera_sim.defaults("bp_poly")(fqs)

# now compare them
print(
    f"""
H1C Defaults:
  Beam Polynomial:
    {h1c_beam}
  Bandpass Polynomial:
    {h1c_bandpass}
"""
)
print(
    f"""
H2C Defaults:
  Beam Polynomial:
    {h2c_beam}
  Bandpass Polynomial:
    {h2c_bandpass}
"""
)

H1C Defaults:
  Beam Polynomial:
    [0.08033913 0.08020828 0.08007931 ... 0.02888961 0.0288856  0.02888128]
  Bandpass Polynomial:
    [-0.00524501 -0.00463106 -0.00402282 ... -0.00538521 -0.00568868
 -0.00599484]


H2C Defaults:
  Beam Polynomial:
    [0.10805394 0.10794104 0.10783932 ... 0.03401091 0.03393028 0.03384834]
  Bandpass Polynomial:
    [10.65096904 10.65504987 10.65915031 ... 22.07424626 22.06666646
 22.05905975]

# another thing
noise = hera_sim.noise.thermal_noise

hera_sim.defaults.set("h1c")
rng = np.random.default_rng(0)
h1c_noise = noise(lsts, fqs, rng=rng)

hera_sim.defaults.set("h2c")
rng = np.random.default_rng(0)
h2c_noise = noise(lsts, fqs, rng=rng)
rng = np.random.default_rng(0)
still_h2c_noise = noise(lsts, fqs, rng=rng)

# passing in a kwarg
other_noise = noise(lsts, fqs, omega_p=np.ones(fqs.size), rng=rng)

# check things
print(f"H1C noise identical to H2C noise? {np.all(h1c_noise == h2c_noise)}")
print(
    f"H2C noise identical to its other version? {np.all(h2c_noise == still_h2c_noise)}"
)

_identical = np.all(h1c_noise == other_noise) or np.all(h2c_noise == other_noise)
print(f"Either noise identical to other noise? {_identical}")

H1C noise identical to H2C noise? False
H2C noise identical to its other version? True
Either noise identical to other noise? False

# check out what the beam and bandpass look like
seasons = ("h1c", "h2c")
beams = {"h1c": h1c_beam, "h2c": h2c_beam}
bps = {"h1c": h1c_bandpass, "h2c": h2c_bandpass}
fig = plt.figure(figsize=(12, 10))
axes = fig.subplots(2, 2)
for j, ax in enumerate(axes[:, 0]):
    seas = seasons[j]
    ax.set_xlabel("Frequency [MHz]", fontsize=12)
    ax.set_ylabel("Beam Size [sr]", fontsize=12)
    ax.set_title(f"HERA {seas.upper()} Beam Model", fontsize=12)
    ax.plot(fqs * 1e3, beams[seas])
for j, ax in enumerate(axes[:, 1]):
    seas = seasons[j]
    ax.set_xlabel("Frequency [MHz]", fontsize=12)
    ax.set_ylabel("Noiseless Bandpass Gain", fontsize=12)
    ax.set_title(f"HERA {seas.upper()} Bandpass Model", fontsize=12)
    ax.plot(fqs * 1e3, bps[seas])
plt.tight_layout()
plt.show()

# let's look at the configuration that's initially loaded in
hera_sim.defaults._raw_config

{'setup': {'frequency_array': {'Nfreqs': 1638,
   'channel_width': 122070.3125,
   'start_freq': 46920776.3671875},
  'time_array': {'Ntimes': 100,
   'integration_time': 8.59,
   'start_time': 2458119.5}},
 'telescope': {'array_layout': {0: array([-14.6       ,  25.28794179,   0.        ]),
   1: array([ 0.        , 25.28794179,  0.        ]),
   2: array([14.6       , 25.28794179,  0.        ]),
   3: array([-21.9      ,  12.6439709,   0.       ]),
   4: array([-7.3      , 12.6439709,  0.       ]),
   5: array([ 7.3      , 12.6439709,  0.       ]),
   6: array([21.9      , 12.6439709,  0.       ]),
   7: array([-29.2,   0. ,   0. ]),
   8: array([-14.6,   0. ,   0. ]),
   9: array([0., 0., 0.]),
   10: array([14.6,  0. ,  0. ]),
   11: array([29.2,  0. ,  0. ]),
   12: array([-21.9      , -12.6439709,   0.       ]),
   13: array([ -7.3      , -12.6439709,   0.       ]),
   14: array([  7.3      , -12.6439709,   0.       ]),
   15: array([ 21.9      , -12.6439709,   0.       ]),
   16: array([-14.6       , -25.28794179,   0.        ]),
   17: array([  0.        , -25.28794179,   0.        ]),
   18: array([ 14.6       , -25.28794179,   0.        ])},
  'bp_poly': <hera_sim.interpolators.Bandpass at 0x7fbbd5548ef0>,
  'omega_p': <hera_sim.interpolators.Beam at 0x7fbbd5549400>,
  'delay_filter_type': 'tophat',
  'fringe_filter_type': 'tophat'},
 'sky': {'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7fbbd550ea50>}}

# and what about the set of defaults actually used?
hera_sim.defaults._config

{'Nfreqs': 1638,
 'channel_width': 122070.3125,
 'start_freq': 46920776.3671875,
 'Ntimes': 100,
 'integration_time': 8.59,
 'start_time': 2458119.5,
 'array_layout': {0: array([-14.6       ,  25.28794179,   0.        ]),
  1: array([ 0.        , 25.28794179,  0.        ]),
  2: array([14.6       , 25.28794179,  0.        ]),
  3: array([-21.9      ,  12.6439709,   0.       ]),
  4: array([-7.3      , 12.6439709,  0.       ]),
  5: array([ 7.3      , 12.6439709,  0.       ]),
  6: array([21.9      , 12.6439709,  0.       ]),
  7: array([-29.2,   0. ,   0. ]),
  8: array([-14.6,   0. ,   0. ]),
  9: array([0., 0., 0.]),
  10: array([14.6,  0. ,  0. ]),
  11: array([29.2,  0. ,  0. ]),
  12: array([-21.9      , -12.6439709,   0.       ]),
  13: array([ -7.3      , -12.6439709,   0.       ]),
  14: array([  7.3      , -12.6439709,   0.       ]),
  15: array([ 21.9      , -12.6439709,   0.       ]),
  16: array([-14.6       , -25.28794179,   0.        ]),
  17: array([  0.        , -25.28794179,   0.        ]),
  18: array([ 14.6       , -25.28794179,   0.        ])},
 'bp_poly': <hera_sim.interpolators.Bandpass at 0x7fbbd5548ef0>,
 'omega_p': <hera_sim.interpolators.Beam at 0x7fbbd5549400>,
 'delay_filter_type': 'tophat',
 'fringe_filter_type': 'tophat',
 'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7fbbd550ea50>}

# let's make two simulator objects
sim1 = Simulator(**init_params)
sim2 = Simulator(**init_params)

The default baseline conjugation convention has changed. In the past it was 'ant2<ant1', it now defaults to 'ant1<ant2'. You can specify the baseline conjugation convention in `obs_param` by setting the obs_param['ordering']['conjugation_convention'] field. This warning will go away in version 1.5.
The use_grid_alg parameter is not set. Defaulting to True to use the new gridding based algorithm (developed by the HERA team) rather than the older clustering based algorithm. This is change to the default, to use the clustering algorithm set use_grid_alg=False.
The default baseline conjugation convention has changed. In the past it was 'ant2<ant1', it now defaults to 'ant1<ant2'. You can specify the baseline conjugation convention in `obs_param` by setting the obs_param['ordering']['conjugation_convention'] field. This warning will go away in version 1.5.
The use_grid_alg parameter is not set. Defaulting to True to use the new gridding based algorithm (developed by the HERA team) rather than the older clustering based algorithm. This is change to the default, to use the clustering algorithm set use_grid_alg=False.

# parameters for two different simulations
hera_sim.defaults.set("h1c")
sim1params = {"pntsrc_foreground": {}, "noiselike_eor": {}, "diffuse_foreground": {}}
sim1.run_sim(**sim1params)

hera_sim.defaults.set("h2c")
sim2params = {"pntsrc_foreground": {}, "noiselike_eor": {}, "diffuse_foreground": {}}
sim2.run_sim(**sim2params)

antpairpol = (0, 1, "xx")
waterfall(sim1, antpairpol)

waterfall(sim2, antpairpol)