hera_sim: A Simple Simulator

hera_sim is a simple simulator that generates instrumental effects and applies them to visibilities.

Contents

Installation

Requirements

Requires:

• numpy

• scipy

• aipy

• pyuvdata

Then, at the command line, navigate to the hera_sim repo/directory, and:

pip install .

If developing, from the top-level directory do:

pip install -e .

Tutorials and FAQs

The following introductory tutorial will help you get started with hera_sim:

Tour of hera_sim

This notebook briefly introduces some of the effects that can be modeled with hera_sim.

%matplotlib notebook
import aipy, uvtools
import numpy as np
import pylab as plt

from hera_sim import foregrounds, noise, sigchain, rfi

/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/visibilities/__init__.py:27: UserWarning: PRISim failed to import.
  warnings.warn("PRISim failed to import.")
/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/visibilities/__init__.py:33: UserWarning: VisGPU failed to import.
  warnings.warn("VisGPU failed to import.")
/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/__init__.py:36: FutureWarning:
In the next major release, all HERA-specific variables will be removed from the codebase. The following variables will need to be accessed through new class-like structures to be introduced in the next major release:

noise.HERA_Tsky_mdl
noise.HERA_BEAM_POLY
sigchain.HERA_NRAO_BANDPASS
rfi.HERA_RFI_STATIONS

Additionally, the next major release will involve modifications to the package's API, which move toward a regularization of the way in which hera_sim methods are interfaced with; in particular, changes will be made such that the Simulator class is the most intuitive way of interfacing with the hera_sim package features.
  FutureWarning)

fqs = np.linspace(.1,.2,1024,endpoint=False)
lsts = np.linspace(0,2*np.pi,10000, endpoint=False)
times = lsts / (2*np.pi) * aipy.const.sidereal_day
bl_len_ns = np.array([30., 0, 0])

Foregrounds

Diffuse Foregrounds

Tsky_mdl = noise.HERA_Tsky_mdl['xx']
vis_fg_diffuse = foregrounds.diffuse_foreground(lsts, fqs, bl_len_ns, Tsky_mdl)

MX, DRNG = 2.5, 3
plt.figure()
plt.subplot(211); uvtools.plot.waterfall(vis_fg_diffuse, mode='log', mx=MX, drng=DRNG); plt.colorbar(); plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(vis_fg_diffuse, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

Point-Source Foregrounds

vis_fg_pntsrc = foregrounds.pntsrc_foreground(lsts, fqs, bl_len_ns, nsrcs=200)

plt.figure()
plt.subplot(211); uvtools.plot.waterfall(vis_fg_pntsrc, mode='log', mx=MX, drng=DRNG); plt.colorbar()#; plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(vis_fg_pntsrc, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

/home/bobby/anaconda3/envs/fix_tutorial/lib/python3.7/site-packages/uvtools/plot.py:40: RuntimeWarning: divide by zero encountered in log10
  data = np.log10(data)

Diffuse and Point-Source Foregrounds

vis_fg = vis_fg_diffuse + vis_fg_pntsrc

plt.figure()
plt.subplot(211); uvtools.plot.waterfall(vis_fg, mode='log', mx=MX, drng=DRNG); plt.colorbar(); plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(vis_fg, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

Noise

tsky = noise.resample_Tsky(fqs,lsts,Tsky_mdl=noise.HERA_Tsky_mdl['xx'])
t_rx = 150.
omega_p = noise.bm_poly_to_omega_p(fqs)
nos_jy = noise.sky_noise_jy(tsky + t_rx, fqs, lsts, omega_p)

plt.figure()
plt.subplot(211); uvtools.plot.waterfall(nos_jy, mode='log', mx=MX, drng=DRNG); plt.colorbar()#; plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(nos_jy, mode='phs'); plt.colorbar()#; plt.ylim(0,4000)
plt.show()

vis_fg_nos = vis_fg + nos_jy

plt.figure()
plt.subplot(211); uvtools.plot.waterfall(vis_fg_nos, mode='log', mx=MX, drng=DRNG); plt.colorbar(); plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(vis_fg_nos, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

RFI

rfi1 = rfi.rfi_stations(fqs, lsts)
rfi2 = rfi.rfi_impulse(fqs, lsts, chance=.02)
rfi3 = rfi.rfi_scatter(fqs, lsts, chance=.001)
rfi_all = rfi1 + rfi2 + rfi3

plt.figure()
plt.subplot(211); uvtools.plot.waterfall(rfi_all, mode='log', mx=MX, drng=DRNG); plt.colorbar(); plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(rfi_all, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

vis_fg_nos_rfi = vis_fg_nos + rfi_all

plt.figure()
plt.subplot(211); uvtools.plot.waterfall(vis_fg_nos_rfi, mode='log', mx=MX, drng=DRNG); plt.colorbar(); plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(vis_fg_nos_rfi, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

Gains

g = sigchain.gen_gains(fqs, [1,2,3])
plt.figure()
for i in g: plt.plot(fqs, np.abs(g[i]), label=str(i))
plt.legend(); plt.show()
gainscale = np.average([np.median(np.abs(g[i])) for i in g])
MXG = MX + np.log10(gainscale)

vis_total = sigchain.apply_gains(vis_fg_nos_rfi, g, (1,2))
plt.figure()
plt.subplot(211); uvtools.plot.waterfall(vis_total, mode='log', mx=MXG, drng=DRNG); plt.colorbar(); plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(vis_total, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

Crosstalk

xtalk = sigchain.gen_whitenoise_xtalk(fqs)
vis_xtalk = sigchain.apply_xtalk(vis_fg_nos_rfi, xtalk)
vis_xtalk = sigchain.apply_gains(vis_xtalk, g, (1,2))
plt.figure()
plt.subplot(211); uvtools.plot.waterfall(vis_xtalk, mode='log', mx=MXG, drng=DRNG); plt.colorbar(); plt.ylim(0,4000)
plt.subplot(212); uvtools.plot.waterfall(vis_xtalk, mode='phs'); plt.colorbar(); plt.ylim(0,4000)
plt.show()

The following tutorial will help you learn how use the Simulator class:

The Simulator Class

The most convenient way to use hera_sim is to use the Simulator class, which builds in all the primary functionality of the hera_sim package in an easy-to-use interface, and adds the ability to consistently write all produced effects into a pyuvdata.UVData object (and to file).

What follows is a quick tour of the main functionality this provides.

Setup

%matplotlib inline
import os
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u

import hera_sim
import uvtools

from hera_sim.simulate import Simulator
from hera_sim.noise import HERA_Tsky_mdl
from hera_sim.data import DATA_PATH
from hera_sim.interpolators import Beam

/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/visibilities/__init__.py:27: UserWarning: PRISim failed to import.
  warnings.warn("PRISim failed to import.")
/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/visibilities/__init__.py:33: UserWarning: VisGPU failed to import.
  warnings.warn("VisGPU failed to import.")
/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/__init__.py:36: FutureWarning:
In the next major release, all HERA-specific variables will be removed from the codebase. The following variables will need to be accessed through new class-like structures to be introduced in the next major release:

noise.HERA_Tsky_mdl
noise.HERA_BEAM_POLY
sigchain.HERA_NRAO_BANDPASS
rfi.HERA_RFI_STATIONS

Additionally, the next major release will involve modifications to the package's API, which move toward a regularization of the way in which hera_sim methods are interfaced with; in particular, changes will be made such that the Simulator class is the most intuitive way of interfacing with the hera_sim package features.
  FutureWarning)

The Simulator Class

The Simulator class holds everything required to generate and do basic analysis on a simulation. It can be instantiated by submitting any one of the following: a filename pointing to an existing simulation in uvfits, uvh5 or miriad format; a UVData object, or a set of keyword arguments which are passed to the hera_sim.io.empty_uvdata function. Any keyword arguments passed to the Simulator initializer are passed directly to the hera_sim.io.empty_uvdata, and those are in turn passed directly to the pyuvsim.simsetup.initialize_uvdata_from_keywords function, so that setting up a simulation object is relatively straightforward. The default parameter values for the empty_uvdata function are all set to None, but are given default values to show a suggested minimal set of parameters necessary to initialize a UVData object (and, in turn, a Simulator object):

hera_sim.io.empty_uvdata(
    Ntimes=None,
    start_time=None,
    integration_time=None,
    Nfreqs=None,
    start_freq=None,
    channel_width=None,
    array_layout=None,
    **kwargs,
)

If you are curious about the full set of parameters that may be used to initialize a UVData object, then please see the documentation for the pyuvsim.simsetup.initialize_uvdata_from_keywords function.

To get us started, let’s make a Simulator object with 100 frequency channels spanning from 100 to 200 MHz, a half-hour of observing time with an H1C-appropriate integration time of 10.7 seconds, and a 7-element hexagonal array.

# first, generate the array. this is returned as a dict
# with antenna numbers as keys and ENU positions as values
ants = hera_sim.antpos.hex_array(2, split_core=False, outriggers=0)

# for clarity about some timing parameters
start_time = 2458115.9 # JD
integration_time = 10.7 # seconds
Ntimes = int(30 * u.min.to("s") / integration_time)

sim = Simulator(
    Nfreqs=100,
    start_freq=1e8,
    bandwidth=1e8,
    Ntimes=Ntimes,
    start_time=start_time,
    integration_time=integration_time,
    array_layout=ants
)

The Simulator class adds some attributes for conveniently accessing metadata:

You can retrieve the frequencies (in GHz) with the freqs attribute. (sim.freqs)

You can retrieve the LSTs (in radians) with the lsts attribute. (sim.lsts)

You can retrieve the antenna array with the antpos attribute. (sim.antpos)

sim.freqs

array([0.1  , 0.101, 0.102, 0.103, 0.104, 0.105, 0.106, 0.107, 0.108,
       0.109, 0.11 , 0.111, 0.112, 0.113, 0.114, 0.115, 0.116, 0.117,
       0.118, 0.119, 0.12 , 0.121, 0.122, 0.123, 0.124, 0.125, 0.126,
       0.127, 0.128, 0.129, 0.13 , 0.131, 0.132, 0.133, 0.134, 0.135,
       0.136, 0.137, 0.138, 0.139, 0.14 , 0.141, 0.142, 0.143, 0.144,
       0.145, 0.146, 0.147, 0.148, 0.149, 0.15 , 0.151, 0.152, 0.153,
       0.154, 0.155, 0.156, 0.157, 0.158, 0.159, 0.16 , 0.161, 0.162,
       0.163, 0.164, 0.165, 0.166, 0.167, 0.168, 0.169, 0.17 , 0.171,
       0.172, 0.173, 0.174, 0.175, 0.176, 0.177, 0.178, 0.179, 0.18 ,
       0.181, 0.182, 0.183, 0.184, 0.185, 0.186, 0.187, 0.188, 0.189,
       0.19 , 0.191, 0.192, 0.193, 0.194, 0.195, 0.196, 0.197, 0.198,
       0.199])

sim.lsts

array([4.58108965, 4.58186991, 4.58265017, 4.58343042, 4.58421068,
       4.58499094, 4.58577119, 4.58655145, 4.58733171, 4.58811196,
       4.58889222, 4.58967247, 4.59045273, 4.59123299, 4.59201324,
       4.5927935 , 4.59357376, 4.59435401, 4.59513427, 4.59591452,
       4.59669478, 4.59747504, 4.59825529, 4.59903555, 4.59981581,
       4.60059606, 4.60137632, 4.60215657, 4.60293683, 4.60371709,
       4.60449734, 4.6052776 , 4.60605786, 4.60683811, 4.60761837,
       4.60839863, 4.60917888, 4.60995914, 4.6107394 , 4.61151965,
       4.61229991, 4.61308017, 4.61386042, 4.61464068, 4.61542093,
       4.61620119, 4.61698145, 4.6177617 , 4.61854196, 4.61932222,
       4.62010247, 4.62088273, 4.62166299, 4.62244324, 4.6232235 ,
       4.62400375, 4.62478401, 4.62556427, 4.62634452, 4.62712478,
       4.62790503, 4.62868529, 4.62946555, 4.6302458 , 4.63102606,
       4.63180632, 4.63258657, 4.63336683, 4.63414709, 4.63492734,
       4.6357076 , 4.63648786, 4.63726811, 4.63804837, 4.63882863,
       4.63960888, 4.64038914, 4.64116939, 4.64194965, 4.64272991,
       4.64351016, 4.64429042, 4.64507068, 4.64585093, 4.64663119,
       4.64741145, 4.6481917 , 4.64897196, 4.64975221, 4.65053247,
       4.65131273, 4.65209298, 4.65287324, 4.65365349, 4.65443375,
       4.65521401, 4.65599426, 4.65677452, 4.65755478, 4.65833503,
       4.65911529, 4.65989555, 4.6606758 , 4.66145606, 4.66223632,
       4.66301657, 4.66379683, 4.66457709, 4.66535734, 4.6661376 ,
       4.66691785, 4.66769811, 4.66847837, 4.66925862, 4.67003888,
       4.67081914, 4.67159939, 4.67237965, 4.67315991, 4.67394016,
       4.67472042, 4.67550068, 4.67628093, 4.67706119, 4.67784144,
       4.6786217 , 4.67940196, 4.68018221, 4.68096247, 4.68174272,
       4.68252298, 4.68330324, 4.68408349, 4.68486375, 4.68564401,
       4.68642426, 4.68720452, 4.68798478, 4.68876503, 4.68954529,
       4.69032555, 4.6911058 , 4.69188606, 4.69266631, 4.69344657,
       4.69422683, 4.69500708, 4.69578734, 4.6965676 , 4.69734785,
       4.69812811, 4.69890837, 4.69968862, 4.70046888, 4.70124914,
       4.70202939, 4.70280965, 4.7035899 , 4.70437016, 4.70515042,
       4.70593067, 4.70671093, 4.70749118, 4.70827144, 4.7090517 ,
       4.70983195, 4.71061221, 4.71139247])

sim.antpos

{0: array([-7.30000000e+00,  1.26439709e+01, -4.36185665e-09]),
 1: array([ 7.30000000e+00,  1.26439709e+01, -3.99203159e-09]),
 2: array([-1.46000000e+01,  6.98581573e-09, -4.65185394e-09]),
 3: array([ 0.00000000e+00,  7.20559015e-09, -4.28202888e-09]),
 4: array([ 1.46000000e+01,  7.42536457e-09, -3.91220382e-09]),
 5: array([-7.30000000e+00, -1.26439709e+01, -4.57202631e-09]),
 6: array([ 7.30000000e+00, -1.26439709e+01, -4.20220125e-09])}

You can also generate a plot of the array layout using the plot_array method:

fig = sim.plot_array()
plt.show()

The data attribute can be used to access the UVData object used to store the simulated data and metadata:

from pyuvdata import UVData
isinstance(sim.data, UVData)

True

We’ll use a standard waterfall plot throughout the notebook to show the progress of the simulation:

def waterfall(
    vis, freqs=sim.freqs * 1e3, lsts=sim.lsts,
    vmax=None, vrange=None, title=None,
):
    """
    A wrapper around the uvtools' waterfall function providing some
    extra labelling and plot adjustment.
    """
    fig, ax = plt.subplots(
        2,1, sharex=True, sharey=True, figsize=(12,10)
    )

    if title is not None:
        ax[0].set_title(title, fontsize=12)
    plt.sca(ax[0])
    uvtools.plot.waterfall(
        vis, mode='log', mx=vmax, drng=vrange,
        extent=(freqs.min(), freqs.max(), lsts.min(), lsts.max())
    )
    plt.colorbar(label=r'log$_{10}$(Vis/Jy)')
    plt.ylabel("LST", fontsize=12)

    plt.sca(ax[1])
    uvtools.plot.waterfall(
        vis,
        mode='phs',
        cmap='twilight',
        extent=(freqs.min(), freqs.max(), lsts.min(), lsts.max())
    )
    plt.colorbar(label='Phase [rad]')
    plt.xlabel("Frequency [MHz]", fontsize=12)
    plt.ylabel("LST", fontsize=12)

Adding Effects

Effects may be added to a simulation by using the add method. This method takes one argument and variable keyword arguments: the required argument component may be either a string identifying the name of a hera_sim class (or an alias thereof, see below), or a callable object that has either lsts or freqs (or both) as a parameter(s) and returns an object with shape (lsts.size, freqs.size) or (freqs.size,) (granted, the Simulator does not do a check on the input parameters or the shape of returned values, but passing a callable object that returns something incorrectly shaped will likely result in an exception being raised). Let’s walk through an example.

A Simple Example

Let’s begin by simulating diffuse foregrounds. To do this, we’ll need to specify a sky temperature model and a beam area interpolator object (or array of values corresponding to a frequency-dependent beam size). Let’s use the sky temperature model and beam polynomial fit for H1C:

from hera_sim.interpolators import Beam
from hera_sim.noise import HERA_Tsky_mdl
from hera_sim.data import DATA_PATH

# HERA_Tsky_mdl is a dictionary with keys 'xx' and 'yy'
Tsky_mdl = HERA_Tsky_mdl['xx']

# there are useful files in the data directory, with this being one
beamfile = os.path.join(DATA_PATH, "HERA_H1C_BEAM_POLY.npy")

# we can make a beam interpolation object from a path to a .npy file
omega_p = Beam(beamfile)

# now to add the foregrounds
sim.add("diffuse_foreground", Tsky_mdl=Tsky_mdl, omega_p=omega_p)

# let's check out the data for the (0,1) baseline
vis = sim.data.get_data((0,1,'xx'))
waterfall(vis, title='Simulated Diffuse Foregrounds')

Now, whenever an effect is simulated using the add method, that effect (and any non-default optional parameters used) is logged: a note of what was simulated is put into the object’s history, and the information is also stored in a hidden attribute _components.

print(sim.data.history)

hera_sim v1.0.0: Added DiffuseForeground using kwargs:
Tsky_mdl = <hera_sim.interpolators.Tsky object at 0x7fb870d28b50>
omega_p = <hera_sim.interpolators.Beam object at 0x7fb86dbd3390>

sim._components

{'diffuse_foreground': {'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7fb870d28b50>,
  'omega_p': <hera_sim.interpolators.Beam at 0x7fb86dbd3390>}}

Note that the actual visibility simulated is not cached, but rather the parameters necessary for recreating are saved and can be used to later re-simulate the effect. Now, the diffuse foreground simulator uses a random number generator, so if we want to ensure repeatability, we have to make sure it’s seeded in some repeatable way…

# refresh the simulation
# this zeros the data array, resets the history, and clears the _components dictionary
sim.refresh()

from hera_sim.foregrounds import DiffuseForeground
# now let's seed by redundant group
sim.add(DiffuseForeground, Tsky_mdl=Tsky_mdl, omega_p=omega_p, seed_mode="redundantly")

# and let's take a look
vis = sim.data.get_data((0,1,'xx'))
waterfall(vis, title='Simulated Diffuse Foregrounds')

# as for the history
print(sim.data.history)

hera_sim v1.0.0: Added DiffuseForeground using kwargs:
Tsky_mdl = <hera_sim.interpolators.Tsky object at 0x7fb870d28b50>
omega_p = <hera_sim.interpolators.Beam object at 0x7fb86dbd3390>
seed_mode = redundantly

# and the _components
sim._components

{hera_sim.foregrounds.DiffuseForeground: {'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7fb870d28b50>,
  'omega_p': <hera_sim.interpolators.Beam at 0x7fb86dbd3390>,
  'seed_mode': 'redundantly'}}

Note that using the DiffuseForeground class to simulate the effect caused the Simulator to log the component using the class itself as the key.

# we can even re-simulate the effect
vis_copy = sim.get(DiffuseForeground, ant1=0, ant2=1, pol='xx')
np.all(np.isclose(vis,vis_copy))

True

# but do baselines within a redundant group agree?
for reds in sim._get_reds():
    red_grp = [sim.data.baseline_to_antnums(red) for red in reds]
    print(list(red_grp))

[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6)]
[(0, 1), (2, 3), (3, 4), (5, 6)]
[(0, 2), (1, 3), (3, 5), (4, 6)]
[(0, 3), (1, 4), (2, 5), (3, 6)]
[(0, 4), (2, 6)]
[(0, 5), (1, 6)]
[(0, 6)]
[(1, 2), (4, 5)]
[(1, 5)]
[(2, 4)]

np.all(
    np.isclose(sim.data.get_data(0,1), sim.data.get_data(2,3))
)

True

A final note on this: when the seed_mode parameter is set to "redundantly", a seed is generated for each redundant group, and that seed is used to ensure that each baseline within a redundant group observes the same realization of the effect; the only other mode that’s currently supported is "once", which just seeds the random number generator once before any effects are simulated.

Filtering Effects

The add method also allows for an optional keyword argument vis_filter. This may be an iterable of arbitrary length, or an iterable of keys of various lengths. This is still under-developed and may have a different implementation in the near future, so this section will not go into any detail about how it works.

As a quick example, let’s choose our filter so that only the 'xx' polarization for the baseline (0,1) (and its conjugate) is simulated. To check that it works, we’ll test that the antennas’ autoc

sim.refresh()
sim.add("noiselike_eor", vis_filter=(0,1,'xx'))
np.all(sim.data.get_data(0,0,'xx') == 0)

True

np.all(sim.data.get_data(0,1,'xx') != 0) and np.all(sim.data.get_data(1,0,'xx') != 0)

True

Ordering of Effects

Some of the effects that can be simulated act on visibilities additively, whereas others act multiplicatively. The Simulator class keeps track of the order in which simulation components have been added, and issues warnings if the simulation components have been added out of order. In particular, it issues warnings if a multiplicative model is added to an empty set of visibilities or if an absolute visibility is added after a multiplicative effect has been added. Let’s check it out.

# first, let's add a multiplicative effect
gains = sim.add("gains", seed_mode="once", ret_vis=True)

Note that you can set the ret_vis parameter to True in order to have the add method return the effect it calculated. In this case, it’s a dictionary whose keys are antenna numbers and whose values are complex gains as a function of frequency.

# take a quick peek at the gains
fig = plt.figure(figsize=(11,7))
ax = fig.add_subplot(111)
ax.plot(sim.freqs * 1e3, np.abs(gains[0]))
ax.set_xlabel("Frequency [MHz]", fontsize=12)
ax.set_ylabel("Gain Amplitude", fontsize=12)
plt.show()

Now let’s try to add an additive effect, like EoR.

sim.add("noiselike_eor")

You are adding visibilities to a data array *after* multiplicative effects have been introduced.

Voila! The Simulator knows that a multiplicative effect (bandpass gains) have already been applied to the data and issues a warning when an additive effect is added afterward. Now let’s refresh the simulation and check what happens if we try to add gains immediately after.

sim.refresh()

# now try adding a multiplicative effect to a fresh simulation
gains = sim.add("gains", seed_mode="once", ret_vis=True)

You are trying to compute a multiplicative effect, but no visibilities have been simulated yet.

# the effect is still simulated, but of course nothing happens
bool(gains) and np.all(sim.data.data_array == 0)

True

# note, however, that this warning is only issued once
sim.add("gains")

Using Custom Components

You can also use custom-defined callable objects as simulation components. Here’s an example using a function:

# first, make up a function
def mock_vis(lsts, freqs):
    return np.ones((lsts.size, freqs.size), dtype=np.complex)

# then refresh the sim and add the effect
sim.refresh()
sim.add(mock_vis)

You are attempting to compute a component but have not specified an ``is_multiplicative`` attribute for the component. The component will be added under the assumption that it is *not* multiplicative.

# and now let's see what it's logged as
print(sim.data.history)

hera_sim v1.0.0: Added mock_vis using kwargs:

# as for the components dictionary
sim._components

{<function __main__.mock_vis(lsts, freqs)>: {}}

Note that a warning is raised if a user-defined function (or callable class) does not have an is_multiplicative attribute—this will be the case for all functions, but class implementations can have such an attribute defined.

Registering Classes

To wrap this section up, let’s briefly cover how the Simulator class knows what to look for when the add method is called. All of the simulation components defined in the hera_sim repo are designed as callable classes which inherit from a registry. A new registry can be created by defining a class and decorating it with the registry decorator:

from hera_sim import registry

@registry
class ExampleRegistry:
    is_multiplicative = False
    pass

class MockVis(ExampleRegistry):
    __aliases__ = ("mock_vis_alias",)
    def __call__(self, lsts, freqs):
        return np.ones((lsts.size, freqs.size), dtype=np.complex)

The above cell defines a new registry called ExampleRegistry and creates a class MockVis that is tracked by ExampleRegistry. Since the registry has its is_multiplicative attribute set to False, any class that inherits from it (and does not change the attribute’s value) will also be considered an additive model.

ExampleRegistry._models

{'MockVis': __main__.MockVis}

MockVis.is_multiplicative

False

Since the MockVis class has "mock_vis_alias" as one of the entries in its __aliases__ attribute, it can be discovered with that string:

sim.refresh()
sim.add("mock_vis_alias")

np.all(sim.data.data_array == 1)

True

print(sim.data.history)

hera_sim v1.0.0: Added MockVis using kwargs:

sim._components

{'mock_vis_alias': {}}

If you try to add a component that is not discoverable, then an exception will be raised and you’ll be shown a list of known aliases:

# show the message this way to let the rest of the kernel run
try:
    sim.add("not_a_known_model")
except AttributeError as err:
    print(err)

The component 'not_a_known_model' wasn't found. The following aliases are known: cross_coupling_xtalk, thermal_noise, noiselikeeor, gen_gains, gen_reflection_gains, impulse, gen_cross_coupling_xtalk, diffuseforeground, rfi_dtv, sigchain_reflections, rfi_impulse, diffuse_foreground, gen_whitenoise_xtalk, white_noise_xtalk, dtv, pointsourceforeground, hexarray, bandpass, gains, stations, rfi_stations, lineararray, reflections, mockvis, thermalnoise, whitenoisecrosstalk, scatter, mock_vis_alias, crosscouplingcrosstalk, bandpass_gain, noiselike_eor, pntsrc_foreground, rfi_scatter

If you would like to see what classes are registered and discoverable by the Simulator class, then you can use the list_discoverable_components function:

hera_sim.list_discoverable_components()

hera_sim.antpos.LinearArray
hera_sim.antpos.HexArray
hera_sim.foregrounds.DiffuseForeground
hera_sim.foregrounds.PointSourceForeground
hera_sim.noise.ThermalNoise
hera_sim.rfi.Stations
hera_sim.rfi.Impulse
hera_sim.rfi.Scatter
hera_sim.rfi.DTV
hera_sim.sigchain.Bandpass
hera_sim.sigchain.Reflections
hera_sim.sigchain.CrossCouplingCrosstalk
hera_sim.sigchain.WhiteNoiseCrosstalk
hera_sim.eor.NoiselikeEoR
__main__.MockVis

One final caveat: currently, the add method is implemented so that an exception is raised when a user-defined function that does not exist in the global namespace is passed as the component. So, if you want to import a simulation component that is implemented as a function from some other module, then be sure to import the function into the global namespace.

The run_sim Method

The Simulator class also features the run_sim method, which allows you to run an entire simulation by either specifying a dictionary whose keys are strings that specify simulation components (or acceptable aliases) compatible with the add method and whose values are the optional parameter settings for the simulation components.

An Example

Let’s keep the example short and sweet: diffuse foregrounds with bandpass gains.

# first, refresh the simulation to reset the history and _components
sim.refresh()

# now make a dictionary of simulation parameters
sim_params = {
    "diffuse_foreground" : {
        "Tsky_mdl" : Tsky_mdl, "omega_p" : omega_p, "seed_mode" : "redundantly"
    },
    "gains" : {"seed_mode" : "once"}
}

sim.run_sim(**sim_params)

# and let's inspect the history
print(sim.data.history)

hera_sim v1.0.0: Added DiffuseForeground using kwargs:
Tsky_mdl = <hera_sim.interpolators.Tsky object at 0x7fb870d28b50>
omega_p = <hera_sim.interpolators.Beam object at 0x7fb86dbd3390>
seed_mode = redundantly
hera_sim v1.0.0: Added Bandpass using kwargs:
seed_mode = once

We can also use a YAML file to control the simulation. Note the use of the YAML tags !Tsky and !Beam; these are defined in the __yaml_constructors module, and they use the Tsky and Beam classes, respectively, from the interpolators module to construct an interpolation object from the provided datafile (and optional interp_kwargs dictionary). Note that currently these classes assume that relative paths are intended to be relative to the hera_sim.data directory, but this behavior may change in the future.

# let's first make a configuration file
import tempfile
tempdir = tempfile.mkdtemp()
config = os.path.join(tempdir, "example.yaml")
with open(config, 'w') as cfg:
    cfg.write(
"""
diffuse_foreground:
    Tsky_mdl: !Tsky
        datafile: HERA_Tsky_Reformatted.npz
    omega_p: !Beam
        datafile: HERA_H1C_BEAM_POLY.npy
    seed_mode: redundantly
gains:
    seed_mode: once
"""
)

# now refresh the simulation and run the new simulation
sim.refresh()
sim.run_sim(config)

# again let's look at the history
print(sim.data.history)

hera_sim v1.0.0: Added DiffuseForeground using kwargs:
Tsky_mdl = <hera_sim.interpolators.Tsky object at 0x7fb8704b5cd0>
omega_p = <hera_sim.interpolators.Beam object at 0x7fb86e8598d0>
seed_mode = redundantly
hera_sim v1.0.0: Added Bandpass using kwargs:
seed_mode = once

# and the components
sim._components

{'diffuse_foreground': {'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7fb8704b5cd0>,
  'omega_p': <hera_sim.interpolators.Beam at 0x7fb86e8598d0>,
  'seed_mode': 'redundantly'},
 'gains': {'seed_mode': 'once'}}

Finally, you can return individual components from a simulation like so:

# return some things, but not all things
sim_params["diffuse_foreground"]["ret_vis"] = True
sim_params["pntsrc_foreground"] = {"seed_mode" : "redundantly", "ret_vis" : True}

# reorder so that multiplicative effects come last
new_sim_params = {"diffuse_foreground" : {}, "pntsrc_foreground" : {}, "gains" : {}}
new_sim_params.update(sim_params)

# refresh the simulation and run
sim.refresh()
results = dict(sim.run_sim(**new_sim_params))

# tada!
results.keys()

dict_keys(['diffuse_foreground', 'pntsrc_foreground'])

np.all(results['diffuse_foreground'].dtype == np.complex)

True

results['diffuse_foreground'].shape == sim.data.data_array.shape

True

np.all(
    np.isclose(results['diffuse_foreground'], sim.get("diffuse_foreground"))
)

True

sim.refresh()
sim.add("noiselike_eor", vis_filter=(0,1,'xx'))

The following tutorial will help you learn how to interface with the defaults module:

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 os
import numpy as np
import matplotlib.pyplot as plt
import yaml

import uvtools
import hera_sim
from hera_sim import Simulator
from hera_sim.data import DATA_PATH
from hera_sim.config import CONFIG_PATH
%matplotlib inline

/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/visibilities/__init__.py:27: UserWarning: PRISim failed to import.
  warnings.warn("PRISim failed to import.")
/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/visibilities/__init__.py:33: UserWarning: VisGPU failed to import.
  warnings.warn("VisGPU failed to import.")
/home/bobby/HERA/dev/fix_tutorial/hera_sim/hera_sim/__init__.py:36: FutureWarning:
In the next major release, all HERA-specific variables will be removed from the codebase. The following variables will need to be accessed through new class-like structures to be introduced in the next major release:

noise.HERA_Tsky_mdl
noise.HERA_BEAM_POLY
sigchain.HERA_NRAO_BANDPASS
rfi.HERA_RFI_STATIONS

Additionally, the next major release will involve modifications to the package's API, which move toward a regularization of the way in which hera_sim methods are interfaced with; in particular, changes will be made such that the Simulator class is the most intuitive way of interfacing with the hera_sim package features.
  FutureWarning)

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

# 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.noise._get_hera_bm_poly()
h1c_bandpass = hera_sim.sigchain._get_hera_bandpass()

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

h2c_beam = hera_sim.noise._get_hera_bm_poly()
h2c_bandpass = hera_sim.sigchain._get_hera_bandpass()

# now compare them
print("H1C Defaults:\n \nBeam Polynomial: \n{}\nBandpass Polynomial: \n{}\n".format(h1c_beam, h1c_bandpass))
print("H2C Defaults:\n \nBeam Polynomial: \n{}\nBandpass Polynomial: \n{}\n".format(h2c_beam, h2c_bandpass))

H1C Defaults:

Beam Polynomial:
[ 8.07774113e+08 -1.02194430e+09  5.59397878e+08 -1.72970713e+08
  3.30317669e+07 -3.98798031e+06  2.97189690e+05 -1.24980700e+04
  2.27220000e+02]
Bandpass Polynomial:
[-2.04689451e+06  1.90683718e+06 -7.41348361e+05  1.53930807e+05
 -1.79976473e+04  1.12270390e+03 -2.91166102e+01]

H2C Defaults:

Beam Polynomial:
[ 1.36744227e+13 -2.55445530e+13  2.14955504e+13 -1.07620674e+13
  3.56602626e+12 -8.22732117e+11  1.35321508e+11 -1.59624378e+10
  1.33794725e+09 -7.75754276e+07  2.94812713e+06 -6.58329699e+04
  6.52944619e+02]
Bandpass Polynomial:
[ 1.56076423e+17 -3.03924841e+17  2.72553042e+17 -1.49206626e+17
  5.56874144e+16 -1.49763003e+16  2.98853436e+15 -4.48609479e+14
  5.07747935e+13 -4.29965657e+12  2.67346077e+11 -1.18007726e+10
  3.48589690e+08 -6.15315646e+06  4.88747021e+04]

# another thing
fqs = np.linspace(0.1,0.2,1024)
lsts = np.linspace(0,2*np.pi,100)

noise = hera_sim.noise.thermal_noise

hera_sim.defaults.set('h1c')
np.random.seed(0)
h1c_noise = noise(fqs,lsts)

hera_sim.defaults.set('h2c')
np.random.seed(0)
h2c_noise = noise(fqs,lsts)
np.random.seed(0)
still_h2c_noise = noise(fqs,lsts)

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

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

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('HERA {} Beam Model'.format(seas.upper()), fontsize=12)
    ax.plot(fqs * 1e3, np.polyval(beams[seas], fqs))
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('HERA {} Bandpass Model'.format(seas.upper()), fontsize=12)
    ax.plot(fqs * 1e3, np.polyval(bps[seas], fqs))
plt.tight_layout()
plt.show()

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

{'foregrounds': {'diffuse_foreground': {'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7f52e274cd90>,
   'omega_p': <hera_sim.interpolators.Beam at 0x7f52e27561d0>}},
 'io': {'empty_uvdata': {'start_freq': 46920776.3671875,
   'channel_width': 122070.3125,
   'integration_time': 8.59}},
 'noise': {'_get_hera_bm_poly': {'bm_poly': 'HERA_H2C_BEAM_POLY.npy'},
  'resample_Tsky': {'Tsky': 180.0, 'mfreq': 0.18, 'index': -2.5},
  'sky_noise_jy': {'inttime': 8.59},
  'thermal_noise': {'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7f52e275e810>,
   'omega_p': <hera_sim.interpolators.Beam at 0x7f52e275eb10>,
   'Trx': 0,
   'inttime': 8.59}},
 'rfi': {'_get_hera_stations': {'rfi_stations': 'HERA_H2C_RFI_STATIONS.npy'},
  'rfi_impulse': {'chance': 0.001, 'strength': 20.0},
  'rfi_scatter': {'chance': 0.0001, 'strength': 10.0, 'std': 10.0},
  'rfi_dtv': {'freq_min': 0.174,
   'freq_max': 0.214,
   'width': 0.008,
   'chance': 0.0001,
   'strength': 10.0,
   'strength_std': 10.0}},
 'sigchain': {'_get_hera_bandpass': {'bandpass': 'HERA_H2C_BANDPASS.npy'},
  'gen_bandpass': {'gain_spread': 0.1},
  'gen_whitenoise_xtalk': {'amplitude': 3.0},
  'gen_cross_coupling_xtalk': {'amp': 0.0, 'dly': 0.0, 'phs': 0.0}}}

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

{'Tsky_mdl': <hera_sim.interpolators.Tsky at 0x7f52e275e810>,
 'omega_p': <hera_sim.interpolators.Beam at 0x7f52e275eb10>,
 'start_freq': 46920776.3671875,
 'channel_width': 122070.3125,
 'integration_time': 8.59,
 'bm_poly': 'HERA_H2C_BEAM_POLY.npy',
 'Tsky': 180.0,
 'mfreq': 0.18,
 'index': -2.5,
 'inttime': 8.59,
 'Trx': 0,
 'rfi_stations': 'HERA_H2C_RFI_STATIONS.npy',
 'chance': 0.0001,
 'strength': 10.0,
 'std': 10.0,
 'freq_min': 0.174,
 'freq_max': 0.214,
 'width': 0.008,
 'strength_std': 10.0,
 'bandpass': 'HERA_H2C_BANDPASS.npy',
 'gain_spread': 0.1,
 'amplitude': 3.0,
 'amp': 0.0,
 'dly': 0.0,
 'phs': 0.0}

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

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

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

waterfall(sim2, antpairpol)

The following tutorial will give you an overview of how to use hera_sim from the command line:

Running hera_sim from the command line

As of v0.2.0 of hera_sim, quick-and-dirty simulations can be run from the command line by creating a configuration file and using hera_sim’s run command to create simulated data in line with the configuration’s specifications. The basic syntax of using hera_sim’s command-line interface is:

$ hera_sim run --help
bash: line 1: hera_sim: command not found

An example configuration file can be found in the config_examples directory of the repo’s top-level directory. Here are its contents:

$ cat -n ../config_examples/template_config.yaml
     1	# This document is intended to serve as a template for constructing new
     2	# configuration YAMLs for use with the command-line interface.
     3	
     4	bda:
     5	    max_decorr: 0
     6	    pre_fs_int_time: !dimensionful
     7	        value: 0.1
     8	        units: 's'
     9	    corr_FoV_angle: !dimensionful
    10	        value: 20
    11	        units: 'deg'
    12	    max_time: !dimensionful
    13	        value: 16
    14	        units: 's'
    15	    corr_int_time: !dimensionful
    16	        value: 2
    17	        units: 's'
    18	filing:
    19	    outdir: '.'
    20	    outfile_name: 'quick_and_dirty_sim.uvh5'
    21	    output_format: 'uvh5'
    22	    clobber: True
    23	# freq and time entries currently configured for hera_sim use
    24	freq:
    25	    n_freq: 100
    26	    channel_width: 122070.3125
    27	    start_freq: 46920776.3671875
    28	time:
    29	    n_times: 10
    30	    integration_time: 8.59
    31	    start_time: 2457458.1738949567
    32	telescope:
    33	    # generate from an antenna layout csv
    34	    # array_layout: 'antenna_layout.csv'
    35	    # generate using hera_sim.antpos
    36	    array_layout: !antpos
    37	        array_type: "hex"
    38	        hex_num: 3
    39	        sep: 14.6
    40	        split_core: False
    41	        outriggers: 0
    42	    omega_p: !Beam
    43	        # non-absolute paths are assumed to be specified relative to the
    44	        # hera_sim data path
    45	        datafile: HERA_H2C_BEAM_MODEL.npz
    46	        interp_kwargs:
    47	            interpolator: interp1d
    48	            fill_value: extrapolate
    49	            # if you want to use a polynomial interpolator instead, then
    50	            # interpolator: poly1d
    51	            # kwargs not accepted for this; see numpy.poly1d documentation
    52	defaults:
    53	    # This must be a string specifying an absolute path to a default
    54	    # configuration file or one of the season default keywords
    55	    'h2c'
    56	systematics:
    57	    rfi:
    58	        # see hera_sim.rfi documentation for details on parameter names
    59	        rfi_stations:
    60	            seed: once
    61	            stations: !!null
    62	        rfi_impulse:
    63	            impulse_chance: 0.001
    64	            impulse_strength: 20.0
    65	        rfi_scatter:
    66	            scatter_chance: 0.0001
    67	            scatter_strength: 10.0
    68	            scatter_std: 10.0
    69	        rfi_dtv:
    70	            seed: once
    71	            dtv_band: 
    72	                - 0.174
    73	                - 0.214
    74	            dtv_channel_width: 0.008
    75	            dtv_chance: 0.0001
    76	            dtv_strength: 10.0
    77	            dtv_std: 10.0
    78	    sigchain:
    79	        gains:
    80	            seed: once
    81	            gain_spread: 0.1
    82	            dly_rng: [-20, 20]
    83	            bp_poly: HERA_H1C_BANDPASS.npy
    84	        sigchain_reflections:
    85	            seed: once
    86	            amp: !!null
    87	            dly: !!null
    88	            phs: !!null
    89	    crosstalk:
    90	        # only one of the two crosstalk methods should be specified
    91	        gen_whitenoise_xtalk:
    92	            amplitude: 3.0
    93	        # gen_cross_coupling_xtalk:
    94	            # seed: initial
    95	            # amp: !!null
    96	            # dly: !!null
    97	            # phs: !!null
    98	    noise:
    99	        thermal_noise:
   100	            seed: initial
   101	            Trx: 0
   102	sky:
   103	    Tsky_mdl: !Tsky
   104	        # non-absolute paths are assumed to be relative to the hera_sim
   105	        # data folder
   106	        datafile: HERA_Tsky_Reformatted.npz
   107	        # interp kwargs are passed to scipy.interp.RectBivariateSpline
   108	        interp_kwargs:
   109	            pol: xx # this is popped when making a Tsky object
   110	    eor:
   111	        noiselike_eor:
   112	            eor_amp: 0.00001
   113	            min_delay: !!null
   114	            max_delay: !!null
   115	            seed: redundant # so redundant baselines see same sky
   116	            fringe_filter_type: tophat
   117	    foregrounds:
   118	        # if using hera_sim.foregrounds
   119	        diffuse_foreground:
   120	            seed: redundant # redundant baselines see same sky
   121	            delay_filter_kwargs:
   122	                standoff: 0
   123	                delay_filter_type: tophat
   124	                normalize: !!null
   125	            fringe_filter_kwargs:
   126	                fringe_filter_type: tophat
   127	        pntsrc_foreground:
   128	            seed: once
   129	            nsrcs: 1000
   130	            Smin: 0.3
   131	            Smax: 300
   132	            beta: -1.5
   133	            spectral_index_mean: -1.0
   134	            spectral_index_std: 0.5
   135	            reference_freq: 0.5
   136	        # Note regarding seed_redundantly:
   137	        # This ensures that baselines within a redundant group see the same sky;
   138	        # however, this does not ensure that the sky is actually consistent. So,
   139	        # while the data produced can be absolutely calibrated, it cannot be
   140	        # used to make sensible images (at least, I don't *think* it can be).
   141	
   142	simulation:
   143	    # specify which components to simulate in desired order
   144	    # this should be a complete list of the things to include if hera_sim
   145	    # is the simulator being used. this will necessarily look different
   146	    # if other simulators are used, but that's not implemented yet
   147	    #
   148	    components: [foregrounds,
   149	                 noise,
   150	                 eor,
   151	                 rfi,
   152	                 sigchain, ]
   153	    # list particular model components to exclude from simulation
   154	    exclude: [sigchain_reflections,
   155	              gen_whitenoise_xtalk,]

The remainder of this tutorial will be spent on exploring each of the items in the above configuration file.

BDA

The following block of text shows all of the options that must be specified if you would like to apply BDA to the simulated data. Note that BDA is applied at the very end of the script, and requires the BDA package to be installed from http://github.com/HERA-Team/baseline_dependent_averaging.

$ sed -n 4,17p ../config_examples/template_config.yaml
bda:
    max_decorr: 0
    pre_fs_int_time: !dimensionful
        value: 0.1
        units: 's'
    corr_FoV_angle: !dimensionful
        value: 20
        units: 'deg'
    max_time: !dimensionful
        value: 16
        units: 's'
    corr_int_time: !dimensionful
        value: 2
        units: 's'

Please refer to the bda.apply_bda documentation for details on what each parameter represents. Note that practically each entry has the tag !dimensionful; this YAML tag converts the entries in value and units to an astropy.units.quantity.Quantity object with the specified value and units.

Filing

The following block of text shows all of the options that may be specified in the filing section; however, not all of these must be specified. In fact, the only parameter that is required to be specified in the config YAML is output_format, and it must be either miriad, uvfits, or uvh5. These are currently the only supported write methods for UVData objects.

$ sed -n 18,24p ../config_examples/template_config.yaml
filing:
    outdir: '.'
    outfile_name: 'quick_and_dirty_sim.uvh5'
    output_format: 'uvh5'
    clobber: True
# freq and time entries currently configured for hera_sim use
freq:

Recall that run can be called with the option --outfile; this specifies the full path to where the simulated data should be saved and overrides the outdir and outfile_name settings from the config YAML. Additionally, one can choose to use the flag -c or --clobber in place of specifying clobber in the config YAML. Finally, the dictionary defined by the kwargs entry has its contents passed to whichever write method is chosen, and the save_seeds option should only be used if the seed_redundantly option is specified for any of the simulation components.

Setup

The following block of text contains three sections: freq, time, and telescope. These sections are used to initialize the Simulator object that is used to perform the simulation. Note that the config YAML shows all of the options that may be specified, but not all options are necessarily required.

$ sed -n 26,53p ../config_examples/template_config.yaml
    channel_width: 122070.3125
    start_freq: 46920776.3671875
time:
    n_times: 10
    integration_time: 8.59
    start_time: 2457458.1738949567
telescope:
    # generate from an antenna layout csv
    # array_layout: 'antenna_layout.csv'
    # generate using hera_sim.antpos
    array_layout: !antpos
        array_type: "hex"
        hex_num: 3
        sep: 14.6
        split_core: False
        outriggers: 0
    omega_p: !Beam
        # non-absolute paths are assumed to be specified relative to the
        # hera_sim data path
        datafile: HERA_H2C_BEAM_MODEL.npz
        interp_kwargs:
            interpolator: interp1d
            fill_value: extrapolate
            # if you want to use a polynomial interpolator instead, then
            # interpolator: poly1d
            # kwargs not accepted for this; see numpy.poly1d documentation
defaults:
    # This must be a string specifying an absolute path to a default

If you are familiar with using configuration files with pyuvsim, then you’ll notice that the sections shown above look very similar to the way config files are constructed for use with pyuvsim. The config files for run were designed as an extension of the pyuvsim config files, with the caveat that some of the naming conventions used in pyuvsim are somewhat different than those used in hera_sim. For information on the parameters listed in the freq and time sections, please refer to the documentation for hera_sim.io.empty_uvdata. As for the telescope section, this is where the antenna array and primary beam are defined. The array_layout entry specifies the array, either by specifying an antenna layout file or by using the !antpos YAML tag and specifying the type of array (currently only linear and hex are supported) and the parameters to be passed to the corresponding function in hera_sim.antpos. The omega_p entry is where the primary beam is specified, and it is currently assumed that the beam is the same for each simulation component (indeed, this simulator is not intended to produce super-realistic simulations, but rather perform simulations quickly and give somewhat realistic results). This entry defines an interpolation object to be used for various hera_sim functions which require such an object; please refer to the documentation for hera_sim.interpolators.Beam for more information. Future versions of hera_sim will provide support for specifying the beam in an antenna layout file, similar to how it is done by pyuvsim.

Defaults

This section of the configuration file is optional to include. This section gives the user the option to use a default configuration to specify different parameters throughout the codebase. Users may define their own default configuration files, or they may use one of the provided season default configurations, located in the config folder. The currently supported season configurations are h1c and h2c. Please see the defaults module/documentation for more information.

$ sed -n 54,57p ../config_examples/template_config.yaml
    # configuration file or one of the season default keywords
    'h2c'
systematics:
    rfi:

Systematics

This is the section where any desired systematic effects can be specified. The block of text shown below details all of the possible options for systematic effects. Note that currently the sigchain_reflections and gen_cross_coupling_xtalk sections cannot easily be worked with; in fact, gen_cross_coupling_xtalk does not work as intended (each baseline has crosstalk show up at the same phase and delay, with the same amplitude, but uses a different autocorrelation visibility). Also note that the rfi section is subject to change, pending a rework of the rfi module.

$ sed -n 58,96p ../config_examples/template_config.yaml
        # see hera_sim.rfi documentation for details on parameter names
        rfi_stations:
            seed: once
            stations: !!null
        rfi_impulse:
            impulse_chance: 0.001
            impulse_strength: 20.0
        rfi_scatter:
            scatter_chance: 0.0001
            scatter_strength: 10.0
            scatter_std: 10.0
        rfi_dtv:
            seed: once
            dtv_band: 
                - 0.174
                - 0.214
            dtv_channel_width: 0.008
            dtv_chance: 0.0001
            dtv_strength: 10.0
            dtv_std: 10.0
    sigchain:
        gains:
            seed: once
            gain_spread: 0.1
            dly_rng: [-20, 20]
            bp_poly: HERA_H1C_BANDPASS.npy
        sigchain_reflections:
            seed: once
            amp: !!null
            dly: !!null
            phs: !!null
    crosstalk:
        # only one of the two crosstalk methods should be specified
        gen_whitenoise_xtalk:
            amplitude: 3.0
        # gen_cross_coupling_xtalk:
            # seed: initial
            # amp: !!null
            # dly: !!null

Note that although these simulation components are listed under systematics, they do not necessarily need to be listed here; the configuration file is formatted as such just for semantic clarity. For information on any particular simulation component listed here, please refer to the corresponding function’s documentation. For those who may not know what it means, !!null is how NoneType objects are specified using pyyaml.

Sky

This section specifies both the sky temperature model to be used throughout the simulation as well as any simulation components which are best interpreted as being associated with the sky (rather than as a systematic effect). Just like the systematics section, these do not necessarily need to exist in the sky section (however, the Tsky_mdl entry must be placed in this section, as that’s where the script looks for it).

$ sed -n 97,130p ../config_examples/template_config.yaml
            # phs: !!null
    noise:
        thermal_noise:
            seed: initial
            Trx: 0
sky:
    Tsky_mdl: !Tsky
        # non-absolute paths are assumed to be relative to the hera_sim
        # data folder
        datafile: HERA_Tsky_Reformatted.npz
        # interp kwargs are passed to scipy.interp.RectBivariateSpline
        interp_kwargs:
            pol: xx # this is popped when making a Tsky object
    eor:
        noiselike_eor:
            eor_amp: 0.00001
            min_delay: !!null
            max_delay: !!null
            seed: redundant # so redundant baselines see same sky
            fringe_filter_type: tophat
    foregrounds:
        # if using hera_sim.foregrounds
        diffuse_foreground:
            seed: redundant # redundant baselines see same sky
            delay_filter_kwargs:
                standoff: 0
                delay_filter_type: tophat
                normalize: !!null
            fringe_filter_kwargs:
                fringe_filter_type: tophat
        pntsrc_foreground:
            seed: once
            nsrcs: 1000
            Smin: 0.3

As of now, run only supports simulating effects using the functions in hera_sim; however, we intend to provide support for using different simulators in the future. If you would like more information regarding the Tsky_mdl entry, please refer to the documentation for the hera_sim.interpolators.Tsky class. Finally, note that the seed_redundantly parameter is specified for each entry in eor and foregrounds; this parameter is used to ensure that baselines within a redundant group all measure the same visibility, which is a necessary feature for data to be absolutely calibrated. Please refer to the documentation for hera_sim.eor and hera_sim.foregrounds for more information on the parameters and functions listed above.

Simulation

This section is used to specify which of the simulation components to include in or exclude from the simulation. There are only two entries in this section: components and exclude. The components entry should be a list specifying which of the groups from the sky and systematics sections should be included in the simulation. The exclude entry should be a list specifying which of the particular models should not be simulated. Here’s an example:

$ sed -n -e 137,138p -e 143,150p ../config_examples/template_config.yaml
        # This ensures that baselines within a redundant group see the same sky;
        # however, this does not ensure that the sky is actually consistent. So,
    # specify which components to simulate in desired order
    # this should be a complete list of the things to include if hera_sim
    # is the simulator being used. this will necessarily look different
    # if other simulators are used, but that's not implemented yet
    #
    components: [foregrounds,
                 noise,
                 eor,

The entries listed above would result in a simulation that includes all models contained in the foregrounds, noise, eor, rfi, and sigchain dictionaries, except for the sigchain_reflections and gen_whitenoise_xtalk models. So the simulation would consist of diffuse and point source foregrounds, thermal noise, noiselike EoR, all types of RFI modeled by hera_sim, and bandpass gains, with the effects simulated in that order. It is important to make sure that effects which enter multiplicatively (i.e. models from sigchain) are simulated after effects that enter additively, since the order that the simulation components are listed in is the same as the order of execution.

API Reference

Modules

hera_sim.vis
hera_sim.antpos	A module for creating antenna array configurations.
hera_sim.defaults	Class for dynamically changing hera_sim parameter defaults.
hera_sim.eor	A module for simulating EoR-like visibilities.
hera_sim.foregrounds	Reimagining of the foregrounds module, using an object-oriented approach.
hera_sim.interpolators	This module provides interfaces to different interpolation classes.
hera_sim.io
hera_sim.noise	Make some noise.
hera_sim.rfi	v1 implementation of RFI.
hera_sim.sigchain	Object-oriented approach to signal chain systematics.
hera_sim.simulate	Re-imagining of the simulation module.
hera_sim.utils	Utility module

Visibility Simulators

hera_sim.visibilities.simulators
hera_sim.visibilities.conversions	A number of mappings which may be useful for visibility simulators.
hera_sim.visibilities.vis_cpu

Contributing

Contributions are welcome, and they are greatly appreciated! Every little bit helps, and credit will always be given.

Bug reports

When reporting a bug please include:

• Your operating system name and version.

• Any details about your local setup that might be helpful in troubleshooting.

• Detailed steps to reproduce the bug.

Documentation improvements

hera_sim could always use more documentation, whether as part of the official hera_sim docs or in docstrings.

Feature requests and feedback

The best way to send feedback is to file an issue at https://github.com/HERA-Team/hera_sim/issues.

If you are proposing a feature:

• Explain in detail how it would work.

• Keep the scope as narrow as possible, to make it easier to implement.

• Remember that this is a volunteer-driven project, and that code contributions are welcome :)

Development

To set up hera_sim for local development:

1. If you are not on the HERA-Team collaboration, make a fork of hera_sim (look for the “Fork” button).

2. Clone the repository locally. If you made a fork in step 1:

   git clone git@github.com:your_name_here/hera_sim.git
   

   Otherwise:

   git clone git@github.com:HERA-Team/hera_sim.git
   

3. Create a branch for local development (you will not be able to push to “master”):

   git checkout -b name-of-your-bugfix-or-feature
   

   Now you can make your changes locally.

4. Make a development environment. We highly recommend using conda for this. With conda, just run:

   conda env create -f travis-environment.yml
   

4. When you’re done making changes, run all the checks, doc builder and spell checker with tox one command:

   tox
   

5. Commit your changes and push your branch to GitHub:

   git add .
   git commit -m "Your detailed description of your changes."
   git push origin name-of-your-bugfix-or-feature
   

6. Submit a pull request through the GitHub website.

Pull Request Guidelines

If you need some code review or feedback while you’re developing the code just make the pull request.

For merging, you should:

1. Include passing tests (run tox) 1.

2. Update documentation when there’s new API, functionality etc.

3. Add a note to CHANGELOG.rst about the changes.

4. Add yourself to AUTHORS.rst.

1

If you don’t have all the necessary python versions available locally you can rely on Travis - it will run the tests for each change you add in the pull request.

It will be slower though …

Developing hera_sim

hera_sim broadly follows the best-practices laid out in XXX.

Todo

where is that best-practices doc?

All docstrings should be written in Google docstring format.

AUTHORS

• HERA-Team - https://github.com/HERA-Team

Changelog

v1.0.0 [???]

Added

Fixed

Changed

• All functions that take frequencies and LSTs as arguments have had their signatures changed to func(lsts, freqs, *args, **kwargs)

• Changes to handling of functions which employ a fringe or delay filtering step with variable keywords for the filters used. Filter keywords are now specified with individual dictionaries called delay_filter_kwargs or fringe_filter_kwargs, depending on the filter used.

  • Functions affected by this change are as follows:

    • diffuse_foreground

    • noiselike_eor

• Changes to parameters that shared the same name but represented conceptually different objects in various functions. Functions affected by this are as follows:

  • utils.rough_delay_filter

  • utils.rough_fringe_filter

  • Most RFI functions

• The io.empty_uvdata function had its default keyword values set to None. The keywords accepted by this function have also been changed to match their names in pyuvsim.simsetup.initialize_uvdata_from_keywords

• Changes to parameters in most RFI models. Optional parameters that are common to many models (but should not share the same name), such as std or chance, have been prefixed with the model name and an underscore (e.g. dtv_chance). Various other parameters have been renamed for consistency/clarity. Note that the freq_min and freq_max parameters for rfi_dtv have been replaced by the single parameter dtv_band, which is a tuple denoting the edges of the DTV band in GHz.

• Functions in utils now use freqs instead of fqs.

v0.3.0 [2019.12.10]

Added

• New sub-package simulators

  • VisibilitySimulators class

    • Provides a common interface to interferometric visibility simulators. Users instantiate one of its subclasses and provide input antenna and sky scenarios.

    • HealVis subclass

    • Provides an interface to the healvis visibility simulator.

  • VisCPU subclass

    • Provides an interface to the viscpu visibility simulator.

  • conversions module

    • Not intended to be interfaced with by the end user; it provides useful coordinate transformations for VisibilitySimulators.

v0.2.0 [2019.11.20]

Added

• Command-line Interface

  • Use anywhere with hera_sim run [options] INPUT

  • Tutorial available on readthedocs

• Enhancement of run_sim method of Simulator class

  • Allows for each simulation component to be returned

    • Components returned as a list of 2-tuples (model_name, visibility)

    • Components returned by specifying ret_vis=True in their kwargs

• Option to seed random number generators for various methods

  • Available via the Simulator.add_ methods by specifying the kwarg seed_redundantly=True

  • Seeds are stored in Simulator object, and may be saved as a npy file when using the Simulator.write_data method

• New YAML tag !antpos

  • Allows for antenna layouts to be constructed using hera_sim.antpos functions by specifying parameters in config file

Fixed

• Changelog formatting for v0.1.0 entry

Changed

• Implementation of defaults module

  • Allows for semantic organization of config files

  • Parameters that have the same name take on the same value

    • e.g. std in various rfi functions only has one value, even if it’s specified multiple times

v0.1.0 [2019.08.28]

Added

• New module interpolators

  • Classes intended to be interfaced with by end-users:

    • Tsky

      • Provides an interface for generating a sky temperature interpolation object when provided with a .npz file and interpolation kwargs.

    • Beam, Bandpass

      • Provides an interface for generating either a poly1d or interp1d interpolation object when provided with an appropriate datafile.

• New module defaults

  • Provides an interface which allows the user to dynamically adjust default parameter settings for various hera_sim functions.

• New module __yaml_constructors

  • Not intended to be interfaced with by the end user; this module just provides a location for defining new YAML tags to be used in conjunction with the defaults module features and the Simulator.run_sim method.

• New directory config

  • Provides a location to store configuration files.

Fixed

Changed

• HERA-specific variables had their definitions removed from the codebase. Objects storing these variables still exist in the codebase, but their definitions now come from loading in data stored in various new files added to the data directory.

v0.0.1

• Initial released version

Indices and tables

• Index

• Module Index

• Search Page